155a74a43SLisandro Dalcin /* file: del2mat.h */
255a74a43SLisandro Dalcin
3*beceaeb6SBarry Smith #if !defined(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 */
Del2Mat_mult(Mat A,Vec x,Vec y)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;
304d86920dSPierre 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 */
Del2Mat_diag(Mat A,Vec D)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