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