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