Actual source code: ex21.c

  1: static const char help[] = "Tests MatGetSchurComplement\n";

  3: #include <petscksp.h>

  5: PetscErrorCode MatNormDifference(Mat A, Mat B, PetscReal *norm)
  6: {
  7:   PetscReal bnorm;

  9:   PetscFunctionBegin;
 10:   PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B));
 11:   PetscCall(MatNorm(B, NORM_FROBENIUS, &bnorm));
 12:   PetscCall(MatAXPY(B, -1.0, A, DIFFERENT_NONZERO_PATTERN));
 13:   PetscCall(MatNorm(B, NORM_FROBENIUS, norm));
 14:   PetscCall(MatDestroy(&B));
 15:   *norm = *norm / bnorm;
 16:   PetscFunctionReturn(PETSC_SUCCESS);
 17: }

 19: PetscErrorCode Create(MPI_Comm comm, Mat *inA, IS *is0, IS *is1)
 20: {
 21:   Mat         A;
 22:   PetscInt    r, rend, M;
 23:   PetscMPIInt rank;

 25:   PetscFunctionBeginUser;
 26:   *inA = 0;
 27:   PetscCall(MatCreate(comm, &A));
 28:   PetscCall(MatSetSizes(A, 4, 4, PETSC_DETERMINE, PETSC_DETERMINE));
 29:   PetscCall(MatSetFromOptions(A));
 30:   PetscCall(MatSetUp(A));
 31:   PetscCall(MatGetOwnershipRange(A, &r, &rend));
 32:   PetscCall(MatGetSize(A, &M, NULL));

 34:   PetscCall(ISCreateStride(comm, 2, r, 1, is0));
 35:   PetscCall(ISCreateStride(comm, 2, r + 2, 1, is1));

 37:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

 39:   {
 40:     PetscInt    rows[4], cols0[5], cols1[5], cols2[3], cols3[3];
 41:     PetscScalar RR = 1000. * rank, vals0[5], vals1[4], vals2[3], vals3[3];

 43:     rows[0] = r;
 44:     rows[1] = r + 1;
 45:     rows[2] = r + 2;
 46:     rows[3] = r + 3;

 48:     cols0[0] = r + 0;
 49:     cols0[1] = r + 1;
 50:     cols0[2] = r + 3;
 51:     cols0[3] = (r + 4) % M;
 52:     cols0[4] = (r + M - 4) % M;

 54:     cols1[0] = r + 1;
 55:     cols1[1] = r + 2;
 56:     cols1[2] = (r + 4 + 1) % M;
 57:     cols1[3] = (r + M - 4 + 1) % M;

 59:     cols2[0] = r;
 60:     cols2[1] = r + 2;
 61:     cols2[2] = (r + 4 + 2) % M;

 63:     cols3[0] = r + 1;
 64:     cols3[1] = r + 3;
 65:     cols3[2] = (r + 4 + 3) % M;

 67:     vals0[0] = RR + 1.;
 68:     vals0[1] = RR + 2.;
 69:     vals0[2] = RR + 3.;
 70:     vals0[3] = RR + 4.;
 71:     vals0[4] = RR + 5.;

 73:     vals1[0] = RR + 6.;
 74:     vals1[1] = RR + 7.;
 75:     vals1[2] = RR + 8.;
 76:     vals1[3] = RR + 9.;

 78:     vals2[0] = RR + 10.;
 79:     vals2[1] = RR + 11.;
 80:     vals2[2] = RR + 12.;

 82:     vals3[0] = RR + 13.;
 83:     vals3[1] = RR + 14.;
 84:     vals3[2] = RR + 15.;
 85:     PetscCall(MatSetValues(A, 1, &rows[0], 5, cols0, vals0, INSERT_VALUES));
 86:     PetscCall(MatSetValues(A, 1, &rows[1], 4, cols1, vals1, INSERT_VALUES));
 87:     PetscCall(MatSetValues(A, 1, &rows[2], 3, cols2, vals2, INSERT_VALUES));
 88:     PetscCall(MatSetValues(A, 1, &rows[3], 3, cols3, vals3, INSERT_VALUES));
 89:   }
 90:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 91:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 92:   *inA = A;
 93:   PetscFunctionReturn(PETSC_SUCCESS);
 94: }

 96: PetscErrorCode Destroy(Mat *A, IS *is0, IS *is1)
 97: {
 98:   PetscFunctionBeginUser;
 99:   PetscCall(MatDestroy(A));
100:   PetscCall(ISDestroy(is0));
101:   PetscCall(ISDestroy(is1));
102:   PetscFunctionReturn(PETSC_SUCCESS);
103: }

