Actual source code: bddcnullspace.c
1: #include <petsc/private/pcbddcimpl.h>
2: #include <petsc/private/pcbddcprivateimpl.h>
3: #include <../src/mat/impls/dense/seq/dense.h>
5: /* E + small_solve */
6: static PetscErrorCode PCBDDCNullSpaceCorrPreSolve(KSP ksp, Vec y, Vec x, void *ctx)
7: {
8: NullSpaceCorrection_ctx corr_ctx = (NullSpaceCorrection_ctx)ctx;
9: Mat K;
11: PetscFunctionBegin;
12: PetscCall(PetscLogEventBegin(corr_ctx->evapply, ksp, 0, 0, 0));
13: PetscCall(MatMultTranspose(corr_ctx->basis_mat, y, corr_ctx->sw[0]));
14: if (corr_ctx->symm) {
15: PetscCall(MatMult(corr_ctx->inv_smat, corr_ctx->sw[0], corr_ctx->sw[1]));
16: } else {
17: PetscCall(MatMultTranspose(corr_ctx->inv_smat, corr_ctx->sw[0], corr_ctx->sw[1]));
18: }
19: PetscCall(VecScale(corr_ctx->sw[1], -1.0));
20: PetscCall(MatMult(corr_ctx->basis_mat, corr_ctx->sw[1], corr_ctx->fw[0]));
21: PetscCall(VecScale(corr_ctx->sw[1], -1.0));
22: PetscCall(KSPGetOperators(ksp, &K, NULL));
23: PetscCall(MatMultAdd(K, corr_ctx->fw[0], y, y));
24: PetscCall(PetscLogEventEnd(corr_ctx->evapply, ksp, 0, 0, 0));
25: PetscFunctionReturn(PETSC_SUCCESS);
26: }
28: /* E^t + small */
29: static PetscErrorCode PCBDDCNullSpaceCorrPostSolve(KSP ksp, Vec y, Vec x, void *ctx)
30: {
31: NullSpaceCorrection_ctx corr_ctx = (NullSpaceCorrection_ctx)ctx;
32: Mat K;
34: PetscFunctionBegin;
35: PetscCall(PetscLogEventBegin(corr_ctx->evapply, ksp, 0, 0, 0));
36: PetscCall(KSPGetOperators(ksp, &K, NULL));
37: if (corr_ctx->symm) {
38: PetscCall(MatMult(K, x, corr_ctx->fw[0]));
39: } else {
40: PetscCall(MatMultTranspose(K, x, corr_ctx->fw[0]));
41: }
42: PetscCall(MatMultTranspose(corr_ctx->basis_mat, corr_ctx->fw[0], corr_ctx->sw[0]));
43: PetscCall(VecScale(corr_ctx->sw[0], -1.0));
44: PetscCall(MatMult(corr_ctx->inv_smat, corr_ctx->sw[0], corr_ctx->sw[2]));
45: PetscCall(MatMultAdd(corr_ctx->basis_mat, corr_ctx->sw[2], x, corr_ctx->fw[0]));
46: PetscCall(VecScale(corr_ctx->fw[0], corr_ctx->scale));
47: /* Sum contributions from approximate solver and projected system */
48: PetscCall(MatMultAdd(corr_ctx->basis_mat, corr_ctx->sw[1], corr_ctx->fw[0], x));
49: PetscCall(PetscLogEventEnd(corr_ctx->evapply, ksp, 0, 0, 0));
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: static PetscErrorCode PCBDDCNullSpaceCorrDestroy(void *ctx)
54: {
55: NullSpaceCorrection_ctx corr_ctx = (NullSpaceCorrection_ctx)ctx;
57: PetscFunctionBegin;
58: PetscCall(VecDestroyVecs(3, &corr_ctx->sw));
59: PetscCall(VecDestroyVecs(1, &corr_ctx->fw));
60: PetscCall(MatDestroy(&corr_ctx->basis_mat));
61: PetscCall(MatDestroy(&corr_ctx->inv_smat));
62: PetscCall(PetscFree(corr_ctx));
63: PetscFunctionReturn(PETSC_SUCCESS);
64: }
66: PetscErrorCode PCBDDCNullSpaceAssembleCorrection(PC pc, PetscBool isdir, PetscBool needscaling)
67: {
68: PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
69: MatNullSpace NullSpace = NULL;
70: KSP local_ksp;
71: NullSpaceCorrection_ctx shell_ctx;
72: Mat local_mat, local_pmat, dmat, Kbasis_mat;
73: Vec v;
74: PetscContainer c;
75: PetscInt basis_size;
76: IS zerorows;
77: PetscBool iscusp;
79: PetscFunctionBegin;
80: if (isdir) local_ksp = pcbddc->ksp_D; /* Dirichlet solver */
81: else local_ksp = pcbddc->ksp_R; /* Neumann solver */
82: PetscCall(KSPGetOperators(local_ksp, &local_mat, &local_pmat));
83: PetscCall(MatGetNearNullSpace(local_pmat, &NullSpace));
84: if (!NullSpace) {
85: if (pcbddc->dbg_flag) PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d doesn't have local (near) nullspace: no need for correction in %s solver \n", PetscGlobalRank, isdir ? "Dirichlet" : "Neumann"));
86: PetscFunctionReturn(PETSC_SUCCESS);
87: }
88: PetscCall(PetscObjectQuery((PetscObject)NullSpace, "_PBDDC_Null_dmat", (PetscObject *)&dmat));
89: PetscCheck(dmat, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing dense matrix");
90: PetscCall(PetscLogEventBegin(PC_BDDC_ApproxSetUp[pcbddc->current_level], pc, 0, 0, 0));
92: PetscCall(PetscNew(&shell_ctx));
93: shell_ctx->scale = 1.0;
94: PetscCall(PetscObjectReference((PetscObject)dmat));
95: shell_ctx->basis_mat = dmat;
96: PetscCall(MatGetSize(dmat, NULL, &basis_size));
97: shell_ctx->evapply = PC_BDDC_ApproxApply[pcbddc->current_level];
99: PetscCall(MatIsSymmetric(local_mat, 0.0, &shell_ctx->symm));
101: /* explicit construct (Phi^T K Phi)^-1 */
102: PetscCall(PetscObjectTypeCompare((PetscObject)local_mat, MATSEQAIJCUSPARSE, &iscusp));
103: if (iscusp) PetscCall(MatConvert(shell_ctx->basis_mat, MATSEQDENSECUDA, MAT_INPLACE_MATRIX, &shell_ctx->basis_mat));
104: PetscCall(MatMatMult(local_mat, shell_ctx->basis_mat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Kbasis_mat));
105: PetscCall(MatTransposeMatMult(Kbasis_mat, shell_ctx->basis_mat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &shell_ctx->inv_smat));
106: PetscCall(MatDestroy(&Kbasis_mat));
107: PetscCall(MatBindToCPU(shell_ctx->inv_smat, PETSC_TRUE));
108: PetscCall(MatFindZeroRows(shell_ctx->inv_smat, &zerorows));
109: if (zerorows) { /* linearly dependent basis */
110: const PetscInt *idxs;
111: PetscInt i, nz;
113: PetscCall(ISGetLocalSize(zerorows, &nz));
114: PetscCall(ISGetIndices(zerorows, &idxs));
115: for (i = 0; i < nz; i++) PetscCall(MatSetValue(shell_ctx->inv_smat, idxs[i], idxs[i], 1.0, INSERT_VALUES));
116: PetscCall(ISRestoreIndices(zerorows, &idxs));
117: PetscCall(MatAssemblyBegin(shell_ctx->inv_smat, MAT_FINAL_ASSEMBLY));
118: PetscCall(MatAssemblyEnd(shell_ctx->inv_smat, MAT_FINAL_ASSEMBLY));
119: }
120: PetscCall(MatLUFactor(shell_ctx->inv_smat, NULL, NULL, NULL));
121: PetscCall(MatSeqDenseInvertFactors_Private(shell_ctx->inv_smat));
122: if (zerorows) { /* linearly dependent basis */
123: const PetscInt *idxs;
124: PetscInt i, nz;
126: PetscCall(ISGetLocalSize(zerorows, &nz));
127: PetscCall(ISGetIndices(zerorows, &idxs));
128: for (i = 0; i < nz; i++) PetscCall(MatSetValue(shell_ctx->inv_smat, idxs[i], idxs[i], 0.0, INSERT_VALUES));
129: PetscCall(ISRestoreIndices(zerorows, &idxs));
130: PetscCall(MatAssemblyBegin(shell_ctx->inv_smat, MAT_FINAL_ASSEMBLY));
131: PetscCall(MatAssemblyEnd(shell_ctx->inv_smat, MAT_FINAL_ASSEMBLY));
132: }
133: PetscCall(ISDestroy(&zerorows));
135: /* Create work vectors in shell context */
136: PetscCall(MatCreateVecs(shell_ctx->inv_smat, &v, NULL));
137: PetscCall(KSPCreateVecs(local_ksp, 1, &shell_ctx->fw, 0, NULL));
138: PetscCall(VecDuplicateVecs(v, 3, &shell_ctx->sw));
139: PetscCall(VecDestroy(&v));
141: /* add special pre/post solve to KSP (see [1], eq. 48) */
142: PetscCall(KSPSetPreSolve(local_ksp, PCBDDCNullSpaceCorrPreSolve, shell_ctx));
143: PetscCall(KSPSetPostSolve(local_ksp, PCBDDCNullSpaceCorrPostSolve, shell_ctx));
144: PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)local_ksp), &c));
145: PetscCall(PetscContainerSetPointer(c, shell_ctx));
146: PetscCall(PetscContainerSetUserDestroy(c, PCBDDCNullSpaceCorrDestroy));
147: PetscCall(PetscObjectCompose((PetscObject)local_ksp, "_PCBDDC_Null_PrePost_ctx", (PetscObject)c));
148: PetscCall(PetscContainerDestroy(&c));
150: /* Create ksp object suitable for extreme eigenvalues' estimation */
151: if (needscaling || pcbddc->dbg_flag) {
152: KSP check_ksp;
153: PC local_pc;
154: Vec work1, work2;
155: const char *prefix;
156: PetscReal test_err, lambda_min, lambda_max;
157: PetscInt k, maxit;
158: PetscBool isspd, isset;
160: PetscCall(VecDuplicate(shell_ctx->fw[0], &work1));
161: PetscCall(VecDuplicate(shell_ctx->fw[0], &work2));
162: PetscCall(KSPCreate(PETSC_COMM_SELF, &check_ksp));
163: PetscCall(KSPSetNestLevel(check_ksp, pc->kspnestlevel));
164: PetscCall(MatIsSPDKnown(local_mat, &isset, &isspd));
165: if (isset && isspd) PetscCall(KSPSetType(check_ksp, KSPCG));
166: PetscCall(PetscObjectIncrementTabLevel((PetscObject)check_ksp, (PetscObject)local_ksp, 0));
167: PetscCall(KSPGetOptionsPrefix(local_ksp, &prefix));
168: PetscCall(KSPSetOptionsPrefix(check_ksp, prefix));
169: PetscCall(KSPAppendOptionsPrefix(check_ksp, "approximate_scale_"));
170: PetscCall(KSPSetErrorIfNotConverged(check_ksp, PETSC_FALSE));
171: PetscCall(KSPSetOperators(check_ksp, local_mat, local_pmat));
172: PetscCall(KSPSetComputeSingularValues(check_ksp, PETSC_TRUE));
173: PetscCall(KSPSetPreSolve(check_ksp, PCBDDCNullSpaceCorrPreSolve, shell_ctx));
174: PetscCall(KSPSetPostSolve(check_ksp, PCBDDCNullSpaceCorrPostSolve, shell_ctx));
175: PetscCall(KSPSetTolerances(check_ksp, PETSC_SMALL, PETSC_SMALL, PETSC_DEFAULT, PETSC_DEFAULT));
176: PetscCall(KSPSetFromOptions(check_ksp));
177: /* setup with default maxit, then set maxit to min(10,any_set_from_command_line) (bug in computing eigenvalues when changing the number of iterations */
178: PetscCall(KSPSetUp(check_ksp));
179: PetscCall(KSPGetPC(local_ksp, &local_pc));
180: PetscCall(KSPSetPC(check_ksp, local_pc));
181: PetscCall(KSPGetTolerances(check_ksp, NULL, NULL, NULL, &maxit));
182: PetscCall(KSPSetTolerances(check_ksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, PetscMin(10, maxit)));
183: PetscCall(VecSetRandom(work2, NULL));
184: PetscCall(MatMult(local_mat, work2, work1));
185: PetscCall(KSPSolve(check_ksp, work1, work1));
186: PetscCall(KSPCheckSolve(check_ksp, pc, work1));
187: PetscCall(VecAXPY(work1, -1., work2));
188: PetscCall(VecNorm(work1, NORM_INFINITY, &test_err));
189: PetscCall(KSPComputeExtremeSingularValues(check_ksp, &lambda_max, &lambda_min));
190: PetscCall(KSPGetIterationNumber(check_ksp, &k));
191: if (pcbddc->dbg_flag) {
192: if (isdir) {
193: PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d infinity error for Dirichlet adapted solver (no scale) %1.14e (it %" PetscInt_FMT ", eigs %1.6e %1.6e)\n", PetscGlobalRank, (double)test_err, k, (double)lambda_min, (double)lambda_max));
194: } else {
195: PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d infinity error for Neumann adapted solver (no scale) %1.14e (it %" PetscInt_FMT ", eigs %1.6e %1.6e)\n", PetscGlobalRank, (double)test_err, k, (double)lambda_min, (double)lambda_max));
196: }
197: }
198: if (needscaling) shell_ctx->scale = 1.0 / lambda_max;
200: if (needscaling && pcbddc->dbg_flag) { /* test for scaling factor */
201: PC new_pc;
203: PetscCall(VecSetRandom(work2, NULL));
204: PetscCall(MatMult(local_mat, work2, work1));
205: PetscCall(PCCreate(PetscObjectComm((PetscObject)check_ksp), &new_pc));
206: PetscCall(PCSetType(new_pc, PCKSP));
207: PetscCall(PCSetOperators(new_pc, local_mat, local_pmat));
208: PetscCall(PCKSPSetKSP(new_pc, local_ksp));
209: PetscCall(KSPSetPC(check_ksp, new_pc));
210: PetscCall(PCDestroy(&new_pc));
211: PetscCall(KSPSetTolerances(check_ksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, maxit));
212: PetscCall(KSPSetPreSolve(check_ksp, NULL, NULL));
213: PetscCall(KSPSetPostSolve(check_ksp, NULL, NULL));
214: PetscCall(KSPSolve(check_ksp, work1, work1));
215: PetscCall(KSPCheckSolve(check_ksp, pc, work1));
216: PetscCall(VecAXPY(work1, -1., work2));
217: PetscCall(VecNorm(work1, NORM_INFINITY, &test_err));
218: PetscCall(KSPComputeExtremeSingularValues(check_ksp, &lambda_max, &lambda_min));
219: PetscCall(KSPGetIterationNumber(check_ksp, &k));
220: if (pcbddc->dbg_flag) {
221: if (isdir) {
222: PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d infinity error for Dirichlet adapted solver (scale %g) %1.14e (it %" PetscInt_FMT ", eigs %1.6e %1.6e)\n", PetscGlobalRank, (double)PetscRealPart(shell_ctx->scale), (double)test_err, k, (double)lambda_min, (double)lambda_max));
223: } else {
224: PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d infinity error for Neumann adapted solver (scale %g) %1.14e (it %" PetscInt_FMT ", eigs %1.6e %1.6e)\n", PetscGlobalRank, (double)PetscRealPart(shell_ctx->scale), (double)test_err, k, (double)lambda_min, (double)lambda_max));
225: }
226: }
227: }
228: PetscCall(KSPDestroy(&check_ksp));
229: PetscCall(VecDestroy(&work1));
230: PetscCall(VecDestroy(&work2));
231: }
232: PetscCall(PetscLogEventEnd(PC_BDDC_ApproxSetUp[pcbddc->current_level], pc, 0, 0, 0));
234: if (pcbddc->dbg_flag) {
235: Vec work1, work2, work3;
236: PetscReal test_err;
238: /* check nullspace basis is solved exactly */
239: PetscCall(VecDuplicate(shell_ctx->fw[0], &work1));
240: PetscCall(VecDuplicate(shell_ctx->fw[0], &work2));
241: PetscCall(VecDuplicate(shell_ctx->fw[0], &work3));
242: PetscCall(VecSetRandom(shell_ctx->sw[0], NULL));
243: PetscCall(MatMult(shell_ctx->basis_mat, shell_ctx->sw[0], work1));
244: PetscCall(VecCopy(work1, work2));
245: PetscCall(MatMult(local_mat, work1, work3));
246: PetscCall(KSPSolve(local_ksp, work3, work1));
247: PetscCall(VecAXPY(work1, -1., work2));
248: PetscCall(VecNorm(work1, NORM_INFINITY, &test_err));
249: if (isdir) {
250: PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d infinity error for Dirichlet nullspace correction solver: %1.14e\n", PetscGlobalRank, (double)test_err));
251: } else {
252: PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d infinity error for Neumann nullspace correction solver: %1.14e\n", PetscGlobalRank, (double)test_err));
253: }
254: PetscCall(VecDestroy(&work1));
255: PetscCall(VecDestroy(&work2));
256: PetscCall(VecDestroy(&work3));
257: }
258: PetscFunctionReturn(PETSC_SUCCESS);
259: }