Actual source code: samgpetsctools.c

      1: #define PETSCKSP_DLL

     3:  #include global.h
     4:  #include petscksp.h
     5:  #include samgfunc.h
     6:  #include petscfunc.h
     7:  #include externc.h

     10: void USER_coo(int * i,int * ndim, double * x, double * y, double * z)
     11: {
     12:   printf("in user_coo");
     13: }

     16: static  double Machine_Precision_Eps = 2.e-16;
     17: /* ------------------------------------------------------------------- */
     20: /*..SamgGetGrid - This routine gets an array of grids
     21:     INPUT:  levels: number of levels created by SAMG
     22:             numnodes: number of nodes on finest grid
     23:             numnonzeros: number of nonzeros on coarsest grid
     24:     OUTPUT: grid  : array of grids                   ..*/
     25: PetscErrorCode PETSCKSP_DLLEXPORT SamgGetGrid(int levels, int numnodes, int numnonzero,
     26:                 GridCtx* grid, void* ctx)
     27: {
     28:    int      k;
     29:    int      ia_shift[MAX_LEVELS], ja_shift[MAX_LEVELS], nnu_cg, nna_cg;
     30:    int      iw_shift[MAX_LEVELS], jw_shift[MAX_LEVELS], rows_weights,
     31:             nna_weights, dummy;
     33:    MatInfo  info;

     35:    /*..Get coarse grid operators..*/
     36:    /*....Initialize ia_shift, ja_shift, nnu_cg and nna_cg....*/
     37:    ia_shift[1] = 1;
     38:    ja_shift[1] = 1;
     39:    nnu_cg = numnodes;
     40:    nna_cg = numnonzero;

     42:    for (k=2;k<=levels;k++){ /*....We do not get the finest level matrix....*/
     43:        /*....Update ia_shift and ja_shift values with nna_cg and nnu_cg
     44:              from previous loop....*/
     45:        ia_shift[k] = ia_shift[k-1] + nna_cg ;
     46:        ja_shift[k] = ja_shift[k-1] + nnu_cg ;

     48:        /*....Get coarse grid matrix on level k....*/
     49:        SamgGetCoarseMat(k, ia_shift[k], ja_shift[k], &(grid[k].A),
     50:                                PETSC_NULL);

     52:        /*....Get size and number of nonzeros of coarse grid matrix on
     53:              level k, i.e. get new nna_cg and nnu_cg values....*/
     54:        MatGetSize(grid[k].A, &nnu_cg, &nnu_cg);
     55:        MatGetInfo(grid[k].A, MAT_LOCAL, &info);
     56:        nna_cg = int(info.nz_used);
     57:    }
     58:
     59:    /*..Get interpolation operators..*/
     60:    /*....Initialize iw_shift, jw_shift and nna_weights....*/
     61:    iw_shift[0] = 1;
     62:    jw_shift[0] = 1;
     63:    nna_weights = 0;
     64:    rows_weights = numnodes;

     66:    for (k=1;k<=levels-1;k++){/*....There's NO interpolation operator
     67:                                    associated to the coarsest level....*/
     68:        /*....Update iw_shift with nna_weights value from
     69:              previous loop....*/
     70:        iw_shift[k] = iw_shift[k-1] + nna_weights ;
     71:        /*....Update jw_shift with rows_weights value from
     72:              current loop....*/
     73:        jw_shift[k] = jw_shift[k-1] + rows_weights ;
     74:
     75:        /*....Get interpolation from level k+1 to level k....*/
     76:        SamgGetInterpolation(k, iw_shift[k], jw_shift[k],
     77:                                    &(grid[k].Interp), PETSC_NULL) ;

     79:        /*....Get number of collumns and number of nonzeros of
     80:              interpolation associated to level k. NOTE: The
     81:              number of collums at this loop equals the number of
     82:              rows at the next loop...*/
     83:        MatGetSize(grid[k].Interp, &dummy, &rows_weights);
     84:        MatGetInfo(grid[k].Interp, MAT_LOCAL, &info);
     85:        nna_weights = int(info.nz_used);
     86:    }

     88:    return 0;
     89: }
     90: /* ------------------------------------------------------------------- */
     93: /*..SamgGetCoarseMat - This routine gets the coarse level matrix on the
     94:     level specified at input.
     95:     WARNING: This routine does not work to get the fine level matrix,
     96:              i.e. the value of k at input should be at least 2
     97:     INPUT:  level:    current grid level
     98:             ia_shift: shift to apply on ia_cg elements
     99:             ja_shift: shift to apply on ja_cg elements
    100:     OUTPUT: coarsemat: coarse level matrix
    101: ..*/

    103: PetscErrorCode PETSCKSP_DLLEXPORT SamgGetCoarseMat(int level, int ia_shift, int ja_shift,
    104:                      Mat* coarsemat, void* ctx)
    105: {
    106:    int      nnu_k, nna_k; /* size and non-zeros of operator on level k     */
    107:    int      *ia_k, *ja_k; /* coarse grid matrix in skyline format          */
    108:    double   *a_k;
    109:    int      *nnz_per_row; /* integer vector to hold the number of nonzeros */
    110:                           /* of each row. This vector will be used to      */
    111:                           /* allocate memory for the matrix, and to store  */
    112:                           /* elements in the matrix                        */
    114:    int      I;

    116:    /*..Get size (nnu_k) and number of non-zeros (nna_k) of operator
    117:      on level k..*/
    118:    SAMGPETSC_get_dim_operator(&level, &nnu_k, &nna_k);

    120:    /*..Now that nnu_cg and nna_cg are known, we can allocate memory for
    121:      coarse level matrix in compresses skyline format..*/
    122:    PetscMalloc(nna_k     * sizeof(double),&a_k);
    123:    PetscMalloc((nnu_k+1) * sizeof(int),&ia_k);
    124:    PetscMalloc(nna_k     * sizeof(int),&ja_k);

    126:    /*..Get coarse grid matrix in skyline format..*/
    127:    SAMGPETSC_get_operator(&level, a_k, ia_k, ja_k);

    129:    /*..Apply shift on each of the ia_cg and ja_cg elements..*/
    130:    SAMGPETSC_apply_shift(ia_k, &nnu_k, &ia_shift,
    131:                          ja_k, &nna_k, &ja_shift);

    133:    PetscMalloc(nnu_k * sizeof(int),&nnz_per_row);

    135:    /*..The numbero f nonzeros entries in row I can be calculated as
    136:        ia[I+1] - 1 - ia[I] + 1 = ia[I+1] - ia[I]                         ..*/
    137:    for (I=0;I<nnu_k;I++)
    138:        nnz_per_row[I] = ia_k[I+1] - ia_k[I];

    140:    /*..Allocate (create) SeqAIJ matrix  for use within PETSc..*/
    141:    MatCreate(PETSC_COMM_WORLD,coarsemat);
    142:    MatSetSizes(*coarsemat,nnu_k,nnu_k,nnu_k,nnu_k);
    143:    MatSetType(*coarsemat,MATSEQAIJ);
    144:    MatSeqAIJSetPreallocation(*coarsemat,0,nnz_per_row);

    146:    /*..Store coarse grid matrix in Petsc Mat object..*/
    147:    for (I=0;I<nnu_k;I++){
    148:       MatSetValues(*coarsemat,
    149:                1,              /* number of rows */
    150:                &I,             /* pointer to global row number */
    151:                nnz_per_row[I], /* number of collums = number of nonzero ... */
    152:                                /* entries in row I                          */
    153:                &(ja_k[ ia_k[I] ]),
    154:                               /* vector global column indices */
    155:                (PetscScalar *) &(a_k[ ia_k[I] ]),
    156:                               /* vector of coefficients */
    157:         INSERT_VALUES);
    158:    }

    160:    MatAssemblyBegin(*coarsemat,MAT_FINAL_ASSEMBLY);
    161:    MatAssemblyEnd(*coarsemat,MAT_FINAL_ASSEMBLY);

    163:    /*..Free memory required by storage in skyline format..*/
    164:    PetscFree(a_k);
    165:    PetscFree(ia_k);
    166:    PetscFree(ja_k);
    167:    /*..Free memory for auxilary array..*/
    168:    PetscFree(nnz_per_row);
    169:
    170:    return 0;
    171: }
    172: /* ------------------------------------------------------------------- */
    175: /*..SamgGetInterpolation - Get interpolation operator that interpolates
    176:     from level k+1 (coarse grid) to k (fine grid), where the input level
    177:     equals k ..*/
    178: /*..WARNING: This routine assumes that the input value for level is strictly
    179:     smaller than the number of levels created..*/
    180: /*..Implementation notes
    181:     o) The interpolation is a rectangular matrix with
    182:        number of rows equal to fine grid dim
    183:                  cols          coarse.
    184: ..*/
    185: PetscErrorCode PETSCKSP_DLLEXPORT SamgGetInterpolation(int level, int iw_shift, int jw_shift,
    186:                          Mat* interpolation, void* ctx)
    187: {
    188:    int     rows_weights, cols_weights, nna_weights;
    189:    int     *iweights, *jweights;
    190:    double  *weights;
    191:    int      *nnz_per_row; /* integer vector to hold the number of nonzeros */
    192:                           /* of each row. This vector will be used to      */
    193:                           /* allocate memory for the matrix, and to store  */
    194:                           /* elements in the matrix                        */
    196:    int      I,  coarser_level=level+1, dummy;

    198:    /*..Get number of rows and number of nonzeros of interpolation operator..*/
    199:    SAMGPETSC_get_dim_interpol(&level, &rows_weights, &nna_weights);

    201:    /*..Get number of cols of interpolation operator. NOTE: The number of
    202:        collums of the interpolation on level k equals the size of
    203:        the coarse grid matrix on the next coarsest grid.
    204:        SAMGPETSC_get_dim_interpol does not allow to get the number of
    205:        collumns of next to coarsest grid..*/
    206:    SAMGPETSC_get_dim_operator(&coarser_level, &cols_weights, &dummy);

    208:    /*..Now that nnu_weights and nna_weights are known, we can allocate
    209:        memory for interpolation operator in compresses skyline format..*/
    210:    PetscMalloc(nna_weights  * sizeof(double),&weights);
    211:
    212:    PetscMalloc((rows_weights+1) * sizeof(int),&iweights);
    213:
    214:    PetscMalloc(nna_weights  * sizeof(int),&jweights);
    215:

    217:    /*..Get interpolation operator in compressed skyline format..*/
    218:    SAMGPETSC_get_interpol(&level, weights, iweights, jweights);

    220:    /*..Apply shift on each of the ia_cg and ja_cg elements..*/
    221:    SAMGPETSC_apply_shift(iweights, &rows_weights, &iw_shift,
    222:                          jweights, &nna_weights, &jw_shift);

    224:    PetscMalloc(rows_weights * sizeof(int),&nnz_per_row);
    225:

    227:    /*..The numbero f nonzeros entries in row I can be calculated as
    228:        ia[I+1] - 1 - ia[I] + 1 = ia[I+1] - ia[I]                         ..*/
    229:    for (I=0;I<rows_weights;I++)
    230:        nnz_per_row[I] = iweights[I+1] - iweights[I];

    232:    /*..Allocate (create) SeqAIJ matrix  for use within PETSc..*/
    233:    MatCreate(PETSC_COMM_WORLD,interpolation);
    234:    MatSetSizes(*interpolation,rows_weights,cols_weights,rows_weights,cols_weights);
    235:    MatSetType(*interpolation,MATSEQAIJ);
    236:    MatSeqAIJSetPreallocation(*interpolation,0,nnz_per_row);

    238:    /*..Store coarse grid matrix in Petsc Mat object..*/
    239:    for (I=0;I<rows_weights;I++){
    240:       MatSetValues(*interpolation,
    241:                1,              /* number of rows */
    242:                &I,             /* pointer to global row number */
    243:                nnz_per_row[I], /* number of collums = number of nonzero ... */
    244:                                /* entries in row I                          */
    245:                &(jweights[ iweights[I] ]),
    246:                               /* vector global column indices */
    247:                (PetscScalar *) &(weights[ iweights[I] ]),
    248:                               /* vector of coefficients */
    249:         INSERT_VALUES);
    250:    }

    252:    MatAssemblyBegin(*interpolation,MAT_FINAL_ASSEMBLY);
    253:    MatAssemblyEnd(*interpolation,MAT_FINAL_ASSEMBLY);

    255:    /*..Free memory required by storage in skyline format..*/
    256:    PetscFree(weights);
    257:    PetscFree(iweights);
    258:    PetscFree(jweights);
    259:    /*..Free memory for auxilary array..*/
    260:    PetscFree(nnz_per_row);

    262:    return 0;
    263: }
    264: /* ------------------------------------------------------------------- */
    267: /*..Write coarser grid operators constructed by SAMG to ASCII file..*/

    269: PetscErrorCode PETSCKSP_DLLEXPORT SamgPetscWriteOperator(const int numnodes, const double* Asky,
    270:                            const int* ia, const int* ja, int extension)
    271: {

    273:      static char filename[80];
    274:      int         I,j,j1,j2;
    275:      FILE        *output;
    276:
    277:      /*..Switch arrays iacopy and jacopy to C conventions..*/
    278:      //     for (j=0;j<=numnodes;j++)
    279:      //         iacopy[j]--;
    280:      //     for (j=0;j<ia[numnodes];j++)
    281:      //         jacopy[j]--;
    282:
    283:      /*....Write matrix to file....*/
    284:       sprintf(filename,"coarsemat.%02u", extension);
    285:       output=fopen(filename,"w");
    286:       fprintf(output, "%% \n");
    287:       fprintf(output, "%% %d %d \n", numnodes, ia[numnodes] );

    289:       for (I=0;I<numnodes;I++){
    290:            j1 = ia[I];
    291:            j2 = ia[I+1] - 1;
    292:            for (j=j1;j<=j2;j++){
    293:                fprintf(output, "%d %d %22.18e\n", I+1, ja[j]+1,
    294:                                 Asky[j] );
    295:                //               printf("%d %d %e \n", I+1, ja[j]+1,
    296:                //                              Asky[j] );
    297:            }
    298:       }
    299:       fclose(output);
    300:
    301:       return 0;
    302: }
    303: /* ------------------------------------------------------------------- */
    306: /*..Write interpolation operators constructed by SAMG to ASCII file..*/

    308: PetscErrorCode PETSCKSP_DLLEXPORT SamgPetscWriteInterpol(const int numrows, const double* weights,
    309:                   const int* iweights, const int* jweights, int extension)
    310: {

    312:      static char filename[80];
    313:      int         I,j,j1,j2,numcols,numnonzero;
    314:      FILE        *output;
    315:
    316:      /*..Set number of nonzeros..*/
    317:      numnonzero = iweights[numrows];

    319:      /*..Determine numcols as the maximum ja value +1..*/
    320:      numcols = jweights[0];
    321:      for (j=0;j<numnonzero;j++){
    322:        if (jweights[j] > numcols) numcols = jweights[j];
    323:      }
    324:      numcols++;

    326:      /*..Write interpolation operator from grid k+1 (coarse grid) grid to k
    327:          (finer grid) to file..*/
    328:       sprintf(filename,"interpol.%02u%02u",
    329:                         extension+1, extension);
    330:       output=fopen(filename,"w");
    331:       fprintf(output, "%% \n%% %d %d %d \n", numrows, numcols, iweights[numrows] );
    332:       for (I=0;I<numrows;I++){
    333:            j1 = iweights[I];
    334:            j2 = iweights[I+1] - 1;
    335:            for (j=j1;j<=j2;j++){
    336:                fprintf(output, "%d %d %22.18e\n", I+1, jweights[j]+1,
    337:                                 weights[j] );
    338:                //               printf("%d %d %e \n", I+1, jweights[j]+1,
    339:                //                                weights[j] );
    340:            }
    341:       }
    342:       fclose(output);

    344: return 0;
    345: }
    346: /* ------------------------------------------------------------------- */
    349: /*..SamgCheckGalerkin - This routine offers a check on the correctness
    350:     of how SAMG interpolation and coarse grid operators are parsed to
    351:     PETSc. This routine computes I^H_h A^h I^h_H by PETSc matrix - matrix
    352:     multiplications, and compares this product with A^H..*/
    353:
    354: PetscErrorCode PETSCKSP_DLLEXPORT SamgCheckGalerkin(int levels, Mat A, GridCtx* grid,
    355:                       void* ctx)
    356: {
    357:    Mat     FineLevelMatrix, Restriction, HalfGalerkin, Galerkin, Diff;
    358:    double  normdiff;
    360:    int  k;

    362:    for (k=1;k<=levels-1;k++){
    363:       if (k==1)
    364:           FineLevelMatrix = A;
    365:           else
    366:           FineLevelMatrix = grid[k].A;
    367:       /*....Compute A^h I^h_H....*/
    368:       MatMatMult(FineLevelMatrix, grid[k].Interp, &HalfGalerkin);
    369:       /*....Get I^h_H....*/
    370:       MatTranspose(grid[k].Interp,&Restriction);
    371:       /*....Compute I^H_h A^h I^h_H....*/
    372:       MatMatMult(Restriction, HalfGalerkin, &Galerkin);
    373:       /*....Compute A^H - I^H_h A^h I^h_H....*/
    374:       MatSubstract(grid[k+1].A, Galerkin, &Diff);
    375:      /*....Compute || A^H - I^H_h A^h I^h_H||_{\infty}....*/
    376:       MatNorm(Diff,NORM_INFINITY,&normdiff);

    378:       printf("SamgCheckGalerkin :: || A^H - I^H_h A^h I^h_H||_{infty} on level %8d = %e\n",
    379:               k+1, normdiff);

    381:       MatDestroy(Restriction);
    382:       MatDestroy(HalfGalerkin);
    383:       MatDestroy(Galerkin);
    384:       MatDestroy(Diff);
    385:    }
    386:    return 0;
    387: }
    388: /* ------------------------------------------------------------------- */
    391: /*..MatSubstract - Computes the difference Term1 - Term2
    392:     INPUT:  Term1, Term2 : The input matrices
    393:     OUTPUT: Prod:          the difference
    394:   NOTE: Memory needed by difference has to freed outside this routine!
    395: ..*/

    397: PetscErrorCode MatSubstract(Mat Term1, Mat Term2, Mat* Diff)
    398: {
    399:    Vec          col_vec1, col_vec2, diff_vec;
    401:    int          rows1, cols1, rows2, cols2, col, row;
    402:    static PetscScalar dminusone = -1.;
    403:    PetscScalar  matrix_element ;
    404:    PetscScalar  *vec_getvalues ;
    405:    double       inf_norm_diff_vec = 0.0 ;
    406:    double       Zero_Element_Treshold = 0.0 ;

    408:    /*..Get sizes of terms..*/
    409:    MatGetSize(Term1, &rows1, &cols1);
    410:    MatGetSize(Term2, &rows2, &cols2);

    412:    /*..Check input..*/
    413:    if ( (cols1 != rows1) || (cols2 != rows2) ){
    414:       SETERRQ(PETSC_ERR_ARG_SIZ,"Error in MatMatMult: cols1 <> rows1 or cols1 <> rows1");
    415:    }

    417:    /*..Create difference of 2 SeqAIJ matrices..*/
    418:    MatCreate(PETSC_COMM_WORLD,Diff);
    419:    MatSetSizes(*Diff,rows1,cols1,rows1,cols1);
    420:    MatSetType(*Diff,MATSEQAIJ);
    421:    MatSeqAIJSetPreallocation(*Diff,0,PETSC_NULL);

    423:    /*..Create vectors..*/
    424:    VecCreate(MPI_COMM_WORLD,&col_vec1);
    425:    VecSetSizes(col_vec1,PETSC_DECIDE,rows1);
    426:    VecSetType(col_vec1,VECSEQ);
    427:    VecDuplicate(col_vec1, &col_vec2);
    428:    VecDuplicate(col_vec1, &diff_vec);

    430:    for (col=0;col<cols1;col++){
    431:        /*..Get collumns..*/
    432:       MatGetColumnVector(Term1,col_vec1,col);
    433:       MatGetColumnVector(Term2,col_vec2,col);
    434:       /*..Substract collumns..*/
    435:       VecWAXPY(diff_vec, dminusone, col_vec2, col_vec1);
    436:
    437:
    438:       /*..Compute norm..*/
    439:       VecNorm( diff_vec, NORM_INFINITY, &inf_norm_diff_vec );
    440:
    441:       /*..Set threshold..*/
    442:       Zero_Element_Treshold = inf_norm_diff_vec * Machine_Precision_Eps ;

    444:       /*..Get Term1(:,col) -  Term2(:,col) values..*/
    445:       VecGetArray(diff_vec, &vec_getvalues);

    447:       for (row=0;row<rows1;row++){
    448:           matrix_element = vec_getvalues[row];
    449:            if ( PetscAbsScalar( matrix_element ) >= Zero_Element_Treshold ) {
    450:                  MatSetValue(*Diff, row, col,matrix_element, INSERT_VALUES );
    451:            }
    452:       }
    453:       VecRestoreArray(diff_vec, &vec_getvalues);
    454:    }

    456:    MatAssemblyBegin(*Diff,MAT_FINAL_ASSEMBLY);
    457:    MatAssemblyEnd(*Diff,MAT_FINAL_ASSEMBLY);

    459:    VecDestroy(col_vec1);
    460:    VecDestroy(col_vec2);
    461:    VecDestroy(diff_vec);

    463:    return 0;
    464: }
    465: /* ------------------------------------------------------------------- */