text:
C-----------------------------------------------------------------------
C AMD:  approximate minimum degree, with aggressive absorption
C-----------------------------------------------------------------------

        SUBROUTINE AMD64
     $          (N, PE, IW, LEN, IWLEN, PFREE, NV, NEXT,
     $          LAST, HEAD, ELEN, DEGREE, NCMPA, W)

        USE,INTRINSIC :: ISO_FORTRAN_ENV
        INTEGER(INT64) N, IWLEN, PFREE, NCMPA, IW (IWLEN), PE (N),
     $              DEGREE (N), NV (N), NEXT (N), LAST (N), HEAD (N),
     $              ELEN (N), W (N), LEN (N)

C Given a representation of the nonzero pattern of a symmetric matrix,
C       A, (excluding the diagonal) perform an approximate minimum
C       (UMFPACK/MA38-style) degree ordering to compute a pivot order
C       such that the introduction of nonzeros (fill-in) in the Cholesky
C       factors A = LL^T are kept low.  At each step, the pivot
C       selected is the one with the minimum UMFPACK/MA38-style
C       upper-bound on the external degree.
C
C       Aggresive absorption is used to tighten the bound on the degree.

C **********************************************************************
C ***** CAUTION:  ARGUMENTS ARE NOT CHECKED FOR ERRORS ON INPUT.  ******
C **********************************************************************

C       References:
C
C       [1] Timothy A. Davis and Iain Duff, "An unsymmetric-pattern
C           multifrontal method for sparse LU factorization", SIAM J.
C           Matrix Analysis and Applications, vol. 18, no. 1, pp.
C           140-158.  Discusses UMFPACK / MA38, which first introduced
C           the approximate minimum degree used by this routine.
C
C       [2] Patrick Amestoy, Timothy A. Davis, and Iain S. Duff, "An
C           approximate degree ordering algorithm," SIAM J. Matrix
C           Analysis and Applications, vol. 17, no. 4, pp. 886-905,
C           1996.  Discusses AMD, AMDBAR, and MC47B.
C
C       [3] Alan George and Joseph Liu, "The evolution of the minimum
C           degree ordering algorithm," SIAM Review, vol. 31, no. 1,
C           pp. 1-19, 1989.  We list below the features mentioned in
C           that paper that this code includes:
C
C       mass elimination:
C               Yes.  MA27 relied on supervariable detection for mass
C               elimination.
C       indistinguishable nodes:
C               Yes (we call these "supervariables").  This was also in
C               the MA27 code - although we modified the method of
C               detecting them (the previous hash was the true degree,
C               which we no longer keep track of).  A supervariable is
C               a set of rows with identical nonzero pattern.  All
C               variables in a supervariable are eliminated together.
C               Each supervariable has as its numerical name that of
C               one of its variables (its principal variable).
C       quotient graph representation:
C               Yes.  We use the term "element" for the cliques formed
C               during elimination.  This was also in the MA27 code.
C               The algorithm can operate in place, but it will work
C               more efficiently if given some "elbow room."
C       element absorption:
C               Yes.  This was also in the MA27 code.
C       external degree:
C               Yes.  The MA27 code was based on the true degree.
C       incomplete degree update and multiple elimination:
C               No.  This was not in MA27, either.  Our method of
C               degree update within MC47B/BD is element-based, not
C               variable-based.  It is thus not well-suited for use
C               with incomplete degree update or multiple elimination.

C-----------------------------------------------------------------------
C Authors, and Copyright (C) 1995 by:
C       Timothy A. Davis, Patrick Amestoy, Iain S. Duff, & John K. Reid.
C
C Acknowledgements:
C       This work (and the UMFPACK package) was supported by the
C       National Science Foundation (ASC-9111263 and DMS-9223088).
C       The UMFPACK/MA38 approximate degree update algorithm, the
C       unsymmetric analog which forms the basis of MC47B/BD, was
C       developed while Tim Davis was supported by CERFACS (Toulouse,
C       France) in a post-doctoral position.
C
C Date:  September, 1995
C-----------------------------------------------------------------------

C-----------------------------------------------------------------------
C INPUT ARGUMENTS (unaltered):
C-----------------------------------------------------------------------

C n:    The matrix order.
C
C       Restriction:  1 .le. n .lt. (iovflo/2)-2, where iovflo is
C       the largest positive integer that your computer can represent.

C iwlen:        The length of iw (1..iwlen).  On input, the matrix is
C       stored in iw (1..pfree-1).  However, iw (1..iwlen) should be
C       slightly larger than what is required to hold the matrix, at
C       least iwlen .ge. pfree + n is recommended.  Otherwise,
C       excessive compressions will take place.
C       *** We do not recommend running this algorithm with ***
C       ***      iwlen .lt. pfree + n.                      ***
C       *** Better performance will be obtained if          ***
C       ***      iwlen .ge. pfree + n                       ***
C       *** or better yet                                   ***
C       ***      iwlen .gt. 1.2 * pfree                     ***
C       *** (where pfree is its value on input).            ***
C       The algorithm will not run at all if iwlen .lt. pfree-1.
C
C       Restriction: iwlen .ge. pfree-1

C-----------------------------------------------------------------------
C INPUT/OUPUT ARGUMENTS:
C-----------------------------------------------------------------------

C pe:   On input, pe (i) is the index in iw of the start of row i, or
C       zero if row i has no off-diagonal non-zeros.
C
C       During execution, it is used for both supervariables and
C       elements:
C
C       * Principal supervariable i:  index into iw of the
C               description of supervariable i.  A supervariable
C               represents one or more rows of the matrix
C               with identical nonzero pattern.
C       * Non-principal supervariable i:  if i has been absorbed
C               into another supervariable j, then pe (i) = -j.
C               That is, j has the same pattern as i.
C               Note that j might later be absorbed into another
C               supervariable j2, in which case pe (i) is still -j,
C               and pe (j) = -j2.
C       * Unabsorbed element e:  the index into iw of the description
C               of element e, if e has not yet been absorbed by a
C               subsequent element.  Element e is created when
C               the supervariable of the same name is selected as
C               the pivot.
C       * Absorbed element e:  if element e is absorbed into element
C               e2, then pe (e) = -e2.  This occurs when the pattern of
C               e (that is, Le) is found to be a subset of the pattern
C               of e2 (that is, Le2).  If element e is "null" (it has
C               no nonzeros outside its pivot block), then pe (e) = 0.
C
C       On output, pe holds the assembly tree/forest, which implicitly
C       represents a pivot order with identical fill-in as the actual
C       order (via a depth-first search of the tree).
C
C       On output:
C       If nv (i) .gt. 0, then i represents a node in the assembly tree,
C       and the parent of i is -pe (i), or zero if i is a root.
C       If nv (i) = 0, then (i,-pe (i)) represents an edge in a
C       subtree, the root of which is a node in the assembly tree.

C pfree:        On input the tail end of the array, iw (pfree..iwlen),
C       is empty, and the matrix is stored in iw (1..pfree-1).
C       During execution, additional data is placed in iw, and pfree
C       is modified so that iw (pfree..iwlen) is always the unused part
C       of iw.  On output, pfree is set equal to the size of iw that
C       would have been needed for no compressions to occur.  If
C       ncmpa is zero, then pfree (on output) is less than or equal to
C       iwlen, and the space iw (pfree+1 ... iwlen) was not used.
C       Otherwise, pfree (on output) is greater than iwlen, and all the
C       memory in iw was used.

C-----------------------------------------------------------------------
C INPUT/MODIFIED (undefined on output):
C-----------------------------------------------------------------------

C len:  On input, len (i) holds the number of entries in row i of the
C       matrix, excluding the diagonal.  The contents of len (1..n)
C       are undefined on output.

C iw:   On input, iw (1..pfree-1) holds the description of each row i
C       in the matrix.  The matrix must be symmetric, and both upper
C       and lower triangular parts must be present.  The diagonal must
C       not be present.  Row i is held as follows:
C
C               len (i):  the length of the row i data structure
C               iw (pe (i) ... pe (i) + len (i) - 1):
C                       the list of column indices for nonzeros
C                       in row i (simple supervariables), excluding
C                       the diagonal.  All supervariables start with
C                       one row/column each (supervariable i is just
C                       row i).
C               if len (i) is zero on input, then pe (i) is ignored
C               on input.
C
C               Note that the rows need not be in any particular order,
C               and there may be empty space between the rows.
C
C       During execution, the supervariable i experiences fill-in.
C       This is represented by placing in i a list of the elements
C       that cause fill-in in supervariable i:
C
C               len (i):  the length of supervariable i
C               iw (pe (i) ... pe (i) + elen (i) - 1):
C                       the list of elements that contain i.  This list
C                       is kept short by removing absorbed elements.
C               iw (pe (i) + elen (i) ... pe (i) + len (i) - 1):
C                       the list of supervariables in i.  This list
C                       is kept short by removing nonprincipal
C                       variables, and any entry j that is also
C                       contained in at least one of the elements
C                       (j in Le) in the list for i (e in row i).
C
C       When supervariable i is selected as pivot, we create an
C       element e of the same name (e=i):
C
C               len (e):  the length of element e
C               iw (pe (e) ... pe (e) + len (e) - 1):
C                       the list of supervariables in element e.
C
C       An element represents the fill-in that occurs when supervariable
C       i is selected as pivot (which represents the selection of row i
C       and all non-principal variables whose principal variable is i).
C       We use the term Le to denote the set of all supervariables
C       in element e.  Absorbed supervariables and elements are pruned
C       from these lists when computationally convenient.
C
C       CAUTION:  THE INPUT MATRIX IS OVERWRITTEN DURING COMPUTATION.
C       The contents of iw are undefined on output.

