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++) {
 47:       MatSetValue(A12, i,j, (PetscScalar)(i+j*n), INSERT_VALUES);
 48:     }
 49:   }
 50:   MatSetValue(A12, 2,1, (PetscScalar)(4), INSERT_VALUES);
 51:   MatAssemblyBegin(A12, MAT_FINAL_ASSEMBLY);
 52:   MatAssemblyEnd(A12, MAT_FINAL_ASSEMBLY);

 54:   /* A21 */
 55:   MatTranspose(A12, MAT_INITIAL_MATRIX, &A21);

 57:   A22 = NULL;

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

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

 70:   /* Tests MatMissingDiagonal_Nest */
 71:   MatMissingDiagonal(A,&flg,NULL);
 72:   if (!flg) {
 73:     PetscPrintf(PETSC_COMM_WORLD,"Unexpected %s\n",flg ? "true" : "false");
 74:   }

 76:   /* Create vectors */
 77:   MatCreateVecs(A12, &h, &f);

 79:   VecSet(f, 1.0);
 80:   VecSet(h, 0.0);

 82:   /* Create block vector */
 83:   tmp_x[0] = f;
 84:   tmp_x[1] = h;

 86:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_x,&b);
 87:   VecAssemblyBegin(b);
 88:   VecAssemblyEnd(b);
 89:   VecDuplicate(b, &x);

 91:   KSPCreate(PETSC_COMM_WORLD, &ksp);
 92:   KSPSetOperators(ksp, A, A);
 93:   KSPSetType(ksp, KSPGMRES);
 94:   KSPGetPC(ksp, &pc);
 95:   PCSetType(pc, PCNONE);
 96:   KSPSetFromOptions(ksp);

 98:   KSPSolve(ksp, b, x);

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

102:   x1 = _tmp_x[0];
103:   x2 = _tmp_x[1];

105:   PetscPrintf(PETSC_COMM_WORLD, "x1 \n");
106:   PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL);
107:   VecView(x1, PETSC_VIEWER_STDOUT_WORLD);
108:   PetscPrintf(PETSC_COMM_WORLD, "x2 \n");
109:   VecView(x2, PETSC_VIEWER_STDOUT_WORLD);
110:   PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);

112:   KSPDestroy(&ksp);
113:   VecDestroy(&x);
114:   VecDestroy(&b);
115:   MatDestroy(&A11);
116:   MatDestroy(&A12);
117:   MatDestroy(&A21);
118:   VecDestroy(&f);
119:   VecDestroy(&h);

121:   MatDestroy(&A);
122:   return 0;
123: }

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

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

138:   n  = 3;
139:   np = 2;
140:   /* Create matrices */
141:   /* A11 */
142:   VecCreate(PETSC_COMM_WORLD, &diag);
143:   VecSetSizes(diag, PETSC_DECIDE, n);
144:   VecSetFromOptions(diag);

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

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

156:   VecDestroy(&diag);

158:   /* A12 */
159:   MatCreate(PETSC_COMM_WORLD, &A12);
160:   MatSetSizes(A12, PETSC_DECIDE, PETSC_DECIDE, n, np);
161:   MatSetType(A12, MATAIJ);
162:   MatSeqAIJSetPreallocation(A12, np, NULL);
163:   MatMPIAIJSetPreallocation(A12, np, NULL,np, NULL);

165:   for (i=0; i<n; i++) {
166:     for (j=0; j<np; j++) {
167:       MatSetValue(A12, i,j, (PetscScalar)(i+j*n), INSERT_VALUES);
168:     }
169:   }
170:   MatSetValue(A12, 2,1, (PetscScalar)(4), INSERT_VALUES);
171:   MatAssemblyBegin(A12, MAT_FINAL_ASSEMBLY);
172:   MatAssemblyEnd(A12, MAT_FINAL_ASSEMBLY);

174:   /* A21 */
175:   MatTranspose(A12, MAT_INITIAL_MATRIX, &A21);

177:   /* Create block matrix */
178:   tmp[0][0] = A11;
179:   tmp[0][1] = A12;
180:   tmp[1][0] = A21;
181:   tmp[1][1] = NULL;

183:   MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,&tmp[0][0],&A);
184:   MatNestSetVecType(A,VECNEST);
185:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
186:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

188:   /* Create vectors */
189:   MatCreateVecs(A, &b, &x);
190:   VecNestGetSubVecs(b,NULL,&tmp_x);
191:   f    = tmp_x[0];
192:   h    = tmp_x[1];

194:   VecSet(f, 1.0);
195:   VecSet(h, 0.0);

197:   KSPCreate(PETSC_COMM_WORLD, &ksp);
198:   KSPSetOperators(ksp, A, A);
199:   KSPGetPC(ksp, &pc);
200:   PCSetType(pc, PCNONE);
201:   KSPSetFromOptions(ksp);

203:   KSPSolve(ksp, b, x);
204:   VecNestGetSubVecs(x,NULL,&tmp_x);
205:   x1   = tmp_x[0];
206:   x2   = tmp_x[1];

208:   PetscPrintf(PETSC_COMM_WORLD, "x1 \n");
209:   PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL);
210:   VecView(x1, PETSC_VIEWER_STDOUT_WORLD);
211:   PetscPrintf(PETSC_COMM_WORLD, "x2 \n");
212:   VecView(x2, PETSC_VIEWER_STDOUT_WORLD);
213:   PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);

215:   KSPDestroy(&ksp);
216:   VecDestroy(&x);
217:   VecDestroy(&b);
218:   MatDestroy(&A11);
219:   MatDestroy(&A12);
220:   MatDestroy(&A21);

222:   MatDestroy(&A);
223:   return 0;
224: }

226: int main(int argc, char **args)
227: {

229:   PetscInitialize(&argc, &args,(char*)0, help);
230:   test_solve();
231:   test_solve_matgetvecs();
232:   PetscFinalize();
233:   return 0;
234: }

236: /*TEST

238:     test:

240:     test:
241:       suffix: 2
242:       nsize: 2

244:     test:
245:       suffix: 3
246:       nsize: 2
247:       args: -ksp_monitor_short -ksp_type bicg
248:       requires: !single

250: TEST*/