1 #include <petsc/private/fortranimpl.h> 2 #include <petscmat.h> 3 4 #if defined(PETSC_HAVE_FORTRAN_CAPS) 5 #define matshellsetoperation_ MATSHELLSETOPERATION 6 #define matcreateshell_ MATCREATESHELL 7 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 8 #define matcreateshell_ matcreateshell 9 #define matshellsetoperation_ matshellsetoperation 10 #endif 11 12 /** 13 * Subset of MatOperation that is supported by the Fortran wrappers. 14 */ 15 enum FortranMatOperation { 16 FORTRAN_MATOP_MULT = 0, 17 FORTRAN_MATOP_MULT_ADD = 1, 18 FORTRAN_MATOP_MULT_TRANSPOSE = 2, 19 FORTRAN_MATOP_MULT_TRANSPOSE_ADD = 3, 20 FORTRAN_MATOP_SOR = 4, 21 FORTRAN_MATOP_TRANSPOSE = 5, 22 FORTRAN_MATOP_GET_DIAGONAL = 6, 23 FORTRAN_MATOP_DIAGONAL_SCALE = 7, 24 FORTRAN_MATOP_ZERO_ENTRIES = 8, 25 FORTRAN_MATOP_AXPY = 9, 26 FORTRAN_MATOP_SHIFT = 10, 27 FORTRAN_MATOP_DIAGONAL_SET = 11, 28 FORTRAN_MATOP_DESTROY = 12, 29 FORTRAN_MATOP_VIEW = 13, 30 FORTRAN_MATOP_CREATE_VECS = 14, 31 FORTRAN_MATOP_GET_DIAGONAL_BLOCK = 15, 32 FORTRAN_MATOP_COPY = 16, 33 FORTRAN_MATOP_SCALE = 17, 34 FORTRAN_MATOP_SET_RANDOM = 18, 35 FORTRAN_MATOP_ASSEMBLY_BEGIN = 19, 36 FORTRAN_MATOP_ASSEMBLY_END = 20, 37 FORTRAN_MATOP_DUPLICATE = 21, 38 FORTRAN_MATOP_MULT_HT = 22, 39 FORTRAN_MATOP_MULT_HT_ADD = 23, 40 FORTRAN_MATOP_SIZE = 24 41 }; 42 43 /* 44 The MatShell Matrix Vector product requires a C routine. 45 This C routine then calls the corresponding Fortran routine that was 46 set by the user. 47 */ 48 PETSC_EXTERN void matcreateshell_(MPI_Comm *comm, PetscInt *m, PetscInt *n, PetscInt *M, PetscInt *N, void *ctx, Mat *mat, PetscErrorCode *ierr) 49 { 50 *ierr = MatCreateShell(MPI_Comm_f2c(*(MPI_Fint *)&*comm), *m, *n, *M, *N, ctx, mat); 51 } 52 53 static PetscErrorCode ourmult(Mat mat, Vec x, Vec y) 54 { 55 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_MULT]))(&mat, &x, &y, &ierr)); 56 return PETSC_SUCCESS; 57 } 58 59 static PetscErrorCode ourmultadd(Mat mat, Vec x, Vec y, Vec z) 60 { 61 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_ADD]))(&mat, &x, &y, &z, &ierr)); 62 return PETSC_SUCCESS; 63 } 64 65 static PetscErrorCode ourmulttranspose(Mat mat, Vec x, Vec y) 66 { 67 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_TRANSPOSE]))(&mat, &x, &y, &ierr)); 68 return PETSC_SUCCESS; 69 } 70 71 static PetscErrorCode ourmulthermitiantranspose(Mat mat, Vec x, Vec y) 72 { 73 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_HT]))(&mat, &x, &y, &ierr)); 74 return PETSC_SUCCESS; 75 } 76 77 static PetscErrorCode ourmulttransposeadd(Mat mat, Vec x, Vec y, Vec z) 78 { 79 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_TRANSPOSE_ADD]))(&mat, &x, &y, &z, &ierr)); 80 return PETSC_SUCCESS; 81 } 82 83 static PetscErrorCode ourmulthermitiantransposeadd(Mat mat, Vec x, Vec y, Vec z) 84 { 85 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_HT_ADD]))(&mat, &x, &y, &z, &ierr)); 86 return PETSC_SUCCESS; 87 } 88 89 static PetscErrorCode oursor(Mat mat, Vec b, PetscReal omega, MatSORType flg, PetscReal shift, PetscInt its, PetscInt lits, Vec x) 90 { 91 PetscErrorCode ierr = PETSC_SUCCESS; 92 93 (*(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); 94 return ierr; 95 } 96 97 static PetscErrorCode ourtranspose(Mat mat, MatReuse reuse, Mat *B) 98 { 99 Mat bb = (Mat)-1; 100 Mat *b = (!B ? &bb : B); 101 102 PetscCallFortranVoidFunction((*(void (*)(Mat *, MatReuse *, Mat *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_TRANSPOSE]))(&mat, &reuse, b, &ierr)); 103 return PETSC_SUCCESS; 104 } 105 106 static PetscErrorCode ourgetdiagonal(Mat mat, Vec x) 107 { 108 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_GET_DIAGONAL]))(&mat, &x, &ierr)); 109 return PETSC_SUCCESS; 110 } 111 112 static PetscErrorCode ourdiagonalscale(Mat mat, Vec l, Vec r) 113 { 114 Vec aa = (Vec)-1; 115 Vec *a = (!l ? &aa : &l); 116 Vec *b = (!r ? &aa : &r); 117 118 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_DIAGONAL_SCALE]))(&mat, a, b, &ierr)); 119 return PETSC_SUCCESS; 120 } 121 122 static PetscErrorCode ourzeroentries(Mat mat) 123 { 124 PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_ZERO_ENTRIES]))(&mat, &ierr)); 125 return PETSC_SUCCESS; 126 } 127 128 static PetscErrorCode ouraxpy(Mat mat, PetscScalar a, Mat X, MatStructure str) 129 { 130 PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscScalar *, Mat *, MatStructure *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_AXPY]))(&mat, &a, &X, &str, &ierr)); 131 return PETSC_SUCCESS; 132 } 133 134 static PetscErrorCode ourshift(Mat mat, PetscScalar a) 135 { 136 PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscScalar *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_SHIFT]))(&mat, &a, &ierr)); 137 return PETSC_SUCCESS; 138 } 139 140 static PetscErrorCode ourdiagonalset(Mat mat, Vec x, InsertMode ins) 141 { 142 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, InsertMode *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_DIAGONAL_SET]))(&mat, &x, &ins, &ierr)); 143 return PETSC_SUCCESS; 144 } 145 146 static PetscErrorCode ourdestroy(Mat mat) 147 { 148 PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_DESTROY]))(&mat, &ierr)); 149 return PETSC_SUCCESS; 150 } 151 152 static PetscErrorCode ourview(Mat mat, PetscViewer v) 153 { 154 PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscViewer *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_VIEW]))(&mat, &v, &ierr)); 155 return PETSC_SUCCESS; 156 } 157 158 static PetscErrorCode ourgetvecs(Mat mat, Vec *l, Vec *r) 159 { 160 Vec aa = (Vec)-1; 161 Vec *a = (!l ? &aa : l); 162 Vec *b = (!r ? &aa : r); 163 164 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_CREATE_VECS]))(&mat, a, b, &ierr)); 165 return PETSC_SUCCESS; 166 } 167 168 static PetscErrorCode ourgetdiagonalblock(Mat mat, Mat *l) 169 { 170 PetscCallFortranVoidFunction((*(void (*)(Mat *, Mat *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_GET_DIAGONAL_BLOCK]))(&mat, l, &ierr)); 171 return PETSC_SUCCESS; 172 } 173 174 static PetscErrorCode ourcopy(Mat mat, Mat B, MatStructure str) 175 { 176 PetscCallFortranVoidFunction((*(void (*)(Mat *, Mat *, MatStructure *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_COPY]))(&mat, &B, &str, &ierr)); 177 return PETSC_SUCCESS; 178 } 179 180 static PetscErrorCode ourscale(Mat mat, PetscScalar a) 181 { 182 PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscScalar *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_SCALE]))(&mat, &a, &ierr)); 183 return PETSC_SUCCESS; 184 } 185 186 static PetscErrorCode oursetrandom(Mat mat, PetscRandom ctx) 187 { 188 PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscRandom *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_SET_RANDOM]))(&mat, &ctx, &ierr)); 189 return PETSC_SUCCESS; 190 } 191 192 static PetscErrorCode ourassemblybegin(Mat mat, MatAssemblyType type) 193 { 194 PetscCallFortranVoidFunction((*(void (*)(Mat *, MatAssemblyType *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_ASSEMBLY_BEGIN]))(&mat, &type, &ierr)); 195 return PETSC_SUCCESS; 196 } 197 198 static PetscErrorCode ourassemblyend(Mat mat, MatAssemblyType type) 199 { 200 PetscCallFortranVoidFunction((*(void (*)(Mat *, MatAssemblyType *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_ASSEMBLY_END]))(&mat, &type, &ierr)); 201 return PETSC_SUCCESS; 202 } 203 204 static PetscErrorCode ourduplicate(Mat mat, MatDuplicateOption op, Mat *M) 205 { 206 PetscCallFortranVoidFunction((*(void (*)(Mat *, MatDuplicateOption *, Mat *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_DUPLICATE]))(&mat, &op, M, &ierr)); 207 return PETSC_SUCCESS; 208 } 209 210 PETSC_EXTERN void matshellsetoperation_(Mat *mat, MatOperation *op, PetscErrorCode (*f)(Mat *, Vec *, Vec *, PetscErrorCode *), PetscErrorCode *ierr) 211 { 212 MPI_Comm comm; 213 214 *ierr = PetscObjectGetComm((PetscObject)*mat, &comm); 215 if (*ierr) return; 216 PetscObjectAllocateFortranPointers(*mat, FORTRAN_MATOP_SIZE); 217 218 switch (*op) { 219 case MATOP_MULT: 220 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFn *)ourmult); 221 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_MULT] = (PetscVoidFn *)f; 222 break; 223 case MATOP_MULT_ADD: 224 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFn *)ourmultadd); 225 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_ADD] = (PetscVoidFn *)f; 226 break; 227 case MATOP_MULT_TRANSPOSE: 228 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFn *)ourmulttranspose); 229 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_TRANSPOSE] = (PetscVoidFn *)f; 230 break; 231 case MATOP_MULT_HERMITIAN_TRANSPOSE: 232 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFn *)ourmulthermitiantranspose); 233 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_HT] = (PetscVoidFn *)f; 234 break; 235 case MATOP_MULT_TRANSPOSE_ADD: 236 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFn *)ourmulttransposeadd); 237 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_TRANSPOSE_ADD] = (PetscVoidFn *)f; 238 break; 239 case MATOP_MULT_HERMITIAN_TRANS_ADD: 240 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFn *)ourmulthermitiantransposeadd); 241 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_HT_ADD] = (PetscVoidFn *)f; 242 break; 243 case MATOP_SOR: 244 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFn *)oursor); 245 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_SOR] = (PetscVoidFn *)f; 246 break; 247 case MATOP_TRANSPOSE: 248 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFn *)ourtranspose); 249 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_TRANSPOSE] = (PetscVoidFn *)f; 250 break; 251 case MATOP_GET_DIAGONAL: 252 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFn *)ourgetdiagonal); 253 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_GET_DIAGONAL] = (PetscVoidFn *)f; 254 break; 255 case MATOP_DIAGONAL_SCALE: 256 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFn *)ourdiagonalscale); 257 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_DIAGONAL_SCALE] = (PetscVoidFn *)f; 258 break; 259 case MATOP_ZERO_ENTRIES: 260 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFn *)ourzeroentries); 261 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_ZERO_ENTRIES] = (PetscVoidFn *)f; 262 break; 263 case MATOP_AXPY: 264 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFn *)ouraxpy); 265 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_AXPY] = (PetscVoidFn *)f; 266 break; 267 case MATOP_SHIFT: 268 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFn *)ourshift); 269 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_SHIFT] = (PetscVoidFn *)f; 270 break; 271 case MATOP_DIAGONAL_SET: 272 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFn *)ourdiagonalset); 273 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_DIAGONAL_SET] = (PetscVoidFn *)f; 274 break; 275 case MATOP_DESTROY: 276 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFn *)ourdestroy); 277 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_DESTROY] = (PetscVoidFn *)f; 278 break; 279 case MATOP_VIEW: 280 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFn *)ourview); 281 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_VIEW] = (PetscVoidFn *)f; 282 break; 283 case MATOP_CREATE_VECS: 284 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFn *)ourgetvecs); 285 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_CREATE_VECS] = (PetscVoidFn *)f; 286 break; 287 case MATOP_GET_DIAGONAL_BLOCK: 288 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFn *)ourgetdiagonalblock); 289 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_GET_DIAGONAL_BLOCK] = (PetscVoidFn *)f; 290 break; 291 case MATOP_COPY: 292 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFn *)ourcopy); 293 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_COPY] = (PetscVoidFn *)f; 294 break; 295 case MATOP_SCALE: 296 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFn *)ourscale); 297 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_SCALE] = (PetscVoidFn *)f; 298 break; 299 case MATOP_SET_RANDOM: 300 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFn *)oursetrandom); 301 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_SET_RANDOM] = (PetscVoidFn *)f; 302 break; 303 case MATOP_ASSEMBLY_BEGIN: 304 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFn *)ourassemblybegin); 305 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_ASSEMBLY_BEGIN] = (PetscVoidFn *)f; 306 break; 307 case MATOP_ASSEMBLY_END: 308 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFn *)ourassemblyend); 309 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_ASSEMBLY_END] = (PetscVoidFn *)f; 310 break; 311 case MATOP_DUPLICATE: 312 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFn *)ourduplicate); 313 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_DUPLICATE] = (PetscVoidFn *)f; 314 break; 315 default: 316 *ierr = PetscError(comm, __LINE__, "MatShellSetOperation_Fortran", __FILE__, PETSC_ERR_ARG_WRONG, PETSC_ERROR_INITIAL, "Cannot set that matrix operation"); 317 *ierr = PETSC_ERR_ARG_WRONG; 318 } 319 } 320