module blacs
! store blacs context and processor data
! Time-stamp: "2003-09-19 10:33:48 cjn"

  public

  integer, save    :: ictxt  ! blacs default context handle

  integer, save    :: ctxt   ! blacs context handle
  integer, save    :: iam    ! processor id
  logical, save    :: io_processor   ! i/o processor flag
  integer, save    :: p, q   ! # processors in grid
  integer, save    :: mycol, myrow ! processor coordinates
  integer, save    :: p0, q0 ! i/o processor coordinates
  integer, save    :: nprocs ! # processors
  integer, save    :: nblock ! blacs blocking factor
  integer, save    :: mynumrows
  integer, save    :: mynumcols
end module blacs

module pdg_ctl
! Time-stamp: "2003-09-17 17:00:13 cjn"
  implicit none
  integer, save  :: call_no = 0

  private
  public pdiag

contains

  subroutine pdiag (imesh, ctxti, myrowi, mycoli, pp, qq)
   use blacs, only: ctxt, p, q, myrow, mycol, iam, ictxt
   integer, intent(in) :: imesh, ctxti, myrowi, mycoli
   integer, intent(in) :: pp, qq
   integer             :: r, c, pnum, BLACS_pnum
   real(8) :: y(10000), x, z(50000)

   call_no = call_no + 1

   ctxt = ctxti  ! set ctxt to default context handle
   p = pp
   q = qq
   myrow = myrowi
   mycol = mycoli

   if (imesh == 2) then ! bulid in a delay
      call random_number (y)
      x = SUM(EXP(-y))
      write (*,*) 'imesh, x (sum y) = ', imesh, x
   else if (imesh == 3) then
      call random_number (z)
      x = SUM(EXP(-z))
      write (*,*) 'imesh, x (sum-z)= ', imesh, x
   end if

   pnum = BLACS_PNUM (ctxt, myrow, mycol)
   call BLACS_PCOORD (ctxt, pnum, r, c)
   write (*,*) '->pdiag: iam, myrow, mycol, r, c = ', iam, myrow, mycol, r, c
 end subroutine pdiag
end module pdg_ctl

program rmx1
! Time-stamp: "2003-09-17 14:39:03 cjn"
  use pdg_ctl, only: pdiag
  use blacs, only: nprocs, iam, ictxt, p, q
  implicit none
  integer, parameter          :: wp = selected_real_kind(12)
  integer, parameter          :: fo = 6
  integer                     :: i, j, ion, ncp, imesh, ibuf(2)
  integer                     :: st1, st2
  real(wp)                    :: t0, t1
  integer                     :: c0, c1, cr
  integer                     :: pin(4), qin(4)
  integer                     :: myrow, mycol, pp, qq, k, iproc, np, nq, npr
  integer                     :: blacs_pnum, ctxt(4)
  type a
     integer, pointer :: map(:,:)
  end type a
  type(a)                    :: maps(4)

  call cpu_time (t0)
  call system_clock (count=c0)

  call BLACS_PINFO (iam, nprocs) ! find process #, total # processors

  if (iam == 0) then ! This is the i/o processor
     write (fo,'(//,15x,a)') '============='
     write (fo,'(15x,a)') 'Program test_blacs'
     write (fo,'(15x,a,//)') '============='
     write (fo,'(a,i6)') 'Number of processors = ', nprocs
  end if

! define blacs grid to include all processes to broadcast information
  call BLACS_GET (-1, 0, ictxt)  ! find default context, ictxt
  call BLACS_GRIDINIT (ictxt, 'Row-major', 1, nprocs)

   if (iam == 0) then
      pin = (/1, 1, 1, 1/)
      qin = (/2, 2, 2, 2/)
      call igebs2d (ictxt, 'all', ' ', 4, 1, pin, 4)
      call igebs2d (ictxt, 'all', ' ', 4, 1, qin, 4)
    else
       call igebr2d (ictxt, 'all', ' ', 4, 1, pin, 4, 0, 0)
       call igebr2d (ictxt, 'all', ' ', 4, 1, qin, 4, 0, 0)
    end if

    call BLACS_BARRIER (ictxt, 'A')

    iproc = -1
    contexts: do i = 1, 4
       np = pin(i)
       nq = qin(i)
       allocate (maps(i)%map(np,nq))
       do k = 1, nq
          do j = 1, np
             iproc = iproc + 1
             if (iproc >= nprocs) call BLACS_EXIT()
             maps(i)%map(j,k) = BLACS_PNUM(ictxt, 0, iproc)
          end do
       end do
       ctxt(i) = ictxt
       call BLACS_GRIDMAP (ctxt(i), maps(i)%map, np, np, nq)
    end do contexts
    write (*,*) 'contexts done, iam = ', iam, ctxt

    mesh_loop: do i = 1, 4
       call BLACS_GRIDINFO (ctxt(i), pp, qq, myrow, mycol)
       if (myrow < pp .and. mycol < qq) then ! in this grid
          write (*,*) 'gridinfo: sdelected iam, myrow, mycol, p, q = ', iam, &
               myrow, mycol, pp, qq
          call pdiag (i, ctxt(i), myrow, mycol, pp, qq)
          call BLACS_GRIDEXIT (ctxt(i))    ! release context
       end if
    end do mesh_loop

    call blacs_barrier (ictxt, 'All')

  if (iam == 0) then
     write (fo,'(/,a,/)') 'end of RMX1'
     call cpu_time (t1)
     write (fo,'(a,f16.4,a)') 'CPU time     = ', t1 - t0, ' secs'
     call system_clock (count=c1, count_rate=cr)
     write (fo,'(a,f16.4,a)') 'Elapsed time = ', REAL(c1-c0,wp) / &
          REAL(cr,wp), ' secs'
  end if

  call BLACS_GRIDEXIT (ictxt)
  call BLACS_EXIT()
  stop
end program rmx1
