Actual source code: samgmgpetsc.c

      1: #define PETSCKSP_DLL

     3:  #include global.h
     4:  #include petscfunc.h
     5:  #include petscksp.h
     6:  #include petscmg.h

      8: static char help[] = "Does PETSc multigrid cycling using hierarchy build by SAMG\n\n";

     10: PetscErrorCode samgmgpetsc(const int numnodes, double* Asky, int* ia,
     11:                 int* ja, double* rhs, double* u_approx,
     12:                 const OPTIONS *options)
     13: {
     14:    /*..Petsc variables..*/
     15:    Vec      x, b;         /* approx solution and RHS */
     16:    Mat      A;            /* linear system matrix */
     17:    KSP     ksp;         /* linear solver context */
     18:    PC       pc;           /* preconditioner context */
     19:    PCType   pctype;       /* preconditioning technique */
     20:    KSP      ksp;          /* KSP context */
     22:    int  its;    /* Error messages and number of iterations */
     23:    /*..Other variables for the PETSc interface..*/
     24:    int      *nnz_per_row; /* integer vector to hold the number of nonzeros */
     25:                           /* of each row. This vector will be used to      */
     26:                           /* allocate memory for the matrix, and to store  */
     27:                           /* elements in the matrix                        */
     28:    int      *cols;        /* cols is a vector of collumn indices used in   */
     29:                           /* assembling the PETSc rhs vector               */
     30:    PetscScalar *sol_array;/* sol_array used to pass the PETSc solution     */
     31:                           /* back to the calling program                   */
     32:    /*..Variables used to customize the convergence criterium to            */
     33:    /*  ||res|| / ||b|| < tol                                               */
     34:    double      bnrm2;
     35:    CONVHIST    *convhist;
     36:    /*..Context for the SAMG preconditioner..*/
     37:    SamgShellPC *samg_ctx;
     38:    /*..Variables to extract SAMG hierarchy..*/
     39:    int         k, levels, numnonzero;
     40:    double      normdiff;
     41:    GridCtx     grid[MAX_LEVELS];
     42:    char        pathfilename[80], basefilename[80];
     43:    /*..Variables for intermediate levels..*/
     44:    KSP         ksp_pre, ksp_post;
     45:    PC          pc_pre, pc_post;
     46:    Mat         FineLevelMatrix;
     47:    int         petsc_level, size;
     48:    /*..Variables for coarse grid solve..*/
     49:    KSP        coarsegridksp;
     50:    PC          coarsegridpc;
     51:    KSP         coarsegridksp;
     52:    int         coarsegrid_n;
     53:    double      coarsegrid_rnorm;
     54:    /*..Variables that determine behaviour of the code..*/
     55:    static PetscErrorCode  debug = *(options->DEBUG);
     56:    /*..Other variables..*/
     57:    int         I;
     58:    PetscTruth  flg, issamg, issamg_print;
     59:    /*..Variables for CPU timings..*/
     60:    PetscLogDouble  v1,v2,t_setup, t_solve;

     62:    /*..Executable statements..*/
     63:    PetscInitialize( (int*) 0, (char ***) 0,(char *) 0, help);

     65:    /*..Get start time of linear system setup..*/
     66:    PetscGetTime(&v1);

     68:    PetscMalloc(numnodes * sizeof(int),&nnz_per_row);

     70:    /*..The numbero f nonzeros entries in row I can be calculated as
     71:        ia[I+1] - 1 - ia[I] + 1 = ia[I+1] - ia[I]                         ..*/
     72:    for (I=0;I<numnodes;I++)
     73:        nnz_per_row[I] = ia[I+1] - ia[I];

     75:    /*..Allocate (create) SeqAIJ matrix  for use within PETSc..*/
     76:    MatCreate(PETSC_COMM_WORLD,&A);
     77:    MatSetSizes(A,numnodes,numnodes,numnodes,numnodes);
     78:    MatSetType(A,MATSEQAIJ);
     79:    MatSeqAIJSetPreallocation(A,0,nnz_per_row);

     81:   /*..Assemble matrix  for use within PETSc..*/
     82:    for (I=0;I<numnodes;I++){
     83:       MatSetValues(A,
     84:                1,              /* number of rows */
     85:                &I,             /* pointer to global row number */
     86:                nnz_per_row[I], /* number of collums = number of nonzero ... */
     87:                                /* entries in row I                          */
     88:                &(ja[ ia[I] ]),
     89:                               /* vector global column indices */
     90:                (PetscScalar *) &(Asky[ ia[I] ]),
     91:                               /* vector of coefficients */
     92:         INSERT_VALUES);
     93:    }

     95:    MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
     96:    MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
     97:    if (debug)
     98:        MatView(A,PETSC_VIEWER_STDOUT_SELF);

    100:    /*..Create solution and rhs vector.  Note that we form  vector from
    101:      scratch and then duplicate as needed..*/
    102:    VecCreate(PETSC_COMM_WORLD,&x);
    103:    VecSetSizes(x,PETSC_DECIDE,numnodes);
    104:    VecSetType(x,VECSEQ);
    105:    VecDuplicate(x,&b);

    107:    PetscMalloc(numnodes * sizeof(int),&cols);
    108:    for (I=0;I<numnodes;I++)
    109:        cols[I] = I;

    111:    /*..Assemble the right-hand side vector for use within PETSc..*/
    112:    VecSetValues(b,numnodes,cols,(PetscScalar*)rhs,INSERT_VALUES);
    113:
    114:    VecAssemblyBegin(b);
    115:    VecAssemblyEnd(b);
    116:    if (debug){
    117:       printf("[PETSc]:The right-hand side \n");
    118:       VecView(b,PETSC_VIEWER_STDOUT_SELF);
    119:       printf("\n");
    120:    }
    121:    VecNorm(b,NORM_2,&bnrm2);

    123:    /*..Assemble the start solution vector for use within PETSc..*/
    124:    VecSetValues(x,numnodes,cols,(PetscScalar*)u_approx,INSERT_VALUES);
    125:
    126:    VecAssemblyBegin(x);
    127:    VecAssemblyEnd(x);

    129:    VecNorm(b,NORM_2,&bnrm2);
    130:    if (debug)
    131:        printf("[PETSc]:The right-hand side norm = %e \n",bnrm2);

    133:    /*..Create linear solver context..*/
    134:    KSPCreate(MPI_COMM_WORLD,&ksp);

    136:    /*..Set operators. Here the matrix that defines the linear system
    137:      also serves as the preconditioning matrix..*/
    138:    KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);

    140:    /*..Extract pc type from context..*/
    141:    KSPGetPC(ksp,&pc);

    143:    /*..Customize tolerances..*/
    144:    KSPSetTolerances(ksp,1e-12,1e-14,PETSC_DEFAULT,
    145:           PETSC_DEFAULT);

    147:     /*..Create user defined context for the shell preconditioner..*/
    148:    SamgShellPCCreate(&samg_ctx);

    150:    /*..Do the setup for the SAMG precondioner..*/
    151:    SamgShellPCSetUp(samg_ctx,A);

    153:    /*..Give user defined preconditioner a name..*/
    154:    PCShellSetName(pc,"SAMG (Scalar mode)");
    155:

    157:    /*..Parse SAMG hierarchy to PETSc variables..*/
    158:    levels = samg_ctx->LEVELS;
    159:    numnonzero = ia[numnodes];
    160:    SamgGetGrid(levels, numnodes, numnonzero, grid, PETSC_NULL);

    162:    /*..Print coarser grid and interpolation operators to file..*/
    163:    PetscOptionsHasName(PETSC_NULL,"-samg_print",&issamg_print);
    164:
    165:    if (issamg_print){
    166:      for (k=2;k<=levels;k++){
    167:          sprintf(pathfilename,"./");
    168:          sprintf(basefilename,"Pcoarsemat.%02u",k);
    169:          PrintMatrix(grid[k].A, pathfilename, basefilename);
    170:      }
    171:      for (k=1;k<=levels-1;k++){
    172:          sprintf(basefilename,"Pinterpol.%02u%02",k, k-1);
    173:          PrintMatrix(grid[k].Interp, pathfilename, basefilename);
    174:      }
    175:    }
    176:
    177:    /*..Perform check on parsing..*/
    178:    PetscOptionsHasName(PETSC_NULL,"-samg_check",&issamg_print);
    179:
    180:    if (issamg_print)
    181:      SamgCheckGalerkin(levels, A, grid, PETSC_NULL);
    182:
    183:    /*..Set KSP solver type..*/
    184:    KSPSetType(ksp,KSPRICHARDSON);
    185:    KSPSetMonitor(ksp,KSPDefaultMonitor,PETSC_NULL, PETSC_NULL);
    186:
    187:    /*..Set MG preconditioner..*/
    188:    PCSetType(pc,PCMG);
    189:    PCMGSetLevels(pc,levels, PETSC_NULL);
    190:    PCMGSetType(pc, PC_MG_MULTIPLICATIVE);
    191:    PCMGSetCycles(pc, 1);
    192:    PCMGSetNumberSmoothUp(pc,1);
    193:    PCMGSetNumberSmoothDown(pc,1);

    195:    /*....Set smoother, work vectors and residual calculation on each
    196:          level....*/
    197:    for (k=1;k<=levels;k++){
    198:        petsc_level = levels - k;
    199:        /*....Get pre-smoothing KSP context....*/
    200:        PCMGGetSmootherDown(pc,petsc_level,&grid[k].ksp_pre);
    201:
    202:        PCMGGetSmootherUp(pc,petsc_level,&grid[k].ksp_post);
    203:
    204:        if (k==1)
    205:           FineLevelMatrix = A;
    206:           else
    207:           FineLevelMatrix = grid[k].A;
    208:        MatGetSize(FineLevelMatrix, &size, &size);
    209:        VecCreate(MPI_COMM_WORLD,&grid[k].x);
    210:        VecSetSizes(grid[k].x,PETSC_DECIDE,size);
    211:        VecSetType(grid[k].x,VECSEQ);
    212:        VecDuplicate(grid[k].x,&grid[k].b);
    213:        VecDuplicate(grid[k].x,&grid[k].r);

    215:        /*....set ksp_pre context....*/
    216:        KSPSetOperators(grid[k].ksp_pre, FineLevelMatrix,
    217:                                FineLevelMatrix, DIFFERENT_NONZERO_PATTERN);
    218:
    219:        KSPGetPC(ksp_pre_pre,&pc_pre);
    220:        KSPSetType(grid[k].ksp_pre, KSPRICHARDSON);
    221:        KSPSetTolerances(grid[k].ksp_pre, 1e-12, 1e-50, 1e7,1);
    222:
    223:        PCSetType(pc_pre, PCSOR);
    224:        PCSORSetSymmetric(pc_pre,SOR_FORWARD_SWEEP);

    226:        /*....set ksp_post context....*/
    227:        KSPSetOperators(grid[k].ksp_post, FineLevelMatrix,
    228:                                FineLevelMatrix, DIFFERENT_NONZERO_PATTERN);
    229:
    230:        KSPGetPC(grid[k].ksp_post,&pc_post);
    231:        KSPSetInitialGuessNonzero(grid[k].ksp_post, PETSC_TRUE);
    232:        KSPSetType(grid[k].ksp_post, KSPRICHARDSON);
    233:        KSPSetTolerances(grid[k].ksp_post, 1e-12, 1e-50, 1e7,1);
    234:
    235:        PCSetType(pc_post, PCSOR);
    236:        PCSORSetSymmetric(pc_post,SOR_BACKWARD_SWEEP);

    238:        PCMGSetX(pc,petsc_level,grid[k].x);
    239:        PCMGSetRhs(pc,petsc_level,grid[k].b);
    240:        PCMGSetR(pc,petsc_level,grid[k].r);
    241:        PCMGSetResidual(pc,petsc_level,PCMGDefaultResidual,FineLevelMatrix);
    242:
    243:    }

    245:    /*....Create interpolation between the levels....*/
    246:    for (k=1;k<=levels-1;k++){
    247:      petsc_level = levels - k;
    248:      PCMGSetInterpolate(pc,petsc_level,grid[k].Interp);
    249:      PCMGSetRestriction(pc,petsc_level,grid[k].Interp);
    250:    }

    252:    /*....Set coarse grid solver....*/
    253:    PCMGGetCoarseSolve(pc,&coarsegridksp);
    254:    KSPSetFromOptions(coarsegridksp);
    255:    KSPSetOperators(coarsegridksp, grid[levels].A, grid[levels].A,
    256:                            DIFFERENT_NONZERO_PATTERN);
    257:    KSPGetPC(coarsegridksp,&coarsegridpc);
    258:    KSPSetType(coarsegridksp, KSPPREONLY);
    259:    PCSetType(coarsegridpc, PCLU);

    261:    /*..Allow above criterea to be overwritten..*/
    262:    KSPSetFromOptions(ksp);

    264:    /*..Indicate that we are going to use a non-zero initial solution..*/
    265:    KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);

    267:    /*..Get end time of linear system setup..*/
    268:    PetscGetTime(&v2);
    269:    t_setup = v2 - v1;

    271:    /*..Get start time of linear solve..*/
    272:    PetscGetTime(&v1);

    274:    /*..Solve linear system..*/
    275:    KSPSolve(ksp,b,x);
    276:    KSPGetIterationNumber(ksp,&its);
    277:
    278:    /*..Print number of iterations..*/
    279:    PetscPrintf(PETSC_COMM_WORLD,"\n** Number of iterations done = %d \n",
    280:                its);

    282:    /*..Get end time of linear solve..*/
    283:    PetscGetTime(&v2);
    284:    t_solve = v2 - v1;

    286:    printf("\n[PETSc]:Time spend in setup = %e \n",t_setup);
    287:    printf("[PETSc]:Time spend in solve = %e \n",t_solve);
    288:    printf("[PETSc]:Total time = %e \n\n", t_setup + t_solve);
    289:
    290:    /*..Copy PETSc solution back..*/
    291:    ierr= VecGetArray(x, &sol_array);
    292:    for (I=0;I<numnodes;I++){
    293:      u_approx[I] = sol_array[I];
    294:    }
    295:    VecRestoreArray(x,&sol_array);

    297:    if (debug){
    298:       printf("[PETSc]:The solution \n");
    299:       VecView(x,PETSC_VIEWER_STDOUT_SELF);
    300:       printf("\n");
    301:    }

    303:    /*..Free work space..*/
    304:    SamgShellPCDestroy(samg_ctx);
    305:    VecDestroy(x);
    306:    VecDestroy(b);
    307:    MatDestroy(A);
    308:    KSPDestroy(ksp);
    309:    for (k=2;k<=levels;k++){
    310:      MatDestroy(grid[k].A);
    311:    }
    312:    for (k=1;k<=levels-1;k++){
    313:      MatDestroy(grid[k].Interp);
    314:    }
    315:    for (k=1;k<=levels-1;k++){
    316:      VecDestroy(grid[k].b);
    317:      VecDestroy(grid[k].x);
    318:      VecDestroy(grid[k].r);
    319:    }

    321:    PetscFinalize();
    322:
    323:    return 0;
    324: }