155a74a43SLisandro Dalcin /* file: del2mat.h */ 255a74a43SLisandro Dalcin 355a74a43SLisandro Dalcin #ifndef DEL2MAT_H 455a74a43SLisandro Dalcin #define DEL2MAT_H 555a74a43SLisandro Dalcin 655a74a43SLisandro Dalcin #include <petsc.h> 755a74a43SLisandro Dalcin #include <petscvec.h> 855a74a43SLisandro Dalcin #include <petscmat.h> 955a74a43SLisandro Dalcin 1055a74a43SLisandro Dalcin /* external Fortran 90 subroutine */ 1155a74a43SLisandro Dalcin #define Del2Apply del2apply_ 1255a74a43SLisandro Dalcin EXTERN_C_BEGIN 1355a74a43SLisandro Dalcin extern void Del2Apply(int*,double*,const double*,double*); 1455a74a43SLisandro Dalcin EXTERN_C_END 1555a74a43SLisandro Dalcin 1655a74a43SLisandro Dalcin /* user data structure and routines 1755a74a43SLisandro Dalcin * defining the matrix-free operator */ 1855a74a43SLisandro Dalcin 1955a74a43SLisandro Dalcin typedef struct { 2055a74a43SLisandro Dalcin PetscInt N; 2155a74a43SLisandro Dalcin PetscScalar *F; 2255a74a43SLisandro Dalcin } Del2Mat; 2355a74a43SLisandro Dalcin 2455a74a43SLisandro Dalcin /* y <- A * x */ 2555a74a43SLisandro Dalcin PetscErrorCode Del2Mat_mult(Mat A, Vec x, Vec y) 2655a74a43SLisandro Dalcin { 2755a74a43SLisandro Dalcin Del2Mat *ctx; 2855a74a43SLisandro Dalcin const PetscScalar *xx; 2955a74a43SLisandro Dalcin PetscScalar *yy; 30*4d86920dSPierre Jolivet 3155a74a43SLisandro Dalcin PetscFunctionBegin; 3255a74a43SLisandro Dalcin PetscCall(MatShellGetContext(A,&ctx)); 3355a74a43SLisandro Dalcin /* get raw vector arrays */ 3455a74a43SLisandro Dalcin PetscCall(VecGetArrayRead(x,&xx)); 3555a74a43SLisandro Dalcin PetscCall(VecGetArray(y,&yy)); 3655a74a43SLisandro Dalcin /* call external Fortran subroutine */ 3755a74a43SLisandro Dalcin Del2Apply(&ctx->N,ctx->F,xx,yy); 3855a74a43SLisandro Dalcin /* restore raw vector arrays */ 3955a74a43SLisandro Dalcin PetscCall(VecRestoreArrayRead(x,&xx)); 4055a74a43SLisandro Dalcin PetscCall(VecRestoreArray(y,&yy)); 4155a74a43SLisandro Dalcin PetscFunctionReturn(PETSC_SUCCESS); 4255a74a43SLisandro Dalcin } 4355a74a43SLisandro Dalcin 4455a74a43SLisandro Dalcin /*D_i <- A_ii */ 4555a74a43SLisandro Dalcin PetscErrorCode Del2Mat_diag(Mat A, Vec D) 4655a74a43SLisandro Dalcin { 4755a74a43SLisandro Dalcin PetscFunctionBegin; 4855a74a43SLisandro Dalcin PetscCall(VecSet(D,6.0)); 4955a74a43SLisandro Dalcin PetscFunctionReturn(PETSC_SUCCESS); 5055a74a43SLisandro Dalcin } 5155a74a43SLisandro Dalcin 5255a74a43SLisandro Dalcin #endif 53