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