Actual source code: ex22.c
petsc-3.4.2 2013-07-02

      1: static const char help[] = "Test MatNest solving a linear system\n\n";

      3: #include <petscksp.h>

      7: PetscErrorCode test_solve(void)
      8: {
      9:   Mat            A11, A12,A21,A22, A, tmp[2][2];
     10:   KSP            ksp;
     11:   PC             pc;
     12:   Vec            b,x, f,h, diag, x1,x2;
     13:   Vec            tmp_x[2],*_tmp_x;
     14:   int            n, np, i,j;

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

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

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

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

     38:   VecDestroy(&diag);

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

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

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

     59:   A22 = NULL;

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

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

     72:   /* Create vectors */
     73:   MatGetVecs(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, SAME_NONZERO_PATTERN);
     89:   KSPSetType(ksp, "gmres");
     90:   KSPGetPC(ksp, &pc);
     91:   PCSetType(pc, "none");
     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:   PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL);
    103:   VecView(x1, PETSC_VIEWER_STDOUT_WORLD);
    104:   PetscPrintf(PETSC_COMM_WORLD, "x2 \n");
    105:   PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL);
    106:   VecView(x2, 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: }

    124: PetscErrorCode test_solve_matgetvecs(void)
    125: {
    126:   Mat            A11, A12,A21, A;
    127:   KSP            ksp;
    128:   PC             pc;
    129:   Vec            b,x, f,h, diag, x1,x2;
    130:   int            n, np, i,j;
    131:   Mat            tmp[2][2];
    132:   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, (double)(i+j*n), INSERT_VALUES);
    168:     }
    169:   }
    170:   MatSetValue(A12, 2,1, (double)(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:   MatGetVecs(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, SAME_NONZERO_PATTERN);
    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:   PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL);
    210:   VecView(x1, PETSC_VIEWER_STDOUT_WORLD);
    211:   PetscPrintf(PETSC_COMM_WORLD, "x2 \n");
    212:   PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL);
    213:   VecView(x2, 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: }

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

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