% Technical report on Global Arrays and High-Performance Fortran
% prepared for Battelle/Pacific Northwest National Laboratory under
% contract 315647-A9E.
%
% $Id: ga-hpf.tex,v 1.1 1997/06/18 19:59:24 bernhold Exp $
%

\documentclass[12pt]{article}

\begin{document}

\title{High-Performance Fortran as a Possible
Successor to Global Arrays in the NWChem Parallel Computational
Chemistry Code\thanks{Prepared for Battelle/Pacific Northwest National
Laboratory under contract 315647-A9E.}}

\author{Yuhong Wen, D. Bryan Carpenter, \\
Erol Akarsu, Tomasz Haupt, and David E. Bernholdt\thanks{Author for
correspondence. E-mail address \emph{bernhold@npac.syr.edu}.}\\
\\
Northeast Parallel Architectures Center\\
Syracuse University\\
111 College Place\\
Syracuse, NY 13422--4100}

\maketitle

\begin{abstract}
This is the abstract
\end{abstract}

\section{Introduction}

When the development of the NWChem parallel computational chemistry
package began at the Pacific Northwest National Laboratory (PNNL) in
1993, the question of what parallel programming environment to use
received careful attention.  It was clear that a straightforward
message-passing model would place a significant burden on the
programmer.  High-Performance Fortran (HPF) provided a higher-level
programming model, but suffered from several other problems.  In the
first place, version 1.0 of the HPF standard addressed only
data-parallel operations, while it was clear that both data- and
task-parallel capabilities would be required for the algorithms
envisioned in NWChem.  Second, at that time, the standard had not yet
been formally approved, and such HPF implementations as did exist were
based on drafts of the standard, were incomplete in their
implementation of the language, and were rather immature.

As a result, the NWChem developers decided to implement their own
parallel programming model, \emph{Global Arrays}, which would be
modeled to a large extent on HPF.  The hope was that a future version
of HPF would support the capabilities required by NWChem, and that
given time, compilers would become more complete and robust, making it
possible to transition NWChem from the global array toolkit (GA) to
HPF, thereby acruing the benefits of an accepted standard and
vendor-supported compilers.

Since the HPF 1.0 standard has now been available for several years,
and since version 2.0 of the language standard (currently in a public
comment period prior to formal acceptance) includes some task
parallelism and other features required by NWChem algorithms, it is
appropriate at this time to evaluate the present state of HPF and the
prospects for using it to replace the GA toolkit in NWChem and other
similar computational chemistry applications.

There are two components to this analysis.  The first is an
examination of current HPF implementations, and their performance in
comparison with GAs.  The second is an examination of whether or not
the HPF 2.0 standard supports the task parallel and other
capabilities required by the algorithms inplemented in NWChem.

\section{Performance of HPF and GAs}

\section{HPF 2.0 Capabilities Analysis}

Version 2.0 of the High-Performance Fortran standard is currently in a
public comment period prior to formal acceptance.  This version of HPF
incorporates some support for task-based parallelism as well as other
features similar to those in the GA toolkit.  In order to understand the
degree to which HPF 2.0 supports the GA programming model, we have
looked at the requirements of two algorithms implemented in NWChem
which are exemplars of the kind of task parallelism used throughout
the code.

We must note that since there are no HPF 2.0 compilers presently
available, the proposed implementations we outline are ``thought
experiments'' and subject to all of the uncertainties of both
interpretation of the standard itself, and actual implementations.

\subsection{NWChem Distributed Data SCF Algorithm}

We investigated the kernel of NWChem Distributed Data SCF Algorithm, and
found that there will be two main problems needed to be considered if we want
to choose High Performance Fortran as a possible successor to global
arrays codes. The first is how to define the arrays and distributed
among the available processors. The second is how to implement the parallel
tasks among multi-processors.

According to the DDSCF codes in Global Arrays, the original density matrices
need to be dynamicly created and irregularly distributed among the
processors using global array function {\em ga\_create\_irreg}. As in
HPF 1.0, there is no support for irregular data distribution, even not in
HPF 2.0. But there are two new distribution formats: {\em GEN\_BLOCK} and
{\em INDIRECT}, in the new approved extensions of HPF 2. We can use the
"generalized" block distribution, {\em GEN\_BLOCK}, to implement this in
DDSCF HPF codes. Because the arrays are dynamicly created and then
irregularly distributed among the processors, though we can define the
arrays as ALLOCATABLE arrays to dynamicly allocated, we still need the
DYNAMIC attribute and the REDISTRIBUTE statements to implement. DYNAMIC
and REDISTRIBUTE features have all been moved from the HPF 2.0 kernel to
approved extensions.

