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

#include <petsc/finclude/petscsys.h>
#include <petsc/finclude/petscvec.h>
#include <petsc/finclude/petscmat.h>
#include <petsc/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,  localRsize
  PetscScalar    :: KArray(1), RArray(1), PetscScalarSize
  PetscOffset    :: KArrayOffset, RArrayOffset,blockShift

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

  call MatDenseGetArray(K,KArray,KArrayOffset,ierr)

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

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

  call MatCreateDense(PETSC_COMM_WORLD,  PETSC_DECIDE, &
      PETSC_DECIDE , nDim, bsize,KArray(KArrayOffset + 1), Km(1), ierr)
  call MatView(Km(1),PETSC_VIEWER_STDOUT_WORLD,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 MatDenseRestoreArray(R,RArray,RArrayOffset,ierr)

  
  do stepIdx= 2, step_k


    blockShift = KArrayOffset + (stepIdx-1) * (local_RRow * local_RCol)
      print*,local_RRow,local_RCol,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 MatView(Km(stepIdx),PETSC_VIEWER_STDOUT_WORLD,ierr)
    
  end do

  call MatDenseRestoreArray(K,KArray,KArrayOffset,ierr)

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


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

!   call MatView(K,PETSC_VIEWER_STDOUT_WORLD,ierr)

end subroutine mpk_monomial