SPAR.f

Go to the documentation of this file.

C                                                                1/15/81

C***********************************************************************

C  ODRV -- DRIVER FOR SPARSE MATRIX REORDERING ROUTINES

C***********************************************************************

        SUBROUTINE  odrvd

     *     (n, ia,ja,a, p,ip, nsp,isp, path, flag)

C

C  DESCRIPTION

C

C    ODRV FINDS A MINIMUM DEGREE ORDERING OF THE ROWS AND COLUMNS OF A

C    SYMMETRIC MATRIX M STORED IN (IA,JA,A) FORMAT (SEE BELOW).  FOR THE

C    REORDERED MATRIX, THE WORK AND STORAGE REQUIRED TO PERFORM GAUSSIAN

C    ELIMINATION IS (USUALLY) SIGNIFICANTLY LESS.

C

C    IF ONLY THE NONZERO ENTRIES IN THE UPPER TRIANGLE OF M ARE BEING

C    STORED, THEN ODRV SYMMETRICALLY REORDERS (IA,JA,A), (OPTIONALLY)

C    WITH THE DIAGONAL ENTRIES PLACED FIRST IN EACH ROW.  THIS IS TO

C    ENSURE THAT IF M(I,J) WILL BE IN THE UPPER TRIANGLE OF M WITH

C    RESPECT TO THE NEW ORDERING, THEN M(I,J) IS STORED IN ROW I (AND

C    THUS M(J,I) IS NOT STORED);  WHEREAS IF M(I,J) WILL BE IN THE

C    STRICT LOWER TRIANGLE OF M, THEN M(J,I) IS STORED IN ROW J (AND

C    THUS M(I,J) IS NOT STORED).

C

C

C  STORAGE OF SPARSE MATRICES

C

C    THE NONZERO ENTRIES OF THE MATRIX M ARE STORED ROW-BY-ROW IN THE

C    ARRAY A.  TO IDENTIFY THE INDIVIDUAL NONZERO ENTRIES IN EACH ROW,

C    WE NEED TO KNOW IN WHICH COLUMN EACH ENTRY LIES.  THESE COLUMN

C    INDICES ARE STORED IN THE ARRAY JA;  I.E., IF  A(K) = M(I,J),  THEN

C    JA(K) = J.  TO IDENTIFY THE INDIVIDUAL ROWS, WE NEED TO KNOW WHERE

C    EACH ROW STARTS.  THESE ROW POINTERS ARE STORED IN THE ARRAY IA;

C    I.E., IF M(I,J) IS THE FIRST NONZERO ENTRY (STORED) IN THE I-TH ROW

C    AND  A(K) = M(I,J),  THEN  IA(I) = K.  MOREOVER, IA(N+1) POINTS TO

C    THE FIRST LOCATION FOLLOWING THE LAST ELEMENT IN THE LAST ROW.

C    THUS, THE NUMBER OF ENTRIES IN THE I-TH ROW IS  IA(I+1) - IA(I),

C    THE NONZERO ENTRIES IN THE I-TH ROW ARE STORED CONSECUTIVELY IN

C

C            A(IA(I)),  A(IA(I)+1),  ..., A(IA(I+1)-1),

C

C    AND THE CORRESPONDING COLUMN INDICES ARE STORED CONSECUTIVELY IN

C

C            JA(IA(I)), JA(IA(I)+1), ..., JA(IA(I+1)-1).

C

C    SINCE THE COEFFICIENT MATRIX IS SYMMETRIC, ONLY THE NONZERO ENTRIES

C    IN THE UPPER TRIANGLE NEED BE STORED.  FOR EXAMPLE, THE MATRIX

C

C             ( 1  0  2  3  0 )

C             ( 0  4  0  0  0 )

C         M = ( 2  0  5  6  0 )

C             ( 3  0  6  7  8 )

C             ( 0  0  0  8  9 )

C

C    COULD BE STORED AS

C

C              1  2  3  4  5  6  7  8  9 10 11 12 13

C         ---+--------------------------------------

C         IA   1  4  5  8 12 14

C         JA   1  3  4  2  1  3  4  1  3  4  5  4  5

C          A   1  2  3  4  2  5  6  3  6  7  8  8  9

C

C    OR (SYMMETRICALLY) AS

C

C              1  2  3  4  5  6  7  8  9

C         ---+--------------------------

C         IA   1  4  5  7  9 10

C         JA   1  3  4  2  3  4  4  5  5

C          A   1  2  3  4  5  6  7  8  9          .

C

C

C  PARAMETERS

C

C    N    - ORDER OF THE MATRIX

C

C    IA   - INTEGER ONE-DIMENSIONAL ARRAY CONTAINING POINTERS TO DELIMIT

C           ROWS IN JA AND A;  DIMENSION = N+1

C

C    JA   - INTEGER ONE-DIMENSIONAL ARRAY CONTAINING THE COLUMN INDICES

C           CORRESPONDING TO THE ELEMENTS OF A;  DIMENSION = NUMBER OF

C           NONZERO ENTRIES IN (THE UPPER TRIANGLE OF) M

C

C    A    - REAL ONE-DIMENSIONAL ARRAY CONTAINING THE NONZERO ENTRIES IN

C           (THE UPPER TRIANGLE OF) M, STORED BY ROWS;  DIMENSION =

C           NUMBER OF NONZERO ENTRIES IN (THE UPPER TRIANGLE OF) M

C

C    P    - INTEGER ONE-DIMENSIONAL ARRAY USED TO RETURN THE PERMUTATION

C           OF THE ROWS AND COLUMNS OF M CORRESPONDING TO THE MINIMUM

C           DEGREE ORDERING;  DIMENSION = N

C

C    IP   - INTEGER ONE-DIMENSIONAL ARRAY USED TO RETURN THE INVERSE OF

C           THE PERMUTATION RETURNED IN P;  DIMENSION = N

C

C    NSP  - DECLARED DIMENSION OF THE ONE-DIMENSIONAL ARRAY ISP;  NSP

C           MUST BE AT LEAST  3N+4K,  WHERE K IS THE NUMBER OF NONZEROES

C           IN THE STRICT UPPER TRIANGLE OF M

C

C    ISP  - INTEGER ONE-DIMENSIONAL ARRAY USED FOR WORKING STORAGE;

C           DIMENSION = NSP

C

C    PATH - INTEGER PATH SPECIFICATION;  VALUES AND THEIR MEANINGS ARE -

C             1  FIND MINIMUM DEGREE ORDERING ONLY

C             2  FIND MINIMUM DEGREE ORDERING AND REORDER SYMMETRICALLY

C                  STORED MATRIX (USED WHEN ONLY THE NONZERO ENTRIES IN

C                  THE UPPER TRIANGLE OF M ARE BEING STORED)

C             3  REORDER SYMMETRICALLY STORED MATRIX AS SPECIFIED BY

C                  INPUT PERMUTATION (USED WHEN AN ORDERING HAS ALREADY

C                  BEEN DETERMINED AND ONLY THE NONZERO ENTRIES IN THE

C                  UPPER TRIANGLE OF M ARE BEING STORED)

C             4  SAME AS 2 BUT PUT DIAGONAL ENTRIES AT START OF EACH ROW

C             5  SAME AS 3 BUT PUT DIAGONAL ENTRIES AT START OF EACH ROW

C

C    FLAG - INTEGER ERROR FLAG;  VALUES AND THEIR MEANINGS ARE -

C               0    NO ERRORS DETECTED

C              9N+K  INSUFFICIENT STORAGE IN MD

C             10N+1  INSUFFICIENT STORAGE IN ODRV

C             11N+1  ILLEGAL PATH SPECIFICATION

C

C

C  CONVERSION FROM REAL TO DOUBLE PRECISION

C

C    CHANGE THE REAL DECLARATIONS IN ODRV AND SRO TO DOUBLE PRECISION

C    DECLARATIONS.

C

C-----------------------------------------------------------------------

C

        INTEGER  ia(n+1), ja(*),  p(n), ip(n),  isp(nsp),  path,  flag,

     *     v, l, head,  tmp, q

C...    REAL  A(1)

        DOUBLE PRECISION  a(*)

        LOGICAL  dflag

C

C----INITIALIZE ERROR FLAG AND VALIDATE PATH SPECIFICATION

        flag = 0

        IF (path.LT.1 .OR. 5.LT.path)  go to 111

C

C----ALLOCATE STORAGE AND FIND MINIMUM DEGREE ORDERING

        IF ((path-1) * (path-2) * (path-4) .NE. 0)  go to 1

          max = (nsp-n)/2

          v    = 1

          l    = v     +  max

          head = l     +  max

          next = head  +  n

          IF (max.LT.n)  go to 110

C

          CALL  md

     *       (n, ia,ja, max,isp(v),isp(l), isp(head),p,ip, isp(v), flag)

          IF (flag.NE.0)  go to 100

C

C----ALLOCATE STORAGE AND SYMMETRICALLY REORDER MATRIX

   1    IF ((path-2) * (path-3) * (path-4) * (path-5) .NE. 0)  go to 2

          tmp = (nsp+1) -      n

          q   = tmp     - (ia(n+1)-1)

          IF (q.LT.1)  go to 110

C

          dflag = path.EQ.4 .OR. path.EQ.5

          CALL sro

     *       (n,  ip,  ia, ja, a,  isp(tmp),  isp(q),  dflag)

C

   2    RETURN

C

C ** ERROR -- ERROR DETECTED IN MD

 100    RETURN

C ** ERROR -- INSUFFICIENT STORAGE

 110    flag = 10*n + 1

        RETURN

C ** ERROR -- ILLEGAL PATH SPECIFIED

 111    flag = 11*n + 1

        RETURN

        END

C***********************************************************************

C***********************************************************************

C  MD -- MINIMUM DEGREE ALGORITHM (BASED ON ELEMENT MODEL)

C***********************************************************************

        SUBROUTINE  md

     *     (n, ia,ja, max, v,l, head,last,next, mark, flag)

C

C  DESCRIPTION

C

C    MD FINDS A MINIMUM DEGREE ORDERING OF THE ROWS AND COLUMNS OF A

C    SYMMETRIC MATRIX M STORED IN (IA,JA,A) FORMAT.

C

C

C  ADDITIONAL PARAMETERS

C

C    MAX  - DECLARED DIMENSION OF THE ONE-DIMENSIONAL ARRAYS V AND L;

C           MAX MUST BE AT LEAST  N+2K,  WHERE K IS THE NUMBER OF

C           NONZEROES IN THE STRICT UPPER TRIANGLE OF M

C

C    V    - INTEGER ONE-DIMENSIONAL WORK ARRAY;  DIMENSION = MAX

C

C    L    - INTEGER ONE-DIMENSIONAL WORK ARRAY;  DIMENSION = MAX

C

C    HEAD - INTEGER ONE-DIMENSIONAL WORK ARRAY;  DIMENSION = N

C

C    LAST - INTEGER ONE-DIMENSIONAL ARRAY USED TO RETURN THE PERMUTATION

C           OF THE ROWS AND COLUMNS OF M CORRESPONDING TO THE MINIMUM

C           DEGREE ORDERING;  DIMENSION = N

C

C    NEXT - INTEGER ONE-DIMENSIONAL ARRAY USED TO RETURN THE INVERSE OF

C           THE PERMUTATION RETURNED IN LAST;  DIMENSION = N

C

C    MARK - INTEGER ONE-DIMENSIONAL WORK ARRAY (MAY BE THE SAME AS V);

C           DIMENSION = N

C

C    FLAG - INTEGER ERROR FLAG;  VALUES AND THEIR MEANINGS ARE -

C             0      NO ERRORS DETECTED

C             11N+1  INSUFFICIENT STORAGE IN MD

C

C

C  DEFINITIONS OF INTERNAL PARAMETERS

C

C    ---------+---------------------------------------------------------

C    V(S)       VALUE FIELD OF LIST ENTRY

C    ---------+---------------------------------------------------------

C    L(S)       LINK FIELD OF LIST ENTRY  (0 => END OF LIST)

C    ---------+---------------------------------------------------------

C    L(VI)      POINTER TO ELEMENT LIST OF UNELIMINATED VERTEX VI

C    ---------+---------------------------------------------------------

C    L(EJ)      POINTER TO BOUNDARY LIST OF ACTIVE ELEMENT EJ

C    ---------+---------------------------------------------------------