C-----------------------------------------------------------------------
C OUTPUT (need not be set on input):
C-----------------------------------------------------------------------

C nv:   During execution, abs (nv (i)) is equal to the number of rows
C       that are represented by the principal supervariable i.  If i is
C       a nonprincipal variable, then nv (i) = 0.  Initially,
C       nv (i) = 1 for all i.  nv (i) .lt. 0 signifies that i is a
C       principal variable in the pattern Lme of the current pivot
C       element me.  On output, nv (e) holds the true degree of element
C       e at the time it was created (including the diagonal part).

C ncmpa:        The number of times iw was compressed.  If this is
C       excessive, then the execution took longer than what could have
C       been.  To reduce ncmpa, try increasing iwlen to be 10% or 20%
C       larger than the value of pfree on input (or at least
C       iwlen .ge. pfree + n).  The fastest performance will be
C       obtained when ncmpa is returned as zero.  If iwlen is set to
C       the value returned by pfree on *output*, then no compressions
C       will occur.

C elen: See the description of iw above.  At the start of execution,
C       elen (i) is set to zero.  During execution, elen (i) is the
C       number of elements in the list for supervariable i.  When e
C       becomes an element, elen (e) = -nel is set, where nel is the
C       current step of factorization.  elen (i) = 0 is done when i
C       becomes nonprincipal.
C
C       For variables, elen (i) .ge. 0 holds until just before the
C       permutation vectors are computed.  For elements,
C       elen (e) .lt. 0 holds.
C
C       On output elen (1..n) holds the inverse permutation (the same
C       as the 'INVP' argument in Sparspak).  That is, if k = elen (i),
C       then row i is the kth pivot row.  Row i of A appears as the
C       (elen(i))-th row in the permuted matrix, PAP^T.

C last: In a degree list, last (i) is the supervariable preceding i,
C       or zero if i is the head of the list.  In a hash bucket,
C       last (i) is the hash key for i.  last (head (hash)) is also
C       used as the head of a hash bucket if head (hash) contains a
C       degree list (see head, below).
C
C       On output, last (1..n) holds the permutation (the same as the
C       'PERM' argument in Sparspak).  That is, if i = last (k), then
C       row i is the kth pivot row.  Row last (k) of A is the k-th row
C       in the permuted matrix, PAP^T.

C-----------------------------------------------------------------------
C LOCAL (not input or output - used only during execution):
C-----------------------------------------------------------------------

C degree:       If i is a supervariable, then degree (i) holds the
C       current approximation of the external degree of row i (an upper
C       bound).  The external degree is the number of nonzeros in row i,
C       minus abs (nv (i)) (the diagonal part).  The bound is equal to
C       the external degree if elen (i) is less than or equal to two.
C
C       We also use the term "external degree" for elements e to refer
C       to |Le \ Lme|.  If e is an element, then degree (e) holds |Le|,
C       which is the degree of the off-diagonal part of the element e
C       (not including the diagonal part).

C head: head is used for degree lists.  head (deg) is the first
C       supervariable in a degree list (all supervariables i in a
C       degree list deg have the same approximate degree, namely,
C       deg = degree (i)).  If the list deg is empty then
C       head (deg) = 0.
C
C       During supervariable detection head (hash) also serves as a
C       pointer to a hash bucket.
C       If head (hash) .gt. 0, there is a degree list of degree hash.
C               The hash bucket head pointer is last (head (hash)).
C       If head (hash) = 0, then the degree list and hash bucket are
C               both empty.
C       If head (hash) .lt. 0, then the degree list is empty, and
C               -head (hash) is the head of the hash bucket.
C       After supervariable detection is complete, all hash buckets
C       are empty, and the (last (head (hash)) = 0) condition is
C       restored for the non-empty degree lists.

C next: next (i) is the supervariable following i in a link list, or
C       zero if i is the last in the list.  Used for two kinds of
C       lists:  degree lists and hash buckets (a supervariable can be
C       in only one kind of list at a time).

C w:    The flag array w determines the status of elements and
C       variables, and the external degree of elements.
C
C       for elements:
C          if w (e) = 0, then the element e is absorbed
C          if w (e) .ge. wflg, then w (e) - wflg is the size of
C               the set |Le \ Lme|, in terms of nonzeros (the
C               sum of abs (nv (i)) for each principal variable i that
C               is both in the pattern of element e and NOT in the
C               pattern of the current pivot element, me).
C          if wflg .gt. w (e) .gt. 0, then e is not absorbed and has
C               not yet been seen in the scan of the element lists in
C               the computation of |Le\Lme| in loop 150 below.
C
C       for variables:
C          during supervariable detection, if w (j) .ne. wflg then j is
C          not in the pattern of variable i
C
C       The w array is initialized by setting w (i) = 1 for all i,
C       and by setting wflg = 2.  It is reinitialized if wflg becomes
C       too large (to ensure that wflg+n does not cause integer
C       overflow).

C-----------------------------------------------------------------------
C LOCAL INTEGERS:
C-----------------------------------------------------------------------

        INTEGER(INT64)
     $          DEG, DEGME, DEXT, DMAX, E, ELENME, ELN, HASH, HMOD, I,
     $          ILAST, INEXT, J, JLAST, JNEXT, K, KNT1, KNT2, KNT3,
     $          LENJ, LN, MAXMEM, ME, MEM, MINDEG, NEL, NEWMEM,
     $          NLEFT, NVI, NVJ, NVPIV, SLENME, WE, WFLG, WNVI, X

C deg:          the degree of a variable or element
C degme:        size, |Lme|, of the current element, me (= degree (me))
C dext:         external degree, |Le \ Lme|, of some element e
C dmax:         largest |Le| seen so far
C e:            an element
C elenme:       the length, elen (me), of element list of pivotal var.
C eln:          the length, elen (...), of an element list
C hash:         the computed value of the hash function
C hmod:         the hash function is computed modulo hmod = max (1,n-1)
C i:            a supervariable
C ilast:        the entry in a link list preceding i
C inext:        the entry in a link list following i
C j:            a supervariable
C jlast:        the entry in a link list preceding j
C jnext:        the entry in a link list, or path, following j
C k:            the pivot order of an element or variable
C knt1:         loop counter used during element construction
C knt2:         loop counter used during element construction
C knt3:         loop counter used during compression
C lenj:         len (j)
C ln:           length of a supervariable list
C maxmem:       amount of memory needed for no compressions
C me:           current supervariable being eliminated, and the
C                       current element created by eliminating that
C                       supervariable
C mem:          memory in use assuming no compressions have occurred
C mindeg:       current minimum degree
C nel:          number of pivots selected so far
C newmem:       amount of new memory needed for current pivot element
C nleft:        n - nel, the number of nonpivotal rows/columns remaining
C nvi:          the number of variables in a supervariable i (= nv (i))
C nvj:          the number of variables in a supervariable j (= nv (j))
C nvpiv:        number of pivots in current element
C slenme:       number of variables in variable list of pivotal variable
C we:           w (e)
C wflg:         used for flagging the w array.  See description of iw.
C wnvi:         wflg - nv (i)
C x:            either a supervariable or an element

C-----------------------------------------------------------------------
C LOCAL POINTERS:
C-----------------------------------------------------------------------

        INTEGER(INT64)
     $          P, P1, P2, P3, PDST, PEND, PJ, PME, PME1, PME2, PN, PSRC

C               Any parameter (pe (...) or pfree) or local variable
C               starting with "p" (for Pointer) is an index into iw,
C               and all indices into iw use variables starting with
C               "p."  The only exception to this rule is the iwlen
C               input argument.

C p:            pointer into lots of things
C p1:           pe (i) for some variable i (start of element list)
C p2:           pe (i) + elen (i) -  1 for some var. i (end of el. list)
C p3:           index of first supervariable in clean list
C pdst:         destination pointer, for compression
C pend:         end of memory to compress
C pj:           pointer into an element or variable
C pme:          pointer into the current element (pme1...pme2)
C pme1:         the current element, me, is stored in iw (pme1...pme2)
C pme2:         the end of the current element
C pn:           pointer into a "clean" variable, also used to compress
C psrc:         source pointer, for compression

C-----------------------------------------------------------------------
C  FUNCTIONS CALLED:
C-----------------------------------------------------------------------

        INTRINSIC MAX, MIN, MOD

C=======================================================================
C  INITIALIZATIONS
C=======================================================================

        WFLG = 2
        MINDEG = 1
        NCMPA = 0
        NEL = 0
        HMOD = MAX (1, N-1)
        DMAX = 0
        MEM = PFREE - 1
        MAXMEM = MEM
	ME = 0

        DO 10 I = 1, N
           LAST (I) = 0
           HEAD (I) = 0
           NV (I) = 1
           W (I) = 1
           ELEN (I) = 0
           DEGREE (I) = LEN (I)
10         CONTINUE

C       ----------------------------------------------------------------
C       initialize degree lists and eliminate rows with no off-diag. nz.
C       ----------------------------------------------------------------

        DO 20 I = 1, N

           DEG = DEGREE (I)

           IF (DEG .GT. 0) THEN

C             ----------------------------------------------------------
C             place i in the degree list corresponding to its degree
C             ----------------------------------------------------------

              INEXT = HEAD (DEG)
              IF (INEXT .NE. 0) LAST (INEXT) = I
              NEXT (I) = INEXT
              HEAD (DEG) = I

           ELSE

C             ----------------------------------------------------------
C             we have a variable that can be eliminated at once because
C             there is no off-diagonal non-zero in its row.
C             ----------------------------------------------------------

              NEL = NEL + 1
              ELEN (I) = -NEL
              PE (I) = 0
              W (I) = 0

              ENDIF

