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