pfarm/cpu/src_mpi_omp_atom/rmx1_in.f90

+module rmX1_in
+! read data defining pdiag run
+! Time-stamp: "2005-07-28 09:55:40 cjn"
+! modified for jpi_coupling april to dec 2006 vmb
+  use precisn, only: wp
+  use io_units, only: nfnml, fo
+  use mpi_params, only: io_processor
+!  use blacs, only: p, q, nblock, ictxt, io_processor
+  implicit none
+
+  integer, save            :: bug1, bug2, bug3, bug4, bug5
+  integer, save            :: bug6, bug7, bug8, bug9
+  integer, save            :: lrgl1, nspn1, npty1
+  integer, save            :: num_timers_dig
+  integer, save            :: p,q, nblock
+  integer, save            :: n_lambda, nhmax, ncmax
+  real(wp), save           :: abvthr, belthr, degeny, rafin
+  integer, save            :: ne(10)
+  real(wp), save           :: esc(21)
+  real(wp), save           :: einc_fine, einc_coarse
+  real(wp), save           :: emax
+  logical, save            :: split_prop
+  real(wp), save           :: tdfn(50)
+  real(wp), pointer, save  :: tconst(:)
+  character(len=20), save  :: flabel_dig
+  character(len=80),save   :: title
+  character(len=120), save :: filh, fltm
+  integer, save            :: dflag(5), pflg(5)
+  integer, save            :: targ_low_sce, targ_high_sce
+  integer, save            :: nz, nelc    
+  logical, save            :: inc_lrp_prop
+  logical, save            :: xdr
+  logical, save            :: uniform_e_grid
+  logical, save            :: jpi_coupling
+  integer, save            :: start_mesh
+
+  private
+  public rd_nmlst
+  public lrgl1, nspn1, npty1, n_lambda, nhmax, ncmax
+  public abvthr, belthr, degeny, rafin, ne, esc, einc_fine
+  public einc_coarse, emax, tdfn, tconst, flabel_dig, title, filh
+  public fltm, dflag, pflg, targ_low_sce, targ_high_sce, start_mesh
+  public inc_lrp_prop, xdr, split_prop, uniform_e_grid, jpi_coupling
+  public bug1, bug2, bug3, bug4, bug5, bug6, bug7, bug8, bug9
+
+contains
+
+  subroutine rd_nmlst
+! parameters required only in the propagation stage:
+    integer                  :: num_rm_gen, num_asy_per_pipe
+    integer                  :: ne_write, nc2, ng, nxsn_store
+    integer                  :: time_limit, num_timers_prop
+    integer                  :: nts_save, nt
+    integer                  :: xsn_list(100)
+    real(wp)                 :: ewron, eps
+    logical                  :: buttle
+    logical                  :: gail_exp
+    character(len=20)        :: flabel_prop
+    integer                  :: ios, ctxt, debug(9)
+
+! Input variables used in diag stage:
+    namelist /phzin/ title, dflag, pflg, einc_fine, nz, nelc,    &
+         einc_coarse, emax, targ_low_sce, targ_high_sce,         &
+         abvthr, belthr, degeny, nhmax, ncmax, filh, start_mesh, &
+         fltm, nspn1, npty1, lrgl1, n_lambda, uniform_e_grid,    &
+         split_prop, inc_lrp_prop, rafin, num_timers_dig,        &
+         ne, esc, flabel_dig, p, q, nblock, xdr, debug,          & 
+         jpi_coupling,                                           &           
+! Input variables used in prop stage                             
+         num_rm_gen, buttle, nxsn_store, xsn_list,               &
+         ne_write, nts_save, time_limit, eps, ng, nt,            &
+         tdfn, gail_exp, num_asy_per_pipe,                       &
+         num_timers_prop, flabel_prop
+
+! default values
+    split_prop = .true.
+    inc_lrp_prop = .true.
+    targ_low_sce = 0
+    targ_high_sce = 0
+    einc_fine = 0.00002_wp
+    einc_coarse = 0.01_wp
+    emax = 1.86_wp
+    nhmax = 8000
+    ncmax = 315
+    ne = 0
+    esc= 0.0_wp
+    title = REPEAT(' ',80)
+    filh  = REPEAT(' ',120)
+    filh = 'H'
+    fltm = 'TARMOM'
+    abvthr = 1.0e-3_wp
+    belthr = 1.0e-3_wp
+    degeny = 1.0e-06_wp
+    npty1 = -1 !L pi must be input
+    rafin = -1.0_wp
+    n_lambda = 2
+    flabel_dig = 'DIG'
+    num_timers_dig = 4
+    xdr = .true.
+    debug = 0
+    uniform_e_grid = .false.
+    start_mesh = 1
+    jpi_coupling = .false. !default to LS coupling
+
+    open (unit=nfnml, file = 'phzin.ctl', status='old', &
+         position = 'rewind', iostat= ios)
+    if (ios /= 0) then
+       write (fo,'(a)') " phzin.ctl file not found - aborting "
+       STOP
+    end if
+    read (nfnml, phzin, iostat=ios)
+    if (ios /= 0) then
+       write (fo,'(a,i6)') " Problem with Namelist file - aborting ",ios
+       write (fo, phzin)
+       STOP
+    end if
+    close (nfnml)
+
+    bug1 = debug(1); bug2 = debug(2); bug3 = debug(3); bug4 = debug(4)
+    bug5 = debug(5); bug6 = debug(6); bug7 = debug(7); bug8 = debug(8)
+    bug9 = debug(9)
+
+! Rules set by the rmatr2/95 interface and enforced in rmx and prm
+!  Exchange runs:
+!   In LS coupling:
+!    lrgl1 = total angular momentum L 
+!    nspn1 = total multiplicity S
+!   In jpi coupling:
+!    lrgl1 = 2*J
+!    nspn1 = 0
+!  No_exchange runs:
+!   In LS coupling:
+!    lrgl1 = total angular momentum L 
+!    nspn1 = -target multiplicity Si
+!   In jpi coupling:
+!    lrgl1 = -2*K
+!    nspn1 = 1 weighting factor for lowest K
+!          = 2 weighting factor for other K
+
+    if (jpi_coupling) nspn1=0 !until specified by interface data
+
+! Test for no-exchange files
+    if (jpi_coupling) then
+       if (MOD(nelc+lrgl1,2)==0) lrgl1=-ABS(lrgl1) !lrgl1 has been input as 2*K 
+    else
+       if (MOD(nelc+nspn1,2)==1) nspn1=-ABS(nspn1) !nspn1 is target multiplicity
+    end if
+! print basic input data or complete namelist phzin if diaflg(1) set
+
+    if(io_processor) then
+      call prt_title
+      write (fo,'(a,2es14.6)') 'Threshold zone, abvthr, belthr = ', &
+         abvthr, belthr
+      write (fo,'(a)') 'H-file: ' // TRIM(filh)
+      write (fo,'(a)') 'Target Moment File: ' // TRIM(fltm)
+      write (fo,'(/)')
+      if (dflag(2) > 0) write (fo,phzin)
+    end if
+  end subroutine rd_nmlst
+
+  subroutine prt_title
+! write a general title indicating the target ion/atom,
+! # continuum electrons, highest electron energy
+    character(len=2), parameter    :: atom(57) = (/ &
+         ' H', 'He', 'Li', 'Be', ' B', ' C', ' N', ' O', ' F', &
+         'Ne', 'Na', 'Mg', 'Al', 'Si', ' P', ' S', 'Cl', 'Ar', &
+         ' K', 'Ca', 'Sc', 'Ti', ' V', 'Cr', 'Mn', 'Fe', 'Co', &
+         'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', &
+         'Rb', 'Sr', ' Y', 'Zr', 'Nb', 'Mo', 'Tc', 'Ru', 'Rh', &
+         'Pd', 'Ag', 'Cd', 'In', 'Sn', 'Sb', 'Te', ' I', 'Xe', &
+         'Cs', 'Ba', 'La'/)
+    character(len=1)                :: chrge
+    character(len=8)                :: today
+    character(len=10)               :: now
+    integer                         :: nzeff
+
+    call date_and_time (today, now)
+! today: ccyymmdd
+    write (fo,'(a,a2,a,a2,a,a4)') 'Job run on ', today(7:8), &
+         '/', today(5:6), '/', today(1:4)
+! now: hhmmss.sss
+    write (fo,'(a,a2,a,a2,a,a5)') 'Time: ', now(1:2), ':', &
+         now(3:4), ':', now(5:10)
+
+    nzeff = nz - nelc
+    if (nzeff > 0) then
+       chrge = '+'
+    else if (nzeff < 0) then
+       chrge = '-'
+       nzeff = - nzeff
+    end if
+
+    if (nzeff > 1) then
+       write (fo,'(5x,37("*"))')
+       write (fo,'(5x,"*",35x,"*")')
+       if (nzeff < 10) then
+          write (fo,'(5x,"*",29x,i1,a1,4x,"*")') nzeff, chrge
+       else
+          write (fo,'(5x,"*",29x,i2,a1,3x,"*")') nzeff, chrge
+       end if
+
+       if (nz > 57) then
+          write (fo,'(5x,"*",4x,a,i4,6x,"*")') 'Electron scatter&
+               &ing by Z = ', nz
+       else
+          write (fo,'(5x,"*",4x,a,a2,6x,"*")') 'Electron scatter&
+               &ing by ', atom(nz)
+       end if
+       write (fo,'(5x,"*",35x,"*")')
+       write (fo,'(5x,37("*"))')
+    else
+       write (fo,'(/,5x,37("*"))')
+       write (fo,'(5x,"*",35x,"*")')
+       if (nzeff == 1) write (fo,'(5x,"*",29x,a1,5x,"*")') chrge
+       if (nz > 57) then
+          write (fo,'(5x,"*",4x,a,i4,6x,"*")') 'Electron scatter&
+               &ing by Z = ', nz
+       else
+          write (fo,'(5x,"*",4x,a,a2,6x,"*")') 'Electron scattering by ',&
+               atom(nz)
+       end if
+       write (fo,'(5x,"*",35x,"*")')
+       write (fo,'(5x,37("*"),/)')
+    end if
+    if (title /= ' ') write (fo,'(/,a,/)') TRIM(title)
+  end subroutine prt_title
+
+end module rmX1_in