20         CONTINUE

C=======================================================================
C  WHILE (selecting pivots) DO
C=======================================================================

30      CONTINUE
        IF (NEL .LT. N) THEN

C=======================================================================
C  GET PIVOT OF MINIMUM DEGREE
C=======================================================================

C          -------------------------------------------------------------
C          find next supervariable for elimination
C          -------------------------------------------------------------

           DO 40 DEG = MINDEG, N
              ME = HEAD (DEG)
              IF (ME .GT. 0) GOTO 50
40            CONTINUE
50         CONTINUE
           MINDEG = DEG

C          -------------------------------------------------------------
C          remove chosen variable from link list
C          -------------------------------------------------------------

           INEXT = NEXT (ME)
           IF (INEXT .NE. 0) LAST (INEXT) = 0
           HEAD (DEG) = INEXT

C          -------------------------------------------------------------
C          me represents the elimination of pivots nel+1 to nel+nv(me).
C          place me itself as the first in this set.  It will be moved
C          to the nel+nv(me) position when the permutation vectors are
C          computed.
C          -------------------------------------------------------------

           ELENME = ELEN (ME)
           ELEN (ME) = - (NEL + 1)
           NVPIV = NV (ME)
           NEL = NEL + NVPIV

C=======================================================================
C  CONSTRUCT NEW ELEMENT
C=======================================================================

C          -------------------------------------------------------------
C          At this point, me is the pivotal supervariable.  It will be
C          converted into the current element.  Scan list of the
C          pivotal supervariable, me, setting tree pointers and
C          constructing new list of supervariables for the new element,
C          me.  p is a pointer to the current position in the old list.
C          -------------------------------------------------------------

C          flag the variable "me" as being in Lme by negating nv (me)
           NV (ME) = -NVPIV
           DEGME = 0

           IF (ELENME .EQ. 0) THEN

C             ----------------------------------------------------------
C             construct the new element in place
C             ----------------------------------------------------------

              PME1 = PE (ME)
              PME2 = PME1 - 1

              DO 60 P = PME1, PME1 + LEN (ME) - 1
                 I = IW (P)
                 NVI = NV (I)
                 IF (NVI .GT. 0) THEN

C                   ----------------------------------------------------
C                   i is a principal variable not yet placed in Lme.
C                   store i in new list
C                   ----------------------------------------------------

                    DEGME = DEGME + NVI
C                   flag i as being in Lme by negating nv (i)
                    NV (I) = -NVI
                    PME2 = PME2 + 1
                    IW (PME2) = I

C                   ----------------------------------------------------
C                   remove variable i from degree list.
C                   ----------------------------------------------------

                    ILAST = LAST (I)
                    INEXT = NEXT (I)
                    IF (INEXT .NE. 0) LAST (INEXT) = ILAST
                    IF (ILAST .NE. 0) THEN
                       NEXT (ILAST) = INEXT
                    ELSE
C                      i is at the head of the degree list
                       HEAD (DEGREE (I)) = INEXT
                       ENDIF

                    ENDIF
60               CONTINUE
C             this element takes no new memory in iw:
              NEWMEM = 0

           ELSE

C             ----------------------------------------------------------
C             construct the new element in empty space, iw (pfree ...)
C             ----------------------------------------------------------

              P = PE (ME)
              PME1 = PFREE
              SLENME = LEN (ME) - ELENME

              DO 120 KNT1 = 1, ELENME + 1

                 IF (KNT1 .GT. ELENME) THEN
C                   search the supervariables in me.
                    E = ME
                    PJ = P
                    LN = SLENME
                 ELSE
C                   search the elements in me.
                    E = IW (P)
                    P = P + 1
                    PJ = PE (E)
                    LN = LEN (E)
                    ENDIF

C                -------------------------------------------------------
C                search for different supervariables and add them to the
C                new list, compressing when necessary. this loop is
C                executed once for each element in the list and once for
C                all the supervariables in the list.
C                -------------------------------------------------------

                 DO 110 KNT2 = 1, LN
                    I = IW (PJ)
                    PJ = PJ + 1
                    NVI = NV (I)
                    IF (NVI .GT. 0) THEN

C                      -------------------------------------------------
C                      compress iw, if necessary
C                      -------------------------------------------------

                       IF (PFREE .GT. IWLEN) THEN
C                         prepare for compressing iw by adjusting
C                         pointers and lengths so that the lists being
C                         searched in the inner and outer loops contain
C                         only the remaining entries.

                          PE (ME) = P
                          LEN (ME) = LEN (ME) - KNT1
                          IF (LEN (ME) .EQ. 0) THEN
C                            nothing left of supervariable me
                             PE (ME) = 0
                             ENDIF
                          PE (E) = PJ
                          LEN (E) = LN - KNT2
                          IF (LEN (E) .EQ. 0) THEN
C                            nothing left of element e
                             PE (E) = 0
                             ENDIF

                          NCMPA = NCMPA + 1
C                         store first item in pe
C                         set first entry to -item
                          DO 70 J = 1, N
                             PN = PE (J)
                             IF (PN .GT. 0) THEN
                                PE (J) = IW (PN)
                                IW (PN) = -J
                                ENDIF
70                           CONTINUE

C                         psrc/pdst point to source/destination
                          PDST = 1
                          PSRC = 1
                          PEND = PME1 - 1

C                         while loop:
80                        CONTINUE
                          IF (PSRC .LE. PEND) THEN
C                            search for next negative entry
                             J = -IW (PSRC)
                             PSRC = PSRC + 1
                             IF (J .GT. 0) THEN
                                IW (PDST) = PE (J)
                                PE (J) = PDST
                                PDST = PDST + 1
C                               copy from source to destination
                                LENJ = LEN (J)
                                DO 90 KNT3 = 0, LENJ - 2
                                   IW (PDST + KNT3) = IW (PSRC + KNT3)
90                                 CONTINUE
                                PDST = PDST + LENJ - 1
                                PSRC = PSRC + LENJ - 1
                                ENDIF
                             GOTO 80
                             ENDIF

C                         move the new partially-constructed element
                          P1 = PDST
                          DO 100 PSRC = PME1, PFREE - 1
                             IW (PDST) = IW (PSRC)
                             PDST = PDST + 1
100                          CONTINUE
                          PME1 = P1
                          PFREE = PDST
                          PJ = PE (E)
                          P = PE (ME)
                          ENDIF

C                      -------------------------------------------------
C                      i is a principal variable not yet placed in Lme
C                      store i in new list
C                      -------------------------------------------------

                       DEGME = DEGME + NVI
C                      flag i as being in Lme by negating nv (i)
                       NV (I) = -NVI
                       IW (PFREE) = I
                       PFREE = PFREE + 1

C                      -------------------------------------------------
C                      remove variable i from degree link list
C                      -------------------------------------------------

                       ILAST = LAST (I)
                       INEXT = NEXT (I)
                       IF (INEXT .NE. 0) LAST (INEXT) = ILAST
                       IF (ILAST .NE. 0) THEN
                          NEXT (ILAST) = INEXT
                       ELSE
C                         i is at the head of the degree list
                          HEAD (DEGREE (I)) = INEXT
                          ENDIF

                       ENDIF
110                 CONTINUE

                 IF (E .NE. ME) THEN
C                   set tree pointer and flag to indicate element e is
C                   absorbed into new element me (the parent of e is me)
                    PE (E) = -ME
                    W (E) = 0
                    ENDIF
120              CONTINUE

              PME2 = PFREE - 1
C             this element takes newmem new memory in iw (possibly zero)
              NEWMEM = PFREE - PME1
              MEM = MEM + NEWMEM
              MAXMEM = MAX (MAXMEM, MEM)
              ENDIF

C          -------------------------------------------------------------
C          me has now been converted into an element in iw (pme1..pme2)
C          -------------------------------------------------------------

C          degme holds the external degree of new element
           DEGREE (ME) = DEGME
           PE (ME) = PME1
           LEN (ME) = PME2 - PME1 + 1

C          -------------------------------------------------------------
C          make sure that wflg is not too large.  With the current
C          value of wflg, wflg+n must not cause integer overflow
C          -------------------------------------------------------------

           IF (WFLG + N .LE. WFLG) THEN
              DO 130 X = 1, N
                 IF (W (X) .NE. 0) W (X) = 1
130              CONTINUE
              WFLG = 2
              ENDIF

C=======================================================================
C  COMPUTE (w (e) - wflg) = |Le\Lme| FOR ALL ELEMENTS
C=======================================================================

C          -------------------------------------------------------------
C          Scan 1:  compute the external degrees of previous elements
C          with respect to the current element.  That is:
C               (w (e) - wflg) = |Le \ Lme|
C          for each element e that appears in any supervariable in Lme.
C          The notation Le refers to the pattern (list of
C          supervariables) of a previous element e, where e is not yet
C          absorbed, stored in iw (pe (e) + 1 ... pe (e) + iw (pe (e))).
C          The notation Lme refers to the pattern of the current element
C          (stored in iw (pme1..pme2)).   If (w (e) - wflg) becomes
C          zero, then the element e will be absorbed in scan 2.
C          -------------------------------------------------------------

           DO 150 PME = PME1, PME2
              I = IW (PME)
              ELN = ELEN (I)
              IF (ELN .GT. 0) THEN
C                note that nv (i) has been negated to denote i in Lme:
                 NVI = -NV (I)
                 WNVI = WFLG - NVI
                 DO 140 P = PE (I), PE (I) + ELN - 1
                    E = IW (P)
                    WE = W (E)
                    IF (WE .GE. WFLG) THEN
