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;

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

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

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

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

 36:   VecDestroy(&diag);

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

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

 52:   /* A21 */
 53:   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:   MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, &tmp[0][0], &A);
 64:   MatNestSetVecType(A, VECNEST);
 65:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 66:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

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

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

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

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

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

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

 94:   KSPSolve(ksp, b, x);

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

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

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

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

117:   MatDestroy(&A);
118:   return 0;
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;

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

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

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

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

152:   VecDestroy(&diag);

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

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

168:   /* A21 */
169:   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:   MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, &tmp[0][0], &A);
178:   MatNestSetVecType(A, VECNEST);
179:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
180:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

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

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

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

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

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

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

216:   MatDestroy(&A);
217:   return 0;
218: }

220: int main(int argc, char **args)
221: {
223:   PetscInitialize(&argc, &args, (char *)0, help);
224:   test_solve();
225:   test_solve_matgetvecs();
226:   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*/