Actual source code: bddcnullspace.c
petsc-3.10.5 2019-03-28
1: #include <../src/ksp/pc/impls/bddc/bddc.h>
2: #include <../src/ksp/pc/impls/bddc/bddcprivate.h>
4: static PetscErrorCode PCBDDCViewNullSpaceCorrectionPC(PC pc,PetscViewer view)
5: {
6: PetscErrorCode ierr;
7: PetscBool isascii;
8: NullSpaceCorrection_ctx pc_ctx;
11: PCShellGetContext(pc,(void**)&pc_ctx);
12: PetscObjectTypeCompare((PetscObject)view,PETSCVIEWERASCII,&isascii);
13: if (isascii) {
14: PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_INFO);
16: PetscViewerASCIIPrintf(view,"L:\n");
17: PetscViewerASCIIPushTab(view);
18: MatView(pc_ctx->Lbasis_mat,view);
19: PetscViewerASCIIPopTab(view);
21: PetscViewerASCIIPrintf(view,"K:\n");
22: PetscViewerASCIIPushTab(view);
23: MatView(pc_ctx->Kbasis_mat,view);
24: PetscViewerASCIIPopTab(view);
26: PetscViewerPopFormat(view);
28: PetscViewerASCIIPrintf(view,"inner preconditioner:\n");
29: PetscViewerASCIIPushTab(view);
30: PCView(pc_ctx->local_pc,view);
31: PetscViewerASCIIPopTab(view);
32: }
33: return(0);
34: }
36: static PetscErrorCode PCBDDCApplyNullSpaceCorrectionPC(PC pc,Vec x,Vec y)
37: {
38: NullSpaceCorrection_ctx pc_ctx;
39: PetscErrorCode ierr;
42: PCShellGetContext(pc,(void**)&pc_ctx);
43: /* E */
44: MatMultTranspose(pc_ctx->Lbasis_mat,x,pc_ctx->work_small_2);
45: MatMultAdd(pc_ctx->Kbasis_mat,pc_ctx->work_small_2,x,pc_ctx->work_full_1);
46: /* P^-1 */
47: PCApply(pc_ctx->local_pc,pc_ctx->work_full_1,pc_ctx->work_full_2);
48: /* E^T */
49: MatMultTranspose(pc_ctx->Kbasis_mat,pc_ctx->work_full_2,pc_ctx->work_small_1);
50: VecScale(pc_ctx->work_small_1,-1.0);
51: MatMultAdd(pc_ctx->Lbasis_mat,pc_ctx->work_small_1,pc_ctx->work_full_2,pc_ctx->work_full_1);
52: /* Sum contributions */
53: MatMultAdd(pc_ctx->basis_mat,pc_ctx->work_small_2,pc_ctx->work_full_1,y);
54: if (pc_ctx->apply_scaling) {
55: VecScale(y,pc_ctx->scale);
56: }
57: return(0);
58: }
60: static PetscErrorCode PCBDDCDestroyNullSpaceCorrectionPC(PC pc)
61: {
62: NullSpaceCorrection_ctx pc_ctx;
63: PetscErrorCode ierr;
66: PCShellGetContext(pc,(void**)&pc_ctx);
67: VecDestroy(&pc_ctx->work_small_1);
68: VecDestroy(&pc_ctx->work_small_2);
69: VecDestroy(&pc_ctx->work_full_1);
70: VecDestroy(&pc_ctx->work_full_2);
71: MatDestroy(&pc_ctx->basis_mat);
72: MatDestroy(&pc_ctx->Lbasis_mat);
73: MatDestroy(&pc_ctx->Kbasis_mat);
74: PCDestroy(&pc_ctx->local_pc);
75: PetscFree(pc_ctx);
76: return(0);
77: }
79: PetscErrorCode PCBDDCNullSpaceAssembleCorrection(PC pc, PetscBool isdir, PetscBool needscaling)
80: {
81: PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
82: PC_IS *pcis = (PC_IS*)pc->data;
83: Mat_IS *matis = (Mat_IS*)pc->pmat->data;
84: MatNullSpace NullSpace = NULL;
85: KSP local_ksp;
86: PC newpc;
87: NullSpaceCorrection_ctx shell_ctx;
88: Mat local_mat,local_pmat,small_mat,inv_small_mat;
89: Vec *nullvecs;
90: VecScatter scatter_ctx;
91: IS is_aux,local_dofs;
92: MatFactorInfo matinfo;
93: PetscScalar *basis_mat,*Kbasis_mat,*array,*array_mat;
94: PetscScalar one = 1.0,zero = 0.0, m_one = -1.0;
95: PetscInt basis_dofs,basis_size,nnsp_size,i,k;
96: PetscBool nnsp_has_cnst;
97: PetscReal test_err,lambda_min,lambda_max;
98: PetscErrorCode ierr;
101: MatGetNullSpace(matis->A,&NullSpace);
102: if (!NullSpace) {
103: MatGetNearNullSpace(matis->A,&NullSpace);
104: }
105: if (!NullSpace) {
106: if (pcbddc->dbg_flag) {
107: PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d doesn't have local (near) nullspace: no need for correction in %s solver \n",PetscGlobalRank,isdir ? "Dirichlet" : "Neumann");
108: }
109: return(0);
110: }
112: /* Infer the local solver */
113: if (isdir) {
114: /* Dirichlet solver */
115: local_ksp = pcbddc->ksp_D;
116: local_dofs = pcis->is_I_local;
117: } else {
118: /* Neumann solver */
119: local_ksp = pcbddc->ksp_R;
120: local_dofs = pcbddc->is_R_local;
121: }
122: ISGetSize(local_dofs,&basis_dofs);
123: KSPGetOperators(local_ksp,&local_mat,&local_pmat);
125: /* Get null space vecs */
126: MatNullSpaceGetVecs(NullSpace,&nnsp_has_cnst,&nnsp_size,(const Vec**)&nullvecs);
127: basis_size = nnsp_size;
128: if (nnsp_has_cnst) basis_size++;
130: /* Create shell ctx */
131: PetscNew(&shell_ctx);
132: shell_ctx->apply_scaling = needscaling;
133: shell_ctx->scale = 1.0;
135: /* Create work vectors in shell context */
136: VecCreate(PETSC_COMM_SELF,&shell_ctx->work_small_1);
137: VecSetSizes(shell_ctx->work_small_1,basis_size,basis_size);
138: VecSetType(shell_ctx->work_small_1,((PetscObject)pcis->vec1_B)->type_name);
139: VecDuplicate(shell_ctx->work_small_1,&shell_ctx->work_small_2);
140: VecCreate(PETSC_COMM_SELF,&shell_ctx->work_full_1);
141: VecSetSizes(shell_ctx->work_full_1,basis_dofs,basis_dofs);
142: VecSetType(shell_ctx->work_full_1,((PetscObject)pcis->vec1_B)->type_name);
143: VecDuplicate(shell_ctx->work_full_1,&shell_ctx->work_full_2);
145: /* Allocate workspace */
146: MatCreateSeqDense(PETSC_COMM_SELF,basis_dofs,basis_size,NULL,&shell_ctx->basis_mat);
147: MatCreateSeqDense(PETSC_COMM_SELF,basis_dofs,basis_size,NULL,&shell_ctx->Kbasis_mat);
148: MatDenseGetArray(shell_ctx->basis_mat,&basis_mat);
149: MatDenseGetArray(shell_ctx->Kbasis_mat,&Kbasis_mat);
151: /* Restrict local null space on selected dofs
152: and compute matrices N and K*N */
153: VecScatterCreate(pcis->vec1_N,local_dofs,shell_ctx->work_full_1,(IS)0,&scatter_ctx);
154: for (k=0;k<nnsp_size;k++) {
155: VecPlaceArray(shell_ctx->work_full_1,(const PetscScalar*)&basis_mat[k*basis_dofs]);
156: VecScatterBegin(scatter_ctx,nullvecs[k],shell_ctx->work_full_1,INSERT_VALUES,SCATTER_FORWARD);
157: VecScatterEnd(scatter_ctx,nullvecs[k],shell_ctx->work_full_1,INSERT_VALUES,SCATTER_FORWARD);
158: VecResetArray(shell_ctx->work_full_1);
159: }
160: if (nnsp_has_cnst) {
161: VecPlaceArray(shell_ctx->work_full_1,(const PetscScalar*)&basis_mat[k*basis_dofs]);
162: VecSet(shell_ctx->work_full_1,one);
163: VecResetArray(shell_ctx->work_full_1);
164: }
166: PetscMalloc1(basis_size,&nullvecs);
167: for (k=0;k<basis_size;k++) {
168: VecCreateSeqWithArray(PETSC_COMM_SELF,1,basis_dofs,basis_mat + k*basis_dofs,&nullvecs[k]);
169: }
170: PCBDDCOrthonormalizeVecs(basis_size,nullvecs);
171: MatNullSpaceCreate(PETSC_COMM_SELF,PETSC_FALSE,basis_size,nullvecs,&NullSpace);
172: MatSetNearNullSpace(local_mat,NullSpace);
173: MatNullSpaceDestroy(&NullSpace);
174: for (k=0;k<basis_size;k++) {
175: VecDestroy(&nullvecs[k]);
176: }
177: PetscFree(nullvecs);
179: for (k=0;k<basis_size;k++) {
180: VecPlaceArray(shell_ctx->work_full_1,(const PetscScalar*)&basis_mat[k*basis_dofs]);
181: VecPlaceArray(shell_ctx->work_full_2,(const PetscScalar*)&Kbasis_mat[k*basis_dofs]);
182: MatMult(local_mat,shell_ctx->work_full_1,shell_ctx->work_full_2);
183: VecResetArray(shell_ctx->work_full_1);
184: VecResetArray(shell_ctx->work_full_2);
185: }
187: VecScatterDestroy(&scatter_ctx);
188: MatDenseRestoreArray(shell_ctx->basis_mat,&basis_mat);
189: MatDenseRestoreArray(shell_ctx->Kbasis_mat,&Kbasis_mat);
190: /* Assemble another Mat object in shell context */
191: MatTransposeMatMult(shell_ctx->basis_mat,shell_ctx->Kbasis_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&small_mat);
192: MatFactorInfoInitialize(&matinfo);
193: ISCreateStride(PETSC_COMM_SELF,basis_size,0,1,&is_aux);
194: MatLUFactor(small_mat,is_aux,is_aux,&matinfo);
195: ISDestroy(&is_aux);
196: PetscMalloc1(basis_size*basis_size,&array_mat);
197: for (k=0;k<basis_size;k++) {
198: VecSet(shell_ctx->work_small_1,zero);
199: VecSetValue(shell_ctx->work_small_1,k,one,INSERT_VALUES);
200: VecAssemblyBegin(shell_ctx->work_small_1);
201: VecAssemblyEnd(shell_ctx->work_small_1);
202: MatSolve(small_mat,shell_ctx->work_small_1,shell_ctx->work_small_2);
203: VecGetArrayRead(shell_ctx->work_small_2,(const PetscScalar**)&array);
204: for (i=0;i<basis_size;i++) {
205: array_mat[i*basis_size+k]=array[i];
206: }
207: VecRestoreArrayRead(shell_ctx->work_small_2,(const PetscScalar**)&array);
208: }
209: MatCreateSeqDense(PETSC_COMM_SELF,basis_size,basis_size,array_mat,&inv_small_mat);
210: MatMatMult(shell_ctx->basis_mat,inv_small_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&shell_ctx->Lbasis_mat);
211: PetscFree(array_mat);
212: MatDestroy(&inv_small_mat);
213: MatDestroy(&small_mat);
214: MatScale(shell_ctx->Kbasis_mat,m_one);
216: /* Rebuild local PC */
217: KSPGetPC(local_ksp,&shell_ctx->local_pc);
218: PetscObjectReference((PetscObject)shell_ctx->local_pc);
219: PCCreate(PETSC_COMM_SELF,&newpc);
220: PCSetOperators(newpc,local_mat,local_mat);
221: PCSetType(newpc,PCSHELL);
222: PCShellSetContext(newpc,shell_ctx);
223: if (isdir) {
224: PCShellSetName(newpc,"Nullspace corrected interior solve");
225: } else {
226: PCShellSetName(newpc,"Nullspace corrected correction solve");
227: }
228: PCShellSetApply(newpc,PCBDDCApplyNullSpaceCorrectionPC);
229: PCShellSetView(newpc,PCBDDCViewNullSpaceCorrectionPC);
230: PCShellSetDestroy(newpc,PCBDDCDestroyNullSpaceCorrectionPC);
231: PCSetUp(newpc);
232: KSPSetPC(local_ksp,newpc);
233: PetscObjectIncrementTabLevel((PetscObject)newpc,(PetscObject)local_ksp,0);
234: KSPSetUp(local_ksp);
236: /* Create ksp object suitable for extreme eigenvalues' estimation */
237: if (needscaling) {
238: KSP check_ksp;
239: Vec *workv;
240: const char* prefix;
242: KSPCreate(PETSC_COMM_SELF,&check_ksp);
243: PetscObjectIncrementTabLevel((PetscObject)check_ksp,(PetscObject)local_ksp,0);
244: KSPGetOptionsPrefix(local_ksp,&prefix);
245: KSPSetOptionsPrefix(check_ksp,prefix);
246: KSPAppendOptionsPrefix(check_ksp,"approxscale_");
247: KSPSetErrorIfNotConverged(check_ksp,PETSC_FALSE);
248: KSPSetOperators(check_ksp,local_mat,local_mat);
249: KSPSetTolerances(check_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,10);
250: KSPSetComputeSingularValues(check_ksp,PETSC_TRUE);
251: VecDuplicateVecs(shell_ctx->work_full_1,2,&workv);
252: KSPSetFromOptions(check_ksp);
253: KSPSetPC(check_ksp,newpc);
254: KSPSetUp(check_ksp);
255: VecSetRandom(workv[1],NULL);
256: MatMult(local_mat,workv[1],workv[0]);
257: KSPSolve(check_ksp,workv[0],workv[0]);
258: VecAXPY(workv[0],-1.,workv[1]);
259: VecNorm(workv[0],NORM_INFINITY,&test_err);
260: KSPComputeExtremeSingularValues(check_ksp,&lambda_max,&lambda_min);
261: KSPGetIterationNumber(check_ksp,&k);
262: if (pcbddc->dbg_flag) {
263: if (isdir) {
264: 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);
265: } else {
266: 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);
267: }
268: }
269: shell_ctx->scale = 1./lambda_max;
270: KSPDestroy(&check_ksp);
271: VecDestroyVecs(2,&workv);
272: }
273: PCDestroy(&newpc);
274: return(0);
275: }
277: PetscErrorCode PCBDDCNullSpaceCheckCorrection(PC pc, PetscBool isdir)
278: {
279: PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
280: Mat_IS *matis = (Mat_IS*)pc->pmat->data;
281: KSP check_ksp,local_ksp;
282: MatNullSpace NullSpace = NULL;
283: NullSpaceCorrection_ctx shell_ctx;
284: PC check_pc;
285: Mat test_mat,local_mat;
286: PetscReal test_err,lambda_min,lambda_max;
287: Vec work1,work2,work3;
288: PetscInt k;
289: const char *prefix;
290: PetscErrorCode ierr;
293: MatGetNullSpace(matis->A,&NullSpace);
294: if (!NullSpace) {
295: MatGetNearNullSpace(matis->A,&NullSpace);
296: }
297: if (!NullSpace) {
298: return(0);
299: }
300: if (!pcbddc->dbg_flag) return(0);
301: if (isdir) local_ksp = pcbddc->ksp_D;
302: else local_ksp = pcbddc->ksp_R;
303: KSPGetOperators(local_ksp,&local_mat,NULL);
304: KSPGetPC(local_ksp,&check_pc);
305: PCShellGetContext(check_pc,(void**)&shell_ctx);
306: VecDuplicate(shell_ctx->work_full_1,&work1);
307: VecDuplicate(shell_ctx->work_full_1,&work2);
308: VecDuplicate(shell_ctx->work_full_1,&work3);
310: VecSetRandom(shell_ctx->work_small_1,NULL);
311: MatMult(shell_ctx->basis_mat,shell_ctx->work_small_1,work1);
312: VecCopy(work1,work2);
313: MatMult(local_mat,work1,work3);
314: PCApply(check_pc,work3,work1);
315: VecAXPY(work1,-1.,work2);
316: VecNorm(work1,NORM_INFINITY,&test_err);
317: if (isdir) {
318: PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Dirichlet nullspace correction solver: %1.14e\n",PetscGlobalRank,test_err);
319: } else {
320: PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Neumann nullspace correction solver: %1.14e\n",PetscGlobalRank,test_err);
321: }
323: MatTransposeMatMult(shell_ctx->Lbasis_mat,shell_ctx->Kbasis_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&test_mat);
324: MatShift(test_mat,1.);
325: MatNorm(test_mat,NORM_INFINITY,&test_err);
326: MatDestroy(&test_mat);
327: if (isdir) {
328: PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Dirichlet nullspace matrices: %1.14e\n",PetscGlobalRank,test_err);
329: } else {
330: PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Neumann nullspace matrices: %1.14e\n",PetscGlobalRank,test_err);
331: }
333: /* Create ksp object suitable for extreme eigenvalues' estimation */
334: KSPCreate(PETSC_COMM_SELF,&check_ksp);
335: PetscObjectIncrementTabLevel((PetscObject)check_ksp,(PetscObject)local_ksp,1);
336: KSPGetOptionsPrefix(local_ksp,&prefix);
337: KSPSetOptionsPrefix(check_ksp,prefix);
338: KSPAppendOptionsPrefix(check_ksp,"approxcheck_");
339: KSPSetErrorIfNotConverged(check_ksp,pc->erroriffailure);
340: KSPSetOperators(check_ksp,local_mat,local_mat);
341: KSPSetTolerances(check_ksp,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,PETSC_DEFAULT);
342: KSPSetComputeSingularValues(check_ksp,PETSC_TRUE);
343: KSPSetFromOptions(check_ksp);
344: KSPSetPC(check_ksp,check_pc);
345: KSPSetUp(check_ksp);
346: VecSetRandom(work1,NULL);
347: MatMult(local_mat,work1,work2);
348: KSPSolve(check_ksp,work2,work2);
349: VecAXPY(work2,-1.,work1);
350: VecNorm(work2,NORM_INFINITY,&test_err);
351: KSPComputeExtremeSingularValues(check_ksp,&lambda_max,&lambda_min);
352: KSPGetIterationNumber(check_ksp,&k);
353: if (isdir) {
354: PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Dirichlet adapted KSP (scale %d) %1.14e (it %d, eigs %1.6e %1.6e)\n",PetscGlobalRank,shell_ctx->apply_scaling,test_err,k,lambda_min,lambda_max);
355: } else {
356: PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Neumann adapted KSP (scale %d) %1.14e (it %d, eigs %1.6e %1.6e)\n",PetscGlobalRank,shell_ctx->apply_scaling,test_err,k,lambda_min,lambda_max);
357: }
358: KSPDestroy(&check_ksp);
359: VecDestroy(&work1);
360: VecDestroy(&work2);
361: VecDestroy(&work3);
362: return(0);
363: }