!
!
!   This program demonstrates use of PETSc dense matrices.
!
      program main
#include <petsc/finclude/petscsys.h>
      use petscsys
      implicit none

      PetscErrorCode ierr

      PetscCallA(PetscInitialize(ierr))

!  Demo of PETSc-allocated dense matrix storage
      call Demo1()

!  Demo of user-allocated dense matrix storage
      call Demo2()

      PetscCallA(PetscFinalize(ierr))
      end

! -----------------------------------------------------------------
!
!  Demo1 -  This subroutine demonstrates the use of PETSc-allocated dense
!  matrix storage.  Here MatDenseGetArrayF90() is used for direct access to the
!  array that stores the dense matrix.
!
!  Note the use of PETSC_NULL_SCALAR in MatCreateSeqDense() to indicate that no
!  storage is being provided by the user. (Do NOT pass a zero in that
!  location.)
!
      subroutine Demo1()
#include <petsc/finclude/petscmat.h>
      use petscmat
      implicit none

      Mat         A
      PetscInt   n,m
      PetscErrorCode ierr
      PetscScalar,pointer :: aa(:,:)

      n = 4
      m = 5

!  Create matrix

      PetscCall(MatCreate(PETSC_COMM_SELF,A,ierr))
      PetscCall(MatSetSizes(A,m,n,m,n,ierr))
      PetscCall(MatSetType(A,MATSEQDENSE,ierr))
      PetscCall(MatSetUp(A,ierr))

!  Access array storage
      PetscCall(MatDenseGetArrayF90(A,aa,ierr))

!  Set matrix values directly
      PetscCall(FillUpMatrix(m,n,aa))

      PetscCall(MatDenseRestoreArrayF90(A,aa,ierr))

!  Finalize matrix assembly
      PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr))
      PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr))

!  View matrix
      PetscCall(MatView(A,PETSC_VIEWER_STDOUT_SELF,ierr))

!  Clean up
      PetscCall(MatDestroy(A,ierr))
      end

! -----------------------------------------------------------------
!
!  Demo2 -  This subroutine demonstrates the use of user-allocated dense
!  matrix storage.
!
      subroutine Demo2()
#include <petsc/finclude/petscmat.h>
      use petscmat
      implicit none

      PetscInt   n,m
      PetscErrorCode ierr
      parameter (m=5,n=4)
      Mat       A
      PetscScalar    aa(m,n)

!  Create matrix
      PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,m,n,aa,A,ierr))

!  Set matrix values directly
      PetscCall(FillUpMatrix(m,n,aa))

!  Finalize matrix assembly
      PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr))
      PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr))

!  View matrix
      PetscCall(MatView(A,PETSC_VIEWER_STDOUT_SELF,ierr))

!  Clean up
      PetscCall(MatDestroy(A,ierr))
      end

! -----------------------------------------------------------------

      subroutine FillUpMatrix(m,n,X)
      PetscInt          m,n,i,j
      PetscScalar      X(m,n)

      do 10, j=1,n
        do 20, i=1,m
          X(i,j) = 1.0/real(i+j-1)
 20     continue
 10   continue
      end