Actual source code: ex5.c

      2: static char help[] = "Solves two linear systems in parallel with KSP.  The code\n\
      3: illustrates repeated solution of linear systems with the same preconditioner\n\
      4: method but different matrices (having the same nonzero structure).  The code\n\
      5: also uses multiple profiling stages.  Input arguments are\n\
      6:   -m <size> : problem size\n\
      7:   -mat_nonsym : use nonsymmetric matrix (default is symmetric)\n\n";

      9: /*T
     10:    Concepts: KSP^repeatedly solving linear systems;
     11:    Concepts: PetscLog^profiling multiple stages of code;
     12:    Processors: n
     13: T*/

     17: /*
     18:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
     19:   automatically includes:
     20:      petscsys.h       - base PETSc routines   petscvec.h - vectors
     21:      petscmat.h - matrices
     22:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
     23:      petscviewer.h - viewers               petscpc.h  - preconditioners
     24: */
     25: #include <petscksp.h>

     27: int main(int argc,char **args)
     28: {
     29:   KSP            ksp;              /* linear solver context */
     30:   Mat            C;                /* matrix */
     31:   Vec            x,u,b;            /* approx solution, RHS, exact solution */
     32:   PetscReal      norm;             /* norm of solution error */
     33:   PetscScalar    v,none = -1.0;
     34:   PetscInt       Ii,J,ldim,low,high,iglobal,Istart,Iend;
     36:   PetscInt       i,j,m = 3,n = 2,its;
     37:   PetscMPIInt    size,rank;
     38:   PetscBool      mat_nonsymmetric = PETSC_FALSE;
     39:   PetscBool      testnewC = PETSC_FALSE,testscaledMat=PETSC_FALSE;
     40: #if defined(PETSC_USE_LOG)
     41:   PetscLogStage stages[2];
     42: #endif

     44:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
     45:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
     46:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
     47:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
     48:   n    = 2*size;

     50:   /*
     51:      Set flag if we are doing a nonsymmetric problem; the default is symmetric.
     52:   */
     53:   PetscOptionsGetBool(NULL,NULL,"-mat_nonsym",&mat_nonsymmetric,NULL);

     55:   PetscOptionsGetBool(NULL,NULL,"-test_scaledMat",&testscaledMat,NULL);

     57:   /*
     58:      Register two stages for separate profiling of the two linear solves.
     59:      Use the runtime option -log_view for a printout of performance
     60:      statistics at the program's conlusion.
     61:   */
     62:   PetscLogStageRegister("Original Solve",&stages[0]);
     63:   PetscLogStageRegister("Second Solve",&stages[1]);

     65:   /* -------------- Stage 0: Solve Original System ---------------------- */
     66:   /*
     67:      Indicate to PETSc profiling that we're beginning the first stage
     68:   */
     69:   PetscLogStagePush(stages[0]);

     71:   /*
     72:      Create parallel matrix, specifying only its global dimensions.
     73:      When using MatCreate(), the matrix format can be specified at
     74:      runtime. Also, the parallel partitioning of the matrix is
     75:      determined by PETSc at runtime.
     76:   */
     77:   MatCreate(PETSC_COMM_WORLD,&C);
     78:   MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
     79:   MatSetFromOptions(C);
     80:   MatSetUp(C);

     82:   /*
     83:      Currently, all PETSc parallel matrix formats are partitioned by
     84:      contiguous chunks of rows across the processors.  Determine which
     85:      rows of the matrix are locally owned.
     86:   */
     87:   MatGetOwnershipRange(C,&Istart,&Iend);

     89:   /*
     90:      Set matrix entries matrix in parallel.
     91:       - Each processor needs to insert only elements that it owns
     92:         locally (but any non-local elements will be sent to the
     93:         appropriate processor during matrix assembly).
     94:       - Always specify global row and columns of matrix entries.
     95:   */
     96:   for (Ii=Istart; Ii<Iend; Ii++) {
     97:     v = -1.0; i = Ii/n; j = Ii - i*n;
     98:     if (i>0)   {J = Ii - n; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
     99:     if (i<m-1) {J = Ii + n; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
    100:     if (j>0)   {J = Ii - 1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
    101:     if (j<n-1) {J = Ii + 1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
    102:     v = 4.0; MatSetValues(C,1,&Ii,1,&Ii,&v,ADD_VALUES);
    103:   }

    105:   /*
    106:      Make the matrix nonsymmetric if desired
    107:   */
    108:   if (mat_nonsymmetric) {
    109:     for (Ii=Istart; Ii<Iend; Ii++) {
    110:       v = -1.5; i = Ii/n;
    111:       if (i>1)   {J = Ii-n-1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
    112:     }
    113:   } else {
    114:     MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE);
    115:     MatSetOption(C,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
    116:   }

    118:   /*
    119:      Assemble matrix, using the 2-step process:
    120:        MatAssemblyBegin(), MatAssemblyEnd()
    121:      Computations can be done while messages are in transition
    122:      by placing code between these two statements.
    123:   */
    124:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
    125:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

    127:   /*
    128:      Create parallel vectors.
    129:       - When using VecSetSizes(), we specify only the vector's global
    130:         dimension; the parallel partitioning is determined at runtime.
    131:       - Note: We form 1 vector from scratch and then duplicate as needed.
    132:   */
    133:   VecCreate(PETSC_COMM_WORLD,&u);
    134:   VecSetSizes(u,PETSC_DECIDE,m*n);
    135:   VecSetFromOptions(u);
    136:   VecDuplicate(u,&b);
    137:   VecDuplicate(b,&x);

    139:   /*
    140:      Currently, all parallel PETSc vectors are partitioned by
    141:      contiguous chunks across the processors.  Determine which
    142:      range of entries are locally owned.
    143:   */
    144:   VecGetOwnershipRange(x,&low,&high);

    146:   /*
    147:     Set elements within the exact solution vector in parallel.
    148:      - Each processor needs to insert only elements that it owns
    149:        locally (but any non-local entries will be sent to the
    150:        appropriate processor during vector assembly).
    151:      - Always specify global locations of vector entries.
    152:   */
    153:   VecGetLocalSize(x,&ldim);
    154:   for (i=0; i<ldim; i++) {
    155:     iglobal = i + low;
    156:     v       = (PetscScalar)(i + 100*rank);
    157:     VecSetValues(u,1,&iglobal,&v,INSERT_VALUES);
    158:   }

    160:   /*
    161:      Assemble vector, using the 2-step process:
    162:        VecAssemblyBegin(), VecAssemblyEnd()
    163:      Computations can be done while messages are in transition,
    164:      by placing code between these two statements.
    165:   */
    166:   VecAssemblyBegin(u);
    167:   VecAssemblyEnd(u);

    169:   /*
    170:      Compute right-hand-side vector
    171:   */
    172:   MatMult(C,u,b);

    174:   /*
    175:     Create linear solver context
    176:   */
    177:   KSPCreate(PETSC_COMM_WORLD,&ksp);

    179:   /*
    180:      Set operators. Here the matrix that defines the linear system
    181:      also serves as the preconditioning matrix.
    182:   */
    183:   KSPSetOperators(ksp,C,C);

    185:   /*
    186:      Set runtime options (e.g., -ksp_type <type> -pc_type <type>)
    187:   */
    188:   KSPSetFromOptions(ksp);

    190:   /*
    191:      Solve linear system.  Here we explicitly call KSPSetUp() for more
    192:      detailed performance monitoring of certain preconditioners, such
    193:      as ICC and ILU.  This call is optional, as KSPSetUp() will
    194:      automatically be called within KSPSolve() if it hasn't been
    195:      called already.
    196:   */
    197:   KSPSetUp(ksp);
    198:   KSPSolve(ksp,b,x);

    200:   /*
    201:      Check the error
    202:   */
    203:   VecAXPY(x,none,u);
    204:   VecNorm(x,NORM_2,&norm);
    205:   KSPGetIterationNumber(ksp,&its);
    206:   if (!testscaledMat || norm > 1.e-7) {
    207:     PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);
    208:   }

    210:   /* -------------- Stage 1: Solve Second System ---------------------- */
    211:   /*
    212:      Solve another linear system with the same method.  We reuse the KSP
    213:      context, matrix and vector data structures, and hence save the
    214:      overhead of creating new ones.

    216:      Indicate to PETSc profiling that we're concluding the first
    217:      stage with PetscLogStagePop(), and beginning the second stage with
    218:      PetscLogStagePush().
    219:   */
    220:   PetscLogStagePop();
    221:   PetscLogStagePush(stages[1]);

    223:   /*
    224:      Initialize all matrix entries to zero.  MatZeroEntries() retains the
    225:      nonzero structure of the matrix for sparse formats.
    226:   */
    227:   MatZeroEntries(C);

    229:   /*
    230:      Assemble matrix again.  Note that we retain the same matrix data
    231:      structure and the same nonzero pattern; we just change the values
    232:      of the matrix entries.
    233:   */
    234:   for (i=0; i<m; i++) {
    235:     for (j=2*rank; j<2*rank+2; j++) {
    236:       v = -1.0;  Ii = j + n*i;
    237:       if (i>0)   {J = Ii - n; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
    238:       if (i<m-1) {J = Ii + n; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
    239:       if (j>0)   {J = Ii - 1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
    240:       if (j<n-1) {J = Ii + 1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
    241:       v = 6.0; MatSetValues(C,1,&Ii,1,&Ii,&v,ADD_VALUES);
    242:     }
    243:   }
    244:   if (mat_nonsymmetric) {
    245:     for (Ii=Istart; Ii<Iend; Ii++) {
    246:       v = -1.5; i = Ii/n;
    247:       if (i>1)   {J = Ii-n-1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
    248:     }
    249:   }
    250:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
    251:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

    253:   if (testscaledMat) {
    254:     PetscRandom    rctx;

    256:     /* Scale a(0,0) and a(M-1,M-1) */
    257:     if (!rank) {
    258:       v = 6.0*0.00001; Ii = 0; J = 0;
    259:       MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);
    260:     } else if (rank == size -1) {
    261:       v = 6.0*0.00001; Ii = m*n-1; J = m*n-1;
    262:       MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);
    263:     }
    264:     MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
    265:     MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

    267:     /* Compute a new right-hand-side vector */
    268:     VecDestroy(&u);
    269:     VecCreate(PETSC_COMM_WORLD,&u);
    270:     VecSetSizes(u,PETSC_DECIDE,m*n);
    271:     VecSetFromOptions(u);

    273:     PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
    274:     PetscRandomSetFromOptions(rctx);
    275:     VecSetRandom(u,rctx);
    276:     PetscRandomDestroy(&rctx);
    277:     VecAssemblyBegin(u);
    278:     VecAssemblyEnd(u);
    279:   }

    281:   PetscOptionsGetBool(NULL,NULL,"-test_newMat",&testnewC,NULL);
    282:   if (testnewC) {
    283:     /*
    284:      User may use a new matrix C with same nonzero pattern, e.g.
    285:       ./ex5 -ksp_monitor -mat_type sbaij -pc_type cholesky -pc_factor_mat_solver_type mumps -test_newMat
    286:     */
    287:     Mat Ctmp;
    288:     MatDuplicate(C,MAT_COPY_VALUES,&Ctmp);
    289:     MatDestroy(&C);
    290:     MatDuplicate(Ctmp,MAT_COPY_VALUES,&C);
    291:     MatDestroy(&Ctmp);
    292:   }

    294:   MatMult(C,u,b);

    296:   /*
    297:      Set operators. Here the matrix that defines the linear system
    298:      also serves as the preconditioning matrix.
    299:   */
    300:   KSPSetOperators(ksp,C,C);

    302:   /*
    303:      Solve linear system
    304:   */
    305:   KSPSetUp(ksp);
    306:   KSPSolve(ksp,b,x);

    308:   /*
    309:      Check the error
    310:   */
    311:   VecAXPY(x,none,u);
    312:   VecNorm(x,NORM_2,&norm);
    313:   KSPGetIterationNumber(ksp,&its);
    314:   if (!testscaledMat || norm > 1.e-7) {
    315:     PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);
    316:   }

    318:   /*
    319:      Free work space.  All PETSc objects should be destroyed when they
    320:      are no longer needed.
    321:   */
    322:   KSPDestroy(&ksp);
    323:   VecDestroy(&u);
    324:   VecDestroy(&x);
    325:   VecDestroy(&b);
    326:   MatDestroy(&C);

    328:   /*
    329:      Indicate to PETSc profiling that we're concluding the second stage
    330:   */
    331:   PetscLogStagePop();

    333:   PetscFinalize();
    334:   return ierr;
    335: }

    340: /*TEST

    342:    test:
    343:       args: -pc_type jacobi -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always

    345:    test:
    346:       suffix: 2
    347:       nsize: 2
    348:       args: -pc_type jacobi -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always -ksp_rtol .000001

    350:    test:
    351:       suffix: 5
    352:       nsize: 2
    353:       args: -ksp_gmres_cgs_refinement_type refine_always -ksp_monitor draw::draw_lg -ksp_monitor_true_residual draw::draw_lg

    355:    test:
    356:       suffix: asm
    357:       nsize: 4
    358:       args: -pc_type asm

    360:    test:
    361:       suffix: asm_baij
    362:       nsize: 4
    363:       args: -pc_type asm -mat_type baij
    364:       output_file: output/ex5_asm.out

    366:    test:
    367:       suffix: redundant_0
    368:       args: -m 1000 -pc_type redundant -pc_redundant_number 1 -redundant_ksp_type gmres -redundant_pc_type jacobi

    370:    test:
    371:       suffix: redundant_1
    372:       nsize: 5
    373:       args: -pc_type redundant -pc_redundant_number 1 -redundant_ksp_type gmres -redundant_pc_type jacobi

    375:    test:
    376:       suffix: redundant_2
    377:       nsize: 5
    378:       args: -pc_type redundant -pc_redundant_number 3 -redundant_ksp_type gmres -redundant_pc_type jacobi

    380:    test:
    381:       suffix: redundant_3
    382:       nsize: 5
    383:       args: -pc_type redundant -pc_redundant_number 5 -redundant_ksp_type gmres -redundant_pc_type jacobi

    385:    test:
    386:       suffix: redundant_4
    387:       nsize: 5
    388:       args: -pc_type redundant -pc_redundant_number 3 -redundant_ksp_type gmres -redundant_pc_type jacobi -psubcomm_type interlaced

    390:    test:
    391:       suffix: superlu_dist
    392:       nsize: 15
    393:       requires: superlu_dist
    394:       args: -pc_type lu -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_equil false -m 150 -mat_superlu_dist_r 3 -mat_superlu_dist_c 5 -test_scaledMat

    396:    test:
    397:       suffix: superlu_dist_2
    398:       nsize: 15
    399:       requires: superlu_dist
    400:       args: -pc_type lu -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_equil false -m 150 -mat_superlu_dist_r 3 -mat_superlu_dist_c 5 -test_scaledMat -mat_superlu_dist_fact SamePattern_SameRowPerm
    401:       output_file: output/ex5_superlu_dist.out

    403:    test:
    404:       suffix: superlu_dist_3
    405:       nsize: 15
    406:       requires: superlu_dist
    407:       args: -pc_type lu -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_equil false -m 500 -mat_superlu_dist_r 3 -mat_superlu_dist_c 5 -test_scaledMat -mat_superlu_dist_fact DOFACT
    408:       output_file: output/ex5_superlu_dist.out

    410:    test:
    411:       suffix: superlu_dist_0
    412:       nsize: 1
    413:       requires: superlu_dist
    414:       args: -pc_type lu -pc_factor_mat_solver_type superlu_dist -test_scaledMat
    415:       output_file: output/ex5_superlu_dist.out

    417: TEST*/