pfarm/cpu/src_mpi_omp_atom/scaling.f90

+module scaling
+! Time-stamp: "2003-09-18 07:51:08 cjn"
+  use precisn, only: wp
+  use io_units, only: fo
+  implicit none
+
+  private
+  public set_charge, charge, get_charge, scale_etarg
+  public scale_cf, scale_radius, get_ionicity
+
+  real(wp), parameter         :: atry = 2.0_wp
+  real(wp), save              :: iz_res = -999
+  logical, save               :: charged_case
+  integer, save               :: ion  ! ionicity
+
+contains
+
+  subroutine set_charge (nz, nelc)
+! define overall Coulomb charge on target
+    integer, intent(in)      :: nz   ! nuclear charge
+    integer, intent(in)      :: nelc ! # target electrons
+
+    iz_res = nz - nelc
+    if (iz_res == 0) then
+       charged_case = .false.
+       ion = 0
+    else
+       charged_case = .true.
+       if (iz_res > 0) then
+          ion = 1
+       else
+          ion = iz_res
+       end if
+    end if
+  end subroutine set_charge
+
+  function get_ionicity () 
+    integer :: get_ionicity
+    get_ionicity = ion
+  end function get_ionicity
+
+  subroutine print_charge
+    write (fo,'(a,i6)') 'Residual charge (iz_res) = ', iz_res
+  end subroutine print_charge
+
+  subroutine get_charge (ires_charge)
+    integer, intent(out)   :: ires_charge
+
+    if (iz_res == -999) then ! not initialized - warn
+       write(fo,'(a)') 'get_charge: attempt to access uninitialized &
+            &charge'
+    else
+      ires_charge = iz_res
+   end if
+  end subroutine get_charge
+
+  subroutine charge (charge_log)
+    logical, intent(inout)  :: charge_log
+
+    charge_log = charged_case
+  end subroutine charge
+
+  subroutine scale_radius (rmatr)
+    real(wp), intent(inout)   :: rmatr
+
+    rmatr = REAL(iz_res) * rmatr
+  end subroutine scale_radius
+
+  subroutine scale_etarg (etarg)
+    real(wp),intent(inout)    :: etarg(:)
+    real(wp)                  :: fac
+
+    if (charged_case) then
+       fac = atry / REAL(iz_res*iz_res,wp)
+       etarg = fac * etarg
+    else
+       fac = atry
+       etarg = etarg * fac
+    end if
+  end subroutine scale_etarg
+
+  subroutine scale_cf (cf,lamd)
+    real(wp), intent(inout)    :: cf(:,:)
+    integer, intent(in)        :: lamd(:,:)
+    real(wp)                   :: flion
+    integer                    :: i,j,iplim,jplim
+
+    iplim = size(cf,dim=1)
+    jplim = size(cf,dim=2)
+    flion = REAL(iz_res,wp)
+    do j = 1, jplim
+       do i = 1, iplim
+          cf(i,j) = atry * cf(i,j) * flion**(lamd(i,j)-1)
+       end do
+    end do
+  end subroutine scale_cf
+
+end module scaling

pfarm/cpu/src_mpi_omp_atom/sec_ham.f90

+module sec_ham
+! generates sector Hamiltonians
+! Time-stamp: "2003-10-03 17:36:23 cjn"
+  use precisn, only: wp
+  use io_units, only: fo
+  use error_out, only: error_check
+
+  implicit none
+
+  integer, save               :: nx       ! quadradure order
+  real(wp), allocatable, save :: xi(:)    ! quadrature nodes
+  real(wp), allocatable, save :: wi(:)    ! quadrature weights
+  real(wp), allocatable, save :: pl(:,:)  ! normalized Legendre polynomials
+  real(wp), allocatable, save :: vbl(:)   ! sector left bdry amplitudes
+  real(wp), allocatable, save :: vbr(:)   ! sector right bdry amplitudesx
+
+  real(wp), save              :: r12, rsq, r1b, r2b
+
+  integer, save               :: ncq
+
+  private
+  public def_rb_vals, h_el, gauleg, h_reset, vbl, vbr, legndr, xi
+  public def_ncq
+
+contains
+
+  subroutine def_ncq (ncp)
+    integer, intent(in)   :: ncp   ! spin partition offset
+    ncq = ncp - 1
+  end subroutine def_ncq
+
+  subroutine def_rb_vals (ral, rar)
+! define sector geometry parameters
+    use hfile_data, only: bbloch
+    real(wp), intent(in)    :: ral   ! left endpoint of sector
+    real(wp), intent(in)    :: rar   ! right endpoint of sector
+
+    r12 = 1.0_wp / (rar - ral)
+    if (bbloch == 0.0_wp) then
+       r1b = 0.0_wp
+       r2b = 0.0_wp
+    else
+       r1b = r12 * bbloch / ral
+       r2b = r12 * bbloch / rar
+    end if
+    rsq = 2.0_wp * r12**2
+  end subroutine def_rb_vals
+
+  function h_el (iglob, jglob)
+! Hamiltonian matrix element for global indices iglob, jglob
+! should be called for iglob <= jglob
+    use sector_potl, only: nc, v
+    use hfile_data, only: ethr
+    use def_sectors, only: nl
+    real(wp)              :: h_el
+    integer, intent(in)   :: iglob, jglob  ! coordinates of element
+    integer               :: i, j, in, jn, ic, jc
+    real(wp)              :: rmn, ss
+    integer               :: ij
+
+    j = jglob
+    jc = 1 + (j-1)/nl     ! channel index
+    jn = j - (jc-1)*nl    ! basis index (inner loop)
+    i = iglob
+    ic = 1 + (i-1)/nl     ! channel index
+    in = i - (ic-1)*nl    ! basis index
+
+    ss = 0.0_wp
+    if (ic == jc) then
+       rmn = SQRT(REAL((2*jn-1)*(2*in-1),wp))
+       ss = - rmn * (r2b - r1b * (-1)**MOD(in+jn,2))
+       if (in == jn) ss = ss + ethr(ncq+jc)
+       if (MOD(in+jn,2) == 0) ss = ss + rsq * rmn * REAL(in*(in-1),wp)
+    end if
+    h_el = ss + SUM(pl(:,in) * wi(:) * pl(:,jn) * v(:,ic,jc))
+  end function h_el
+
+  subroutine legndr
+! Legndr: normalized legendre polynomials at quadrature points
+    use def_sectors, only: nl, nx, w_sect
+    real(wp)               :: sgn
+    integer                :: k, status
+
+    allocate (pl(nx,nl), vbl(nl), vbr(nl), stat=status)
+    call error_check (status, 'legndr: allocation')
+
+    pl(:,1) = SQRT(0.5_wp)
+    vbr(1) = SQRT(1.0_wp/w_sect)
+    vbl(1) = vbr(1)
+    if (nl == 1) return
+    pl(:,2) = SQRT(1.5_wp) * xi
+    vbr(2) = SQRT(3.0_wp/w_sect)
+    vbl(2) = - vbr(2)
+    sgn = -1.0_wp
+    do k = 3, nl
+       pl(:,k) = (SQRT(REAL(2*k-3,wp)) * xi * pl(:,k-1) - REAL(k-2,wp) * &
+            pl(:,k-2) / SQRT(REAL(2*k-5,wp))) * SQRT(REAL(2*k-1,wp)) / &
+            REAL(k-1,wp)
+       vbr(k) = SQRT(REAL(2*k-1,wp)/w_sect)
+       sgn = - sgn
+       vbl(k) = sgn * vbr(k)
+    end do
+  end subroutine legndr
+
+  subroutine gauleg
+!  Gauss-Legendre quadrature abscissas and weights
+    use def_sectors, only: nx
+    integer                 :: i, ierr
+    real(wp)                :: b(nx-1), d(nx), e(nx)     ! automatic
+    real(wp)                :: z(nx,nx)                 ! automatic
+    real(wp)                :: work(MAX(1,2*nx-2))     ! automatic
+    integer                 :: status
+
+    allocate (xi(nx), wi(nx), stat=status)
+    call error_check (status, 'gauleg: allocation')
+
+    do i = 1, nx-1
+       b(i) = REAL(i*i, wp) / REAL((2*i+1) * (2*i-1), wp)
+    end do
+    d = 0.0_wp                  ! diagonal elements
+    e(:nx-1) = - SQRT(b(1:nx-1))   ! sub-diagonal elements
+    call dsteqr ('i', nx, d, e(:nx-1), z, nx, work, ierr) ! diagonalization
+    if (ierr /= 0) then
+       write (fo,'(a,i6)') 'GAULEG: error code from dsteqr:', ierr
+    end if
+    e = 2.0_wp * z(1,:)**2
+    xi = d
+    wi = e
+  end subroutine gauleg
+
+  subroutine h_reset
+! deallocate memory
+    integer    :: status
+    deallocate (xi, wi, vbr, vbl, pl, stat=status)
+    call error_check (status, 'h_reset: deallocation')
+  end subroutine h_reset
+end module sec_ham

