Actual source code: ex21.c
1: static const char help[] = "Tests MatGetSchurComplement\n";
3: #include <petscksp.h>
5: PetscErrorCode Create(MPI_Comm comm, Mat *inA, IS *is0, IS *is1)
6: {
7: Mat A;
8: PetscInt r, rend, M;
9: PetscMPIInt rank;
11: PetscFunctionBeginUser;
12: *inA = 0;
13: PetscCall(MatCreate(comm, &A));
14: PetscCall(MatSetSizes(A, 4, 4, PETSC_DETERMINE, PETSC_DETERMINE));
15: PetscCall(MatSetFromOptions(A));
16: PetscCall(MatSetUp(A));
17: PetscCall(MatGetOwnershipRange(A, &r, &rend));
18: PetscCall(MatGetSize(A, &M, NULL));
20: PetscCall(ISCreateStride(comm, 2, r, 1, is0));
21: PetscCall(ISCreateStride(comm, 2, r + 2, 1, is1));
23: PetscCallMPI(MPI_Comm_rank(comm, &rank));
25: {
26: PetscInt rows[4], cols0[5], cols1[5], cols2[3], cols3[3];
27: PetscScalar RR = 1000. * rank, vals0[5], vals1[4], vals2[3], vals3[3];
29: rows[0] = r;
30: rows[1] = r + 1;
31: rows[2] = r + 2;
32: rows[3] = r + 3;
34: cols0[0] = r + 0;
35: cols0[1] = r + 1;
36: cols0[2] = r + 3;
37: cols0[3] = (r + 4) % M;
38: cols0[4] = (r + M - 4) % M;
40: cols1[0] = r + 1;
41: cols1[1] = r + 2;
42: cols1[2] = (r + 4 + 1) % M;
43: cols1[3] = (r + M - 4 + 1) % M;
45: cols2[0] = r;
46: cols2[1] = r + 2;
47: cols2[2] = (r + 4 + 2) % M;
49: cols3[0] = r + 1;
50: cols3[1] = r + 3;
51: cols3[2] = (r + 4 + 3) % M;
53: vals0[0] = RR + 1.;
54: vals0[1] = RR + 2.;
55: vals0[2] = RR + 3.;
56: vals0[3] = RR + 4.;
57: vals0[4] = RR + 5.;
59: vals1[0] = RR + 6.;
60: vals1[1] = RR + 7.;
61: vals1[2] = RR + 8.;
62: vals1[3] = RR + 9.;
64: vals2[0] = RR + 10.;
65: vals2[1] = RR + 11.;
66: vals2[2] = RR + 12.;
68: vals3[0] = RR + 13.;
69: vals3[1] = RR + 14.;
70: vals3[2] = RR + 15.;
71: PetscCall(MatSetValues(A, 1, &rows[0], 5, cols0, vals0, INSERT_VALUES));
72: PetscCall(MatSetValues(A, 1, &rows[1], 4, cols1, vals1, INSERT_VALUES));
73: PetscCall(MatSetValues(A, 1, &rows[2], 3, cols2, vals2, INSERT_VALUES));
74: PetscCall(MatSetValues(A, 1, &rows[3], 3, cols3, vals3, INSERT_VALUES));
75: }
76: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
77: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
78: *inA = A;
79: PetscFunctionReturn(PETSC_SUCCESS);
80: }
82: PetscErrorCode Destroy(Mat *A, IS *is0, IS *is1)
83: {
84: PetscFunctionBeginUser;
85: PetscCall(MatDestroy(A));
86: PetscCall(ISDestroy(is0));
87: PetscCall(ISDestroy(is1));
88: PetscFunctionReturn(PETSC_SUCCESS);
89: }
91: int main(int argc, char *argv[])
92: {
93: Mat A, S = NULL, Sexplicit = NULL, Sp, B, C;
94: MatSchurComplementAinvType ainv_type = MAT_SCHUR_COMPLEMENT_AINV_DIAG;
95: IS is0, is1;
96: PetscBool flg;
97: PetscInt m, N = 10;
99: PetscFunctionBeginUser;
100: PetscCall(PetscInitialize(&argc, &argv, 0, help));
101: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex21", "KSP");
102: 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));
103: PetscOptionsEnd();
105: /* Test the Schur complement one way */
106: PetscCall(Create(PETSC_COMM_WORLD, &A, &is0, &is1));
107: PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
108: PetscCall(ISView(is0, PETSC_VIEWER_STDOUT_WORLD));
109: PetscCall(ISView(is1, PETSC_VIEWER_STDOUT_WORLD));
110: PetscCall(MatGetSchurComplement(A, is0, is0, is1, is1, MAT_INITIAL_MATRIX, &S, ainv_type, MAT_IGNORE_MATRIX, NULL));
111: PetscCall(MatSetFromOptions(S));
112: PetscCall(MatComputeOperator(S, MATAIJ, &Sexplicit));
113: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nExplicit Schur complement of (0,0) in (1,1)\n"));
114: PetscCall(MatView(Sexplicit, PETSC_VIEWER_STDOUT_WORLD));
115: if (ainv_type == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
116: PetscCall(MatSchurComplementSetAinvType(S, MAT_SCHUR_COMPLEMENT_AINV_FULL));
117: PetscCall(MatSchurComplementGetPmat(S, MAT_INITIAL_MATRIX, &Sp));
118: PetscCall(MatMultEqual(Sp, Sexplicit, 10, &flg));
119: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Sp != S");
120: PetscCall(MatSchurComplementSetAinvType(S, MAT_SCHUR_COMPLEMENT_AINV_DIAG));
121: PetscCall(MatDestroy(&Sp));
122: }
123: PetscCall(Destroy(&A, &is0, &is1));
124: if (ainv_type == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
125: PetscCall(MatGetLocalSize(Sexplicit, &m, NULL));
126: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)Sexplicit), m, PETSC_DECIDE, PETSC_DECIDE, N, NULL, &B));
127: PetscCall(MatSetRandom(B, NULL));
128: PetscCall(MatMatMult(S, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C));
129: PetscCall(MatMatMultEqual(Sexplicit, B, C, 10, &flg));
130: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "S*B != C");
131: PetscCall(MatDestroy(&C));
132: PetscCall(MatDestroy(&B));
133: }
134: PetscCall(MatDestroy(&S));
135: PetscCall(MatDestroy(&Sexplicit));
137: /* And the other */
138: PetscCall(Create(PETSC_COMM_WORLD, &A, &is0, &is1));
139: PetscCall(MatGetSchurComplement(A, is1, is1, is0, is0, MAT_INITIAL_MATRIX, &S, ainv_type, MAT_IGNORE_MATRIX, NULL));
140: PetscCall(MatSetFromOptions(S));
141: PetscCall(MatComputeOperator(S, MATAIJ, &Sexplicit));
142: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nExplicit Schur complement of (1,1) in (0,0)\n"));
143: PetscCall(MatView(Sexplicit, PETSC_VIEWER_STDOUT_WORLD));
144: if (ainv_type == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
145: PetscCall(MatSchurComplementSetAinvType(S, MAT_SCHUR_COMPLEMENT_AINV_FULL));
146: PetscCall(MatSchurComplementGetPmat(S, MAT_INITIAL_MATRIX, &Sp));
147: PetscCall(MatMultEqual(Sp, Sexplicit, 10, &flg));
148: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Sp != S");
149: PetscCall(MatSchurComplementSetAinvType(S, MAT_SCHUR_COMPLEMENT_AINV_DIAG));
150: PetscCall(MatDestroy(&Sp));
151: }
152: PetscCall(Destroy(&A, &is0, &is1));
153: if (ainv_type == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
154: PetscCall(MatGetLocalSize(Sexplicit, &m, NULL));
155: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)Sexplicit), m, PETSC_DECIDE, PETSC_DECIDE, N, NULL, &B));
156: PetscCall(MatSetRandom(B, NULL));
157: PetscCall(MatMatMult(S, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C));
158: PetscCall(MatMatMultEqual(Sexplicit, B, C, 10, &flg));
159: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "S*B != C");
160: PetscCall(MatDestroy(&C));
161: PetscCall(MatDestroy(&B));
162: }
163: PetscCall(MatDestroy(&S));
164: PetscCall(MatDestroy(&Sexplicit));
166: /* This time just the preconditioning matrix. */
167: PetscCall(Create(PETSC_COMM_WORLD, &A, &is0, &is1));
168: PetscCall(MatGetSchurComplement(A, is0, is0, is1, is1, MAT_IGNORE_MATRIX, NULL, ainv_type, MAT_INITIAL_MATRIX, &S));
169: PetscCall(MatSetFromOptions(S));
170: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nPreconditioning Schur complement of (0,0) in (1,1)\n"));
171: PetscCall(MatView(S, PETSC_VIEWER_STDOUT_WORLD));
172: /* Modify and refresh */
173: PetscCall(MatShift(A, 1.));
174: PetscCall(MatGetSchurComplement(A, is0, is0, is1, is1, MAT_IGNORE_MATRIX, NULL, ainv_type, MAT_REUSE_MATRIX, &S));
175: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nAfter update\n"));
176: PetscCall(MatView(S, PETSC_VIEWER_STDOUT_WORLD));
177: PetscCall(Destroy(&A, &is0, &is1));
178: PetscCall(MatDestroy(&S));
180: PetscCall(PetscFinalize());
181: return 0;
182: }
184: /*TEST
185: test:
186: suffix: diag_1
187: args: -mat_schur_complement_ainv_type diag
188: nsize: 1
189: test:
190: suffix: blockdiag_1
191: args: -mat_schur_complement_ainv_type blockdiag
192: nsize: 1
193: test:
194: suffix: diag_2
195: args: -mat_schur_complement_ainv_type diag
196: nsize: 2
197: test:
198: suffix: blockdiag_2
199: args: -mat_schur_complement_ainv_type blockdiag
200: nsize: 2
201: test:
202: # does not work with single because residual norm computed by GMRES recurrence formula becomes invalid
203: requires: !single
204: suffix: diag_3
205: args: -mat_schur_complement_ainv_type diag -ksp_rtol 1e-12
206: nsize: 3
207: test:
208: # does not work with single because residual norm computed by GMRES recurrence formula becomes invalid
209: requires: !single
210: suffix: blockdiag_3
211: args: -mat_schur_complement_ainv_type blockdiag
212: nsize: 3
213: TEST*/