1 #include <petsc/private/ftnimpl.h> 2 #include <petscmat.h> 3 4 #if defined(PETSC_HAVE_FORTRAN_CAPS) 5 #define matshellsetoperation_ MATSHELLSETOPERATION 6 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 7 #define matshellsetoperation_ matshellsetoperation 8 #endif 9 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 }; 40 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 } 51 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 } 57 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 } 63 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 } 69 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 } 75 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 } 81 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; 85 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 } 89 90 static PetscErrorCode ourtranspose(Mat mat, MatReuse reuse, Mat *B) 91 { 92 Mat bb = (Mat)-1; 93 Mat *b = (!B ? &bb : B); 94 95 PetscCallFortranVoidFunction((*(void (*)(Mat *, MatReuse *, Mat *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_TRANSPOSE]))(&mat, &reuse, b, &ierr)); 96 return PETSC_SUCCESS; 97 } 98 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 } 104 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); 110 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 } 114 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 } 120 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 } 126 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 } 132 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 } 138 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 } 144 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 } 150 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); 156 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 } 160 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 } 166 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 } 172 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 } 178 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 } 184 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 } 190 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 } 196 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 } 203 204 PETSC_EXTERN void matshellsetoperation_(Mat *mat, MatOperation *op, PetscErrorCode (*f)(Mat *, Vec *, Vec *, PetscErrorCode *), PetscErrorCode *ierr) 205 { 206 MPI_Comm comm; 207 208 *ierr = PetscObjectGetComm((PetscObject)*mat, &comm); 209 if (*ierr) return; 210 PetscObjectAllocateFortranPointers(*mat, FORTRAN_MATOP_SIZE); 211 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 } 314