subroutine mpk_monomial (K, A, R, step_k, rank, sizeMPI)
	implicit none

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

	Mat            :: K, Km(step_k)
	Mat            :: A, R
	PetscMPIInt    :: sizeMPI, rank
	PetscInt       :: nDim, bsize, step_k, local_RRow, local_RCol, genIdx
	PetscInt       :: ierr
	PetscInt       :: stepIdx, blockShift, localRsize
  PetscScalar    :: KArray(1), RArray(1), PetscScalarSize
  PetscOffset    :: KArrayOffset, RArrayOffset

	call MatGetSize(R, nDim, bsize, ierr)
 	if (rank == 0) then
  	  print*,'Mat Size = ', nDim, bsize
  	end if

  call MatGetArray(K,KArray,KArrayOffset,ierr)

  call MatGetLocalSize(R,local_RRow,local_RCol)
!   print *, "local_RRow,local_RCol", local_RRow,local_RCol

  ! get arry from R to add values to K(1)
  call MatGetArray(R,RArray,RArrayOffset,ierr)

  call MatCreateDense(PETSC_COMM_WORLD,  PETSC_DECIDE, &
                        PETSC_DECIDE , nDim, bsize,KArray(KArrayOffset + 1), Km(1), ierr)


!   call PetscMemmove(KArray(KArrayOffset + 1),RArray(RArrayOffset + 1) &
!                  ,local_RRow * local_RCol * STORAGE_SIZE(PetscScalarSize), ierr)

  localRsize = local_RRow * local_RCol
  do genIdx= 1, localRsize
    KArray(KArrayOffset + genIdx) = RArray(RArrayOffset + genIdx)
  end do


  call MatRestoreArray(R,RArray,RArrayOffset,ierr)

  call MatAssemblyBegin(Km(1), MAT_FINAL_ASSEMBLY, ierr)
  call MatAssemblyEnd  (Km(1), MAT_FINAL_ASSEMBLY, ierr)
  
  do stepIdx= 2, step_k

    blockShift = KArrayOffset + (stepIdx-1) * (local_RRow * local_RCol)

    call MatCreateDense(PETSC_COMM_WORLD,  PETSC_DECIDE, &
                        PETSC_DECIDE , nDim, bsize,KArray(blockShift+1), Km(stepIdx), ierr)
    call MatAssemblyBegin(Km(stepIdx), MAT_FINAL_ASSEMBLY, ierr)
    call MatAssemblyEnd  (Km(stepIdx), MAT_FINAL_ASSEMBLY, ierr)
    
  end do

  call MatRestoreArray(K,KArray,KArrayOffset,ierr)

!   do stepIdx= 2, step_k
  do stepIdx= 2,2 

    call MatMatMult(A,Km(stepIdx-1),MAT_REUSE_MATRIX,PETSC_DEFAULT_INTEGER,Km(stepIdx), ierr)
!     call MatMatMult(A,Km(stepIdx-1),MAT_INITIAL_MATRIX,PETSC_DEFAULT_INTEGER,Km(stepIdx), ierr)
  end do

!   call MatView(K,PETSC_VIEWER_STDOUT_WORLD,ierr)

end subroutine mpk_monomial