Actual source code: ex22.c

  1: static const char help[] = "Test MatNest solving a linear system\n\n";

  3: #include <petscksp.h>

  5: PetscErrorCode test_solve(void)
  6: {
  7:   Mat       A11, A12, A21, A22, A, tmp[2][2];
  8:   KSP       ksp;
  9:   PC        pc;
 10:   Vec       b, x, f, h, diag, x1, x2;
 11:   Vec       tmp_x[2], *_tmp_x;
 12:   PetscInt  n, np, i, j;
 13:   PetscBool flg;

 15:   PetscFunctionBeginUser;
 16:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s \n", PETSC_FUNCTION_NAME));

 18:   n  = 3;
 19:   np = 2;
 20:   /* Create matrices */
 21:   /* A11 */
 22:   PetscCall(VecCreate(PETSC_COMM_WORLD, &diag));
 23:   PetscCall(VecSetSizes(diag, PETSC_DECIDE, n));
 24:   PetscCall(VecSetFromOptions(diag));

 26:   PetscCall(VecSet(diag, 1.0 / 10.0)); /* so inverse = diag(10) */

 28:   /* As a test, create a diagonal matrix for A11 */
 29:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A11));
 30:   PetscCall(MatSetSizes(A11, PETSC_DECIDE, PETSC_DECIDE, n, n));
 31:   PetscCall(MatSetType(A11, MATAIJ));
 32:   PetscCall(MatSeqAIJSetPreallocation(A11, n, NULL));
 33:   PetscCall(MatMPIAIJSetPreallocation(A11, np, NULL, np, NULL));
 34:   PetscCall(MatDiagonalSet(A11, diag, INSERT_VALUES));

 36:   PetscCall(VecDestroy(&diag));

 38:   /* A12 */
 39:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A12));
 40:   PetscCall(MatSetSizes(A12, PETSC_DECIDE, PETSC_DECIDE, n, np));
 41:   PetscCall(MatSetType(A12, MATAIJ));
 42:   PetscCall(MatSeqAIJSetPreallocation(A12, np, NULL));
 43:   PetscCall(MatMPIAIJSetPreallocation(A12, np, NULL, np, NULL));

 45:   for (i = 0; i < n; i++) {
 46:     for (j = 0; j < np; j++) PetscCall(MatSetValue(A12, i, j, i + j * n, INSERT_VALUES));
 47:   }
 48:   PetscCall(MatSetValue(A12, 2, 1, 4, INSERT_VALUES));
 49:   PetscCall(MatAssemblyBegin(A12, MAT_FINAL_ASSEMBLY));
 50:   PetscCall(MatAssemblyEnd(A12, MAT_FINAL_ASSEMBLY));

 52:   /* A21 */
 53:   PetscCall(MatTranspose(A12, MAT_INITIAL_MATRIX, &A21));

 55:   A22 = NULL;

 57:   /* Create block matrix */
 58:   tmp[0][0] = A11;
 59:   tmp[0][1] = A12;
 60:   tmp[1][0] = A21;
 61:   tmp[1][1] = A22;

 63:   PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, &tmp[0][0], &A));
 64:   PetscCall(MatNestSetVecType(A, VECNEST));
 65:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 66:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 68:   /* Tests MatMissingDiagonal_Nest */
 69:   PetscCall(MatMissingDiagonal(A, &flg, NULL));
 70:   if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Unexpected %s\n", flg ? "true" : "false"));

 72:   /* Create vectors */
 73:   PetscCall(MatCreateVecs(A12, &h, &f));

 75:   PetscCall(VecSet(f, 1.0));
 76:   PetscCall(VecSet(h, 0.0));

 78:   /* Create block vector */
 79:   tmp_x[0] = f;
 80:   tmp_x[1] = h;

 82:   PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, tmp_x, &b));
 83:   PetscCall(VecAssemblyBegin(b));
 84:   PetscCall(VecAssemblyEnd(b));
 85:   PetscCall(VecDuplicate(b, &x));

 87:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
 88:   PetscCall(KSPSetOperators(ksp, A, A));
 89:   PetscCall(KSPSetType(ksp, KSPGMRES));
 90:   PetscCall(KSPGetPC(ksp, &pc));
 91:   PetscCall(PCSetType(pc, PCNONE));
 92:   PetscCall(KSPSetFromOptions(ksp));

 94:   PetscCall(KSPSolve(ksp, b, x));

 96:   PetscCall(VecNestGetSubVecs(x, NULL, &_tmp_x));

 98:   x1 = _tmp_x[0];
 99:   x2 = _tmp_x[1];

101:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "x1 \n"));
102:   PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL));
103:   PetscCall(VecView(x1, PETSC_VIEWER_STDOUT_WORLD));
104:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "x2 \n"));
105:   PetscCall(VecView(x2, PETSC_VIEWER_STDOUT_WORLD));
106:   PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));

108:   PetscCall(KSPDestroy(&ksp));
109:   PetscCall(VecDestroy(&x));
110:   PetscCall(VecDestroy(&b));
111:   PetscCall(MatDestroy(&A11));
112:   PetscCall(MatDestroy(&A12));
113:   PetscCall(MatDestroy(&A21));
114:   PetscCall(VecDestroy(&f));
115:   PetscCall(VecDestroy(&h));