C                      unabsorbed element e has been seen in this loop
                       WE = WE - NVI
                    ELSE IF (WE .NE. 0) THEN
C                      e is an unabsorbed element
C                      this is the first we have seen e in all of Scan 1
                       WE = DEGREE (E) + WNVI
                       ENDIF
                    W (E) = WE
140                 CONTINUE
                 ENDIF
150           CONTINUE

C=======================================================================
C  DEGREE UPDATE AND ELEMENT ABSORPTION
C=======================================================================

C          -------------------------------------------------------------
C          Scan 2:  for each i in Lme, sum up the degree of Lme (which
C          is degme), plus the sum of the external degrees of each Le
C          for the elements e appearing within i, plus the
C          supervariables in i.  Place i in hash list.
C          -------------------------------------------------------------

           DO 180 PME = PME1, PME2
              I = IW (PME)
              P1 = PE (I)
              P2 = P1 + ELEN (I) - 1
              PN = P1
              HASH = 0
              DEG = 0

C             ----------------------------------------------------------
C             scan the element list associated with supervariable i
C             ----------------------------------------------------------

              DO 160 P = P1, P2
                 E = IW (P)
C                dext = | Le \ Lme |
                 DEXT = W (E) - WFLG
                 IF (DEXT .GT. 0) THEN
                    DEG = DEG + DEXT
                    IW (PN) = E
                    PN = PN + 1
                    HASH = HASH + E
                 ELSE IF (DEXT .EQ. 0) THEN
C                   aggressive absorption: e is not adjacent to me, but
C                   the |Le \ Lme| is 0, so absorb it into me
                    PE (E) = -ME
                    W (E) = 0
                 ELSE
C                   element e has already been absorbed, due to
C                   regular absorption, in do loop 120 above. Ignore it.
                    CONTINUE
                    ENDIF
160              CONTINUE

C             count the number of elements in i (including me):
              ELEN (I) = PN - P1 + 1

C             ----------------------------------------------------------
C             scan the supervariables in the list associated with i
C             ----------------------------------------------------------

              P3 = PN
              DO 170 P = P2 + 1, P1 + LEN (I) - 1
                 J = IW (P)
                 NVJ = NV (J)
                 IF (NVJ .GT. 0) THEN
C                   j is unabsorbed, and not in Lme.
C                   add to degree and add to new list
                    DEG = DEG + NVJ
                    IW (PN) = J
                    PN = PN + 1
                    HASH = HASH + J
                    ENDIF
170              CONTINUE

C             ----------------------------------------------------------
C             update the degree and check for mass elimination
C             ----------------------------------------------------------

              IF (DEG .EQ. 0) THEN

C                -------------------------------------------------------
C                mass elimination
C                -------------------------------------------------------

C                There is nothing left of this node except for an
C                edge to the current pivot element.  elen (i) is 1,
C                and there are no variables adjacent to node i.
C                Absorb i into the current pivot element, me.

                 PE (I) = -ME
                 NVI = -NV (I)
                 DEGME = DEGME - NVI
                 NVPIV = NVPIV + NVI
                 NEL = NEL + NVI
                 NV (I) = 0
                 ELEN (I) = 0

              ELSE

C                -------------------------------------------------------
C                update the upper-bound degree of i
C                -------------------------------------------------------

C                the following degree does not yet include the size
C                of the current element, which is added later:
                 DEGREE (I) = MIN (DEGREE (I), DEG)

C                -------------------------------------------------------
C                add me to the list for i
C                -------------------------------------------------------

C                move first supervariable to end of list
                 IW (PN) = IW (P3)
C                move first element to end of element part of list
                 IW (P3) = IW (P1)
C                add new element to front of list.
                 IW (P1) = ME
C                store the new length of the list in len (i)
                 LEN (I) = PN - P1 + 1

C                -------------------------------------------------------
C                place in hash bucket.  Save hash key of i in last (i).
C                -------------------------------------------------------

                 HASH = MOD (HASH, HMOD) + 1
                 J = HEAD (HASH)
                 IF (J .LE. 0) THEN
C                   the degree list is empty, hash head is -j
                    NEXT (I) = -J
                    HEAD (HASH) = -I
                 ELSE
C                   degree list is not empty
C                   use last (head (hash)) as hash head
                    NEXT (I) = LAST (J)
                    LAST (J) = I
                    ENDIF
                 LAST (I) = HASH
                 ENDIF
180           CONTINUE

           DEGREE (ME) = DEGME

C          -------------------------------------------------------------
C          Clear the counter array, w (...), by incrementing wflg.
C          -------------------------------------------------------------

           DMAX = MAX (DMAX, DEGME)
           WFLG = WFLG + DMAX

C          make sure that wflg+n does not cause integer overflow
           IF (WFLG + N .LE. WFLG) THEN
              DO 190 X = 1, N
                 IF (W (X) .NE. 0) W (X) = 1
190              CONTINUE
              WFLG = 2
              ENDIF
C          at this point, w (1..n) .lt. wflg holds

C=======================================================================
C  SUPERVARIABLE DETECTION
C=======================================================================

           DO 250 PME = PME1, PME2
              I = IW (PME)
              IF (NV (I) .LT. 0) THEN
C                i is a principal variable in Lme

C                -------------------------------------------------------
C                examine all hash buckets with 2 or more variables.  We
C                do this by examing all unique hash keys for super-
C                variables in the pattern Lme of the current element, me
C                -------------------------------------------------------

                 HASH = LAST (I)
C                let i = head of hash bucket, and empty the hash bucket
                 J = HEAD (HASH)
                 IF (J .EQ. 0) GOTO 250
                 IF (J .LT. 0) THEN
C                   degree list is empty
                    I = -J
                    HEAD (HASH) = 0
                 ELSE
C                   degree list is not empty, restore last () of head
                    I = LAST (J)
                    LAST (J) = 0
                    ENDIF
                 IF (I .EQ. 0) GOTO 250

C                while loop:
200              CONTINUE
                 IF (NEXT (I) .NE. 0) THEN

C                   ----------------------------------------------------
C                   this bucket has one or more variables following i.
C                   scan all of them to see if i can absorb any entries
C                   that follow i in hash bucket.  Scatter i into w.
C                   ----------------------------------------------------

                    LN = LEN (I)
                    ELN = ELEN (I)
C                   do not flag the first element in the list (me)
                    DO 210 P = PE (I) + 1, PE (I) + LN - 1
                       W (IW (P)) = WFLG
210                    CONTINUE

C                   ----------------------------------------------------
C                   scan every other entry j following i in bucket
C                   ----------------------------------------------------

                    JLAST = I
                    J = NEXT (I)

C                   while loop:
220                 CONTINUE
                    IF (J .NE. 0) THEN

C                      -------------------------------------------------
C                      check if j and i have identical nonzero pattern
C                      -------------------------------------------------

                       IF (LEN (J) .NE. LN) THEN
C                         i and j do not have same size data structure
                          GOTO 240
                          ENDIF
                       IF (ELEN (J) .NE. ELN) THEN
C                         i and j do not have same number of adjacent el
                          GOTO 240
                          ENDIF
C                      do not flag the first element in the list (me)
                       DO 230 P = PE (J) + 1, PE (J) + LN - 1
                          IF (W (IW (P)) .NE. WFLG) THEN
C                            an entry (iw(p)) is in j but not in i
                             GOTO 240
                             ENDIF
230                       CONTINUE

C                      -------------------------------------------------
C                      found it!  j can be absorbed into i
C                      -------------------------------------------------

                       PE (J) = -I
C                      both nv (i) and nv (j) are negated since they
C                      are in Lme, and the absolute values of each
C                      are the number of variables in i and j:
                       NV (I) = NV (I) + NV (J)
                       NV (J) = 0
                       ELEN (J) = 0
C                      delete j from hash bucket
                       J = NEXT (J)
                       NEXT (JLAST) = J
                       GOTO 220

C                      -------------------------------------------------
240                    CONTINUE
C                      j cannot be absorbed into i
C                      -------------------------------------------------

                       JLAST = J
                       J = NEXT (J)
                       GOTO 220
                       ENDIF

C                   ----------------------------------------------------
C                   no more variables can be absorbed into i
C                   go to next i in bucket and clear flag array
C                   ----------------------------------------------------

                    WFLG = WFLG + 1
                    I = NEXT (I)
                    IF (I .NE. 0) GOTO 200
                    ENDIF
                 ENDIF
250           CONTINUE

C=======================================================================
C  RESTORE DEGREE LISTS AND REMOVE NONPRINCIPAL SUPERVAR. FROM ELEMENT
C=======================================================================

           P = PME1
           NLEFT = N - NEL
           DO 260 PME = PME1, PME2
              I = IW (PME)
              NVI = -NV (I)
              IF (NVI .GT. 0) THEN
C                i is a principal variable in Lme
C                restore nv (i) to signify that i is principal
                 NV (I) = NVI

C                -------------------------------------------------------
C                compute the external degree (add size of current elem)
C                -------------------------------------------------------

                 DEG = MIN (DEGREE (I) + DEGME - NVI, NLEFT - NVI)

C                -------------------------------------------------------
C                place the supervariable at the head of the degree list
C                -------------------------------------------------------

                 INEXT = HEAD (DEG)
                 IF (INEXT .NE. 0) LAST (INEXT) = I
                 NEXT (I) = INEXT
                 LAST (I) = 0
                 HEAD (DEG) = I

C                -------------------------------------------------------
C                save the new degree, and find the minimum degree
C                -------------------------------------------------------

                 MINDEG = MIN (MINDEG, DEG)
                 DEGREE (I) = DEG

