Actual source code: zshellf.c

      1: #include <petsc/private/ftnimpl.h>
      2: #include <petscmat.h>

      4: #if defined(PETSC_HAVE_FORTRAN_CAPS)
      5:   #define matshellsetoperation_ MATSHELLSETOPERATION
      6: #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE)
      7:   #define matshellsetoperation_ matshellsetoperation
      8: #endif

     10: /**
     11:  * Subset of MatOperation that is supported by the Fortran wrappers.
     12:  */
     13: enum FortranMatOperation {
     14:   FORTRAN_MATOP_MULT               = 0,
     15:   FORTRAN_MATOP_MULT_ADD           = 1,
     16:   FORTRAN_MATOP_MULT_TRANSPOSE     = 2,
     17:   FORTRAN_MATOP_MULT_TRANSPOSE_ADD = 3,
     18:   FORTRAN_MATOP_SOR                = 4,
     19:   FORTRAN_MATOP_TRANSPOSE          = 5,
     20:   FORTRAN_MATOP_GET_DIAGONAL       = 6,
     21:   FORTRAN_MATOP_DIAGONAL_SCALE     = 7,
     22:   FORTRAN_MATOP_ZERO_ENTRIES       = 8,
     23:   FORTRAN_MATOP_AXPY               = 9,
     24:   FORTRAN_MATOP_SHIFT              = 10,
     25:   FORTRAN_MATOP_DIAGONAL_SET       = 11,
     26:   FORTRAN_MATOP_DESTROY            = 12,
     27:   FORTRAN_MATOP_VIEW               = 13,
     28:   FORTRAN_MATOP_CREATE_VECS        = 14,
     29:   FORTRAN_MATOP_GET_DIAGONAL_BLOCK = 15,
     30:   FORTRAN_MATOP_COPY               = 16,
     31:   FORTRAN_MATOP_SCALE              = 17,
     32:   FORTRAN_MATOP_SET_RANDOM         = 18,
     33:   FORTRAN_MATOP_ASSEMBLY_BEGIN     = 19,
     34:   FORTRAN_MATOP_ASSEMBLY_END       = 20,
     35:   FORTRAN_MATOP_DUPLICATE          = 21,
     36:   FORTRAN_MATOP_MULT_HT            = 22,
     37:   FORTRAN_MATOP_MULT_HT_ADD        = 23,
     38:   FORTRAN_MATOP_SIZE               = 24
     39: };

     41: /*
     42:   The MatShell Matrix Vector product requires a C routine.
     43:   This C routine then calls the corresponding Fortran routine that was
     44:   set by the user.
     45: */
     46: static PetscErrorCode ourmult(Mat mat, Vec x, Vec y)
     47: {
     48:   PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_MULT]))(&mat, &x, &y, &ierr));
     49:   return PETSC_SUCCESS;
     50: }

     52: static PetscErrorCode ourmultadd(Mat mat, Vec x, Vec y, Vec z)
     53: {
     54:   PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_ADD]))(&mat, &x, &y, &z, &ierr));
     55:   return PETSC_SUCCESS;
     56: }

     58: static PetscErrorCode ourmulttranspose(Mat mat, Vec x, Vec y)
     59: {
     60:   PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_TRANSPOSE]))(&mat, &x, &y, &ierr));
     61:   return PETSC_SUCCESS;
     62: }

     64: static PetscErrorCode ourmulthermitiantranspose(Mat mat, Vec x, Vec y)
     65: {
     66:   PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_HT]))(&mat, &x, &y, &ierr));
     67:   return PETSC_SUCCESS;
     68: }

     70: static PetscErrorCode ourmulttransposeadd(Mat mat, Vec x, Vec y, Vec z)
     71: {
     72:   PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_TRANSPOSE_ADD]))(&mat, &x, &y, &z, &ierr));
     73:   return PETSC_SUCCESS;
     74: }

     76: static PetscErrorCode ourmulthermitiantransposeadd(Mat mat, Vec x, Vec y, Vec z)
     77: {
     78:   PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_HT_ADD]))(&mat, &x, &y, &z, &ierr));
     79:   return PETSC_SUCCESS;
     80: }

     82: static PetscErrorCode oursor(Mat mat, Vec b, PetscReal omega, MatSORType flg, PetscReal shift, PetscInt its, PetscInt lits, Vec x)
     83: {
     84:   PetscErrorCode ierr = PETSC_SUCCESS;

     86:   (*(void (*)(Mat *, Vec *, PetscReal *, MatSORType *, PetscReal *, PetscInt *, PetscInt *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_SOR]))(&mat, &b, &omega, &flg, &shift, &its, &lits, &x, &ierr);
     87:   return ierr;
     88: }

     90: static PetscErrorCode ourtranspose(Mat mat, MatReuse reuse, Mat *B)
     91: {
     92:   Mat  bb = (Mat)-1;
     93:   Mat *b  = (!B ? &bb : B);

     95:   PetscCallFortranVoidFunction((*(void (*)(Mat *, MatReuse *, Mat *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_TRANSPOSE]))(&mat, &reuse, b, &ierr));
     96:   return PETSC_SUCCESS;
     97: }

     99: static PetscErrorCode ourgetdiagonal(Mat mat, Vec x)
    100: {
    101:   PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_GET_DIAGONAL]))(&mat, &x, &ierr));
    102:   return PETSC_SUCCESS;
    103: }

    105: static PetscErrorCode ourdiagonalscale(Mat mat, Vec l, Vec r)
    106: {
    107:   Vec  aa = (Vec)-1;
    108:   Vec *a  = (!l ? &aa : &l);
    109:   Vec *b  = (!r ? &aa : &r);

    111:   PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_DIAGONAL_SCALE]))(&mat, a, b, &ierr));
    112:   return PETSC_SUCCESS;
    113: }

    115: static PetscErrorCode ourzeroentries(Mat mat)
    116: {
    117:   PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_ZERO_ENTRIES]))(&mat, &ierr));
    118:   return PETSC_SUCCESS;
    119: }

    121: static PetscErrorCode ouraxpy(Mat mat, PetscScalar a, Mat X, MatStructure str)
    122: {
    123:   PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscScalar *, Mat *, MatStructure *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_AXPY]))(&mat, &a, &X, &str, &ierr));
    124:   return PETSC_SUCCESS;
    125: }

    127: static PetscErrorCode ourshift(Mat mat, PetscScalar a)
    128: {
    129:   PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscScalar *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_SHIFT]))(&mat, &a, &ierr));
    130:   return PETSC_SUCCESS;
    131: }

    133: static PetscErrorCode ourdiagonalset(Mat mat, Vec x, InsertMode ins)
    134: {
    135:   PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, InsertMode *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_DIAGONAL_SET]))(&mat, &x, &ins, &ierr));
    136:   return PETSC_SUCCESS;
    137: }

    139: static PetscErrorCode ourdestroy(Mat mat)
    140: {
    141:   PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_DESTROY]))(&mat, &ierr));
    142:   return PETSC_SUCCESS;
    143: }

    145: static PetscErrorCode ourview(Mat mat, PetscViewer v)
    146: {
    147:   PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscViewer *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_VIEW]))(&mat, &v, &ierr));
    148:   return PETSC_SUCCESS;
    149: }

    151: static PetscErrorCode ourgetvecs(Mat mat, Vec *l, Vec *r)
    152: {
    153:   Vec  aa = (Vec)-1;
    154:   Vec *a  = (!l ? &aa : l);
    155:   Vec *b  = (!r ? &aa : r);

    157:   PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_CREATE_VECS]))(&mat, a, b, &ierr));
    158:   return PETSC_SUCCESS;
    159: }

    161: static PetscErrorCode ourgetdiagonalblock(Mat mat, Mat *l)
    162: {
    163:   PetscCallFortranVoidFunction((*(void (*)(Mat *, Mat *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_GET_DIAGONAL_BLOCK]))(&mat, l, &ierr));
    164:   return PETSC_SUCCESS;
    165: }

    167: static PetscErrorCode ourcopy(Mat mat, Mat B, MatStructure str)
    168: {
    169:   PetscCallFortranVoidFunction((*(void (*)(Mat *, Mat *, MatStructure *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_COPY]))(&mat, &B, &str, &ierr));
    170:   return PETSC_SUCCESS;
    171: }

    173: static PetscErrorCode ourscale(Mat mat, PetscScalar a)
    174: {
    175:   PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscScalar *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_SCALE]))(&mat, &a, &ierr));
    176:   return PETSC_SUCCESS;
    177: }

    179: static PetscErrorCode oursetrandom(Mat mat, PetscRandom ctx)
    180: {
    181:   PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscRandom *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_SET_RANDOM]))(&mat, &ctx, &ierr));
    182:   return PETSC_SUCCESS;
    183: }

    185: static PetscErrorCode ourassemblybegin(Mat mat, MatAssemblyType type)
    186: {
    187:   PetscCallFortranVoidFunction((*(void (*)(Mat *, MatAssemblyType *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_ASSEMBLY_BEGIN]))(&mat, &type, &ierr));
    188:   return PETSC_SUCCESS;
    189: }

    191: static PetscErrorCode ourassemblyend(Mat mat, MatAssemblyType type)
    192: {
    193:   PetscCallFortranVoidFunction((*(void (*)(Mat *, MatAssemblyType *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_ASSEMBLY_END]))(&mat, &type, &ierr));
    194:   return PETSC_SUCCESS;
    195: }

    197: static PetscErrorCode ourduplicate(Mat mat, MatDuplicateOption op, Mat *M)
    198: {
    199:   *((void **)(M)) = (void *)-2; // Initialize matrix since it will be passed to Fortran
    200:   PetscCallFortranVoidFunction((*(void (*)(Mat *, MatDuplicateOption *, Mat *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_DUPLICATE]))(&mat, &op, M, &ierr));
    201:   return PETSC_SUCCESS;
    202: }

    204: PETSC_EXTERN void matshellsetoperation_(Mat *mat, MatOperation *op, PetscErrorCode (*f)(Mat *, Vec *, Vec *, PetscErrorCode *), PetscErrorCode *ierr)
    205: {
    206:   MPI_Comm comm;

    208:   *ierr = PetscObjectGetComm((PetscObject)*mat, &comm);
    209:   if (*ierr) return;
    210:   PetscObjectAllocateFortranPointers(*mat, FORTRAN_MATOP_SIZE);

    212:   switch (*op) {
    213:   case MATOP_MULT:
    214:     *ierr                                                          = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourmult);
    215:     ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_MULT] = (PetscFortranCallbackFn *)f;
    216:     break;
    217:   case MATOP_MULT_ADD:
    218:     *ierr                                                              = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourmultadd);
    219:     ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_ADD] = (PetscFortranCallbackFn *)f;
    220:     break;
    221:   case MATOP_MULT_TRANSPOSE:
    222:     *ierr                                                                    = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourmulttranspose);
    223:     ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_TRANSPOSE] = (PetscFortranCallbackFn *)f;
    224:     break;
    225:   case MATOP_MULT_HERMITIAN_TRANSPOSE:
    226:     *ierr                                                             = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourmulthermitiantranspose);
    227:     ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_HT] = (PetscFortranCallbackFn *)f;
    228:     break;
    229:   case MATOP_MULT_TRANSPOSE_ADD:
    230:     *ierr                                                                        = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourmulttransposeadd);
    231:     ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_TRANSPOSE_ADD] = (PetscFortranCallbackFn *)f;
    232:     break;
    233:   case MATOP_MULT_HERMITIAN_TRANS_ADD:
    234:     *ierr                                                                 = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourmulthermitiantransposeadd);
    235:     ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_HT_ADD] = (PetscFortranCallbackFn *)f;
    236:     break;
    237:   case MATOP_SOR:
    238:     *ierr                                                         = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)oursor);
    239:     ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_SOR] = (PetscFortranCallbackFn *)f;
    240:     break;
    241:   case MATOP_TRANSPOSE:
    242:     *ierr                                                               = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourtranspose);
    243:     ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_TRANSPOSE] = (PetscFortranCallbackFn *)f;
    244:     break;
    245:   case MATOP_GET_DIAGONAL:
    246:     *ierr                                                                  = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourgetdiagonal);
    247:     ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_GET_DIAGONAL] = (PetscFortranCallbackFn *)f;
    248:     break;
    249:   case MATOP_DIAGONAL_SCALE:
    250:     *ierr                                                                    = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourdiagonalscale);
    251:     ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_DIAGONAL_SCALE] = (PetscFortranCallbackFn *)f;
    252:     break;
    253:   case MATOP_ZERO_ENTRIES:
    254:     *ierr                                                                  = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourzeroentries);
    255:     ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_ZERO_ENTRIES] = (PetscFortranCallbackFn *)f;
    256:     break;
    257:   case MATOP_AXPY:
    258:     *ierr                                                          = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ouraxpy);
    259:     ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_AXPY] = (PetscFortranCallbackFn *)f;
    260:     break;
    261:   case MATOP_SHIFT:
    262:     *ierr                                                           = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourshift);
    263:     ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_SHIFT] = (PetscFortranCallbackFn *)f;
    264:     break;
    265:   case MATOP_DIAGONAL_SET:
    266:     *ierr                                                                  = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourdiagonalset);
    267:     ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_DIAGONAL_SET] = (PetscFortranCallbackFn *)f;
    268:     break;
    269:   case MATOP_DESTROY:
    270:     *ierr                                                             = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourdestroy);
    271:     ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_DESTROY] = (PetscFortranCallbackFn *)f;
    272:     break;
    273:   case MATOP_VIEW:
    274:     *ierr                                                          = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourview);
    275:     ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_VIEW] = (PetscFortranCallbackFn *)f;
    276:     break;
    277:   case MATOP_CREATE_VECS:
    278:     *ierr                                                                 = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourgetvecs);
    279:     ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_CREATE_VECS] = (PetscFortranCallbackFn *)f;
    280:     break;
    281:   case MATOP_GET_DIAGONAL_BLOCK:
    282:     *ierr                                                                        = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourgetdiagonalblock);
    283:     ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_GET_DIAGONAL_BLOCK] = (PetscFortranCallbackFn *)f;
    284:     break;
    285:   case MATOP_COPY:
    286:     *ierr                                                          = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourcopy);
    287:     ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_COPY] = (PetscFortranCallbackFn *)f;
    288:     break;
    289:   case MATOP_SCALE:
    290:     *ierr                                                           = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourscale);
    291:     ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_SCALE] = (PetscFortranCallbackFn *)f;
    292:     break;
    293:   case MATOP_SET_RANDOM:
    294:     *ierr                                                                = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)oursetrandom);
    295:     ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_SET_RANDOM] = (PetscFortranCallbackFn *)f;
    296:     break;
    297:   case MATOP_ASSEMBLY_BEGIN:
    298:     *ierr                                                                    = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourassemblybegin);
    299:     ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_ASSEMBLY_BEGIN] = (PetscFortranCallbackFn *)f;
    300:     break;
    301:   case MATOP_ASSEMBLY_END:
    302:     *ierr                                                                  = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourassemblyend);
    303:     ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_ASSEMBLY_END] = (PetscFortranCallbackFn *)f;
    304:     break;
    305:   case MATOP_DUPLICATE:
    306:     *ierr                                                               = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourduplicate);
    307:     ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_DUPLICATE] = (PetscFortranCallbackFn *)f;
    308:     break;
    309:   default:
    310:     *ierr = PetscError(comm, __LINE__, "MatShellSetOperation_Fortran", __FILE__, PETSC_ERR_ARG_WRONG, PETSC_ERROR_INITIAL, "Cannot set that matrix operation");
    311:     *ierr = PETSC_ERR_ARG_WRONG;
    312:   }
    313: }