105: int main(int argc, char *argv[])
106: {
107:   Mat                        A, S = NULL, Sexplicit = NULL, Sp, B, C;
108:   MatSchurComplementAinvType ainv_type = MAT_SCHUR_COMPLEMENT_AINV_DIAG;
109:   IS                         is0, is1;
110:   PetscBool                  flg;
111:   PetscInt                   m, N = 10, M;

113:   PetscFunctionBeginUser;
114:   PetscCall(PetscInitialize(&argc, &argv, 0, help));
115:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex21", "KSP");
116:   PetscCall(PetscOptionsEnum("-mat_schur_complement_ainv_type", "Type of approximation for inv(A00) used when assembling Sp = A11 - A10 inv(A00) A01", "MatSchurComplementAinvType", MatSchurComplementAinvTypes, (PetscEnum)ainv_type, (PetscEnum *)&ainv_type, NULL));
117:   PetscOptionsEnd();

119:   /* Test the Schur complement one way */
120:   PetscCall(Create(PETSC_COMM_WORLD, &A, &is0, &is1));
121:   PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
122:   PetscCall(ISView(is0, PETSC_VIEWER_STDOUT_WORLD));
123:   PetscCall(ISView(is1, PETSC_VIEWER_STDOUT_WORLD));
124:   PetscCall(MatGetSchurComplement(A, is0, is0, is1, is1, MAT_INITIAL_MATRIX, &S, ainv_type, MAT_IGNORE_MATRIX, NULL));
125:   PetscCall(MatSetFromOptions(S));
126:   PetscCall(MatComputeOperator(S, MATAIJ, &Sexplicit));
127:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nExplicit Schur complement of (0,0) in (1,1)\n"));
128:   PetscCall(MatView(Sexplicit, PETSC_VIEWER_STDOUT_WORLD));
129:   if (ainv_type == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
130:     PetscCall(MatSchurComplementSetAinvType(S, MAT_SCHUR_COMPLEMENT_AINV_FULL));
131:     PetscCall(MatSchurComplementGetPmat(S, MAT_INITIAL_MATRIX, &Sp));
132:     PetscCall(MatMultEqual(Sp, Sexplicit, 10, &flg));
133:     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Sp != S");
134:     PetscCall(MatSchurComplementSetAinvType(S, MAT_SCHUR_COMPLEMENT_AINV_DIAG));
135:     PetscCall(MatDestroy(&Sp));
136:   }
137:   PetscCall(Destroy(&A, &is0, &is1));
138:   if (ainv_type == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
139:     PetscCall(MatGetLocalSize(Sexplicit, &m, NULL));
140:     PetscCall(MatCreateDense(PetscObjectComm((PetscObject)Sexplicit), m, PETSC_DECIDE, PETSC_DECIDE, N, NULL, &B));
141:     PetscCall(MatSetRandom(B, NULL));
142:     PetscCall(MatMatMult(S, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C));
143:     PetscCall(MatMatMultEqual(Sexplicit, B, C, 10, &flg));
144:     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "S*B != C");
145:     PetscCall(MatDestroy(&C));
146:     PetscCall(MatDestroy(&B));
147:   }
148:   PetscCall(MatDestroy(&S));
149:   PetscCall(MatDestroy(&Sexplicit));

151:   /* And the other */
152:   PetscCall(Create(PETSC_COMM_WORLD, &A, &is0, &is1));
153:   PetscCall(MatGetSchurComplement(A, is1, is1, is0, is0, MAT_INITIAL_MATRIX, &S, ainv_type, MAT_IGNORE_MATRIX, NULL));
154:   PetscCall(MatSetFromOptions(S));
155:   PetscCall(MatComputeOperator(S, MATAIJ, &Sexplicit));
156:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nExplicit Schur complement of (1,1) in (0,0)\n"));
157:   PetscCall(MatView(Sexplicit, PETSC_VIEWER_STDOUT_WORLD));

