Actual source code: ex22.c

petsc-3.12.2 2019-11-22
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:   int            n, np, i,j;

 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:   /* Create vectors */
 71:   MatCreateVecs(A12, &h, &f);

 73:   VecSet(f, 1.0);
 74:   VecSet(h, 0.0);

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

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

 85:   KSPCreate(PETSC_COMM_WORLD, &ksp);
 86:   KSPSetOperators(ksp, A, A);
 87:   KSPSetType(ksp, "gmres");
 88:   KSPGetPC(ksp, &pc);
 89:   PCSetType(pc, "none");
 90:   KSPSetFromOptions(ksp);

 92:   KSPSolve(ksp, b, x);

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

 96:   x1 = _tmp_x[0];
 97:   x2 = _tmp_x[1];

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

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

115:   MatDestroy(&A);
116:   return(0);
117: }


120: PetscErrorCode test_solve_matgetvecs(void)
121: {
122:   Mat            A11, A12,A21, A;
123:   KSP            ksp;
124:   PC             pc;
125:   Vec            b,x, f,h, diag, x1,x2;
126:   int            n, np, i,j;
127:   Mat            tmp[2][2];
128:   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++) {
163:       MatSetValue(A12, i,j, (PetscScalar)(i+j*n), INSERT_VALUES);
164:     }
165:   }
166:   MatSetValue(A12, 2,1, (PetscScalar)(4), INSERT_VALUES);
167:   MatAssemblyBegin(A12, MAT_FINAL_ASSEMBLY);
168:   MatAssemblyEnd(A12, MAT_FINAL_ASSEMBLY);

170:   /* A21 */
171:   MatTranspose(A12, MAT_INITIAL_MATRIX, &A21);

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

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

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

190:   VecSet(f, 1.0);
191:   VecSet(h, 0.0);

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

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

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

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

218:   MatDestroy(&A);
219:   return(0);
220: }

222: int main(int argc, char **args)
223: {

226:   PetscInitialize(&argc, &args,(char*)0, help);if (ierr) return ierr;
227:   test_solve();
228:   test_solve_matgetvecs();
229:   PetscFinalize();
230:   return ierr;
231: }

233: /*TEST

235:     test:

237:     test:
238:       suffix: 2
239:       nsize: 2

241:     test:
242:       suffix: 3
243:       nsize: 2
244:       args: -ksp_monitor_short -ksp_type bicg
245:       requires: !single

247: TEST*/