C    HEAD(D)    VJ => VJ HEAD OF D-LIST D

C                0 => NO VERTEX IN D-LIST D

C

C

C                                VI UNELIMINATED VERTEX

C                        VI IN EK                   VI NOT IN EK

C    ---------+-----------------------------+---------------------------

C    NEXT(VI)   UNDEFINED BUT NONNEGATIVE     VJ => VJ NEXT IN D-LIST

C                                              0 => VI TAIL OF D-LIST

C    ---------+-----------------------------+---------------------------

C    LAST(VI)   (NOT SET UNTIL MDP)           -D => VI HEAD OF D-LIST D

C              -VK => COMPUTE DEGREE          VJ => VJ LAST IN D-LIST

C               EJ => VI PROTOTYPE OF EJ       0 => VI NOT IN ANY D-LIST

C                0 => DO NOT COMPUTE DEGREE

C    ---------+-----------------------------+---------------------------

C    MARK(VI)   MARK(VK)                      NONNEGATIVE TAG < MARK(VK)

C

C

C                                 VI ELIMINATED VERTEX

C                    EI ACTIVE ELEMENT                  OTHERWISE

C    ---------+-----------------------------+---------------------------

C    NEXT(VI)   -J => VI WAS J-TH VERTEX      -J => VI WAS J-TH VERTEX

C                     TO BE ELIMINATED              TO BE ELIMINATED

C    ---------+-----------------------------+---------------------------

C    LAST(VI)    M => SIZE OF EI = M          UNDEFINED

C    ---------+-----------------------------+---------------------------

C    MARK(VI)   -M => OVERLAP COUNT OF EI     UNDEFINED

C                     WITH EK = M

C               OTHERWISE NONNEGATIVE TAG

C                     < MARK(VK)

C

C-----------------------------------------------------------------------

C

        INTEGER  ia(*), ja(*),  v(*), l(*),  head(*), last(*), next(*),

     *     mark(*),  flag,  tag, dmin, vk,ek, tail

C

C----INITIALIZATION

        tag = 0

        CALL  mdi

     *     (n, ia,ja, max,v,l, head,last,next, mark,tag, flag)

        IF (flag.NE.0)  RETURN

C

        k = 0

        dmin = 1

C

C----WHILE  K < N  DO

   1    IF (k.GE.n)  go to 4

C

C------SEARCH FOR VERTEX OF MINIMUM DEGREE

   2      IF (head(dmin).GT.0)  go to 3

            dmin = dmin + 1

            go to 2

C

C------REMOVE VERTEX VK OF MINIMUM DEGREE FROM DEGREE LIST

   3      vk = head(dmin)

          head(dmin) = next(vk)

          IF (head(dmin).GT.0)  last(head(dmin)) = -dmin

C

C------NUMBER VERTEX VK, ADJUST TAG, AND TAG VK

          k = k+1

          next(vk) = -k

          last(vk) = dmin - 1

          tag = tag + last(vk)

          mark(vk) = tag

C

C------FORM ELEMENT EK FROM UNELIMINATED NEIGHBORS OF VK

          CALL  mdm

     *       (vk,tail, v,l, last,next, mark)

C

C------PURGE INACTIVE ELEMENTS AND DO MASS ELIMINATION

          CALL  mdp

     *       (k,vk,tail, v,l, head,last,next, mark)

C

C------UPDATE DEGREES OF UNELIMINATED VERTICES IN EK

          CALL  mdu

     *       (vk,dmin, v,l, head,last,next, mark)

C

          go to 1

C

C----GENERATE INVERSE PERMUTATION FROM PERMUTATION

   4    DO 5 k=1,n

          next(k) = -next(k)

   5      last(next(k)) = k

C

        RETURN

        END

C

C***********************************************************************

C  MDI -- INITIALIZATION

C***********************************************************************

        SUBROUTINE  mdi

     *     (n, ia,ja, max,v,l, head,last,next, mark,tag, flag)

        INTEGER  ia(*), ja(*),  v(*), l(*),  head(*), last(*), next(*),

     *     mark(*), tag,  flag,  sfs, vi,dvi, vj

C

C----INITIALIZE DEGREES, ELEMENT LISTS, AND DEGREE LISTS

        DO 1 vi=1,n

          mark(vi) = 1

          l(vi) = 0

   1      head(vi) = 0

        sfs = n+1

C

C----CREATE NONZERO STRUCTURE

C----FOR EACH NONZERO ENTRY A(VI,VJ) IN STRICT UPPER TRIANGLE

        DO 3 vi=1,n

          jmin = ia(vi)

          jmax = ia(vi+1) - 1

          IF (jmin.GT.jmax)  go to 3

          DO 2 j=jmin,jmax

            vj = ja(j)

            IF (vi.GE.vj)  go to 2

              IF (sfs.GE.max)  go to 101

C

C------ENTER VJ IN ELEMENT LIST FOR VI

              mark(vi) = mark(vi) + 1

              v(sfs) = vj

              l(sfs) = l(vi)

              l(vi) = sfs

              sfs = sfs+1

C

C------ENTER VI IN ELEMENT LIST FOR VJ

              mark(vj) = mark(vj) + 1

              v(sfs) = vi

              l(sfs) = l(vj)

              l(vj) = sfs

              sfs = sfs+1

   2        CONTINUE

   3      CONTINUE

C

C----CREATE DEGREE LISTS AND INITIALIZE MARK VECTOR

        DO 4 vi=1,n

          dvi = mark(vi)

          next(vi) = head(dvi)

          head(dvi) = vi

          last(vi) = -dvi

          IF (next(vi).GT.0)  last(next(vi)) = vi

   4      mark(vi) = tag

C

        RETURN

C

C ** ERROR -- INSUFFICIENT STORAGE

 101    flag = 9*n + vi

        RETURN

        END

C

C***********************************************************************

C  MDM -- FORM ELEMENT FROM UNELIMINATED NEIGHBORS OF VK

C***********************************************************************

        SUBROUTINE  mdm

     *     (vk,tail, v,l, last,next, mark)

        INTEGER  vk, tail,  v(*), l(*),   last(*), next(*),   mark(*),

     *     tag, s,ls,vs,es, b,lb,vb, blp,blpmax

C

C----INITIALIZE TAG AND LIST OF UNELIMINATED NEIGHBORS

        tag = mark(vk)

        tail = vk

C

C----FOR EACH VERTEX/ELEMENT VS/ES IN ELEMENT LIST OF VK

        ls = l(vk)

   1    s = ls

        IF (s.EQ.0)  go to 5

          ls = l(s)

          vs = v(s)

          IF (next(vs).LT.0)  go to 2

C

C------IF VS IS UNELIMINATED VERTEX, THEN TAG AND APPEND TO LIST OF

C------UNELIMINATED NEIGHBORS

            mark(vs) = tag

            l(tail) = s

            tail = s

            go to 4

C

C------IF ES IS ACTIVE ELEMENT, THEN ...

C--------FOR EACH VERTEX VB IN BOUNDARY LIST OF ELEMENT ES

   2        lb = l(vs)

            blpmax = last(vs)

            DO 3 blp=1,blpmax

              b = lb

              lb = l(b)

              vb = v(b)

C

C----------IF VB IS UNTAGGED VERTEX, THEN TAG AND APPEND TO LIST OF

C----------UNELIMINATED NEIGHBORS

              IF (mark(vb).GE.tag)  go to 3

                mark(vb) = tag

                l(tail) = b

                tail = b

   3          CONTINUE

C

C--------MARK ES INACTIVE

            mark(vs) = tag

C

   4      go to 1

C

C----TERMINATE LIST OF UNELIMINATED NEIGHBORS

   5    l(tail) = 0

C

        RETURN

        END

C

C***********************************************************************

C  MDP -- PURGE INACTIVE ELEMENTS AND DO MASS ELIMINATION

C***********************************************************************

        SUBROUTINE  mdp

     *     (k,ek,tail, v,l, head,last,next, mark)

        INTEGER  ek, tail,  v(*), l(*),  head(*), last(*), next(*),

     *     mark(*),  tag, free, li,vi,lvi,evi, s,ls,vs, ilp,ilpmax

C

C----INITIALIZE TAG

        tag = mark(ek)

C

C----FOR EACH VERTEX VI IN EK

        li = ek

        ilpmax = last(ek)

        IF (ilpmax.LE.0)  go to 12

        DO 11 ilp=1,ilpmax

          i = li

          li = l(i)

          vi = v(li)

C

C------REMOVE VI FROM DEGREE LIST

          IF (last(vi).EQ.0)  go to 3

            IF (last(vi).GT.0)  go to 1

              head(-last(vi)) = next(vi)

              go to 2

   1          next(last(vi)) = next(vi)

   2        IF (next(vi).GT.0)  last(next(vi)) = last(vi)

C

C------REMOVE INACTIVE ITEMS FROM ELEMENT LIST OF VI

   3      ls = vi

   4      s = ls

          ls = l(s)

          IF (ls.EQ.0)  go to 6

            vs = v(ls)

            IF (mark(vs).LT.tag)  go to 5

              free = ls

              l(s) = l(ls)

              ls = s

   5        go to 4

C

C------IF VI IS INTERIOR VERTEX, THEN REMOVE FROM LIST AND ELIMINATE

   6      lvi = l(vi)

          IF (lvi.NE.0)  go to 7

            l(i) = l(li)

            li = i

C

            k = k+1

            next(vi) = -k

            last(ek) = last(ek) - 1

            go to 11

C

C------ELSE ...

C--------CLASSIFY VERTEX VI

   7        IF (l(lvi).NE.0)  go to 9

              evi = v(lvi)

              IF (next(evi).GE.0)  go to 9

                IF (mark(evi).LT.0)  go to 8

C

C----------IF VI IS PROTOTYPE VERTEX, THEN MARK AS SUCH, INITIALIZE

C----------OVERLAP COUNT FOR CORRESPONDING ELEMENT, AND MOVE VI TO END

C----------OF BOUNDARY LIST

                  last(vi) = evi

                  mark(evi) = -1

                  l(tail) = li

                  tail = li

                  l(i) = l(li)

                  li = i

                  go to 10

C

C----------ELSE IF VI IS DUPLICATE VERTEX, THEN MARK AS SUCH AND ADJUST

C----------OVERLAP COUNT FOR CORRESPONDING ELEMENT

   8              last(vi) = 0

                  mark(evi) = mark(evi) - 1

                  go to 10

C

C----------ELSE MARK VI TO COMPUTE DEGREE

   9              last(vi) = -ek

C

C--------INSERT EK IN ELEMENT LIST OF VI

  10        v(free) = ek

            l(free) = l(vi)

            l(vi) = free

  11      CONTINUE

C

C----TERMINATE BOUNDARY LIST

  12    l(tail) = 0

C

        RETURN

        END

C

C***********************************************************************

C  MDU -- UPDATE DEGREES OF UNELIMINATED VERTICES IN EK

C***********************************************************************

        SUBROUTINE  mdu

     *     (ek,dmin, v,l, head,last,next, mark)

        INTEGER  ek, dmin,  v(*), l(*),  head(*), last(*), next(*),

     *     mark(*),  tag, vi,evi,dvi, s,vs,es, b,vb, ilp,ilpmax,

     *     blp,blpmax

C

C----INITIALIZE TAG

        tag = mark(ek) - last(ek)

C

C----FOR EACH VERTEX VI IN EK

        i = ek

        ilpmax = last(ek)

        IF (ilpmax.LE.0)  go to 11

        DO 10 ilp=1,ilpmax

          i = l(i)

          vi = v(i)

          IF (last(vi))  1, 10, 8

C

C------IF VI NEITHER PROTOTYPE NOR DUPLICATE VERTEX, THEN MERGE ELEMENTS

C------TO COMPUTE DEGREE

   1        tag = tag + 1

            dvi = last(ek)

C

C--------FOR EACH VERTEX/ELEMENT VS/ES IN ELEMENT LIST OF VI

            s = l(vi)

   2        s = l(s)

            IF (s.EQ.0)  go to 9

              vs = v(s)

              IF (next(vs).LT.0)  go to 3

C

C----------IF VS IS UNELIMINATED VERTEX, THEN TAG AND ADJUST DEGREE

                mark(vs) = tag

                dvi = dvi + 1

                go to 5

C

C----------IF ES IS ACTIVE ELEMENT, THEN EXPAND

