xref: /petsc/src/ksp/pc/impls/bddc/bddcnullspace.c (revision 4e8208cbcbc709572b8abe32f33c78b69c819375) 
1 #include <petsc/private/pcbddcimpl.h>
2 #include <petsc/private/pcbddcprivateimpl.h>
3 #include <../src/mat/impls/dense/seq/dense.h>
4 
5 /* E + small_solve */
PCBDDCNullSpaceCorrPreSolve(KSP ksp,Vec y,Vec x,PetscCtx ctx)6 static PetscErrorCode PCBDDCNullSpaceCorrPreSolve(KSP ksp, Vec y, Vec x, PetscCtx ctx)
7 {
8   NullSpaceCorrection_ctx corr_ctx = (NullSpaceCorrection_ctx)ctx;
9   Mat                     K;
10 
11   PetscFunctionBegin;
12   PetscCall(PetscLogEventBegin(corr_ctx->evapply, ksp, 0, 0, 0));
13   PetscCall(MatMultTranspose(corr_ctx->basis_mat, y, corr_ctx->sw[0]));
14   if (corr_ctx->symm) {
15     PetscCall(MatMult(corr_ctx->inv_smat, corr_ctx->sw[0], corr_ctx->sw[1]));
16   } else {
17     PetscCall(MatMultTranspose(corr_ctx->inv_smat, corr_ctx->sw[0], corr_ctx->sw[1]));
18   }
19   PetscCall(VecScale(corr_ctx->sw[1], -1.0));
20   PetscCall(MatMult(corr_ctx->basis_mat, corr_ctx->sw[1], corr_ctx->fw[0]));
21   PetscCall(VecScale(corr_ctx->sw[1], -1.0));
22   PetscCall(KSPGetOperators(ksp, &K, NULL));
23   PetscCall(MatMultAdd(K, corr_ctx->fw[0], y, y));
24   PetscCall(PetscLogEventEnd(corr_ctx->evapply, ksp, 0, 0, 0));
25   PetscFunctionReturn(PETSC_SUCCESS);
26 }
27 
28 /* E^t + small */
PCBDDCNullSpaceCorrPostSolve(KSP ksp,Vec y,Vec x,PetscCtx ctx)29 static PetscErrorCode PCBDDCNullSpaceCorrPostSolve(KSP ksp, Vec y, Vec x, PetscCtx ctx)
30 {
31   NullSpaceCorrection_ctx corr_ctx = (NullSpaceCorrection_ctx)ctx;
32   Mat                     K;
33 
34   PetscFunctionBegin;
35   PetscCall(PetscLogEventBegin(corr_ctx->evapply, ksp, 0, 0, 0));
36   PetscCall(KSPGetOperators(ksp, &K, NULL));
37   if (corr_ctx->symm) {
38     PetscCall(MatMult(K, x, corr_ctx->fw[0]));
39   } else {
40     PetscCall(MatMultTranspose(K, x, corr_ctx->fw[0]));
41   }
42   PetscCall(MatMultTranspose(corr_ctx->basis_mat, corr_ctx->fw[0], corr_ctx->sw[0]));
43   PetscCall(VecScale(corr_ctx->sw[0], -1.0));
44   PetscCall(MatMult(corr_ctx->inv_smat, corr_ctx->sw[0], corr_ctx->sw[2]));
45   PetscCall(MatMultAdd(corr_ctx->basis_mat, corr_ctx->sw[2], x, corr_ctx->fw[0]));
46   PetscCall(VecScale(corr_ctx->fw[0], corr_ctx->scale));
47   /* Sum contributions from approximate solver and projected system */
48   PetscCall(MatMultAdd(corr_ctx->basis_mat, corr_ctx->sw[1], corr_ctx->fw[0], x));
49   PetscCall(PetscLogEventEnd(corr_ctx->evapply, ksp, 0, 0, 0));
50   PetscFunctionReturn(PETSC_SUCCESS);
51 }
52 
PCBDDCNullSpaceCorrDestroy(PetscCtxRt ctx)53 static PetscErrorCode PCBDDCNullSpaceCorrDestroy(PetscCtxRt ctx)
54 {
55   NullSpaceCorrection_ctx corr_ctx = *(NullSpaceCorrection_ctx *)ctx;
56 
57   PetscFunctionBegin;
58   PetscCall(VecDestroyVecs(3, &corr_ctx->sw));
59   PetscCall(VecDestroyVecs(1, &corr_ctx->fw));
60   PetscCall(MatDestroy(&corr_ctx->basis_mat));
61   PetscCall(MatDestroy(&corr_ctx->inv_smat));
62   PetscCall(PetscFree(corr_ctx));
63   PetscFunctionReturn(PETSC_SUCCESS);
64 }
65 
PCBDDCNullSpaceAssembleCorrection(PC pc,PetscBool isdir,PetscBool needscaling)66 PetscErrorCode PCBDDCNullSpaceAssembleCorrection(PC pc, PetscBool isdir, PetscBool needscaling)
67 {
68   PC_BDDC                *pcbddc    = (PC_BDDC *)pc->data;
69   MatNullSpace            NullSpace = NULL;
70   KSP                     local_ksp;
71   NullSpaceCorrection_ctx shell_ctx;
72   Mat                     local_mat, local_pmat, dmat, Kbasis_mat;
73   Vec                     v;
74   PetscInt                basis_size;
75   IS                      zerorows;
76   PetscBool               iscusp;
77 
78   PetscFunctionBegin;
79   if (isdir) local_ksp = pcbddc->ksp_D; /* Dirichlet solver */
80   else local_ksp = pcbddc->ksp_R;       /* Neumann solver */
81   PetscCall(KSPGetOperators(local_ksp, &local_mat, &local_pmat));
82   PetscCall(MatGetNearNullSpace(local_pmat, &NullSpace));
83   if (!NullSpace) {
84     if (pcbddc->dbg_flag) PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d doesn't have local (near) nullspace: no need for correction in %s solver \n", PetscGlobalRank, isdir ? "Dirichlet" : "Neumann"));
85     PetscFunctionReturn(PETSC_SUCCESS);
86   }
87   PetscCall(PetscObjectQuery((PetscObject)NullSpace, "_PBDDC_Null_dmat", (PetscObject *)&dmat));
88   PetscCheck(dmat, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing dense matrix");
89   PetscCall(PetscLogEventBegin(PC_BDDC_ApproxSetUp[pcbddc->current_level], pc, 0, 0, 0));
90 
91   PetscCall(PetscNew(&shell_ctx));
92   shell_ctx->scale = 1.0;
93   PetscCall(PetscObjectReference((PetscObject)dmat));
94   shell_ctx->basis_mat = dmat;
95   PetscCall(MatGetSize(dmat, NULL, &basis_size));
96   shell_ctx->evapply = PC_BDDC_ApproxApply[pcbddc->current_level];
97 
98   PetscCall(MatIsSymmetric(local_mat, 0.0, &shell_ctx->symm));
99 
100   /* explicit construct (Phi^T K Phi)^-1 */
101   PetscCall(PetscObjectTypeCompare((PetscObject)local_mat, MATSEQAIJCUSPARSE, &iscusp));
102   if (iscusp) PetscCall(MatConvert(shell_ctx->basis_mat, MATSEQDENSECUDA, MAT_INPLACE_MATRIX, &shell_ctx->basis_mat));
103   PetscCall(MatMatMult(local_mat, shell_ctx->basis_mat, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Kbasis_mat));
104   PetscCall(MatTransposeMatMult(Kbasis_mat, shell_ctx->basis_mat, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &shell_ctx->inv_smat));
105   PetscCall(MatDestroy(&Kbasis_mat));
106   PetscCall(MatBindToCPU(shell_ctx->inv_smat, PETSC_TRUE));
107   PetscCall(MatFindZeroRows(shell_ctx->inv_smat, &zerorows));
108   if (zerorows) { /* linearly dependent basis */
109     const PetscInt *idxs;
110     PetscInt        i, nz;
111 
112     PetscCall(ISGetLocalSize(zerorows, &nz));
113     PetscCall(ISGetIndices(zerorows, &idxs));
114     for (i = 0; i < nz; i++) PetscCall(MatSetValue(shell_ctx->inv_smat, idxs[i], idxs[i], 1.0, INSERT_VALUES));
115     PetscCall(ISRestoreIndices(zerorows, &idxs));
116     PetscCall(MatAssemblyBegin(shell_ctx->inv_smat, MAT_FINAL_ASSEMBLY));
117     PetscCall(MatAssemblyEnd(shell_ctx->inv_smat, MAT_FINAL_ASSEMBLY));
118   }
119   PetscCall(MatLUFactor(shell_ctx->inv_smat, NULL, NULL, NULL));
120   PetscCall(MatSeqDenseInvertFactors_Private(shell_ctx->inv_smat));
121   if (zerorows) { /* linearly dependent basis */
122     const PetscInt *idxs;
123     PetscInt        i, nz;
124 
125     PetscCall(ISGetLocalSize(zerorows, &nz));
126     PetscCall(ISGetIndices(zerorows, &idxs));
127     for (i = 0; i < nz; i++) PetscCall(MatSetValue(shell_ctx->inv_smat, idxs[i], idxs[i], 0.0, INSERT_VALUES));
128     PetscCall(ISRestoreIndices(zerorows, &idxs));
129     PetscCall(MatAssemblyBegin(shell_ctx->inv_smat, MAT_FINAL_ASSEMBLY));
130     PetscCall(MatAssemblyEnd(shell_ctx->inv_smat, MAT_FINAL_ASSEMBLY));
131   }
132   PetscCall(ISDestroy(&zerorows));
133 
134   /* Create work vectors in shell context */
135   PetscCall(MatCreateVecs(shell_ctx->inv_smat, &v, NULL));
136   PetscCall(KSPCreateVecs(local_ksp, 1, &shell_ctx->fw, 0, NULL));
137   PetscCall(VecDuplicateVecs(v, 3, &shell_ctx->sw));
138   PetscCall(VecDestroy(&v));
139 
140   /* add special pre/post solve to KSP (see [1], eq. 48) */
141   PetscCall(KSPSetPreSolve(local_ksp, PCBDDCNullSpaceCorrPreSolve, shell_ctx));
142   PetscCall(KSPSetPostSolve(local_ksp, PCBDDCNullSpaceCorrPostSolve, shell_ctx));
143   PetscCall(PetscObjectContainerCompose((PetscObject)local_ksp, "_PCBDDC_Null_PrePost_ctx", shell_ctx, PCBDDCNullSpaceCorrDestroy));
144 
145   /* Create ksp object suitable for extreme eigenvalues' estimation */
146   if (needscaling || pcbddc->dbg_flag) {
147     KSP         check_ksp;
148     PC          local_pc;
149     Vec         work1, work2;
150     const char *prefix;
151     PetscReal   test_err, lambda_min, lambda_max;
152     PetscInt    k, maxit;
153     PetscBool   isspd, isset;
154 
155     PetscCall(VecDuplicate(shell_ctx->fw[0], &work1));
156     PetscCall(VecDuplicate(shell_ctx->fw[0], &work2));
157     PetscCall(KSPCreate(PETSC_COMM_SELF, &check_ksp));
158     PetscCall(KSPSetNestLevel(check_ksp, pc->kspnestlevel));
159     PetscCall(MatIsSPDKnown(local_mat, &isset, &isspd));
160     if (isset && isspd) PetscCall(KSPSetType(check_ksp, KSPCG));
161     PetscCall(PetscObjectIncrementTabLevel((PetscObject)check_ksp, (PetscObject)local_ksp, 0));
162     PetscCall(KSPGetOptionsPrefix(local_ksp, &prefix));
163     PetscCall(KSPSetOptionsPrefix(check_ksp, prefix));
164     PetscCall(KSPAppendOptionsPrefix(check_ksp, "approximate_scale_"));
165     PetscCall(KSPSetErrorIfNotConverged(check_ksp, PETSC_FALSE));
166     PetscCall(KSPSetOperators(check_ksp, local_mat, local_pmat));
167     PetscCall(KSPSetComputeSingularValues(check_ksp, PETSC_TRUE));
168     PetscCall(KSPSetPreSolve(check_ksp, PCBDDCNullSpaceCorrPreSolve, shell_ctx));
169     PetscCall(KSPSetPostSolve(check_ksp, PCBDDCNullSpaceCorrPostSolve, shell_ctx));
170     PetscCall(KSPSetTolerances(check_ksp, PETSC_SMALL, PETSC_SMALL, PETSC_CURRENT, PETSC_CURRENT));
171     PetscCall(KSPSetFromOptions(check_ksp));
172     /* setup with default maxit, then set maxit to min(10,any_set_from_command_line) (bug in computing eigenvalues when changing the number of iterations */
173     PetscCall(KSPSetUp(check_ksp));
174     PetscCall(KSPGetPC(local_ksp, &local_pc));
175     PetscCall(KSPSetPC(check_ksp, local_pc));
176     PetscCall(KSPGetTolerances(check_ksp, NULL, NULL, NULL, &maxit));
177     PetscCall(KSPSetTolerances(check_ksp, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, PetscMin(10, maxit)));
178     PetscCall(VecSetRandom(work2, NULL));
179     PetscCall(MatMult(local_mat, work2, work1));
180     PetscCall(KSPSolve(check_ksp, work1, work1));
181     PetscCall(KSPCheckSolve(check_ksp, pc, work1));
182     PetscCall(VecAXPY(work1, -1., work2));
183     PetscCall(VecNorm(work1, NORM_INFINITY, &test_err));
184     PetscCall(KSPComputeExtremeSingularValues(check_ksp, &lambda_max, &lambda_min));
185     PetscCall(KSPGetIterationNumber(check_ksp, &k));
186     if (pcbddc->dbg_flag) {
187       if (isdir) {
188         PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d infinity error for Dirichlet adapted solver (no scale) %1.14e (it %" PetscInt_FMT ", eigs %1.6e %1.6e)\n", PetscGlobalRank, (double)test_err, k, (double)lambda_min, (double)lambda_max));
189       } else {
190         PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d infinity error for Neumann adapted solver (no scale) %1.14e (it %" PetscInt_FMT ", eigs %1.6e %1.6e)\n", PetscGlobalRank, (double)test_err, k, (double)lambda_min, (double)lambda_max));
191       }
192     }
193     if (needscaling) shell_ctx->scale = 1.0 / lambda_max;
194 
195     if (needscaling && pcbddc->dbg_flag) { /* test for scaling factor */
196       PC new_pc;
197 
198       PetscCall(VecSetRandom(work2, NULL));
199       PetscCall(MatMult(local_mat, work2, work1));
200       PetscCall(PCCreate(PetscObjectComm((PetscObject)check_ksp), &new_pc));
201       PetscCall(PCSetType(new_pc, PCKSP));
202       PetscCall(PCSetOperators(new_pc, local_mat, local_pmat));
203       PetscCall(PCKSPSetKSP(new_pc, local_ksp));
204       PetscCall(KSPSetPC(check_ksp, new_pc));
205       PetscCall(PCDestroy(&new_pc));
206       PetscCall(KSPSetTolerances(check_ksp, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, maxit));
207       PetscCall(KSPSetPreSolve(check_ksp, NULL, NULL));
208       PetscCall(KSPSetPostSolve(check_ksp, NULL, NULL));
209       PetscCall(KSPSolve(check_ksp, work1, work1));
210       PetscCall(KSPCheckSolve(check_ksp, pc, work1));
211       PetscCall(VecAXPY(work1, -1., work2));
212       PetscCall(VecNorm(work1, NORM_INFINITY, &test_err));
213       PetscCall(KSPComputeExtremeSingularValues(check_ksp, &lambda_max, &lambda_min));
214       PetscCall(KSPGetIterationNumber(check_ksp, &k));
215       if (pcbddc->dbg_flag) {
216         if (isdir) {
217           PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d infinity error for Dirichlet adapted solver (scale %g) %1.14e (it %" PetscInt_FMT ", eigs %1.6e %1.6e)\n", PetscGlobalRank, (double)PetscRealPart(shell_ctx->scale), (double)test_err, k, (double)lambda_min, (double)lambda_max));
218         } else {
219           PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d infinity error for Neumann adapted solver (scale %g) %1.14e (it %" PetscInt_FMT ", eigs %1.6e %1.6e)\n", PetscGlobalRank, (double)PetscRealPart(shell_ctx->scale), (double)test_err, k, (double)lambda_min, (double)lambda_max));
220         }
221       }
222     }
223     PetscCall(KSPDestroy(&check_ksp));
224     PetscCall(VecDestroy(&work1));
225     PetscCall(VecDestroy(&work2));
226   }
227   PetscCall(PetscLogEventEnd(PC_BDDC_ApproxSetUp[pcbddc->current_level], pc, 0, 0, 0));
228 
229   if (pcbddc->dbg_flag) {
230     Vec       work1, work2, work3;
231     PetscReal test_err;
232 
233     /* check nullspace basis is solved exactly */
234     PetscCall(VecDuplicate(shell_ctx->fw[0], &work1));
235     PetscCall(VecDuplicate(shell_ctx->fw[0], &work2));
236     PetscCall(VecDuplicate(shell_ctx->fw[0], &work3));
237     PetscCall(VecSetRandom(shell_ctx->sw[0], NULL));
238     PetscCall(MatMult(shell_ctx->basis_mat, shell_ctx->sw[0], work1));
239     PetscCall(VecCopy(work1, work2));
240     PetscCall(MatMult(local_mat, work1, work3));
241     PetscCall(KSPSolve(local_ksp, work3, work1));
242     PetscCall(VecAXPY(work1, -1., work2));
243     PetscCall(VecNorm(work1, NORM_INFINITY, &test_err));
244     if (isdir) {
245       PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d infinity error for Dirichlet nullspace correction solver: %1.14e\n", PetscGlobalRank, (double)test_err));
246     } else {
247       PetscCall(PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer, "Subdomain %04d infinity error for Neumann nullspace correction solver: %1.14e\n", PetscGlobalRank, (double)test_err));
248     }
249     PetscCall(VecDestroy(&work1));
250     PetscCall(VecDestroy(&work2));
251     PetscCall(VecDestroy(&work3));
252   }
253   PetscFunctionReturn(PETSC_SUCCESS);
254 }
255