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:   PetscLogEventBegin(corr_ctx->evapply, ksp, 0, 0, 0);
 12:   MatMultTranspose(corr_ctx->basis_mat, y, corr_ctx->sw[0]);
 13:   if (corr_ctx->symm) {
 14:     MatMult(corr_ctx->inv_smat, corr_ctx->sw[0], corr_ctx->sw[1]);
 15:   } else {
 16:     MatMultTranspose(corr_ctx->inv_smat, corr_ctx->sw[0], corr_ctx->sw[1]);
 17:   }
 18:   VecScale(corr_ctx->sw[1], -1.0);
 19:   MatMult(corr_ctx->basis_mat, corr_ctx->sw[1], corr_ctx->fw[0]);
 20:   VecScale(corr_ctx->sw[1], -1.0);
 21:   KSPGetOperators(ksp, &K, NULL);
 22:   MatMultAdd(K, corr_ctx->fw[0], y, y);
 23:   PetscLogEventEnd(corr_ctx->evapply, ksp, 0, 0, 0);
 24:   return 0;
 25: }

 27: /* E^t + small */
 28: static PetscErrorCode PCBDDCNullSpaceCorrPostSolve(KSP ksp, Vec y, Vec x, void *ctx)
 29: {
 30:   NullSpaceCorrection_ctx corr_ctx = (NullSpaceCorrection_ctx)ctx;
 31:   Mat                     K;

 33:   PetscLogEventBegin(corr_ctx->evapply, ksp, 0, 0, 0);
 34:   KSPGetOperators(ksp, &K, NULL);
 35:   if (corr_ctx->symm) {
 36:     MatMult(K, x, corr_ctx->fw[0]);
 37:   } else {
 38:     MatMultTranspose(K, x, corr_ctx->fw[0]);
 39:   }
 40:   MatMultTranspose(corr_ctx->basis_mat, corr_ctx->fw[0], corr_ctx->sw[0]);
 41:   VecScale(corr_ctx->sw[0], -1.0);
 42:   MatMult(corr_ctx->inv_smat, corr_ctx->sw[0], corr_ctx->sw[2]);
 43:   MatMultAdd(corr_ctx->basis_mat, corr_ctx->sw[2], x, corr_ctx->fw[0]);
 44:   VecScale(corr_ctx->fw[0], corr_ctx->scale);
 45:   /* Sum contributions from approximate solver and projected system */
 46:   MatMultAdd(corr_ctx->basis_mat, corr_ctx->sw[1], corr_ctx->fw[0], x);
 47:   PetscLogEventEnd(corr_ctx->evapply, ksp, 0, 0, 0);
 48:   return 0;
 49: }

 51: static PetscErrorCode PCBDDCNullSpaceCorrDestroy(void *ctx)
 52: {
 53:   NullSpaceCorrection_ctx corr_ctx = (NullSpaceCorrection_ctx)ctx;

 55:   VecDestroyVecs(3, &corr_ctx->sw);
 56:   VecDestroyVecs(1, &corr_ctx->fw);
 57:   MatDestroy(&corr_ctx->basis_mat);
 58:   MatDestroy(&corr_ctx->inv_smat);
 59:   PetscFree(corr_ctx);
 60:   return 0;
 61: }

 63: PetscErrorCode PCBDDCNullSpaceAssembleCorrection(PC pc, PetscBool isdir, PetscBool needscaling)
 64: {
 65:   PC_BDDC                *pcbddc    = (PC_BDDC *)pc->data;
 66:   MatNullSpace            NullSpace = NULL;
 67:   KSP                     local_ksp;
 68:   NullSpaceCorrection_ctx shell_ctx;
 69:   Mat                     local_mat, local_pmat, dmat, Kbasis_mat;
 70:   Vec                     v;
 71:   PetscContainer          c;
 72:   PetscInt                basis_size;
 73:   IS                      zerorows;
 74:   PetscBool               iscusp;

 76:   if (isdir) local_ksp = pcbddc->ksp_D; /* Dirichlet solver */
 77:   else local_ksp = pcbddc->ksp_R;       /* Neumann solver */
 78:   KSPGetOperators(local_ksp, &local_mat, &local_pmat);
 79:   MatGetNearNullSpace(local_pmat, &NullSpace);
 80:   if (!NullSpace) {
 81:     if (pcbddc->dbg_flag) PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d doesn't have local (near) nullspace: no need for correction in %s solver \n", PetscGlobalRank, isdir ? "Dirichlet" : "Neumann");
 82:     return 0;
 83:   }
 84:   PetscObjectQuery((PetscObject)NullSpace, "_PBDDC_Null_dmat", (PetscObject *)&dmat);
 86:   PetscLogEventBegin(PC_BDDC_ApproxSetUp[pcbddc->current_level], pc, 0, 0, 0);

 88:   PetscNew(&shell_ctx);
 89:   shell_ctx->scale = 1.0;
 90:   PetscObjectReference((PetscObject)dmat);
 91:   shell_ctx->basis_mat = dmat;
 92:   MatGetSize(dmat, NULL, &basis_size);
 93:   shell_ctx->evapply = PC_BDDC_ApproxApply[pcbddc->current_level];

 95:   MatIsSymmetric(local_mat, 0.0, &shell_ctx->symm);

 97:   /* explicit construct (Phi^T K Phi)^-1 */
 98:   PetscObjectTypeCompare((PetscObject)local_mat, MATSEQAIJCUSPARSE, &iscusp);
 99:   if (iscusp) MatConvert(shell_ctx->basis_mat, MATSEQDENSECUDA, MAT_INPLACE_MATRIX, &shell_ctx->basis_mat);
100:   MatMatMult(local_mat, shell_ctx->basis_mat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Kbasis_mat);
101:   MatTransposeMatMult(Kbasis_mat, shell_ctx->basis_mat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &shell_ctx->inv_smat);
102:   MatDestroy(&Kbasis_mat);
103:   MatBindToCPU(shell_ctx->inv_smat, PETSC_TRUE);
104:   MatFindZeroRows(shell_ctx->inv_smat, &zerorows);
105:   if (zerorows) { /* linearly dependent basis */
106:     const PetscInt *idxs;
107:     PetscInt        i, nz;

109:     ISGetLocalSize(zerorows, &nz);
110:     ISGetIndices(zerorows, &idxs);
111:     for (i = 0; i < nz; i++) MatSetValue(shell_ctx->inv_smat, idxs[i], idxs[i], 1.0, INSERT_VALUES);
112:     ISRestoreIndices(zerorows, &idxs);
113:     MatAssemblyBegin(shell_ctx->inv_smat, MAT_FINAL_ASSEMBLY);
114:     MatAssemblyEnd(shell_ctx->inv_smat, MAT_FINAL_ASSEMBLY);
115:   }
116:   MatLUFactor(shell_ctx->inv_smat, NULL, NULL, NULL);
117:   MatSeqDenseInvertFactors_Private(shell_ctx->inv_smat);
118:   if (zerorows) { /* linearly dependent basis */
119:     const PetscInt *idxs;
120:     PetscInt        i, nz;

122:     ISGetLocalSize(zerorows, &nz);
123:     ISGetIndices(zerorows, &idxs);
124:     for (i = 0; i < nz; i++) MatSetValue(shell_ctx->inv_smat, idxs[i], idxs[i], 0.0, INSERT_VALUES);
125:     ISRestoreIndices(zerorows, &idxs);
126:     MatAssemblyBegin(shell_ctx->inv_smat, MAT_FINAL_ASSEMBLY);
127:     MatAssemblyEnd(shell_ctx->inv_smat, MAT_FINAL_ASSEMBLY);
128:   }
129:   ISDestroy(&zerorows);

131:   /* Create work vectors in shell context */
132:   MatCreateVecs(shell_ctx->inv_smat, &v, NULL);
133:   KSPCreateVecs(local_ksp, 1, &shell_ctx->fw, 0, NULL);
134:   VecDuplicateVecs(v, 3, &shell_ctx->sw);
135:   VecDestroy(&v);

137:   /* add special pre/post solve to KSP (see [1], eq. 48) */
138:   KSPSetPreSolve(local_ksp, PCBDDCNullSpaceCorrPreSolve, shell_ctx);
139:   KSPSetPostSolve(local_ksp, PCBDDCNullSpaceCorrPostSolve, shell_ctx);
140:   PetscContainerCreate(PetscObjectComm((PetscObject)local_ksp), &c);
141:   PetscContainerSetPointer(c, shell_ctx);
142:   PetscContainerSetUserDestroy(c, PCBDDCNullSpaceCorrDestroy);
143:   PetscObjectCompose((PetscObject)local_ksp, "_PCBDDC_Null_PrePost_ctx", (PetscObject)c);
144:   PetscContainerDestroy(&c);

146:   /* Create ksp object suitable for extreme eigenvalues' estimation */
147:   if (needscaling || pcbddc->dbg_flag) {
148:     KSP         check_ksp;
149:     PC          local_pc;
150:     Vec         work1, work2;
151:     const char *prefix;
152:     PetscReal   test_err, lambda_min, lambda_max;
153:     PetscInt    k, maxit;
154:     PetscBool   isspd, isset;

156:     VecDuplicate(shell_ctx->fw[0], &work1);
157:     VecDuplicate(shell_ctx->fw[0], &work2);
158:     KSPCreate(PETSC_COMM_SELF, &check_ksp);
159:     MatIsSPDKnown(local_mat, &isset, &isspd);
160:     if (isset && isspd) KSPSetType(check_ksp, KSPCG);
161:     PetscObjectIncrementTabLevel((PetscObject)check_ksp, (PetscObject)local_ksp, 0);
162:     KSPGetOptionsPrefix(local_ksp, &prefix);
163:     KSPSetOptionsPrefix(check_ksp, prefix);
164:     KSPAppendOptionsPrefix(check_ksp, "approximate_scale_");
165:     KSPSetErrorIfNotConverged(check_ksp, PETSC_FALSE);
166:     KSPSetOperators(check_ksp, local_mat, local_pmat);
167:     KSPSetComputeSingularValues(check_ksp, PETSC_TRUE);
168:     KSPSetPreSolve(check_ksp, PCBDDCNullSpaceCorrPreSolve, shell_ctx);
169:     KSPSetPostSolve(check_ksp, PCBDDCNullSpaceCorrPostSolve, shell_ctx);
170:     KSPSetTolerances(check_ksp, PETSC_SMALL, PETSC_SMALL, PETSC_DEFAULT, PETSC_DEFAULT);
171:     KSPSetFromOptions(check_ksp);
172:     /* 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 */
173:     KSPSetUp(check_ksp);
174:     KSPGetPC(local_ksp, &local_pc);
175:     KSPSetPC(check_ksp, local_pc);
176:     KSPGetTolerances(check_ksp, NULL, NULL, NULL, &maxit);
177:     KSPSetTolerances(check_ksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, PetscMin(10, maxit));
178:     VecSetRandom(work2, NULL);
179:     MatMult(local_mat, work2, work1);
180:     KSPSolve(check_ksp, work1, work1);
181:     KSPCheckSolve(check_ksp, pc, work1);
182:     VecAXPY(work1, -1., work2);
183:     VecNorm(work1, NORM_INFINITY, &test_err);
184:     KSPComputeExtremeSingularValues(check_ksp, &lambda_max, &lambda_min);
185:     KSPGetIterationNumber(check_ksp, &k);
186:     if (pcbddc->dbg_flag) {
187:       if (isdir) {
188:         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);
189:       } else {
190:         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);
191:       }
192:     }
193:     if (needscaling) shell_ctx->scale = 1.0 / lambda_max;

