ma27ad.f

  710   CONTINUE
C WE COPY IN REVERSE DIRECTION TO AVOID OVERWRITE PROBLEMS.
        LAELL = ((LIELL+1)*LIELL)/2
        KK = POSFAC + LAELL
        IF (KK.NE.ASTK) GO TO 720
        ASTK = ASTK - LAELL
        GO TO 740
C THE MOVE AND ZEROING OF ARRAY A IS PERFORMED WITH TWO LOOPS SO
C THAT THE CRAY-1 WILL VECTORIZE THEM SAFELY.
  720   KMAX = KK - 1
C THE FOLLOWING SPECIAL COMMENT FORCES VECTORIZATION ON THE CRAY-1.
CDIR$ IVDEP
        DO 730 K = 1,LAELL
          KK = KK - 1
          ASTK = ASTK - 1
          A(ASTK) = A(KK)
  730   CONTINUE
        KMAX = MIN(KMAX,ASTK-1)
        DO 735 K = KK,KMAX
          A(K) = ZERO
  735   CONTINUE
  740   AZERO = MIN(AZERO,ASTK-1)
  750   IF (NPIV.EQ.0) IWPOS = IOLDPS
  760 CONTINUE
C
C END OF LOOP ON TREE NODES.
C
      IW(1) = NBLK
      IF (NTWO.GT.0) IW(1) = -NBLK
      NRLBDU = POSFAC - 1
      NIRBDU = IWPOS - 1
      IF (NTOTPV.EQ.N) GO TO 810
      INFO(1) = 3
      INFO(2) = NTOTPV
      GO TO 810
C **** ERROR RETURNS ****
  770 INFO(1) = -3
      GO TO 810

  780 INFO(1) = -4
      INFO(2) = LA + MAX(POSFAC+LNASS,APOS2-LTOPST+2) - ASTK
      GO TO 810

  790 INFO(1) = -5
      INFO(2) = NTOTPV + 1
      GO TO 810

  800 INFO(1) = -6
      INFO(2) = NTOTPV + 1
  810 CONTINUE
      INFO(9) = NRLBDU
      INFO(10) = NIRBDU
      INFO(12) = NCMPBR
      INFO(13) = NCMPBI
      INFO(14) = NTWO
      INFO(15) = NEIG

      RETURN
      END
      SUBROUTINE MA27PD(A,IW,J1,J2,ITOP,IREAL,NCMPBR,NCMPBI)
