program mattrans

implicit none

#include <finclude/petscsys.h>
#include <finclude/petscis.h>
#include <finclude/petscvec.h>
#include <finclude/petscmat.h>
#include <finclude/petscpc.h>
#include <finclude/petscksp.h>
#include <finclude/petscdmda.h>

  MPI_Comm    :: comm
  PetscMPIInt :: myrank, mysize
  PetscErrorCode :: ierr
  PetscInt    :: nrows
  PetscInt    :: Mcols
  PetscInt    :: indices(4,4)
  PetscInt    :: ii,jj,offset
  PetscReal :: real_array4(4)
  Mat         :: A, Atrans
  PetscViewer :: viewer

  call PetscInitialize(PETSC_NULL_CHARACTER,ierr)

  comm = PETSC_COMM_WORLD
  call MPI_Comm_rank(comm,myrank,ierr)
  call MPI_Comm_rank(comm,mysize,ierr)

  if (myrank == 0) then
    nrows = 1
    indices(1,1) = 36
    indices(2,1) = 37
    indices(3,1) = 42
    indices(4,1) = 43
  else
    nrows = 4
    indices(1,1) = 37
    indices(2,1) = 38
    indices(3,1) = 43
    indices(4,1) = 44

    indices(1,2) = 38
    indices(2,2) = 39
    indices(3,2) = 44
    indices(4,2) = 45

    indices(1,3) = 39
    indices(2,3) = 40
    indices(3,3) = 45
    indices(4,3) = 46

    indices(1,4) = 40
    indices(2,4) = 41
    indices(3,4) = 46
    indices(4,4) = 47
  endif
  
  Mcols = 4

  call MatCreateMPIAIJ(comm, &
                       nrows, &
                       PETSC_DETERMINE, &
                       PETSC_DETERMINE, &
                       48, &
                       4, &
                       PETSC_NULL_INTEGER, &
                       4, &
                       PETSC_NULL_INTEGER, &
                       A, &
                       ierr)


  offset=0
  call MPI_Exscan(nrows,offset,1,MPIU_INTEGER,MPI_SUM,comm,ierr)
  real_array4 = 1.d0

  do ii = 1,nrows
      call MatSetValues(A,1,ii-1+offset, &
                        4,indices(:,ii),real_array4, &
                        INSERT_VALUES,ierr)
  enddo

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

  call PetscViewerASCIIOpen(comm,'A.out',viewer,ierr)
  call MatView(A,viewer,ierr)
  call PetscViewerDestroy(viewer,ierr)

  call MatTranspose(A,MAT_INITIAL_MATRIX,Atrans,ierr)

  call PetscViewerASCIIOpen(comm,'Atrans.out',viewer,ierr)
  call MatView(Atrans,viewer,ierr)
  call PetscViewerDestroy(viewer,ierr)

  call PetscFinalize(PETSC_NULL_CHARACTER,ierr)
end program mattrans