!
!
!   This program demonstrates use of MatShellSetOperation()
!
      subroutine mymatmult(A, x, y, ierr)
#include <petsc/finclude/petscmat.h>
        use petscmat
        implicit none

        Mat A
        Vec x, y
        PetscErrorCode ierr

        print *, 'Called MatMult'
      end

      subroutine mymatmultadd(A, x, y, z, ierr)
        use petscmat
        implicit none
        Mat A
        Vec x, y, z
        PetscErrorCode ierr

        print *, 'Called MatMultAdd'
      end

      subroutine mymatmulttranspose(A, x, y, ierr)
        use petscmat
        implicit none
        Mat A
        Vec x, y
        PetscErrorCode ierr

        print *, 'Called MatMultTranspose'
      end

      subroutine mymatmulthermitiantranspose(A, x, y, ierr)
        use petscmat
        implicit none
        Mat A
        Vec x, y
        PetscErrorCode ierr

        print *, 'Called MatMultHermitianTranspose'
      end

      subroutine mymatmulttransposeadd(A, x, y, z, ierr)
        use petscmat
        implicit none
        Mat A
        Vec x, y, z
        PetscErrorCode ierr

        print *, 'Called MatMultTransposeAdd'
      end

      subroutine mymatmulthermitiantransposeadd(A, x, y, z, ierr)
        use petscmat
        implicit none
        Mat A
        Vec x, y, z
        PetscErrorCode ierr

        print *, 'Called MatMultHermitianTransposeAdd'
      end

      subroutine mymattranspose(A, reuse, B, ierr)
        use petscmat
        implicit none
        Mat A, B
        MatReuse reuse
        PetscErrorCode ierr
        PetscInt i12, i0

        i12 = 12
        i0 = 0
        PetscCallA(MatCreateShell(PETSC_COMM_SELF, i12, i12, i12, i12, i0, B, ierr))
        PetscCallA(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY, ierr))
        PetscCallA(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY, ierr))

        print *, 'Called MatTranspose'
      end

      subroutine mymatgetdiagonal(A, x, ierr)
        use petscmat
        implicit none
        Mat A
        Vec x
        PetscErrorCode ierr

        print *, 'Called MatGetDiagonal'
      end

      subroutine mymatdiagonalscale(A, x, y, ierr)
        use petscmat
        implicit none
        Mat A
        Vec x, y
        PetscErrorCode ierr

        print *, 'Called MatDiagonalScale'
      end

      subroutine mymatzeroentries(A, ierr)
        use petscmat
        implicit none
        Mat A
        PetscErrorCode ierr

        print *, 'Called MatZeroEntries'
      end

      subroutine mymataxpy(A, alpha, B, str, ierr)
        use petscmat
        implicit none
        Mat A, B
        PetscScalar alpha
        MatStructure str
        PetscErrorCode ierr

        print *, 'Called MatAXPY'
      end

      subroutine mymatshift(A, alpha, ierr)
        use petscmat
        implicit none
        Mat A
        PetscScalar alpha
        PetscErrorCode ierr

        print *, 'Called MatShift'
      end

      subroutine mymatdiagonalset(A, x, ins, ierr)
        use petscmat
        implicit none
        Mat A
        Vec x
        InsertMode ins
        PetscErrorCode ierr

        print *, 'Called MatDiagonalSet'
      end

      subroutine mymatdestroy(A, ierr)
        use petscmat
        implicit none
        Mat A
        PetscErrorCode ierr

        print *, 'Called MatDestroy'
      end

      subroutine mymatview(A, viewer, ierr)
        use petscmat
        implicit none
        Mat A
        PetscViewer viewer
        PetscErrorCode ierr

        print *, 'Called MatView'
      end

      subroutine mymatgetvecs(A, x, y, ierr)
        use petscmat
        implicit none
        Mat A
        Vec x, y
        PetscErrorCode ierr

        print *, 'Called MatCreateVecs'
      end

      program main
        use petscmat
        implicit none

        Mat m, mt
        Vec x, y, z
        PetscScalar a
        PetscViewer viewer
        MatOperation op
        PetscErrorCode ierr
        PetscInt i12, i0
        external mymatmult
        external mymatmultadd
        external mymatmulttranspose
        external mymatmulthermitiantranspose
        external mymatmulttransposeadd
        external mymatmulthermitiantransposeadd
        external mymattranspose
        external mymatgetdiagonal
        external mymatdiagonalscale
        external mymatzeroentries
        external mymataxpy
        external mymatshift
        external mymatdiagonalset
        external mymatdestroy
        external mymatview
        external mymatgetvecs

        PetscCallA(PetscInitialize(ierr))

        viewer = PETSC_VIEWER_STDOUT_SELF
        i12 = 12
        i0 = 0
        PetscCallA(VecCreateSeq(PETSC_COMM_SELF, i12, x, ierr))
        PetscCallA(VecCreateSeq(PETSC_COMM_SELF, i12, y, ierr))
        PetscCallA(VecCreateSeq(PETSC_COMM_SELF, i12, z, ierr))
        PetscCallA(MatCreateShell(PETSC_COMM_SELF, i12, i12, i12, i12, i0, m, ierr))
        PetscCallA(MatShellSetManageScalingShifts(m, ierr))
        PetscCallA(MatAssemblyBegin(m, MAT_FINAL_ASSEMBLY, ierr))
        PetscCallA(MatAssemblyEnd(m, MAT_FINAL_ASSEMBLY, ierr))

        op = MATOP_MULT
        PetscCallA(MatShellSetOperation(m, op, mymatmult, ierr))
        op = MATOP_MULT_ADD
        PetscCallA(MatShellSetOperation(m, op, mymatmultadd, ierr))
        op = MATOP_MULT_TRANSPOSE
        PetscCallA(MatShellSetOperation(m, op, mymatmulttranspose, ierr))
        op = MATOP_MULT_HERMITIAN_TRANSPOSE
        PetscCallA(MatShellSetOperation(m, op, mymatmulthermitiantranspose, ierr))
        op = MATOP_MULT_TRANSPOSE_ADD
        PetscCallA(MatShellSetOperation(m, op, mymatmulttransposeadd, ierr))
        op = MATOP_MULT_HERMITIAN_TRANS_ADD
        PetscCallA(MatShellSetOperation(m, op, mymatmulthermitiantransposeadd, ierr))
        op = MATOP_TRANSPOSE
        PetscCallA(MatShellSetOperation(m, op, mymattranspose, ierr))
        op = MATOP_GET_DIAGONAL
        PetscCallA(MatShellSetOperation(m, op, mymatgetdiagonal, ierr))
        op = MATOP_DIAGONAL_SCALE
        PetscCallA(MatShellSetOperation(m, op, mymatdiagonalscale, ierr))
        op = MATOP_ZERO_ENTRIES
        PetscCallA(MatShellSetOperation(m, op, mymatzeroentries, ierr))
        op = MATOP_AXPY
        PetscCallA(MatShellSetOperation(m, op, mymataxpy, ierr))
        op = MATOP_SHIFT
        PetscCallA(MatShellSetOperation(m, op, mymatshift, ierr))
        op = MATOP_DIAGONAL_SET
        PetscCallA(MatShellSetOperation(m, op, mymatdiagonalset, ierr))
        op = MATOP_DESTROY
        PetscCallA(MatShellSetOperation(m, op, mymatdestroy, ierr))
        op = MATOP_VIEW
        PetscCallA(MatShellSetOperation(m, op, mymatview, ierr))
        op = MATOP_CREATE_VECS
        PetscCallA(MatShellSetOperation(m, op, mymatgetvecs, ierr))

        PetscCallA(MatMult(m, x, y, ierr))
        PetscCallA(MatMultAdd(m, x, y, z, ierr))
        PetscCallA(MatMultTranspose(m, x, y, ierr))
        PetscCallA(MatMultHermitianTranspose(m, x, y, ierr))
        PetscCallA(MatMultTransposeAdd(m, x, y, z, ierr))
        PetscCallA(MatMultHermitianTransposeAdd(m, x, y, z, ierr))
        PetscCallA(MatTranspose(m, MAT_INITIAL_MATRIX, mt, ierr))
        PetscCallA(MatGetDiagonal(m, x, ierr))
        PetscCallA(MatDiagonalScale(m, x, y, ierr))
        PetscCallA(MatZeroEntries(m, ierr))
        a = 102.
        PetscCallA(MatAXPY(m, a, mt, SAME_NONZERO_PATTERN, ierr))
        PetscCallA(MatShift(m, a, ierr))
        PetscCallA(MatDiagonalSet(m, x, INSERT_VALUES, ierr))
        PetscCallA(MatView(m, viewer, ierr))
        PetscCallA(MatCreateVecs(m, x, y, ierr))
        PetscCallA(MatDestroy(m, ierr))
        PetscCallA(MatDestroy(mt, ierr))
        PetscCallA(VecDestroy(x, ierr))
        PetscCallA(VecDestroy(y, ierr))
        PetscCallA(VecDestroy(z, ierr))

        PetscCallA(PetscFinalize(ierr))
      end

!/*TEST
!
!   testset:
!     args: -malloc_dump
!     filter: sort -b
!     filter_output: sort -b
!     test:
!       suffix: 1
!       requires: !complex
!     test:
!       suffix: 2
!       requires: complex
!
!TEST*/