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