C                -------------------------------------------------------
C                place the supervariable in the element pattern
C                -------------------------------------------------------

                 IW (P) = I
                 P = P + 1
                 ENDIF
260           CONTINUE

C=======================================================================
C  FINALIZE THE NEW ELEMENT
C=======================================================================

           NV (ME) = NVPIV + DEGME
C          nv (me) is now the degree of pivot (including diagonal part)
C          save the length of the list for the new element me
           LEN (ME) = P - PME1
           IF (LEN (ME) .EQ. 0) THEN
C             there is nothing left of the current pivot element
              PE (ME) = 0
              W (ME) = 0
              ENDIF
           IF (NEWMEM .NE. 0) THEN
C             element was not constructed in place: deallocate part
C             of it (final size is less than or equal to newmem,
C             since newly nonprincipal variables have been removed).
              PFREE = P
              MEM = MEM - NEWMEM + LEN (ME)
              ENDIF

C=======================================================================
C          END WHILE (selecting pivots)
           GOTO 30
           ENDIF
C=======================================================================

C=======================================================================
C  COMPUTE THE PERMUTATION VECTORS
C=======================================================================

C       ----------------------------------------------------------------
C       The time taken by the following code is O(n).  At this
C       point, elen (e) = -k has been done for all elements e,
C       and elen (i) = 0 has been done for all nonprincipal
C       variables i.  At this point, there are no principal
C       supervariables left, and all elements are absorbed.
C       ----------------------------------------------------------------

C       ----------------------------------------------------------------
C       compute the ordering of unordered nonprincipal variables
C       ----------------------------------------------------------------

        DO 290 I = 1, N
           IF (ELEN (I) .EQ. 0) THEN

C             ----------------------------------------------------------
C             i is an un-ordered row.  Traverse the tree from i until
C             reaching an element, e.  The element, e, was the
C             principal supervariable of i and all nodes in the path
C             from i to when e was selected as pivot.
C             ----------------------------------------------------------

              J = -PE (I)
C             while (j is a variable) do:
270           CONTINUE
              IF (ELEN (J) .GE. 0) THEN
                 J = -PE (J)
                 GOTO 270
                 ENDIF
              E = J

C             ----------------------------------------------------------
C             get the current pivot ordering of e
C             ----------------------------------------------------------

              K = -ELEN (E)

C             ----------------------------------------------------------
C             traverse the path again from i to e, and compress the
C             path (all nodes point to e).  Path compression allows
C             this code to compute in O(n) time.  Order the unordered
C             nodes in the path, and place the element e at the end.
C             ----------------------------------------------------------

              J = I
C             while (j is a variable) do:
280           CONTINUE
              IF (ELEN (J) .GE. 0) THEN
                 JNEXT = -PE (J)
                 PE (J) = -E
                 IF (ELEN (J) .EQ. 0) THEN
C                   j is an unordered row
                    ELEN (J) = K
                    K = K + 1
                    ENDIF
                 J = JNEXT
                 GOTO 280
                 ENDIF
C             leave elen (e) negative, so we know it is an element
              ELEN (E) = -K
              ENDIF
290        CONTINUE

C       ----------------------------------------------------------------
C       reset the inverse permutation (elen (1..n)) to be positive,
C       and compute the permutation (last (1..n)).
C       ----------------------------------------------------------------

        DO 300 I = 1, N
           K = ABS (ELEN (I))
           LAST (K) = I
           ELEN (I) = K
300        CONTINUE

C=======================================================================
C  RETURN THE MEMORY USAGE IN IW
C=======================================================================

C       If maxmem is less than or equal to iwlen, then no compressions
C       occurred, and iw (maxmem+1 ... iwlen) was unused.  Otherwise
C       compressions did occur, and iwlen would have had to have been
C       greater than or equal to maxmem for no compressions to occur.
C       Return the value of maxmem in the pfree argument.

        PFREE = MAXMEM

        RETURN
        END

max_stars_repo_name: masuday/fsparse

max_issues_repo_name: masuday/fsparse

max_forks_repo_name: masuday/fsparse


text:
      SUBROUTINE FB01VD( N, M, L, P, LDP, A, LDA, B, LDB, C, LDC, Q,
     $                   LDQ, R, LDR, K, LDK, TOL, IWORK, DWORK, LDWORK,
     $                   INFO )
C
C     SLICOT RELEASE 5.5.
C
C     Copyright (c) 2002-2012 NICONET e.V.
C
C     PURPOSE
C
C     To compute one recursion of the conventional Kalman filter
C     equations. This is one update of the Riccati difference equation
C     and the Kalman filter gain.
C
C     ARGUMENTS
C
C     Input/Output Parameters
C
C     N       (input) INTEGER
C             The actual state dimension, i.e., the order of the
C             matrices P      and A .  N >= 0.
C                       i|i-1      i
C
C     M       (input) INTEGER
C             The actual input dimension, i.e., the order of the matrix
C             Q .  M >= 0.
C              i
C
C     L       (input) INTEGER
C             The actual output dimension, i.e., the order of the matrix
C             R .  L >= 0.
C              i
C
C     P       (input/output) DOUBLE PRECISION array, dimension (LDP,N)
C             On entry, the leading N-by-N part of this array must
C             contain P     , the state covariance matrix at instant
C                      i|i-1
C             (i-1). The upper triangular part only is needed.
C             On exit, if INFO = 0, the leading N-by-N part of this
C             array contains P     , the state covariance matrix at
C                             i+1|i
C             instant i. The strictly lower triangular part is not set.
C             Otherwise, the leading N-by-N part of this array contains
C             P     , its input value.
C              i|i-1
C
C     LDP     INTEGER
C             The leading dimension of array P.  LDP >= MAX(1,N).
C
C     A       (input) DOUBLE PRECISION array, dimension (LDA,N)
C             The leading N-by-N part of this array must contain A ,
C                                                                 i
C             the state transition matrix of the discrete system at
C             instant i.
C
C     LDA     INTEGER
C             The leading dimension of array A.  LDA >= MAX(1,N).
C
C     B       (input) DOUBLE PRECISION array, dimension (LDB,M)
C             The leading N-by-M part of this array must contain B ,
C                                                                 i
C             the input weight matrix of the discrete system at
C             instant i.
C
C     LDB     INTEGER
C             The leading dimension of array B.  LDB >= MAX(1,N).
C
C     C       (input) DOUBLE PRECISION array, dimension (LDC,N)
C             The leading L-by-N part of this array must contain C ,
C                                                                 i
C             the output weight matrix of the discrete system at
C             instant i.
C
C     LDC     INTEGER
C             The leading dimension of array C.  LDC >= MAX(1,L).
C
C     Q       (input) DOUBLE PRECISION array, dimension (LDQ,M)
C             The leading M-by-M part of this array must contain Q ,
C                                                                 i
C             the input (process) noise covariance matrix at instant i.
C             The diagonal elements of this array are modified by the
C             routine, but are restored on exit.
C
C     LDQ     INTEGER
C             The leading dimension of array Q.  LDQ >= MAX(1,M).
C
C     R       (input/output) DOUBLE PRECISION array, dimension (LDR,L)
C             On entry, the leading L-by-L part of this array must
C             contain R , the output (measurement) noise covariance
C                      i
C             matrix at instant i.
C             On exit, if INFO = 0, or INFO = L+1, the leading L-by-L
C                                                                  1/2
C             upper triangular part of this array contains (RINOV )   ,
C                                                                i
C             the square root (left Cholesky factor) of the covariance
C             matrix of the innovations at instant i.
C
C     LDR     INTEGER
C             The leading dimension of array R.  LDR >= MAX(1,L).
C
C     K       (output) DOUBLE PRECISION array, dimension (LDK,L)
C             If INFO = 0, the leading N-by-L part of this array
C             contains K , the Kalman filter gain matrix at instant i.
C                       i
C             If INFO > 0, the leading N-by-L part of this array
C             contains the matrix product P     C'.
C                                          i|i-1 i
C
C     LDK     INTEGER
C             The leading dimension of array K.  LDK >= MAX(1,N).
C
C     Tolerances
C
C     TOL     DOUBLE PRECISION
C             The tolerance to be used to test for near singularity of
C             the matrix RINOV . If the user sets TOL > 0, then the
C                             i
C             given value of TOL is used as a lower bound for the
C             reciprocal condition number of that matrix; a matrix whose
C             estimated condition number is less than 1/TOL is
C             considered to be nonsingular. If the user sets TOL <= 0,
C             then an implicitly computed, default tolerance, defined by
C             TOLDEF = L*L*EPS, is used instead, where EPS is the
C             machine precision (see LAPACK Library routine DLAMCH).
C
C     Workspace
C
C     IWORK   INTEGER array, dimension (L)
C
C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
C             On exit, if INFO = 0, or INFO = L+1, DWORK(1) returns an
C             estimate of the reciprocal of the condition number (in the
C             1-norm) of the matrix RINOV .
C                                        i
C
C     LDWORK  The length of the array DWORK.
C             LDWORK >= MAX(1,L*N+3*L,N*N,N*M).
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -k, the k-th argument had an illegal
C                   value;
C             = k:  if INFO = k, 1 <= k <= L, the leading minor of order
C                   k of the matrix RINOV  is not positive-definite, and
C                                        i
C                   its Cholesky factorization could not be completed;
C             = L+1: the matrix RINOV  is singular, i.e., the condition
C                                    i
C                   number estimate of RINOV  (in the 1-norm) exceeds
C                                           i
C                   1/TOL.
C
C     METHOD
C
C     The conventional Kalman filter gain used at the i-th recursion
C     step is of the form
C
C                            -1
C        K  = P     C'  RINOV  ,
C         i    i|i-1 i       i
C
C     where RINOV  = C P     C' + R , and the state covariance matrix
C                i    i i|i-1 i    i
C
C     P      is updated by the discrete-time difference Riccati equation
C      i|i-1
C
C        P      = A  (P      - K C P     ) A'  + B Q B'.
C         i+1|i    i   i|i-1    i i i|i-1   i     i i i
C
C     Using these two updates, the combined time and measurement update
C     of the state X      is given by
C                   i|i-1
C
C        X      = A X      + A K (Y  - C X     ),
C         i+1|i    i i|i-1    i i  i    i i|i-1
C
C     where Y  is the new observation at step i.
C            i
C
C     REFERENCES
C
C     [1] Anderson, B.D.O. and Moore, J.B.
C         Optimal Filtering,
C         Prentice Hall, Englewood Cliffs, New Jersey, 1979.
C
C     [2] Verhaegen, M.H.G. and Van Dooren, P.
C         Numerical Aspects of Different Kalman Filter Implementations.
C         IEEE Trans. Auto. Contr., AC-31, pp. 907-917, 1986.
C
C     NUMERICAL ASPECTS
C
C     The algorithm requires approximately
C
C             3   2
C      3/2 x N + N  x (3 x L + M/2)
C
C     operations.
C
C     CONTRIBUTORS
C
C     Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Feb. 1997.
C     Supersedes Release 2.0 routine FB01JD by M.H.G. Verhaegen,
C     M. Vanbegin, and P. Van Dooren.
C
C     REVISIONS
C
C     February 20, 1998, November 20, 2003, April 20, 2004.
C
C     KEYWORDS
C
C     Kalman filtering, optimal filtering, recursive estimation.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION  ZERO, ONE, TWO
      PARAMETER         ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
