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_SIZE = 22 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 PETSC_EXTERN void matcreateshell_(MPI_Comm *comm, PetscInt *m, PetscInt *n, PetscInt *M, PetscInt *N, void *ctx, Mat *mat, PetscErrorCode *ierr) 47 { 48 *ierr = MatCreateShell(MPI_Comm_f2c(*(MPI_Fint *)&*comm), *m, *n, *M, *N, ctx, mat); 49 } 50 51 static PetscErrorCode ourmult(Mat mat, Vec x, Vec y) 52 { 53 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_MULT]))(&mat, &x, &y, &ierr)); 54 return PETSC_SUCCESS; 55 } 56 57 static PetscErrorCode ourmultadd(Mat mat, Vec x, Vec y, Vec z) 58 { 59 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_ADD]))(&mat, &x, &y, &z, &ierr)); 60 return PETSC_SUCCESS; 61 } 62 63 static PetscErrorCode ourmulttranspose(Mat mat, Vec x, Vec y) 64 { 65 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_TRANSPOSE]))(&mat, &x, &y, &ierr)); 66 return PETSC_SUCCESS; 67 } 68 69 static PetscErrorCode ourmulttransposeadd(Mat mat, Vec x, Vec y, Vec z) 70 { 71 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_TRANSPOSE_ADD]))(&mat, &x, &y, &z, &ierr)); 72 return PETSC_SUCCESS; 73 } 74 75 static PetscErrorCode oursor(Mat mat, Vec b, PetscReal omega, MatSORType flg, PetscReal shift, PetscInt its, PetscInt lits, Vec x) 76 { 77 PetscErrorCode ierr = PETSC_SUCCESS; 78 79 (*(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); 80 return ierr; 81 } 82 83 static PetscErrorCode ourtranspose(Mat mat, MatReuse reuse, Mat *B) 84 { 85 Mat bb = (Mat)-1; 86 Mat *b = (!B ? &bb : B); 87 88 PetscCallFortranVoidFunction((*(void (*)(Mat *, MatReuse *, Mat *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_TRANSPOSE]))(&mat, &reuse, b, &ierr)); 89 return PETSC_SUCCESS; 90 } 91 92 static PetscErrorCode ourgetdiagonal(Mat mat, Vec x) 93 { 94 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_GET_DIAGONAL]))(&mat, &x, &ierr)); 95 return PETSC_SUCCESS; 96 } 97 98 static PetscErrorCode ourdiagonalscale(Mat mat, Vec l, Vec r) 99 { 100 Vec aa = (Vec)-1; 101 Vec *a = (!l ? &aa : &l); 102 Vec *b = (!r ? &aa : &r); 103 104 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_DIAGONAL_SCALE]))(&mat, a, b, &ierr)); 105 return PETSC_SUCCESS; 106 } 107 108 static PetscErrorCode ourzeroentries(Mat mat) 109 { 110 PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_ZERO_ENTRIES]))(&mat, &ierr)); 111 return PETSC_SUCCESS; 112 } 113 114 static PetscErrorCode ouraxpy(Mat mat, PetscScalar a, Mat X, MatStructure str) 115 { 116 PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscScalar *, Mat *, MatStructure *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_AXPY]))(&mat, &a, &X, &str, &ierr)); 117 return PETSC_SUCCESS; 118 } 119 120 static PetscErrorCode ourshift(Mat mat, PetscScalar a) 121 { 122 PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscScalar *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_SHIFT]))(&mat, &a, &ierr)); 123 return PETSC_SUCCESS; 124 } 125 126 static PetscErrorCode ourdiagonalset(Mat mat, Vec x, InsertMode ins) 127 { 128 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, InsertMode *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_DIAGONAL_SET]))(&mat, &x, &ins, &ierr)); 129 return PETSC_SUCCESS; 130 } 131 132 static PetscErrorCode ourdestroy(Mat mat) 133 { 134 PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_DESTROY]))(&mat, &ierr)); 135 return PETSC_SUCCESS; 136 } 137 138 static PetscErrorCode ourview(Mat mat, PetscViewer v) 139 { 140 PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscViewer *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_VIEW]))(&mat, &v, &ierr)); 141 return PETSC_SUCCESS; 142 } 143 144 static PetscErrorCode ourgetvecs(Mat mat, Vec *l, Vec *r) 145 { 146 Vec aa = (Vec)-1; 147 Vec *a = (!l ? &aa : l); 148 Vec *b = (!r ? &aa : r); 149 150 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_CREATE_VECS]))(&mat, a, b, &ierr)); 151 return PETSC_SUCCESS; 152 } 153 154 static PetscErrorCode ourgetdiagonalblock(Mat mat, Mat *l) 155 { 156 PetscCallFortranVoidFunction((*(void (*)(Mat *, Mat *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_GET_DIAGONAL_BLOCK]))(&mat, l, &ierr)); 157 return PETSC_SUCCESS; 158 } 159 160 static PetscErrorCode ourcopy(Mat mat, Mat B, MatStructure str) 161 { 162 PetscCallFortranVoidFunction((*(void (*)(Mat *, Mat *, MatStructure *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_COPY]))(&mat, &B, &str, &ierr)); 163 return PETSC_SUCCESS; 164 } 165 166 static PetscErrorCode ourscale(Mat mat, PetscScalar a) 167 { 168 PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscScalar *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_SCALE]))(&mat, &a, &ierr)); 169 return PETSC_SUCCESS; 170 } 171 172 static PetscErrorCode oursetrandom(Mat mat, PetscRandom ctx) 173 { 174 PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscRandom *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_SET_RANDOM]))(&mat, &ctx, &ierr)); 175 return PETSC_SUCCESS; 176 } 177 178 static PetscErrorCode ourassemblybegin(Mat mat, MatAssemblyType type) 179 { 180 PetscCallFortranVoidFunction((*(void (*)(Mat *, MatAssemblyType *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_ASSEMBLY_BEGIN]))(&mat, &type, &ierr)); 181 return PETSC_SUCCESS; 182 } 183 184 static PetscErrorCode ourassemblyend(Mat mat, MatAssemblyType type) 185 { 186 PetscCallFortranVoidFunction((*(void (*)(Mat *, MatAssemblyType *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_ASSEMBLY_END]))(&mat, &type, &ierr)); 187 return PETSC_SUCCESS; 188 } 189 190 static PetscErrorCode ourduplicate(Mat mat, MatDuplicateOption op, Mat *M) 191 { 192 PetscCallFortranVoidFunction((*(void (*)(Mat *, MatDuplicateOption *, Mat *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_DUPLICATE]))(&mat, &op, M, &ierr)); 193 return PETSC_SUCCESS; 194 } 195 196 PETSC_EXTERN void matshellsetoperation_(Mat *mat, MatOperation *op, PetscErrorCode (*f)(Mat *, Vec *, Vec *, PetscErrorCode *), PetscErrorCode *ierr) 197 { 198 MPI_Comm comm; 199 200 *ierr = PetscObjectGetComm((PetscObject)*mat, &comm); 201 if (*ierr) return; 202 PetscObjectAllocateFortranPointers(*mat, FORTRAN_MATOP_SIZE); 203 204 switch (*op) { 205 case MATOP_MULT: 206 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFunction)ourmult); 207 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_MULT] = (PetscVoidFunction)f; 208 break; 209 case MATOP_MULT_ADD: 210 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFunction)ourmultadd); 211 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_ADD] = (PetscVoidFunction)f; 212 break; 213 case MATOP_MULT_TRANSPOSE: 214 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFunction)ourmulttranspose); 215 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_TRANSPOSE] = (PetscVoidFunction)f; 216 break; 217 case MATOP_MULT_TRANSPOSE_ADD: 218 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFunction)ourmulttransposeadd); 219 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_TRANSPOSE_ADD] = (PetscVoidFunction)f; 220 break; 221 case MATOP_SOR: 222 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFunction)oursor); 223 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_SOR] = (PetscVoidFunction)f; 224 break; 225 case MATOP_TRANSPOSE: 226 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFunction)ourtranspose); 227 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_TRANSPOSE] = (PetscVoidFunction)f; 228 break; 229 case MATOP_GET_DIAGONAL: 230 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFunction)ourgetdiagonal); 231 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_GET_DIAGONAL] = (PetscVoidFunction)f; 232 break; 233 case MATOP_DIAGONAL_SCALE: 234 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFunction)ourdiagonalscale); 235 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_DIAGONAL_SCALE] = (PetscVoidFunction)f; 236 break; 237 case MATOP_ZERO_ENTRIES: 238 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFunction)ourzeroentries); 239 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_ZERO_ENTRIES] = (PetscVoidFunction)f; 240 break; 241 case MATOP_AXPY: 242 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFunction)ouraxpy); 243 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_AXPY] = (PetscVoidFunction)f; 244 break; 245 case MATOP_SHIFT: 246 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFunction)ourshift); 247 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_SHIFT] = (PetscVoidFunction)f; 248 break; 249 case MATOP_DIAGONAL_SET: 250 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFunction)ourdiagonalset); 251 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_DIAGONAL_SET] = (PetscVoidFunction)f; 252 break; 253 case MATOP_DESTROY: 254 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFunction)ourdestroy); 255 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_DESTROY] = (PetscVoidFunction)f; 256 break; 257 case MATOP_VIEW: 258 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFunction)ourview); 259 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_VIEW] = (PetscVoidFunction)f; 260 break; 261 case MATOP_CREATE_VECS: 262 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFunction)ourgetvecs); 263 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_CREATE_VECS] = (PetscVoidFunction)f; 264 break; 265 case MATOP_GET_DIAGONAL_BLOCK: 266 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFunction)ourgetdiagonalblock); 267 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_GET_DIAGONAL_BLOCK] = (PetscVoidFunction)f; 268 break; 269 case MATOP_COPY: 270 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFunction)ourcopy); 271 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_COPY] = (PetscVoidFunction)f; 272 break; 273 case MATOP_SCALE: 274 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFunction)ourscale); 275 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_SCALE] = (PetscVoidFunction)f; 276 break; 277 case MATOP_SET_RANDOM: 278 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFunction)oursetrandom); 279 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_SET_RANDOM] = (PetscVoidFunction)f; 280 break; 281 case MATOP_ASSEMBLY_BEGIN: 282 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFunction)ourassemblybegin); 283 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_ASSEMBLY_BEGIN] = (PetscVoidFunction)f; 284 break; 285 case MATOP_ASSEMBLY_END: 286 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFunction)ourassemblyend); 287 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_ASSEMBLY_END] = (PetscVoidFunction)f; 288 break; 289 case MATOP_DUPLICATE: 290 *ierr = MatShellSetOperation(*mat, *op, (PetscVoidFunction)ourduplicate); 291 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_DUPLICATE] = (PetscVoidFunction)f; 292 break; 293 default: 294 *ierr = PetscError(comm, __LINE__, "MatShellSetOperation_Fortran", __FILE__, PETSC_ERR_ARG_WRONG, PETSC_ERROR_INITIAL, "Cannot set that matrix operation"); 295 *ierr = PETSC_ERR_ARG_WRONG; 296 } 297 } 298