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