C     .. Scalar Arguments ..
      INTEGER           INFO, L, LDA, LDB, LDC, LDK, LDP, LDQ, LDR,
     $                  LDWORK, M, N
      DOUBLE PRECISION  TOL
C     .. Array Arguments ..
      INTEGER           IWORK(*)
      DOUBLE PRECISION  A(LDA,*), B(LDB,*), C(LDC,*), DWORK(*),
     $                  K(LDK,*), P(LDP,*), Q(LDQ,*), R(LDR,*)
C     .. Local Scalars ..
      INTEGER           J, JWORK, LDW, N1
      DOUBLE PRECISION  RCOND, RNORM, TOLDEF
C     .. External Functions ..
      DOUBLE PRECISION  DLAMCH, DLANSY
      EXTERNAL          DLAMCH, DLANSY
C     .. External Subroutines ..
      EXTERNAL          DAXPY, DCOPY, DGEMV, DLACPY, DLASET, DPOCON,
     $                  DPOTRF, DSCAL, DTRMM, DTRSM, MB01RD, XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC         DBLE, MAX
C     .. Executable Statements ..
C
C     Test the input scalar arguments.
C
      INFO = 0
      N1 = MAX( 1, N )
      IF( N.LT.0 ) THEN
         INFO = -1
      ELSE IF( M.LT.0 ) THEN
         INFO = -2
      ELSE IF( L.LT.0 ) THEN
         INFO = -3
      ELSE IF( LDP.LT.N1 ) THEN
         INFO = -5
      ELSE IF( LDA.LT.N1 ) THEN
         INFO = -7
      ELSE IF( LDB.LT.N1 ) THEN
         INFO = -9
      ELSE IF( LDC.LT.MAX( 1, L ) ) THEN
         INFO = -11
      ELSE IF( LDQ.LT.MAX( 1, M ) ) THEN
         INFO = -13
      ELSE IF( LDR.LT.MAX( 1, L ) ) THEN
         INFO = -15
      ELSE IF( LDK.LT.N1 ) THEN
         INFO = -17
      ELSE IF( LDWORK.LT.MAX( 1, L*N + 3*L, N*N, N*M ) ) THEN
         INFO = -21
      END IF
C
      IF ( INFO.NE.0 ) THEN
C
C        Error return.
C
         CALL XERBLA( 'FB01VD', -INFO )
         RETURN
      END IF
C
C     Quick return if possible.
C
      IF ( MAX( N, L ).EQ.0 ) THEN
         DWORK(1) = ONE
         RETURN
      END IF
C
C     Efficiently compute RINOV = CPC' + R in R and put CP in DWORK and
C     PC' in K. (The content of DWORK on exit from MB01RD is used.)
C     Workspace: need L*N.
C
C     (Note: Comments in the code beginning "Workspace:" describe the
C     minimal amount of real workspace needed at that point in the
C     code.)
C
      CALL MB01RD( 'Upper', 'No transpose', L, N, ONE, ONE, R, LDR, C,
     $             LDC, P, LDP, DWORK, LDWORK, INFO )
      LDW = MAX( 1, L )
C
      DO 10 J = 1, L
         CALL DCOPY( N, DWORK(J), LDW, K(1,J), 1 )
   10 CONTINUE
C
      CALL DLACPY( 'Full', L, N, C, LDC, DWORK, LDW )
      CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Non-unit', L, N, ONE,
     $            P, LDP, DWORK, LDW )
      CALL DSCAL( N, TWO, P, LDP+1 )
C
      DO 20 J = 1, L
         CALL DAXPY( N, ONE, K(1,J), 1, DWORK(J), LDW )
         CALL DCOPY( N, DWORK(J), LDW, K(1,J), 1 )
   20 CONTINUE
C
C     Calculate the Cholesky decomposition U'U of the innovation
C     covariance matrix RINOV, and its reciprocal condition number.
C     Workspace: need L*N + 3*L.
C
      JWORK = L*N + 1
      RNORM = DLANSY( '1-norm', 'Upper', L, R, LDR, DWORK(JWORK) )
C
      TOLDEF = TOL
      IF ( TOLDEF.LE.ZERO )
     $   TOLDEF = DBLE( L*L )*DLAMCH( 'Epsilon' )
      CALL DPOTRF( 'Upper', L, R, LDR, INFO )
      IF ( INFO.NE.0 )
     $   RETURN
C
      CALL DPOCON( 'Upper', L, R, LDR, RNORM, RCOND, DWORK(JWORK),
     $             IWORK, INFO )
C
      IF ( RCOND.LT.TOLDEF ) THEN
C
C        Error return: RINOV is numerically singular.
C
         INFO = L+1
         DWORK(1) = RCOND
         RETURN
      END IF
C
      IF ( L.GT.1 )
     $   CALL DLASET( 'Lower', L-1, L-1, ZERO, ZERO, R(2,1),LDR )
C                                                          -1
C     Calculate the Kalman filter gain matrix  K = PC'RINOV .
C     Workspace: need L*N.
C
      CALL DTRSM( 'Right', 'Upper', 'No transpose', 'Non-unit', N, L,
     $            ONE, R, LDR, K, LDK )
      CALL DTRSM( 'Right', 'Upper', 'Transpose', 'Non-unit', N, L,
     $            ONE, R, LDR, K, LDK )
C
C     First part of the Riccati equation update: compute A(P-KCP)A'.
C     The upper triangular part of the symmetric matrix P-KCP is formed.
C     Workspace: need max(L*N,N*N).
C
      JWORK = 1
C
      DO 30 J = 1, N
         CALL DGEMV( 'No transpose', J, L, -ONE, K, LDK, DWORK(JWORK),
     $               1, ONE, P(1,J), 1 )
         JWORK = JWORK + L
   30 CONTINUE
C
      CALL MB01RD( 'Upper', 'No transpose', N, N, ZERO, ONE, P, LDP, A,
     $             LDA, P, LDP, DWORK, LDWORK, INFO )
C
C     Second part of the Riccati equation update: add BQB'.
C     Workspace: need N*M.
C
      CALL MB01RD( 'Upper', 'No transpose', N, M, ONE, ONE, P, LDP, B,
     $             LDB, Q, LDQ, DWORK, LDWORK, INFO )
      CALL DSCAL( M, TWO, Q, LDQ+1 )
C
C     Set the reciprocal of the condition number estimate.
C
      DWORK(1) = RCOND
C
      RETURN
C *** Last line of FB01VD ***
      END

ext: f

max_stars_repo_name: bgin/MissileSimulation

max_issues_repo_name: bgin/MissileSimulation

max_forks_repo_name: bgin/MissileSimulation


text:
module openmp
use types, only: dp
implicit none
private
public omp_get_thread_num, omp_get_num_threads, omp_get_max_threads, &
    omp_get_wtime, with_openmp

contains

    integer function omp_get_thread_num()
    omp_get_thread_num = 0
    end function

    integer function omp_get_num_threads()
    omp_get_num_threads = 1
    end function

    integer function omp_get_max_threads()
    omp_get_max_threads = 1
    end function

    real(dp) function omp_get_wtime()
    omp_get_wtime = 0
    end function

    logical function with_openmp()
    with_openmp = .false.
    end function

end module

max_stars_repo_name: certik/hfsolver

max_issues_repo_name: certik/hfsolver

max_forks_repo_name: certik/hfsolver