pfarm/cpu/src_mpi_omp_atom/slp.f90

+module slp
+! LS quantum numbers
+! Time-stamp: "03/07/04 14:33:29 cjn"
+  implicit none
+
+  private
+  public qno, eq_qno, def_qno, val_qno
+
+  type qno
+     integer     lrgl     ! orbital am q. #
+     integer     nspn     ! spin multiplicity
+     integer     npty     ! parity q. #
+  end type qno
+
+  interface operator (==)
+     module procedure eq_qno
+  end interface
+
+contains
+
+  function eq_qno (q1, q2)
+! define equivalence for q # type
+    logical               :: eq_qno
+    type(qno), intent(in) :: q1
+    type(qno), intent(in) :: q2
+    eq_qno = .true.
+    if (q1%lrgl /= q2%lrgl .or. q1%nspn /= q2%nspn .or. &
+         q1%npty /= q2%npty) eq_qno = .false.
+  end function eq_qno
+
+  function def_qno (l, s, p)
+! define q. # type
+    type(qno)              :: def_qno
+    integer, intent(in)    :: l   ! orbital am
+    integer, intent(in)    :: s   ! spin multiplicity
+    integer, intent(in)    :: p   ! parity
+
+    def_qno%lrgl = l
+    def_qno%nspn = s
+    def_qno%npty = p
+  end function def_qno
+
+  function val_qno (q)
+! return values of the q. #s in q
+    integer                :: val_qno(3)
+    type(qno), intent(in)  :: q
+
+    val_qno(1) = q%lrgl
+    val_qno(2) = q%nspn
+    val_qno(3) = q%npty
+  end function val_qno
+
+end module slp

pfarm/cpu/src_mpi_omp_atom/test_blacs.f90

+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

pfarm/cpu/src_mpi_omp_atom/xdr_dummy_files.f90

+module xdr_files
+! f90 dummy interface to XDR routines
+! used when library libfxdr.a is not available
+! all calls result in noops except open_xdr which terminates
+! Time-stamp: "03/03/18 08:04:42 cjn"
+  use io_units, only: fo
+  implicit none
+
+  integer, save   :: hbufl = 2048   ! xdr buffer length
+
+  private
+  public xdr_io, open_xdr, close_xdr, rewind_xdr
+  public hbufl
+
+  interface xdr_io
+     module procedure xdr_int, xdr_imt, xdr_real
+     module procedure xdr_rlmt, xdr_chrs
+  end interface
+
+contains
+   function open_xdr (file, action)
+! open xdr file, fortran open syntax style
+    integer                          :: open_xdr
+    character(len=*), intent(in)     :: file
+    character(len=*), intent(in)     :: action
+
+    open_xdr = 0
+    if (open_xdr == 0) then
+       write (fo,'(a)') 'open_xdr called: libxdr.a is not available'
+       stop
+    end if
+  end function open_xdr
+
+  subroutine close_xdr (xdr)
+    integer, intent(in)     :: xdr
+  end subroutine close_xdr
+
+  subroutine rewind_xdr (xdr)
+    integer, intent(in)     :: xdr
+  end subroutine rewind_xdr
+
+  subroutine xdr_int (xdr, ival)
+    integer, intent(in)      ::  xdr
+    integer, intent(inout)   ::  ival
+  end subroutine xdr_int
+
+  subroutine xdr_imt (xdr, ival, nels)
+    integer, intent(in)            :: xdr
+    integer, intent(in)            :: nels
+    integer, intent(inout)         :: ival(:)
+  end subroutine xdr_imt
+
+  subroutine xdr_real (xdr, dval)
+    use precisn, only: wp
+    integer, intent(in)      :: xdr
+    real(wp), intent(inout)   :: dval
+  end subroutine xdr_real
+
+  subroutine xdr_rlmt (xdr, dval, nels)
+    use precisn, only: wp
+    integer, intent(in)           :: xdr
+    integer, intent(in)           :: nels
+    real(wp), intent(inout)      :: dval(nels)
+  end subroutine xdr_rlmt
+
+  subroutine xdr_chrs (xdr, string)
+    integer, intent(in)            :: xdr
+    character(len=*), intent(inout) :: string
+  end subroutine xdr_chrs
+
+end module xdr_files

pfarm/cpu/src_mpi_omp_atom/xdr_files.f90

+module xdr_files
+! f90 interface to XDR routines
+! uses library libfxdr.a
+! Time-stamp: "03/03/18 08:06:30 cjn"
+  use io_units, only: fo
+  implicit none
+
+  interface xdr_io
+     module procedure xdr_int, xdr_imt, xdr_real
+     module procedure xdr_rlmt, xdr_chrs
+  end interface
+
+  logical, save :: real64=.false.
+  integer, save :: hbufl = 2048 ! xdr file buffer length
+
+  private
+  public xdr_io, open_xdr, close_xdr, rewind_xdr, hbufl
+
+contains
+   function open_xdr (file, action)
+! open xdr file, fortran open syntax style
+    integer                          :: open_xdr
+    character(len=*), intent(in)     :: file
+    character(len=*), intent(in)     :: action
+    integer                          :: xdr_open, imode, lenf
+
+    if (action(:1) == 'r' .or. action(:1) == 'R') then
+       imode = 1
+    elseif (action(:1) == 'w' .or. action(:1) == 'W') then
+       imode = 2
+    else
+       write (fo,'(a,a5)') 'open_xdr error: unknown action =', &
+            TRIM(action)
+       stop
+    endif
+
+    lenf = LEN(TRIM(file))
+    open_xdr = xdr_open (lenf, file, imode) + 1
+  end function open_xdr
+
+  subroutine close_xdr (xdr)
+    integer, intent(in)     :: xdr
+
+    call xdr_close(xdr-1)
+  end subroutine close_xdr
+
+  subroutine rewind_xdr (xdr)
+    integer, intent(in)     :: xdr
+
+    call xdr_rewind (xdr-1)
+  end subroutine rewind_xdr
+
+  subroutine xdr_int (xdr, ival)
+    integer, intent(in)      ::  xdr
+    integer, intent(inout)   ::  ival
+
+    call xdr_fint (xdr-1, ival)
+  end subroutine xdr_int
+
+  subroutine xdr_imt (xdr, ival, nels)
+    integer, intent(in)            :: xdr
+    integer, intent(in)            :: nels
+    integer, intent(inout)         :: ival(:)
+
+    call xdr_imat (xdr-1, nels, ival)
+  end subroutine xdr_imt
+
+  subroutine xdr_real (xdr, dval)
+    use precisn, only: wp
+    integer, intent(in)      :: xdr
+    real(wp), intent(inout)   :: dval
+
+    if (real64) then
+       call xdr_real64 (xdr-1, dval)
+    else
+       call xdr_fdouble (xdr-1, dval)
+    endif
+  end subroutine xdr_real
+
+  subroutine xdr_rlmt (xdr, dval, nels)
+    use precisn, only: wp
+    integer, intent(in)           :: xdr
+    integer, intent(in)           :: nels
+    real(wp), intent(inout)        :: dval(nels)
+
+    if (real64) then
+       call xdr_rmat64 (xdr-1, nels, dval)
+    else
+       call xdr_dmat (xdr-1, nels, dval)
+    endif
+  end subroutine xdr_rlmt
+
+  subroutine xdr_chrs (xdr, string)
+    integer, intent(in)         :: xdr
+    character(*), intent(inout) :: string
+
+    call xdr_fstring (xdr-1, LEN(string), string)
+  end subroutine xdr_chrs
+
+end module xdr_files

pfarm/cpu/src_mpi_omp_mol/.#pdg_ctl.f90

+ags@eslogin006.20921:1512035368
\ No newline at end of file

pfarm/cpu/src_mpi_omp_mol/DEFS

+DEFS_cray
\ No newline at end of file

pfarm/cpu/src_mpi_omp_mol/DEFS_cray

+LIBS = 
+STATIC_LIBS = -L../../lib -lfxdr
+
+# Default Archer2 environment (Cray) 
+F90 = ftn -Ofast
+
+# Skylake options
+#F90FLAGS = -mtune=skylake -mkl=parallel -O3
+
+# AMD options
+F90FLAGS = -O3
+
+F90FLAGS1 = $(F90FLAGS)
+LDFLAGS=$(F90FLAGS)

pfarm/cpu/src_mpi_omp_mol/DEFS_cray~

+LIBS = 
+STATIC_LIBS = -L../../lib -lfxdr
+
+# Default Archer2 environment (Cray) 
+F90 = ftn
+
+# Skylake options
+#F90FLAGS = -mtune=skylake -mkl=parallel -O3
+
+# AMD options
+F90FLAGS = -O3
+
+F90FLAGS1 = $(F90FLAGS)
+LDFLAGS=$(F90FLAGS)

pfarm/cpu/src_mpi_omp_mol/DEFS_skl

+LIBS = 
+STATIC_LIBS = -L../lib -lfxdr
+
+F90 = mpif90
+
+# KNL options
+#F90FLAGS = -mtune=knl -mkl=parallel -O3
+
+# Skylake options
+F90FLAGS = -mtune=skylake -mkl=parallel -Ofast
+
+F90FLAGS1 = $(F90FLAGS)
+LDFLAGS=$(F90FLAGS)

pfarm/cpu/src_mpi_omp_mol/Makefile