195:     if (needscaling && pcbddc->dbg_flag) { /* test for scaling factor */
196:       PC new_pc;

198:       VecSetRandom(work2, NULL);
199:       MatMult(local_mat, work2, work1);
200:       PCCreate(PetscObjectComm((PetscObject)check_ksp), &new_pc);
201:       PCSetType(new_pc, PCKSP);
202:       PCSetOperators(new_pc, local_mat, local_pmat);
203:       PCKSPSetKSP(new_pc, local_ksp);
204:       KSPSetPC(check_ksp, new_pc);
205:       PCDestroy(&new_pc);
206:       KSPSetTolerances(check_ksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, maxit);
207:       KSPSetPreSolve(check_ksp, NULL, NULL);
208:       KSPSetPostSolve(check_ksp, NULL, NULL);
209:       KSPSolve(check_ksp, work1, work1);
210:       KSPCheckSolve(check_ksp, pc, work1);
211:       VecAXPY(work1, -1., work2);
212:       VecNorm(work1, NORM_INFINITY, &test_err);
213:       KSPComputeExtremeSingularValues(check_ksp, &lambda_max, &lambda_min);
214:       KSPGetIterationNumber(check_ksp, &k);
215:       if (pcbddc->dbg_flag) {
216:         if (isdir) {
217:           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);
218:         } else {
219:           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);
220:         }
221:       }
222:     }
223:     KSPDestroy(&check_ksp);
224:     VecDestroy(&work1);
225:     VecDestroy(&work2);
226:   }
227:   PetscLogEventEnd(PC_BDDC_ApproxSetUp[pcbddc->current_level], pc, 0, 0, 0);