text:
subroutine amrex_probinit (init,name,namlen,problo,probhi) bind(c)

  use parallel
  use probdata_module
  use model_parser_module
  use bl_error_module
  use bl_constants_module
  use amrex_fort_module, only : rt => amrex_real
  use meth_params_module, only : const_grav

  implicit none

  integer init, namlen
  integer name(namlen)
  real(rt) :: problo(3), probhi(3)

  real(rt) :: offset
  integer untin,i

  namelist /fortin/ model_name, apply_vel_field, &
       velpert_scale, velpert_amplitude, velpert_height_loc, num_vortices, &
       shear_height_loc, shear_amplitude, &
       cutoff_density, interp_BC, zero_vels

  integer, parameter :: maxlen = 256
  character probin*(maxlen)

  real(rt) :: max_hse_err, dpdr, rhog, hse_err, dr_model

  ! Build "probin" filename from C++ land --
  ! the name of file containing fortin namelist.

  if (namlen .gt. maxlen) call bl_error("probin file name too long")

  do i = 1, namlen
     probin(i:i) = char(name(i))
  end do

  ! Namelist defaults
  apply_vel_field = .false.
  velpert_scale = 1.0e2_rt
  velpert_amplitude = 1.0e2_rt
  velpert_height_loc = 6.5e3_rt
  num_vortices = 1
  cutoff_density = 50.e0_rt
  interp_BC = .false.
  zero_vels = .false.
  shear_height_loc = 0.0d0

  ! Read namelists
  untin = 9
  open(untin,file=probin(1:namlen),form='formatted',status='old')
  read(untin,fortin)
  close(unit=untin)

  ! Read initial model
  call read_model_file(model_name)

  ! HSE check
  max_hse_err = -1.e30
  dr_model = model_r(2) - model_r(1)
  do i = 2, npts_model-2
     dpdr = (model_state(i,ipres_model) - model_state(i-1,ipres_model))/dr_model
     rhog = HALF*(model_state(i,idens_model) + model_state(i-1,idens_model))*const_grav
     hse_err = abs(dpdr - rhog)

     ! only count the error if we have a non-zero gradient in the next cell
     if (model_state(i+1,ipres_model) /= model_state(i,ipres_model)) then
        max_hse_err = max(max_hse_err, hse_err)
     endif
  enddo

  print *, 'maximum hse error = ', max_hse_err

  if (parallel_IOProcessor()) then
     do i = 1, npts_model
        print *, i, model_r(i), model_state(i,idens_model)
     enddo
  endif

  ! set local variable defaults
  center(1) = HALF*(problo(1)+probhi(1))
  center(2) = HALF*(problo(2)+probhi(2))
  center(3) = HALF*(problo(3)+probhi(3))

  ! velocity perturbation stuff
  offset = (probhi(1) - problo(1)) / (num_vortices)

  allocate(xloc_vortices(num_vortices))

  do i = 1, num_vortices
     xloc_vortices(i) = (dble(i-1) + HALF) * offset + problo(1)
  enddo

end subroutine amrex_probinit

! ::: -----------------------------------------------------------
! ::: This routine is called at problem setup time and is used
! ::: to initialize data on each grid.
! :::
! ::: NOTE:  all arrays have one cell of ghost zones surrounding
! :::        the grid interior.  Values in these cells need not
! :::        be set here.
! :::
! ::: INPUTS/OUTPUTS:
! :::
! ::: level     => amr level of grid
! ::: time      => time at which to init data
! ::: lo,hi     => index limits of grid interior (cell centered)
! ::: nstate    => number of state components.  You should know
! :::		   this already!
! ::: state     <=  Scalar array
! ::: delta     => cell size
! ::: xlo,xhi   => physical locations of lower left and upper
! :::              right hand corner of grid.  (does not include
! :::		   ghost region).
! ::: -----------------------------------------------------------
subroutine ca_initdata(level,time,lo,hi,nscal, &
     state,state_l1,state_l2,state_l3,state_h1,state_h2, &
     state_h3,delta,xlo,xhi)

  use bl_constants_module
  use probdata_module
  use interpolate_module
  use eos_module
  use meth_params_module, only : NVAR, URHO, UMX, UMY,UMZ, UEDEN, UEINT, UFS, UTEMP
  use network, only: nspec
  use model_parser_module
  use eos_type_module
  use amrex_fort_module, only : rt => amrex_real

  implicit none

  integer level, nscal
  integer lo(3), hi(3)
  integer state_l1,state_l2,state_l3,state_h1,state_h2,state_h3
  real(rt) :: xlo(3), xhi(3), time, delta(3)
  real(rt) :: state(state_l1:state_h1,state_l2:state_h2,state_l3:state_h3,NVAR)

  real(rt) :: xdist,ydist,zdist,x,y,z,r,upert(3)
  integer i,j,k,n,vortex

  type (eos_t) :: eos_state

  do k = lo(3), hi(3)
     z = xlo(3) + delta(3)*(dble(k-lo(3)) + HALF)
     do j = lo(2), hi(2)
        do i = lo(1), hi(1)

           state(i,j,k,URHO)  = interpolate(z,npts_model,model_r, &
                model_state(:,idens_model))
           state(i,j,k,UTEMP) = interpolate(z,npts_model,model_r, &
                model_state(:,itemp_model))

           do n = 1, nspec
              state(i,j,k,UFS-1+n) = interpolate(z,npts_model,model_r, &
                   model_state(:,ispec_model-1+n))
           enddo

           eos_state%rho = state(i,j,k,URHO)
           eos_state%T = state(i,j,k,UTEMP)
           eos_state%xn(:) = state(i,j,k,UFS:)

           call eos(eos_input_rt, eos_state)
           state(i,j,k,UEINT) = eos_state%e

        end do
     end do
  end do

  ! switch to conserved quantities
  do k = lo(3), hi(3)
     do j = lo(2), hi(2)
        do i = lo(1), hi(1)

           state(i,j,k,UEDEN) = state(i,j,k,URHO) * state(i,j,k,UEINT)
           state(i,j,k,UEINT) = state(i,j,k,URHO) * state(i,j,k,UEINT)

           do n = 1,nspec
              state(i,j,k,UFS+n-1) = state(i,j,k,URHO) * state(i,j,k,UFS+n-1)
           end do

        enddo
     end do
  enddo

  ! Initial velocities
  state(:,:,:,UMX:UMZ) = ZERO

  ! Now add the velocity perturbation (update the kinetic energy too)
  if (apply_vel_field) then
     do k = lo(3), hi(3)
        z = xlo(3) + delta(3)*(dble(j-lo(3)) + HALF)
        zdist = z - velpert_height_loc
        do i = lo(1), hi(1)
           x = xlo(1) + delta(1)*(dble(i-lo(1)) + HALF)

           if (z >= shear_height_loc) then
              state(:,:,:,UMX) = state(:,:,:,URHO)*shear_amplitude
              state(:,:,:,UMY) = ZERO
              state(:,:,:,UMZ) = ZERO
           endif

           upert = ZERO

           ! loop over each vortex
           do vortex = 1, num_vortices

              xdist = x - xloc_vortices(vortex)
              r = sqrt(xdist**2 + zdist**2)

              upert(1) = upert(1) - (zdist/velpert_scale) * &
                   velpert_amplitude * exp( -r**2/(TWO*velpert_scale**2)) &
                   * (-ONE)**vortex

              upert(3) = upert(3) + (xdist/velpert_scale) * &
                   velpert_amplitude * exp(-r**2/(TWO*velpert_scale**2)) &
                   * (-ONE)**vortex
           enddo
           do j = lo(1), hi(1)
              if (mod(j, int(abs(xhi(2)-xlo(2)/delta(2))/4.0_rt)) == 0) then
                 state(i,j,k,UMX) = state(i,j,k,UMX) + state(i,j,k,URHO) * upert(1)
                 state(i,j,k,UMZ) = state(i,j,k,UMZ) + state(i,j,k,URHO) * upert(3)
                 state(i,j,k,UEDEN) = state(i,j,k,UEDEN) + HALF*sum(state(i,j,k,UMX:UMZ)**2)/state(i,j,k,URHO)
              end if
           end do

        end do
     end do
  endif

end subroutine ca_initdata

subroutine ca_initrad(level,time,lo,hi,nrad, &
     rad_state,rad_state_l1,rad_state_l2, rad_state_l3, &
     rad_state_h1,rad_state_h2, rad_state_h3, &
     delta,xlo,xhi)

  use bl_constants_module
  use probdata_module
  use amrex_fort_module, only : rt => amrex_real

  integer level, nrad
  integer lo(3), hi(3)
  integer rad_state_l1,rad_state_l2,rad_state_l3
  integer rad_state_h1,rad_state_h2,rad_state_h3
  real(rt) :: xlo(3), xhi(3), time, delta(3)
  real(rt) :: rad_state(rad_state_l1:rad_state_h1, &
       rad_state_l2:rad_state_h2,rad_state_l3:rad_state_h3, nrad)

  integer i,j,k
  do k = lo(3), hi(3)
     do j = lo(2), hi(2)
        do i = lo(1), hi(1)
           rad_state(i,j,k,:) = ZERO
        end do
     end do
  end do
end subroutine ca_initrad

max_stars_repo_name: yingtchen/Castro

max_issues_repo_name: yingtchen/Castro

max_forks_repo_name: yingtchen/Castro


text:
c     Subroutine to output variables: relh/wspd/airtmax/airtmin
c     AJ_Kettle, Oct3/2017

      SUBROUTINE export_datafile2(s_directory_output,s_date,l_zipfile,
     +  s_arch_geog_stnid,
     +  s_arch_geog_hgt_m,s_arch_geog_lat_deg,s_arch_geog_lon_deg,
     +  s_arch_geog_stdate,s_arch_geog_endate,i_arch_lprod,
     +  i_arch_dayavg_relh_ngood,f_arch_dayavg_relh_min_pc,
     +  f_arch_dayavg_relh_max_pc,f_arch_dayavg_relh_avg_pc,
     +  i_arch_dayavg_wspd_ngood,f_arch_dayavg_wspd_min_ms,
     +  f_arch_dayavg_wspd_max_ms,f_arch_dayavg_wspd_avg_ms,
     +  i_arch_daymax_airt_ngood,f_arch_daymax_airt_min_c,
     +  f_arch_daymax_airt_max_c,f_arch_daymax_airt_avg_c,
     +  i_arch_daymin_airt_ngood,f_arch_daymin_airt_min_c,
     +  f_arch_daymin_airt_max_c,f_arch_daymin_airt_avg_c)

      IMPLICIT NONE