+PROG = ../bin/exdig_mpi_omp_mol
+
+include DEFS
+
+XDR = xdr_files
+
+SRCS =	amp.f90 def_sectors_bench.f90 pdg_ctl.f90 dist_mat.f90 \
+	energy_grid.f90 rmx1.f90 potl.f90 scaling.f90 \
+	potl_coefs.f90 precisn.f90 rd_Hfile.f90 rmx1_in.f90 \
+	sec_ham.f90 slp.f90 error_out.f90 io_units.f90 mpi_params.f90
+SRC1 = $(XDR).f90
+
+
+OBJS = $(SRCS:.f90=.o)
+OBJ1 = $(SRC1:.f90=.o)
+
+
+all: $(PROG)
+
+$(PROG): $(OBJ1) $(OBJS) $(ULIBS)
+	$(F90) $(LDFLAGS) -o $@ $(OBJS) $(OBJ1) $(LIBS) $(STATIC_LIBS)
+
+$(OBJ1): $(SRC1)
+	$(F90) $(F90FLAGS) -c  $<
+
+clean:
+	rm -f $(PROG) $(OBJS) $(OBJ1) *.mod *~
+very_clean:
+	rm -f $(PROG) $(OBJS) $(OBJ1) *.mod *~
+
+.SUFFIXES: $(SUFFIXES) .f90
+
+.f90.o:
+	$(F90) $(F90FLAGS) -c $<
+
+amp.o: amp.f90 def_sectors_bench.o dist_mat.o energy_grid.o \
+	io_units.o precisn.o slp.o $(XDR).o error_out.o mpi_params.o
+def_sectors_bench.o: def_sectors_bench.f90 energy_grid.o io_units.o precisn.o \
+	rd_Hfile.o rmx1_in.o error_out.o mpi_params.o
+pdg_ctl.o: pdg_ctl.f90 amp.o def_sectors_bench.o dist_mat.o io_units.o potl.o \
+	potl_coefs.o precisn.o rd_Hfile.o rmx1_in.o slp.o mpi_params.o
+dist_mat.o: dist_mat.f90 io_units.o precisn.o sec_ham.o error_out.o
+energy_grid.o: energy_grid.f90 io_units.o precisn.o rd_Hfile.o \
+	rmx1_in.o error_out.o mpi_params.o
+error_prt.o: error_prt.f90 io_units.o
+rmx1.o: rmx1.f90 pdg_ctl.o energy_grid.o io_units.o precisn.o rd_Hfile.o \
+	rmx1_in.o slp.o error_out.o
+potl.o: potl.f90 io_units.o potl_coefs.o precisn.o rd_Hfile.o rmx1_in.o
+potl_coefs.o: potl_coefs.f90 io_units.o precisn.o scaling.o slp.o $(XDR).o
+rd_Hfile.o: rd_Hfile.f90 io_units.o precisn.o rmx1_in.o scaling.o $(XDR).o error_out.o mpi_params.o
+rmx1_in.o: rmx1_in.f90 io_units.o precisn.o error_out.o mpi_params.o
+scaling.o: scaling.f90 io_units.o precisn.o
+sec_ham.o: sec_ham.f90 precisn.o io_units.o potl.o error_out.o
+$(XDR).o: $(XDR).f90 io_units.o precisn.o
+slp.o: slp.f90
+io_units.o: io_units.f90
+precisn.o: precisn.f90
+error_out.o: error_out.f90 io_units.o
+

pfarm/cpu/src_mpi_omp_mol/a_descriptor.f90

+module array_descriptor
+! define parameters used to descibe BLACS distributed arrays:
+! Time-stamp: "03/07/04 14:31:08 cjn"
+  implicit none
+
+  public
+
+  integer, parameter :: dlen_ = 9
+  integer, parameter :: block_cyclic_2d = 1
+  integer, parameter :: ctxt_ = 2
+  integer, parameter :: m_ = 3
+  integer, parameter :: n_ = 4
+  integer, parameter :: mb_ = 5
+  integer, parameter :: nb_ = 6
+  integer, parameter :: rsrc_ = 7
+  integer, parameter :: csrc_ = 8
+  integer, parameter :: lld_ = 9
+end module array_descriptor

pfarm/cpu/src_mpi_omp_mol/amp.f90