In order to rewrite the dynamic global array allocation in Global Array
codes into HPF 2 codes, we mainly need to rewrite the function wihch will
call subroutine {\em ga\_create\_irreg} into HPF statements.
That's in function {\em ga\_create\_atom\_blocked} in file
{\em ga\_atom\_blk.F}. Global Array codes as:

\begin{verbatim}
      ......
      status =  ga_create_irreg(MT_DBL, nbf, nbf, name,
     $     map1, nblock1, map2, nblock2, g_a)
      ......
\end{verbatim}

They need to be rewritten in HPF as :

\begin{verbatim}
      integer hpf_map1(max_nproc), hpf_map2(max_nproc)
      REAL, ALLOCATABLE (:, :) :: g_array
!HPF$ DYNAMIC :: g_array
      ......
      call map_array_format_change(map1, hpf_map1)
      call map_array_format_change(map2, hpf_map2)
      ALLOCATE (g_array(nbf, nbf))
!HPF$ REDISTRIBUTE (g_array(GEN_BLOCK(hpf_map1, hpf_map2)))
      ......
\end{verbatim}

In Global Array codes, {\em map(i)} is the low bound of the
global array in processor {\em i}, and subroutine
{\em map\_array\_format\_change}
is aimed to change it to {\em hpf\_map(i)} corresponding the element
number of the global array in processor {\em i}, which is the HPF
"generalized" BLOCK distribution mode needed. And another thing we
need to pay attention is that in the funcation call to
{\em ga\_create\_atom\_blocked}, there's only a name parameter, in HPF we
need to change it to a global array, given an array name as a parameter.

The second thing we need to focus on is how to implement the parallel
tasks of DDSCF in HPF. There are two ways to express the parallel
execution in HPF, one is to use FORALL, another is to use INDEPENDENT
directive.

FORALL statement is designed to implement the prallel array assignment
among the processors. There's strict limitations in FORALL statement
and constructor. We can't call a subroutine in the FORALL loop body, and
the function calls in the loop body should be pure functions. In
the DDSCF codes, we truely
need the subroutine calls in the implementation of the parallel tasks.
So we have to use INDEPENDENT directive to implement.

The HPF INDEPENDENT directive allows the programmer to give information
to the compiler concerning oppertunities for parallel execution. HPF
allows a user to make assertions that no data object is defined by one
iteration of a DO loop and used (read or written) by another, and the
compiler may rely on them in its translation process to execute the
program in parallel.

In contrast to HPF 1.0, the INDEPENDENT assertion of HPF 2.0 allows
reductions to be performed in INDEPENDENT loops, provided the reduction
operator is a built-in, associative and commutative FORTRAN operator
(such as .AND.) or function (such as MAX). In this DDSCF kernel algorithm,
the data parallel computation cannot be expressed in HPF 1.0 as an
INDEPENDENT loop because several loop iterations update one or more
varibles. In such cases parallelism maybe possible and desirable
because the order of updates is immaterial to the final result.

In the implementation of the kernel algorithm of DDSCF, the Global
Array codes to implement the parallel tasks in subroutine
{\em fock\_2e\_a} in the file {\em fock\_2e\_a.F} is like:

\begin{verbatim}
      do ib = nblock, 1, -1
         iatlo = blocks(1,ib)
         iathi = blocks(2,ib)
         if (util_print('report loop', print_debug))
     $        call fock_2e_report(iathi, natoms)
         do jb = 1, ib
            jatlo = blocks(1,jb)
            jathi = blocks(2,jb)
            sij = schwarz_atom_block(iatlo,iathi,jatlo,jathi)
*
*     Cannot use symmetry of pairs since the blocking of atoms means
*     that we may have (kl) > (ij).
*
            if (sij*smax .ge. dentol) then
               do kb = ib, 1, -1
                  katlo = blocks(1,kb)
                  kathi = blocks(2,kb)
                  lbhi = kb
                  if (ib .eq. kb) lbhi = jb
                  do lb = 1, lbhi
                     latlo = blocks(1,lb)
                     lathi = blocks(2,lb)
                     skl = schwarz_atom_block(katlo,kathi,latlo,lathi)
                     otest = (sij*skl) .gt. dentol
                     if (otest .and. oskel) otest=sym_atom_block_quartet
     $                    (geom, iatlo, iathi, jatlo, jathi,
     $                       katlo, kathi, latlo, lathi)
