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;

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

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

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

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

 37:   VecDestroy(&diag);

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

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

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

 58:   A22 = NULL;

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

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

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

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

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

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

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

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

 99:   KSPSolve(ksp, b, x);

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

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

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

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

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

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

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

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

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

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

158:   VecDestroy(&diag);

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

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

176:   /* A21 */
177:   MatTranspose(A12, MAT_INITIAL_MATRIX, &A21);

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

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

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

196:   VecSet(f, 1.0);
197:   VecSet(h, 0.0);

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

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

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

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

224:   MatDestroy(&A);
225:   return(0);
226: }

228: int main(int argc, char **args)
229: {

232:   PetscInitialize(&argc, &args,(char*)0, help);if (ierr) return ierr;
233:   test_solve();
234:   test_solve_matgetvecs();
235:   PetscFinalize();
236:   return ierr;
237: }

239: /*TEST

241:     test:

243:     test:
244:       suffix: 2
245:       nsize: 2

247:     test:
248:       suffix: 3
249:       nsize: 2
250:       args: -ksp_monitor_short -ksp_type bicg
251:       requires: !single

253: TEST*/