+module rw_amplitudes
+! calculated sector reduced-width amplitudes
+! Time-stamp: "2003-10-06 14:23:21 cjn"
+  use io_units, only: fo, fa
+  use precisn, only: wp
+  use def_sectors, only: nh, nsect
+  use mpi_params, only: taskid, numtasks, io_processor
+  implicit none
+
+  integer, save     :: ixdro   ! XDR unit number for rwa file
+
+  private
+  public start_sa_file, start_sa_file_id, close_sa_file, ampltd
+
+contains
+
+  subroutine start_sa_file (sqno, st, sect)
+! open XDR file to hold sector reduced width amplitudes
+    use xdr_files, only: open_xdr, xdr_io
+    use energy_grid, only: fine_mesh
+    use slp, only: qno, val_qno
+    use rmx1_in, only: xdr_amp_out, split_prop
+
+    type(qno), intent(in)   :: sqno        ! scattering q#s
+    integer, intent(in)     :: st         ! second target spin
+    integer, intent(in)     :: sect         ! sector ID
+    character(len=11)       :: hname
+    character(len=2)        :: sp, l, isect
+    character(len=1)        :: p
+    character(len=4)        :: stem
+    integer                 :: ios, qn(3), status
+    integer                 :: nspn ! scattering spin.  
+
+! File naming convention:
+!  AMPCSSLLP: coarse mesh, 1st partition (or all channels)
+!  AMPFSSLLP: fine mesh, 1st partition (or all channels)
+!  AMPCTTLLP: coarse mesh, 2nd partition. TT = target spin multiplicity
+!  AMPFTTLLP: fine mesh, 2nd partition, TT = target spin multiplicity
+! For no-exchange LS files target spin multiplicity is used
+! For J files L is replaced by 2*J and total spin multiplicity by 0 
+!                                  or target spin multiplicity by 2*K
+! For no-exchange files L is replaced by 2*K 
+!                   and total spin multiplicity by 0
+
+    ixdro = fa
+    stem = 'AMPC'
+    if (fine_mesh) stem = 'AMPF'
+    qn = val_qno(sqno)
+    if (split_prop) then   ! channel splitting in use
+       write (sp,'(i2.2)') st  ! use target spin (splitting parameter)
+    else
+       if (qn(1) < 0) qn(2) = 0 ! for jpi_coupling
+       write (sp,'(i2.2)') abs(qn(2)) ! scattering spin (0 for jpi_coupling)
+    end if
+    write (l,'(i2.2)') abs(qn(1))
+    write (p, '(i1.1)') qn(3)
+    write (isect, '(i2.2)') sect
+    hname = stem // sp // l // p // isect
+
+    if (xdr_amp_out) then   ! open XDR output file for Hamiltonian
+       ixdro = open_xdr (file=TRIM(hname), action='write')
+       call xdr_io (ixdro, nh)
+       call xdr_io (ixdro, nsect)
+    else
+       open (unit=ixdro, file=TRIM(hname), access='sequential',  &
+            status='replace', form='unformatted', action='write', &
+            iostat=ios)
+       if (ios /= 0) then
+          write (fo,'(a,i6)') 'start_sa_file: error opening ' // &
+               TRIM(hname) // ' iostat = ', ios
+          stop
+       end if
+       write (ixdro) nh, nsect
+    end if
+    if (xdr_amp_out) then
+       write (fo,'(a)') 'XDR Representation Amplitude file: ' // &
+            &TRIM(hname)
+    else
+       write (fo,'(a)') 'Native Representation Amplitude file: ' // &
+            &TRIM(hname)
+    end if
+  end subroutine start_sa_file
+
+  subroutine start_sa_file_id (sqno, st, sect)
+! open XDR file to hold sector reduced width amplitudes
+    use xdr_files, only: open_xdr, xdr_io
+    use energy_grid, only: fine_mesh
+    use slp, only: qno, val_qno
+    use rmx1_in, only: xdr_amp_out, split_prop
+
+    type(qno), intent(in)   :: sqno        ! scattering q#s
+    integer, intent(in)     :: st         ! second target spin
+    integer, intent(in)     :: sect        ! sector ID
+    character(len=13)       :: hname
+    character(len=2)        :: sp, l, isect, s_mpi
+    character(len=1)        :: p
+    character(len=4)        :: stem
+    integer                 :: ios, qn(3), status
+    integer                 :: nspn ! scattering spin.  
+
+! File naming convention:
+!  AMPCSSLLP: coarse mesh, 1st partition (or all channels)
+!  AMPFSSLLP: fine mesh, 1st partition (or all channels)
+!  AMPCTTLLP: coarse mesh, 2nd partition. TT = target spin multiplicity
+!  AMPFTTLLP: fine mesh, 2nd partition, TT = target spin multiplicity
+! For no-exchange LS files target spin multiplicity is used
+! For J files L is replaced by 2*J and total spin multiplicity by 0 
+!                                  or target spin multiplicity by 2*K
+! For no-exchange files L is replaced by 2*K 
+!                   and total spin multiplicity by 0
+
+    ixdro = fa
+    stem = 'AMPC'
+    if (fine_mesh) stem = 'AMPF'
+    qn = val_qno(sqno)
+    if (split_prop) then   ! channel splitting in use
+       write (sp,'(i2.2)') st  ! use target spin (splitting parameter)
+    else
+       if (qn(1) < 0) qn(2) = 0 ! for jpi_coupling
+       write (sp,'(i2.2)') abs(qn(2)) ! scattering spin (0 for jpi_coupling)
+    end if
+    write (l,'(i2.2)') abs(qn(1))
+    write (p, '(i1.1)') qn(3)
+    write (isect, '(i2.2)') sect
+    write (s_mpi, '(i2.2)') taskid
+    hname = stem // sp // l // p // isect // s_mpi
+
+    if (xdr_amp_out) then   ! open XDR output file for Hamiltonian
+       ixdro = open_xdr (file=TRIM(hname), action='write')
+       call xdr_io (ixdro, nh)
+       call xdr_io (ixdro, nsect)
+    else
+       open (unit=ixdro, file=TRIM(hname), access='sequential',  &
+            status='replace', form='unformatted', action='write', &
+            iostat=ios)
+       if (ios /= 0) then
+          write (fo,'(a,i6)') 'start_sa_file: error opening ' // &
+               TRIM(hname) // ' iostat = ', ios
+          stop
+       end if
+       write (ixdro) nh, nsect
+    end if
+    if (xdr_amp_out) then
+       write (fo,'(a)') 'XDR Representation Amplitude file: ' // &
+            &TRIM(hname)
+    else
+       write (fo,'(a)') 'Native Representation Amplitude file: ' // &
+            &TRIM(hname)
+    end if
+  end subroutine start_sa_file_id
+
+  subroutine close_sa_file ()
+   use xdr_files, only: close_xdr
+   call flush(ixdro)
+   call close_xdr(ixdro)
+   
+  end subroutine close_sa_file
+
+
+  subroutine ampltd (ra1, ra2, nc, nl, nh)
+! ampltd : computes surface amplitudes and writes to file
+    use xdr_files, only: xdr_io
+    use sec_ham, only: vbl, vbr
+    use dist_mat, only: evals
+    use rmx1_in, only: bug8, xdr_amp_out
+!    use blacs, only: p0, q0, ctxt, io_processor, p_error, dlen_, &
+!         ctxt_, mb_
+    use error_out, only: error_check
+    use dist_mat, only: z
+    integer, intent(in)      :: nc         ! # channels
+    integer, intent(in)      :: nl         ! # Legendre functions
+    integer, intent(in)      :: nh         ! sec Hamiltonian dimension
+    real(wp), intent(in)     :: ra1, ra2   ! end pt radii of subrange
+    real(wp), allocatable    :: ampal(:), ampar(:) ! local amps
+    real(wp), allocatable    :: evec(:)
+    integer                  :: k, j, ik, i, jj, status
+    integer                  :: mm, nn, lwork, mb, nb, rsrc, csrc
+    integer                  :: ldd, isize, istart, iend
+    integer                  :: jsize, jstart, jend
+    real(wp)                 :: alpha, beta, sml, smr
+!    integer                  :: descx(dlen_)
+
+!    if (io_processor) then
+       allocate (ampal(nc), ampar(nc), stat=status)
+       if (status /= 0) then
+          write (fo,'(a,i6)') 'ampltd: allocation error = ', status
+          STOP
+       end if
+    
+       if (bug8 > 0) then
+          write (fo,'(a,e20.12)') 'Sector eigenvalue sum = ', SUM(evals)
+       else if (bug8 > 1) then
+          write (fo,'(a)') 'R-Matrix Eigenvalues (Ryd)'
+          write (fo,'(5f12.5)') evals
+       end if
+
+! write eigenvalues to disk
+       if (xdr_amp_out) then
+          call xdr_io (ixdro, evals, nh)  ! write R-matrix eigenvalues
+       else
+          write (ixdro) evals
+       end if
+!    end if
+
+    allocate (evec(nh), stat=status)
+    call error_check (status, 'ampltd: evec allocation')
+    lwork = 96
+!    lwork = descz(mb_)   ! eigenvectgor block size
+!    ctxt = descz(ctxt_)  ! context
+
+! form descriptor for evec
+    mm = MAX(1, MIN(nh, lwork))
+    nn = MAX(1, INT(lwork / mm))
+    mb = mm
+    nb = nn
+!    rsrc = p0
+!    csrc = q0
+!    ldd = MAX(1, mm)
+!    call descset (descx, mm, nn, mb, nb, rsrc, csrc, ctxt, ldd)
+
+    alpha = 1.0_wp
+    beta = 0.0_wp
+    sml = 0.0_wp
+    smr = 0.0_wp
+
+    evec = z
+    jst: do jstart = 1, nh, nn
+       jend = MIN(nh, jstart+nn-1)
+       jsize = jend - jstart + 1
+
+ !      row_blocks: do istart = 1, nh, mm
+ !         iend = MIN(nh, istart+mm-1)
+ !         isize = iend - istart + 1
+ !         call pdgeadd ('notrans', isize, jsize, alpha, z, istart, &
+ !              jstart, descz, beta, evec(istart:iend), 1, 1, descx)
+ !      end do row_blocks
+ !      call BLACS_BARRIER (ctxt, 'A')
+! now have the full eigenvector on the i/o processor:
+
+!       if (io_processor) then
+          ampal = 0.0_wp
+          ampar = 0.0_wp
+          ik = 0
+          do i = 1, nc
+             do j = 1, nl
+                ik = ik + 1
+                ampal(i) = ampal(i) + evec(ik) * vbl(j)
+                ampar(i) = ampar(i) + evec(ik) * vbr(j)
+             end do
+          end do
+          
+          if (bug8 > 0) then
+             sml = sml + SUM(ampal)
+             smr = smr + SUM(ampar)
+          else if (bug8 > 2) then
+             write (fo,'(a,f10.5,a,i4)') 'Amplitudes (a.u.) at ', ra1,&
+                  ' Column: ', jstart
+             write (fo,'(5e14.6)') ampal
+             write (fo,'(a,f10.5,a,i4)') 'Amplitudes (a.u.) at ', ra2,&
+                  ' Column: ', jstart
+             write (fo,'(5e14.6)') ampar
+          end if
+
+! write sector reduced width amplitudes to disk
+          if (xdr_amp_out) then
+             call xdr_io (ixdro, ampal, nc) ! write left amplitudes to disk
+             call xdr_io (ixdro, ampar, nc) ! write right amplitudes to disk
+          else
+             write (ixdro) ampal(1:nc)
+             write (ixdro) ampar(1:nc)
+          end if
+ !      end if
+    end do jst
+!    if (io_processor) then
+       if (bug8 > 0) write (fo,'(a,2e20.12)') 'RWAmplitude sums = ', &
+            sml, smr
+       deallocate (ampal, ampar, evec, stat=status)
+!    else
+!       deallocate (evec, stat=status)
+!    end if
+    call error_check (status, 'ampltd: deallocation')
+  end subroutine ampltd
+end module rw_amplitudes

pfarm/cpu/src_mpi_omp_mol/amp.f90.ori

