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