!--------------------------------------------
!CAF=/usr/lib/x86_64-linux-gnu/open-coarrays/openmpi/lib    (locate libcaf_mpi) 
!mpifort -fcoarray=lib 03_array.f90  -lscalapack-openmpi  -L$CAF     -lcaf_mpi -o m ; mpirun -n 4 ./m ; head -n 100 out0*[0-9]
!
! If too many processes, then:   
!  n=4;  mpirun -n $n --host $(hostname):$n ./m | sort 
!--------------------------------------------
program     test_scalapack
 use, intrinsic :: iso_fortran_env
 implicit none
!--------------------------------------------
 integer, parameter       :: dp = real64
 integer                  :: context
 integer                  ::  nprow,  npcol, nprocs
 integer                  :: myprow, mypcol, myrank
 integer                  :: myimg , nimg  , myi   , myj
!--------------------------------------------
 integer, parameter       :: DLEN = 9        ! matrix descriptor size
 integer, dimension(DLEN) :: desca           ! matrix descriptors
 integer, external        :: numroc          ! must provide definitions or interfaces to BLACS functions
 integer                  :: info
!--------------------------------------------
 integer                  :: N,MB,NB,Mloc,Nloc
 integer                  :: i,j ,iloc,jloc
 real(dp), allocatable    :: A    (:,:)[:,:]
 real(dp), allocatable    :: Aglob(:,:)
!--------------------------------------------
 character(100)           :: fname
 integer                  :: fout
!--------------------------------------------
 call init
!--------------------------------------------
 N    = 16
 MB   = N/nprow ; NB = N/npcol           ! Fit one block per image/process
 Mloc = numroc(N,MB,myprow,0,nprow) 
 Nloc = numroc(N,NB,mypcol,0,npcol) 
 if(Mloc /= N/nprow) error stop
 if(Nloc /= N/npcol) error stop
 if(Nloc /= Mloc   ) error stop
!Array descriptor:
 call descinit(desca,N,N,MB,NB,0,0,context,Mloc,info)           !DESCINIT( DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT,LLD, INFO )
!--------------------------------------------
 allocate(A(Nloc,Nloc)[nprow,*])
 myi = this_image(A,1); myj = this_image(A,2)
 write(fout,'(A    )') '00  myimg myrank myprow mypcol     MB  nprow'
 write(fout,'(A,6I7)') '00',myimg,myrank,myprow,mypcol,    MB, nprow
!--------------------------------------------
!Define the local arrays: The global array is Aglob(i,j) = i + j /100
 A    = 0.0_dp
 do i = 1, N ; do j = 1, N                        
  call pdelset(A,i,j,desca, i + j / 100.0_dp)       ! Parallel Double ELement SET = PDELSET  (expensive...)
 end do      ; end do
!--------------------------------------------
!Print local arrays:
 do iloc = 1, Nloc
  write(fout,'(A,1000F7.2)') '01',(A(iloc,jloc),jloc=1,Nloc)
 end do
!--------------------------------------------
!Construct global array:
 allocate(Aglob(N,N))
 Aglob=0.0_dp
 do i = 1, N; do j = 1, N                           ! Parallel Double ELement GET
  call pdelget('A',' ',Aglob(i,j),A,i,j,desca)      ! SUBROUTINE PDELGET( SCOPE, TOP, ALPHA, A, IA, JA, DESCA )
 end do     ; end do

 write(fout,'(A)')'------------------------------'
 do i = 1, N
  write(fout,'(A,1000F7.2)') '02',(Aglob(i,j),j=1,N)
 end do

contains
!--------------------------------------------
 subroutine     init
  
 call sca_init

 write(fname,'(A,I0.3)') 'out',myimg
 open(newunit=fout,file=fname)

 end subroutine init
!--------------------------------------------
 subroutine     sca_init
  call blacs_pinfo(myrank,nprocs)
  nprow = sqrt(nprocs+1.0e-6_dp); npcol = nprow
  if(nprocs /= nprow * npcol) error stop
  call blacs_get(-1,0, context)
  call blacs_gridinit(context,'Col-major', nprow, npcol)
  call blacs_gridinfo(context,nprow,npcol,myprow,mypcol)
  myimg = this_image()
  nimg  = num_images()
 end subroutine sca_init
!--------------------------------------------
end program test_scalapack
!--------------------------------------------