117:   PetscCall(MatDestroy(&A));
118:   PetscFunctionReturn(PETSC_SUCCESS);
119: }

121: PetscErrorCode test_solve_matgetvecs(void)
122: {
123:   Mat      A11, A12, A21, A;
124:   KSP      ksp;
125:   PC       pc;
126:   Vec      b, x, f, h, diag, x1, x2;
127:   PetscInt n, np, i, j;
128:   Mat      tmp[2][2];
129:   Vec     *tmp_x;

131:   PetscFunctionBeginUser;
132:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s \n", PETSC_FUNCTION_NAME));

134:   n  = 3;
135:   np = 2;
136:   /* Create matrices */
137:   /* A11 */
138:   PetscCall(VecCreate(PETSC_COMM_WORLD, &diag));
139:   PetscCall(VecSetSizes(diag, PETSC_DECIDE, n));
140:   PetscCall(VecSetFromOptions(diag));

142:   PetscCall(VecSet(diag, 1.0 / 10.0)); /* so inverse = diag(10) */

144:   /* As a test, create a diagonal matrix for A11 */
145:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A11));
146:   PetscCall(MatSetSizes(A11, PETSC_DECIDE, PETSC_DECIDE, n, n));
147:   PetscCall(MatSetType(A11, MATAIJ));
148:   PetscCall(MatSeqAIJSetPreallocation(A11, n, NULL));
149:   PetscCall(MatMPIAIJSetPreallocation(A11, np, NULL, np, NULL));
150:   PetscCall(MatDiagonalSet(A11, diag, INSERT_VALUES));

152:   PetscCall(VecDestroy(&diag));

154:   /* A12 */
155:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A12));
156:   PetscCall(MatSetSizes(A12, PETSC_DECIDE, PETSC_DECIDE, n, np));
157:   PetscCall(MatSetType(A12, MATAIJ));
158:   PetscCall(MatSeqAIJSetPreallocation(A12, np, NULL));
159:   PetscCall(MatMPIAIJSetPreallocation(A12, np, NULL, np, NULL));

161:   for (i = 0; i < n; i++) {
162:     for (j = 0; j < np; j++) PetscCall(MatSetValue(A12, i, j, i + j * n, INSERT_VALUES));
163:   }
164:   PetscCall(MatSetValue(A12, 2, 1, 4, INSERT_VALUES));
165:   PetscCall(MatAssemblyBegin(A12, MAT_FINAL_ASSEMBLY));
166:   PetscCall(MatAssemblyEnd(A12, MAT_FINAL_ASSEMBLY));

168:   /* A21 */
169:   PetscCall(MatTranspose(A12, MAT_INITIAL_MATRIX, &A21));

171:   /* Create block matrix */
172:   tmp[0][0] = A11;
173:   tmp[0][1] = A12;
174:   tmp[1][0] = A21;
175:   tmp[1][1] = NULL;

177:   PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, &tmp[0][0], &A));
178:   PetscCall(MatNestSetVecType(A, VECNEST));
179:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
180:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

182:   /* Create vectors */
183:   PetscCall(MatCreateVecs(A, &b, &x));
184:   PetscCall(VecNestGetSubVecs(b, NULL, &tmp_x));
185:   f = tmp_x[0];
186:   h = tmp_x[1];

188:   PetscCall(VecSet(f, 1.0));
189:   PetscCall(VecSet(h, 0.0));

191:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
192:   PetscCall(KSPSetOperators(ksp, A, A));
193:   PetscCall(KSPGetPC(ksp, &pc));
194:   PetscCall(PCSetType(pc, PCNONE));
195:   PetscCall(KSPSetFromOptions(ksp));

197:   PetscCall(KSPSolve(ksp, b, x));
198:   PetscCall(VecNestGetSubVecs(x, NULL, &tmp_x));
199:   x1 = tmp_x[0];
200:   x2 = tmp_x[1];

202:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "x1 \n"));
203:   PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL));
204:   PetscCall(VecView(x1, PETSC_VIEWER_STDOUT_WORLD));
205:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "x2 \n"));
206:   PetscCall(VecView(x2, PETSC_VIEWER_STDOUT_WORLD));
207:   PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));

209:   PetscCall(KSPDestroy(&ksp));
210:   PetscCall(VecDestroy(&x));
211:   PetscCall(VecDestroy(&b));
212:   PetscCall(MatDestroy(&A11));
213:   PetscCall(MatDestroy(&A12));
214:   PetscCall(MatDestroy(&A21));

216:   PetscCall(MatDestroy(&A));
217:   PetscFunctionReturn(PETSC_SUCCESS);
218: }

220: int main(int argc, char **args)
221: {
222:   PetscFunctionBeginUser;
223:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
224:   PetscCall(test_solve());
225:   PetscCall(test_solve_matgetvecs());
226:   PetscCall(PetscFinalize());
227:   return 0;
228: }

230: /*TEST

232:     test:

234:     test:
235:       suffix: 2
236:       nsize: 2

238:     test:
239:       suffix: 3
240:       nsize: 2
241:       args: -ksp_monitor_short -ksp_type bicg
242:       requires: !single

244: TEST*/