1 #include <petsc/private/ftnimpl.h>
2 #include <petscmat.h>
3
4 #if defined(PETSC_HAVE_FORTRAN_CAPS)
5 #define matshellsetoperation_ MATSHELLSETOPERATION
6 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE)
7 #define matshellsetoperation_ matshellsetoperation
8 #endif
9
10 /**
11 * Subset of MatOperation that is supported by the Fortran wrappers.
12 */
13 enum FortranMatOperation {
14 FORTRAN_MATOP_MULT = 0,
15 FORTRAN_MATOP_MULT_ADD = 1,
16 FORTRAN_MATOP_MULT_TRANSPOSE = 2,
17 FORTRAN_MATOP_MULT_TRANSPOSE_ADD = 3,
18 FORTRAN_MATOP_SOR = 4,
19 FORTRAN_MATOP_TRANSPOSE = 5,
20 FORTRAN_MATOP_GET_DIAGONAL = 6,
21 FORTRAN_MATOP_DIAGONAL_SCALE = 7,
22 FORTRAN_MATOP_ZERO_ENTRIES = 8,
23 FORTRAN_MATOP_AXPY = 9,
24 FORTRAN_MATOP_SHIFT = 10,
25 FORTRAN_MATOP_DIAGONAL_SET = 11,
26 FORTRAN_MATOP_DESTROY = 12,
27 FORTRAN_MATOP_VIEW = 13,
28 FORTRAN_MATOP_CREATE_VECS = 14,
29 FORTRAN_MATOP_GET_DIAGONAL_BLOCK = 15,
30 FORTRAN_MATOP_COPY = 16,
31 FORTRAN_MATOP_SCALE = 17,
32 FORTRAN_MATOP_SET_RANDOM = 18,
33 FORTRAN_MATOP_ASSEMBLY_BEGIN = 19,
34 FORTRAN_MATOP_ASSEMBLY_END = 20,
35 FORTRAN_MATOP_DUPLICATE = 21,
36 FORTRAN_MATOP_MULT_HT = 22,
37 FORTRAN_MATOP_MULT_HT_ADD = 23,
38 FORTRAN_MATOP_SIZE = 24
39 };
40
41 /*
42 The MatShell Matrix Vector product requires a C routine.
43 This C routine then calls the corresponding Fortran routine that was
44 set by the user.
45 */
ourmult(Mat mat,Vec x,Vec y)46 static PetscErrorCode ourmult(Mat mat, Vec x, Vec y)
47 {
48 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_MULT]))(&mat, &x, &y, &ierr));
49 return PETSC_SUCCESS;
50 }
51
ourmultadd(Mat mat,Vec x,Vec y,Vec z)52 static PetscErrorCode ourmultadd(Mat mat, Vec x, Vec y, Vec z)
53 {
54 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_ADD]))(&mat, &x, &y, &z, &ierr));
55 return PETSC_SUCCESS;
56 }
57
ourmulttranspose(Mat mat,Vec x,Vec y)58 static PetscErrorCode ourmulttranspose(Mat mat, Vec x, Vec y)
59 {
60 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_TRANSPOSE]))(&mat, &x, &y, &ierr));
61 return PETSC_SUCCESS;
62 }
63
ourmulthermitiantranspose(Mat mat,Vec x,Vec y)64 static PetscErrorCode ourmulthermitiantranspose(Mat mat, Vec x, Vec y)
65 {
66 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_HT]))(&mat, &x, &y, &ierr));
67 return PETSC_SUCCESS;
68 }
69
ourmulttransposeadd(Mat mat,Vec x,Vec y,Vec z)70 static PetscErrorCode ourmulttransposeadd(Mat mat, Vec x, Vec y, Vec z)
71 {
72 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_TRANSPOSE_ADD]))(&mat, &x, &y, &z, &ierr));
73 return PETSC_SUCCESS;
74 }
75
ourmulthermitiantransposeadd(Mat mat,Vec x,Vec y,Vec z)76 static PetscErrorCode ourmulthermitiantransposeadd(Mat mat, Vec x, Vec y, Vec z)
77 {
78 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_HT_ADD]))(&mat, &x, &y, &z, &ierr));
79 return PETSC_SUCCESS;
80 }
81
oursor(Mat mat,Vec b,PetscReal omega,MatSORType flg,PetscReal shift,PetscInt its,PetscInt lits,Vec x)82 static PetscErrorCode oursor(Mat mat, Vec b, PetscReal omega, MatSORType flg, PetscReal shift, PetscInt its, PetscInt lits, Vec x)
83 {
84 PetscErrorCode ierr = PETSC_SUCCESS;
85
86 (*(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);
87 return ierr;
88 }
89
ourtranspose(Mat mat,MatReuse reuse,Mat * B)90 static PetscErrorCode ourtranspose(Mat mat, MatReuse reuse, Mat *B)
91 {
92 Mat bb = (Mat)-1;
93 Mat *b = (!B ? &bb : B);
94
95 PetscCallFortranVoidFunction((*(void (*)(Mat *, MatReuse *, Mat *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_TRANSPOSE]))(&mat, &reuse, b, &ierr));
96 return PETSC_SUCCESS;
97 }
98
ourgetdiagonal(Mat mat,Vec x)99 static PetscErrorCode ourgetdiagonal(Mat mat, Vec x)
100 {
101 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_GET_DIAGONAL]))(&mat, &x, &ierr));
102 return PETSC_SUCCESS;
103 }
104
ourdiagonalscale(Mat mat,Vec l,Vec r)105 static PetscErrorCode ourdiagonalscale(Mat mat, Vec l, Vec r)
106 {
107 Vec aa = (Vec)-1;
108 Vec *a = (!l ? &aa : &l);
109 Vec *b = (!r ? &aa : &r);
110
111 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_DIAGONAL_SCALE]))(&mat, a, b, &ierr));
112 return PETSC_SUCCESS;
113 }
114
ourzeroentries(Mat mat)115 static PetscErrorCode ourzeroentries(Mat mat)
116 {
117 PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_ZERO_ENTRIES]))(&mat, &ierr));
118 return PETSC_SUCCESS;
119 }
120
ouraxpy(Mat mat,PetscScalar a,Mat X,MatStructure str)121 static PetscErrorCode ouraxpy(Mat mat, PetscScalar a, Mat X, MatStructure str)
122 {
123 PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscScalar *, Mat *, MatStructure *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_AXPY]))(&mat, &a, &X, &str, &ierr));
124 return PETSC_SUCCESS;
125 }
126
ourshift(Mat mat,PetscScalar a)127 static PetscErrorCode ourshift(Mat mat, PetscScalar a)
128 {
129 PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscScalar *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_SHIFT]))(&mat, &a, &ierr));
130 return PETSC_SUCCESS;
131 }
132
ourdiagonalset(Mat mat,Vec x,InsertMode ins)133 static PetscErrorCode ourdiagonalset(Mat mat, Vec x, InsertMode ins)
134 {
135 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, InsertMode *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_DIAGONAL_SET]))(&mat, &x, &ins, &ierr));
136 return PETSC_SUCCESS;
137 }
138
ourdestroy(Mat mat)139 static PetscErrorCode ourdestroy(Mat mat)
140 {
141 PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_DESTROY]))(&mat, &ierr));
142 return PETSC_SUCCESS;
143 }
144
ourview(Mat mat,PetscViewer v)145 static PetscErrorCode ourview(Mat mat, PetscViewer v)
146 {
147 PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscViewer *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_VIEW]))(&mat, &v, &ierr));
148 return PETSC_SUCCESS;
149 }
150
ourgetvecs(Mat mat,Vec * l,Vec * r)151 static PetscErrorCode ourgetvecs(Mat mat, Vec *l, Vec *r)
152 {
153 Vec aa = (Vec)-1;
154 Vec *a = (!l ? &aa : l);
155 Vec *b = (!r ? &aa : r);
156
157 PetscCallFortranVoidFunction((*(void (*)(Mat *, Vec *, Vec *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_CREATE_VECS]))(&mat, a, b, &ierr));
158 return PETSC_SUCCESS;
159 }
160
ourgetdiagonalblock(Mat mat,Mat * l)161 static PetscErrorCode ourgetdiagonalblock(Mat mat, Mat *l)
162 {
163 PetscCallFortranVoidFunction((*(void (*)(Mat *, Mat *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_GET_DIAGONAL_BLOCK]))(&mat, l, &ierr));
164 return PETSC_SUCCESS;
165 }
166
ourcopy(Mat mat,Mat B,MatStructure str)167 static PetscErrorCode ourcopy(Mat mat, Mat B, MatStructure str)
168 {
169 PetscCallFortranVoidFunction((*(void (*)(Mat *, Mat *, MatStructure *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_COPY]))(&mat, &B, &str, &ierr));
170 return PETSC_SUCCESS;
171 }
172
ourscale(Mat mat,PetscScalar a)173 static PetscErrorCode ourscale(Mat mat, PetscScalar a)
174 {
175 PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscScalar *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_SCALE]))(&mat, &a, &ierr));
176 return PETSC_SUCCESS;
177 }
178
oursetrandom(Mat mat,PetscRandom ctx)179 static PetscErrorCode oursetrandom(Mat mat, PetscRandom ctx)
180 {
181 PetscCallFortranVoidFunction((*(void (*)(Mat *, PetscRandom *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_SET_RANDOM]))(&mat, &ctx, &ierr));
182 return PETSC_SUCCESS;
183 }
184
ourassemblybegin(Mat mat,MatAssemblyType type)185 static PetscErrorCode ourassemblybegin(Mat mat, MatAssemblyType type)
186 {
187 PetscCallFortranVoidFunction((*(void (*)(Mat *, MatAssemblyType *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_ASSEMBLY_BEGIN]))(&mat, &type, &ierr));
188 return PETSC_SUCCESS;
189 }
190
ourassemblyend(Mat mat,MatAssemblyType type)191 static PetscErrorCode ourassemblyend(Mat mat, MatAssemblyType type)
192 {
193 PetscCallFortranVoidFunction((*(void (*)(Mat *, MatAssemblyType *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_ASSEMBLY_END]))(&mat, &type, &ierr));
194 return PETSC_SUCCESS;
195 }
196
ourduplicate(Mat mat,MatDuplicateOption op,Mat * M)197 static PetscErrorCode ourduplicate(Mat mat, MatDuplicateOption op, Mat *M)
198 {
199 *((void **)(M)) = (void *)-2; // Initialize matrix since it will be passed to Fortran
200 PetscCallFortranVoidFunction((*(void (*)(Mat *, MatDuplicateOption *, Mat *, PetscErrorCode *))(((PetscObject)mat)->fortran_func_pointers[FORTRAN_MATOP_DUPLICATE]))(&mat, &op, M, &ierr));
201 return PETSC_SUCCESS;
202 }
203
matshellsetoperation_(Mat * mat,MatOperation * op,PetscErrorCode (* f)(Mat *,Vec *,Vec *,PetscErrorCode *),PetscErrorCode * ierr)204 PETSC_EXTERN void matshellsetoperation_(Mat *mat, MatOperation *op, PetscErrorCode (*f)(Mat *, Vec *, Vec *, PetscErrorCode *), PetscErrorCode *ierr)
205 {
206 MPI_Comm comm;
207
208 *ierr = PetscObjectGetComm((PetscObject)*mat, &comm);
209 if (*ierr) return;
210 PetscObjectAllocateFortranPointers(*mat, FORTRAN_MATOP_SIZE);
211
212 switch (*op) {
213 case MATOP_MULT:
214 *ierr = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourmult);
215 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_MULT] = (PetscFortranCallbackFn *)f;
216 break;
217 case MATOP_MULT_ADD:
218 *ierr = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourmultadd);
219 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_ADD] = (PetscFortranCallbackFn *)f;
220 break;
221 case MATOP_MULT_TRANSPOSE:
222 *ierr = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourmulttranspose);
223 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_TRANSPOSE] = (PetscFortranCallbackFn *)f;
224 break;
225 case MATOP_MULT_HERMITIAN_TRANSPOSE:
226 *ierr = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourmulthermitiantranspose);
227 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_HT] = (PetscFortranCallbackFn *)f;
228 break;
229 case MATOP_MULT_TRANSPOSE_ADD:
230 *ierr = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourmulttransposeadd);
231 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_TRANSPOSE_ADD] = (PetscFortranCallbackFn *)f;
232 break;
233 case MATOP_MULT_HERMITIAN_TRANS_ADD:
234 *ierr = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourmulthermitiantransposeadd);
235 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_MULT_HT_ADD] = (PetscFortranCallbackFn *)f;
236 break;
237 case MATOP_SOR:
238 *ierr = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)oursor);
239 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_SOR] = (PetscFortranCallbackFn *)f;
240 break;
241 case MATOP_TRANSPOSE:
242 *ierr = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourtranspose);
243 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_TRANSPOSE] = (PetscFortranCallbackFn *)f;
244 break;
245 case MATOP_GET_DIAGONAL:
246 *ierr = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourgetdiagonal);
247 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_GET_DIAGONAL] = (PetscFortranCallbackFn *)f;
248 break;
249 case MATOP_DIAGONAL_SCALE:
250 *ierr = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourdiagonalscale);
251 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_DIAGONAL_SCALE] = (PetscFortranCallbackFn *)f;
252 break;
253 case MATOP_ZERO_ENTRIES:
254 *ierr = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourzeroentries);
255 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_ZERO_ENTRIES] = (PetscFortranCallbackFn *)f;
256 break;
257 case MATOP_AXPY:
258 *ierr = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ouraxpy);
259 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_AXPY] = (PetscFortranCallbackFn *)f;
260 break;
261 case MATOP_SHIFT:
262 *ierr = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourshift);
263 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_SHIFT] = (PetscFortranCallbackFn *)f;
264 break;
265 case MATOP_DIAGONAL_SET:
266 *ierr = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourdiagonalset);
267 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_DIAGONAL_SET] = (PetscFortranCallbackFn *)f;
268 break;
269 case MATOP_DESTROY:
270 *ierr = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourdestroy);
271 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_DESTROY] = (PetscFortranCallbackFn *)f;
272 break;
273 case MATOP_VIEW:
274 *ierr = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourview);
275 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_VIEW] = (PetscFortranCallbackFn *)f;
276 break;
277 case MATOP_CREATE_VECS:
278 *ierr = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourgetvecs);
279 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_CREATE_VECS] = (PetscFortranCallbackFn *)f;
280 break;
281 case MATOP_GET_DIAGONAL_BLOCK:
282 *ierr = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourgetdiagonalblock);
283 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_GET_DIAGONAL_BLOCK] = (PetscFortranCallbackFn *)f;
284 break;
285 case MATOP_COPY:
286 *ierr = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourcopy);
287 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_COPY] = (PetscFortranCallbackFn *)f;
288 break;
289 case MATOP_SCALE:
290 *ierr = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourscale);
291 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_SCALE] = (PetscFortranCallbackFn *)f;
292 break;
293 case MATOP_SET_RANDOM:
294 *ierr = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)oursetrandom);
295 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_SET_RANDOM] = (PetscFortranCallbackFn *)f;
296 break;
297 case MATOP_ASSEMBLY_BEGIN:
298 *ierr = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourassemblybegin);
299 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_ASSEMBLY_BEGIN] = (PetscFortranCallbackFn *)f;
300 break;
301 case MATOP_ASSEMBLY_END:
302 *ierr = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourassemblyend);
303 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_ASSEMBLY_END] = (PetscFortranCallbackFn *)f;
304 break;
305 case MATOP_DUPLICATE:
306 *ierr = MatShellSetOperation(*mat, *op, (PetscErrorCodeFn *)ourduplicate);
307 ((PetscObject)*mat)->fortran_func_pointers[FORTRAN_MATOP_DUPLICATE] = (PetscFortranCallbackFn *)f;
308 break;
309 default:
310 *ierr = PetscError(comm, __LINE__, "MatShellSetOperation_Fortran", __FILE__, PETSC_ERR_ARG_WRONG, PETSC_ERROR_INITIAL, "Cannot set that matrix operation");
311 *ierr = PETSC_ERR_ARG_WRONG;
312 }
313 }
314