+module rw_amplitudes
+! calculated sector reduced-width amplitudes
+! Time-stamp: "2003-10-06 14:23:21 cjn"
+  use io_units, only: fo, fa
+  use precisn, only: wp
+  use def_sectors, only: nh, nsect
+  implicit none
+
+  integer, save     :: ixdro   ! XDR unit number for rwa file
+
+  private
+  public start_sa_file, ampltd
+
+contains
+
+  subroutine start_sa_file (sqno, st)
+! open XDR file to hold sector reduced width amplitudes
+    use xdr_files, only: open_xdr, xdr_io
+    use energy_grid, only: fine_mesh
+    use slp, only: qno, val_qno
+    use rmx1_in, only: xdr, split_prop
+    type(qno), intent(in)   :: sqno        ! scattering q#s
+    integer, intent(in)     :: st         ! second target spin
+    character(len=9)        :: hname
+    character(len=2)        :: sp, l
+    character(len=1)        :: p
+    character(len=4)        :: stem
+    integer                 :: ios, qn(3), status
+    integer                 :: nspn ! scattering spin.  
+
+! File naming convention:
+!  AMPCSSLLP: coarse mesh, 1st partition (or all channels)
+!  AMPFSSLLP: fine mesh, 1st partition (or all channels)
+!  AMPCTTLLP: coarse mesh, 2nd partition. TT = target spin multiplicity
+!  AMPFTTLLP: fine mesh, 2nd partition, TT = target spin multiplicity
+! For no-exchange LS files target spin multiplicity is used
+! For J files L is replaced by 2*J and total spin multiplicity by 0 
+!                                  or target spin multiplicity by 2*K
+! For no-exchange files L is replaced by 2*K 
+!                   and total spin multiplicity by 0
+
+    ixdro = fa
+    stem = 'AMPC'
+    if (fine_mesh) stem = 'AMPF'
+    qn = val_qno(sqno)
+    if (split_prop) then   ! channel splitting in use
+       write (sp,'(i2.2)') st  ! use target spin (splitting parameter)
+    else
+       if (qn(1) < 0) qn(2) = 0 ! for jpi_coupling
+       write (sp,'(i2.2)') abs(qn(2)) ! scattering spin (0 for jpi_coupling)
+    end if
+    write (l,'(i2.2)') abs(qn(1))
+    write (p, '(i1.1)') qn(3)
+    hname = stem // sp // l // p
+
+    if (xdr) then   ! open XDR output file for Hamiltonian
+       ixdro = open_xdr (file=TRIM(hname), action='write')
+       call xdr_io (ixdro, nh)
+       call xdr_io (ixdro, nsect)
+    else
+       open (unit=ixdro, file=TRIM(hname), access='sequential',  &
+            status='replace', form='unformatted', action='write', &
+            iostat=ios)
+       if (ios /= 0) then
+          write (fo,'(a,i6)') 'start_sa_file: error opening ' // &
+               TRIM(hname) // ' iostat = ', ios
+          stop
+       end if
+       write (ixdro) nh, nsect
+    end if
+    if (xdr) then
+       write (fo,'(a)') 'XDR Representation Amplitude file: ' // &
+            &TRIM(hname)
+    else
+       write (fo,'(a)') 'Native Representation Amplitude file: ' // &
+            &TRIM(hname)
+    end if
+  end subroutine start_sa_file
+
+  subroutine ampltd (ra1, ra2, nc, nl, nh)
+! ampltd : computes surface amplitudes and writes to file
+    use xdr_files, only: xdr_io
+    use sec_ham, only: vbl, vbr
+    use dist_mat, only: evals
+    use rmx1_in, only: bug8, xdr
+    use blacs, only: p0, q0, ctxt, io_processor, p_error, dlen_, &
+         ctxt_, mb_
+    use dist_mat, only: z, descz
+    integer, intent(in)      :: nc         ! # channels
+    integer, intent(in)      :: nl         ! # Legendre functions
+    integer, intent(in)      :: nh         ! sec Hamiltonian dimension
+    real(wp), intent(in)     :: ra1, ra2   ! end pt radii of subrange
+    real(wp), allocatable    :: ampal(:), ampar(:) ! local amps
+    real(wp), allocatable    :: evec(:)
+    integer                  :: k, j, ik, i, jj, status
+    integer                  :: mm, nn, lwork, mb, nb, rsrc, csrc
+    integer                  :: ldd, isize, istart, iend
+    integer                  :: jsize, jstart, jend
+    real(wp)                 :: alpha, beta, sml, smr
+    integer                  :: descx(dlen_)
+
+    if (io_processor) then
+       allocate (ampal(nc), ampar(nc), stat=status)
+       if (status /= 0) then
+          write (fo,'(a,i6)') 'ampltd: allocation error = ', status
+          call BLACS_ABORT (ctxt, status)
+       end if
+    
+       if (bug8 > 0) then
+          write (fo,'(a,e20.12)') 'Sector eigenvalue sum = ', SUM(evals)
+       else if (bug8 > 1) then
+          write (fo,'(a)') 'R-Matrix Eigenvalues (Ryd)'
+          write (fo,'(5f12.5)') evals
+       end if
+
+! write eigenvalues to disk
+       if (xdr) then
+          call xdr_io (ixdro, evals, nh)  ! write R-matrix eigenvalues
+       else
+          write (ixdro) evals
+       end if
+    end if
+
+    allocate (evec(nh), stat=status)
+    call p_error (status, 'ampltd: evec allocation')
+    lwork = descz(mb_)   ! eigenvectgor block size
+    ctxt = descz(ctxt_)  ! context
+
+! form descriptor for evec
+    mm = MAX(1, MIN(nh, lwork))
+    nn = MAX(1, INT(lwork / mm))
+    mb = mm
+    nb = nn
+    rsrc = p0
+    csrc = q0
+    ldd = MAX(1, mm)
+    call descset (descx, mm, nn, mb, nb, rsrc, csrc, ctxt, ldd)
+
+    alpha = 1.0_wp
+    beta = 0.0_wp
+    sml = 0.0_wp
+    smr = 0.0_wp
+    jst: do jstart = 1, nh, nn
+       jend = MIN(nh, jstart+nn-1)
+       jsize = jend - jstart + 1
+       evec = 0.0_wp
+
+       row_blocks: do istart = 1, nh, mm
+          iend = MIN(nh, istart+mm-1)
+          isize = iend - istart + 1
+          call pdgeadd ('notrans', isize, jsize, alpha, z, istart, &
+               jstart, descz, beta, evec(istart:iend), 1, 1, descx)
+       end do row_blocks
+       call BLACS_BARRIER (ctxt, 'A')
+! now have the full eigenvector on the i/o processor:
+
+       if (io_processor) then
+          ampal = 0.0_wp
+          ampar = 0.0_wp
+          ik = 0
+          do i = 1, nc
+             do j = 1, nl
+                ik = ik + 1
+                ampal(i) = ampal(i) + evec(ik) * vbl(j)
+                ampar(i) = ampar(i) + evec(ik) * vbr(j)
+             end do
+          end do
+          
+          if (bug8 > 0) then
+             sml = sml + SUM(ampal)
+             smr = smr + SUM(ampar)
+          else if (bug8 > 2) then
+             write (fo,'(a,f10.5,a,i4)') 'Amplitudes (a.u.) at ', ra1,&
+                  ' Column: ', jstart
+             write (fo,'(5e14.6)') ampal
+             write (fo,'(a,f10.5,a,i4)') 'Amplitudes (a.u.) at ', ra2,&
+                  ' Column: ', jstart
+             write (fo,'(5e14.6)') ampar
+          end if
+
+! write sector reduced width amplitudes to disk
+          if (xdr) then
+             call xdr_io (ixdro, ampal, nc) ! write left amplitudes to disk
+             call xdr_io (ixdro, ampar, nc) ! write right amplitudes to disk
+          else
+             write (ixdro) ampal(1:nc)
+             write (ixdro) ampar(1:nc)
+          end if
+       end if
+    end do jst
+    if (io_processor) then
+       if (bug8 > 0) write (fo,'(a,2e20.12)') 'RWAmplitude sums = ', &
+            sml, smr
+       deallocate (ampal, ampar, evec, stat=status)
+    else
+       deallocate (evec, stat=status)
+    end if
+    call p_error (status, 'ampltd: deallocation')
+  end subroutine ampltd
+end module rw_amplitudes

pfarm/cpu/src_mpi_omp_mol/blacs.f90

+module blacs
+! store blacs context and processor data
+! Time-stamp: "2003-09-23 13:41:59 cjn"
+  implicit none
+
+! define parameters used to descibe BLACS distributed arrays:
+  integer, parameter :: dlen_ = 9
+  integer, parameter :: block_cyclic_2d = 1
+  integer, parameter :: ctxt_ = 2
+  integer, parameter :: m_ = 3
+  integer, parameter :: n_ = 4
+  integer, parameter :: mb_ = 5
+  integer, parameter :: nb_ = 6
+  integer, parameter :: rsrc_ = 7
+  integer, parameter :: csrc_ = 8
+  integer, parameter :: lld_ = 9
+
+! initial linear context:
+  integer, save      :: ictxt  ! initial blacs context handle
+  integer, save      :: iam    ! processor id in initial context
+  integer, save      :: nprocs ! total # processors
+
+! context for main calculation: p * q <= nprocs
+  integer, save      :: p, q   ! # processors in grid
+  integer, save      :: ctxt   ! blacs context handle
+  integer, save      :: mycol, myrow ! processor coordinates
+  logical, save      :: io_processor   ! i/o processor flag
+  integer, save      :: p0, q0 ! i/o processor coordinates
+  integer, save      :: nblock ! blacs blocking factor
+  integer, save      :: mynumrows   ! Hamiltonian local # rows
+  integer, save      :: mynumcols   ! Hamiltonian local # cols
+
+  integer, save      :: cur_ctxt=-5 ! current contxt for p_error
+
+  private
+  public p_error, set_cur_ctxt
+  public dlen_, block_cyclic_2d, ctxt_, m_, n_, mb_, nb_
+  public rsrc_, csrc_, lld_
+  public ictxt, nprocs, iam
+  public ctxt, io_processor, p, q, mycol, myrow, p0, q0
+  public nblock, mynumrows, mynumcols
+
+contains
+
+  subroutine set_cur_ctxt (n)
+! Define the Blacs context handle to be used by p_error
+    integer, intent(in)    :: n
+    select case (n)
+    case (1)
+       cur_ctxt = ictxt
+    case (2)
+       cur_ctxt = ctxt
+    case default
+       cur_ctxt = ictxt
+    end select
+  end subroutine set_cur_ctxt
+
+  subroutine p_error (err, msg)
+    use io_units, only: fo
+    integer, intent(in)          :: err ! flag
+    character(len=*), intent(in) :: msg   ! error message
+    integer                      :: no
+
+    no = 1
+    if (cur_ctxt /= -5) then
+       call igsum2d (cur_ctxt, 'all', ' ', 1, 1, err, 1, -1, -1)
+       if (err /= 0) then
+          if (iam == 0) write (fo, fmt='(a)') msg
+          call BLACS_ABORT (cur_ctxt, no)
+       end if
+    else
+       if (err /= 0) then
+          write (fo, fmt='(a)') msg
+          call BLACS_GET (-1, 0, cur_ctxt) ! find default context
+          call BLACS_ABORT (cur_ctxt, no)
+       end if
+    end if
+  end subroutine p_error
+
+end module blacs

pfarm/cpu/src_mpi_omp_mol/def_sectors.f90

