Actual source code: bddcnullspace.c
1: #include <petsc/private/pcbddcimpl.h>
2: #include <petsc/private/pcbddcprivateimpl.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;
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: }
28: /* E^t + small */
29: static PetscErrorCode PCBDDCNullSpaceCorrPostSolve(KSP ksp, Vec y, Vec x, void *ctx)
30: {
31: NullSpaceCorrection_ctx corr_ctx = (NullSpaceCorrection_ctx)ctx;
32: Mat K;
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: }
53: static PetscErrorCode PCBDDCNullSpaceCorrDestroy(void *ctx)
54: {
55: NullSpaceCorrection_ctx corr_ctx = (NullSpaceCorrection_ctx)ctx;
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: }
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;
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));
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];
98: PetscCall(MatIsSymmetric(local_mat, 0.0, &shell_ctx->symm));
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;
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;
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));
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));
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));
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;
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;
195: if (needscaling && pcbddc->dbg_flag) { /* test for scaling factor */
196: PC new_pc;
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));
229: if (pcbddc->dbg_flag) {
230: Vec work1, work2, work3;
231: PetscReal test_err;
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: }