subroutine ddistr(ictxt, M, N, blockSize, A0, A, descA) 
! Distribute double precision (selected_real_kind(8)) matrix A0 (MxN) from 
! root node to all processes in context ictxt. 
!
! A must be allocated like this:
! allocate(A(numroc(M, blockSize, myrow, 0, nprow), &
!            numroc(N, blockSize, mycol, 0, npcol)))
!
! Initial version 6.6.2002 Antti.Vanne@iki.fi
!


implicit none
integer, parameter :: dp = selected_real_kind(8), debug = 0
integer, intent(in) :: ictxt, M, N, blockSize 
integer :: myrow, mycol, rootNodeContext, nprow, npcol, isRootNode, &
     NRU, info, numroc
real(kind=dp), dimension(:,:) :: A
real(kind=dp), dimension(M,N), intent(in) :: A0
integer, dimension(9) :: descA, descA0

call sl_init (rootNodeContext, 1, 1) ! create 1 node context 
                                     ! for loading matrices
call blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol)

isRootNode = 0
if (myrow == 0 .and. mycol == 0) isRootNode = 1

!
! LOAD MATRIX ON ROOT NODE AND CREATE DESC FOR IT
!
if (isRootNode) then 
   NRU = numroc( M, M, MYROW, 0, NPROW )
   CALL DESCINIT( descA0, M, N, M, N, 0, 0, rootNodeContext, max(1, NRU), INFO )
else
   descA0(1:9) = 0
   descA0(2) = -1
end if 
!
!  CREATE DESC FOR DISTRIBUTED MATRIX
!
NRU = numroc( M, blockSize, MYROW, 0, NPROW )
CALL DESCINIT( descA, M, N, blockSize, blockSize, 0, 0, &
     ictxt, max(1, NRU), INFO )

!
!  DISTRIBUTE DATA
!
if (debug) write(*,*) "node r=", myrow, "c=", mycol, "M=", M, "N=", N
call PDGEMR2D( M, N, A0, 1, 1, descA0, A, 1, 1, descA, descA( 2 ) )

if (isRootNode) then 
   call blacs_gridexit( rootNodeContext ) 
end if

end subroutine ddistr