+module def_sectors
+! Time-stamp: "2003-10-14 12:51:24 cjn"
+  use precisn, only: wp
+  use io_units, only: fo
+  use error_out, only: error_check
+  implicit none
+
+  private
+  public get_nsect_nl, hack_get_nsect_nl, reset_sec, sectors
+  public nsect, nh, nl, nx, asect, w_sect
+
+  integer, save               :: nsect   ! # sectors
+  integer, save               :: nh      ! Hamiltonian dimension (all channels)
+  integer, save               :: nl      ! Legendre order
+  integer, save               :: nx      ! # Gauss points
+  real(wp), allocatable, save :: asect(:)
+  real(wp)                    :: w_sect
+
+contains
+
+  subroutine reset_sec
+! deallocate memory
+    integer            :: status
+    deallocate (asect, stat=status)
+    call error_check(status, " Error in reset_sec ")
+  end subroutine reset_sec
+
+  subroutine sectors (nc)
+ ! calculate sector end; rlast lies on the end of last sector
+    use energy_grid, only: fine_mesh, ebig
+    use rmx1_in, only: nhmax, ncmax, rafin
+    use hfile_data, only: rmatr
+!    use blacs, only: ctxt, io_processor
+    integer, intent(in) :: nc    ! # channels
+    real(wp)            :: rfirst, rlast ! propagation range
+    integer             :: ibuf(4), status, i
+    real(wp)            :: dum(1)
+
+
+       rfirst = rmatr
+       rlast = rafin
+       if (rafin < rmatr) then
+          write (fo,'(a)') 'sectors: propagation range error'
+          write (fo,'(a,2f12.6)') 'rmatr, rafin = ', rmatr, rafin
+          STOP
+       end if
+       call get_nsect_nl (nhmax, ncmax, rfirst, rlast, ebig, fine_mesh)
+       ! AGS - The below call can be used to hack Ham matrix size for benchmarking purposes
+       ! call hack_get_nsect_nl (nhmax, ncmax, rfirst, rlast, ebig, fine_mesh)
+       w_sect = (rlast - rfirst) / REAL(nsect,wp)   ! actual sector width
+       allocate (asect(nsect+1), stat=status)
+       asect(1) = rfirst
+       do i = 1, nsect
+          asect(i+1) = asect(1) + i * w_sect
+       end do
+       nh = nl * nc
+
+! output sector radius
+       write (fo,'(/a,i4)') 'Number of basis functions *** HACKED for Fine & Coarse Regions ***, nl = ', nl
+       write (fo,'(a,i4)') 'Number of abscissae, nx       = ', nx
+       write (fo,'(/a,i8/)') 'Hamiltonian dimension = ', nh
+       write (fo,'(a,f12.6)') 'Sector width   = ', w_sect
+       do i = 1, nsect+1
+          write (fo,'(i4,5x,a,f12.6)') i, 'asect = ', asect(i)
+       end do
+
+   end subroutine sectors
+
+  subroutine get_nsect_nl (nhmax, ncmax, rfirst, rlast, ebig, &
+       fine_mesh)
+! use PGB formula to determine sector sizes and # of basis functions
+! formula:   kmax =  nl * pi   ~=   nl
+!                   ---------       --       delta(a) = rlast-rfirst
+!                   3*delta(a)    delta(a)
+    integer, intent(in)  :: nhmax   ! max dimension of sector Hamiltonian
+    integer, intent(in)  :: ncmax   ! max # channels for this L
+    real(wp), intent(in) :: rfirst ! initial radius
+    real(wp), intent(in) :: rlast  ! final radius
+    real(wp), intent(in) :: ebig   ! maximum energy (for this mesh)
+    logical, intent(in)  :: fine_mesh ! fine mesh flag
+    real(wp)             :: ansect, kmax , anlmax
+
+    kmax = SQRT(ebig)
+    if (fine_mesh) then
+       anlmax = nhmax / ncmax
+       ansect = (rlast - rfirst) * kmax / anlmax
+       nsect = MAX(INT(ansect+1.0_wp), 2)
+       nl = MAX(INT(kmax*(rlast-rfirst)/REAL(nsect,wp)+1.0_wp), 10)
+    else ! coarse mesh, therefore set nl to 10
+       nl = 10
+       ansect = (rlast - rfirst) * kmax / REAL(nl,wp)
+       nsect = MAX(INT(ansect+1.0_wp), 2)
+    end if
+    nx = nl + 5
+  end subroutine get_nsect_nl
+
+  subroutine hack_get_nsect_nl (nhmax, ncmax, rfirst, rlast, ebig, &
+       fine_mesh)
+! **** AGS Hacked to provide larger Hamiltonian matrices ****
+! use PGB formula to determine sector sizes and # of basis functions
+! formula:   kmax =  nl * pi   ~=   nl
+!                   ---------       --       delta(a) = rlast-rfirst
+!                   3*delta(a)    delta(a)
+    integer, intent(in)  :: nhmax   ! max dimension of sector Hamiltonian
+    integer, intent(in)  :: ncmax   ! max # channels for this L
+    real(wp), intent(in) :: rfirst ! initial radius
+    real(wp), intent(in) :: rlast  ! final radius
+    real(wp), intent(in) :: ebig   ! maximum energy (for this mesh)
+    logical, intent(in)  :: fine_mesh ! fine mesh flag
+    real(wp)             :: ansect, kmax , anlmax
+
+    kmax = SQRT(ebig)
+    if (fine_mesh) then
+       anlmax = nhmax / ncmax
+       ansect = (rlast - rfirst) * kmax / anlmax
+       nsect = MAX(INT(ansect+1.0_wp), 2)
+       nl = MAX(INT(kmax*(rlast-rfirst)/REAL(nsect,wp)+1.0_wp), 10)
+!       nl = MAX(INT(kmax*(rlast-rfirst)/REAL(nsect,wp)+1.0_wp), 4)
+       nl = 6
+    else ! coarse mesh, therefore set nl to 10
+       nl = 4
+!       nl = 10
+       ansect = (rlast - rfirst) * kmax / REAL(nl,wp)
+       nsect = MAX(INT(ansect+1.0_wp), 2)
+    end if
+    nx = nl + 5
+  end subroutine hack_get_nsect_nl
+
+end module def_sectors

pfarm/cpu/src_mpi_omp_mol/def_sectors_bench.f90

+module def_sectors
+! Time-stamp: "2003-10-14 12:51:24 cjn"
+  use precisn, only: wp
+  use io_units, only: fo
+  use error_out, only: error_check
+  use mpi_params, only: io_processor
+  implicit none
+
+  private
+  public reset_sec, sectors
+  public nsect, nh, nl, nx, asect, w_sect
+
+  integer, save               :: nsect   ! # sectors
+  integer, save               :: nh      ! Hamiltonian dimension (all channels)
+  integer, save               :: nl      ! Legendre order
+  integer, save               :: nx      ! # Gauss points
+  real(wp), allocatable, save :: asect(:)
+  real(wp)                    :: w_sect
+
+contains
+
+  subroutine reset_sec
+! deallocate memory
+    integer            :: status
+    deallocate (asect, stat=status)
+    call error_check(status, " Error in reset_sec ")
+  end subroutine reset_sec
+
+  subroutine sectors (nc)
+ ! calculate sector end; rlast lies on the end of last sector
+    use energy_grid, only: fine_mesh, ebig
+    use rmx1_in, only: nhmax, ncmax, rafin
+    use hfile_data, only: rmatr
+    integer, intent(in) :: nc    ! # channels
+    real(wp)            :: rfirst, rlast ! propagation range
+    integer             :: ibuf(4), status, i
+    real(wp)            :: dum(1)
+
+
+       rfirst = rmatr
+       rlast = rafin
+       if (rafin < rmatr) then
+          write (fo,'(a)') 'sectors: propagation range error'
+          write (fo,'(a,2f12.6)') 'rmatr, rafin = ', rmatr, rafin
+          STOP
+       end if
+
+       ! Benchmark version of code allowing easier problem size variation by user
+       call bench_get_nsect_nl (nhmax, ncmax, rfirst, rlast, ebig, fine_mesh)
+
+       w_sect = (rlast - rfirst) / REAL(nsect,wp)   ! actual sector width
+       allocate (asect(nsect+1), stat=status)
+       asect(1) = rfirst
+       do i = 1, nsect
+          asect(i+1) = asect(1) + i * w_sect
+       end do
+       nh = nl * nc
+
+! output sector radius
+       if(io_processor) then
+         write (fo,'(/a,i4)') 'Number of basis functions (Benchmark Case) = ', nl
+         write (fo,'(a,i4)') 'Number of abscissae, nx = ', nx
+         write (fo,'(a,i8)') 'Hamiltonian dimension = ', nh
+         write (fo,'(a,i4)') 'Total Number of Sectors = ', nsect
+         write (fo,'(a,f12.6)') 'Sector width   = ', w_sect
+         do i = 1, nsect+1
+           write (fo,'(i4,5x,a,f12.6)') i, 'asect = ', asect(i)
+         end do
+         write (fo,'(/)')
+         call flush(fo)
+       end if
+   end subroutine sectors
+
+  subroutine get_nsect_nl (nhmax, ncmax, rfirst, rlast, ebig, &
+       fine_mesh)
+! use PGB formula to determine sector sizes and # of basis functions
+! formula:   kmax =  nl * pi   ~=   nl
+!                   ---------       --       delta(a) = rlast-rfirst
+!                   3*delta(a)    delta(a)
+    integer, intent(in)  :: nhmax   ! max dimension of sector Hamiltonian
+    integer, intent(in)  :: ncmax   ! max # channels for this L
+    real(wp), intent(in) :: rfirst ! initial radius
+    real(wp), intent(in) :: rlast  ! final radius
+    real(wp), intent(in) :: ebig   ! maximum energy (for this mesh)
+    logical, intent(in)  :: fine_mesh ! fine mesh flag
+    real(wp)             :: ansect, kmax , anlmax
+
+    kmax = SQRT(ebig)
+    if (fine_mesh) then
+       anlmax = nhmax / ncmax
+       ansect = (rlast - rfirst) * kmax / anlmax
+       nsect = MAX(INT(ansect+1.0_wp), 2)
+       nl = MAX(INT(kmax*(rlast-rfirst)/REAL(nsect,wp)+1.0_wp), 10)
+    else ! coarse mesh, therefore set nl to 10
+       nl = 10
+       ansect = (rlast - rfirst) * kmax / REAL(nl,wp)
+       nsect = MAX(INT(ansect+1.0_wp), 2)
+    end if
+    nx = nl + 5
+  end subroutine get_nsect_nl
+
+  subroutine bench_get_nsect_nl (nhmax, ncmax, rfirst, rlast, ebig, &
+       fine_mesh)
+    ! **** Benchmark version of routine to allow greater user control of problem size ****
+    integer, intent(in)  :: nhmax   ! max dimension of sector Hamiltonian
+    integer, intent(in)  :: ncmax   ! max # channels for this L
+    real(wp), intent(in) :: rfirst ! initial radius
+    real(wp), intent(in) :: rlast  ! final radius
+    real(wp), intent(in) :: ebig   ! maximum energy (for this mesh)
+    logical, intent(in)  :: fine_mesh ! fine mesh flag
+    real(wp)             :: ansect, kmax , anlmax
+    CHARACTER(len=255) :: c_nl, c_nsect
+    integer :: i_nl, i_nsect, status_val
+
+    kmax = SQRT(ebig)
+    if (fine_mesh) then
+
+       nsect = -1
+       i_nl = -1
+       CALL get_environment_variable("RMX_NL_FINE", c_nl, STATUS=status_val)
+       read( c_nl, '(i10)' ) i_nl
+       if (i_nl > 0 ) then
+          nl = i_nl
+       else
+          nl = 12
+       end if
+
+       CALL get_environment_variable("RMX_NSECT_FINE", c_nsect, STATUS=status_val)
+       read( c_nsect, '(i10)' ) i_nsect
+       if (i_nsect > 0 ) then
+          nsect = i_nsect
+       else
+          nsect = 5
+       end if
+
+       if(io_processor) then
+         write (*,'(a,i0)') "BENCHMARK CODE: Number of Sectors in Fine Region = ", nsect
+         write (*,'(a,i0)') "BENCHMARK CODE: Number of Basis Functions used in Fine Region = ", nl
+       end if
+
+    else ! coarse mesh
+       i_nsect=-1
+       i_nl=-1
+       CALL get_environment_variable("RMX_NL_COARSE", c_nl, STATUS=status_val)
+       read( c_nl, '(i10)' ) i_nl
+       if (i_nl > 0 ) then
+          nl = i_nl
+       else
+          nl = 6
+       end if
+
+       CALL get_environment_variable("RMX_NSECT_COARSE", c_nsect, STATUS=status_val)
+       read( c_nsect, '(i10)' ) i_nsect
+       if (i_nsect > 0 ) then
+          nsect = i_nsect
+       else
+          nsect = 20
+       end if
+
+       if(io_processor) then
+         write (*,'(a,i0)') "BENCHMARK CODE: Number of Sectors in Coarse Region = ", nsect
+         write (*,'(a,i0)') "BENCHMARK CODE: Number of Basis Functions used in Coarse Region = ", nl
+       end if
+
+    end if
+
+    ! number of Gauss points
+    nx = nl + 5
+
+  end subroutine bench_get_nsect_nl
+
+end module def_sectors

