1 /* file: del2mat.h */ 2 3 #ifndef DEL2MAT_H 4 #define DEL2MAT_H 5 6 #include <petsc.h> 7 #include <petscvec.h> 8 #include <petscmat.h> 9 10 /* external Fortran 90 subroutine */ 11 #define Del2Apply del2apply_ 12 EXTERN_C_BEGIN 13 extern void Del2Apply(int*,double*,const double*,double*); 14 EXTERN_C_END 15 16 /* user data structure and routines 17 * defining the matrix-free operator */ 18 19 typedef struct { 20 PetscInt N; 21 PetscScalar *F; 22 } Del2Mat; 23 24 /* y <- A * x */ 25 PetscErrorCode Del2Mat_mult(Mat A, Vec x, Vec y) 26 { 27 Del2Mat *ctx; 28 const PetscScalar *xx; 29 PetscScalar *yy; 30 PetscFunctionBegin; 31 PetscCall(MatShellGetContext(A,&ctx)); 32 /* get raw vector arrays */ 33 PetscCall(VecGetArrayRead(x,&xx)); 34 PetscCall(VecGetArray(y,&yy)); 35 /* call external Fortran subroutine */ 36 Del2Apply(&ctx->N,ctx->F,xx,yy); 37 /* restore raw vector arrays */ 38 PetscCall(VecRestoreArrayRead(x,&xx)); 39 PetscCall(VecRestoreArray(y,&yy)); 40 PetscFunctionReturn(PETSC_SUCCESS); 41 } 42 43 /*D_i <- A_ii */ 44 PetscErrorCode Del2Mat_diag(Mat A, Vec D) 45 { 46 PetscFunctionBegin; 47 PetscCall(VecSet(D,6.0)); 48 PetscFunctionReturn(PETSC_SUCCESS); 49 } 50 51 #endif 52