159:   /* Test Mat-Mat operations with B AIJ */
160:   {
161:     Mat       B, C, Ce, Cee, Cer;
162:     PetscReal err, tol = 10 * PETSC_SMALL;
163:     PetscErrorCode (*funcs[])(Mat, Mat, MatReuse, PetscReal, Mat *) = {MatMatMult, MatMatTransposeMult, MatPtAP, MatRARt};
164:     const char *names[]                                             = {"MatMatMult", "MatMatTransposeMult", "MatPtAP", "MatRARt"};
165:     PetscBool   browsacols[]                                        = {PETSC_TRUE, PETSC_FALSE, PETSC_TRUE, PETSC_FALSE};

167:     for (PetscInt i = 0; i < 4; i++) {
168:       PetscCall(MatGetLocalSize(S, NULL, &m));
169:       PetscCall(MatGetSize(S, NULL, &M));
170:       PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
171:       if (browsacols[i]) PetscCall(MatSetSizes(B, m, PETSC_DECIDE, M, 11));
172:       else PetscCall(MatSetSizes(B, PETSC_DECIDE, m, 11, M));
173:       PetscCall(MatSetType(B, MATAIJ));
174:       PetscCall(MatSeqAIJSetPreallocation(B, PETSC_DEFAULT, NULL));
175:       PetscCall(MatMPIAIJSetPreallocation(B, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL));
176:       PetscCall(MatSetUp(B));
177:       PetscCall(MatSetRandom(B, NULL));
178:       PetscCall((*funcs[i])(S, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C));
179:       PetscCall(MatComputeOperator(C, MATAIJ, &Ce));
180:       PetscCall(MatMatMult(S, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &C));
181:       PetscCall(MatComputeOperator(C, MATAIJ, &Cer));
182:       PetscCall((*funcs[i])(Sexplicit, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Cee));
183:       PetscCall(MatNormDifference(Ce, Cee, &err));
184:       PetscCheck(err <= tol, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in initial %s(): %g", names[i], (double)err);
185:       PetscCall(MatNormDifference(Cer, Cee, &err));
186:       PetscCheck(err <= tol, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in reuse %s(): %g", names[i], (double)err);
187:       PetscCall(MatDestroy(&C));
188:       PetscCall(MatDestroy(&Ce));
189:       PetscCall(MatDestroy(&Cer));
190:       PetscCall(MatDestroy(&Cee));
191:       PetscCall(MatDestroy(&B));
192:     }
193:   }
194:   if (ainv_type == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
195:     PetscCall(MatSchurComplementSetAinvType(S, MAT_SCHUR_COMPLEMENT_AINV_FULL));
196:     PetscCall(MatSchurComplementGetPmat(S, MAT_INITIAL_MATRIX, &Sp));
197:     PetscCall(MatMultEqual(Sp, Sexplicit, 10, &flg));
198:     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Sp != S");
199:     PetscCall(MatSchurComplementSetAinvType(S, MAT_SCHUR_COMPLEMENT_AINV_DIAG));
200:     PetscCall(MatDestroy(&Sp));
201:   }
202:   PetscCall(Destroy(&A, &is0, &is1));
203:   if (ainv_type == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
204:     PetscCall(MatGetLocalSize(Sexplicit, &m, NULL));
205:     PetscCall(MatCreateDense(PetscObjectComm((PetscObject)Sexplicit), m, PETSC_DECIDE, PETSC_DECIDE, N, NULL, &B));
206:     PetscCall(MatSetRandom(B, NULL));
207:     PetscCall(MatMatMult(S, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C));
208:     PetscCall(MatMatMultEqual(Sexplicit, B, C, 10, &flg));
209:     PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "S*B != C");
210:     PetscCall(MatDestroy(&C));
211:     PetscCall(MatDestroy(&B));
212:   }
213:   PetscCall(MatDestroy(&S));
214:   PetscCall(MatDestroy(&Sexplicit));

216:   /* This time just the matrix used to construct the preconditioner. */
217:   PetscCall(Create(PETSC_COMM_WORLD, &A, &is0, &is1));
218:   PetscCall(MatGetSchurComplement(A, is0, is0, is1, is1, MAT_IGNORE_MATRIX, NULL, ainv_type, MAT_INITIAL_MATRIX, &S));
219:   PetscCall(MatSetFromOptions(S));
220:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nPreconditioning Schur complement of (0,0) in (1,1)\n"));
221:   PetscCall(MatView(S, PETSC_VIEWER_STDOUT_WORLD));
222:   /* Modify and refresh */
223:   PetscCall(MatShift(A, 1.));
224:   PetscCall(MatGetSchurComplement(A, is0, is0, is1, is1, MAT_IGNORE_MATRIX, NULL, ainv_type, MAT_REUSE_MATRIX, &S));
225:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nAfter update\n"));
226:   PetscCall(MatView(S, PETSC_VIEWER_STDOUT_WORLD));
227:   PetscCall(Destroy(&A, &is0, &is1));
228:   PetscCall(MatDestroy(&S));

230:   PetscCall(PetscFinalize());
231:   return 0;
232: }

234: /*TEST
235:   test:
236:     suffix: diag_1
237:     args: -mat_schur_complement_ainv_type diag
238:     nsize: 1
239:   test:
240:     suffix: blockdiag_1
241:     args: -mat_schur_complement_ainv_type blockdiag
242:     nsize: 1
243:   test:
244:     suffix: diag_2
245:     args: -mat_schur_complement_ainv_type diag
246:     nsize: 2
247:   test:
248:     suffix: blockdiag_2
249:     args: -mat_schur_complement_ainv_type blockdiag
250:     nsize: 2
251:   test:
252:     # does not work with single because residual norm computed by GMRES recurrence formula becomes invalid
253:     requires: !single
254:     suffix: diag_3
255:     args: -mat_schur_complement_ainv_type diag -ksp_rtol 1e-12
256:     nsize: 3
257:   test:
258:     # does not work with single because residual norm computed by GMRES recurrence formula becomes invalid
259:     requires: !single
260:     suffix: blockdiag_3
261:     args: -mat_schur_complement_ainv_type blockdiag
262:     nsize: 3
263: TEST*/