Actual source code: samgpetsc.c

      1: #define PETSCKSP_DLL
     2:  #include petscksp.h
     3:  #include src/mat/impls/aij/seq/aij.h
     4:  #include global.h
     5:  #include externc.h
     6:  #include samgfunc.h
     7:  #include petscfunc.h

      9: /*MC
     10:      PCSAMG -   SAMG + PETSc interface
     11:
     12:     This interface allows e.g. to call samg as a shell preconditioner.
     13:     samg is the new algebraic multigrid code by Klaus Stueben [1].
     14:     Reference for SAMG
     15:     [1] K. St\"{u}ben,"Algebraic Multigrid: An Introduction for Positive
     16:         Definite Problems with Applications", in [2], pp. 413--532,
     17:         Academic Press, 2001.
     18:     [2] U. Trottenberg, C. Oosterlee and A. Sch\"{u}ller, "Multigrid
     19:         Methods", Academic Press, 2001.
     20:     [1] is also available as
     21:     [3] K. St\"{u}ben,"Algebraic Multigrid: An Introduction for Positive
     22:         Definite Problems with Applications", Tech. Rep. 53, German
     23:         National Research Center for Information Technology (GMD),
     24:         Schloss Birlinhoven, D-53754 Sankt-Augustin, Germany, March 1999
     25:     For more information on the SAMG-PETSc interface and examples of it's
     26:     use, see
     27:     [4] D. Lahaye "Algebraic Multigrid for Time-Harmonic Magnetic Field
     28:         Computations", PhD Thesis, KU Leuven, Belgium, December 2001.
     29:         (www.cs.kuleuven.ac.be/~domenico)
     30:
     31:    Notes on PETSc part of this interface
     32:    Information on how to set up shell preconditioners in PETSc can be
     33:    in the PETSc documentation under preconditioners.

     35:    This preconditioner has not been completely organized to match the PETSc style,
     36:    see src/ksp/pc/impls/samgpetsc.c for the PC shell routines.

     38:     SAMG is a commercial product available from Klaus Stueben.
     39:
     40:    Level: developer

     42:    Contributed by Domenico Lahaye (domenico.lahaye@cs.kuleuven.ac.be)   January 2001
     43:
     44: M*/

     46: /* ------------------------------------------------------------------- */
     49: /*.. SamgShellPCCreate - This routine creates a user-defined
     50:      preconditioner context.

     52:      Output Parameter:
     53:      shell - user-defined preconditioner context..*/

     55: PetscErrorCode PETSCKSP_DLLEXPORT SamgShellPCCreate(SamgShellPC **shell)
     56: {
     57:    SamgShellPC *newctx;

     60:    PetscNew(SamgShellPC,&newctx);
     61:    *shell = newctx;
     62:    return 0;
     63: }

     65: /* ------------------------------------------------------------------- */
     68: /*..SamgShellPCSetUp - This routine sets up a user-defined
     69:     ramg preconditioner context.

     71:     Input Parameters:
     72:     shell - user-defined preconditioner context
     73:     pmat  - preconditioner matrix

     75:     Output Parameter:
     76:     shell - fully set up user-defined preconditioner context

     78:     This routine calls the setup phase of RAMG..*/
     79:
     80: PetscErrorCode PETSCKSP_DLLEXPORT SamgShellPCSetUp(SamgShellPC *shell, Mat pmat)
     81: {
     83:    int  numnodes, numnonzero, nnz_count;
     84:    int              j, I, J, ncols_getrow, *cols_getrow, *diag;
     85:    PetscScalar      *vals_getrow;
     86:    MatInfo          info;
     87:    Mat_SeqAIJ       *aij = (Mat_SeqAIJ*)pmat->data;
     88:    /*..SAMG variables..*/
     89:    SAMG_PARAM       *samg_param;
     90:    double           *u_approx, *rhs, *Asky;
     91:    int              *ia, *ja;
     92:    /*..Primary SAMG parameters..*/
     93:    /*....System of equations....*/
     94:    int      matrix;
     95:    /*....Start solution and stopping criterium....*/
     96:    int      ifirst;
     97:    double   eps;
     98:    /*....Scalar and coupled system..*/
     99:    int      nsys=1, ndiu=1, ndip=1;
    100:    int      iscale[1], iu[1], ip[1];
    101:    /*....Approach and smoother....*/
    102:    int      nsolve;
    103:    /*....Cycling process....*/
    104:    int      ncyc;
    105:    /*....Repeated calls....*/
    106:    int      iswtch;
    107:    /*....Initial dimensioning..*/
    108:    double   a_cmplx, g_cmplx, p_cmplx, w_avrge;
    109:    /*....Class 1 input parameters controlling input/output..*/
    110:    int      idump, iout;
    111:    double   chktol;
    112:    /*....Output parameters....*/
    113:    double   res_in, res_out;
    114:    int      ncyc_done;
    115:    /*..Numbers of levels created by SAMG..*/
    116:    int      levels;
    117:    /*..Auxilary integer to set secondary parameters..*/
    118:    int      intin;

    120:    /*..Get size and number of unknowns of preconditioner matrix..*/
    121:    MatGetSize(pmat, &numnodes, &numnodes);
    122:    MatGetInfo(pmat,MAT_LOCAL,&info);
    123:    numnonzero = int(info.nz_used);

    125:    /*..Allocate memory for RAMG variables..*/
    126:    PetscMalloc(numnonzero   * sizeof(double),&Asky);
    127:    PetscMalloc((numnodes+1) * sizeof(int),&ia);
    128:    PetscMalloc(numnonzero   * sizeof(int),&ja);
    129:    PetscMalloc(numnodes     * sizeof(double),&u_approx);
    130:    PetscMalloc(numnodes     * sizeof(double),&rhs);

    132:    /*..Store PETSc matrix in compressed skyline format required by SAMG..*/
    133:    nnz_count = 0;
    134:    MatMarkDiagonal_SeqAIJ(pmat);
    135:    diag = aij->diag;

    137:    for (I=0;I<numnodes;I++){
    138:      ia[I]           = nnz_count;

    140:      /*....put in diagonal entry first....*/
    141:      ja[nnz_count]   = I;
    142:      Asky[nnz_count] = aij->a[diag[I]];
    143:      nnz_count++;

    145:      /*....put in off-diagonals....*/
    146:      ncols_getrow = aij->i[I+1] - aij->i[I];
    147:      vals_getrow  = aij->a + aij->i[I];
    148:      cols_getrow  = aij->j + aij->i[I];
    149:      for (j=0;j<ncols_getrow;j++){
    150:        J = cols_getrow[j];
    151:        if (J != I) {
    152:          Asky[nnz_count] = vals_getrow[j];
    153:          ja[nnz_count]   = J;
    154:          nnz_count++;
    155:        }
    156:      }
    157:    }
    158:    ia[numnodes] = nnz_count;

    160:    /*..Allocate memory for SAMG parameters..*/
    161:    PetscNew(SAMG_PARAM,&samg_param);

    163:    /*..Set SAMG parameters..*/
    164:    SamgGetParam(samg_param);

    166:    /*..Set Primary parameters..*/
    167:    matrix  = 12;
    168:    ifirst  = (*samg_param).IFIRST;
    169:    eps     = (*samg_param).EPS;
    170:    nsolve  = (*samg_param).NSOLVE;
    171:    ncyc    = (*samg_param).NCYC;
    172:    iswtch  = (*samg_param).ISWTCH;
    173:    a_cmplx = (*samg_param).A_CMPLX;
    174:    g_cmplx = (*samg_param).G_CMPLX;
    175:    p_cmplx = (*samg_param).P_CMPLX;
    176:    w_avrge = (*samg_param).W_AVRGE;
    177:    chktol  = (*samg_param).CHKTOL;
    178:    idump   = (*samg_param).IDUMP;
    179:    iout    = (*samg_param).IOUT;
    180:    /*..Set secondary parameters..*/
    181:    SAMG_set_levelx(&(samg_param->LEVELX));
    182:    SAMG_set_nptmn(&(samg_param->NPTMN));
    183:    SAMG_set_ecg(&(samg_param->ECG));
    184:    SAMG_set_ewt(&(samg_param->EWT));
    185:    SAMG_set_ncg(&(samg_param->NCG));
    186:    SAMG_set_nwt(&(samg_param->NWT));
    187:    SAMG_set_etr(&(samg_param->ETR));
    188:    SAMG_set_ntr(&(samg_param->NTR));
    189:    SAMG_set_nrd(&(samg_param->NRD));
    190:    SAMG_set_nru(&(samg_param->NRU));
    191:    SAMG_set_nrc(&(samg_param->NRC));
    192:    intin = 6; SAMG_set_logio(&intin);
    193:    intin = 1; SAMG_set_mode_debug(&intin);
    194:    intin = 0; SAMG_set_iter_pre(&intin);
    195:    intin = 0; SAMG_set_np_opt(&intin);
    196:    intin = 0; SAMG_set_ncgrad_default(&intin);

    198:    /*..Reset ncyc such that only setup is performed. This is done by setting
    199:      the last digit of ncyc (the number of cycles performed) equal to zero.
    200:      The first digits of ncyc (related to the solve phase) become
    201:      irrelevant.
    202:    ..*/
    203:    ncyc   = 1000;

    205:    /*..Switch arrays ia and ja to Fortran conventions..*/
    206:    for (j=0;j<=numnodes;j++)
    207:        ia[j]++;
    208:    for (j=0;j<numnonzero;j++)
    209:        ja[j]++;

    211:    PetscPrintf(MPI_COMM_WORLD,"\n\n");
    212:    PetscPrintf(MPI_COMM_WORLD,"******************************************\n");
    213:    PetscPrintf(MPI_COMM_WORLD,"*** Start Setup SAMG code (scal. mode) ***\n");
    214:    PetscPrintf(MPI_COMM_WORLD,"******************************************\n");
    215:    PetscPrintf(MPI_COMM_WORLD,"\n\n");

    217:    /*..Call SAMG..*/
    218:    SAMG(&numnodes, &numnonzero, &nsys,
    219:         ia, ja, Asky, rhs, u_approx, iu, &ndiu, ip, &ndip, &matrix,
    220:         iscale, &res_in, &res_out, &ncyc_done, &ierr, &nsolve,
    221:         &ifirst, &eps, &ncyc, &iswtch, &a_cmplx, &g_cmplx,
    222:          &p_cmplx, &w_avrge, &chktol, &idump, &iout);

    224:    PetscPrintf(MPI_COMM_WORLD,"\n\n");
    225:    PetscPrintf(MPI_COMM_WORLD,"******************************************\n");
    226:    PetscPrintf(MPI_COMM_WORLD,"*** End Setup SAMG code (scal. mode)   ***\n");
    227:    PetscPrintf(MPI_COMM_WORLD,"******************************************\n");
    228:    PetscPrintf(MPI_COMM_WORLD,"\n\n");

    230:    /*..Get number of levels created..*/
    231:    SAMGPETSC_get_levels(&levels);

    233:    /*..Store RAMG output in PETSc context..*/
    234:    shell->A        = Asky;
    235:    shell->IA       = ia;
    236:    shell->JA       = ja;
    237:    shell->PARAM    = samg_param;
    238:    (*shell).LEVELS = levels;

    240:    return 0;
    241: }

    243: /* ------------------------------------------------------------------- */
    246: /*..SamgShellPCApply - This routine applies the AMG code as preconditioner

    248:     Input Parameters:
    249:     ctx - user-defined context, as set by SamgShellSetApply()
    250:     x - input vector

    252:     Output Parameter:
    253:     y - preconditioned vector

    255:    Notes:
    256:    Note that the PCSHELL preconditioner passes a void pointer as the
    257:    first input argument.  This can be cast to be the whatever the user
    258:    has set (via PCSetShellApply()) the application-defined context to be.

    260:    ..*/
    261: /*..To apply AMG as a preconditioner we set:                        */
    262: /*        i) rhs-vector equal to the residual                       */
    263: /*       ii) start solution equal to zero                           */
    264: /*  Implementation notes:                                           */
    265: /*  For the residual (vector r) we take the values from the vector  */
    266: /*  using VecGetArray. No explicit memory allocation for            */
    267: /*  vals_getarray is thus  needed as VecGetArray takes care of it.  */
    268: /*  The allocated memory is freed again using a call to             */
    269: /*  VecRestoreArray. The values in vals_getarray are then copy to   */
    270: /*  rhs of the AMG code using memcpy.                               */

    272: PetscErrorCode PETSCKSP_DLLEXPORT SamgShellPCApply(void *ctx, Vec r, Vec z)
    273: {
    275:    int  I, numnodes, numnonzero, *cols;
    276:    SamgShellPC      *shell = (SamgShellPC *) ctx;
    277:    double           *u_approx, *rhs, *Asky, *vals_getarray;
    278:    int              *ia, *ja;
    279:    SAMG_PARAM       *samg_param;
    280:    /*..Primary SAMG parameters..*/
    281:    /*....System of equations....*/
    282:    int              matrix;
    283:    /*....Start solution and stopping criterium....*/
    284:    int              ifirst;
    285:    double           eps;
    286:    /*....Scalar and coupled system..*/
    287:    int              nsys=1, ndiu=1, ndip=1;
    288:    int              iscale[1], iu[1], ip[1];
    289:    /*....Approach and smoother....*/
    290:    int              nsolve;
    291:    /*....Cycling process....*/
    292:    int              ncyc;
    293:    /*....Repeated calls....*/
    294:    int              iswtch;
    295:    /*....Initial dimensioning..*/
    296:    double           a_cmplx, g_cmplx, p_cmplx, w_avrge;
    297:    /*....Class 1 input parameters controlling input/output..*/
    298:    int              idump, iout;
    299:    double           chktol;
    300:    /*....Output parameters....*/
    301:    double           res_in, res_out;
    302:    int              ncyc_done;

    304:    /*..Get values from context..*/
    305:    Asky       = shell->A;
    306:    ia         = shell->IA;
    307:    ja         = shell->JA;
    308:    samg_param = shell->PARAM;

    310:    /*..Get numnodes and numnonzeros..*/
    311:    /*....numnodes can be determined as the size of the input vector r....*/
    312:    VecGetSize(r,&numnodes);
    313:    /*....numnonzero is determined from the pointer ia....*/
    314:    /*....Remember that ia following Fortran conventions....*/
    315:    numnonzero = ia[numnodes]-1;

    317:    /*..Set the rhs of the call to ramg equal to the residual..*/
    318:    VecGetArray(r,&vals_getarray);

    320:    /*..Allocate memory for rhs and initial solution of call to samg..*/
    321:    PetscMalloc(numnodes     * sizeof(double),&u_approx);
    322:    PetscMalloc(numnodes     * sizeof(double),&rhs);

    324:    /*..Set rhs of call to ramg..*/
    325:    memcpy(rhs, vals_getarray, numnodes * sizeof(*rhs));
    326:
    327:    /*..Set initial solution of call to ramg to zero..*/
    328:    for (I=0;I<numnodes;I++){
    329:        u_approx[I] = 0.;
    330:    }

    332:    /*..Set Primary parameters..*/
    333:    matrix  = (*samg_param).MATRIX;
    334:    ifirst  = (*samg_param).IFIRST;
    335:    eps     = (*samg_param).EPS;
    336:    nsolve  = (*samg_param).NSOLVE;
    337:    ncyc    = (*samg_param).NCYC;
    338:    iswtch  = (*samg_param).ISWTCH;
    339:    a_cmplx = (*samg_param).A_CMPLX;
    340:    g_cmplx = (*samg_param).G_CMPLX;
    341:    p_cmplx = (*samg_param).P_CMPLX;
    342:    w_avrge = (*samg_param).W_AVRGE;
    343:    chktol  = (*samg_param).CHKTOL;
    344:    idump   = (*samg_param).IDUMP;
    345:    iout    = (*samg_param).IOUT;

    347:    /*..Redefine iswtch to bypass setup..*/
    348:    /*....First digit of iswtch = 2: bypass setup and do not release memory
    349:          upon return....*/
    350:    /*....Second digit of iswtch = 1: memory extension switch....*/
    351:    /*....Third and fourth digit of iswtch: n_default. If n_default = 0,
    352:      the user has to set secondary parameters....*/
    353:    iswtch = 210;

    355:    /*..Call SAMG..*/
    356:    SAMG(&numnodes, &numnonzero, &nsys,
    357:         ia, ja, Asky, rhs, u_approx, iu, &ndiu, ip, &ndip, &matrix,
    358:         iscale, &res_in, &res_out, &ncyc_done, &ierr, &nsolve,
    359:         &ifirst, &eps, &ncyc, &iswtch, &a_cmplx, &g_cmplx,
    360:         &p_cmplx, &w_avrge, &chktol, &idump, &iout);

    362:    /*..Create auxilary integer array..*/
    363:    PetscMalloc(numnodes * sizeof(double),&cols);

    365:    for (I=0;I<numnodes;I++)
    366:        cols[I] = I;

    368:    /*..Store values computed by SAMG into the PETSc vector z..*/
    369:    VecSetValues(z,numnodes,cols,u_approx,INSERT_VALUES);
    370:

    372:    /*..Restore PETSc rhs vector..*/
    373:    VecRestoreArray(r, &vals_getarray);

    375:    PetscFree(cols);
    376:    PetscFree(rhs);
    377:    PetscFree(u_approx);

    379:    return 0;
    380: }

    382: /* ------------------------------------------------------------------- */
    385: /*..RamgShellPCDestroy - This routine destroys a user-defined
    386:     preconditioner context.

    388:     Input Parameter:
    389:     shell - user-defined preconditioner context..*/

    391: PetscErrorCode PETSCKSP_DLLEXPORT SamgShellPCDestroy(SamgShellPC *shell)
    392: {
    393:   /*..Free memory allocated by samg..*/
    394:   SAMG_cleanup();
    395:   /*..Free PCShell context..*/
    396:   PetscFree(shell->A);
    397:   PetscFree(shell->IA);
    398:   PetscFree(shell->JA);
    399:   PetscFree(shell->PARAM);
    400:   PetscFree(shell);

    402:   return 0;
    403: }
    404: /* ------------------------------------------------------------------- */
    407: /*..SamgGetParam - Gets SAMG parameters specified at runtime
    408:     OUTPUT: The parameters set in the SAMG_PARAM context
    409: ..*/

    411: PetscErrorCode PETSCKSP_DLLEXPORT SamgGetParam(SAMG_PARAM *samg_param)
    412: {

    415:   /*..Set default SAMG paramets..*/
    416:   /*....Class 0 SAMG parameters....*/
    417:   (*samg_param).MATRIX    = 12;
    418:   /*....Primary SAMG parameters....*/
    419:   /*......If ifirst=0, the vector u, as passed to SAMG, is taken as first
    420:           approximation......*/
    421:   (*samg_param).IFIRST    = 0;
    422:   (*samg_param).EPS       = 1e-12;
    423:   /*......nsolve =1 denotes scalar approach......*/
    424:   (*samg_param).NSOLVE    = 1;
    425:   /*......note: in the AMG-PETSc interface the number of cycles is required
    426:           to equal one to assure that in the PCApply routine AMG only performs
    427:           one cycle......*/
    428:   (*samg_param).NCYC      = 1001;
    429:   (*samg_param).ISWTCH    = 410;
    430:   (*samg_param).A_CMPLX   = 2.5;
    431:   (*samg_param).G_CMPLX   = 1.9;
    432:   (*samg_param).P_CMPLX   = 1.9;
    433:   (*samg_param).W_AVRGE   = 2.5;
    434:   (*samg_param).CHKTOL    = 1e-8;
    435:   (*samg_param).IDUMP     = 0;
    436:   (*samg_param).IOUT      = 2;
    437:   /*....Secundary SAMG parameters....*/
    438:   (*samg_param).LEVELX    = 25;
    439:   (*samg_param).NPTMN     = 100;
    440:   (*samg_param).ECG       = 0.25;
    441:   (*samg_param).EWT       = 0.5;
    442:   (*samg_param).NCG       = 1000;
    443:   (*samg_param).NWT       = 3000;
    444:   (*samg_param).ETR       = 12.2;
    445:   (*samg_param).NTR       = 2;
    446:   (*samg_param).NRD       = 131;
    447:   (*samg_param).NRC       = 0;
    448:   (*samg_param).NRU       = -131;

    450:   /*..Overwrite default values by values specified at runtime..*/
    451:   /*....Primary SAMG parameters....*/
    452:   PetscOptionsGetInt(PETSC_NULL,"-pc_samg_iswtch",&(*samg_param).ISWTCH,
    453:                        PETSC_NULL);

    455:   PetscOptionsGetInt(PETSC_NULL,"-pc_samg_ncyc",&(*samg_param).NCYC,
    456:                        PETSC_NULL);

    458:   PetscOptionsGetInt(PETSC_NULL,"-pc_samg_iout",&(*samg_param).IOUT,
    459:                        PETSC_NULL);

    461:  /*....Secundary SAMG parameters....*/
    462:   PetscOptionsGetInt(PETSC_NULL,"-pc_samg_levelx",&(*samg_param).LEVELX,
    463:                        PETSC_NULL);

    465:   PetscOptionsGetInt(PETSC_NULL,"-pc_samg_nptmn",&(*samg_param).NPTMN,
    466:                        PETSC_NULL);
    467:   PetscOptionsGetInt(PETSC_NULL,"-pc_samg_nrd",&(*samg_param).NRD,
    468:                        PETSC_NULL);

    470:   PetscOptionsGetInt(PETSC_NULL,"-pc_samg_nrc",&(*samg_param).NRC,
    471:                        PETSC_NULL);

    473:   PetscOptionsGetInt(PETSC_NULL,"-pc_samg_nru",&(*samg_param).NRU,
    474:                        PETSC_NULL);

    476:   return 0;
    477: }