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