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