c
c     Load balance over non-zero interactions
c
                     if (otest .and. (ijkl .eq. next)) then
                        call fock_2e_task(
     $                       geom, basis, oskel,
     $                       iatlo, jatlo, katlo, latlo,
     $                       iathi, jathi, kathi, lathi,
     $                       ijk_prev,
     $                       nfock, vg_dens, vg_fock,
     $                       jfac, kfac,
     $                       dij, dik, dli, djk, dlj, dlk,
     $                       fij, fik, fli, fjk, flj, flk,
     $                       tmp, tol2e, dentol)
c
                        next = nxtask(nproc, 1)
                     end if
                     if (otest) ijkl = ijkl + 1
                  end do
               enddo
            endif
         end do
      end do
\end{verbatim}

This kernel algorithm should be rewritten in HPF 2.0 like:

\begin{verbatim}
!HPF$ INDEPENDENT, NEW(jb, iatlo, iathi),REDUCTION(fij, fik, fli, fjk, flj,
flk)

      do ib = nblock, 1, -1
         iatlo = blocks(1,ib)
         iathi = blocks(2,ib)
         if (util_print('report loop', print_debug))
     $        call fock_2e_report(iathi, natoms)

!HPF$ INDEPENDENT, NEW(kb, jatlo, jathi, sij)

         do jb = 1, ib
            jatlo = blocks(1,jb)
            jathi = blocks(2,jb)
            sij = schwarz_atom_block(iatlo,iathi,jatlo,jathi)
*
*     Cannot use symmetry of pairs since the blocking of atoms means
*     that we may have (kl) > (ij).
*
            if (sij*smax .ge. dentol) then

!HPF$ INDEPENDENT, NEW(lb, katlo, kathi, lbhi)

               do kb = ib, 1, -1
                  katlo = blocks(1,kb)
                  kathi = blocks(2,kb)
                  lbhi = kb
                  if (ib .eq. kb) lbhi = jb

!HPF$ INDEPENDENT, NEW(latlo, lathi, skl, otest, next, ijkl)

                  do lb = 1, lbhi
                     latlo = blocks(1,lb)
                     lathi = blocks(2,lb)
                     skl = schwarz_atom_block(katlo,kathi,latlo,lathi)
                     otest = (sij*skl) .gt. dentol
                     if (otest .and. oskel) otest=sym_atom_block_quartet
     $                    (geom, iatlo, iathi, jatlo, jathi,
     $                       katlo, kathi, latlo, lathi)
c
c     Load balance over non-zero interactions
c
                     if (otest .and. (ijkl .eq. next)) then
                        call fock_2e_task(
     $                       geom, basis, oskel,
     $                       iatlo, jatlo, katlo, latlo,
     $                       iathi, jathi, kathi, lathi,
     $                       ijk_prev,
     $                       nfock, vg_dens, vg_fock,
     $                       jfac, kfac,
     $                       dij, dik, dli, djk, dlj, dlk,
     $                       fij, fik, fli, fjk, flj, flk,
     $                       tmp, tol2e, dentol)
c
                        next = nxtask(nproc, 1)
                     end if
                     if (otest) ijkl = ijkl + 1
                  end do
               enddo
            endif
         end do
      end do
\end{verbatim}

The remaining thing mainly to do is to replace the Global Array
functions in the DDSCF codes in HPF way. According to the kernel
algorithm, the codes needed to be rewritten are mainly in function
{\em uhn\_energy}  in file {\em uhf.F}. These Global Array functions
can be replaced in HPF in the following, respectively.

\begin{verbatim}
ga_dgemm(transa, transb, m, n, k, alpha, g_a, g_b, beta, g_c )
\end{verbatim}

can be replaced like:

\begin{verbatim}
FORALL (i=1:k, j=1:k) d(i,j) = alpha * OPA(1,j) * OPB(i,j)
FORALL (i=1:m, j=1:n) c(i,j) = d(i,j) + beta * C(i,j)
\end{verbatim}

where OPA = A or TRANSPOSE(A), OPB = B or TRANSPOSE(B).

\begin{verbatim}
ga_zero(A)
\end{verbatim}

can be replaced like:

