Actual source code: bddcnullspace.c

  1: #include <../src/ksp/pc/impls/bddc/bddc.h>
  2: #include <../src/ksp/pc/impls/bddc/bddcprivate.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) {
 82:       PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d doesn't have local (near) nullspace: no need for correction in %s solver \n",PetscGlobalRank,isdir ? "Dirichlet" : "Neumann");
 83:     }
 84:     return 0;
 85:   }
 86:   PetscObjectQuery((PetscObject)NullSpace,"_PBDDC_Null_dmat",(PetscObject*)&dmat);
 88:   PetscLogEventBegin(PC_BDDC_ApproxSetUp[pcbddc->current_level],pc,0,0,0);

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

 97:   MatGetOption(local_mat,MAT_SYMMETRIC,&shell_ctx->symm);

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

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

128:     ISGetLocalSize(zerorows,&nz);
129:     ISGetIndices(zerorows,&idxs);
130:     for (i=0;i<nz;i++) {
131:       MatSetValue(shell_ctx->inv_smat,idxs[i],idxs[i],0.0,INSERT_VALUES);
132:     }
133:     ISRestoreIndices(zerorows,&idxs);
134:     MatAssemblyBegin(shell_ctx->inv_smat,MAT_FINAL_ASSEMBLY);
135:     MatAssemblyEnd(shell_ctx->inv_smat,MAT_FINAL_ASSEMBLY);
136:   }
137:   ISDestroy(&zerorows);

139:   /* Create work vectors in shell context */
140:   MatCreateVecs(shell_ctx->inv_smat,&v,NULL);
141:   KSPCreateVecs(local_ksp,1,&shell_ctx->fw,0,NULL);
142:   VecDuplicateVecs(v,3,&shell_ctx->sw);
143:   VecDestroy(&v);

145:   /* add special pre/post solve to KSP (see [1], eq. 48) */
146:   KSPSetPreSolve(local_ksp,PCBDDCNullSpaceCorrPreSolve,shell_ctx);
147:   KSPSetPostSolve(local_ksp,PCBDDCNullSpaceCorrPostSolve,shell_ctx);
148:   PetscContainerCreate(PetscObjectComm((PetscObject)local_ksp),&c);
149:   PetscContainerSetPointer(c,shell_ctx);
150:   PetscContainerSetUserDestroy(c,PCBDDCNullSpaceCorrDestroy);
151:   PetscObjectCompose((PetscObject)local_ksp,"_PCBDDC_Null_PrePost_ctx",(PetscObject)c);
152:   PetscContainerDestroy(&c);

