!
!
!   This program tests storage of PETSc Dense matrix.
!   It Creates a Dense matrix, and stores it into a file,
!   and then reads it back in as a SeqDense and MPIDense
!   matrix, and prints out the contents.
!
#include <petsc/finclude/petscmat.h>
program main
  use petscmat
  implicit none

  PetscErrorCode ierr
  PetscInt row, col, ten
  PetscMPIInt rank
  PetscScalar v
  Mat A
  PetscViewer view

  PetscCallA(PetscInitialize(ierr))

  PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
!
!     Proc-0 Create a seq-dense matrix and write it to a file
!
  if (rank == 0) then
    ten = 10
    PetscCallA(MatCreateSeqDense(PETSC_COMM_SELF, ten, ten, PETSC_NULL_SCALAR_ARRAY, A, ierr))
    v = 1.0
    do row = 0, 9
      do col = 0, 9
        PetscCallA(MatSetValue(A, row, col, v, INSERT_VALUES, ierr))
        v = v + 1.0
      end do
    end do

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

    PetscCallA(PetscObjectSetName(A, 'Original Matrix', ierr))
    PetscCallA(MatView(A, PETSC_VIEWER_STDOUT_SELF, ierr))
!
!        Now Write this matrix to a binary file
!
    PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_SELF, 'dense.mat', FILE_MODE_WRITE, view, ierr))
    PetscCallA(MatView(A, view, ierr))
    PetscCallA(PetscViewerDestroy(view, ierr))
    PetscCallA(MatDestroy(A, ierr))
!
!        Read this matrix into a SeqDense matrix

    PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_SELF, 'dense.mat', FILE_MODE_READ, view, ierr))
    PetscCallA(MatCreate(PETSC_COMM_SELF, A, ierr))
    PetscCallA(MatSetType(A, MATSEQDENSE, ierr))
    PetscCallA(MatLoad(A, view, ierr))

    PetscCallA(PetscObjectSetName(A, 'SeqDense Matrix read in from file', ierr))
    PetscCallA(MatView(A, PETSC_VIEWER_STDOUT_SELF, ierr))
    PetscCallA(MatDestroy(A, ierr))
    PetscCallA(PetscViewerDestroy(view, ierr))
  end if

!
!     Use a barrier, so that the procs do not try opening the file before
!     it is created.
!
  PetscCallMPIA(MPI_Barrier(PETSC_COMM_WORLD, ierr))

  PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_WORLD, 'dense.mat', FILE_MODE_READ, view, ierr))
  PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
  PetscCallA(MatSetType(A, MATMPIDENSE, ierr))
  PetscCallA(MatLoad(A, view, ierr))

  PetscCallA(PetscObjectSetName(A, 'MPIDense Matrix read in from file', ierr))
  PetscCallA(MatView(A, PETSC_VIEWER_STDOUT_WORLD, ierr))
  PetscCallA(MatDestroy(A, ierr))
  PetscCallA(PetscViewerDestroy(view, ierr))
  PetscCallA(PetscFinalize(ierr))
end

!/*TEST
!
!   test:
!      nsize: 2
!      output_file: output/ex63_1.out
!
!TEST*/