\begin{verbatim}
A = 0
\end{verbatim}

\begin{verbatim}
ga_copy(A, B)
\end{verbatim}

can be replaced like:

\begin{verbatim}
A = B
\end{verbatim}

There are two main functions needed to be replaced when implementing
this DDSCF kernel algorithm. One is function
{\em fock\_get\_blk}  and {\em get\_atom\_block}
in file {\em fock\_get\_blk.F}, where to call the Global
Array function {\em ga\_get} to fetch the original computing data from
global arrays into local arrays before the computing.

Function {\em ga\_get} can be rewritten in HPF as:

\begin{verbatim}
      dimension g_a(1:dim1,1:dim2)
      dimension buf(1:ld, 1:*)

      buf(1:ihi-ilo+1, 1:jhi-jlo+1) = g_a(ilo:ihi, jlo:jhi)
\end{verbatim}

Another is function {\em fock\_upd\_blk} and {\em upd\_atom\_block}
in file {\em fock\_upd\_blk.F}, where  to call the Global Array
function {\em ga\_acc} to store the result from local arrays to the
global arrays.

Function {\em ga\_acc} can be rewritten in HPF as:

\begin{verbatim}
      dimension g_a(1:dim1,1:dim2)
      dimension buf(1:ld, 1:*)

      g_a(ilo:ihi, jlo:jhi) = g_a(ilo:ihi, jlo:jhi) +
                              alpha*buf(1:ihi-ilo+1, 1:jhi-jlo+1)
\end{verbatim}

\emph{Say something about why HPF 1 can't handle this}

\emph{Say something about new HPF 2 capabilities that you considered}

\emph{Say something about what you think is the best approach and
why.  Sketch code required to implement.}

\emph{If ``standard extensions'' are required, please discuss how
likely they are to be implemented (say something about utility to
other application domains if you know).}

\emph{Is there anything relevant in the standard which you fell is
poen to interpretation and may effect the ability or performance of
your proposed implementation? Talk about any modifications or
additional features you might propose adding to the language to
improve performance, ease of programming, etc.  (Anything we should
comment on in the public comment period?)}

\subsection{NWChem MP2 Back-Transformation Algorithm}

We investigated the Global Array codes of NWChem MP2
Back-Transformation Algorithm and found that this algorithm can
also be rewritten in HPF 2.0. And it will be easier than that in
DDSCF algorithm, for the simple calling structure in the kernel loop.
We can just use FORALL statemants and constructors to implement the
parallel tasks in this algorithm, instead of using INDEPENDENT
directive.

The main things we need to focus on when rewriting the
codes will also be in the two area, dynamic array allocation and
distribution, and the implementation of parallel tasks. First is the
dynamic array allcation. What we need to replace is the function call
to {\em ga\_create\_irreg} in subroutine {\em mp2\_nonsep\_uhf} in file
{\em mp2\_back\_transform.F}. We can use the same strategy we used in
the DDSCF codes to deal with this irregular distribution, and dynamic
array allocation problem.

After the set up of parallel tasks of {\em mp2\_backt\_info} function
call and the dynamic array allocation, the remaining thing is how to
implement them in parallel way in HPF. Then we need to rewrite the
function {\em mp2\_back\_transform\_uhf} in file
{\em mp2\_back\_transform.F} in to HPF codes. The main points we need
to focus on is to rewrite the following Global Array functions into
HPF codes with FORALL statements.

\begin{verbatim}
      call ga_fill_patch(g_buf, 1, nva, 1, nbfpair, 99.0d0)
\end{verbatim}

will be changed to :

\begin{verbatim}
      FORALL (i=1:nva, j=1:nbfpair) ga_buf(i,j) = 99.0d0
\end{verbatim}

Function {\em dgemm} can be treated as the same way we discussed in the
last section. And another Global Array function

\begin{verbatim}
            call ga_put(g_buf, a-nva_lo+1, a-nva_lo+1,
     $           1, nbfpair, tmp, 1)
\end{verbatim}

to update
the global result array after the computing can be replaced with:
will be replaced with HPF codes like:

\begin{verbatim}
      dimension g_buf(1:dim1,1:dim2)
      dimension tmp(1:1, 1:*)

      g_a(a-nva_lo+1:a-nva_lo+1, 1:nbfpair) = buf(1:1, 1:nbfpair)
\end{verbatim}

\section{Conclusions and Recommendations}

\end{document}