C------------CHECK FOR OUTMATCHED VERTEX

   3            IF (mark(vs).LT.0)  go to 6

C

C------------FOR EACH VERTEX VB IN ES

                b = vs

                blpmax = last(vs)

                DO 4 blp=1,blpmax

                  b = l(b)

                  vb = v(b)

C

C--------------IF VB IS UNTAGGED, THEN TAG AND ADJUST DEGREE

                  IF (mark(vb).GE.tag)  go to 4

                    mark(vb) = tag

                    dvi = dvi + 1

   4              CONTINUE

C

   5          go to 2

C

C------ELSE IF VI IS OUTMATCHED VERTEX, THEN ADJUST OVERLAPS BUT DO NOT

C------COMPUTE DEGREE

   6        last(vi) = 0

            mark(vs) = mark(vs) - 1

   7        s = l(s)

            IF (s.EQ.0)  go to 10

              vs = v(s)

              IF (mark(vs).LT.0)  mark(vs) = mark(vs) - 1

              go to 7

C

C------ELSE IF VI IS PROTOTYPE VERTEX, THEN CALCULATE DEGREE BY

C------INCLUSION/EXCLUSION AND RESET OVERLAP COUNT

   8        evi = last(vi)

            dvi = last(ek) + last(evi) + mark(evi)

            mark(evi) = 0

C

C------INSERT VI IN APPROPRIATE DEGREE LIST

   9      next(vi) = head(dvi)

          head(dvi) = vi

          last(vi) = -dvi

          IF (next(vi).GT.0)  last(next(vi)) = vi

          IF (dvi.LT.dmin)  dmin = dvi

C

  10      CONTINUE

C

  11    RETURN

        END

C***********************************************************************

C***********************************************************************

C  SRO -- SYMMETRIC REORDERING OF SPARSE SYMMETRIC MATRIX

C***********************************************************************

        SUBROUTINE  sro

     *     (n, ip, ia,ja,a, q, r, dflag)

C

C  DESCRIPTION

C

C    THE NONZERO ENTRIES OF THE MATRIX M ARE ASSUMED TO BE STORED

C    SYMMETRICALLY IN (IA,JA,A) FORMAT (I.E., NOT BOTH M(I,J) AND M(J,I)

C    ARE STORED IF I NE J).

C

C    SRO DOES NOT REARRANGE THE ORDER OF THE ROWS, BUT DOES MOVE

C    NONZEROES FROM ONE ROW TO ANOTHER TO ENSURE THAT IF M(I,J) WILL BE

C    IN THE UPPER TRIANGLE OF M WITH RESPECT TO THE NEW ORDERING, THEN

C    M(I,J) IS STORED IN ROW I (AND THUS M(J,I) IS NOT STORED);  WHEREAS

C    IF M(I,J) WILL BE IN THE STRICT LOWER TRIANGLE OF M, THEN M(J,I) IS

C    STORED IN ROW J (AND THUS M(I,J) IS NOT STORED).

C

C

C  ADDITIONAL PARAMETERS

C

C    Q     - INTEGER ONE-DIMENSIONAL WORK ARRAY;  DIMENSION = N

C

C    R     - INTEGER ONE-DIMENSIONAL WORK ARRAY;  DIMENSION = NUMBER OF

C            NONZERO ENTRIES IN THE UPPER TRIANGLE OF M

C

C    DFLAG - LOGICAL VARIABLE;  IF DFLAG = .TRUE., THEN STORE NONZERO

C            DIAGONAL ELEMENTS AT THE BEGINNING OF THE ROW

C

C-----------------------------------------------------------------------

C

        INTEGER  ip(*),  ia(*), ja(*),  q(*), r(*)

C...    REAL  A(1),  AK

        DOUBLE PRECISION  a(*),  ak

        LOGICAL  dflag

C

C

C--PHASE 1 -- FIND ROW IN WHICH TO STORE EACH NONZERO

C----INITIALIZE COUNT OF NONZEROES TO BE STORED IN EACH ROW

        DO 1 i=1,n

  1       q(i) = 0

C

C----FOR EACH NONZERO ELEMENT A(J)

        DO 3 i=1,n

          jmin = ia(i)

          jmax = ia(i+1) - 1

          IF (jmin.GT.jmax)  go to 3

          DO 2 j=jmin,jmax

C

C--------FIND ROW (=R(J)) AND COLUMN (=JA(J)) IN WHICH TO STORE A(J) ...

            k = ja(j)

            IF (ip(k).LT.ip(i))  ja(j) = i

            IF (ip(k).GE.ip(i))  k = i

            r(j) = k

C

