Actual source code: ex52.c

      2: static char help[] = "Solves a linear system in parallel with KSP. Modified from ex2.c \n\
      3:                       Illustrate how to use external packages MUMPS, SUPERLU and STRUMPACK \n\
      4: Input parameters include:\n\
      5:   -random_exact_sol : use a random exact solution vector\n\
      6:   -view_exact_sol   : write exact solution vector to stdout\n\
      7:   -m <mesh_x>       : number of mesh points in x-direction\n\
      8:   -n <mesh_y>       : number of mesh points in y-direction\n\n";

     10: #include <petscksp.h>

     12: int main(int argc,char **args)
     13: {
     14:   Vec            x,b,u;    /* approx solution, RHS, exact solution */
     15:   Mat            A,F;
     16:   KSP            ksp;      /* linear solver context */
     17:   PC             pc;
     18:   PetscRandom    rctx;     /* random number generator context */
     19:   PetscReal      norm;     /* norm of solution error */
     20:   PetscInt       i,j,Ii,J,Istart,Iend,m = 8,n = 7,its;
     22:   PetscBool      flg=PETSC_FALSE,flg_ilu=PETSC_FALSE,flg_ch=PETSC_FALSE;
     23: #if defined(PETSC_HAVE_MUMPS)
     24:   PetscBool      flg_mumps=PETSC_FALSE,flg_mumps_ch=PETSC_FALSE;
     25: #endif
     26: #if defined(PETSC_HAVE_SUPERLU) || defined(PETSC_HAVE_SUPERLU_DIST)
     27:   PetscBool      flg_superlu=PETSC_FALSE;
     28: #endif
     29: #if defined(PETSC_HAVE_STRUMPACK)
     30:   PetscBool      flg_strumpack=PETSC_FALSE;
     31: #endif
     32:   PetscScalar    v;
     33:   PetscMPIInt    rank,size;
     34: #if defined(PETSC_USE_LOG)
     35:   PetscLogStage  stage;
     36: #endif

     38:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
     39:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
     40:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
     41:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
     42:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
     43:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     44:          Compute the matrix and right-hand-side vector that define
     45:          the linear system, Ax = b.
     46:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     47:   MatCreate(PETSC_COMM_WORLD,&A);
     48:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
     49:   MatSetFromOptions(A);
     50:   MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);
     51:   MatSeqAIJSetPreallocation(A,5,NULL);
     52:   MatSetUp(A);

     54:   /*
     55:      Currently, all PETSc parallel matrix formats are partitioned by
     56:      contiguous chunks of rows across the processors.  Determine which
     57:      rows of the matrix are locally owned.
     58:   */
     59:   MatGetOwnershipRange(A,&Istart,&Iend);

     61:   /*
     62:      Set matrix elements for the 2-D, five-point stencil in parallel.
     63:       - Each processor needs to insert only elements that it owns
     64:         locally (but any non-local elements will be sent to the
     65:         appropriate processor during matrix assembly).
     66:       - Always specify global rows and columns of matrix entries.

     68:      Note: this uses the less common natural ordering that orders first
     69:      all the unknowns for x = h then for x = 2h etc; Hence you see J = Ii +- n
     70:      instead of J = I +- m as you might expect. The more standard ordering
     71:      would first do all variables for y = h, then y = 2h etc.

     73:    */
     74:   PetscLogStageRegister("Assembly", &stage);
     75:   PetscLogStagePush(stage);
     76:   for (Ii=Istart; Ii<Iend; Ii++) {
     77:     v = -1.0; i = Ii/n; j = Ii - i*n;
     78:     if (i>0)   {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
     79:     if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
     80:     if (j>0)   {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
     81:     if (j<n-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
     82:     v = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);
     83:   }

     85:   /*
     86:      Assemble matrix, using the 2-step process:
     87:        MatAssemblyBegin(), MatAssemblyEnd()
     88:      Computations can be done while messages are in transition
     89:      by placing code between these two statements.
     90:   */
     91:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
     92:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
     93:   PetscLogStagePop();

     95:   /* A is symmetric. Set symmetric flag to enable ICC/Cholesky preconditioner */
     96:   MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);

     98:   /*
     99:      Create parallel vectors.
    100:       - We form 1 vector from scratch and then duplicate as needed.
    101:       - When using VecCreate(), VecSetSizes and VecSetFromOptions()
    102:         in this example, we specify only the
    103:         vector's global dimension; the parallel partitioning is determined
    104:         at runtime.
    105:       - When solving a linear system, the vectors and matrices MUST
    106:         be partitioned accordingly.  PETSc automatically generates
    107:         appropriately partitioned matrices and vectors when MatCreate()
    108:         and VecCreate() are used with the same communicator.
    109:       - The user can alternatively specify the local vector and matrix
    110:         dimensions when more sophisticated partitioning is needed
    111:         (replacing the PETSC_DECIDE argument in the VecSetSizes() statement
    112:         below).
    113:   */
    114:   VecCreate(PETSC_COMM_WORLD,&u);
    115:   VecSetSizes(u,PETSC_DECIDE,m*n);
    116:   VecSetFromOptions(u);
    117:   VecDuplicate(u,&b);
    118:   VecDuplicate(b,&x);

    120:   /*
    121:      Set exact solution; then compute right-hand-side vector.
    122:      By default we use an exact solution of a vector with all
    123:      elements of 1.0;  Alternatively, using the runtime option
    124:      -random_sol forms a solution vector with random components.
    125:   */
    126:   PetscOptionsGetBool(NULL,NULL,"-random_exact_sol",&flg,NULL);
    127:   if (flg) {
    128:     PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
    129:     PetscRandomSetFromOptions(rctx);
    130:     VecSetRandom(u,rctx);
    131:     PetscRandomDestroy(&rctx);
    132:   } else {
    133:     VecSet(u,1.0);
    134:   }
    135:   MatMult(A,u,b);

    137:   /*
    138:      View the exact solution vector if desired
    139:   */
    140:   flg  = PETSC_FALSE;
    141:   PetscOptionsGetBool(NULL,NULL,"-view_exact_sol",&flg,NULL);
    142:   if (flg) {VecView(u,PETSC_VIEWER_STDOUT_WORLD);}

    144:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    145:                 Create the linear solver and set various options
    146:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    148:   /*
    149:      Create linear solver context
    150:   */
    151:   KSPCreate(PETSC_COMM_WORLD,&ksp);
    152:   KSPSetOperators(ksp,A,A);

    154:   /*
    155:     Example of how to use external package MUMPS
    156:     Note: runtime options
    157:           '-ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps -mat_mumps_icntl_7 2 -mat_mumps_icntl_1 0.0'
    158:           are equivalent to these procedural calls
    159:   */
    160: #if defined(PETSC_HAVE_MUMPS)
    161:   flg_mumps    = PETSC_FALSE;
    162:   flg_mumps_ch = PETSC_FALSE;
    163:   PetscOptionsGetBool(NULL,NULL,"-use_mumps_lu",&flg_mumps,NULL);
    164:   PetscOptionsGetBool(NULL,NULL,"-use_mumps_ch",&flg_mumps_ch,NULL);
    165:   if (flg_mumps || flg_mumps_ch) {
    166:     KSPSetType(ksp,KSPPREONLY);
    167:     PetscInt  ival,icntl;
    168:     PetscReal val;
    169:     KSPGetPC(ksp,&pc);
    170:     if (flg_mumps) {
    171:       PCSetType(pc,PCLU);
    172:     } else if (flg_mumps_ch) {
    173:       MatSetOption(A,MAT_SPD,PETSC_TRUE); /* set MUMPS id%SYM=1 */
    174:       PCSetType(pc,PCCHOLESKY);
    175:     }
    176:     PCFactorSetMatSolverType(pc,MATSOLVERMUMPS);
    177:     PCFactorSetUpMatSolverType(pc); /* call MatGetFactor() to create F */
    178:     PCFactorGetMatrix(pc,&F);

    180:     /* sequential ordering */
    181:     icntl = 7; ival = 2;
    182:     MatMumpsSetIcntl(F,icntl,ival);

    184:     /* threshold for row pivot detection */
    185:     MatMumpsSetIcntl(F,24,1);
    186:     icntl = 3; val = 1.e-6;
    187:     MatMumpsSetCntl(F,icntl,val);

    189:     /* compute determinant of A */
    190:     MatMumpsSetIcntl(F,33,1);
    191:   }
    192: #endif

    194:   /*
    195:     Example of how to use external package SuperLU
    196:     Note: runtime options
    197:           '-ksp_type preonly -pc_type ilu -pc_factor_mat_solver_type superlu -mat_superlu_ilu_droptol 1.e-8'
    198:           are equivalent to these procedual calls
    199:   */
    200: #if defined(PETSC_HAVE_SUPERLU) || defined(PETSC_HAVE_SUPERLU_DIST)
    201:   flg_ilu     = PETSC_FALSE;
    202:   flg_superlu = PETSC_FALSE;
    203:   PetscOptionsGetBool(NULL,NULL,"-use_superlu_lu",&flg_superlu,NULL);
    204:   PetscOptionsGetBool(NULL,NULL,"-use_superlu_ilu",&flg_ilu,NULL);
    205:   if (flg_superlu || flg_ilu) {
    206:     KSPSetType(ksp,KSPPREONLY);
    207:     KSPGetPC(ksp,&pc);
    208:     if (flg_superlu) {
    209:       PCSetType(pc,PCLU);
    210:     } else if (flg_ilu) {
    211:       PCSetType(pc,PCILU);
    212:     }
    213:     if (size == 1) {
    214: #if !defined(PETSC_HAVE_SUPERLU)
    215:       SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This test requires SUPERLU");
    216: #else
    217:       PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU);
    218: #endif
    219:     } else {
    220: #if !defined(PETSC_HAVE_SUPERLU_DIST)
    221:       SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This test requires SUPERLU_DIST");
    222: #else
    223:       PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU_DIST);
    224: #endif
    225:     }
    226:     PCFactorSetUpMatSolverType(pc); /* call MatGetFactor() to create F */
    227:     PCFactorGetMatrix(pc,&F);
    228: #if defined(PETSC_HAVE_SUPERLU)
    229:     if (size == 1) {
    230:       MatSuperluSetILUDropTol(F,1.e-8);
    231:     }
    232: #endif
    233:   }
    234: #endif

    237:   /*
    238:     Example of how to use external package STRUMPACK
    239:     Note: runtime options
    240:           '-pc_type lu/ilu \
    241:            -pc_factor_mat_solver_type strumpack \
    242:            -mat_strumpack_reordering METIS \
    243:            -mat_strumpack_colperm 0 \
    244:            -mat_strumpack_hss_rel_tol 1.e-3 \
    245:            -mat_strumpack_hss_min_sep_size 50 \
    246:            -mat_strumpack_max_rank 100 \
    247:            -mat_strumpack_leaf_size 4'
    248:        are equivalent to these procedural calls

    250:     We refer to the STRUMPACK-sparse manual, section 5, for more info on
    251:     how to tune the preconditioner.
    252:   */
    253: #if defined(PETSC_HAVE_STRUMPACK)
    254:   flg_ilu       = PETSC_FALSE;
    255:   flg_strumpack = PETSC_FALSE;
    256:   PetscOptionsGetBool(NULL,NULL,"-use_strumpack_lu",&flg_strumpack,NULL);
    257:   PetscOptionsGetBool(NULL,NULL,"-use_strumpack_ilu",&flg_ilu,NULL);
    258:   if (flg_strumpack || flg_ilu) {
    259:     KSPSetType(ksp,KSPPREONLY);
    260:     KSPGetPC(ksp,&pc);
    261:     if (flg_strumpack) {
    262:       PCSetType(pc,PCLU);
    263:     } else if (flg_ilu) {
    264:       PCSetType(pc,PCILU);
    265:     }
    266: #if !defined(PETSC_HAVE_STRUMPACK)
    267:     SETERRQ(PETSC_COMM_WORLD,PETSC_,"This test requires STRUMPACK");
    268: #endif
    269:     PCFactorSetMatSolverType(pc,MATSOLVERSTRUMPACK);
    270:     PCFactorSetUpMatSolverType(pc); /* call MatGetFactor() to create F */
    271:     PCFactorGetMatrix(pc,&F);
    272: #if defined(PETSC_HAVE_STRUMPACK)
    273:     /* Set the fill-reducing reordering.                              */
    274:     MatSTRUMPACKSetReordering(F,MAT_STRUMPACK_METIS);
    275:     /* Since this is a simple discretization, the diagonal is always  */
    276:     /* nonzero, and there is no need for the extra MC64 permutation.  */
    277:     MatSTRUMPACKSetColPerm(F,PETSC_FALSE);
    278:     /* The compression tolerance used when doing low-rank compression */
    279:     /* in the preconditioner. This is problem specific!               */
    280:     MatSTRUMPACKSetHSSRelTol(F,1.e-3);
    281:     /* Set minimum matrix size for HSS compression to 15 in order to  */
    282:     /* demonstrate preconditioner on small problems. For performance  */
    283:     /* a value of say 500 is better.                                  */
    284:     MatSTRUMPACKSetHSSMinSepSize(F,15);
    285:     /* You can further limit the fill in the preconditioner by        */
    286:     /* setting a maximum rank                                         */
    287:     MatSTRUMPACKSetHSSMaxRank(F,100);
    288:     /* Set the size of the diagonal blocks (the leafs) in the HSS     */
    289:     /* approximation. The default value should be better for real     */
    290:     /* problems. This is mostly for illustration on a small problem.  */
    291:     MatSTRUMPACKSetHSSLeafSize(F,4);
    292: #endif
    293:   }
    294: #endif

    296:   /*
    297:     Example of how to use procedural calls that are equivalent to
    298:           '-ksp_type preonly -pc_type lu/ilu -pc_factor_mat_solver_type petsc'
    299:   */
    300:   flg     = PETSC_FALSE;
    301:   flg_ilu = PETSC_FALSE;
    302:   flg_ch  = PETSC_FALSE;
    303:   PetscOptionsGetBool(NULL,NULL,"-use_petsc_lu",&flg,NULL);
    304:   PetscOptionsGetBool(NULL,NULL,"-use_petsc_ilu",&flg_ilu,NULL);
    305:   PetscOptionsGetBool(NULL,NULL,"-use_petsc_ch",&flg_ch,NULL);
    306:   if (flg || flg_ilu || flg_ch) {
    307:     Vec diag;

    309:     KSPSetType(ksp,KSPPREONLY);
    310:     KSPGetPC(ksp,&pc);
    311:     if (flg) {
    312:       PCSetType(pc,PCLU);
    313:     } else if (flg_ilu) {
    314:       PCSetType(pc,PCILU);
    315:     } else if (flg_ch) {
    316:       PCSetType(pc,PCCHOLESKY);
    317:     }
    318:     PCFactorSetMatSolverType(pc,MATSOLVERPETSC);
    319:     PCFactorSetUpMatSolverType(pc); /* call MatGetFactor() to create F */
    320:     PCFactorGetMatrix(pc,&F);

    322:     /* Test MatGetDiagonal() */
    323:     KSPSetUp(ksp);
    324:     VecDuplicate(x,&diag);
    325:     MatGetDiagonal(F,diag);
    326:     /* VecView(diag,PETSC_VIEWER_STDOUT_WORLD); */
    327:     VecDestroy(&diag);
    328:   }

    330:   KSPSetFromOptions(ksp);

    332:   /* Get info from matrix factors */
    333:   KSPSetUp(ksp);

    335: #if defined(PETSC_HAVE_MUMPS)
    336:   if (flg_mumps || flg_mumps_ch) {
    337:     PetscInt  icntl,infog34;
    338:     PetscReal cntl,rinfo12,rinfo13;
    339:     icntl = 3;
    340:     MatMumpsGetCntl(F,icntl,&cntl);

    342:     /* compute determinant */
    343:     if (!rank) {
    344:       MatMumpsGetInfog(F,34,&infog34);
    345:       MatMumpsGetRinfog(F,12,&rinfo12);
    346:       MatMumpsGetRinfog(F,13,&rinfo13);
    347:       PetscPrintf(PETSC_COMM_SELF,"  Mumps row pivot threshold = %g\n",cntl);
    348:       PetscPrintf(PETSC_COMM_SELF,"  Mumps determinant = (%g, %g) * 2^%D \n",(double)rinfo12,(double)rinfo13,infog34);
    349:     }
    350:   }
    351: #endif

    353:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    354:                       Solve the linear system
    355:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    356:   KSPSolve(ksp,b,x);

    358:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    359:                       Check solution and clean up
    360:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    361:   VecAXPY(x,-1.0,u);
    362:   VecNorm(x,NORM_2,&norm);
    363:   KSPGetIterationNumber(ksp,&its);

    365:   /*
    366:      Print convergence information.  PetscPrintf() produces a single
    367:      print statement from all processes that share a communicator.
    368:      An alternative is PetscFPrintf(), which prints to a file.
    369:   */
    370:   if (norm < 1.e-12) {
    371:     PetscPrintf(PETSC_COMM_WORLD,"Norm of error < 1.e-12 iterations %D\n",its);
    372:   } else {
    373:     PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);
    374:  }

    376:   /*
    377:      Free work space.  All PETSc objects should be destroyed when they
    378:      are no longer needed.
    379:   */
    380:   KSPDestroy(&ksp);
    381:   VecDestroy(&u);  VecDestroy(&x);
    382:   VecDestroy(&b);  MatDestroy(&A);

    384:   /*
    385:      Always call PetscFinalize() before exiting a program.  This routine
    386:        - finalizes the PETSc libraries as well as MPI
    387:        - provides summary and diagnostic information if certain runtime
    388:          options are chosen (e.g., -log_view).
    389:   */
    390:   PetscFinalize();
    391:   return ierr;
    392: }

    395: /*TEST

    397:    test:
    398:       args: -use_petsc_lu
    399:       output_file: output/ex52_2.out

    401:    test:
    402:       suffix: mumps
    403:       nsize: 3
    404:       requires: mumps
    405:       args: -use_mumps_lu
    406:       output_file: output/ex52_1.out

    408:    test:
    409:       suffix: mumps_2
    410:       nsize: 3
    411:       requires: mumps
    412:       args: -use_mumps_ch
    413:       output_file: output/ex52_1.out

    415:    test:
    416:       suffix: mumps_3
    417:       nsize: 3
    418:       requires: mumps
    419:       args: -use_mumps_ch -mat_type sbaij
    420:       output_file: output/ex52_1.out

    422:    test:
    423:       suffix: mumps_omp_2
    424:       nsize: 4
    425:       requires: mumps hwloc openmp pthread define(PETSC_HAVE_MPI_PROCESS_SHARED_MEMORY)
    426:       args: -use_mumps_lu -mat_mumps_use_omp_threads 2
    427:       output_file: output/ex52_1.out

    429:    test:
    430:       suffix: mumps_omp_3
    431:       nsize: 4
    432:       requires: mumps hwloc openmp pthread define(PETSC_HAVE_MPI_PROCESS_SHARED_MEMORY)
    433:       args: -use_mumps_ch -mat_mumps_use_omp_threads 3
    434:       # Ignore the warning since we are intentionally testing the imbalanced case
    435:       filter: grep -v "Warning: number of OpenMP threads"
    436:       output_file: output/ex52_1.out

    438:    test:
    439:       suffix: mumps_omp_4
    440:       nsize: 4
    441:       requires: mumps hwloc openmp pthread define(PETSC_HAVE_MPI_PROCESS_SHARED_MEMORY)
    442:       # let petsc guess a proper number for threads
    443:       args: -use_mumps_ch -mat_type sbaij -mat_mumps_use_omp_threads
    444:       output_file: output/ex52_1.out

    446:    test:
    447:       suffix: strumpack
    448:       requires: strumpack
    449:       args: -use_strumpack_lu
    450:       output_file: output/ex52_3.out

    452:    test:
    453:       suffix: strumpack_2
    454:       nsize: 2
    455:       requires: strumpack
    456:       args: -use_strumpack_lu
    457:       output_file: output/ex52_3.out

    459:    test:
    460:       suffix: strumpack_ilu
    461:       requires: strumpack
    462:       args: -use_strumpack_ilu
    463:       output_file: output/ex52_3.out

    465:    test:
    466:       suffix: strumpack_ilu_2
    467:       nsize: 2
    468:       requires: strumpack
    469:       args: -use_strumpack_ilu
    470:       output_file: output/ex52_3.out

    472:    test:
    473:       suffix: superlu
    474:       requires: superlu superlu_dist
    475:       args: -use_superlu_lu
    476:       output_file: output/ex52_2.out

    478:    test:
    479:       suffix: superlu_dist
    480:       nsize: 2
    481:       requires: superlu superlu_dist
    482:       args: -use_superlu_lu
    483:       output_file: output/ex52_2.out

    485:    test:
    486:       suffix: superlu_ilu
    487:       requires: superlu superlu_dist
    488:       args: -use_superlu_ilu
    489:       output_file: output/ex52_2.out

    491: TEST*/