program test

  implicit none

#include <petsc/finclude/petscsys.h>
#include <petsc/finclude/petscvec.h>
#include <petsc/finclude/petscmat.h>
#include <petsc/finclude/petscviewer.h>


  PetscViewer    :: view
  ! sparse matrix
  Mat            :: A
  ! distributed dense matrix of size n x m
  Mat            :: B, X, R, QDlt, AQDlt
  ! distributed dense matrix of size n x (m x k)
  Mat            :: Q, K, AQ_p, AQ
  ! local dense matrix (every process keep the identical copies), (m x k) x (m x k)
  Mat            :: AConjPara, QtAQ, QtAQ_p, Dlt

  PetscInt       :: nDim, mDim, rhsNDim,rhsMDim,ierr, maxIter, iter, step_k,bsize
  PetscInt       :: ownRowS,ownRowE
  PetscScalar, allocatable :: XInit(:,:)
  PetscInt       :: XInitI, XInitJ
  PetscScalar    :: v=1.0
  PetscBool      :: flg
  PetscMPIInt    :: size, rank

  character(128) ::  fin, rhsfin


  call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
  call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
  call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)

  ! read binary matrix file
  call PetscOptionsGetString(PETSC_NULL_CHARACTER,'-f',fin,flg,ierr)
  call PetscOptionsGetString(PETSC_NULL_CHARACTER,'-r',rhsfin,flg,ierr)

  call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-i',maxIter,flg,ierr)
  call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-k',step_k,flg,ierr)
  call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-w',bsize,flg,ierr)


  call PetscViewerBinaryOpen(PETSC_COMM_WORLD,fin,FILE_MODE_READ,view,ierr)
  call MatCreate(PETSC_COMM_WORLD,A,ierr)
  call MatSetType(A,MATAIJ,ierr)
  call MatLoad(A,view,ierr)
  call PetscViewerDestroy(view,ierr)
  ! for the time being, assume mDim == nDim is true
  call MatGetSize(A, nDim, mDim, ierr)
  call MatView(A,PETSC_VIEWER_STDOUT_WORLD,ierr)    
  if (rank == 0) then
    print*,'Mat Size = ', nDim, mDim
  end if


  call MatGetOwnershipRange(A,ownRowS,ownRowE, ierr)

  ! create right-and-side matrix
  ! for the time being, choose row-wise decomposition
  ! for the time being, assume nDim%size = 0
  call MatCreateDense(PETSC_COMM_WORLD, (ownRowE - ownRowS), &
                      PETSC_DECIDE, nDim, bsize,PETSC_NULL_SCALAR, B, ierr)
  call PetscViewerBinaryOpen(PETSC_COMM_WORLD,rhsfin,FILE_MODE_READ,view, ierr)
  call MatLoad(B,view,ierr)
  call PetscViewerDestroy(view,ierr)
  call MatGetSize(B, rhsMDim, rhsNDim, ierr)
  if (rank == 0) then
    print*,'MRHS Size actually are:', rhsMDim, rhsNDim
    print*,'MRHS Size should be:', nDim, bsize
  end if

  call MatView(B,PETSC_VIEWER_STDOUT_WORLD,ierr)    
  ! inintial value guses X
  allocate(XInit(nDim,bsize))
  do XInitI=1, nDim
    do XInitJ=1, bsize
      XInit(XInitI,XInitJ) = 1.0
    end do
  end do

  call MatCreateDense(PETSC_COMM_WORLD, (ownRowE - ownRowS), &
                      PETSC_DECIDE, nDim, bsize,XInit, X, ierr)



  !  B, X, R, QDlt, AQDlt
  call MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, R, ierr)


  call MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, QDlt, ierr)


  call MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, AQDlt, ierr)


! Q, K, AQ_p, AQ of size n x (m x k)
  call MatCreateDense(PETSC_COMM_WORLD, (ownRowE - ownRowS), &
                      PETSC_DECIDE, nDim, (bsize*step_k),PETSC_NULL_SCALAR, Q, ierr)
  call MatAssemblyBegin(Q, MAT_FINAL_ASSEMBLY, ierr)
  call MatAssemblyEnd(Q, MAT_FINAL_ASSEMBLY, ierr)
 
  call MatDuplicate(Q, MAT_DO_NOT_COPY_VALUES, K, ierr)

  call MatDuplicate(Q, MAT_DO_NOT_COPY_VALUES, AQ_p, ierr)


  call MatDuplicate(Q, MAT_DO_NOT_COPY_VALUES, AQ, ierr)


! QtAQ, QtAQ_p, Dlt of size (m x k) x (m x k)
  call MatCreateSeqDense(PETSC_COMM_SELF,(bsize*step_k),(bsize*step_k),&
                         PETSC_NULL_SCALAR, QtAQ, ierr)

  call MatDuplicate(QtAQ, MAT_DO_NOT_COPY_VALUES, QtAQ_p    , ierr)

  call MatDuplicate(QtAQ, MAT_DO_NOT_COPY_VALUES, Dlt       , ierr)

  call MatDuplicate(QtAQ, MAT_DO_NOT_COPY_VALUES, AConjPara , ierr)

! calculation for R

! call matrix powers kernel
  call mpk_monomial (K, A, R, step_k, rank,size)

! destory matrices
  deallocate(XInit)

  call MatDestroy(B, ierr)
  call MatDestroy(X, ierr)
  call MatDestroy(R, ierr)
  call MatDestroy(QDlt, ierr)
  call MatDestroy(AQDlt, ierr)
  call MatDestroy(Q, ierr)
  call MatDestroy(K, ierr)
  call MatDestroy(AQ_p, ierr)
  call MatDestroy(AQ, ierr)
  call MatDestroy(QtAQ, ierr)
  call MatDestroy(QtAQ_p, ierr)
  call MatDestroy(Dlt, ierr)


  call PetscFinalize(ierr)

  stop

end program test