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