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