Actual source code: ex22.c

petsc-3.13.4 2020-08-01
Report Typos and Errors
  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: }


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

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

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

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

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

159:   VecDestroy(&diag);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

240: /*TEST

242:     test:

244:     test:
245:       suffix: 2
246:       nsize: 2

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

254: TEST*/