C--------... AND INCREMENT COUNT OF NONZEROES (=Q(R(J)) IN THAT ROW

  2         q(k) = q(k) + 1

  3       CONTINUE

C

C

C--PHASE 2 -- FIND NEW IA AND PERMUTATION TO APPLY TO (JA,A)

C----DETERMINE POINTERS TO DELIMIT ROWS IN PERMUTED (JA,A)

        DO 4 i=1,n

          ia(i+1) = ia(i) + q(i)

  4       q(i) = ia(i+1)

C

C----DETERMINE WHERE EACH (JA(J),A(J)) IS STORED IN PERMUTED (JA,A)

C----FOR EACH NONZERO ELEMENT (IN REVERSE ORDER)

        ilast = 0

        jmin = ia(1)

        jmax = ia(n+1) - 1

        j = jmax

        DO 6 jdummy=jmin,jmax

          i = r(j)

          IF (.NOT.dflag .OR. ja(j).NE.i .OR. i.EQ.ilast)  go to 5

C

C------IF DFLAG, THEN PUT DIAGONAL NONZERO AT BEGINNING OF ROW

            r(j) = ia(i)

            ilast = i

            go to 6

C

C------PUT (OFF-DIAGONAL) NONZERO IN LAST UNUSED LOCATION IN ROW

  5         q(i) = q(i) - 1

            r(j) = q(i)

C

  6       j = j-1

C

C

C--PHASE 3 -- PERMUTE (JA,A) TO UPPER TRIANGULAR FORM (WRT NEW ORDERING)

        DO 8 j=jmin,jmax

  7       IF (r(j).EQ.j)  go to 8

            k = r(j)

            r(j) = r(k)

            r(k) = k

            jak = ja(k)

            ja(k) = ja(j)

            ja(j) = jak

            ak = a(k)

            a(k) = a(j)

            a(j) = ak

            go to 7

  8       CONTINUE

C

        RETURN

        END

C                                                                1/15/81

C***********************************************************************

C  SDRV -- DRIVER FOR SPARSE SYMMETRIC POSITIVE DEFINITE MATRIX ROUTINES

C***********************************************************************

        SUBROUTINE sdrvd

     *     (n, p,ip, ia,ja,a, b, z, nsp,isp,rsp,esp, path, flag)

C

C  DESCRIPTION

C

C    SDRV SOLVES SPARSE SYMMETRIC POSITIVE DEFINITE SYSTEMS OF LINEAR

C    EQUATIONS.  THE SOLUTION PROCESS IS DIVIDED INTO THREE STAGES --

C

C      SSF - THE COEFFICIENT MATRIX M IS FACTORED SYMBOLICALLY TO

C            DETERMINE WHERE FILLIN WILL OCCUR DURING THE NUMERIC

C            FACTORIZATION.

C

C      SNF - M IS FACTORED NUMERICALLY INTO THE PRODUCT UT-D-U, WHERE

C            D IS DIAGONAL AND U IS UNIT UPPER TRIANGULAR.

C

C      SNS - THE LINEAR SYSTEM  MX = B  IS SOLVED USING THE UT-D-U

C            FACTORIZATION FROM SNF.

C

C    FOR SEVERAL SYSTEMS WITH THE SAME COEFFICIENT MATRIX, SSF AND SNF

C    NEED BE DONE ONLY ONCE (FOR THE FIRST SYSTEM);  THEN SNS IS DONE

C    ONCE FOR EACH ADDITIONAL RIGHT-HAND SIDE.  FOR SEVERAL SYSTEMS

C    WHOSE COEFFICIENT MATRICES HAVE THE SAME NONZERO STRUCTURE, SSF

C    NEED BE DONE ONLY ONCE (FOR THE FIRST SYSTEM);  THEN SNF AND SNS

C    ARE DONE ONCE FOR EACH ADDITIONAL SYSTEM.

C

C

C  STORAGE OF SPARSE MATRICES

C

C    THE NONZERO ENTRIES OF THE MATRIX M ARE STORED ROW-BY-ROW IN THE

C    ARRAY A.  TO IDENTIFY THE INDIVIDUAL NONZERO ENTRIES IN EACH ROW,

C    WE NEED TO KNOW IN WHICH COLUMN EACH ENTRY LIES.  THESE COLUMN

C    INDICES ARE STORED IN THE ARRAY JA;  I.E., IF  A(K) = M(I,J),  THEN

C    JA(K) = J.  TO IDENTIFY THE INDIVIDUAL ROWS, WE NEED TO KNOW WHERE

C    EACH ROW STARTS.  THESE ROW POINTERS ARE STORED IN THE ARRAY IA;

C    I.E., IF M(I,J) IS THE FIRST NONZERO ENTRY (STORED) IN THE I-TH ROW

C    AND  A(K) = M(I,J),  THEN  IA(I) = K.  MOREOVER, IA(N+1) POINTS TO

C    THE FIRST LOCATION FOLLOWING THE LAST ELEMENT IN THE LAST ROW.

C    THUS, THE NUMBER OF ENTRIES IN THE I-TH ROW IS  IA(I+1) - IA(I),

C    THE NONZERO ENTRIES IN THE I-TH ROW ARE STORED CONSECUTIVELY IN

C

C            A(IA(I)),  A(IA(I)+1),  ..., A(IA(I+1)-1),

C

C    AND THE CORRESPONDING COLUMN INDICES ARE STORED CONSECUTIVELY IN

C

C            JA(IA(I)), JA(IA(I)+1), ..., JA(IA(I+1)-1).

C

C    SINCE THE COEFFICIENT MATRIX IS SYMMETRIC, ONLY THE NONZERO ENTRIES

C    IN THE UPPER TRIANGLE NEED BE STORED, FOR EXAMPLE, THE MATRIX

C

C             ( 1  0  2  3  0 )

C             ( 0  4  0  0  0 )

C         M = ( 2  0  5  6  0 )

C             ( 3  0  6  7  8 )

C             ( 0  0  0  8  9 )

C

C    COULD BE STORED AS

C

C              1  2  3  4  5  6  7  8  9 10 11 12 13

C         ---+--------------------------------------

C         IA   1  4  5  8 12 14

C         JA   1  3  4  2  1  3  4  1  3  4  5  4  5

C          A   1  2  3  4  2  5  6  3  6  7  8  8  9

C

C    OR (SYMMETRICALLY) AS

C

C              1  2  3  4  5  6  7  8  9

C         ---+--------------------------

C         IA   1  4  5  7  9 10

C         JA   1  3  4  2  3  4  4  5  5

C          A   1  2  3  4  5  6  7  8  9          .

C

C

C  REORDERING THE ROWS AND COLUMNS OF M

C

C    A SYMMETRIC PERMUTATION OF THE ROWS AND COLUMNS OF THE COEFFICIENT

C    MATRIX M (E.G., WHICH REDUCES FILLIN OR ENHANCES NUMERICAL

C    STABILITY) MUST BE SPECIFIED.  THE SOLUTION Z IS RETURNED IN THE

C    ORIGINAL ORDER.

C

C    TO SPECIFY THE TRIVIAL ORDERING (I.E., THE IDENTITY PERMUTATION),

C    SET  P(I) = IP(I) = I,  I=1,...,N.  IN THIS CASE, P AND IP CAN BE

C    THE SAME ARRAY.

C

C    IF A NONTRIVIAL ORDERING (I.E., NOT THE IDENTITY PERMUTATION) IS

C    SPECIFIED AND M IS STORED SYMMETRICALLY (I.E., NOT BOTH M(I,J) AND

C    M(J,I) ARE STORED FOR I NE J), THEN ODRV SHOULD BE CALLED (WITH

C    PATH = 3 OR 5) TO SYMMETRICALLY REORDER (IA,JA,A) BEFORE CALLING

C    SDRV.  THIS IS TO ENSURE THAT IF M(I,J) WILL BE IN THE UPPER

C    TRIANGLE OF M WITH RESPECT TO THE NEW ORDERING, THEN M(I,J) IS

C    STORED IN ROW I (AND THUS M(J,I) IS NOT STORED);  WHEREAS IF M(I,J)

C    WILL BE IN THE STRICT LOWER TRIANGLE OF M, THEN M(J,I) IS STORED IN

C    ROW J (AND THUS M(I,J) IS NOT STORED).

C

C

C  PARAMETERS

C

C    N    - NUMBER OF VARIABLES/EQUATIONS

C

C    P    - INTEGER ONE-DIMENSIONAL ARRAY SPECIFYING A PERMUTATION OF

C           THE ROWS AND COLUMNS OF M;  DIMENSION = N

C

C    IP   - INTEGER ONE-DIMENSIONAL ARRAY CONTAINING THE INVERSE OF THE

C           PERMUTATION SPECIFIED IN P;  I.E., IP(P(I)) = I, I=1,...,N;

C           DIMENSION = N

C

C    IA   - INTEGER ONE-DIMENSIONAL ARRAY CONTAINING POINTERS TO DELIMIT

C           ROWS IN JA AND A;  DIMENSION = N+1

C

C    JA   - INTEGER ONE-DIMENSIONAL ARRAY CONTAINING THE COLUMN INDICES

C           CORRESPONDING TO THE ELEMENTS OF A;  DIMENSION = NUMBER OF

C           NONZERO ENTRIES IN M STORED

C

C    A    - REAL ONE-DIMENSIONAL ARRAY CONTAINING THE NONZERO ENTRIES IN

C           THE COEFFICIENT MATRIX M, STORED BY ROWS;  DIMENSION =

C           NUMBER OF NONZERO ENTRIES IN M STORED

C

C    B    - REAL ONE-DIMENSIONAL ARRAY CONTAINING THE RIGHT-HAND SIDE B;

C           B AND Z CAN BE THE SAME ARRAY;  DIMENSION = N

C

C    Z    - REAL ONE-DIMENSIONAL ARRAY CONTAINING THE SOLUTION X;  Z AND

C           B CAN BE THE SAME ARRAY;  DIMENSION = N

C

C    NSP  - DECLARED DIMENSION OF THE ONE-DIMENSIONAL ARRAYS ISP AND

C           RSP;  NSP MUST BE (SUBSTANTIALLY) LARGER THAN  3N+2K,  WHERE

C           K = NUMBER OF NONZERO ENTRIES IN THE UPPER TRIANGLE OF M

C

C    ISP  - INTEGER ONE-DIMENSIONAL ARRAY USED FOR WORKING STORAGE;  ISP

C           AND RSP SHOULD BE EQUIVALENCED;  DIMENSION = NSP

C

C    RSP  - REAL ONE-DIMENSIONAL ARRAY USED FOR WORKING STORAGE;  RSP

C           AND ISP SHOULD BE EQUIVALENCED;  DIMENSION = NSP

C

C    ESP  - INTEGER VARIABLE;  IF SUFFICIENT STORAGE WAS AVAILABLE TO

C           PERFORM THE SYMBOLIC FACTORIZATION (SSF), THEN ESP IS SET TO

C           THE AMOUNT OF EXCESS STORAGE PROVIDED (NEGATIVE IF

C           INSUFFICIENT STORAGE WAS AVAILABLE TO PERFORM THE NUMERIC

C           FACTORIZATION (SNF))

C

C    PATH - INTEGER PATH SPECIFICATION;  VALUES AND THEIR MEANINGS ARE -

C             1  PERFORM SSF, SNF, AND SNS

C             2  PERFORM SNF AND SNS (ISP/RSP IS ASSUMED TO HAVE BEEN

C                  SET UP IN AN EARLIER CALL TO SDRV (FOR SSF))

C             3  PERFORM SNS ONLY (ISP/RSP IS ASSUMED TO HAVE BEEN SET

C                  UP IN AN EARLIER CALL TO SDRV (FOR SSF AND SNF))

C             4  PERFORM SSF

C             5  PERFORM SSF AND SNF

C             6  PERFORM SNF ONLY (ISP/RSP IS ASSUMED TO HAVE BEEN SET

C                  UP IN AN EARLIER CALL TO SDRV (FOR SSF))

C

C    FLAG - INTEGER ERROR FLAG;  VALUES AND THEIR MEANINGS ARE -

C               0     NO ERRORS DETECTED

C              2N+K   DUPLICATE ENTRY IN A  --  ROW = K

C              6N+K   INSUFFICIENT STORAGE IN SSF  --  ROW = K

C              7N+1   INSUFFICIENT STORAGE IN SNF

C              8N+K   ZERO PIVOT  --  ROW = K

C             10N+1   INSUFFICIENT STORAGE IN SDRV

C             11N+1   ILLEGAL PATH SPECIFICATION

C

C

C  CONVERSION FROM REAL TO DOUBLE PRECISION

C

C    CHANGE THE REAL DECLARATIONS IN SDRV, SNF, AND SNS TO DOUBLE

C    PRECISION DECLARATIONS;  AND CHANGE THE VALUE IN THE DATA STATEMENT

C    FOR THE INTEGER VARIABLE RATIO (IN SDRV) FROM 1 TO 2.

C

C-----------------------------------------------------------------------

C

        INTEGER  p(*), ip(*),  ia(*), ja(*),  isp(*), esp,  path,  flag,

     *     ratio,  q, mark, d, u, tmp,  umax

C...    REAL  A(1),  B(1),  Z(1),  RSP(1)

C...    DATA  RATIO/1/

        DOUBLE PRECISION  a(*),  b(*),  z(*),  rsp(*)

        DATA  ratio/1/

C

C----VALIDATE PATH SPECIFICATION

        IF (path.LT.1 .OR. 6.LT.path)  go to 111

C

C----ALLOCATE STORAGE AND FACTOR M SYMBOLICALLY TO DETERMINE FILL-IN

        iju   = 1

        iu    = iju     +  n

        jl    = iu      + n+1

        ju    = jl      +  n

        q     = (nsp+1) -  n

        mark  = q       -  n

        jumax = mark    - ju

C

        IF ((path-1) * (path-4) * (path-5) .NE. 0)  go to 1

          IF (jumax.LE.0)  go to 110

          CALL ssf

     *       (n,  p, ip,  ia, ja,  isp(iju), isp(ju), isp(iu), jumax,

     *        isp(q), isp(mark), isp(jl), flag)

          IF (flag.NE.0)  go to 100

C

C----ALLOCATE STORAGE AND FACTOR M NUMERICALLY

   1    il   = ju      + isp(iju+(n-1))

        tmp  = ((il-1)+(ratio-1)) / ratio  +  1

        d    = tmp     + n

        u    = d       + n

        umax = (nsp+1) - u

        esp  = umax    - (isp(iu+n)-1)

C

        IF ((path-1) * (path-2) * (path-5) * (path-6) .NE. 0)  go to 2

          IF (umax.LE.0)  go to 110

          CALL snf

     *       (n,  p, ip,  ia, ja, a,

     *        rsp(d),  isp(iju), isp(ju), isp(iu), rsp(u), umax,

     *        isp(il),  isp(jl),  flag)

          IF (flag.NE.0)  go to 100

C

C----SOLVE SYSTEM OF LINEAR EQUATIONS  MX = B

   2    IF ((path-1) * (path-2) * (path-3) .NE. 0)  go to 3

          IF (umax.LE.0)  go to 110

          CALL sns

     *       (n,  p,  rsp(d), isp(iju), isp(ju), isp(iu), rsp(u),  z, b,

     *        rsp(tmp))

C

   3    RETURN

C

C ** ERROR -- ERROR DETECTED IN SSF, SNF, OR SNS

 100    RETURN

C ** ERROR -- INSUFFICIENT STORAGE

 110    flag = 10*n + 1

        RETURN

C ** ERROR -- ILLEGAL PATH SPECIFICATION

 111    flag = 11*n + 1

        RETURN

        END

C

C

C***********************************************************************

C  INTERNAL DOCUMENTATION FOR SSF, SNF, AND SNS

C***********************************************************************

C

C  COMPRESSED STORAGE OF SPARSE MATRICES

C

C    THE STRICT UPPER TRIANGULAR PORTION OF THE MATRIX U IS STORED IN

C    (IA,JA,A) FORMAT USING THE ARRAYS IU, JU, AND U, EXCEPT THAT AN

C    ADDITIONAL ARRAY IJU IS USED TO REDUCE THE STORAGE REQUIRED FOR JU

C    BY ALLOWING SOME SEQUENCES OF COLUMN INDICES TO CORRESPOND TO MORE

C    THAN ONE ROW.  FOR I < N, IJU(I) IS THE INDEX IN JU OF THE FIRST

C    ENTRY FOR THE I-TH ROW;  IJU(N) IS THE NUMBER OF ENTRIES IN JU.

C    THUS, THE NUMBER OF ENTRIES IN THE I-TH ROW IS  IU(I+1) - IU(I),

C    THE NONZERO ENTRIES OF THE I-TH ROW ARE STORED CONSECUTIVELY IN

C

C        U(IU(I)),   U(IU(I)+1),   ..., U(IU(I+1)-1),

C

C    AND THE CORRESPONDING COLUMN INDICES ARE STORED CONSECUTIVELY IN

C

C        JU(IJU(I)), JU(IJU(I)+1), ..., JU(IJU(I)+IU(I+1)-IU(I)+1).

C

C    COMPRESSION IN JU OCCURS IN TWO WAYS.  FIRST, IF A ROW I WAS MERGED

C    INTO ROW K, AND THE NUMBER OF ELEMENTS MERGED IN FROM (THE TAIL

C    PORTION OF) ROW I IS THE SAME AS THE FINAL LENGTH OF ROW K, THEN

C    THE KTH ROW AND THE TAIL OF ROW I ARE IDENTICAL AND IJU(K) POINTS

C    TO THE START OF THE TAIL.  SECOND, IF SOME TAIL PORTION OF THE

C    (K-1)ST ROW IS IDENTICAL TO THE HEAD OF THE KTH ROW, THEN IJU(K)

C    POINTS TO THE START OF THAT TAIL PORTION.  FOR EXAMPLE, THE NONZERO

C    STRUCTURE OF THE STRICT UPPER TRIANGULAR PART OF THE MATRIX

C

C             ( D 0 0 0 X X X )

C             ( 0 D 0 X X 0 0 )

C             ( 0 0 D 0 X X 0 )

C         U = ( 0 0 0 D X X 0 )

C             ( 0 0 0 0 D X X )

C             ( 0 0 0 0 0 D X )

C             ( 0 0 0 0 0 0 D )

C

C    WOULD BE STORED AS

C

C               1  2  3  4  5  6  7  8

C         ----+------------------------

C          IU   1  4  6  8 10 12 12 12

C          JU   5  6  7  4  5  6

C         IJU   1  4  5  5  2  3  6           .

C

C    THE DIAGONAL ENTRIES OF U ARE EQUAL TO ONE AND ARE NOT STORED.

C

C

C  ADDITIONAL PARAMETERS

C

C    D     - REAL ONE-DIMENSIONAL ARRAY CONTAINING THE RECIPROCALS OF

C            THE DIAGONAL ENTRIES OF THE MATRIX D;  DIMENSION = N

C

C    IJU   - INTEGER ONE-DIMENSIONAL ARRAY CONTAINING POINTERS TO THE

C            START OF EACH ROW IN JU;  DIMENSION = N

C

C    IU    - INTEGER ONE-DIMENSIONAL ARRAY CONTAINING POINTERS TO

C            DELIMIT ROWS IN U;  DIMENSION = N+1

C

C    JU    - INTEGER ONE-DIMENSIONAL ARRAY CONTAINING THE COLUMN INDICES

C            CORRESPONDING TO THE ELEMENTS OF U;  DIMENSION = JUMAX

C

C    JUMAX - DECLARED DIMENSION OF THE ONE-DIMENSIONAL ARRAY JU;  JUMAX

C            MUST BE AT LEAST THE SIZE OF U MINUS COMPRESSION (IJU(N)

C            AFTER THE CALL TO SSF)

C

C    U     - REAL ONE-DIMENSIONAL ARRAY CONTAINING THE NONZERO ENTRIES

C            IN THE STRICT UPPER TRIANGLE OF U, STORED BY ROWS;

C            DIMENSION = UMAX

C

C    UMAX  - DECLARED DIMENSION OF THE ONE-DIMENSIONAL ARRAY U;  UMAX

C            MUST BE AT LEAST THE NUMBER OF NONZERO ENTRIES IN THE

C            STRICT UPPER TRIANGLE OF M PLUS FILLIN (IU(N+1)-1 AFTER THE

C            CALL TO SSF)

C

C

C***********************************************************************

C  SSF --  SYMBOLIC UT-D-U FACTORIZATION OF SPARSE SYMMETRIC MATRIX

C***********************************************************************

        SUBROUTINE  ssf

     *     (n, p,ip, ia,ja, iju,ju,iu,jumax, q, mark, jl, flag)

C

C  ADDITIONAL PARAMETERS

C

C    Q     - INTEGER ONE-DIMENSIONAL WORK ARRAY;  DIMENSION = N

C

C    MARK  - INTEGER ONE-DIMENSIONAL WORK ARRAY;  DIMENSION = N

C

C    JL    - INTEGER ONE-DIMENSIONAL WORK ARRAY;  DIMENSION = N

C

C

C  DEFINITIONS OF INTERNAL PARAMETERS (DURING K-TH STAGE OF ELIMINATION)

C

C    Q CONTAINS AN ORDERED LINKED LIST REPRESENTATION OF THE NONZERO

C      STRUCTURE OF THE K-TH ROW OF U --

C        Q(K) IS THE FIRST COLUMN WITH A NONZERO ENTRY

C        Q(I) IS THE NEXT COLUMN WITH A NONZERO ENTRY AFTER COLUMN I

C      IN EITHER CASE, Q(I) = N+1 INDICATES THE END OF THE LIST

C

C    JL CONTAINS LISTS OF ROWS TO BE MERGED INTO UNELIMINATED ROWS --

C        I GE K => JL(I) IS THE FIRST ROW TO BE MERGED INTO ROW I

C        I LT K => JL(I) IS THE ROW FOLLOWING ROW I IN SOME LIST OF ROWS

C      IN EITHER CASE, JL(I) = 0 INDICATES THE END OF A LIST

C

C    MARK(I) IS THE LAST ROW STORED IN JU FOR WHICH U(MARK(I),I) NE 0

C

C    JUMIN AND JUPTR ARE THE INDICES IN JU OF THE FIRST AND LAST

C      ELEMENTS IN THE LAST ROW SAVED IN JU

C

C    LUK IS THE NUMBER OF NONZERO ENTRIES IN THE K-TH ROW

C

C-----------------------------------------------------------------------

C

        INTEGER  p(*), ip(*),  ia(*), ja(*),  iju(*), ju(*), iu(*),

     *     q(*),  mark(*),  jl(*),  flag,  tag, vj, qm

        LOGICAL  clique

C

C----INITIALIZATION

        jumin = 1

        juptr = 0

        iu(1) = 1

        DO 1 k=1,n

          mark(k) = 0

   1      jl(k) = 0

C

C----FOR EACH ROW K

        DO 18 k=1,n

          luk = 0

          q(k) = n+1

C

          tag = mark(k)

          clique = .false.

          IF (jl(k).NE.0)  clique = jl(jl(k)).EQ.0

C

C------INITIALIZE NONZERO STRUCTURE OF K-TH ROW TO ROW P(K) OF M

          jmin = ia(p(k))

          jmax = ia(p(k)+1) - 1

          IF (jmin.GT.jmax)  go to 4

          DO 3 j=jmin,jmax

            vj = ip(ja(j))

            IF (vj.LE.k)  go to 3

C

              qm = k

   2          m = qm

              qm = q(m)

              IF (qm.LT.vj)  go to 2

              IF (qm.EQ.vj)  go to 102

                luk = luk+1

                q(m) = vj

                q(vj) = qm

                IF (mark(vj).NE.tag)  clique = .false.

C

   3        CONTINUE

C

C------IF EXACTLY ONE ROW IS TO BE MERGED INTO THE K-TH ROW AND THERE IS

C------A NONZERO ENTRY IN EVERY COLUMN IN THAT ROW IN WHICH THERE IS A

C------NONZERO ENTRY IN ROW P(K) OF M, THEN DO NOT COMPUTE FILL-IN, JUST

C------USE THE COLUMN INDICES FOR THE ROW WHICH WAS TO HAVE BEEN MERGED

   4      IF (.NOT.clique)  go to 5

            iju(k) = iju(jl(k)) + 1

            luk = iu(jl(k)+1) - (iu(jl(k))+1)

            go to 17

C

C------MODIFY NONZERO STRUCTURE OF K-TH ROW BY COMPUTING FILL-IN

C------FOR EACH ROW I TO BE MERGED IN

   5      lmax = 0

          iju(k) = juptr

C

          i = k

   6      i = jl(i)

          IF (i.EQ.0)  go to 10

C

C--------MERGE ROW I INTO K-TH ROW

            lui = iu(i+1) - (iu(i)+1)

            jmin = iju(i) +  1

            jmax = iju(i) + lui

            qm = k

C

            DO 8 j=jmin,jmax

              vj = ju(j)

   7          m = qm

              qm = q(m)

              IF (qm.LT.vj)  go to 7

              IF (qm.EQ.vj)  go to 8

                luk = luk+1

                q(m) = vj

                q(vj) = qm

                qm = vj

   8          CONTINUE

C

C--------REMEMBER LENGTH AND POSITION IN JU OF LONGEST ROW MERGED

            IF (lui.LE.lmax)  go to 9

              lmax = lui

              iju(k) = jmin

C

   9        go to 6

C

C------IF THE K-TH ROW IS THE SAME LENGTH AS THE LONGEST ROW MERGED,

C------THEN USE THE COLUMN INDICES FOR THAT ROW

  10      IF (luk.EQ.lmax)  go to 17

C

C------IF THE TAIL OF THE LAST ROW SAVED IN JU IS THE SAME AS THE HEAD

C------OF THE K-TH ROW, THEN OVERLAP THE TWO SETS OF COLUMN INDICES --

C--------SEARCH LAST ROW SAVED FOR FIRST NONZERO ENTRY IN K-TH ROW ...

            i = q(k)

            IF (jumin.GT.juptr)  go to 12

            DO 11 jmin=jumin,juptr

              IF (ju(jmin)-i)  11, 13, 12

  11          CONTINUE

  12        go to 15

C

C--------... AND THEN TEST WHETHER TAIL MATCHES HEAD OF K-TH ROW

  13        iju(k) = jmin

            DO 14 j=jmin,juptr

              IF (ju(j).NE.i)  go to 15

              i = q(i)

              IF (i.GT.n)  go to 17

  14          CONTINUE

            juptr = jmin - 1

C

C------SAVE NONZERO STRUCTURE OF K-TH ROW IN JU

  15      i = k

          jumin = juptr +  1

          juptr = juptr + luk

          IF (juptr.GT.jumax)  go to 106

          DO 16 j=jumin,juptr

            i = q(i)

            ju(j) = i

  16        mark(i) = k

          iju(k) = jumin

C

C------ADD K TO ROW LIST FOR FIRST NONZERO ELEMENT IN K-TH ROW

  17      IF (luk.LE.1)  go to 18

            i = ju(iju(k))

            jl(k) = jl(i)

            jl(i) = k

C

  18      iu(k+1) = iu(k) + luk

C

        flag = 0

        RETURN

C

C ** ERROR -- DUPLICATE ENTRY IN A

 102    flag = 2*n + p(k)

        RETURN

C ** ERROR -- INSUFFICIENT STORAGE FOR JU

 106    flag = 6*n + k

        RETURN

        END

C

C***********************************************************************

C  SNF -- NUMERICAL UT-D-U FACTORIZATION OF SPARSE SYMMETRIC POSITIVE

C         DEFINITE MATRIX

C***********************************************************************

        SUBROUTINE  snf

     *     (n, p,ip, ia,ja,a, d, iju,ju,iu,u,umax, il, jl, flag)

C

C  ADDITIONAL PARAMETERS

C

C    IL    - INTEGER ONE-DIMENSIONAL WORK ARRAY;  DIMENSION = N

C

C    JL    - INTEGER ONE-DIMENSIONAL WORK ARRAY;  DIMENSION = N

C

C

C  DEFINITIONS OF INTERNAL PARAMETERS (DURING K-TH STAGE OF ELIMINATION)

C

C    (D(I),I=K,N) CONTAINS THE K-TH ROW OF U (EXPANDED)

C

C    IL(I) POINTS TO THE FIRST NONZERO ELEMENT IN COLUMNS K,...,N OF

C      ROW I OF U

C

C    JL CONTAINS LISTS OF ROWS TO BE ADDED TO UNELIMINATED ROWS --

C      I GE K => JL(I) IS THE FIRST ROW TO BE ADDED TO ROW I

C      I LT K => JL(I) IS THE ROW FOLLOWING ROW I IN SOME LIST OF ROWS

C      IN EITHER CASE, JL(I) = 0 INDICATES THE END OF A LIST

C

C-----------------------------------------------------------------------

C

        INTEGER  p(*), ip(*),  ia(*), ja(*),  iju(*), ju(*), iu(*),

     *     umax,  il(*),  jl(*),  flag,  vj

C...    REAL  A(1),  D(1), U(1),  DK, UKIDI

        DOUBLE PRECISION  a(*),  d(*), u(*),  dk, ukidi

C

C----CHECK FOR SUFFICIENT STORAGE FOR U

        IF (iu(n+1)-1 .GT. umax)  go to 107

C

C----INITIALIZATION

        DO 1 k=1,n

          kk=p(k)

          ip(kk)=k

          d(k) = 0

   1      jl(k) = 0

C

C----FOR EACH ROW K

        DO 11 k=1,n

C

C------INITIALIZE K-TH ROW WITH ELEMENTS NONZERO IN ROW P(K) OF M

   3      jmin = ia(p(k))

          jmax = ia(p(k)+1) - 1

          IF (jmin.GT.jmax) go to 5

          DO 4 j=jmin,jmax

            vj = ip(ja(j))

            IF (k.LE.vj)  d(vj) = a(j)

   4        CONTINUE

C

C------MODIFY K-TH ROW BY ADDING IN THOSE ROWS I WITH U(I,K) NE 0

C------FOR EACH ROW I TO BE ADDED IN

   5      dk = d(k)

          i = jl(k)

   6      IF (i.EQ.0)  go to 9

            nexti = jl(i)

C

C--------COMPUTE MULTIPLIER AND UPDATE DIAGONAL ELEMENT

            ili = il(i)

            ukidi = - u(ili) * d(i)

            dk = dk + ukidi * u(ili)

            u(ili) = ukidi

C

C--------ADD MULTIPLE OF ROW I TO K-TH ROW ...

            jmin = ili     + 1

            jmax = iu(i+1) - 1

            IF (jmin.GT.jmax)  go to 8

              mu = iju(i) - iu(i)

              DO 7 j=jmin,jmax

   7            d(ju(mu+j)) = d(ju(mu+j)) + ukidi * u(j)

C

C--------... AND ADD I TO ROW LIST FOR NEXT NONZERO ENTRY

              il(i) = jmin

              j = ju(mu+jmin)

              jl(i) = jl(j)

              jl(j) = i

C

   8        i = nexti

            go to 6

C

C------CHECK FOR ZERO PIVOT AND SAVE DIAGONAL ELEMENT

   9      IF (dk.EQ.0)  go to 108

          d(k) = 1 / dk

C

C------SAVE NONZERO ENTRIES IN K-TH ROW OF U ...

          jmin = iu(k)

          jmax = iu(k+1) - 1

          IF (jmin.GT.jmax)  go to 11

            mu = iju(k) - jmin

            DO 10 j=jmin,jmax

              jumuj = ju(mu+j)

              u(j) = d(jumuj)

  10          d(jumuj) = 0

C

C------... AND ADD K TO ROW LIST FOR FIRST NONZERO ENTRY IN K-TH ROW

            il(k) = jmin

            i = ju(mu+jmin)

            jl(k) = jl(i)

            jl(i) = k

  11      CONTINUE

C

        flag = 0

        RETURN

C

C ** ERROR -- INSUFFICIENT STORAGE FOR U

 107    flag = 7*n + 1

        RETURN

C ** ERROR -- ZERO PIVOT

 108    flag = 8*n + k

        RETURN

        END

C

C***********************************************************************

C  SNS -- SOLUTION OF SPARSE SYMMETRIC POSITIVE DEFINITE SYSTEM OF

C         LINEAR EQUATIONS  MX = B  GIVEN UT-D-U FACTORIZATION OF M

C***********************************************************************

        SUBROUTINE  sns

     *     (n, p, d, iju,ju,iu,u, z, b, tmp)

        INTEGER  p(*),  iju(*), ju(*), iu(*)

C...    REAL  D(1), U(1),  Z(1), B(1),  TMP(1),  TMPK, SUM

        DOUBLE PRECISION  d(*), u(*),  z(*), b(*),  tmp(*),  tmpk, sum

C

C  ADDITIONAL PARAMETERS

C

C    TMP   - REAL ONE-DIMENSIONAL WORK ARRAY;  DIMENSION = N

C

C-----------------------------------------------------------------------

C

C----SET TMP TO PERMUTED B

        DO 1 k=1,n

   1      tmp(k) = b(p(k))

C

C----SOLVE  UT D Y = B  BY FORWARD SUBSTITUTION

        DO 3 k=1,n

          tmpk = tmp(k)

          jmin = iu(k)

          jmax = iu(k+1) - 1

          IF (jmin.GT.jmax)  go to 3

          mu = iju(k) - jmin

          DO 2 j=jmin,jmax

   2        tmp(ju(mu+j)) = tmp(ju(mu+j)) + u(j) * tmpk

   3      tmp(k) = tmpk * d(k)

C

C----SOLVE  U X = Y  BY BACK SUBSTITUTION

        k = n

        DO 6 i=1,n

          sum = tmp(k)

          jmin = iu(k)

          jmax = iu(k+1) - 1

          IF (jmin.GT.jmax)  go to 5

          mu = iju(k) - jmin

          DO 4 j=jmin,jmax

   4        sum = sum + u(j) * tmp(ju(mu+j))

   5      tmp(k) = sum

          z(p(k)) = sum

   6      k = k-1

C

        RETURN

        END

      SUBROUTINE cdrvd

     *     (n, r,c,ic, ia,ja,a, b, z, nsp,isp,rsp,esp, path, flag)

C-----------------------------------------------------------------------

C*** DRIVER FOR SUBROUTINES FOR SOLVING SPARSE NONSYMMETRIC SYSTEMS OF

C       LINEAR EQUATIONS (COMPRESSED POINTER STORAGE)

C-----------------------------------------------------------------------

C

C

      INTEGER  r(*), c(*), ic(*),  ia(*), ja(*),  isp(*), esp,  path,

     *   flag,  d, u, q, row, tmp, ar,  umax

      DOUBLE PRECISION  a(*), b(*), z(*), rsp(*)

C

C

      DATA lratio/2/

C

      IF (path.LT.1 .OR. 5.LT.path)  go to 111

C

      il   = 1

      ijl  = il  + (n+1)

      iu   = ijl +   n

      iju  = iu  + (n+1)

      irl  = iju +   n

      jrl  = irl +   n

      jl   = jrl +   n

C

C

      IF ((path-1) * (path-5) .NE. 0)  go to 5

        max = (lratio*nsp + 1 - jl) - (n+1) - 5*n

        jlmax = max/2

        q     = jl   + jlmax

        ira   = q    + (n+1)

        jra   = ira  +   n

        irac  = jra  +   n

        iru   = irac +   n

        jru   = iru  +   n

        jutmp = jru  +   n

        jumax = lratio*nsp  + 1 - jutmp

        esp = max/lratio

        IF (jlmax.LE.0 .OR. jumax.LE.0)  go to 110

C

        DO 1 i=1,n

          IF (c(i).NE.i)  go to 2

 1        CONTINUE

        go to 3

 2      ar = nsp + 1 - n

        CALL nroc

     *     (n, ic, ia,ja,a, isp(il), rsp(ar), isp(iu), flag)

        IF (flag.NE.0)  go to 100

C

 3      CALL nsfc

     *     (n, r, ic, ia,ja,

     *      jlmax, isp(il), isp(jl), isp(ijl),

     *      jumax, isp(iu), isp(jutmp), isp(iju),

     *      isp(q), isp(ira), isp(jra), isp(irac),

     *      isp(irl), isp(jrl), isp(iru), isp(jru),  flag)

        IF(flag .NE. 0)  go to 100

C

        jlmax = isp(ijl+n-1)

        ju    = jl + jlmax

        jumax = isp(iju+n-1)

        IF (jumax.LE.0)  go to 5

        DO 4 j=1,jumax

 4        isp(ju+j-1) = isp(jutmp+j-1)

C

C

 5    jlmax = isp(ijl+n-1)

      ju    = jl  + jlmax

      jumax = isp(iju+n-1)

      l     = (ju + jumax - 2 + lratio)  /  lratio    +    1

      lmax  = isp(il+n) - 1

      d     = l   + lmax

      u     = d   + n

      row   = nsp + 1 - n

      tmp   = row - n

      umax  = tmp - u

      esp   = umax - (isp(iu+n) - 1)

C

      IF ((path-1) * (path-2) .NE. 0)  go to 6

        IF (umax.LE.0)  go to 110

        CALL nnfc

     *     (n,  r, c, ic,  ia, ja, a, z, b,

     *      lmax, isp(il), isp(jl), isp(ijl), rsp(l),  rsp(d),

     *      umax, isp(iu), isp(ju), isp(iju), rsp(u),

     *      rsp(row), rsp(tmp),  isp(irl), isp(jrl),  flag)

        IF(flag .NE. 0)  go to 100

C

 6    IF ((path-3) .NE. 0)  go to 7

        CALL nnsc

     *     (n,  r, c,  isp(il), isp(jl), isp(ijl), rsp(l),

     *      rsp(d),    isp(iu), isp(ju), isp(iju), rsp(u),

     *      z, b,  rsp(tmp))

C

 7    IF ((path-4) .NE. 0)  go to 8

        CALL nntc

     *     (n,  r, c,  isp(il), isp(jl), isp(ijl), rsp(l),

     *      rsp(d),    isp(iu), isp(ju), isp(iju), rsp(u),

     *      z, b,  rsp(tmp))

 8    RETURN

C

C

 100  RETURN

C

 110  flag = 10*n + 1

      RETURN

C

 111  flag = 11*n + 1

      RETURN

      END

C-----------------------------------------------------------------------

      SUBROUTINE nnfc

     *     (n, r,c,ic, ia,ja,a, z, b,

     *      lmax,il,jl,ijl,l, d, umax,iu,ju,iju,u,

     *      row, tmp, irl,jrl, flag)

C-----------------------------------------------------------------------

C*** SUBROUTINE NNFC

C*** NUMERICAL LDU-FACTORIZATION OF SPARSE NONSYMMETRIC MATRIX AND

C      SOLUTION OF SYSTEM OF LINEAR EQUATIONS (COMPRESSED POINTER

C      STORAGE)

C-----------------------------------------------------------------------

C

C

      INTEGER rk,umax

      INTEGER  r(*), c(*), ic(*), ia(*), ja(*), il(*), jl(*), ijl(*)

      INTEGER  iu(*), ju(*), iju(*), irl(*), jrl(*), flag

      DOUBLE PRECISION  a(*), l(*), d(*), u(*), z(*), b(*), row(*)

      DOUBLE PRECISION  tmp(*), lki, sum, dk

C

C

      IF(il(n+1)-1 .GT. lmax) go to 104

      IF(iu(n+1)-1 .GT. umax) go to 107

      DO 1 k=1,n

        irl(k) = il(k)

        jrl(k) = 0

 1      CONTINUE

C

C

      DO 19 k=1,n

C

        row(k) = 0

        i1 = 0

        IF (jrl(k) .EQ. 0) go to 3

        i = jrl(k)

 2      i2 = jrl(i)

        jrl(i) = i1

        i1 = i

        row(i) = 0

        i = i2

        IF (i .NE. 0) go to 2

C

 3      jmin = iju(k)

        jmax = jmin + iu(k+1) - iu(k) - 1

        IF (jmin .GT. jmax) go to 5

        DO 4 j=jmin,jmax

 4        row(ju(j)) = 0

C

 5      rk = r(k)

        jmin = ia(rk)

        jmax = ia(rk+1) - 1

        DO 6 j=jmin,jmax

          row(ic(ja(j))) = a(j)

 6        CONTINUE

C

        sum = b(rk)

        i = i1

        IF (i .EQ. 0) go to 10

C

 7        lki = -row(i)

C

          l(irl(i)) = -lki

          sum = sum + lki * tmp(i)

          jmin = iu(i)

          jmax = iu(i+1) - 1

          IF (jmin .GT. jmax) go to 9

          mu = iju(i) - jmin

          DO 8 j=jmin,jmax

 8          row(ju(mu+j)) = row(ju(mu+j)) + lki * u(j)

 9        i = jrl(i)

          IF (i .NE. 0) go to 7

C

C

 10     IF (row(k) .EQ. 0.0d0) go to 108

        dk = 1.0d0 / row(k)

        d(k) = dk

        tmp(k) = sum * dk

        IF (k .EQ. n) go to 19

        jmin = iu(k)

        jmax = iu(k+1) - 1

        IF (jmin .GT. jmax) go to 12

        mu = iju(k) - jmin

        DO 11 j=jmin,jmax

 11       u(j) = row(ju(mu+j)) * dk

 12     CONTINUE

C

C

        i = i1

        IF (i .EQ. 0) go to 18

 14     irl(i) = irl(i) + 1

        i1 = jrl(i)

        IF (irl(i) .GE. il(i+1)) go to 17

        ijlb = irl(i) - il(i) + ijl(i)

        j = jl(ijlb)

 15     IF (i .GT. jrl(j)) go to 16

          j = jrl(j)

          go to 15

 16     jrl(i) = jrl(j)

        jrl(j) = i

 17     i = i1

        IF (i .NE. 0) go to 14

 18     IF (irl(k) .GE. il(k+1)) go to 19

        j = jl(ijl(k))

        jrl(k) =jrl(j)

        jrl(j) = k

 19     CONTINUE

C

C

      k = n

      DO 22 i=1,n

        sum = tmp(k)

        jmin = iu(k)

        jmax = iu(k+1) - 1

        IF (jmin .GT. jmax)  go to 21

        mu = iju(k) - jmin

        DO 20 j=jmin,jmax

 20       sum = sum - u(j) * tmp(ju(mu+j))

 21     tmp(k) = sum

        z(c(k)) = sum

 22     k = k-1

      flag = 0

      RETURN

C

C

 104  flag = 4*n + 1

      RETURN

C

 107  flag = 7*n + 1

      RETURN

C

 108  flag = 8*n + k

      RETURN

      END

C-----------------------------------------------------------------------

      SUBROUTINE nnsc

     *     (n, r, c, il, jl, ijl, l, d, iu, ju, iju, u, z, b, tmp)

C-----------------------------------------------------------------------

C*** SUBROUTINE NNSC

C*** NUMERICAL SOLUTION OF SPARSE NONSYMMETRIC SYSTEM OF LINEAR

C      EQUATIONS GIVEN LDU-FACTORIZATION (COMPRESSED POINTER STORAGE)

C-----------------------------------------------------------------------

C

C

      INTEGER r(*), c(*), il(*), jl(*), ijl(*), iu(*), ju(*), iju(*)

      DOUBLE PRECISION  l(*), d(*), u(*), b(*), z(*), tmp(*), tmpk,sum

C

C

      DO 1 k=1,n

 1      tmp(k) = b(r(k))

C

      DO 3 k = 1,n

        jmin = il(k)

        jmax = il(k+1) - 1

        tmpk = -d(k) * tmp(k)

        tmp(k) = -tmpk

        IF (jmin .GT. jmax) go to 3

        ml = ijl(k) - jmin

        DO 2 j=jmin,jmax

 2        tmp(jl(ml+j)) = tmp(jl(ml+j)) + tmpk * l(j)

 3      CONTINUE

C

      k = n

      DO 6 i=1,n

        sum = -tmp(k)

        jmin = iu(k)

        jmax = iu(k+1) - 1

        IF (jmin .GT. jmax) go to 5

        mu = iju(k) - jmin

        DO 4 j=jmin,jmax

 4        sum = sum + u(j) * tmp(ju(mu+j))

 5      tmp(k) = -sum

        z(c(k)) = -sum

        k = k - 1

 6      CONTINUE

      RETURN

      END

C-----------------------------------------------------------------------

      SUBROUTINE nntc

     *     (n, r, c, il, jl, ijl, l, d, iu, ju, iju, u, z, b, tmp)

C-----------------------------------------------------------------------

C*** SUBROUTINE NNTC

C*** NUMERIC SOLUTION OF THE TRANSPOSE OF A SPARSE NONSYMMETRIC SYSTEM

C      OF LINEAR EQUATIONS GIVEN LU-FACTORIZATION (COMPRESSED POINTER

C      STORAGE)

C-----------------------------------------------------------------------

C

C

      INTEGER r(*), c(*), il(*), jl(*), ijl(*), iu(*), ju(*), iju(*)

      DOUBLE PRECISION l(*), d(*), u(*), b(*), z(*), tmp(*), tmpk,sum

C

C

      DO 1 k=1,n

 1      tmp(k) = b(c(k))

C

      DO 3 k=1,n

        jmin = iu(k)

        jmax = iu(k+1) - 1

        tmpk = -tmp(k)

        IF (jmin .GT. jmax) go to 3

        mu = iju(k) - jmin

        DO 2 j=jmin,jmax

 2        tmp(ju(mu+j)) = tmp(ju(mu+j)) + tmpk * u(j)

 3      CONTINUE

C

      k = n

      DO 6 i=1,n

        sum = -tmp(k)

        jmin = il(k)

        jmax = il(k+1) - 1

        IF (jmin .GT. jmax) go to 5

        ml = ijl(k) - jmin

        DO 4 j=jmin,jmax

 4        sum = sum + l(j) * tmp(jl(ml+j))

 5      tmp(k) = -sum * d(k)

        z(r(k)) = tmp(k)

        k = k - 1

 6      CONTINUE

      RETURN

      END

C-----------------------------------------------------------------------

      SUBROUTINE nroc (N, IC, IA, JA, A, JAR, AR, P, FLAG)

C-----------------------------------------------------------------------

C    YALE SPARSE MATRIX PACKAGE - NONSYMMETRIC CODES

C-----------------------------------------------------------------------

C*** SUBROUTINE NROC

C*** REORDERS ROWS OF A, LEAVING ROW ORDER UNCHANGED

C-----------------------------------------------------------------------

C

      INTEGER  ic(*), ia(*), ja(*), jar(*), p(*), flag

      DOUBLE PRECISION  a(*), ar(*)

C

C

      DO 5 k=1,n

        jmin = ia(k)

        jmax = ia(k+1) - 1

        IF(jmin .GT. jmax) go to 5

        p(n+1) = n + 1

C

        DO 3 j=jmin,jmax

          newj = ic(ja(j))

          i = n + 1

 1        IF(p(i) .GE. newj) go to 2

            i = p(i)

            go to 1

 2        IF(p(i) .EQ. newj) go to 102

          p(newj) = p(i)

          p(i) = newj

          jar(newj) = ja(j)

          ar(newj) = a(j)

 3        CONTINUE

C

        i = n + 1

        DO 4 j=jmin,jmax

          i = p(i)

          ja(j) = jar(i)

 4        a(j) = ar(i)

 5      CONTINUE

      flag = 0

      RETURN

C

C

 102  flag = n + k

      RETURN

      END

C-----------------------------------------------------------------------

      SUBROUTINE nsfc

     *      (n, r, ic, ia,ja, jlmax,il,jl,ijl, jumax,iu,ju,iju,

     *       q, ira,jra, irac, irl,jrl, iru,jru, flag)

C-----------------------------------------------------------------------

C*** SUBROUTINE NSFC

C*** SYMBOLIC LDU-FACTORIZATION OF NONSYMMETRIC SPARSE MATRIX

C      (COMPRESSED POINTER STORAGE)

C-----------------------------------------------------------------------

C

C

      INTEGER cend, qm, rend, rk, vj

      INTEGER ia(*), ja(*), ira(*), jra(*), il(*), jl(*), ijl(*)

      INTEGER iu(*), ju(*), iju(*), irl(*), jrl(*), iru(*), jru(*)

      INTEGER r(*), ic(*), q(*), irac(*), flag

C

C

      np1 = n + 1

      jlmin = 1

      jlptr = 0

      il(1) = 1

      jumin = 1

      juptr = 0

      iu(1) = 1

      DO 1 k=1,n

        irac(k) = 0

        jra(k) = 0

        jrl(k) = 0

 1      jru(k) = 0

C

      DO 2 k=1,n

        rk = r(k)

        iak = ia(rk)

        IF (iak .GE. ia(rk+1))  go to 101

        jaiak = ic(ja(iak))

        IF (jaiak .GT. k)  go to 105

        jra(k) = irac(jaiak)

        irac(jaiak) = k

 2      ira(k) = iak

C

C

      DO 41 k=1,n

C

C

        q(np1) = np1

        luk = -1

C

        vj = irac(k)

        IF (vj .EQ. 0)  go to 5

 3        qm = np1

 4        m = qm

          qm = q(m)

          IF (qm .LT. vj)  go to 4

          IF (qm .EQ. vj)  go to 102

            luk = luk + 1

            q(m) = vj

            q(vj) = qm

            vj = jra(vj)

            IF (vj .NE. 0)  go to 3

C

 5      lastid = 0

        lasti = 0

        ijl(k) = jlptr

        i = k

 6        i = jru(i)

          IF (i .EQ. 0)  go to 10

          qm = np1

          jmin = irl(i)

          jmax = ijl(i) + il(i+1) - il(i) - 1

          long = jmax - jmin

          IF (long .LT. 0)  go to 6

          jtmp = jl(jmin)

          IF (jtmp .NE. k)  long = long + 1

          IF (jtmp .EQ. k)  r(i) = -r(i)

          IF (lastid .GE. long)  go to 7

            lasti = i

            lastid = long

C

 7        DO 9 j=jmin,jmax

            vj = jl(j)

 8          m = qm

            qm = q(m)

            IF (qm .LT. vj)  go to 8

            IF (qm .EQ. vj)  go to 9

              luk = luk + 1

              q(m) = vj

              q(vj) = qm

              qm = vj

 9          CONTINUE

            go to 6

C

C

 10     qm = q(np1)

        IF (qm .NE. k)  go to 105

        IF (luk .EQ. 0)  go to 17

        IF (lastid .NE. luk)  go to 11

C

        irll = irl(lasti)

        ijl(k) = irll + 1

        IF (jl(irll) .NE. k)  ijl(k) = ijl(k) - 1

        go to 17

C

 11     IF (jlmin .GT. jlptr)  go to 15

        qm = q(qm)

        DO 12 j=jlmin,jlptr

          IF (jl(j) - qm)  12, 13, 15

 12       CONTINUE

        go to 15

 13     ijl(k) = j

        DO 14 i=j,jlptr

          IF (jl(i) .NE. qm)  go to 15

          qm = q(qm)

          IF (qm .GT. n)  go to 17

 14       CONTINUE

        jlptr = j - 1

C

 15     jlmin = jlptr + 1

        ijl(k) = jlmin

        IF (luk .EQ. 0)  go to 17

        jlptr = jlptr + luk

        IF (jlptr .GT. jlmax)  go to 103

          qm = q(np1)

          DO 16 j=jlmin,jlptr

            qm = q(qm)

 16         jl(j) = qm

 17     irl(k) = ijl(k)

        il(k+1) = il(k) + luk

C

C

        q(np1) = np1

        luk = -1

C

        rk = r(k)

        jmin = ira(k)

        jmax = ia(rk+1) - 1

        IF (jmin .GT. jmax)  go to 20

        DO 19 j=jmin,jmax

          vj = ic(ja(j))

          qm = np1

 18       m = qm

          qm = q(m)

          IF (qm .LT. vj)  go to 18

          IF (qm .EQ. vj)  go to 102

            luk = luk + 1

            q(m) = vj

            q(vj) = qm

 19       CONTINUE

C

 20     lastid = 0

        lasti = 0

        iju(k) = juptr

        i = k

        i1 = jrl(k)

 21       i = i1

          IF (i .EQ. 0)  go to 26

          i1 = jrl(i)

          qm = np1

          jmin = iru(i)

          jmax = iju(i) + iu(i+1) - iu(i) - 1

          long = jmax - jmin

          IF (long .LT. 0)  go to 21

          jtmp = ju(jmin)

          IF (jtmp .EQ. k)  go to 22

C

            long = long + 1

            cend = ijl(i) + il(i+1) - il(i)

            irl(i) = irl(i) + 1

            IF (irl(i) .GE. cend)  go to 22

              j = jl(irl(i))

              jrl(i) = jrl(j)

              jrl(j) = i

 22       IF (lastid .GE. long)  go to 23

            lasti = i

            lastid = long

C

 23       DO 25 j=jmin,jmax

            vj = ju(j)

 24         m = qm

            qm = q(m)

            IF (qm .LT. vj)  go to 24

            IF (qm .EQ. vj)  go to 25

              luk = luk + 1

              q(m) = vj

              q(vj) = qm

              qm = vj

 25         CONTINUE

          go to 21

C

 26     IF (il(k+1) .LE. il(k))  go to 27

          j = jl(irl(k))

          jrl(k) = jrl(j)

          jrl(j) = k

C

C

 27     qm = q(np1)

        IF (qm .NE. k)  go to 105

        IF (luk .EQ. 0)  go to 34

        IF (lastid .NE. luk)  go to 28

C

        irul = iru(lasti)

        iju(k) = irul + 1

        IF (ju(irul) .NE. k)  iju(k) = iju(k) - 1

        go to 34

C

 28     IF (jumin .GT. juptr)  go to 32

        qm = q(qm)

        DO 29 j=jumin,juptr

          IF (ju(j) - qm)  29, 30, 32

 29       CONTINUE

        go to 32

 30     iju(k) = j

        DO 31 i=j,juptr

          IF (ju(i) .NE. qm)  go to 32

          qm = q(qm)

          IF (qm .GT. n)  go to 34

 31       CONTINUE

        juptr = j - 1

C

 32     jumin = juptr + 1

        iju(k) = jumin

        IF (luk .EQ. 0)  go to 34

        juptr = juptr + luk

        IF (juptr .GT. jumax)  go to 106

          qm = q(np1)

          DO 33 j=jumin,juptr

            qm = q(qm)

 33         ju(j) = qm

 34     iru(k) = iju(k)

        iu(k+1) = iu(k) + luk

C

C

        i = k

 35       i1 = jru(i)

          IF (r(i) .LT. 0)  go to 36

          rend = iju(i) + iu(i+1) - iu(i)

          IF (iru(i) .GE. rend)  go to 37

            j = ju(iru(i))

            jru(i) = jru(j)

            jru(j) = i

            go to 37

 36       r(i) = -r(i)

 37       i = i1

          IF (i .EQ. 0)  go to 38

          iru(i) = iru(i) + 1

          go to 35

C

C

 38     i = irac(k)

        IF (i .EQ. 0)  go to 41

 39       i1 = jra(i)

          ira(i) = ira(i) + 1

          IF (ira(i) .GE. ia(r(i)+1))  go to 40

          irai = ira(i)

          jairai = ic(ja(irai))

          IF (jairai .GT. i)  go to 40

          jra(i) = irac(jairai)

          irac(jairai) = i

 40       i = i1

          IF (i .NE. 0)  go to 39

 41     CONTINUE

C

      ijl(n) = jlptr

      iju(n) = juptr

      flag = 0

      RETURN

C

C

 101  flag = n + rk

      RETURN

C

 102  flag = 2*n + rk

      RETURN

C

 103  flag = 3*n + k

      RETURN

C

 105  flag = 5*n + k

      RETURN

C

 106  flag = 6*n + k

      RETURN

      END

C-----------------------------------------------------------------------

      SUBROUTINE nnfcad

     *     (n, r,c,ic, ia,ja,a, z, b,

     *      lmax,il,jl,ijl, umax,iu,ju,iju,u,

     *      row, tmp, irl,jrl, flag)

C-----------------------------------------------------------------------

C*** SUBROUTINE NNFC

C*** NUMERICAL LDU-FACTORIZATION OF SPARSE NONSYMMETRIC MATRIX AND

C      SOLUTION OF SYSTEM OF LINEAR EQUATIONS (COMPRESSED POINTER

C      STORAGE)

C-----------------------------------------------------------------------

C

C

      INTEGER rk,umax

      INTEGER  r(*), c(*), ic(*), ia(*), ja(*), il(*), jl(*), ijl(*)

      INTEGER  iu(*), ju(*), iju(*), irl(*), jrl(*), flag

      REAL*8  a(*),  u(*), z(*), b(*), row(*)

      REAL*8 tmp(*), lki, sum, dk

C

C

      IF(il(n+1)-1 .GT. lmax) go to 104

      IF(iu(n+1)-1 .GT. umax) go to 107

      DO 1 k=1,n

        irl(k) = il(k)

        jrl(k) = 0

 1      CONTINUE

C

C

      DO 19 k=1,n

C

        row(k) = 0

        i1 = 0

        IF (jrl(k) .EQ. 0) go to 3

        i = jrl(k)

 2      i2 = jrl(i)

        jrl(i) = i1

        i1 = i

        row(i) = 0

        i = i2

        IF (i .NE. 0) go to 2

C

 3      jmin = iju(k)

        jmax = jmin + iu(k+1) - iu(k) - 1

        IF (jmin .GT. jmax) go to 5

        DO 4 j=jmin,jmax

 4        row(ju(j)) = 0

C

 5      rk = r(k)

        jmin = ia(rk)

        jmax = ia(rk+1) - 1

        DO 6 j=jmin,jmax

          row(ic(ja(j))) = a(j)

 6        CONTINUE

C

        sum = b(rk)

        i = i1

        IF (i .EQ. 0) go to 10

C

 7        lki = -row(i)

C

          sum = sum + lki * tmp(i)

          jmin = iu(i)

          jmax = iu(i+1) - 1

          IF (jmin .GT. jmax) go to 9

          mu = iju(i) - jmin

          DO 8 j=jmin,jmax

 8          row(ju(mu+j)) = row(ju(mu+j)) + lki * u(j)

 9        i = jrl(i)

          IF (i .NE. 0) go to 7

C

C

 10     IF (row(k) .EQ. 0.0d0) go to 108

        dk = 1.0d0 / row(k)

        tmp(k) = sum * dk

        IF (k .EQ. n) go to 19

        jmin = iu(k)

        jmax = iu(k+1) - 1

        IF (jmin .GT. jmax) go to 12

        mu = iju(k) - jmin

        DO 11 j=jmin,jmax

 11       u(j) = row(ju(mu+j)) * dk

 12     CONTINUE

C

C

        i = i1

        IF (i .EQ. 0) go to 18

 14     irl(i) = irl(i) + 1

        i1 = jrl(i)

        IF (irl(i) .GE. il(i+1)) go to 17

        ijlb = irl(i) - il(i) + ijl(i)

        j = jl(ijlb)

 15     IF (i .GT. jrl(j)) go to 16

          j = jrl(j)

          go to 15

 16     jrl(i) = jrl(j)

        jrl(j) = i

 17     i = i1

        IF (i .NE. 0) go to 14

 18     IF (irl(k) .GE. il(k+1)) go to 19

        j = jl(ijl(k))

        jrl(k) =jrl(j)

        jrl(j) = k

 19     CONTINUE

C

C

      k = n

      DO 22 i=1,n

        sum = tmp(k)

        jmin = iu(k)

        jmax = iu(k+1) - 1

        IF (jmin .GT. jmax)  go to 21

        mu = iju(k) - jmin

        DO 20 j=jmin,jmax

 20       sum = sum - u(j) * tmp(ju(mu+j))

 21     tmp(k) = sum

        z(c(k)) = sum

 22     k = k-1

      flag = 0

      RETURN

C

C

 104  flag = 4*n + 1

      RETURN

C

 107  flag = 7*n + 1

      RETURN

C

 108  flag = 8*n + k

      RETURN

      END

      SUBROUTINE adrvd

     *     (n, r,c,ic, ia,ja,a, b, z, nsp,isp,rsp,esp, path, flag)

C-----------------------------------------------------------------------

C*** DRIVER FOR SUBROUTINES FOR SOLVING SPARSE NONSYMMETRIC SYSTEMS OF

C       LINEAR EQUATIONS (COMPRESSED POINTER STORAGE)

C-----------------------------------------------------------------------

C

C

      INTEGER  r(*), c(*), ic(*),  ia(*), ja(*),  isp(*), esp,  path,

     *   flag, u,  q, row, tmp, ar,  umax

      REAL*8 a(*), b(*), z(*), rsp(*)

C

C

      DATA lratio/1/

C

C

      il   = 1

      ijl  = il  + (n+1)

      iu   = ijl +   n

      iju  = iu  + (n+1)

      irl  = iju +   n

      jrl  = irl +   n

      jl   = jrl +   n

C

C

        max = (lratio*nsp + 1 - jl) - (n+1) - 5*n

        jlmax = max/2

        q     = jl   + jlmax

        ira   = q    + (n+1)

        jra   = ira  +   n

        irac  = jra  +   n

        iru   = irac +   n

        jru   = iru  +   n

        jutmp = jru  +   n

        jumax = lratio*nsp  + 1 - jutmp

        esp = max/lratio

        IF (jlmax.LE.0 .OR. jumax.LE.0)  go to 110

C

        DO 1 i=1,n

          IF (c(i).NE.i)  go to 2

 1        CONTINUE

        go to 3

 2      ar = nsp + 1 - n

        CALL nroc

     *     (n, ic, ia,ja,a, isp(il), rsp(ar), isp(iu), flag)

        IF (flag.NE.0)  go to 100

C

 3      CALL nsfc

     *     (n, r, ic, ia,ja,

     *      jlmax, isp(il), isp(jl), isp(ijl),

     *      jumax, isp(iu), isp(jutmp), isp(iju),

     *      isp(q), isp(ira), isp(jra), isp(irac),

     *      isp(irl), isp(jrl), isp(iru), isp(jru),  flag)

        IF(flag .NE. 0)  go to 100

C

        jlmax = isp(ijl+n-1)

        ju    = jl + jlmax

        jumax = isp(iju+n-1)

        IF (jumax.LE.0)  go to 5

        DO 4 j=1,jumax

 4        isp(ju+j-1) = isp(jutmp+j-1)

C

C

 5    jlmax = isp(ijl+n-1)

      ju    = jl  + jlmax

      jumax = isp(iju+n-1)

      lmax = isp(il+n) - 1

      u     = (ju + jumax - 2 + lratio ) / lratio + 1

      umax = isp(iu+n) - 1

      row   = u + umax

      tmp   = row + n

      esp = nsp - (tmp + n)

      IF(esp.LT.0) goto 110

C

        CALL nnfcad

     *     (n,  r, c, ic,  ia, ja, a, z, b,

     *      lmax, isp(il), isp(jl), isp(ijl),

     *      umax, isp(iu), isp(ju), isp(iju), rsp(u),

     *      rsp(row), rsp(tmp),  isp(irl), isp(jrl),  flag)

        IF(flag .NE. 0)  go to 100

 8    RETURN

C

C

 100  RETURN

C

 110  flag = 10*n + 1

      RETURN

C

      END

      SUBROUTINE pdrvd

     *     (n, ia,ja,a, p,ip, nsp,isp, path, flag)

C-----------------------------------------------------------------------

C  ODRV -- DRIVER FOR SPARSE MATRIX REORDERING ROUTINES

C-----------------------------------------------------------------------

      INTEGER  ia(*), ja(*), p(*), ip(*), isp(*), path, flag,

     *   v, l, head,  tmp, q

      DOUBLE PRECISION  a(*)

      LOGICAL  dflag

C

C

      flag = 0

      IF (path.LT.1 .OR. 5.LT.path) go to 111

C

C

      IF ((path-1) * (path-2) * (path-4) .NE. 0) go to 1

        max = (nsp-n)/2

        v    = 1

        l    = v     +  max

        head = l     +  max

        next = head  +  n

        IF (max.LT.n)  go to 110

C

        CALL  mdn

     *     (n, ia,ja, max,isp(v),isp(l), isp(head),p,ip, isp(v), flag)

        IF (flag.NE.0)  go to 100

C

C

 1    IF ((path-2) * (path-3) * (path-4) * (path-5) .NE. 0)  go to 2

        tmp = (nsp+1) -      n

        q   = tmp     - (ia(n+1)-1)

        IF (q.LT.1)  go to 110

C

        dflag = path.EQ.4 .OR. path.EQ.5

        CALL sro

     *     (n,  ip,  ia, ja, a,  isp(tmp),  isp(q),  dflag)

C

 2    RETURN

C

C ** ERROR -- ERROR DETECTED IN MD

 100  RETURN

C ** ERROR -- INSUFFICIENT STORAGE

 110  flag = 10*n + 1

      RETURN

C ** ERROR -- ILLEGAL PATH SPECIFIED

 111  flag = 11*n + 1

      RETURN

      END

      SUBROUTINE mdn

     *     (n, ia,ja, max, v,l, head,last,next, mark, flag)

C-----------------------------------------------------------------------

C  MD -- MINIMUM DEGREE ALGORITHM (BASED ON ELEMENT MODEL)

C-----------------------------------------------------------------------

      INTEGER  ia(*), ja(*),  v(*), l(*),  head(*), last(*), next(*),

     *   mark(*),  flag,  tag, dmin, vk,ek, tail

      equivalence(vk,ek)

C

C

      tag = 0

      CALL  mdin

     *   (n, ia,ja, max,v,l, head,last,next, mark,tag, flag)

      IF (flag.NE.0)  RETURN

C

      k = 0

      dmin = 1

C

C

 1    IF (k.GE.n)  go to 4

C

C

 2      IF (head(dmin).GT.0)  go to 3

          dmin = dmin + 1

          go to 2

C

C

 3      vk = head(dmin)

        head(dmin) = next(vk)

        IF (head(dmin).GT.0)  last(head(dmin)) = -dmin

C

C

        k = k+1

        next(vk) = -k

        last(ek) = dmin - 1

        tag = tag + last(ek)

        mark(vk) = tag

C

C

        CALL  mdm

     *     (vk,tail, v,l, last,next, mark)

C

C

        CALL  mdp

     *     (k,ek,tail, v,l, head,last,next, mark)

C

C

        CALL  mdu

     *     (ek,dmin, v,l, head,last,next, mark)

C

        go to 1

C

C

 4    DO 5 k=1,n

        next(k) = -next(k)

 5      last(next(k)) = k

C

      RETURN

      END

C-----------------------------------------------------------------------

      SUBROUTINE mdin

     *     (n, ia,ja, max,v,l, head,last,next, mark,tag, flag)

C-----------------------------------------------------------------------

C  MDI -- INITIALIZATION

C-----------------------------------------------------------------------

      INTEGER  ia(*), ja(*),  v(*), l(*),  head(*), last(*), next(*),

     *   mark(*), tag,  flag,  sfs, vi,dvi, vj

C

C

      DO 1 vi=1,n

        mark(vi) = 1

        l(vi) = 0

 1      head(vi) = 0

      sfs = n+1

C

C

      DO 6 vi=1,n

        jmin = ia(vi)

        jmax = ia(vi+1) - 1

        IF (jmin.GT.jmax)  go to 6

        DO 5 j=jmin,jmax

          vj = ja(j)

          IF (vj-vi) 2, 5, 4

C

C

 2        lvk = vi

          kmax = mark(vi) - 1

          IF (kmax .EQ. 0) go to 4

          DO 3 k=1,kmax

            lvk = l(lvk)

            IF (v(lvk).EQ.vj) go to 5

 3          CONTINUE

C

 4          IF (sfs.GE.max)  go to 101

C

C

            mark(vi) = mark(vi) + 1

            v(sfs) = vj

            l(sfs) =l(vi)

            l(vi) = sfs

            sfs = sfs+1

C

C

            mark(vj) = mark(vj) + 1

            v(sfs) = vi

            l(sfs) = l(vj)

            l(vj) = sfs

            sfs = sfs+1

 5        CONTINUE

 6      CONTINUE

C

C

      DO 7 vi=1,n

        dvi = mark(vi)

        next(vi) = head(dvi)

        head(dvi) =vi

        last(vi) = -dvi

        nextvi = next(vi)

        IF (nextvi.GT.0)  last(nextvi) = vi

 7      mark(vi) = tag

C

      RETURN

C

C

 101  flag = 9*n + vi

      RETURN

      END