!
!
!  Tests MatGetArray() on a dense matrix
!

      program main
      implicit none


! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!                    Include files
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
!  The following include statements are required for Fortran programs
!  that use PETSc vectors:
!     petsc.h  - base PETSc routines
!     petscvec.h - defines (INSERT_VALUES)
!     petscmat.h    - matrices
!     petscmat.h90  - to allow access to Fortran 90 features of matrices

#include "include/finclude/petsc.h"
#include "include/finclude/petscvec.h"
#include "include/finclude/petscmat.h"
#include "include/finclude/petscmat.h90"

      Mat A
      integer i,j,m,n,ierr
      integer rstart,rend
      PetscScalar  v
      PetscScalar, pointer :: array(:,:)


      call PetscInitialize(PETSC_NULL_CHARACTER,ierr)

      m = 3
      n = 2
!
!      Create a parallel dense matrix shared by all processors 
!
      call MatCreateMPIDense(PETSC_COMM_WORLD,PETSC_DECIDE,                &
     &                       PETSC_DECIDE,m,n,PETSC_NULL_SCALAR,A,ierr)

!
!     Set values into the matrix. All processors set all values.
!
      do 10, i=0,m-1
        do 20, j=0,n-1
          v = 9.d0/(i+j+1)
          call MatSetValues(A,1,i,1,j,v,INSERT_VALUES,ierr)
 20     continue
 10   continue

      call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr) 
      call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr) 

!
!       Print the matrix to the screen 
!
      call MatView(A,PETSC_VIEWER_STDOUT_WORLD,ierr) 


!
!      Print the local portion of the matrix to the screen
!
      call MatGetArrayF90(A,array,ierr) 
      call MatGetOwnershipRange(A,rstart,rend,ierr)
      call PetscSequentialPhaseBegin(PETSC_COMM_WORLD,1,ierr)
      do 30 i=1,rend-rstart
         write(6,100) (PetscRealPart(array(i,j)),j=1,n)
 30   continue
 100  format(2F6.2)

      call PetscSequentialPhaseEnd(PETSC_COMM_WORLD,1,ierr)

      call MatRestoreArrayF90(A,array,ierr)   
!
!      Free the space used by the matrix
!
      call MatDestroy(A,ierr) 
      call PetscFinalize(ierr)
      end

