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 31 PetscFunctionBegin; 32 PetscCall(MatShellGetContext(A,&ctx)); 33 /* get raw vector arrays */ 34 PetscCall(VecGetArrayRead(x,&xx)); 35 PetscCall(VecGetArray(y,&yy)); 36 /* call external Fortran subroutine */ 37 Del2Apply(&ctx->N,ctx->F,xx,yy); 38 /* restore raw vector arrays */ 39 PetscCall(VecRestoreArrayRead(x,&xx)); 40 PetscCall(VecRestoreArray(y,&yy)); 41 PetscFunctionReturn(PETSC_SUCCESS); 42 } 43 44 /*D_i <- A_ii */ 45 PetscErrorCode Del2Mat_diag(Mat A, Vec D) 46 { 47 PetscFunctionBegin; 48 PetscCall(VecSet(D,6.0)); 49 PetscFunctionReturn(PETSC_SUCCESS); 50 } 51 52 #endif 53