C THIS SUBROUTINE PERFORMS A VERY SIMPLE COMPRESS (BLOCK MOVE).
C     ENTRIES J1 TO J2 (INCL.) IN A OR IW AS APPROPRIATE ARE MOVED TO
C     OCCUPY THE POSITIONS IMMEDIATELY PRIOR TO POSITION ITOP.
C A/IW HOLD THE ARRAY BEING COMPRESSED.
C J1/J2 DEFINE THE ENTRIES BEING MOVED.
C ITOP DEFINES THE POSITION IMMEDIATELY AFTER THE POSITIONS TO WHICH
C     J1 TO J2 ARE MOVED.
C IREAL MUST BE SET BY THE USER TO 2 IF THE MOVE IS ON ARRAY IW,
C     ANY OTHER VALUE WILL PERFORM THE MOVE ON A.
C NCMPBR and NCMPBI, see INFO(12) and INFO(13) in MA27A/AD (ACCUMULATE
C     THE NUMBER OF COMPRESSES OF THE REALS AND INTEGERS PERFORMED BY
C     MA27B/BD.
C
C     .. Scalar Arguments ..
      INTEGER IREAL,ITOP,J1,J2,NCMPBR,NCMPBI
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION A(*)
      INTEGER IW(*)
C     ..
C     .. Local Scalars ..
      INTEGER IPOS,JJ,JJJ
C     ..
C     .. Executable Statements ..
      IPOS = ITOP - 1
      IF (J2.EQ.IPOS) GO TO 50
      IF (IREAL.EQ.2) GO TO 20
      NCMPBR = NCMPBR + 1
      IF (J1.GT.J2) GO TO 40
      DO 10 JJJ = J1,J2
        JJ = J2 - JJJ + J1
        A(IPOS) = A(JJ)
        IPOS = IPOS - 1
   10 CONTINUE
      GO TO 40

   20 NCMPBI = NCMPBI + 1
      IF (J1.GT.J2) GO TO 40
      DO 30 JJJ = J1,J2
        JJ = J2 - JJJ + J1
        IW(IPOS) = IW(JJ)
        IPOS = IPOS - 1
   30 CONTINUE
   40 J2 = ITOP - 1
      J1 = IPOS + 1
   50 RETURN

      END
      SUBROUTINE MA27QD(N,A,LA,IW,LIW,W,MAXFNT,RHS,IW2,NBLK,LATOP,ICNTL)
C THIS SUBROUTINE PERFORMS FORWARD ELIMINATION
C     USING THE FACTOR U TRANSPOSE STORED IN A/IA AFTER MA27B/BD.
C
C N      - MUST BE SET TO THE ORDER OF THE MATRIX. NOT ALTERED
C          BY MA27Q/QD.
C A      - MUST BE SET TO HOLD THE REAL VALUES
C          CORRESPONDING TO THE FACTORS OF DINVERSE AND U.  THIS MUST BE
C          UNCHANGED SINCE THE PRECEDING CALL TO MA27B/BD.  NOT ALTERED
C          BY MA27Q/QD.
C LA     - LENGTH OF ARRAY A.  NOT ALTERED BY MA27Q/QD.
C IW     - HOLDS THE INTEGER INDEXING
C          INFORMATION FOR THE MATRIX FACTORS IN A.  THIS MUST BE
C          UNCHANGED SINCE THE PRECEDING CALL TO MA27B/BD.  NOT ALTERED
C          BY MA27Q/QD.
C LIW    - LENGTH OF ARRAY IW.  NOT ALTERED BY MA27Q/QD.
C W      - USED
C          AS WORKSPACE BY MA27Q/QD TO HOLD THE COMPONENTS OF THE RIGHT
C          HAND SIDES CORRESPONDING TO CURRENT BLOCK PIVOTAL ROWS.
C MAXFNT - MUST BE SET TO THE LARGEST NUMBER OF
C          VARIABLES IN ANY BLOCK PIVOT ROW.  THIS VALUE WILL HAVE
C          BEEN OUTPUT BY MA27B/BD.  NOT ALTERED BY MA27Q/QD.
C RHS    - ON INPUT,
C          MUST BE SET TO HOLD THE RIGHT HAND SIDES FOR THE EQUATIONS
C          WHICH THE USER DESIRES TO SOLVE.  ON OUTPUT, RHS WILL HOLD
C          THE MODIFIED VECTORS CORRESPONDING TO PERFORMING FORWARD
C          ELIMINATION ON THE RIGHT HAND SIDES.
C IW2    - NEED NOT BE SET ON ENTRY. ON EXIT IW2(I) (I = 1,NBLK)
C          WILL HOLD POINTERS TO THE
C          BEGINNING OF EACH BLOCK PIVOT IN ARRAY IW.
C NBLK   - NUMBER OF BLOCK PIVOT ROWS. NOT ALTERED BY MA27Q/QD.
C LATOP  - NEED NOT BE SET ON ENTRY. ON EXIT, IT IS THE POSITION IN
C          A OF THE LAST ENTRY IN THE  FACTORS. IT MUST BE PASSED
C          UNCHANGED TO MA27R/RD.
C ICNTL is an INTEGER array of length 30, see MA27A/AD.
C   ICNTL(IFRLVL+I) I=1,20 IS USED TO CONTROL WHETHER DIRECT OR INDIRECT
C     ACCESS IS USED BY MA27C/CD.  INDIRECT ACCESS IS EMPLOYED
C     IN FORWARD AND BACK SUBSTITUTION RESPECTIVELY IF THE SIZE OF
C     A BLOCK IS LESS THAN ICNTL(IFRLVL+MIN(10,NPIV)) AND
C     ICNTL(IFRLVL+10+MIN(10,NPIV)) RESPECTIVELY, WHERE NPIV IS THE
C     NUMBER OF PIVOTS IN THE BLOCK.
C
      INTEGER IFRLVL
      PARAMETER ( IFRLVL=5 )
C     .. Scalar Arguments ..
      INTEGER LA,LATOP,LIW,MAXFNT,N,NBLK
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION A(LA),RHS(N),W(MAXFNT)
      INTEGER IW(LIW),IW2(NBLK),ICNTL(30)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION W1,W2
      INTEGER APOS,IBLK,IFR,ILVL,IPIV,IPOS,IRHS,IROW,IST,J,J1,J2,J3,JJ,
     +        JPIV,K,K1,K2,K3,LIELL,NPIV
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,MIN
C     ..
C     .. Executable Statements ..
C APOS. RUNNING POINTER TO CURRENT PIVOT POSITION IN ARRAY A.
C IPOS. RUNNING POINTER TO BEGINNING OF BLOCK PIVOT ROW IN IW.
      APOS = 1
      IPOS = 1
      J2 = 0
      IBLK = 0
      NPIV = 0
      DO 140 IROW = 1,N
        IF (NPIV.GT.0) GO TO 90
        IBLK = IBLK + 1
        IF (IBLK.GT.NBLK) GO TO 150
        IPOS = J2 + 1
C SET UP POINTER IN PREPARATION FOR BACK SUBSTITUTION.
        IW2(IBLK) = IPOS
C ABS(LIELL) IS NUMBER OF VARIABLES (COLUMNS) IN BLOCK PIVOT ROW.
        LIELL = -IW(IPOS)
C NPIV IS NUMBER OF PIVOTS (ROWS) IN BLOCK PIVOT.
        NPIV = 1
        IF (LIELL.GT.0) GO TO 10
        LIELL = -LIELL
        IPOS = IPOS + 1
        NPIV = IW(IPOS)
   10   J1 = IPOS + 1
        J2 = IPOS + LIELL
        ILVL = MIN(NPIV,10)
        IF (LIELL.LT.ICNTL(IFRLVL+ILVL)) GO TO 90
C
C PERFORM OPERATIONS USING DIRECT ADDRESSING.
C
C LOAD APPROPRIATE COMPONENTS OF RIGHT HAND SIDES INTO ARRAY W.
        IFR = 0
        DO 20 JJ = J1,J2
          J = ABS(IW(JJ)+0)
          IFR = IFR + 1
          W(IFR) = RHS(J)
   20   CONTINUE
C JPIV IS USED AS A FLAG SO THAT IPIV IS INCREMENTED CORRECTLY AFTER
C THE USE OF A 2 BY 2 PIVOT.
        JPIV = 1
        J3 = J1
C PERFORM OPERATIONS.
        DO 70 IPIV = 1,NPIV
          JPIV = JPIV - 1
          IF (JPIV.EQ.1) GO TO 70
C JUMP IF WE HAVE A 2 BY 2 PIVOT.
          IF (IW(J3).LT.0) GO TO 40
C PERFORM FORWARD SUBSTITUTION USING 1 BY 1 PIVOT.
          JPIV = 1
          J3 = J3 + 1
          APOS = APOS + 1
          IST = IPIV + 1
          IF (LIELL.LT.IST) GO TO 70
          W1 = W(IPIV)
          K = APOS
          DO 30 J = IST,LIELL
            W(J) = W(J) + A(K)*W1
            K = K + 1
   30     CONTINUE
          APOS = APOS + LIELL - IST + 1
          GO TO 70
C PERFORM OPERATIONS WITH 2 BY 2 PIVOT.
   40     JPIV = 2
          J3 = J3 + 2
          APOS = APOS + 2
          IST = IPIV + 2
          IF (LIELL.LT.IST) GO TO 60
          W1 = W(IPIV)
          W2 = W(IPIV+1)
          K1 = APOS
          K2 = APOS + LIELL - IPIV
          DO 50 J = IST,LIELL
            W(J) = W(J) + W1*A(K1) + W2*A(K2)
            K1 = K1 + 1
            K2 = K2 + 1
   50     CONTINUE
   60     APOS = APOS + 2* (LIELL-IST+1) + 1
   70   CONTINUE
C RELOAD W BACK INTO RHS.
        IFR = 0
        DO 80 JJ = J1,J2
          J = ABS(IW(JJ)+0)
          IFR = IFR + 1
          RHS(J) = W(IFR)
   80   CONTINUE
        NPIV = 0
        GO TO 140
C
C PERFORM OPERATIONS USING INDIRECT ADDRESSING.
C
C JUMP IF WE HAVE A 2 BY 2 PIVOT.
   90   IF (IW(J1).LT.0) GO TO 110
C PERFORM FORWARD SUBSTITUTION USING 1 BY 1 PIVOT.
        NPIV = NPIV - 1
        APOS = APOS + 1
        J1 = J1 + 1
        IF (J1.GT.J2) GO TO 140
        IRHS = IW(J1-1)
        W1 = RHS(IRHS)
        K = APOS
        DO 100 J = J1,J2
          IRHS = ABS(IW(J)+0)
          RHS(IRHS) = RHS(IRHS) + A(K)*W1
          K = K + 1
  100   CONTINUE
        APOS = APOS + J2 - J1 + 1
        GO TO 140
C PERFORM OPERATIONS WITH 2 BY 2 PIVOT
  110   NPIV = NPIV - 2
        J1 = J1 + 2
        APOS = APOS + 2
        IF (J1.GT.J2) GO TO 130
        IRHS = -IW(J1-2)
        W1 = RHS(IRHS)
        IRHS = IW(J1-1)
        W2 = RHS(IRHS)
        K1 = APOS
        K3 = APOS + J2 - J1 + 2
        DO 120 J = J1,J2
          IRHS = ABS(IW(J)+0)
          RHS(IRHS) = RHS(IRHS) + W1*A(K1) + W2*A(K3)
          K1 = K1 + 1
          K3 = K3 + 1
  120   CONTINUE
  130   APOS = APOS + 2* (J2-J1+1) + 1
  140 CONTINUE
  150 LATOP = APOS - 1
      RETURN

      END
      SUBROUTINE MA27RD(N,A,LA,IW,LIW,W,MAXFNT,RHS,IW2,NBLK,LATOP,ICNTL)
C THIS SUBROUTINE PERFORMS BACKWARD ELIMINATION OPERATIONS
C     USING THE FACTORS DINVERSE AND U
C     STORED IN A/IW AFTER MA27B/BD.
C
C N      - MUST BE SET TO THE ORDER OF THE MATRIX. NOT ALTERED
C          BY MA27R/RD.
C A      - MUST BE SET TO HOLD THE REAL VALUES CORRESPONDING
C          TO THE FACTORS OF DINVERSE AND U.  THIS MUST BE
C          UNCHANGED SINCE THE PRECEDING CALL TO MA27B/BD.  NOT ALTERED
C          BY MA27R/RD.
C LA     - LENGTH OF ARRAY A. NOT ALTERED BY MA27R/RD.
C IW     - HOLDS THE INTEGER INDEXING
C          INFORMATION FOR THE MATRIX FACTORS IN A.  THIS MUST BE
C          UNCHANGED SINCE THE PRECEDING CALL TO MA27B/BD.  NOT ALTERED
C          BY MA27R/RD.
C LIW    - LENGTH OF ARRAY IW.  NOT ALTERED BY MA27R/RD.
C W      - USED
C          AS WORKSPACE BY MA27R/RD TO HOLD THE COMPONENTS OF THE RIGHT
C          HAND SIDES CORRESPONDING TO CURRENT BLOCK PIVOTAL ROWS.
C MAXFNT - INTEGER VARIABLE.  MUST BE SET TO THE LARGEST NUMBER OF
C          VARIABLES IN ANY BLOCK PIVOT ROW.  THIS VALUE WAS GIVEN
C          ON OUTPUT FROM MA27B/BD.  NOT ALTERED BY MA27R/RD.
C RHS    - ON INPUT,
C          MUST BE SET TO HOLD THE RIGHT HAND SIDE MODIFIED BY THE
C          FORWARD SUBSTITUTION OPERATIONS.  ON OUTPUT, RHS WILL HOLD
C          THE SOLUTION VECTOR.
C IW2    - ON ENTRY IW2(I) (I = 1,NBLK)
C          MUST HOLD POINTERS TO THE
C          BEGINNING OF EACH BLOCK PIVOT IN ARRAY IW, AS SET BY
C          MA27Q/QD.
C NBLK   - NUMBER OF BLOCK PIVOT ROWS. NOT ALTERED BY MA27R/RD.
C LATOP  - IT IS THE POSITION IN
C          A OF THE LAST ENTRY IN THE  FACTORS. IT MUST BE UNCHANGED
C          SINCE THE CALL TO MA27Q/QD.  IT IS NOT ALTERED BY MA27R/RD.
C ICNTL is an INTEGER array of length 30, see MA27A/AD.
C   ICNTL(IFRLVL+I) I=1,20 IS USED TO CONTROL WHETHER DIRECT OR INDIRECT
C     ACCESS IS USED BY MA27C/CD.  INDIRECT ACCESS IS EMPLOYED
C     IN FORWARD AND BACK SUBSTITUTION RESPECTIVELY IF THE SIZE OF
C     A BLOCK IS LESS THAN ICNTL(IFRLVL+MIN(10,NPIV)) AND
C     ICNTL(IFRLVL+10+MIN(10,NPIV)) RESPECTIVELY, WHERE NPIV IS THE
C     NUMBER OF PIVOTS IN THE BLOCK.
C
      INTEGER IFRLVL
      PARAMETER ( IFRLVL=5 )
C
C     .. Scalar Arguments ..
      INTEGER LA,LATOP,LIW,MAXFNT,N,NBLK
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION A(LA),RHS(N),W(MAXFNT)
      INTEGER IW(LIW),IW2(NBLK),ICNTL(30)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION W1,W2
      INTEGER APOS,APOS2,I1RHS,I2RHS,IBLK,IFR,IIPIV,IIRHS,ILVL,IPIV,
     +        IPOS,IRHS,IST,J,J1,J2,JJ,JJ1,JJ2,JPIV,JPOS,K,LIELL,LOOP,
     +        NPIV
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,MIN
C     ..
C     .. Executable Statements ..
C APOS. RUNNING POINTER TO CURRENT PIVOT POSITION IN ARRAY A.
C IPOS. RUNNING POINTER TO BEGINNING OF CURRENT BLOCK PIVOT ROW.
      APOS = LATOP + 1
      NPIV = 0
      IBLK = NBLK + 1
C RUN THROUGH BLOCK PIVOT ROWS IN THE REVERSE ORDER.
      DO 180 LOOP = 1,N
        IF (NPIV.GT.0) GO TO 110
        IBLK = IBLK - 1
        IF (IBLK.LT.1) GO TO 190
        IPOS = IW2(IBLK)
C ABS(LIELL) IS NUMBER OF VARIABLES (COLUMNS) IN BLOCK PIVOT ROW.
        LIELL = -IW(IPOS)
C NPIV IS NUMBER OF PIVOTS (ROWS) IN BLOCK PIVOT.
        NPIV = 1
        IF (LIELL.GT.0) GO TO 10
        LIELL = -LIELL
        IPOS = IPOS + 1
        NPIV = IW(IPOS)
   10   JPOS = IPOS + NPIV
        J2 = IPOS + LIELL
        ILVL = MIN(10,NPIV) + 10
        IF (LIELL.LT.ICNTL(IFRLVL+ILVL)) GO TO 110
C
C PERFORM OPERATIONS USING DIRECT ADDRESSING.
C
        J1 = IPOS + 1
C LOAD APPROPRIATE COMPONENTS OF RHS INTO W.
        IFR = 0
        DO 20 JJ = J1,J2
          J = ABS(IW(JJ)+0)
          IFR = IFR + 1
          W(IFR) = RHS(J)
   20   CONTINUE
C JPIV IS USED AS A FLAG SO THAT IPIV IS INCREMENTED CORRECTLY AFTER
C     THE USE OF A 2 BY 2 PIVOT.
        JPIV = 1
C PERFORM ELIMINATIONS.
        DO 90 IIPIV = 1,NPIV
          JPIV = JPIV - 1
          IF (JPIV.EQ.1) GO TO 90
          IPIV = NPIV - IIPIV + 1
          IF (IPIV.EQ.1) GO TO 30
C JUMP IF WE HAVE A 2 BY 2 PIVOT.
          IF (IW(JPOS-1).LT.0) GO TO 60
C PERFORM BACK-SUBSTITUTION USING 1 BY 1 PIVOT.
   30     JPIV = 1
          APOS = APOS - (LIELL+1-IPIV)
          IST = IPIV + 1
          W1 = W(IPIV)*A(APOS)
          IF (LIELL.LT.IST) GO TO 50
          JJ1 = APOS + 1
          DO 40 J = IST,LIELL
            W1 = W1 + A(JJ1)*W(J)
            JJ1 = JJ1 + 1
   40     CONTINUE
   50     W(IPIV) = W1
          JPOS = JPOS - 1
          GO TO 90
C PERFORM BACK-SUBSTITUTION OPERATIONS WITH 2 BY 2 PIVOT
   60     JPIV = 2
          APOS2 = APOS - (LIELL+1-IPIV)
          APOS = APOS2 - (LIELL+2-IPIV)
          IST = IPIV + 1
          W1 = W(IPIV-1)*A(APOS) + W(IPIV)*A(APOS+1)
          W2 = W(IPIV-1)*A(APOS+1) + W(IPIV)*A(APOS2)
          IF (LIELL.LT.IST) GO TO 80
          JJ1 = APOS + 2
          JJ2 = APOS2 + 1
          DO 70 J = IST,LIELL
            W1 = W1 + W(J)*A(JJ1)
            W2 = W2 + W(J)*A(JJ2)
            JJ1 = JJ1 + 1
            JJ2 = JJ2 + 1
   70     CONTINUE
   80     W(IPIV-1) = W1
          W(IPIV) = W2
          JPOS = JPOS - 2
   90   CONTINUE
C RELOAD WORKING VECTOR INTO SOLUTION VECTOR.
        IFR = 0
        DO 100 JJ = J1,J2
          J = ABS(IW(JJ)+0)
          IFR = IFR + 1
          RHS(J) = W(IFR)
  100   CONTINUE
        NPIV = 0
        GO TO 180
C
C PERFORM OPERATIONS USING INDIRECT ADDRESSING.
C
  110   IF (NPIV.EQ.1) GO TO 120
C JUMP IF WE HAVE A 2 BY 2 PIVOT.
        IF (IW(JPOS-1).LT.0) GO TO 150
C PERFORM BACK-SUBSTITUTION USING 1 BY 1 PIVOT.
  120   NPIV = NPIV - 1
        APOS = APOS - (J2-JPOS+1)
        IIRHS = IW(JPOS)
        W1 = RHS(IIRHS)*A(APOS)
        J1 = JPOS + 1
        IF (J1.GT.J2) GO TO 140
        K = APOS + 1
        DO 130 J = J1,J2
          IRHS = ABS(IW(J)+0)
          W1 = W1 + A(K)*RHS(IRHS)
          K = K + 1
  130   CONTINUE
  140   RHS(IIRHS) = W1
        JPOS = JPOS - 1
        GO TO 180
C PERFORM OPERATIONS WITH 2 BY 2 PIVOT
  150   NPIV = NPIV - 2
        APOS2 = APOS - (J2-JPOS+1)
        APOS = APOS2 - (J2-JPOS+2)
        I1RHS = -IW(JPOS-1)
        I2RHS = IW(JPOS)
        W1 = RHS(I1RHS)*A(APOS) + RHS(I2RHS)*A(APOS+1)
        W2 = RHS(I1RHS)*A(APOS+1) + RHS(I2RHS)*A(APOS2)
        J1 = JPOS + 1
        IF (J1.GT.J2) GO TO 170
        JJ1 = APOS + 2
        JJ2 = APOS2 + 1
        DO 160 J = J1,J2
          IRHS = ABS(IW(J)+0)
          W1 = W1 + RHS(IRHS)*A(JJ1)
          W2 = W2 + RHS(IRHS)*A(JJ2)
          JJ1 = JJ1 + 1
          JJ2 = JJ2 + 1
  160   CONTINUE
  170   RHS(I1RHS) = W1
        RHS(I2RHS) = W2
        JPOS = JPOS - 2
  180 CONTINUE
  190 RETURN

      END