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*/