Actual source code: bddcnullspace.c
petsc-3.6.4 2016-04-12
1: #include <../src/ksp/pc/impls/bddc/bddc.h>
2: #include <../src/ksp/pc/impls/bddc/bddcprivate.h>
6: PetscErrorCode PCBDDCNullSpaceAssembleCoarse(PC pc, Mat coarse_mat, MatNullSpace* CoarseNullSpace)
7: {
8: PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
9: Mat_IS *matis = (Mat_IS*)pc->pmat->data;
10: MatNullSpace tempCoarseNullSpace=NULL;
11: const Vec *nsp_vecs;
12: Vec *coarse_nsp_vecs,local_vec,local_primal_vec,wcoarse_vec,wcoarse_rhs;
13: PetscInt nsp_size,coarse_nsp_size,i;
14: PetscBool nsp_has_cnst;
15: PetscReal test_null;
19: tempCoarseNullSpace = 0;
20: coarse_nsp_size = 0;
21: coarse_nsp_vecs = 0;
22: MatNullSpaceGetVecs(pcbddc->NullSpace,&nsp_has_cnst,&nsp_size,&nsp_vecs);
23: if (coarse_mat) {
24: PetscMalloc1(nsp_size+1,&coarse_nsp_vecs);
25: for (i=0;i<nsp_size+1;i++) {
26: MatCreateVecs(coarse_mat,&coarse_nsp_vecs[i],NULL);
27: }
28: if (pcbddc->dbg_flag) {
29: MatCreateVecs(coarse_mat,&wcoarse_vec,&wcoarse_rhs);
30: }
31: }
32: MatCreateVecs(pcbddc->ConstraintMatrix,&local_vec,&local_primal_vec);
33: if (nsp_has_cnst) {
34: VecSet(local_vec,1.0);
35: MatMult(pcbddc->ConstraintMatrix,local_vec,local_primal_vec);
36: VecScatterBegin(pcbddc->coarse_loc_to_glob,local_primal_vec,pcbddc->coarse_vec,INSERT_VALUES,SCATTER_FORWARD);
37: VecScatterEnd(pcbddc->coarse_loc_to_glob,local_primal_vec,pcbddc->coarse_vec,INSERT_VALUES,SCATTER_FORWARD);
38: if (coarse_mat) {
39: PetscScalar *array_out;
40: const PetscScalar *array_in;
41: PetscInt lsize;
42: if (pcbddc->dbg_flag) {
43: PetscViewer dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)coarse_mat));
44: MatMult(coarse_mat,wcoarse_vec,wcoarse_rhs);
45: VecNorm(wcoarse_rhs,NORM_INFINITY,&test_null);
46: PetscViewerASCIIPrintf(dbg_viewer,"Constant coarse null space error % 1.14e\n",test_null);
47: PetscViewerFlush(dbg_viewer);
48: }
49: VecGetLocalSize(pcbddc->coarse_vec,&lsize);
50: VecGetArrayRead(pcbddc->coarse_vec,&array_in);
51: VecGetArray(coarse_nsp_vecs[coarse_nsp_size],&array_out);
52: PetscMemcpy(array_out,array_in,lsize*sizeof(PetscScalar));
53: VecRestoreArray(coarse_nsp_vecs[coarse_nsp_size],&array_out);
54: VecRestoreArrayRead(pcbddc->coarse_vec,&array_in);
55: coarse_nsp_size++;
56: }
57: }
58: for (i=0;i<nsp_size;i++) {
59: VecScatterBegin(matis->ctx,nsp_vecs[i],local_vec,INSERT_VALUES,SCATTER_FORWARD);
60: VecScatterEnd(matis->ctx,nsp_vecs[i],local_vec,INSERT_VALUES,SCATTER_FORWARD);
61: MatMult(pcbddc->ConstraintMatrix,local_vec,local_primal_vec);
62: VecScatterBegin(pcbddc->coarse_loc_to_glob,local_primal_vec,pcbddc->coarse_vec,INSERT_VALUES,SCATTER_FORWARD);
63: VecScatterEnd(pcbddc->coarse_loc_to_glob,local_primal_vec,pcbddc->coarse_vec,INSERT_VALUES,SCATTER_FORWARD);
64: if (coarse_mat) {
65: PetscScalar *array_out;
66: const PetscScalar *array_in;
67: PetscInt lsize;
68: if (pcbddc->dbg_flag) {
69: PetscViewer dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)coarse_mat));
70: MatMult(coarse_mat,wcoarse_vec,wcoarse_rhs);
71: VecNorm(wcoarse_rhs,NORM_2,&test_null);
72: PetscViewerASCIIPrintf(dbg_viewer,"Vec %d coarse null space error % 1.14e\n",i,test_null);
73: PetscViewerFlush(dbg_viewer);
74: }
75: VecGetLocalSize(pcbddc->coarse_vec,&lsize);
76: VecGetArrayRead(pcbddc->coarse_vec,&array_in);
77: VecGetArray(coarse_nsp_vecs[coarse_nsp_size],&array_out);
78: PetscMemcpy(array_out,array_in,lsize*sizeof(PetscScalar));
79: VecRestoreArray(coarse_nsp_vecs[coarse_nsp_size],&array_out);
80: VecRestoreArrayRead(pcbddc->coarse_vec,&array_in);
81: coarse_nsp_size++;
82: }
83: }
84: if (coarse_nsp_size > 0) {
85: PCBDDCOrthonormalizeVecs(coarse_nsp_size,coarse_nsp_vecs);
86: MatNullSpaceCreate(PetscObjectComm((PetscObject)coarse_mat),PETSC_FALSE,coarse_nsp_size,coarse_nsp_vecs,&tempCoarseNullSpace);
87: for (i=0;i<nsp_size+1;i++) {
88: VecDestroy(&coarse_nsp_vecs[i]);
89: }
90: }
91: if (coarse_mat) {
92: PetscFree(coarse_nsp_vecs);
93: if (pcbddc->dbg_flag) {
94: VecDestroy(&wcoarse_vec);
95: VecDestroy(&wcoarse_rhs);
96: }
97: }
98: VecDestroy(&local_vec);
99: VecDestroy(&local_primal_vec);
100: *CoarseNullSpace = tempCoarseNullSpace;
101: return(0);
102: }
107: static PetscErrorCode PCBDDCApplyNullSpaceCorrectionPC(PC pc,Vec x,Vec y)
108: {
109: NullSpaceCorrection_ctx pc_ctx;
110: PetscErrorCode ierr;
113: PCShellGetContext(pc,(void**)&pc_ctx);
114: /* E */
115: MatMultTranspose(pc_ctx->Lbasis_mat,x,pc_ctx->work_small_2);
116: MatMultAdd(pc_ctx->Kbasis_mat,pc_ctx->work_small_2,x,pc_ctx->work_full_1);
117: /* P^-1 */
118: PCApply(pc_ctx->local_pc,pc_ctx->work_full_1,pc_ctx->work_full_2);
119: /* E^T */
120: MatMultTranspose(pc_ctx->Kbasis_mat,pc_ctx->work_full_2,pc_ctx->work_small_1);
121: VecScale(pc_ctx->work_small_1,-1.0);
122: MatMultAdd(pc_ctx->Lbasis_mat,pc_ctx->work_small_1,pc_ctx->work_full_2,pc_ctx->work_full_1);
123: /* Sum contributions */
124: MatMultAdd(pc_ctx->basis_mat,pc_ctx->work_small_2,pc_ctx->work_full_1,y);
125: return(0);
126: }
130: static PetscErrorCode PCBDDCDestroyNullSpaceCorrectionPC(PC pc)
131: {
132: NullSpaceCorrection_ctx pc_ctx;
133: PetscErrorCode ierr;
136: PCShellGetContext(pc,(void**)&pc_ctx);
137: VecDestroy(&pc_ctx->work_small_1);
138: VecDestroy(&pc_ctx->work_small_2);
139: VecDestroy(&pc_ctx->work_full_1);
140: VecDestroy(&pc_ctx->work_full_2);
141: MatDestroy(&pc_ctx->basis_mat);
142: MatDestroy(&pc_ctx->Lbasis_mat);
143: MatDestroy(&pc_ctx->Kbasis_mat);
144: PCDestroy(&pc_ctx->local_pc);
145: PetscFree(pc_ctx);
146: return(0);
147: }
149: /*
150: PETSC_EXTERN PetscErrorCode PCBDDCApplyNullSpaceCorrectionPC(PC,Vec,Vec);
151: PETSC_EXTERN PetscErrorCode PCBDDCDestroyNullSpaceCorrectionPC(PC);
152: */
156: PetscErrorCode PCBDDCNullSpaceAssembleCorrection(PC pc, PetscBool isdir, IS local_dofs)
157: {
158: PC_BDDC *pcbddc = (PC_BDDC*)pc->data;
159: PC_IS *pcis = (PC_IS*)pc->data;
160: Mat_IS* matis = (Mat_IS*)pc->pmat->data;
161: KSP local_ksp;
162: PC newpc;
163: NullSpaceCorrection_ctx shell_ctx;
164: Mat local_mat,local_pmat,small_mat,inv_small_mat;
165: Vec work1,work2;
166: const Vec *nullvecs;
167: VecScatter scatter_ctx;
168: IS is_aux;
169: MatFactorInfo matinfo;
170: PetscScalar *basis_mat,*Kbasis_mat,*array,*array_mat;
171: PetscScalar one = 1.0,zero = 0.0, m_one = -1.0;
172: PetscInt basis_dofs,basis_size,nnsp_size,i,k;
173: PetscBool nnsp_has_cnst;
177: /* Infer the local solver */
178: ISGetSize(local_dofs,&basis_dofs);
179: if (isdir) {
180: /* Dirichlet solver */
181: local_ksp = pcbddc->ksp_D;
182: } else {
183: /* Neumann solver */
184: local_ksp = pcbddc->ksp_R;
185: }
186: KSPGetOperators(local_ksp,&local_mat,&local_pmat);
188: /* Get null space vecs */
189: MatNullSpaceGetVecs(pcbddc->NullSpace,&nnsp_has_cnst,&nnsp_size,&nullvecs);
190: basis_size = nnsp_size;
191: if (nnsp_has_cnst) {
192: basis_size++;
193: }
195: if (basis_dofs) {
196: /* Create shell ctx */
197: PetscNew(&shell_ctx);
199: /* Create work vectors in shell context */
200: VecCreate(PETSC_COMM_SELF,&shell_ctx->work_small_1);
201: VecSetSizes(shell_ctx->work_small_1,basis_size,basis_size);
202: VecSetType(shell_ctx->work_small_1,VECSEQ);
203: VecDuplicate(shell_ctx->work_small_1,&shell_ctx->work_small_2);
204: VecCreate(PETSC_COMM_SELF,&shell_ctx->work_full_1);
205: VecSetSizes(shell_ctx->work_full_1,basis_dofs,basis_dofs);
206: VecSetType(shell_ctx->work_full_1,VECSEQ);
207: VecDuplicate(shell_ctx->work_full_1,&shell_ctx->work_full_2);
209: /* Allocate workspace */
210: MatCreateSeqDense(PETSC_COMM_SELF,basis_dofs,basis_size,NULL,&shell_ctx->basis_mat );
211: MatCreateSeqDense(PETSC_COMM_SELF,basis_dofs,basis_size,NULL,&shell_ctx->Kbasis_mat);
212: MatDenseGetArray(shell_ctx->basis_mat,&basis_mat);
213: MatDenseGetArray(shell_ctx->Kbasis_mat,&Kbasis_mat);
215: /* Restrict local null space on selected dofs (Dirichlet or Neumann)
216: and compute matrices N and K*N */
217: VecDuplicate(shell_ctx->work_full_1,&work1);
218: VecDuplicate(shell_ctx->work_full_1,&work2);
219: VecScatterCreate(pcis->vec1_N,local_dofs,work1,(IS)0,&scatter_ctx);
220: }
222: for (k=0;k<nnsp_size;k++) {
223: VecScatterBegin(matis->ctx,nullvecs[k],pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);
224: VecScatterEnd(matis->ctx,nullvecs[k],pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);
225: if (basis_dofs) {
226: VecPlaceArray(work1,(const PetscScalar*)&basis_mat[k*basis_dofs]);
227: VecScatterBegin(scatter_ctx,pcis->vec1_N,work1,INSERT_VALUES,SCATTER_FORWARD);
228: VecScatterEnd(scatter_ctx,pcis->vec1_N,work1,INSERT_VALUES,SCATTER_FORWARD);
229: VecPlaceArray(work2,(const PetscScalar*)&Kbasis_mat[k*basis_dofs]);
230: MatMult(local_mat,work1,work2);
231: VecResetArray(work1);
232: VecResetArray(work2);
233: }
234: }
236: if (basis_dofs) {
237: if (nnsp_has_cnst) {
238: VecPlaceArray(work1,(const PetscScalar*)&basis_mat[k*basis_dofs]);
239: VecSet(work1,one);
240: VecPlaceArray(work2,(const PetscScalar*)&Kbasis_mat[k*basis_dofs]);
241: MatMult(local_mat,work1,work2);
242: VecResetArray(work1);
243: VecResetArray(work2);
244: }
245: VecDestroy(&work1);
246: VecDestroy(&work2);
247: VecScatterDestroy(&scatter_ctx);
248: MatDenseRestoreArray(shell_ctx->basis_mat,&basis_mat);
249: MatDenseRestoreArray(shell_ctx->Kbasis_mat,&Kbasis_mat);
251: /* Assemble another Mat object in shell context */
252: MatTransposeMatMult(shell_ctx->basis_mat,shell_ctx->Kbasis_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&small_mat);
253: MatFactorInfoInitialize(&matinfo);
254: ISCreateStride(PETSC_COMM_SELF,basis_size,0,1,&is_aux);
255: MatLUFactor(small_mat,is_aux,is_aux,&matinfo);
256: ISDestroy(&is_aux);
257: PetscMalloc1(basis_size*basis_size,&array_mat);
258: for (k=0;k<basis_size;k++) {
259: VecSet(shell_ctx->work_small_1,zero);
260: VecSetValue(shell_ctx->work_small_1,k,one,INSERT_VALUES);
261: VecAssemblyBegin(shell_ctx->work_small_1);
262: VecAssemblyEnd(shell_ctx->work_small_1);
263: MatSolve(small_mat,shell_ctx->work_small_1,shell_ctx->work_small_2);
264: VecGetArrayRead(shell_ctx->work_small_2,(const PetscScalar**)&array);
265: for (i=0;i<basis_size;i++) {
266: array_mat[i*basis_size+k]=array[i];
267: }
268: VecRestoreArrayRead(shell_ctx->work_small_2,(const PetscScalar**)&array);
269: }
270: MatCreateSeqDense(PETSC_COMM_SELF,basis_size,basis_size,array_mat,&inv_small_mat);
271: MatMatMult(shell_ctx->basis_mat,inv_small_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&shell_ctx->Lbasis_mat);
272: PetscFree(array_mat);
273: MatDestroy(&inv_small_mat);
274: MatDestroy(&small_mat);
275: MatScale(shell_ctx->Kbasis_mat,m_one);
277: /* Rebuild local PC */
278: KSPGetPC(local_ksp,&shell_ctx->local_pc);
279: PetscObjectReference((PetscObject)shell_ctx->local_pc);
280: PCCreate(PETSC_COMM_SELF,&newpc);
281: PCSetOperators(newpc,local_mat,local_mat);
282: PCSetType(newpc,PCSHELL);
283: PCShellSetContext(newpc,shell_ctx);
284: PCShellSetApply(newpc,PCBDDCApplyNullSpaceCorrectionPC);
285: PCShellSetDestroy(newpc,PCBDDCDestroyNullSpaceCorrectionPC);
286: PCSetUp(newpc);
287: KSPSetPC(local_ksp,newpc);
288: PCDestroy(&newpc);
289: KSPSetUp(local_ksp);
290: }
291: /* test */
292: if (pcbddc->dbg_flag && basis_dofs) {
293: KSP check_ksp;
294: PC check_pc;
295: Mat test_mat;
296: Vec work3;
297: PetscReal test_err,lambda_min,lambda_max;
298: PetscBool setsym,issym=PETSC_FALSE;
299: PetscInt tabs;
301: PetscViewerASCIIGetTab(pcbddc->dbg_viewer,&tabs);
302: KSPGetPC(local_ksp,&check_pc);
303: VecDuplicate(shell_ctx->work_full_1,&work1);
304: VecDuplicate(shell_ctx->work_full_1,&work2);
305: VecDuplicate(shell_ctx->work_full_1,&work3);
306: VecSetRandom(shell_ctx->work_small_1,NULL);
307: MatMult(shell_ctx->basis_mat,shell_ctx->work_small_1,work1);
308: VecCopy(work1,work2);
309: MatMult(local_mat,work1,work3);
310: PCApply(check_pc,work3,work1);
311: VecAXPY(work1,m_one,work2);
312: VecNorm(work1,NORM_INFINITY,&test_err);
313: PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d error for nullspace correction for ",PetscGlobalRank);
314: PetscViewerASCIIUseTabs(pcbddc->dbg_viewer,PETSC_FALSE);
315: if (isdir) {
316: PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Dirichlet ");
317: } else {
318: PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Neumann ");
319: }
320: PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"solver is :%1.14e\n",test_err);
321: PetscViewerASCIISetTab(pcbddc->dbg_viewer,tabs);
322: PetscViewerASCIIUseTabs(pcbddc->dbg_viewer,PETSC_TRUE);
324: MatTransposeMatMult(shell_ctx->Lbasis_mat,shell_ctx->Kbasis_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&test_mat);
325: MatShift(test_mat,one);
326: MatNorm(test_mat,NORM_INFINITY,&test_err);
327: MatDestroy(&test_mat);
328: PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d error for nullspace matrices is :%1.14e\n",PetscGlobalRank,test_err);
330: /* Create ksp object suitable for extreme eigenvalues' estimation */
331: KSPCreate(PETSC_COMM_SELF,&check_ksp);
332: KSPSetErrorIfNotConverged(check_ksp,pc->erroriffailure);
333: KSPSetOperators(check_ksp,local_mat,local_mat);
334: KSPSetTolerances(check_ksp,1.e-8,1.e-8,PETSC_DEFAULT,basis_dofs);
335: KSPSetComputeSingularValues(check_ksp,PETSC_TRUE);
336: MatIsSymmetricKnown(pc->pmat,&setsym,&issym);
337: if (issym) {
338: KSPSetType(check_ksp,KSPCG);
339: }
340: KSPSetPC(check_ksp,check_pc);
341: KSPSetUp(check_ksp);
342: VecSetRandom(work1,NULL);
343: MatMult(local_mat,work1,work2);
344: KSPSolve(check_ksp,work2,work2);
345: VecAXPY(work2,m_one,work1);
346: VecNorm(work2,NORM_INFINITY,&test_err);
347: KSPComputeExtremeSingularValues(check_ksp,&lambda_max,&lambda_min);
348: KSPGetIterationNumber(check_ksp,&k);
349: PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d error for adapted KSP %1.14e (it %d, eigs %1.6e %1.6e)\n",PetscGlobalRank,test_err,k,lambda_min,lambda_max);
350: KSPDestroy(&check_ksp);
351: VecDestroy(&work1);
352: VecDestroy(&work2);
353: VecDestroy(&work3);
354: }
355: /* all processes shoud call this, even the void ones */
356: if (pcbddc->dbg_flag) {
357: PetscViewerFlush(pcbddc->dbg_viewer);
358: }
359: return(0);
360: }
364: PetscErrorCode PCBDDCNullSpaceAdaptGlobal(PC pc)
365: {
366: PC_IS* pcis = (PC_IS*)(pc->data);
367: PC_BDDC* pcbddc = (PC_BDDC*)(pc->data);
368: KSP inv_change;
369: const Vec *nsp_vecs;
370: Vec *new_nsp_vecs;
371: PetscInt i,nsp_size,new_nsp_size,start_new;
372: PetscBool nsp_has_cnst;
373: MatNullSpace new_nsp;
377: /* create KSP for change of basis */
378: MatGetSize(pcbddc->ChangeOfBasisMatrix,&i,NULL);
379: KSPCreate(PetscObjectComm((PetscObject)pc),&inv_change);
380: KSPSetErrorIfNotConverged(inv_change,pc->erroriffailure);
381: KSPSetOperators(inv_change,pcbddc->ChangeOfBasisMatrix,pcbddc->ChangeOfBasisMatrix);
382: KSPSetTolerances(inv_change,1.e-8,1.e-8,PETSC_DEFAULT,2*i);
383: if (pcbddc->dbg_flag) {
384: KSPMonitorSet(inv_change,KSPMonitorDefault,pcbddc->dbg_viewer,NULL);
385: }
386: KSPSetUp(inv_change);
388: /* get nullspace and transform it */
389: MatNullSpaceGetVecs(pcbddc->NullSpace,&nsp_has_cnst,&nsp_size,&nsp_vecs);
390: new_nsp_size = nsp_size;
391: if (nsp_has_cnst) {
392: new_nsp_size++;
393: }
394: VecDuplicateVecs(pcis->vec1_global,new_nsp_size,&new_nsp_vecs);
396: start_new = 0;
397: if (nsp_has_cnst) {
398: start_new = 1;
399: VecSet(new_nsp_vecs[0],1.0);
400: if (pcbddc->dbg_flag) {
401: PetscViewerFlush(pcbddc->dbg_viewer);
402: PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"Mapping constant in nullspace\n");
403: }
404: KSPSolve(inv_change,new_nsp_vecs[0],new_nsp_vecs[0]);
405: }
406: for (i=0;i<nsp_size;i++) {
407: PetscViewerFlush(pcbddc->dbg_viewer);
408: PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"Mapping %dth vector in nullspace\n",i);
409: KSPSolve(inv_change,nsp_vecs[i],new_nsp_vecs[i+start_new]);
410: }
411: PCBDDCOrthonormalizeVecs(new_nsp_size,new_nsp_vecs);
412: MatNullSpaceCreate(PetscObjectComm((PetscObject)pc),PETSC_FALSE,new_nsp_size,new_nsp_vecs,&new_nsp);
413: PCBDDCSetNullSpace(pc,new_nsp);
415: /* free */
416: KSPDestroy(&inv_change);
417: MatNullSpaceDestroy(&new_nsp);
418: VecDestroyVecs(new_nsp_size,&new_nsp_vecs);
420: /* check */
421: if (pcbddc->dbg_flag) {
422: PetscBool nsp_t=PETSC_FALSE;
423: Mat temp_mat;
424: Mat_IS* matis = (Mat_IS*)pc->pmat->data;
426: temp_mat = matis->A;
427: matis->A = pcbddc->local_mat;
428: pcbddc->local_mat = temp_mat;
429: MatNullSpaceTest(pcbddc->NullSpace,pc->pmat,&nsp_t);
430: PetscPrintf(PetscObjectComm((PetscObject)(pc->pmat)),"Check nullspace with change of basis: %d\n",nsp_t);
431: temp_mat = matis->A;
432: matis->A = pcbddc->local_mat;
433: pcbddc->local_mat = temp_mat;
434: }
435: return(0);
436: }