pfarm/cpu/src_mpi_omp_mol/dist_mat.f90

+module dist_mat
+! Form diestributed sector Hamiltonian matrix and diagonalize
+! Time-stamp: "2003-09-23 17:08:23 cjn"
+  use io_units, only: fo
+  use precisn, only: wp, ip_long
+!  use blacs, only: ctxt, myrow, mycol, p, q, nblock, mynumrows,&
+!       mynumcols, p_error, dlen_, nb_, csrc_
+  use error_out, only: error_check
+  implicit none
+
+  integer, save               :: ny           ! a dimension
+  real(wp), allocatable, save :: a(:)         ! distributed matrix
+  real(wp), allocatable, save :: z(:)         ! eigenvectors
+  real(wp), allocatable, save :: evals(:)     ! eigenvalues
+!  integer, save               :: desca(dlen_) ! Hamiltonian descriptor
+!  integer, save               :: descz(dlen_) ! Eigenvector descriptor
+
+  private
+  public A_fill, sh_diag, kill_az
+  public evals, z
+
+contains
+
+  subroutine A_fill (n)
+! fill lower triangle of symmetric matrix, A
+    use sec_ham, only: h_el
+    integer, intent(in)   :: n       ! matrix dimension
+    integer               :: j, jj, jloc, jglob, i, ii, iloc, iglob
+    integer               :: status, info, numroc
+    integer :: ind
+
+    ny = n
+
+! Compute local dimensions of matrix A
+!    mynumrows = NUMROC(n, nblock, myrow, 0, p)
+!    mynumcols = NUMROC(n, nblock, mycol, 0, q)
+
+! Allocate local part of A
+    allocate (a(ny*ny), z(ny*ny), stat=status)
+    call error_check (status, ' Allocation Error, A_fill')
+
+    a = 0.0_wp        ! initialize A
+
+! Form array descriptor od destributed matrix A
+!    call descinit (desca, n, n, nblock, nblock, 0, 0, ctxt, mynumrows,&
+!         info)
+!    descz = desca
+
+!    do j = 1, mynumcols, nblock                   ! local column blocks
+!       do jj = 1, MIN(nblock, mynumcols-j+1)      ! all columns
+!          jloc = j - 1 + jj                       ! local column index
+!          jglob = (j-1) * q + mycol * nblock + jj ! global column index
+
+!          do i = 1, mynumrows, nblock
+!             do ii = 1, MIN(nblock, mynumrows-i+1)
+!                iloc = i - 1 + ii                       ! local index
+!                iglob = (i-1) * p + myrow * nblock + ii ! global index
+!                if (iglob > jglob) cycle
+                
+!                a(iloc+(jloc-1)*mynumrows) = h_el(iglob, jglob)
+
+!             end do
+!          end do
+!       end do
+!    end do
+
+! sequential version - full matrix
+  ind = 0
+  do j = 1, ny
+    do i = 1, j
+      ind = ind + 1
+      a(ind) = h_el(i,j)
+    end do
+  end do
+
+  end subroutine A_fill
+
+  subroutine sh_diag(imesh, sect)
+    ! diagonalize the matrix a
+    integer, intent(in)          :: imesh, sect
+    real(wp), allocatable        :: work(:)
+    integer, allocatable         :: iwork(:)
+    integer                      :: lwork, liwork
+    real(wp)                     :: alpha, beta
+    character(len=1)             :: jobz, uplo
+    logical                      :: ismycol
+    real(wp)                     :: dnorm2
+    integer                      :: status, ii, nstrt, trilwmin
+    integer                      :: np, nq, lda, iproc, info, numroc
+    integer                      :: indxg2p, lwork1, ip, iq, imyrow, imycol
+    integer :: rows, cols
+    real(wp)                     :: t0, t1
+    integer(ip_long)             :: c0, c1, cr
+    real(wp), external :: dnrm2
+
+! initialize eigenvector distributed array, z.
+     alpha = 0.0_wp
+     beta = 1.0_wp
+!     call pdlaset ('all', ny, ny, alpha, beta, z, 1, 1, descz)
+!     np = mynumrows
+!     nq = mynumcols5
+
+! perform diagonalization:
+     jobz = 'V'
+     uplo = 'U'
+     liwork = 3 + 5*ny
+     lwork = 1 + 6*ny + 2*ny*ny
+
+     allocate (work(lwork), iwork(liwork), evals(ny), stat=status)
+     call error_check (status, 'sh_diag: allocation error')
+
+     write (fo,'(/,a,i0,a,i0,a)') ' Mesh ',imesh,', Sector ',sect,': Starting DSYEVD '
+     call flush(fo)
+
+     call cpu_time (t0)
+     call system_clock (count=c0)
+
+     call dsyevd (jobz, uplo, ny, a, ny, evals, work, lwork, iwork, liwork, info)
+
+     write (fo,'(/,a,i0,a,i0,a)') ' Mesh ',imesh,', Sector ',sect,': Ended DSYEVD '
+
+     call cpu_time (t1)
+     call system_clock (count=c1, count_rate=cr)
+     
+     write (fo,'(/,a,i0,a,i0,a,f14.3,a)') ' Mesh ',imesh,', Sector ',sect,&
+          ': DSYEVD CPU time = ', t1 - t0, ' secs'
+     write (fo,'(a,i0,a,i0,a,f14.3,a,/)') ' Mesh ',imesh,', Sector ',sect,&
+          ': DSYEVD Elapsed time = ', REAL(c1-c0,wp) / REAL(cr,wp), ' secs'
+
+     call error_check (info, 'sh_diag: pdsyevd error')
+
+     deallocate (work, iwork, stat=status)
+     call error_check (status, 'sh_diag: deallocation error')
+
+     write (fo,'(/,a,i0,a,i0,a,5f14.4)') ' Mesh ',imesh,', Sector ',sect,': first five eigenvalues = ', evals(1:5)
+     write (fo,'(a,i0,a,i0,a,5f14.4,/)') ' Mesh ',imesh,', Sector ',sect,': final five eigenvalues = ', evals(ny-4:ny)
+
+     call flush(fo)
+
+     z = a
+     
+     evecs: do ii = 1, ny ! normalize eigenvectors
+! indxg2p computes process coord which posseses entry of a
+! distributed matrix specified by a global index INDXGLOB.
+!        iproc = indxg2p(ii, descz(nb_), mycol, descz(csrc_), q)
+!        ismycol = (iproc == mycol)
+!        if (ismycol) then
+!           call pdnrm2 (ny, dnorm2, z, 1, ii, descz, 1)
+            
+           dnorm2 = dnrm2(ny,z,1)    
+           alpha = 1.0_wp / REAL(dnorm2,wp)
+           call dscal(ny, alpha, z, 1)
+!        end if
+     end do evecs
+   end subroutine sh_diag
+
+   subroutine kill_az
+     integer       :: status
+     deallocate (a, z, evals, stat=status)
+     call error_check (status, 'kill_az: deallocation error')
+   end subroutine kill_az
+
+end module dist_mat