154:   /* Create ksp object suitable for extreme eigenvalues' estimation */
155:   if (needscaling || pcbddc->dbg_flag) {
156:     KSP         check_ksp;
157:     PC          local_pc;
158:     Vec         work1,work2;
159:     const char* prefix;
160:     PetscReal   test_err,lambda_min,lambda_max;
161:     PetscInt    k,maxit;

163:     VecDuplicate(shell_ctx->fw[0],&work1);
164:     VecDuplicate(shell_ctx->fw[0],&work2);
165:     KSPCreate(PETSC_COMM_SELF,&check_ksp);
166:     if (local_mat->spd) {
167:       KSPSetType(check_ksp,KSPCG);
168:     }
169:     PetscObjectIncrementTabLevel((PetscObject)check_ksp,(PetscObject)local_ksp,0);
170:     KSPGetOptionsPrefix(local_ksp,&prefix);
171:     KSPSetOptionsPrefix(check_ksp,prefix);
172:     KSPAppendOptionsPrefix(check_ksp,"approximate_scale_");
173:     KSPSetErrorIfNotConverged(check_ksp,PETSC_FALSE);
174:     KSPSetOperators(check_ksp,local_mat,local_pmat);
175:     KSPSetComputeSingularValues(check_ksp,PETSC_TRUE);
176:     KSPSetPreSolve(check_ksp,PCBDDCNullSpaceCorrPreSolve,shell_ctx);
177:     KSPSetPostSolve(check_ksp,PCBDDCNullSpaceCorrPostSolve,shell_ctx);
178:     KSPSetTolerances(check_ksp,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,PETSC_DEFAULT);
179:     KSPSetFromOptions(check_ksp);
180:     /* setup with default maxit, then set maxit to min(10,any_set_from_command_line) (bug in computing eigenvalues when chaning the number of iterations */
181:     KSPSetUp(check_ksp);
182:     KSPGetPC(local_ksp,&local_pc);
183:     KSPSetPC(check_ksp,local_pc);
184:     KSPGetTolerances(check_ksp,NULL,NULL,NULL,&maxit);
185:     KSPSetTolerances(check_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,PetscMin(10,maxit));
186:     VecSetRandom(work2,NULL);
187:     MatMult(local_mat,work2,work1);
188:     KSPSolve(check_ksp,work1,work1);
189:     KSPCheckSolve(check_ksp,pc,work1);
190:     VecAXPY(work1,-1.,work2);
191:     VecNorm(work1,NORM_INFINITY,&test_err);
192:     KSPComputeExtremeSingularValues(check_ksp,&lambda_max,&lambda_min);
193:     KSPGetIterationNumber(check_ksp,&k);
194:     if (pcbddc->dbg_flag) {
195:       if (isdir) {
196:         PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Dirichlet adapted solver (no scale) %1.14e (it %D, eigs %1.6e %1.6e)\n",PetscGlobalRank,test_err,k,lambda_min,lambda_max);
197:       } else {
198:         PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Neumann adapted solver (no scale) %1.14e (it %D, eigs %1.6e %1.6e)\n",PetscGlobalRank,test_err,k,lambda_min,lambda_max);
199:       }
200:     }
201:     if (needscaling) shell_ctx->scale = 1.0/lambda_max;

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

206:       VecSetRandom(work2,NULL);
207:       MatMult(local_mat,work2,work1);
208:       PCCreate(PetscObjectComm((PetscObject)check_ksp),&new_pc);
209:       PCSetType(new_pc,PCKSP);
210:       PCSetOperators(new_pc,local_mat,local_pmat);
211:       PCKSPSetKSP(new_pc,local_ksp);
212:       KSPSetPC(check_ksp,new_pc);
213:       PCDestroy(&new_pc);
214:       KSPSetTolerances(check_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,maxit);
215:       KSPSetPreSolve(check_ksp,NULL,NULL);
216:       KSPSetPostSolve(check_ksp,NULL,NULL);
217:       KSPSolve(check_ksp,work1,work1);
218:       KSPCheckSolve(check_ksp,pc,work1);
219:       VecAXPY(work1,-1.,work2);
220:       VecNorm(work1,NORM_INFINITY,&test_err);
221:       KSPComputeExtremeSingularValues(check_ksp,&lambda_max,&lambda_min);
222:       KSPGetIterationNumber(check_ksp,&k);
223:       if (pcbddc->dbg_flag) {
224:         if (isdir) {
225:           PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Dirichlet adapted solver (scale %g) %1.14e (it %D, eigs %1.6e %1.6e)\n",PetscGlobalRank,(double)PetscRealPart(shell_ctx->scale),test_err,k,lambda_min,lambda_max);
226:         } else {
227:           PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Neumann adapted solver (scale %g) %1.14e (it %D, eigs %1.6e %1.6e)\n",PetscGlobalRank,(double)PetscRealPart(shell_ctx->scale),test_err,k,lambda_min,lambda_max);
228:         }
229:       }
230:     }
231:     KSPDestroy(&check_ksp);
232:     VecDestroy(&work1);
233:     VecDestroy(&work2);
234:   }
235:   PetscLogEventEnd(PC_BDDC_ApproxSetUp[pcbddc->current_level],pc,0,0,0);

237:   if (pcbddc->dbg_flag) {
238:     Vec       work1,work2,work3;
239:     PetscReal test_err;

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