c************************************************************************
      CHARACTER(LEN=200)  :: s_directory_output
      CHARACTER(LEN=8)    :: s_date
      INTEGER             :: l_zipfile

      INTEGER             :: i_arch_lprod(1100)
      CHARACTER(LEN=5)    :: s_arch_geog_stnid(1100)
      CHARACTER(LEN=8)    :: s_arch_geog_hgt_m(1100)
      CHARACTER(LEN=8)    :: s_arch_geog_lat_deg(1100)
      CHARACTER(LEN=8)    :: s_arch_geog_lon_deg(1100)
      CHARACTER(LEN=8)    :: s_arch_geog_stdate(1100)
      CHARACTER(LEN=8)    :: s_arch_geog_endate(1100)

      INTEGER             :: i_arch_dayavg_relh_ngood(1100)
      REAL                :: f_arch_dayavg_relh_ngood(1100)
      REAL                :: f_arch_dayavg_relh_nbad(1100)
      REAL                :: f_arch_dayavg_relh_min_pc(1100)
      REAL                :: f_arch_dayavg_relh_max_pc(1100)
      REAL                :: f_arch_dayavg_relh_avg_pc(1100)

      INTEGER             :: i_arch_dayavg_wspd_ngood(1100)
      REAL                :: f_arch_dayavg_wspd_ngood(1100)
      REAL                :: f_arch_dayavg_wspd_nbad(1100)
      REAL                :: f_arch_dayavg_wspd_min_ms(1100)
      REAL                :: f_arch_dayavg_wspd_max_ms(1100)
      REAL                :: f_arch_dayavg_wspd_avg_ms(1100)

      INTEGER             :: i_arch_daymax_airt_ngood(1100)
      REAL                :: f_arch_daymax_airt_ngood(1100)
      REAL                :: f_arch_daymax_airt_nbad(1100)
      REAL                :: f_arch_daymax_airt_min_c(1100)
      REAL                :: f_arch_daymax_airt_max_c(1100)
      REAL                :: f_arch_daymax_airt_avg_c(1100)

      INTEGER             :: i_arch_daymin_airt_ngood(1100)
      REAL                :: f_arch_daymin_airt_ngood(1100)
      REAL                :: f_arch_daymin_airt_nbad(1100)
      REAL                :: f_arch_daymin_airt_min_c(1100)
      REAL                :: f_arch_daymin_airt_max_c(1100)
      REAL                :: f_arch_daymin_airt_avg_c(1100)

      INTEGER             :: i,j,k
c************************************************************************

      OPEN(UNIT=1,
     +  FILE=TRIM(s_directory_output)//'stats_day_data2.dat',
     +  FORM='formatted',
     +  STATUS='NEW',ACTION='WRITE')

      WRITE(UNIT=1,FMT=3009) 'DWD - Day Historical Data'
      WRITE(UNIT=1,FMT=3009) 'relh/wspd/airtmax/airtmin'
      WRITE(UNIT=1,FMT=3009) '                         '
      WRITE(UNIT=1,FMT=3009) '                         '
3009  FORMAT(a25)

      WRITE(unit=1,FMT=3029) 'AJ_Kettle ',s_date
3029  FORMAT(t1,a10,t12,a8)
      WRITE(unit=1,FMT=3023) '                                        '
3023  FORMAT(t1,a40)

      WRITE(unit=1,FMT=3025) 'Stnid',
     +  'Height  ','Latitude','Longitud',
     +  'St_date ','En_date ','N_TOTAL ',
     +  'RELH  ','RELH  ','RELH  ','RELH  ',
     +  'WSPD  ','WSPD  ','WSPD  ','WSPD  ',
     +  'AIRTH ','AIRTH ','AIRTH ','AIRTH ',
     +  'AIRTL ','AIRTL ','AIRTL ','AIRTL '
      WRITE(unit=1,FMT=3025) 'Stnid',
     +  'Height  ','Latitude','Longitud',
     +  'St_date ','En_date ','N_TOTAL ',
     +  'NGOOD ','MIN   ','MAX   ','AVG   ',
     +  'NGOOD ','MIN   ','MAX   ','AVG   ',
     +  'NGOOD ','MIN   ','MAX   ','AVG   ',
     +  'NGOOD ','MIN   ','MAX   ','AVG   '
      WRITE(unit=1,FMT=3025) 'Stnid',
     +  'Height  ','Latitude','Longitud',
     +  'St_date ','En_date ','N_TOTAL ',
     +  '      ','%     ','%     ','%     ',
     +  '      ','m/s   ','m/s   ','m/s   ',
     +  '      ','C     ','C     ','C     ',
     +  '      ','C     ','C     ','C     '
      WRITE(unit=1,FMT=3025) '+---+',
     +  '+------+','+------+','+------+',
     +  '+------+','+------+','+------+',
     +  '+----+','+----+','+----+','+----+',
     +  '+----+','+----+','+----+','+----+',
     +  '+----+','+----+','+----+','+----+',
     +  '+----+','+----+','+----+','+----+'
3025  FORMAT(t1,a5,t7,
     +  a8,t16,a8,t25,a8,t34,
     +  a8,t43,a8,t52,a8,t61,
     +  a6,t67,a6,t73,a6,t79,a6,t85,
     +  a6,t91,a6,t97,a6,t103,a6,t109,
     +  a6,t115,a6,t121,a6,t127,a6,t133,
     +  a6,t139,a6,t145,a6,t151,a6)

      DO i=1,l_zipfile

       WRITE(unit=1,FMT=3027)
     +  s_arch_geog_stnid(i),
     +  s_arch_geog_hgt_m(i),
     +    s_arch_geog_lat_deg(i),s_arch_geog_lon_deg(i),
     +  s_arch_geog_stdate(i),s_arch_geog_endate(i),i_arch_lprod(i),
     +  i_arch_dayavg_relh_ngood(i),f_arch_dayavg_relh_min_pc(i),
     +  f_arch_dayavg_relh_max_pc(i),f_arch_dayavg_relh_avg_pc(i),
     +  i_arch_dayavg_wspd_ngood(i),f_arch_dayavg_wspd_min_ms(i),
     +  f_arch_dayavg_wspd_max_ms(i),f_arch_dayavg_wspd_avg_ms(i),
     +  i_arch_daymax_airt_ngood(i),f_arch_daymax_airt_min_c(i),
     +  f_arch_daymax_airt_max_c(i),f_arch_daymax_airt_avg_c(i),
     +  i_arch_daymin_airt_ngood(i),f_arch_daymin_airt_min_c(i),
     +  f_arch_daymin_airt_max_c(i),f_arch_daymin_airt_avg_c(i)
3027   FORMAT(t1,a5,t7,
     +  a8,t16,a8,t25,a8,t34,
     +  a8,t43,a8,t52,i8,t61,
     +  i6,t67,f6.1,t73,f6.1,t79,f6.2,t85,
     +  i6,t91,f6.2,t97,f6.2,t103,f6.2,t109,
     +  i6,t115,f6.1,t121,f6.1,t127,f6.2,t133,
     +  i6,t139,f6.1,t145,f6.1,t151,f6.2)
      ENDDO

      CLOSE(UNIT=1)

      RETURN
      END

max_stars_repo_name: ajkettle/glamod-nuim

max_issues_repo_name: ajkettle/glamod-nuim

max_forks_repo_name: ajkettle/glamod-nuim


text:
      program master

      include 'p4f.h'

      integer*4 i,slaves,type,from,retcde,recvlen,buflen,loops,time
      character*500 buffer
      integer*4 TAGCNT, TAGDAT, TAGEND
      parameter (TAGCNT = 10)
      parameter (TAGDAT = 20)
      parameter (TAGEND = 30)

      call p4init()
      call p4crpg()

      slaves = p4ntotids() - 1
      length = 0
      buflen = 0

      do 10 i = 1,slaves
         call p4sendr(TAGCNT,i,buffer,length,retcde)
10    continue

20    print *,'How long is the message: '
      read (*,*,end=50) buflen
      buflen = MIN(500,buflen)

25    print *,'How many times around the loop? '
      read (*,*,end=50) loops

      length = buflen
      time = p4clock()
      do 45 j=1,loops

C send in reverse order of the slaves
      call p4send(TAGDAT,slaves,buffer,length,retcde)
      buffer = ' '
      type = TAGDAT
      from = -1
      call p4recv(type,from,buffer,buflen,recvlen,type,retcde)

      print *,'Loop number = ',j
      print 101,'MASTER receives from = ',from,' length = ',recvlen
101   format(a,i2,a,i6,/)
      length = recvlen
45    continue
      oldtime = time
      time = p4clock()
      print *,'Time used: ',time-oldtime, ' milliseconds'
      length = 0
50    continue

      do 60 i = 1,slaves
         call p4sendr(TAGEND,i,buffer,buflen,retcde)
60    continue

      call p4cleanup()
      print *,'Master exiting normally'
      end

max_stars_repo_name: RWTH-OS/MP-MPICH

max_issues_repo_name: RWTH-OS/MP-MPICH

max_forks_repo_name: RWTH-OS/MP-MPICH