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