229:   if (pcbddc->dbg_flag) {
230:     Vec       work1, work2, work3;
231:     PetscReal test_err;

233:     /* check nullspace basis is solved exactly */
234:     VecDuplicate(shell_ctx->fw[0], &work1);
235:     VecDuplicate(shell_ctx->fw[0], &work2);
236:     VecDuplicate(shell_ctx->fw[0], &work3);
237:     VecSetRandom(shell_ctx->sw[0], NULL);
238:     MatMult(shell_ctx->basis_mat, shell_ctx->sw[0], work1);
239:     VecCopy(work1, work2);
240:     MatMult(local_mat, work1, work3);
241:     KSPSolve(local_ksp, work3, work1);
242:     VecAXPY(work1, -1., work2);
243:     VecNorm(work1, NORM_INFINITY, &test_err);
244:     if (isdir) {
245:       PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d infinity error for Dirichlet nullspace correction solver: %1.14e\n", PetscGlobalRank, (double)test_err);
246:     } else {
247:       PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d infinity error for Neumann nullspace correction solver: %1.14e\n", PetscGlobalRank, (double)test_err);
248:     }
249:     VecDestroy(&work1);
250:     VecDestroy(&work2);
251:     VecDestroy(&work3);
252:   }
253:   return 0;
254: }