C
c######date   01 jan 1984     copyright ukaea, harwell.
c######alias ma28ad ma28bd ma28cd
c###### calls   ma30    mc20    mc22    mc23    mc24
      SUBROUTINE MA28AD(N, NZ, A, LICN, IRN, LIRN, ICN, U, IKEEP, IW, W,
     * IFLAG)
c this subroutine performs the lu factorization of a.
c
c the parameters are as follows.....
c n     order of matrix  not altered by subroutine.
c nz    number of non-zeros in input matrix  not altered by subroutine.
c a is a  real array  length licn.  holds non-zeros of matrix on entry
c     and non-zeros of factors on exit.  reordered by mc20a/ad and
c     mc23a/ad and altered by ma30a/ad.
c licn  integer  length of arrays a and icn.  not altered by subroutine.
c irn   integer array of length lirn.  holds row indices on input.
c     used as workspace by ma30a/ad to hold column orientation of
c     matrix.
c lirn  integer  length of array irn. not altered by the subroutine.
c icn   integer array of length licn.  holds column indices on entry
c     and column indices of decomposed matrix on exit. reordered by
c     mc20a/ad and mc23a/ad and altered by ma30a/ad.
c u     real variable  set by user to control bias towards numeric or
c     sparsity pivoting.  u=1.0 gives partial pivoting while u=0. does
c     not check multipliers at all.  values of u greater than one are
c     treated as one while negative values are treated as zero.  not
c     altered by subroutine.
c ikeep  integer array of length 5*n  used as workspace by ma28a/ad
c     (see later comments).  it is not required to be set on entry
c     and, on exit, it contains information about the decomposition.
c     it should be preserved between this call and subsequent calls
c     to ma28b/bd or ma28c/cd.
c     ikeep(i,1),i=1,n  holds the total length of the part of row i
c     in the diagonal block.
c     row ikeep(i,2),i=1,n  of the input matrix is the ith row in
c     pivot order.
c     column ikeep(i,3),i=1,n  of the input matrix is the ith column
c     in pivot order.
c     ikeep(i,4),i=1,n  holds the length of the part of row i in
c     the l part of the l/u decomposition.
c     ikeep(i,5),i=1,n  holds the length of the part of row i in the
c     off-diagonal blocks.  if there is only one diagonal block,
c     ikeep(1,5) will be set to -1.
c iw    integer array of length 8*n.  if the option nsrch.le.n is
c     used, then the length of array iw can be reduced to 7*n.
c w real array  length n.  used by mc24a/ad both as workspace and to
c     return growth estimate in w(1).  the use of this array by ma28a/ad
c     is thus optional depending on common block logical variable grow.
c iflag  integer variable  used as error flag by routine.  a positive
c     or zero value on exit indicates success.  possible negative
c     values are -1 through -14.
c
      INTEGER N, NZ, LICN, LIRN, IFLAG
      INTEGER IRN(LIRN), ICN(LICN), IKEEP(N,5), IW(N,8)
      DOUBLE PRECISION A(LICN), U, W(N)
c
c common and private variables.
c     common block ma28f/fd is used merely
c     to communicate with common block ma30f/fd  so that the user
c     need not declare this common block in his main program.
c the common block variables are as follows ...
c lp,mp  integer  default value 6 (line printer).  unit number
c     for error messages and duplicate element warning resp.
c nlp,mlp  integer  unit number for messages from ma30a/ad and
c     mc23a/ad resp.  set by ma28a/ad to value of lp.
c lblock  logical  default value true.  if true mc23a/ad is used
c     to first permute the matrix to block lower triangular form.
c grow    logical  default value true.  if true then an estimate
c     of the increase in size of matrix elements during l/u
c     decomposition is given by mc24a/ad.
c eps,rmin,resid  real/double precision variables not referenced
c     by ma28a/ad.
c irncp,icncp  integer  set to number of compresses on arrays irn and
c     icn/a respectively.
c minirn,minicn  integer  minimum length of arrays irn and icn/a
c     respectively, for success on future runs.
c irank  integer   estimated rank of matrix.
c mirncp,micncp,mirank,mirn,micn integer variables.  used to
c     communicate between ma30f/fd and ma28f/fd values of abovenamed
c     variables with somewhat similar names.
c abort1,abort2  logical variables with default value true.  if false
c     then decomposition will be performed even if the matrix is
c     structurally or numerically singular respectively.
c aborta,abortb  logical variables used to communicate values of
c     abort1 and abort2 to ma30a/ad.
c abort  logical  used to communicate value of abort1 to mc23a/ad.
c abort3  logical variable not referenced by ma28a/ad.
c idisp   integer array  length 2.  used to communicate information
c     on decomposition between this call to ma28a/ad and subsequent
c     calls to ma28b/bd and ma28c/cd.  on exit, idisp(1) and
c     idisp(2) indicate position in arrays a and icn of the
c     first and last elements in the l/u decomposition of the
c     diagonal blocks, respectively.
c numnz  integer  structural rank of matrix.
c num    integer  number of diagonal blocks.
c large  integer  size of largest diagonal block.
c
c see block data for further comments on common block variables.
c see code for comments on private variables.
c
      DOUBLE PRECISION TOL, THEMAX, BIG, DXMAX, ERRMAX, DRES, CGCE,
     * TOL1, BIG1, UPRIV, RMIN, EPS, RESID, ZERO
      INTEGER IDISP(2)
      LOGICAL GROW, LBLOCK, ABORT, ABORT1, ABORT2, ABORT3, ABORTA,
     * ABORTB, LBIG, LBIG1
      COMMON /MA28ED/ LP, MP, LBLOCK, GROW
      COMMON /MA28FD/ EPS, RMIN, RESID, IRNCP, ICNCP, MINIRN, MINICN,
     * IRANK, ABORT1, ABORT2
      COMMON /MA28GD/ IDISP
      COMMON /MA28HD/ TOL, THEMAX, BIG, DXMAX, ERRMAX, DRES, CGCE,
     * NDROP, MAXIT, NOITER, NSRCH, ISTART, LBIG
      COMMON /MA30ID/ TOL1, BIG1, NDROP1, NSRCH1, LBIG1
      COMMON /MA30ED/ NLP, ABORTA, ABORTB, ABORT3
      COMMON /MA30FD/ MIRNCP, MICNCP, MIRANK, MIRN, MICN
      COMMON /MC23BD/ MLP, NUMNZ, NUM, LARGE, ABORT
      COMMON /LPIVOT/ LPIV(10),LNPIV(10),MAPIV,MANPIV,IAVPIV,
     *                IANPIV,KOUNTL
c
c some  initialization and transfer of information between
c     common blocks (see earlier comments).
      DATA ZERO /0.0D0/
      IFLAG = 0
      ABORTA = ABORT1
      ABORTB = ABORT2
      ABORT = ABORT1
      MLP = LP
      NLP = LP
      TOL1 = TOL
      LBIG1 = LBIG
      NSRCH1 = NSRCH
c upriv private copy of u is used in case it is outside
c     range  zero to one  and  is thus altered by ma30a/ad.
      UPRIV = U
c simple data check on input variables and array dimensions.
      IF (N.GT.0) GO TO 10
      IFLAG = -8
      IF (LP.NE.0) WRITE (LP,99999) N
      GO TO 210
   10 IF (NZ.GT.0) GO TO 20
      IFLAG = -9
      IF (LP.NE.0) WRITE (LP,99998) NZ
      GO TO 210
   20 IF (LICN.GE.NZ) GO TO 30
      IFLAG = -10
      IF (LP.NE.0) WRITE (LP,99997) LICN
      GO TO 210
   30 IF (LIRN.GE.NZ) GO TO 40
      IFLAG = -11
      IF (LP.NE.0) WRITE (LP,99996) LIRN
      GO TO 210
c
c data check to see if all indices lie between 1 and n.
   40 DO 50 I=1,NZ
        IF (IRN(I).GT.0 .AND. IRN(I).LE.N .AND. ICN(I).GT.0 .AND.
     *   ICN(I).LE.N) GO TO 50
        IF (IFLAG.EQ.0 .AND. LP.NE.0) WRITE (LP,99995)
        IFLAG = -12
        IF (LP.NE.0) WRITE (LP,99994) I, A(I), IRN(I), ICN(I)
   50 CONTINUE
      IF (IFLAG.LT.0) GO TO 220
c
c sort matrix into row order.
      CALL MC20AD(N, NZ, A, ICN, IW, IRN, 0)
c part of ikeep is used here as a work-array.  ikeep(i,2) is
c     the last row to have a non-zero in column i.  ikeep(i,3)
c     is the off-set of column i from the start of the row.
      DO 60 I=1,N
        IKEEP(I,2) = 0
        IKEEP(I,1) = 0
   60 CONTINUE
c
c check for duplicate elements .. summing any such entries and
c     printing a warning message on unit mp.
c move is equal to the number of duplicate elements found.
      MOVE = 0
c the loop also calculates the largest element in the matrix, themax.
      THEMAX = ZERO
c j1 is position in arrays of first non-zero in row.
      J1 = IW(1,1)
      DO 130 I=1,N
        IEND = NZ + 1
        IF (I.NE.N) IEND = IW(I+1,1)
        LENGTH = IEND - J1
        IF (LENGTH.EQ.0) GO TO 130
        J2 = IEND - 1
        NEWJ1 = J1 - MOVE
        DO 120 JJ=J1,J2
          J = ICN(JJ)
          THEMAX = DMAX1(THEMAX,DABS(A(JJ)))
          IF (IKEEP(J,2).EQ.I) GO TO 110
c first time column has ocurred in current row.
          IKEEP(J,2) = I
          IKEEP(J,3) = JJ - MOVE - NEWJ1
          IF (MOVE.EQ.0) GO TO 120
c shift necessary because of  previous duplicate element.
          NEWPOS = JJ - MOVE
          A(NEWPOS) = A(JJ)
          ICN(NEWPOS) = ICN(JJ)
          GO TO 120
c duplicate element.
  110     MOVE = MOVE + 1
          LENGTH = LENGTH - 1
          JAY = IKEEP(J,3) + NEWJ1
          IF (MP.NE.0) WRITE (MP,99993) I, J, A(JJ)
          A(JAY) = A(JAY) + A(JJ)
          THEMAX = DMAX1(THEMAX,DABS(A(JAY)))
  120   CONTINUE
        IKEEP(I,1) = LENGTH
        J1 = IEND
  130 CONTINUE
c
c knum is actual number of non-zeros in matrix with any multiple
c     entries counted only once.
      KNUM = NZ - MOVE
      IF (.NOT.LBLOCK) GO TO 140
c
c perform block triangularisation.
      CALL MC23AD(N, ICN, A, LICN, IKEEP, IDISP, IKEEP(1,2),
     *IKEEP(1,3), IKEEP(1,5), IW(1,3), IW)
      IF (IDISP(1).GT.0) GO TO 170
      IFLAG = -7
      IF (IDISP(1).EQ.-1) IFLAG = -1
      IF (LP.NE.0) WRITE (LP,99992)
      GO TO 210
c
c block triangularization not requested.
c move structure to end of data arrays in preparation for
c     ma30a/ad.
c also set lenoff(1) to -1 and set permutation arrays.
  140 DO 150 I=1,KNUM
        II = KNUM - I + 1
        NEWPOS = LICN - I + 1
        ICN(NEWPOS) = ICN(II)
        A(NEWPOS) = A(II)
  150 CONTINUE
      IDISP(1) = 1
      IDISP(2) = LICN - KNUM + 1
      DO 160 I=1,N
        IKEEP(I,2) = I
        IKEEP(I,3) = I
  160 CONTINUE
      IKEEP(1,5) = -1
  170 IF (LBIG) BIG1 = THEMAX
      IF (NSRCH.LE.N) GO TO 180
c
c perform l/u decomosition on diagonal blocks.
      CALL MA30AD(N, ICN, A, LICN, IKEEP, IKEEP(1,4), IDISP,
     *IKEEP(1,2), IKEEP(1,3), IRN, LIRN, IW(1,2), IW(1,3), IW(1,4),
     *IW(1,5), IW(1,6), IW(1,7), IW(1,8), IW, UPRIV, IFLAG)
      GO TO 190
c this call if used if nsrch has been set less than or equal n.
c     in this case, two integer work arrays of length can be saved.
  180 CALL MA30AD(N, ICN, A, LICN, IKEEP, IKEEP(1,4), IDISP,
     * IKEEP(1,2), IKEEP(1,3), IRN, LIRN, IW(1,2), IW(1,3), IW(1,4),
     * IW(1,5), IW, IW, IW(1,6), IW, UPRIV, IFLAG)
c
c transfer common block information.
  190 MINIRN = MAX0(MIRN,NZ)
      MINICN = MAX0(MICN,NZ)
      IRNCP = MIRNCP
      ICNCP = MICNCP
      IRANK = MIRANK
      NDROP = NDROP1
      IF (LBIG) BIG = BIG1
      IF (IFLAG.GE.0) GO TO 200
      IF (LP.NE.0) WRITE (LP,99991)
      GO TO 210
c
c reorder off-diagonal blocks according to pivot permutation.
  200 I1 = IDISP(1) - 1
      IF (I1.NE.0) CALL MC22AD(N, ICN, A, I1, IKEEP(1,5), IKEEP(1,2),
     * IKEEP(1,3), IW, IRN)
      I1 = IDISP(1)
      IEND = LICN - I1 + 1
c
c optionally calculate element growth estimate.
      IF (GROW) CALL MC24AD(N, ICN, A(I1), IEND, IKEEP, IKEEP(1,4), W)
c increment growth estimate by original maximum element.
      IF (GROW) W(1) = W(1) + THEMAX
      IF (GROW .AND. N.GT.1) W(2) = THEMAX
c set flag if the only error is due to duplicate elements.
      IF (IFLAG.GE.0 .AND. MOVE.NE.0) IFLAG = -14
      GO TO 220
  210 IF (LP.NE.0) WRITE (LP,99990)
  220 RETURN
99999 FORMAT (36X, 17HN OUT OF RANGE = , I10)
99998 FORMAT (36X, 18HNZ NON POSITIVE = , I10)
99997 FORMAT (36X, 17HLICN TOO SMALL = , I10)
99996 FORMAT (36X, 17HLIRN TOO SMALL = , I10)
99995 FORMAT (54H ERROR RETURN FROM MA28A/AD BECAUSE INDICES FOUND OUT ,
     * 8HOF RANGE)
99994 FORMAT (1X, I6, 22HTH ELEMENT WITH VALUE , 1PD22.14, 9H IS OUT O,
     * 21HF RANGE WITH INDICES , I8, 2H ,, I8)
99993 FORMAT (31H DUPLICATE ELEMENT IN POSITION , I8, 2H ,, I8,
     * 12H WITH VALUE , 1PD22.14)
99992 FORMAT (36X, 26HERROR RETURN FROM MC23A/AD)
99991 FORMAT (36X, 26HERROR RETURN FROM MA30A/AD)
99990 FORMAT (36H+ERROR RETURN FROM MA28A/AD BECAUSE )
      END
c
c==============================================================
c
      SUBROUTINE MA28BD(N, NZ, A, LICN, IVECT, JVECT, ICN, IKEEP, IW, W,
     * IFLAG)
c this subroutine factorizes a matrix of a similar sparsity
c     pattern to that previously factorized by ma28a/ad.
c the parameters are as follows ...
c n      integer  order of matrix  not altered by subroutine.
c nz     integer  number of non-zeros in input matrix  not altered
c     by subroutine.
c a      real/double precision array  length licn.  holds non-zeros of
c     matrix on entry and non-zeros of factors on exit.  reordered by
c     ma28d/dd and altered by subroutine ma30b/bd.
c licn   integer  length of arrays a and icn.  not altered by
c     subroutine.
c ivect,jvect  integer arrays of length nz.  hold row and column
c     indices of non-zeros respectively.  not altered by subroutine.
c icn    integer array of length licn.  same array as output from
c     ma28a/ad.  unchanged by ma28b/bd.
c ikeep  integer array of length 5*n.  same array as output from
c     ma28a/ad.  unchanged by ma28b/bd.
c iw     integer array  length 5*n.  used as workspace by ma28d/dd and
c     ma30b/bd.
c w      real/double precision array  length n.  used as workspace
c     by ma28d/dd,ma30b/bd and (optionally) mc24a/ad.
c iflag  integer  used as error flag with positive or zero value
c     indicating success.
c
      INTEGER N, NZ, LICN, IW(N,5), IFLAG
      INTEGER IKEEP(N,5), IVECT(NZ), JVECT(NZ), ICN(LICN)
      DOUBLE PRECISION A(LICN), W(N)
c
c private and common variables.
c unless otherwise stated common block variables are as in ma28a/ad.
c     those variables referenced by ma28b/bd are mentioned below.
c lp,mp  integers  used as in ma28a/ad as unit number for error and
c     warning messages, respectively.
c nlp    integer variable used to give value of lp to ma30e/ed.
c eps    real/double precision  ma30b/bd will output a positive value
c     for iflag if any modulus of the ratio of pivot element to the
c     largest element in its row (u part only) is less than eps (unless
c     eps is greater than 1.0 when no action takes place).
c rmin   real/double precision  variable equal to the value of this
c     minimum ratio in cases where eps is less than or equal to 1.0.
c meps,mrmin  real/double precision variables used by the subroutine
c     to communicate between common blocks ma28f/fd and ma30g/gd.
c idisp  integer array  length 2  the same as that used by ma28a/ad.
c     it is unchanged by ma28b/bd.
c
c see block data or ma28a/ad for further comments on variables
c     in common.
c see code for comments on private variables.
c
      LOGICAL GROW, LBLOCK, ABORTA, ABORTB, ABORT1, ABORT2, ABORT3,
     * LBIG, LBIG1
      INTEGER IDISP(2)
      DOUBLE PRECISION EPS, MEPS, RMIN, MRMIN, RESID, TOL,
     * THEMAX, BIG, DXMAX, ERRMAX, DRES, CGCE, TOL1, BIG1
c
      COMMON /MA28ED/ MP, LP, LBLOCK, GROW
      COMMON /MA28FD/ EPS, RMIN, RESID, IRNCP, ICNCP, MINIRN, MINICN,
     * IRANK, ABORT1, ABORT2
      COMMON /MA28GD/ IDISP
      COMMON /MA28HD/ TOL, THEMAX, BIG, DXMAX, ERRMAX, DRES, CGCE,
     * NDROP, MAXIT, NOITER, NSRCH, ISTART, LBIG
      COMMON /MA30ED/ NLP, ABORTA, ABORTB, ABORT3
      COMMON /MA30GD/ MEPS, MRMIN
      COMMON /MA30ID/ TOL1, BIG1, NDROP1, NSRCH1, LBIG1
c
c check to see if elements were dropped in previous ma28a/ad call.
      IF (NDROP.EQ.0) GO TO 10
      IFLAG = -15
      WRITE (6,99999) IFLAG, NDROP
      GO TO 70
   10 IFLAG = 0
      MEPS = EPS
      NLP = LP
c simple data check on variables.
      IF (N.GT.0) GO TO 20
      IFLAG = -11
      IF (LP.NE.0) WRITE (LP,99998) N
      GO TO 60
   20 IF (NZ.GT.0) GO TO 30
      IFLAG = -10
      IF (LP.NE.0) WRITE (LP,99997) NZ
      GO TO 60
   30 IF (LICN.GE.NZ) GO TO 40
      IFLAG = -9
      IF (LP.NE.0) WRITE (LP,99996) LICN
      GO TO 60
c
   40 CALL MA28DD(N, A, LICN, IVECT, JVECT, NZ, ICN, IKEEP, IKEEP(1,4),
     * IKEEP(1,5), IKEEP(1,2), IKEEP(1,3), IW(1,3), IW, W(1), IFLAG)
c themax is largest element in matrix.
      THEMAX = W(1)
      IF (LBIG) BIG1 = THEMAX
c idup equals one if there were duplicate elements, zero otherwise.
      IDUP = 0
      IF (IFLAG.EQ.(N+1)) IDUP = 1
      IF (IFLAG.LT.0) GO TO 60
c
c perform row-gauss elimination on the structure received from ma28d/dd
      CALL MA30BD(N, ICN, A, LICN, IKEEP, IKEEP(1,4), IDISP,
     * IKEEP(1,2), IKEEP(1,3), W, IW, IFLAG)
c
c transfer common block information.
      IF (LBIG) BIG1 = BIG
      RMIN = MRMIN
      IF (IFLAG.GE.0) GO TO 50
      IFLAG = -2
      IF (LP.NE.0) WRITE (LP,99995)
      GO TO 60
c
c optionally calculate the growth parameter.
   50 I1 = IDISP(1)
      IEND = LICN - I1 + 1
      IF (GROW) CALL MC24AD(N, ICN, A(I1), IEND, IKEEP, IKEEP(1,4), W)
c increment estimate by largest element in input matrix.
      IF (GROW) W(1) = W(1) + THEMAX
      IF (GROW .AND. N.GT.1) W(2) = THEMAX
c set flag if the only error is due to duplicate elements.
      IF (IDUP.EQ.1 .AND. IFLAG.GE.0) IFLAG = -14
      GO TO 70
   60 IF (LP.NE.0) WRITE (LP,99994)
   70 RETURN
99999 FORMAT (39H ERROR RETURN FROM MA28B/BD WITH IFLAG=, I4/I7, 4H ENT,
     * 39HRIES DROPPED FROM STRUCTURE BY MA28A/AD)
99998 FORMAT (36X, 17HN OUT OF RANGE = , I10)
99997 FORMAT (36X, 18HNZ NON POSITIVE = , I10)
99996 FORMAT (36X, 17HLICN TOO SMALL = , I10)
99995 FORMAT (36X, 26HERROR RETURN FROM MA30B/BD)
99994 FORMAT (36H+ERROR RETURN FROM MA28B/BD BECAUSE )
      END
c
c============================================================
c
      SUBROUTINE MA28CD(N, A, LICN, ICN, IKEEP, RHS, W, MTYPE)
c
c this subroutine uses the factors from ma28a/ad or ma28b/bd to
c     solve a system of equations without iterative refinement.
c the parameters are ...
c n   integer  order of matrix  not altered by subroutine.
c a   real/double precision array  length licn.  the same array as
c     was used in the most recent call to ma28a/ad or ma28b/bd.
c licn  integer  length of arrays a and icn.  not altered by
c     subroutine.
c icn    integer array of length licn.  same array as output from
c     ma28a/ad.  unchanged by ma28c/cd.
c ikeep  integer array of length 5*n.  same array as output from
c     ma28a/ad.  unchanged by ma28c/cd.
c rhs    real/double precision array  length n.  on entry, it holds the
c     right hand side.  on exit, the solution vector.
c w      real/double precision array  length n. used as workspace by
c     ma30c/cd.
c mtype  integer  used to tell ma30c/cd to solve the direct equation
c     (mtype=1) or its transpose (mtype.ne.1).
c
      DOUBLE PRECISION A(LICN), RHS(N), W(N), RESID, MRESID, EPS, RMIN
      INTEGER IDISP(2)
      INTEGER ICN(LICN), IKEEP(N,5)
      LOGICAL ABORT1, ABORT2
c common block variables.
c unless otherwise stated common block variables are as in ma28a/ad.
c     those variables referenced by ma28c/cd are mentioned below.
c resid  real/double precision  variable returns maximum residual of
c     equations where pivot was zero.
c mresid  real/double precision variable used by ma28c/cd to
c     communicate between ma28f/fd and ma30h/hd.
c idisp  integer array  length 2  the same as that used by ma28a/ad.
c     it is unchanged by ma28b/bd.
c
c further information on common block variables can be found in block
c     data or ma28a/ad.
      COMMON /MA28FD/ EPS, RMIN, RESID, IRNCP, ICNCP, MINIRN, MINICN,
     * IRANK, ABORT1, ABORT2
      COMMON /MA28GD/ IDISP
      COMMON /MA30HD/ MRESID
c
c this call performs the solution of the set of equations.
      CALL MA30CD(N, ICN, A, LICN, IKEEP, IKEEP(1,4), IKEEP(1,5),
     * IDISP, IKEEP(1,2), IKEEP(1,3), RHS, W, MTYPE)
c transfer common block information.
      RESID = MRESID
      RETURN
      END
      SUBROUTINE MA28DD(N, A, LICN, IVECT, JVECT, NZ, ICN, LENR, LENRL,
     * LENOFF, IP, IQ, IW1, IW, W1, IFLAG)
c this subroutine need never be called by the user directly.
c     it sorts the user's matrix into the structure of the decomposed
c     form and checks for the presence of duplicate entries or
c     non-zeros lying outside the sparsity pattern of the decomposition
c     it also calculates the largest element in the input matrix.
      DOUBLE PRECISION A(LICN), ZERO, W1, AA
      INTEGER IW(N,2), IDISP(2)
      INTEGER ICN(LICN), IVECT(NZ), JVECT(NZ), IP(N), IQ(N),
     * LENR(N), IW1(N,3), LENRL(N), LENOFF(N)
      LOGICAL LBLOCK, GROW, BLOCKL
      COMMON /MA28ED/ LP, MP, LBLOCK, GROW
      COMMON /MA28GD/ IDISP
      DATA ZERO /0.0D0/
      BLOCKL = LENOFF(1).GE.0
c iw1(i,3)  is set to the block in which row i lies and the
c     inverse permutations to ip and iq are set in iw1(.,1) and
c     iw1(.,2) resp.
c pointers to beginning of the part of row i in diagonal and
c   off-diagonal blocks are set in iw(i,2) and iw(i,1) resp.
      IBLOCK = 1
      IW(1,1) = 1
      IW(1,2) = IDISP(1)
      DO 10 I=1,N
        IW1(I,3) = IBLOCK
        IF (IP(I).LT.0) IBLOCK = IBLOCK + 1
        II = IABS(IP(I)+0)
        IW1(II,1) = I
        JJ = IQ(I)
        JJ = IABS(JJ)
        IW1(JJ,2) = I
        IF (I.EQ.1) GO TO 10
        IF (BLOCKL) IW(I,1) = IW(I-1,1) + LENOFF(I-1)
        IW(I,2) = IW(I-1,2) + LENR(I-1)
   10 CONTINUE
c place each non-zero in turn into its correct location
c    in the a/icn array.
      IDISP2 = IDISP(2)
      DO 170 I=1,NZ
c necessary to avoid reference to unassigned element of icn.
        IF (I.GT.IDISP2) GO TO 20
        IF (ICN(I).LT.0) GO TO 170
   20   IOLD = IVECT(I)
        JOLD = JVECT(I)
        AA = A(I)
c this is a dummy loop for following a chain of interchanges.
c   it will be executed nz times in total.
        DO 140 IDUMMY=1,NZ
c perform some validity checks on iold and jold.
          IF (IOLD.LE.N .AND. IOLD.GT.0 .AND. JOLD.LE.N .AND.
     *     JOLD.GT.0) GO TO 30
          IF (LP.NE.0) WRITE (LP,99999) I, A(I), IOLD, JOLD
          IFLAG = -12
          GO TO 180
   30     INEW = IW1(IOLD,1)
          JNEW = IW1(JOLD,2)
c are we in a valid block and is it diagonal or off-diagonal?
          IF (IW1(INEW,3)-IW1(JNEW,3)) 40, 60, 50
   40     IFLAG = -13
          IF (LP.NE.0) WRITE (LP,99998) IOLD, JOLD
          GO TO 180
   50     J1 = IW(INEW,1)
          J2 = J1 + LENOFF(INEW) - 1
          GO TO 110
c element is in diagonal block.
   60     J1 = IW(INEW,2)
          IF (INEW.GT.JNEW) GO TO 70
          J2 = J1 + LENR(INEW) - 1
          J1 = J1 + LENRL(INEW)
          GO TO 110
   70     J2 = J1 + LENRL(INEW)
c binary search of ordered list  .. element in l part of row.
          DO 100 JDUMMY=1,N
            MIDPT = (J1+J2)/2
            JCOMP = IABS(ICN(MIDPT)+0)
            IF (JNEW-JCOMP) 80, 130, 90
   80       J2 = MIDPT
            GO TO 100
   90       J1 = MIDPT
  100     CONTINUE
          IFLAG = -13
          IF (LP.NE.0) WRITE (LP,99997) IOLD, JOLD
          GO TO 180
c linear search ... element in l part of row or off-diagonal blocks.
  110     DO 120 MIDPT=J1,J2
            IF (IABS(ICN(MIDPT)+0).EQ.JNEW) GO TO 130
  120     CONTINUE
          IFLAG = -13
          IF (LP.NE.0) WRITE (LP,99997) IOLD, JOLD
          GO TO 180
c equivalent element of icn is in position midpt.
  130     IF (ICN(MIDPT).LT.0) GO TO 160
          IF (MIDPT.GT.NZ .OR. MIDPT.LE.I) GO TO 150
          W1 = A(MIDPT)
          A(MIDPT) = AA
          AA = W1
          IOLD = IVECT(MIDPT)
          JOLD = JVECT(MIDPT)
          ICN(MIDPT) = -ICN(MIDPT)
  140   CONTINUE
  150   A(MIDPT) = AA
        ICN(MIDPT) = -ICN(MIDPT)
        GO TO 170
  160   A(MIDPT) = A(MIDPT) + AA
c set flag for duplicate elements.
        IFLAG = N + 1
  170 CONTINUE
c reset icn array  and zero elements in l/u but not in a.
c also calculate maximum element of a.
  180 W1 = ZERO
      DO 200 I=1,IDISP2
        IF (ICN(I).LT.0) GO TO 190
        A(I) = ZERO
        GO TO 200
  190   ICN(I) = -ICN(I)
        W1 = DMAX1(W1,DABS(A(I)))
  200 CONTINUE
      RETURN
99999 FORMAT (9H ELEMENT , I6, 12H WITH VALUE , 1PD22.14, 10H HAS INDIC,
     * 3HES , I8, 2H ,, I8/36X, 20HINDICES OUT OF RANGE)
99998 FORMAT (36X, 8HNON-ZERO, I7, 2H ,, I6, 23H IN ZERO OFF-DIAGONAL B,
     * 4HLOCK)
99997 FORMAT (36X, 8H ELEMENT, I6, 2H ,, I6, 23H WAS NOT IN L/U PATTERN)
      END
c######date   01 jan 1984     copyright ukaea, harwell.
c######alias ma30ad
c
c===================================================================
c
      SUBROUTINE MA30AD(NN, ICN, A, LICN, LENR, LENRL, IDISP, IP, IQ,
     * IRN, LIRN, LENC, IFIRST, LASTR, NEXTR, LASTC, NEXTC, IPTR, IPC,
     * U, IFLAG)
c if  the user requires a more convenient data interface then the ma28
c     package should be used.  the ma28 subroutines call the ma30
c     subroutines after checking the user's input data and optionally
c     using mc23a/ad to permute the matrix to block triangular form.
c this package of subroutines (ma30a/ad, ma30b/bd, ma30c/cd and
c     ma30d/dd) performs operations pertinent to the solution of a
c     general sparse n by n system of linear equations (i.e. solve
c     ax=b). structually singular matrices are permitted including
c     those with row or columns consisting entirely of zeros (i.e.
c     including rectangular matrices).  it is assumed that the
c     non-zeros of the matrix a do not differ widely in size.  if
c     necessary a prior call of the scaling subroutine mc19a/ad may be
c     made.
c a discussion of the design of these subroutines is given by duff and
c     reid (acm trans math software 5 pp 18-35,1979 (css 48)) while
c     fuller details of the implementation are given in duff (harwell
c     report aere-r 8730,1977).  the additional pivoting option in
c     ma30a/ad and the use of drop tolerances (see common block
c     ma30i/id) were added to the package after joint work with reid,
c     schaumburg, wasniewski and zlatev (duff, reid, schaumburg,
c     wasniewski and zlatev, harwell report css 135, 1983).
c
c ma30a/ad performs the lu decomposition of the diagonal blocks of the
c     permutation paq of a sparse matrix a, where input permutations
c     p1 and q1 are used to define the diagonal blocks.  there may be
c     non-zeros in the off-diagonal blocks but they are unaffected by
c     ma30a/ad. p and p1 differ only within blocks as do q and q1. the
c     permutations p1 and q1 may be found by calling mc23a/ad or the
c     matrix may be treated as a single block by using p1=q1=i. the
c     matrix non-zeros should be held compactly by rows, although it
c     should be noted that the user can supply the matrix by columns
c     to get the lu decomposition of a transpose.
c
c the parameters are...
c this description should also be consulted for further information on
c     most of the parameters of ma30b/bd and ma30c/cd.
c
c n  is an integer variable which must be set by the user to the order
c     of the matrix.  it is not altered by ma30a/ad.
c icn is an integer array of length licn. positions idisp(2) to
c     licn must be set by the user to contain the column indices of
c     the non-zeros in the diagonal blocks of p1*a*q1. those belonging
c     to a single row must be contiguous but the ordering of column
c     indices with each row is unimportant. the non-zeros of row i
c     precede those of row i+1,i=1,...,n-1 and no wasted space is
c     allowed between the rows.  on output the column indices of the
c     lu decomposition of paq are held in positions idisp(1) to
c     idisp(2), the rows are in pivotal order, and the column indices
c     of the l part of each row are in pivotal order and precede those
c     of u. again there is no wasted space either within a row or
c     between the rows. icn(1) to icn(idisp(1)-1), are neither
c     required nor altered. if mc23a/ad been called, these will hold
c     information about the off-diagonal blocks.
c a is a real/double precision array of length licn whose entries
c     idisp(2) to licn must be set by the user to the  values of the
c     non-zero entries of the matrix in the order indicated by  icn.
c     on output a will hold the lu factors of the matrix where again
c     the position in the matrix is determined by the corresponding
c     values in icn. a(1) to a(idisp(1)-1) are neither required nor
c     altered.
c licn  is an integer variable which must be set by the user to the
c     length of arrays icn and a. it must be big enough for a and icn
c     to hold all the non-zeros of l and u and leave some "elbow
c     room".  it is possible to calculate a minimum value for licn by
c     a preliminary run of ma30a/ad. the adequacy of the elbow room
c     can be judged by the size of the common block variable icncp. it
c     is not altered by ma30a/ad.
c lenr  is an integer array of length n.  on input, lenr(i) should
c     equal the number of non-zeros in row i, i=1,...,n of the
c     diagonal blocks of p1*a*q1. on output, lenr(i) will equal the
c     total number of non-zeros in row i of l and row i of u.
c lenrl  is an integer array of length n. on output from ma30a/ad,
c     lenrl(i) will hold the number of non-zeros in row i of l.
c idisp  is an integer array of length 2. the user should set idisp(1)
c     to be the first available position in a/icn for the lu
c     decomposition while idisp(2) is set to the position in a/icn of
c     the first non-zero in the diagonal blocks of p1*a*q1. on output,
c     idisp(1) will be unaltered while idisp(2) will be set to the
c     position in a/icn of the last non-zero of the lu decomposition.
c ip  is an integer array of length n which holds a permutation of
c     the integers 1 to n.  on input to ma30a/ad, the absolute value of
c     ip(i) must be set to the row of a which is row i of p1*a*q1. a
c     negative value for ip(i) indicates that row i is at the end of a
c     diagonal block.  on output from ma30a/ad, ip(i) indicates the row
c     of a which is the i th row in paq. ip(i) will still be negative
c     for the last row of each block (except the last).
c iq is an integer array of length n which again holds a
c     permutation of the integers 1 to n.  on input to ma30a/ad, iq(j)
c     must be set to the column of a which is column j of p1*a*q1. on
c     output from ma30a/ad, the absolute value of iq(j) indicates the
c     column of a which is the j th in paq.  for rows, i say, in which
c     structural or numerical singularity is detected iq(i) is
c     negated.
c irn  is an integer array of length lirn used as workspace by
c     ma30a/ad.
c lirn  is an integer variable. it should be greater than the
c     largest number of non-zeros in a diagonal block of p1*a*q1 but
c     need not be as large as licn. it is the length of array irn and
c     should be large enough to hold the active part of any block,
c     plus some "elbow room", the  a posteriori  adequacy of which can
c     be estimated by examining the size of common block variable
c     irncp.
c lenc,ifirst,lastr,nextr,lastc,nextc are all integer arrays of
c     length n which are used as workspace by ma30a/ad.  if nsrch is
c     set to a value less than or equal to n, then arrays lastc and
c     nextc are not referenced by ma30a/ad and so can be dummied in
c     the call to ma30a/ad.
c iptr,ipc are integer arrays of length n which are used as workspace
c     by ma30a/ad.
c u  is a real/double precision variable which should be set by the
c     user to a value between 0. and 1.0. if less than zero it is
c     reset to zero and if its value is 1.0 or greater it is reset to
c     0.9999 (0.999999999 in d version).  it determines the balance
c     between pivoting for sparsity and for stability, values near
c     zero emphasizing sparsity and values near one emphasizing
c     stability. we recommend u=0.1 as a posible first trial value.
c     the stability can be judged by a later call to mc24a/ad or by
c     setting lbig to .true.
c iflag  is an integer variable. it will have a non-negative value if
c     ma30a/ad is successful. negative values indicate error
c     conditions while positive values indicate that the matrix has
c     been successfully decomposed but is singular. for each non-zero
c     value, an appropriate message is output on unit lp.  possible
c     non-zero values for iflag are ...
c
c -1  the matrix is structually singular with rank given by irank in
c     common block ma30f/fd.
c +1  if, however, the user wants the lu decomposition of a
c     structurally singular matrix and sets the common block variable
c     abort1 to .false., then, in the event of singularity and a
c     successful decomposition, iflag is returned with the value +1
c     and no message is output.
c -2  the matrix is numerically singular (it may also be structually
c     singular) with estimated rank given by irank in common block
c     ma30f/fd.
c +2  the  user can choose to continue the decomposition even when a
c     zero pivot is encountered by setting common block variable
c     abort2 to .false.  if a singularity is encountered, iflag will
c     then return with a value of +2, and no message is output if the
c     decomposition has been completed successfully.
c -3  lirn has not been large enough to continue with the
c     decomposition.  if the stage was zero then common block variable
c     minirn gives the length sufficient to start the decomposition on
c     this block.  for a successful decomposition on this block the
c     user should make lirn slightly (say about n/2) greater than this
c     value.
c -4  licn not large enough to continue with the decomposition.
c -5  the decomposition has been completed but some of the lu factors
c     have been discarded to create enough room in a/icn to continue
c     the decomposition. the variable minicn in common block ma30f/fd
c     then gives the size that licn should be to enable the
c     factorization to be successful.  if the user sets common block
c     variable abort3 to .true., then the subroutine will exit
c     immediately instead of destroying any factors and continuing.
c -6  both licn and lirn are too small. termination has been caused by
c     lack of space in irn (see error iflag= -3), but already some of
c     the lu factors in a/icn have been lost (see error iflag= -5).
c     minicn gives the minimum amount of space required in a/icn for
c     decomposition up to this point.
c
      DOUBLE PRECISION A(LICN), U, AU, UMAX, AMAX, ZERO, PIVRAT, PIVR,
     * TOL, BIG, ANEW, AANEW, SCALE
      INTEGER IPTR(NN), PIVOT, PIVEND, DISPC, OLDPIV, OLDEND, PIVROW,
     * ROWI, IPC(NN), IDISP(2), COLUPD
      INTEGER ICN(LICN), LENR(NN), LENRL(NN), IP(NN), IQ(NN),
     * LENC(NN), IRN(LIRN), IFIRST(NN), LASTR(NN), NEXTR(NN),
     * LASTC(NN), NEXTC(NN)
      LOGICAL ABORT1, ABORT2, ABORT3, LBIG
c for comments of common block variables see block data subprogram.
      COMMON /MA30ED/ LP, ABORT1, ABORT2, ABORT3
      COMMON /MA30FD/ IRNCP, ICNCP, IRANK, MINIRN, MINICN
      COMMON /MA30ID/ TOL, BIG, NDROP, NSRCH, LBIG
      COMMON /LPIVOT/ LPIV(10),LNPIV(10),MAPIV,MANPIV,IAVPIV,
     *                IANPIV,KOUNTL
c
      DATA UMAX/.999999999D0/
      DATA ZERO /0.0D0/
      MSRCH = NSRCH
      NDROP = 0
      DO 1272 KK=1,10
        LNPIV(KK)=0
        LPIV(KK)=0
 1272 CONTINUE
      MAPIV = 0
      MANPIV = 0
      IAVPIV = 0
      IANPIV = 0
      KOUNTL = 0
      MINIRN = 0
      MINICN = IDISP(1) - 1
      MOREI = 0
      IRANK = NN
      IRNCP = 0
      ICNCP = 0
      IFLAG = 0
c reset u if necessary.
      U = DMIN1(U,UMAX)
c ibeg is the position of the next pivot row after elimination step
c     using it.
      U = DMAX1(U,ZERO)
      IBEG = IDISP(1)
c iactiv is the position of the first entry in the active part of a/icn.
      IACTIV = IDISP(2)
c nzrow is current number of non-zeros in active and unprocessed part
c     of row file icn.
      NZROW = LICN - IACTIV + 1
      MINICN = NZROW + MINICN
c
c count the number of diagonal blocks and set up pointers to the
c     beginnings of the rows.
c num is the number of diagonal blocks.
      NUM = 1
      IPTR(1) = IACTIV
      IF (NN.EQ.1) GO TO 20
      NNM1 = NN - 1
      DO 10 I=1,NNM1
        IF (IP(I).LT.0) NUM = NUM + 1
        IPTR(I+1) = IPTR(I) + LENR(I)
   10 CONTINUE
c ilast is the last row in the previous block.
   20 ILAST = 0
c
c ***********************************************
c ****    lu decomposition of block nblock   ****
c ***********************************************
c
c each pass through this loop performs lu decomposition on one
c     of the diagonal blocks.
      DO 1000 NBLOCK=1,NUM
        ISTART = ILAST + 1
        DO 30 IROWS=ISTART,NN
          IF (IP(IROWS).LT.0) GO TO 40
   30   CONTINUE
        IROWS = NN
   40   ILAST = IROWS
c n is the number of rows in the current block.
c istart is the index of the first row in the current block.
c ilast is the index of the last row in the current block.
c iactiv is the position of the first entry in the block.
c itop is the position of the last entry in the block.
        N = ILAST - ISTART + 1
        IF (N.NE.1) GO TO 90
c
c code for dealing with 1x1 block.
        LENRL(ILAST) = 0
        ISING = ISTART
        IF (LENR(ILAST).NE.0) GO TO 50
c block is structurally singular.
        IRANK = IRANK - 1
        ISING = -ISING
        IF (IFLAG.NE.2 .AND. IFLAG.NE.-5) IFLAG = 1
        IF (.NOT.ABORT1) GO TO 80
        IDISP(2) = IACTIV
        IFLAG = -1
        IF (LP.NE.0) WRITE (LP,99999)
c     return
        GO TO 1120
   50   SCALE = DABS(A(IACTIV))
        IF (SCALE.EQ.ZERO) GO TO 60
        IF (LBIG) BIG = DMAX1(BIG,SCALE)
        GO TO 70
   60   ISING = -ISING
        IRANK = IRANK - 1
        IPTR(ILAST) = 0
        IF (IFLAG.NE.-5) IFLAG = 2
        IF (.NOT.ABORT2) GO TO 70
        IDISP(2) = IACTIV
        IFLAG = -2
        IF (LP.NE.0) WRITE (LP,99998)
        GO TO 1120
   70   A(IBEG) = A(IACTIV)
        ICN(IBEG) = ICN(IACTIV)
        IACTIV = IACTIV + 1
        IPTR(ISTART) = 0
        IBEG = IBEG + 1
        NZROW = NZROW - 1
   80   LASTR(ISTART) = ISTART
        IPC(ISTART) = -ISING
        GO TO 1000
c
c non-trivial block.
   90   ITOP = LICN
        IF (ILAST.NE.NN) ITOP = IPTR(ILAST+1) - 1
c
c set up column oriented storage.
        DO 100 I=ISTART,ILAST
          LENRL(I) = 0
          LENC(I) = 0
  100   CONTINUE
        IF (ITOP-IACTIV.LT.LIRN) GO TO 110
        MINIRN = ITOP - IACTIV + 1
        PIVOT = ISTART - 1
        GO TO 1100
c
c calculate column counts.
  110   DO 120 II=IACTIV,ITOP
          I = ICN(II)
          LENC(I) = LENC(I) + 1
  120   CONTINUE
c set up column pointers so that ipc(j) points to position after end
c     of column j in column file.
        IPC(ILAST) = LIRN + 1
        J1 = ISTART + 1
        DO 130 JJ=J1,ILAST
          J = ILAST - JJ + J1 - 1
          IPC(J) = IPC(J+1) - LENC(J+1)
  130   CONTINUE
        DO 150 INDROW=ISTART,ILAST
          J1 = IPTR(INDROW)
          J2 = J1 + LENR(INDROW) - 1
          IF (J1.GT.J2) GO TO 150
          DO 140 JJ=J1,J2
            J = ICN(JJ)
            IPOS = IPC(J) - 1
            IRN(IPOS) = INDROW
            IPC(J) = IPOS
  140     CONTINUE
  150   CONTINUE
c dispc is the lowest indexed active location in the column file.
        DISPC = IPC(ISTART)
        NZCOL = LIRN - DISPC + 1
        MINIRN = MAX0(NZCOL,MINIRN)
        NZMIN = 1
c
c initialize array ifirst.  ifirst(i) = +/- k indicates that row/col
c     k has i non-zeros.  if ifirst(i) = 0, there is no row or column
c     with i non zeros.
        DO 160 I=1,N
          IFIRST(I) = 0
  160   CONTINUE
c
c compute ordering of row and column counts.
c first run through columns (from column n to column 1).
        DO 180 JJ=ISTART,ILAST
          J = ILAST - JJ + ISTART
          NZ = LENC(J)
          IF (NZ.NE.0) GO TO 170
          IPC(J) = 0
          GO TO 180
  170     IF (NSRCH.LE.NN) GO TO 180
          ISW = IFIRST(NZ)
          IFIRST(NZ) = -J
          LASTC(J) = 0
          NEXTC(J) = -ISW
          ISW1 = IABS(ISW)
          IF (ISW.NE.0) LASTC(ISW1) = J
  180   CONTINUE
c now run through rows (again from n to 1).
        DO 210 II=ISTART,ILAST
          I = ILAST - II + ISTART
          NZ = LENR(I)
          IF (NZ.NE.0) GO TO 190
          IPTR(I) = 0
          LASTR(I) = 0
          GO TO 210
  190     ISW = IFIRST(NZ)
          IFIRST(NZ) = I
          IF (ISW.GT.0) GO TO 200
          NEXTR(I) = 0
          LASTR(I) = ISW
          GO TO 210
  200     NEXTR(I) = ISW
          LASTR(I) = LASTR(ISW)
          LASTR(ISW) = I
  210   CONTINUE
c
c
c **********************************************
c ****    start of main elimination loop    ****
c **********************************************
        DO 980 PIVOT=ISTART,ILAST
c
c first find the pivot using markowitz criterion with stability
c     control.
c jcost is the markowitz cost of the best pivot so far,.. this
c     pivot is in row ipiv and column jpiv.
          NZ2 = NZMIN
          JCOST = N*N
c
c examine rows/columns in order of ascending count.
          DO 340 L=1,2
            PIVRAT = ZERO
            ISRCH = 1
            LL = L
c a pass with l equal to 2 is only performed in the case of singularity.
            DO 330 NZ=NZ2,N
              IF (JCOST.LE.(NZ-1)**2) GO TO 420
              IJFIR = IFIRST(NZ)
              IF (IJFIR) 230, 220, 240
  220         IF (LL.EQ.1) NZMIN = NZ + 1
              GO TO 330
  230         LL = 2
              IJFIR = -IJFIR
              GO TO 290
  240         LL = 2
c scan rows with nz non-zeros.
              DO 270 IDUMMY=1,N
                IF (JCOST.LE.(NZ-1)**2) GO TO 420
                IF (ISRCH.GT.MSRCH) GO TO 420
                IF (IJFIR.EQ.0) GO TO 280
c row ijfir is now examined.
                I = IJFIR
                IJFIR = NEXTR(I)
c first calculate multiplier threshold level.
                AMAX = ZERO
                J1 = IPTR(I) + LENRL(I)
                J2 = IPTR(I) + LENR(I) - 1
                DO 250 JJ=J1,J2
                  AMAX = DMAX1(AMAX,DABS(A(JJ)))
  250           CONTINUE
                AU = AMAX*U
                ISRCH = ISRCH + 1
c scan row for possible pivots
                DO 260 JJ=J1,J2
                  IF (DABS(A(JJ)).LE.AU .AND. L.EQ.1) GO TO 260
                  J = ICN(JJ)
                  KCOST = (NZ-1)*(LENC(J)-1)
                  IF (KCOST.GT.JCOST) GO TO 260
                  PIVR = ZERO
                  IF (AMAX.NE.ZERO) PIVR = DABS(A(JJ))/AMAX
                  IF (KCOST.EQ.JCOST .AND. (PIVR.LE.PIVRAT .OR.
     *             NSRCH.GT.NN+1)) GO TO 260
c best pivot so far is found.
                  JCOST = KCOST
                  IJPOS = JJ
                  IPIV = I
                  JPIV = J
                  IF (MSRCH.GT.NN+1 .AND. JCOST.LE.(NZ-1)**2) GO TO 420
                  PIVRAT = PIVR
  260           CONTINUE
  270         CONTINUE
c
c columns with nz non-zeros now examined.
  280         IJFIR = IFIRST(NZ)
              IJFIR = -LASTR(IJFIR)
  290         IF (JCOST.LE.NZ*(NZ-1)) GO TO 420
              IF (MSRCH.LE.NN) GO TO 330
              DO 320 IDUMMY=1,N
                IF (IJFIR.EQ.0) GO TO 330
                J = IJFIR
                IJFIR = NEXTC(IJFIR)
                I1 = IPC(J)
                I2 = I1 + NZ - 1
c scan column j.
                DO 310 II=I1,I2
                  I = IRN(II)
                  KCOST = (NZ-1)*(LENR(I)-LENRL(I)-1)
                  IF (KCOST.GE.JCOST) GO TO 310
c pivot has best markowitz count so far ... now check its
c     suitability on numeric grounds by examining the other non-zeros
c     in its row.
                  J1 = IPTR(I) + LENRL(I)
                  J2 = IPTR(I) + LENR(I) - 1
c we need a stability check on singleton columns because of possible
c     problems with underdetermined systems.
                  AMAX = ZERO
                  DO 300 JJ=J1,J2
                    AMAX = DMAX1(AMAX,DABS(A(JJ)))
                    IF (ICN(JJ).EQ.J) JPOS = JJ
  300             CONTINUE
                  IF (DABS(A(JPOS)).LE.AMAX*U .AND. L.EQ.1) GO TO 310
                  JCOST = KCOST
                  IPIV = I
                  JPIV = J
                  IJPOS = JPOS
                  IF (AMAX.NE.ZERO) PIVRAT = DABS(A(JPOS))/AMAX
                  IF (JCOST.LE.NZ*(NZ-1)) GO TO 420
  310           CONTINUE
c
  320         CONTINUE
c
  330       CONTINUE
c in the event of singularity, we must make sure all rows and columns
c are tested.
            MSRCH = N
c
c matrix is numerically or structurally singular  ... which it is will
c     be diagnosed later.
            IRANK = IRANK - 1
  340     CONTINUE
c assign rest of rows and columns to ordering array.
c matrix is structurally singular.
          IF (IFLAG.NE.2 .AND. IFLAG.NE.-5) IFLAG = 1
          IRANK = IRANK - ILAST + PIVOT + 1
          IF (.NOT.ABORT1) GO TO 350
          IDISP(2) = IACTIV
          IFLAG = -1
          IF (LP.NE.0) WRITE (LP,99999)
          GO TO 1120
  350     K = PIVOT - 1
          DO 390 I=ISTART,ILAST
            IF (LASTR(I).NE.0) GO TO 390
            K = K + 1
            LASTR(I) = K
            IF (LENRL(I).EQ.0) GO TO 380
            MINICN = MAX0(MINICN,NZROW+IBEG-1+MOREI+LENRL(I))
            IF (IACTIV-IBEG.GE.LENRL(I)) GO TO 360
            CALL MA30DD(A, ICN, IPTR(ISTART), N, IACTIV, ITOP, .TRUE.)
c check now to see if ma30d/dd has created enough available space.
            IF (IACTIV-IBEG.GE.LENRL(I)) GO TO 360
c create more space by destroying previously created lu factors.
            MOREI = MOREI + IBEG - IDISP(1)
            IBEG = IDISP(1)
            IF (LP.NE.0) WRITE (LP,99997)
            IFLAG = -5
            IF (ABORT3) GO TO 1090
  360       J1 = IPTR(I)
            J2 = J1 + LENRL(I) - 1
            IPTR(I) = 0
            DO 370 JJ=J1,J2
              A(IBEG) = A(JJ)
              ICN(IBEG) = ICN(JJ)
              ICN(JJ) = 0
              IBEG = IBEG + 1
  370       CONTINUE
            NZROW = NZROW - LENRL(I)
  380       IF (K.EQ.ILAST) GO TO 400
  390     CONTINUE
  400     K = PIVOT - 1
          DO 410 I=ISTART,ILAST
            IF (IPC(I).NE.0) GO TO 410
            K = K + 1
            IPC(I) = K
            IF (K.EQ.ILAST) GO TO 990
  410     CONTINUE
c
c the pivot has now been found in position (ipiv,jpiv) in location
c     ijpos in row file.
c update column and row ordering arrays to correspond with removal
c     of the active part of the matrix.
  420     ISING = PIVOT
          IF (A(IJPOS).NE.ZERO) GO TO 430
c numerical singularity is recorded here.
          ISING = -ISING
          IF (IFLAG.NE.-5) IFLAG = 2
          IF (.NOT.ABORT2) GO TO 430
          IDISP(2) = IACTIV
          IFLAG = -2
          IF (LP.NE.0) WRITE (LP,99998)
          GO TO 1120
  430     OLDPIV = IPTR(IPIV) + LENRL(IPIV)
          OLDEND = IPTR(IPIV) + LENR(IPIV) - 1
c changes to column ordering.
          IF (NSRCH.LE.NN) GO TO 460
          COLUPD = NN + 1
            LENPP = OLDEND-OLDPIV+1
            IF (LENPP.LT.4) LPIV(1) = LPIV(1) + 1
            IF (LENPP.GE.4 .AND. LENPP.LE.6) LPIV(2) = LPIV(2) + 1
            IF (LENPP.GE.7 .AND. LENPP.LE.10) LPIV(3) = LPIV(3) + 1
            IF (LENPP.GE.11 .AND. LENPP.LE.15) LPIV(4) = LPIV(4) + 1
            IF (LENPP.GE.16 .AND. LENPP.LE.20) LPIV(5) = LPIV(5) + 1
            IF (LENPP.GE.21 .AND. LENPP.LE.30) LPIV(6) = LPIV(6) + 1
            IF (LENPP.GE.31 .AND. LENPP.LE.50) LPIV(7) = LPIV(7) + 1
            IF (LENPP.GE.51 .AND. LENPP.LE.70) LPIV(8) = LPIV(8) + 1
            IF (LENPP.GE.71 .AND. LENPP.LE.100) LPIV(9) = LPIV(9) + 1
            IF (LENPP.GE.101) LPIV(10) = LPIV(10) + 1
            MAPIV = MAX0(MAPIV,LENPP)
            IAVPIV = IAVPIV + LENPP
          DO 450 JJ=OLDPIV,OLDEND
            J = ICN(JJ)
            LC = LASTC(J)
            NC = NEXTC(J)
            NEXTC(J) = -COLUPD
            IF (JJ.NE.IJPOS) COLUPD = J
            IF (NC.NE.0) LASTC(NC) = LC
            IF (LC.EQ.0) GO TO 440
            NEXTC(LC) = NC
            GO TO 450
  440       NZ = LENC(J)
            ISW = IFIRST(NZ)
            IF (ISW.GT.0) LASTR(ISW) = -NC
            IF (ISW.LT.0) IFIRST(NZ) = -NC
  450     CONTINUE
c changes to row ordering.
  460     I1 = IPC(JPIV)
          I2 = I1 + LENC(JPIV) - 1
          DO 480 II=I1,I2
            I = IRN(II)
            LR = LASTR(I)
            NR = NEXTR(I)
            IF (NR.NE.0) LASTR(NR) = LR
            IF (LR.LE.0) GO TO 470
            NEXTR(LR) = NR
            GO TO 480
  470       NZ = LENR(I) - LENRL(I)
            IF (NR.NE.0) IFIRST(NZ) = NR
            IF (NR.EQ.0) IFIRST(NZ) = LR
  480     CONTINUE
c
c move pivot to position lenrl+1 in pivot row and move pivot row
c     to the beginning of the available storage.
c the l part and the pivot in the old copy of the pivot row is
c     nullified while, in the strictly upper triangular part, the
c     column indices, j say, are overwritten by the corresponding
c     entry of iq (iq(j)) and iq(j) is set to the negative of the
c     displacement of the column index from the pivot entry.
          IF (OLDPIV.EQ.IJPOS) GO TO 490
          AU = A(OLDPIV)
          A(OLDPIV) = A(IJPOS)
          A(IJPOS) = AU
          ICN(IJPOS) = ICN(OLDPIV)
          ICN(OLDPIV) = JPIV
c check to see if there is space immediately available in a/icn to
c     hold new copy of pivot row.
  490     MINICN = MAX0(MINICN,NZROW+IBEG-1+MOREI+LENR(IPIV))
          IF (IACTIV-IBEG.GE.LENR(IPIV)) GO TO 500
          CALL MA30DD(A, ICN, IPTR(ISTART), N, IACTIV, ITOP, .TRUE.)
          OLDPIV = IPTR(IPIV) + LENRL(IPIV)
          OLDEND = IPTR(IPIV) + LENR(IPIV) - 1
c check now to see if ma30d/dd has created enough available space.
          IF (IACTIV-IBEG.GE.LENR(IPIV)) GO TO 500
c create more space by destroying previously created lu factors.
          MOREI = MOREI + IBEG - IDISP(1)
          IBEG = IDISP(1)
          IF (LP.NE.0) WRITE (LP,99997)
          IFLAG = -5
          IF (ABORT3) GO TO 1090
          IF (IACTIV-IBEG.GE.LENR(IPIV)) GO TO 500
c there is still not enough room in a/icn.
          IFLAG = -4
          GO TO 1090
c copy pivot row and set up iq array.
  500     IJPOS = 0
          J1 = IPTR(IPIV)
c
          DO 530 JJ=J1,OLDEND
            A(IBEG) = A(JJ)
            ICN(IBEG) = ICN(JJ)
            IF (IJPOS.NE.0) GO TO 510
            IF (ICN(JJ).EQ.JPIV) IJPOS = IBEG
            ICN(JJ) = 0
            GO TO 520
  510       K = IBEG - IJPOS
            J = ICN(JJ)
            ICN(JJ) = IQ(J)
            IQ(J) = -K
  520       IBEG = IBEG + 1
  530     CONTINUE
c
          IJP1 = IJPOS + 1
          PIVEND = IBEG - 1
          LENPIV = PIVEND - IJPOS
          NZROW = NZROW - LENRL(IPIV) - 1
          IPTR(IPIV) = OLDPIV + 1
          IF (LENPIV.EQ.0) IPTR(IPIV) = 0
c
c remove pivot row (including pivot) from column oriented file.
          DO 560 JJ=IJPOS,PIVEND
            J = ICN(JJ)
            I1 = IPC(J)
            LENC(J) = LENC(J) - 1
c i2 is last position in new column.
            I2 = IPC(J) + LENC(J) - 1
            IF (I2.LT.I1) GO TO 550
            DO 540 II=I1,I2
              IF (IRN(II).NE.IPIV) GO TO 540
              IRN(II) = IRN(I2+1)
              GO TO 550
  540       CONTINUE
  550       IRN(I2+1) = 0
  560     CONTINUE
          NZCOL = NZCOL - LENPIV - 1
c
c go down the pivot column and for each row with a non-zero add
c     the appropriate multiple of the pivot row to it.
c we loop on the number of non-zeros in the pivot column since
c     ma30d/dd may change its actual position.
c
          NZPC = LENC(JPIV)
          IF (NZPC.EQ.0) GO TO 900
          DO 840 III=1,NZPC
            II = IPC(JPIV) + III - 1
            I = IRN(II)
c search row i for non-zero to be eliminated, calculate multiplier,
c     and place it in position lenrl+1 in its row.
c  idrop is the number of non-zero entries dropped from row    i
c        because these fall beneath tolerance level.
c
            IDROP = 0
            J1 = IPTR(I) + LENRL(I)
            IEND = IPTR(I) + LENR(I) - 1
            DO 570 JJ=J1,IEND
              IF (ICN(JJ).NE.JPIV) GO TO 570
c if pivot is zero, rest of column is and so multiplier is zero.
              AU = ZERO
              IF (A(IJPOS).NE.ZERO) AU = -A(JJ)/A(IJPOS)
              IF (LBIG) BIG = DMAX1(BIG,DABS(AU))
              A(JJ) = A(J1)
              A(J1) = AU
              ICN(JJ) = ICN(J1)
              ICN(J1) = JPIV
              LENRL(I) = LENRL(I) + 1
              GO TO 580
  570       CONTINUE
c jump if pivot row is a singleton.
  580       IF (LENPIV.EQ.0) GO TO 840
c now perform necessary operations on rest of non-pivot row i.
            ROWI = J1 + 1
            IOP = 0
c jump if all the pivot row causes fill-in.
            IF (ROWI.GT.IEND) GO TO 650
c perform operations on current non-zeros in row i.
c innermost loop.
            LENPP = IEND-ROWI+1
            IF (LENPP.LT.4) LNPIV(1) = LNPIV(1) + 1
            IF (LENPP.GE.4 .AND. LENPP.LE.6) LNPIV(2) = LNPIV(2) + 1
            IF (LENPP.GE.7 .AND. LENPP.LE.10) LNPIV(3) = LNPIV(3) + 1
            IF (LENPP.GE.11 .AND. LENPP.LE.15) LNPIV(4) = LNPIV(4) + 1
            IF (LENPP.GE.16 .AND. LENPP.LE.20) LNPIV(5) = LNPIV(5) + 1
            IF (LENPP.GE.21 .AND. LENPP.LE.30) LNPIV(6) = LNPIV(6) + 1
            IF (LENPP.GE.31 .AND. LENPP.LE.50) LNPIV(7) = LNPIV(7) + 1
            IF (LENPP.GE.51 .AND. LENPP.LE.70) LNPIV(8) = LNPIV(8) + 1
            IF (LENPP.GE.71 .AND. LENPP.LE.100) LNPIV(9) = LNPIV(9) + 1
            IF (LENPP.GE.101) LNPIV(10) = LNPIV(10) + 1
            MANPIV = MAX0(MANPIV,LENPP)
            IANPIV = IANPIV + LENPP
            KOUNTL = KOUNTL + 1
            DO 590 JJ=ROWI,IEND
              J = ICN(JJ)
              IF (IQ(J).GT.0) GO TO 590
              IOP = IOP + 1
              PIVROW = IJPOS - IQ(J)
              A(JJ) = A(JJ) + AU*A(PIVROW)
              IF (LBIG) BIG = DMAX1(DABS(A(JJ)),BIG)
              ICN(PIVROW) = -ICN(PIVROW)
              IF (DABS(A(JJ)).LT.TOL) IDROP = IDROP + 1
  590       CONTINUE
c
c  jump if no non-zeros in non-pivot row have been removed
c       because these are beneath the drop-tolerance  tol.
c
            IF (IDROP.EQ.0) GO TO 650
c
c  run through non-pivot row compressing row so that only
c      non-zeros greater than   tol   are stored.  all non-zeros
c      less than   tol   are also removed from the column structure.
c
            JNEW = ROWI
            DO 630 JJ=ROWI,IEND
              IF (DABS(A(JJ)).LT.TOL) GO TO 600
              A(JNEW) = A(JJ)
              ICN(JNEW) = ICN(JJ)
              JNEW = JNEW + 1
              GO TO 630
c
c  remove non-zero entry from column structure.
c
  600         J = ICN(JJ)
              I1 = IPC(J)
              I2 = I1 + LENC(J) - 1
              DO 610 II=I1,I2
                IF (IRN(II).EQ.I) GO TO 620
  610         CONTINUE
  620         IRN(II) = IRN(I2)
              IRN(I2) = 0
              LENC(J) = LENC(J) - 1
              IF (NSRCH.LE.NN) GO TO 630
c remove column from column chain and place in update chain.
              IF (NEXTC(J).LT.0) GO TO 630
c jump if column already in update chain.
              LC = LASTC(J)
              NC = NEXTC(J)
              NEXTC(J) = -COLUPD
              COLUPD = J
              IF (NC.NE.0) LASTC(NC) = LC
              IF (LC.EQ.0) GO TO 622
              NEXTC(LC) = NC
              GO TO 630
  622         NZ = LENC(J) + 1
              ISW = IFIRST(NZ)
              IF (ISW.GT.0) LASTR(ISW) = -NC
              IF (ISW.LT.0) IFIRST(NZ) = -NC
  630       CONTINUE
            DO 640 JJ=JNEW,IEND
              ICN(JJ) = 0
  640       CONTINUE
c the value of idrop might be different from that calculated earlier
c     because, we may now have dropped some non-zeros which were not
c     modified by the pivot row.
            IDROP = IEND + 1 - JNEW
            IEND = JNEW - 1
            LENR(I) = LENR(I) - IDROP
            NZROW = NZROW - IDROP
            NZCOL = NZCOL - IDROP
            NDROP = NDROP + IDROP
  650       IFILL = LENPIV - IOP
c jump is if there is no fill-in.
            IF (IFILL.EQ.0) GO TO 750
c now for the fill-in.
            MINICN = MAX0(MINICN,MOREI+IBEG-1+NZROW+IFILL+LENR(I))
c see if there is room for fill-in.
c get maximum space for row i in situ.
            DO 660 JDIFF=1,IFILL
              JNPOS = IEND + JDIFF
              IF (JNPOS.GT.LICN) GO TO 670
              IF (ICN(JNPOS).NE.0) GO TO 670
  660       CONTINUE
c there is room for all the fill-in after the end of the row so it
c     can be left in situ.
c next available space for fill-in.
            IEND = IEND + 1
            GO TO 750
c jmore spaces for fill-in are required in front of row.
  670       JMORE = IFILL - JDIFF + 1
            I1 = IPTR(I)
c we now look in front of the row to see if there is space for
c     the rest of the fill-in.
            DO 680 JDIFF=1,JMORE
              JNPOS = I1 - JDIFF
              IF (JNPOS.LT.IACTIV) GO TO 690
              IF (ICN(JNPOS).NE.0) GO TO 700
  680       CONTINUE
  690       JNPOS = I1 - JMORE
            GO TO 710
c whole row must be moved to the beginning of available storage.
  700       JNPOS = IACTIV - LENR(I) - IFILL
c jump if there is space immediately available for the shifted row.
  710       IF (JNPOS.GE.IBEG) GO TO 730
            CALL MA30DD(A, ICN, IPTR(ISTART), N, IACTIV, ITOP, .TRUE.)
            I1 = IPTR(I)
            IEND = I1 + LENR(I) - 1
            JNPOS = IACTIV - LENR(I) - IFILL
            IF (JNPOS.GE.IBEG) GO TO 730
c no space available so try to create some by throwing away previous
c     lu decomposition.
            MOREI = MOREI + IBEG - IDISP(1) - LENPIV - 1
            IF (LP.NE.0) WRITE (LP,99997)
            IFLAG = -5
            IF (ABORT3) GO TO 1090
c keep record of current pivot row.
            IBEG = IDISP(1)
            ICN(IBEG) = JPIV
            A(IBEG) = A(IJPOS)
            IJPOS = IBEG
            DO 720 JJ=IJP1,PIVEND
              IBEG = IBEG + 1
              A(IBEG) = A(JJ)
              ICN(IBEG) = ICN(JJ)
  720       CONTINUE
            IJP1 = IJPOS + 1
            PIVEND = IBEG
            IBEG = IBEG + 1
            IF (JNPOS.GE.IBEG) GO TO 730
c this still does not give enough room.
            IFLAG = -4
            GO TO 1090
  730       IACTIV = MIN0(IACTIV,JNPOS)
c move non-pivot row i.
            IPTR(I) = JNPOS
            DO 740 JJ=I1,IEND
              A(JNPOS) = A(JJ)
              ICN(JNPOS) = ICN(JJ)
              JNPOS = JNPOS + 1
              ICN(JJ) = 0
  740       CONTINUE
c first new available space.
            IEND = JNPOS
  750       NZROW = NZROW + IFILL
c innermost fill-in loop which also resets icn.
            IDROP = 0
            DO 830 JJ=IJP1,PIVEND
              J = ICN(JJ)
              IF (J.LT.0) GO TO 820
              ANEW = AU*A(JJ)
              AANEW = DABS(ANEW)
              IF (AANEW.GE.TOL) GO TO 760
              IDROP = IDROP + 1
              NDROP = NDROP + 1
              NZROW = NZROW - 1
              MINICN = MINICN - 1
              IFILL = IFILL - 1
              GO TO 830
  760         IF (LBIG) BIG = DMAX1(AANEW,BIG)
              A(IEND) = ANEW
              ICN(IEND) = J
              IEND = IEND + 1
c
c put new entry in column file.
              MINIRN = MAX0(MINIRN,NZCOL+LENC(J)+1)
              JEND = IPC(J) + LENC(J)
              JROOM = NZPC - III + 1 + LENC(J)
              IF (JEND.GT.LIRN) GO TO 770
              IF (IRN(JEND).EQ.0) GO TO 810
  770         IF (JROOM.LT.DISPC) GO TO 780
c compress column file to obtain space for new copy of column.
              CALL MA30DD(A, IRN, IPC(ISTART), N, DISPC, LIRN, .FALSE.)
              IF (JROOM.LT.DISPC) GO TO 780
              JROOM = DISPC - 1
              IF (JROOM.GE.LENC(J)+1) GO TO 780
c column file is not large enough.
              GO TO 1100
c copy column to beginning of file.
  780         JBEG = IPC(J)
              JEND = IPC(J) + LENC(J) - 1
              JZERO = DISPC - 1
              DISPC = DISPC - JROOM
              IDISPC = DISPC
              DO 790 II=JBEG,JEND
                IRN(IDISPC) = IRN(II)
                IRN(II) = 0
                IDISPC = IDISPC + 1
  790         CONTINUE
              IPC(J) = DISPC
              JEND = IDISPC
              DO 800 II=JEND,JZERO
                IRN(II) = 0
  800         CONTINUE
  810         IRN(JEND) = I
              NZCOL = NZCOL + 1
              LENC(J) = LENC(J) + 1
c end of adjustment to column file.
              GO TO 830
c
  820         ICN(JJ) = -J
  830       CONTINUE
            IF (IDROP.EQ.0) GO TO 834
            DO 832 KDROP=1,IDROP
            ICN(IEND) = 0
            IEND = IEND + 1
  832       CONTINUE
  834       LENR(I) = LENR(I) + IFILL
c end of scan of pivot column.
  840     CONTINUE
c
c
c remove pivot column from column oriented storage and update row
c     ordering arrays.
          I1 = IPC(JPIV)
          I2 = IPC(JPIV) + LENC(JPIV) - 1
          NZCOL = NZCOL - LENC(JPIV)
          DO 890 II=I1,I2
            I = IRN(II)
            IRN(II) = 0
            NZ = LENR(I) - LENRL(I)
            IF (NZ.NE.0) GO TO 850
            LASTR(I) = 0
            GO TO 890
  850       IFIR = IFIRST(NZ)
            IFIRST(NZ) = I
            IF (IFIR) 860, 880, 870
  860       LASTR(I) = IFIR
            NEXTR(I) = 0
            GO TO 890
  870       LASTR(I) = LASTR(IFIR)
            NEXTR(I) = IFIR
            LASTR(IFIR) = I
            GO TO 890
  880       LASTR(I) = 0
            NEXTR(I) = 0
            NZMIN = MIN0(NZMIN,NZ)
  890     CONTINUE
c restore iq and nullify u part of old pivot row.
c    record the column permutation in lastc(jpiv) and the row
c    permutation in lastr(ipiv).
  900     IPC(JPIV) = -ISING
          LASTR(IPIV) = PIVOT
          IF (LENPIV.EQ.0) GO TO 980
          NZROW = NZROW - LENPIV
          JVAL = IJP1
          JZER = IPTR(IPIV)
          IPTR(IPIV) = 0
          DO 910 JCOUNT=1,LENPIV
            J = ICN(JVAL)
            IQ(J) = ICN(JZER)
            ICN(JZER) = 0
            JVAL = JVAL + 1
            JZER = JZER + 1
  910     CONTINUE
c adjust column ordering arrays.
          IF (NSRCH.GT.NN) GO TO 920
          DO 916 JJ=IJP1,PIVEND
            J = ICN(JJ)
            NZ = LENC(J)
            IF (NZ.NE.0) GO TO 914
            IPC(J) = 0
            GO TO 916
  914       NZMIN = MIN0(NZMIN,NZ)
  916     CONTINUE
          GO TO 980
  920     JJ = COLUPD
          DO 970 JDUMMY=1,NN
            J = JJ
            IF (J.EQ.NN+1) GO TO 980
            JJ = -NEXTC(J)
            NZ = LENC(J)
            IF (NZ.NE.0) GO TO 924
            IPC(J) = 0
            GO TO 970
  924       IFIR = IFIRST(NZ)
            LASTC(J) = 0
            IF (IFIR) 930, 940, 950
  930       IFIRST(NZ) = -J
            IFIR = -IFIR
            LASTC(IFIR) = J
            NEXTC(J) = IFIR
            GO TO 970
  940       IFIRST(NZ) = -J
            NEXTC(J) = 0
            GO TO 960
  950       LC = -LASTR(IFIR)
            LASTR(IFIR) = -J
            NEXTC(J) = LC
            IF (LC.NE.0) LASTC(LC) = J
  960       NZMIN = MIN0(NZMIN,NZ)
  970     CONTINUE
  980   CONTINUE
c ********************************************
c ****    end of main elimination loop    ****
c ********************************************
c
c reset iactiv to point to the beginning of the next block.
  990   IF (ILAST.NE.NN) IACTIV = IPTR(ILAST+1)
 1000 CONTINUE
c
c ********************************************
c ****    end of deomposition of block    ****
c ********************************************
c
c record singularity (if any) in iq array.
      IF (IRANK.EQ.NN) GO TO 1020
      DO 1010 I=1,NN
        IF (IPC(I).LT.0) GO TO 1010
        ISING = IPC(I)
        IQ(ISING) = -IQ(ISING)
        IPC(I) = -ISING
 1010 CONTINUE
c
c run through lu decomposition changing column indices to that of new
c     order and permuting lenr and lenrl arrays according to pivot
c     permutations.
 1020 ISTART = IDISP(1)
      IEND = IBEG - 1
      IF (IEND.LT.ISTART) GO TO 1040
      DO 1030 JJ=ISTART,IEND
        JOLD = ICN(JJ)
        ICN(JJ) = -IPC(JOLD)
 1030 CONTINUE
 1040 DO 1050 II=1,NN
        I = LASTR(II)
        NEXTR(I) = LENR(II)
        IPTR(I) = LENRL(II)
 1050 CONTINUE
      DO 1060 I=1,NN
        LENRL(I) = IPTR(I)
        LENR(I) = NEXTR(I)
 1060 CONTINUE
c
c update permutation arrays ip and iq.
      DO 1070 II=1,NN
        I = LASTR(II)
        J = -IPC(II)
        NEXTR(I) = IABS(IP(II)+0)
        IPTR(J) = IABS(IQ(II)+0)
 1070 CONTINUE
      DO 1080 I=1,NN
        IF (IP(I).LT.0) NEXTR(I) = -NEXTR(I)
        IP(I) = NEXTR(I)
        IF (IQ(I).LT.0) IPTR(I) = -IPTR(I)
        IQ(I) = IPTR(I)
 1080 CONTINUE
      IP(NN) = IABS(IP(NN)+0)
      IDISP(2) = IEND
      GO TO 1120
c
c   ***    error returns    ***
 1090 IDISP(2) = IACTIV
      IF (LP.EQ.0) GO TO 1120
      WRITE (LP,99996)
      GO TO 1110
 1100 IF (IFLAG.EQ.-5) IFLAG = -6
      IF (IFLAG.NE.-6) IFLAG = -3
      IDISP(2) = IACTIV
      IF (LP.EQ.0) GO TO 1120
      IF (IFLAG.EQ.-3) WRITE (LP,99995)
      IF (IFLAG.EQ.-6) WRITE (LP,99994)
 1110 PIVOT = PIVOT - ISTART + 1
      WRITE (LP,99993) PIVOT, NBLOCK, ISTART, ILAST
      IF (PIVOT.EQ.0) WRITE (LP,99992) MINIRN
c
c
 1120 RETURN
99999 FORMAT (54H ERROR RETURN FROM MA30A/AD BECAUSE MATRIX IS STRUCTUR,
     * 13HALLY SINGULAR)
99998 FORMAT (54H ERROR RETURN FROM MA30A/AD BECAUSE MATRIX IS NUMERICA,
     * 12HLLY SINGULAR)
99997 FORMAT (48H LU DECOMPOSITION DESTROYED TO CREATE MORE SPACE)
99996 FORMAT (54H ERROR RETURN FROM MA30A/AD BECAUSE LICN NOT BIG ENOUG,
     * 1HH)
99995 FORMAT (54H ERROR RETURN FROM MA30A/AD BECAUSE LIRN NOT BIG ENOUG,
     * 1HH)
99994 FORMAT (51H ERROR RETURN FROM MA30A/AD LIRN AND LICN TOO SMALL)
99993 FORMAT (10H AT STAGE , I5, 10H IN BLOCK , I5, 16H WITH FIRST ROW ,
     * I5, 14H AND LAST ROW , I5)
99992 FORMAT (34H TO CONTINUE SET LIRN TO AT LEAST , I8)
      END
c
c=====================================================================
c
      SUBROUTINE MA30BD(N, ICN, A, LICN, LENR, LENRL, IDISP, IP, IQ, W,
     * IW, IFLAG)
c ma30b/bd performs the lu decomposition of the diagonal blocks of a
c     new matrix paq of the same sparsity pattern, using information
c     from a previous call to ma30a/ad. the entries of the input
c     matrix  must already be in their final positions in the lu
c     decomposition structure.  this routine executes about five times
c     faster than ma30a/ad.
c
c we now describe the argument list for ma30b/bd. consult ma30a/ad for
c     further information on these parameters.
c n  is an integer variable set to the order of the matrix.
c icn is an integer array of length licn. it should be unchanged
c     since the last call to ma30a/ad. it is not altered by ma30b/bd.
c a  is a real/double precision array of length licn the user must set
c     entries idisp(1) to idisp(2) to contain the entries in the
c     diagonal blocks of the matrix paq whose column numbers are held
c     in icn, using corresponding positions. note that some zeros may
c     need to be held explicitly. on output entries idisp(1) to
c     idisp(2) of array a contain the lu decomposition of the diagonal
c     blocks of paq. entries a(1) to a(idisp(1)-1) are neither
c     required nor altered by ma30b/bd.
c licn  is an integer variable which must be set by the user to the
c     length of arrays a and icn. it is not altered by ma30b/bd.
c lenr,lenrl are integer arrays of length n. they should be
c     unchanged since the last call to ma30a/ad. they are not altered
c     by ma30b/bd.
c idisp  is an integer array of length 2. it should be unchanged since
c     the last call to ma30a/ad. it is not altered by ma30b/bd.
c ip,iq  are integer arrays of length n. they should be unchanged
c     since the last call to ma30a/ad. they are not altered by
c     ma30b/bd.
c w  is a real/double precision array of length n which is used as
c     workspace by ma30b/bd.
c iw  is an integer array of length n which is used as workspace by
c     ma30b/bd.
c iflag  is an integer variable. on output from ma30b/bd, iflag has
c     the value zero if the factorization was successful, has the
c     value i if pivot i was very small and has the value -i if an
c     unexpected singularity was detected at stage i of the
c     decomposition.
c
      DOUBLE PRECISION A(LICN), W(N), AU, EPS, ROWMAX, ZERO, ONE, RMIN,
     * TOL, BIG
      LOGICAL ABORT1, ABORT2, ABORT3, STAB, LBIG
      INTEGER IW(N), IDISP(2), PIVPOS
      INTEGER ICN(LICN), LENR(N), LENRL(N), IP(N), IQ(N)
c see block data for comments on variables in common.
      COMMON /MA30ED/ LP, ABORT1, ABORT2, ABORT3
      COMMON /MA30ID/ TOL, BIG, NDROP, NSRCH, LBIG
      COMMON /MA30GD/ EPS, RMIN
      DATA ZERO /0.0D0/, ONE /1.0D0/
      STAB = EPS.LE.ONE
      RMIN = EPS
      ISING = 0
      IFLAG = 0
      DO 10 I=1,N
        W(I) = ZERO
   10 CONTINUE
c set up pointers to the beginning of the rows.
      IW(1) = IDISP(1)
      IF (N.EQ.1) GO TO 25
      DO 20 I=2,N
        IW(I) = IW(I-1) + LENR(I-1)
   20 CONTINUE
c
c   ****   start  of main loop    ****
c at step i, row i of a is transformed to row i of l/u by adding
c     appropriate multiples of rows 1 to i-1.
c     .... using row-gauss elimination.
   25 DO 160 I=1,N
c istart is beginning of row i of a and row i of l.
        ISTART = IW(I)
c ifin is end of row i of a and row i of u.
        IFIN = ISTART + LENR(I) - 1
c ilend is end of row i of l.
        ILEND = ISTART + LENRL(I) - 1
        IF (ISTART.GT.ILEND) GO TO 90
c load row i of a into vector w.
        DO 30 JJ=ISTART,IFIN
          J = ICN(JJ)
          W(J) = A(JJ)
   30   CONTINUE
c
c add multiples of appropriate rows of  i to i-1  to row i.
        DO 70 JJ=ISTART,ILEND
          J = ICN(JJ)
c ipivj is position of pivot in row j.
          IPIVJ = IW(J) + LENRL(J)
c form multiplier au.
          AU = -W(J)/A(IPIVJ)
          IF (LBIG) BIG = DMAX1(DABS(AU),BIG)
          W(J) = AU
c au * row j (u part) is added to row i.
          IPIVJ = IPIVJ + 1
          JFIN = IW(J) + LENR(J) - 1
          IF (IPIVJ.GT.JFIN) GO TO 70
c innermost loop.
          IF (LBIG) GO TO 50
          DO 40 JAYJAY=IPIVJ,JFIN
            JAY = ICN(JAYJAY)
            W(JAY) = W(JAY) + AU*A(JAYJAY)
   40     CONTINUE
          GO TO 70
   50     DO 60 JAYJAY=IPIVJ,JFIN
            JAY = ICN(JAYJAY)
            W(JAY) = W(JAY) + AU*A(JAYJAY)
            BIG = DMAX1(DABS(W(JAY)),BIG)
   60     CONTINUE
   70   CONTINUE
c
c reload w back into a (now l/u)
        DO 80 JJ=ISTART,IFIN
          J = ICN(JJ)
          A(JJ) = W(J)
          W(J) = ZERO
   80   CONTINUE
c we now perform the stability checks.
   90   PIVPOS = ILEND + 1
        IF (IQ(I).GT.0) GO TO 140
c matrix had singularity at this point in ma30a/ad.
c is it the first such pivot in current block ?
        IF (ISING.EQ.0) ISING = I
c does current matrix have a singularity in the same place ?
        IF (PIVPOS.GT.IFIN) GO TO 100
        IF (A(PIVPOS).NE.ZERO) GO TO 170
c it does .. so set ising if it is not the end of the current block
c check to see that appropriate part of l/u is zero or null.
  100   IF (ISTART.GT.IFIN) GO TO 120
        DO 110 JJ=ISTART,IFIN
          IF (ICN(JJ).LT.ISING) GO TO 110
          IF (A(JJ).NE.ZERO) GO TO 170
  110   CONTINUE
  120   IF (PIVPOS.LE.IFIN) A(PIVPOS) = ONE
        IF (IP(I).GT.0 .AND. I.NE.N) GO TO 160
c end of current block ... reset zero pivots and ising.
        DO 130 J=ISING,I
          IF ((LENR(J)-LENRL(J)).EQ.0) GO TO 130
          JJ = IW(J) + LENRL(J)
          A(JJ) = ZERO
  130   CONTINUE
        ISING = 0
        GO TO 160
c matrix had non-zero pivot in ma30a/ad at this stage.
  140   IF (PIVPOS.GT.IFIN) GO TO 170
        IF (A(PIVPOS).EQ.ZERO) GO TO 170
        IF (.NOT.STAB) GO TO 160
        ROWMAX = ZERO
        DO 150 JJ=PIVPOS,IFIN
          ROWMAX = DMAX1(ROWMAX,DABS(A(JJ)))
  150   CONTINUE
        IF (DABS(A(PIVPOS))/ROWMAX.GE.RMIN) GO TO 160
        IFLAG = I
        RMIN = DABS(A(PIVPOS))/ROWMAX
c   ****    end of main loop    ****
  160 CONTINUE
c
      GO TO 180
c   ***   error return   ***
  170 IF (LP.NE.0) WRITE (LP,99999) I
      IFLAG = -I
c
  180 RETURN
99999 FORMAT (54H ERROR RETURN FROM MA30B/BD SINGULARITY DETECTED IN RO,
     * 1HW, I8)
      END
c
c======================================================================
c
      SUBROUTINE MA30CD(N, ICN, A, LICN, LENR, LENRL, LENOFF, IDISP, IP,
     * IQ, X, W, MTYPE)
c ma30c/cd uses the factors produced by ma30a/ad or ma30b/bd to solve
c     ax=b or a transpose x=b when the matrix p1*a*q1 (paq) is block
c     lower triangular (including the case of only one diagonal
c     block).
c
c we now describe the argument list for ma30c/cd.
c n  is an integer variable set to the order of the matrix. it is not
c     altered by the subroutine.
c icn is an integer array of length licn. entries idisp(1) to
c     idisp(2) should be unchanged since the last call to ma30a/ad. if
c     the matrix has more than one diagonal block, then column indices
c     corresponding to non-zeros in sub-diagonal blocks of paq must
c     appear in positions 1 to idisp(1)-1. for the same row those
c     entries must be contiguous, with those in row i preceding those
c     in row i+1 (i=1,...,n-1) and no wasted space between rows.
c     entries may be in any order within each row. it is not altered
c     by ma30c/cd.
c a  is a real/double precision array of length licn.  entries
c     idisp(1) to idisp(2) should be unchanged since the last call to
c     ma30a/ad or ma30b/bd.  if the matrix has more than one diagonal
c     block, then the values of the non-zeros in sub-diagonal blocks
c     must be in positions 1 to idisp(1)-1 in the order given by icn.
c     it is not altered by ma30c/cd.
c licn  is an integer variable set to the size of arrays icn and a.
c     it is not altered by ma30c/cd.
c lenr,lenrl are integer arrays of length n which should be
c     unchanged since the last call to ma30a/ad. they are not altered
c     by ma30c/cd.
c lenoff  is an integer array of length n. if the matrix paq (or
c     p1*a*q1) has more than one diagonal block, then lenoff(i),
c     i=1,...,n should be set to the number of non-zeros in row i of
c     the matrix paq which are in sub-diagonal blocks.  if there is
c     only one diagonal block then lenoff(1) may be set to -1, in
c     which case the other entries of lenoff are never accessed. it is
c     not altered by ma30c/cd.
c idisp  is an integer array of length 2 which should be unchanged
c     since the last call to ma30a/ad. it is not altered by ma30c/cd.
c ip,iq are integer arrays of length n which should be unchanged
c     since the last call to ma30a/ad. they are not altered by
c     ma30c/cd.
c x is a real/double precision array of length n. it must be set by
c     the user to the values of the right hand side vector b for the
c     equations being solved.  on exit from ma30c/cd it will be equal
c     to the solution x required.
c w  is a real/double precision array of length n which is used as
c     workspace by ma30c/cd.
c mtype is an integer variable which must be set by the user. if
c     mtype=1, then the solution to the system ax=b is returned; any
c     other value for mtype will return the solution to the system a
c     transpose x=b. it is not altered by ma30c/cd.
c
      DOUBLE PRECISION A(LICN), X(N), W(N), WII, WI, RESID, ZERO
      LOGICAL NEG, NOBLOC
      INTEGER IDISP(2)
      INTEGER ICN(LICN), LENR(N), LENRL(N), LENOFF(N), IP(N), IQ(N)
c see block data for comments on variables in common.
      COMMON /MA30HD/ RESID
      DATA ZERO /0.0D0/
c
c the final value of resid is the maximum residual for an inconsistent
c     set of equations.
      RESID = ZERO
c nobloc is .true. if subroutine block has been used previously and
c     is .false. otherwise.  the value .false. means that lenoff
c     will not be subsequently accessed.
      NOBLOC = LENOFF(1).LT.0
      IF (MTYPE.NE.1) GO TO 140
c
c we now solve   a * x = b.
c neg is used to indicate when the last row in a block has been
c     reached.  it is then set to true whereafter backsubstitution is
c     performed on the block.
      NEG = .FALSE.
c ip(n) is negated so that the last row of the last block can be
c     recognised.  it is reset to its positive value on exit.
      IP(N) = -IP(N)
c preorder vector ... w(i) = x(ip(i))
      DO 10 II=1,N
        I = IP(II)
        I = IABS(I)
        W(II) = X(I)
   10 CONTINUE
c lt holds the position of the first non-zero in the current row of the
c     off-diagonal blocks.
      LT = 1
c ifirst holds the index of the first row in the current block.
      IFIRST = 1
c iblock holds the position of the first non-zero in the current row
c     of the lu decomposition of the diagonal blocks.
      IBLOCK = IDISP(1)
c if i is not the last row of a block, then a pass through this loop
c     adds the inner product of row i of the off-diagonal blocks and w
c     to w and performs forward elimination using row i of the lu
c     decomposition.   if i is the last row of a block then, after
c     performing these aforementioned operations, backsubstitution is
c     performed using the rows of the block.
      DO 120 I=1,N
        WI = W(I)
        IF (NOBLOC) GO TO 30
        IF (LENOFF(I).EQ.0) GO TO 30
c operations using lower triangular blocks.
c ltend is the end of row i in the off-diagonal blocks.
        LTEND = LT + LENOFF(I) - 1
        DO 20 JJ=LT,LTEND
          J = ICN(JJ)
          WI = WI - A(JJ)*W(J)
   20   CONTINUE
c lt is set the beginning of the next off-diagonal row.
        LT = LTEND + 1
c set neg to .true. if we are on the last row of the block.
   30   IF (IP(I).LT.0) NEG = .TRUE.
        IF (LENRL(I).EQ.0) GO TO 50
c forward elimination phase.
c iend is the end of the l part of row i in the lu decomposition.
        IEND = IBLOCK + LENRL(I) - 1
        DO 40 JJ=IBLOCK,IEND
          J = ICN(JJ)
          WI = WI + A(JJ)*W(J)
   40   CONTINUE
c iblock is adjusted to point to the start of the next row.
   50   IBLOCK = IBLOCK + LENR(I)
        W(I) = WI
        IF (.NOT.NEG) GO TO 120
c back substitution phase.
c j1 is position in a/icn after end of block beginning in row ifirst
c     and ending in row i.
        J1 = IBLOCK
c are there any singularities in this block?  if not, continue with
c     the backsubstitution.
        IB = I
        IF (IQ(I).GT.0) GO TO 70
        DO 60 III=IFIRST,I
          IB = I - III + IFIRST
          IF (IQ(IB).GT.0) GO TO 70
          J1 = J1 - LENR(IB)
          RESID = DMAX1(RESID,DABS(W(IB)))
          W(IB) = ZERO
   60   CONTINUE
c entire block is singular.
        GO TO 110
c each pass through this loop performs the back-substitution
c     operations for a single row, starting at the end of the block and
c     working through it in reverse order.
   70   DO 100 III=IFIRST,IB
          II = IB - III + IFIRST
c j2 is end of row ii.
          J2 = J1 - 1
c j1 is beginning of row ii.
          J1 = J1 - LENR(II)
c jpiv is the position of the pivot in row ii.
          JPIV = J1 + LENRL(II)
          JPIVP1 = JPIV + 1
c jump if row  ii of u has no non-zeros.
          IF (J2.LT.JPIVP1) GO TO 90
          WII = W(II)
          DO 80 JJ=JPIVP1,J2
            J = ICN(JJ)
            WII = WII - A(JJ)*W(J)
   80     CONTINUE
          W(II) = WII
   90     W(II) = W(II)/A(JPIV)
  100   CONTINUE
  110   IFIRST = I + 1
        NEG = .FALSE.
  120 CONTINUE
c
c reorder solution vector ... x(i) = w(iqinverse(i))
      DO 130 II=1,N
        I = IQ(II)
        I = IABS(I)
        X(I) = W(II)
  130 CONTINUE
      IP(N) = -IP(N)
      GO TO 320
c
c
c we now solve   atranspose * x = b.
c preorder vector ... w(i)=x(iq(i))
  140 DO 150 II=1,N
        I = IQ(II)
        I = IABS(I)
        W(II) = X(I)
  150 CONTINUE
c lj1 points to the beginning the current row in the off-diagonal
c     blocks.
      LJ1 = IDISP(1)
c iblock is initialized to point to the beginning of the block after
c     the last one ]
      IBLOCK = IDISP(2) + 1
c ilast is the last row in the current block.
      ILAST = N
c iblend points to the position after the last non-zero in the
c     current block.
      IBLEND = IBLOCK
c each pass through this loop operates with one diagonal block and
c     the off-diagonal part of the matrix corresponding to the rows
c     of this block.  the blocks are taken in reverse order and the
c     number of times the loop is entered is min(n,no. blocks+1).
      DO 290 NUMBLK=1,N
        IF (ILAST.EQ.0) GO TO 300
        IBLOCK = IBLOCK - LENR(ILAST)
c this loop finds the index of the first row in the current block..
c     it is first and iblock is set to the position of the beginning
c     of this first row.
        DO 160 K=1,N
          II = ILAST - K
          IF (II.EQ.0) GO TO 170
          IF (IP(II).LT.0) GO TO 170
          IBLOCK = IBLOCK - LENR(II)
  160   CONTINUE
  170   IFIRST = II + 1
c j1 points to the position of the beginning of row i (lt part) or pivot
        J1 = IBLOCK
c forward elimination.
c each pass through this loop performs the operations for one row of the
c     block.  if the corresponding entry of w is zero then the
c     operations can be avoided.
        DO 210 I=IFIRST,ILAST
          IF (W(I).EQ.ZERO) GO TO 200
c jump if row i singular.
          IF (IQ(I).LT.0) GO TO 220
c j2 first points to the pivot in row i and then is made to point to the
c     first non-zero in the u transpose part of the row.
          J2 = J1 + LENRL(I)
          WI = W(I)/A(J2)
          IF (LENR(I)-LENRL(I).EQ.1) GO TO 190
          J2 = J2 + 1
c j3 points to the end of row i.
          J3 = J1 + LENR(I) - 1
          DO 180 JJ=J2,J3
            J = ICN(JJ)
            W(J) = W(J) - A(JJ)*WI
  180     CONTINUE
  190     W(I) = WI
  200     J1 = J1 + LENR(I)
  210   CONTINUE
        GO TO 240
c deals with rest of block which is singular.
  220   DO 230 II=I,ILAST
          RESID = DMAX1(RESID,DABS(W(II)))
          W(II) = ZERO
  230   CONTINUE
c back substitution.
c this loop does the back substitution on the rows of the block in
c     the reverse order doing it simultaneously on the l transpose part
c     of the diagonal blocks and the off-diagonal blocks.
  240   J1 = IBLEND
        DO 280 IBACK=IFIRST,ILAST
          I = ILAST - IBACK + IFIRST
c j1 points to the beginning of row i.
          J1 = J1 - LENR(I)
          IF (LENRL(I).EQ.0) GO TO 260
c j2 points to the end of the l transpose part of row i.
          J2 = J1 + LENRL(I) - 1
          DO 250 JJ=J1,J2
            J = ICN(JJ)
            W(J) = W(J) + A(JJ)*W(I)
  250     CONTINUE
  260     IF (NOBLOC) GO TO 280
c operations using lower triangular blocks.
          IF (LENOFF(I).EQ.0) GO TO 280
c lj2 points to the end of row i of the off-diagonal blocks.
          LJ2 = LJ1 - 1
c lj1 points to the beginning of row i of the off-diagonal blocks.
          LJ1 = LJ1 - LENOFF(I)
          DO 270 JJ=LJ1,LJ2
            J = ICN(JJ)
            W(J) = W(J) - A(JJ)*W(I)
  270     CONTINUE
  280   CONTINUE
        IBLEND = J1
        ILAST = IFIRST - 1
  290 CONTINUE
c reorder solution vector ... x(i)=w(ipinverse(i))
  300 DO 310 II=1,N
        I = IP(II)
        I = IABS(I)
        X(I) = W(II)
  310 CONTINUE
c
  320 RETURN
      END
      SUBROUTINE MA30DD(A, ICN, IPTR, N, IACTIV, ITOP, REALS)
c this subroutine performs garbage collection operations on the
c     arrays a, icn and irn.
c iactiv is the first position in arrays a/icn from which the compress
c     starts.  on exit, iactiv equals the position of the first entry
c     in the compressed part of a/icn
c
      DOUBLE PRECISION A(ITOP)
      LOGICAL REALS
      INTEGER IPTR(N)
      INTEGER ICN(ITOP)
c see block data for comments on variables in common.
      COMMON /MA30FD/ IRNCP, ICNCP, IRANK, MINIRN, MINICN
c
      IF (REALS) ICNCP = ICNCP + 1
      IF (.NOT.REALS) IRNCP = IRNCP + 1
c set the first non-zero entry in each row to the negative of the
c     row/col number and hold this row/col index in the row/col
c     pointer.  this is so that the beginning of each row/col can
c     be recognized in the subsequent scan.
      DO 10 J=1,N
        K = IPTR(J)
        IF (K.LT.IACTIV) GO TO 10
        IPTR(J) = ICN(K)
        ICN(K) = -J
   10 CONTINUE
      KN = ITOP + 1
      KL = ITOP - IACTIV + 1
c go through arrays in reverse order compressing to the back so
c     that there are no zeros held in positions iactiv to itop in icn.
c     reset first entry of each row/col and pointer array iptr.
      DO 30 K=1,KL
        JPOS = ITOP - K + 1
        IF (ICN(JPOS).EQ.0) GO TO 30
        KN = KN - 1
        IF (REALS) A(KN) = A(JPOS)
        IF (ICN(JPOS).GE.0) GO TO 20
c first non-zero of row/col has been located
        J = -ICN(JPOS)
        ICN(JPOS) = IPTR(J)
        IPTR(J) = KN
   20   ICN(KN) = ICN(JPOS)
   30 CONTINUE
      IACTIV = KN
      RETURN
      END
c######date   01 jan 1984     copyright ukaea, harwell.
c######alias mc13d
      SUBROUTINE MC13D(N,ICN,LICN,IP,LENR,IOR,IB,NUM,IW)
      INTEGER IP(N)
      INTEGER ICN(LICN),LENR(N),IOR(N),IB(N),IW(N,3)
      CALL MC13E(N,ICN,LICN,IP,LENR,IOR,IB,NUM,IW(1,1),IW(1,2),IW(1,3))
      RETURN
      END
      SUBROUTINE MC13E(N,ICN,LICN,IP,LENR,ARP,IB,NUM,LOWL,NUMB,PREV)
      INTEGER STP,DUMMY
      INTEGER IP(N)
c
c arp(i) is one less than the number of unsearched edges leaving
c     node i.  at the end of the algorithm it is set to a
c     permutation which puts the matrix in block lower
c     triangular form.
c ib(i) is the position in the ordering of the start of the ith
c     block.  ib(n+1-i) holds the node number of the ith node
c     on the stack.
c lowl(i) is the smallest stack position of any node to which a path
c     from node i has been found.  it is set to n+1 when node i
c     is removed from the stack.
c numb(i) is the position of node i in the stack if it is on
c     it, is the permuted order of node i for those nodes
c     whose final position has been found and is otherwise zero.
c prev(i) is the node at the end of the path when node i was
c     placed on the stack.
      INTEGER ICN(LICN),LENR(N),ARP(N),IB(N),LOWL(N),NUMB(N),
     1PREV(N)
c
c
c   icnt is the number of nodes whose positions in final ordering have
c     been found.
      ICNT=0
c num is the number of blocks that have been found.
      NUM=0
      NNM1=N+N-1
c
c initialization of arrays.
      DO 20 J=1,N
      NUMB(J)=0
      ARP(J)=LENR(J)-1
   20 CONTINUE
c
c
      DO 120 ISN=1,N
c look for a starting node
      IF (NUMB(ISN).NE.0) GO TO 120
      IV=ISN
c ist is the number of nodes on the stack ... it is the stack pointer.
      IST=1
c put node iv at beginning of stack.
      LOWL(IV)=1
      NUMB(IV)=1
      IB(N)=IV
c
c the body of this loop puts a new node on the stack or backtracks.
      DO 110 DUMMY=1,NNM1
      I1=ARP(IV)
c have all edges leaving node iv been searched.
      IF (I1.LT.0) GO TO 60
      I2=IP(IV)+LENR(IV)-1
      I1=I2-I1
c
c look at edges leaving node iv until one enters a new node or
c     all edges are exhausted.
      DO 50 II=I1,I2
      IW=ICN(II)
c has node iw been on stack already.
      IF (NUMB(IW).EQ.0) GO TO 100
c update value of lowl(iv) if necessary.
  50  LOWL(IV)=MIN0(LOWL(IV),LOWL(IW))
c
c there are no more edges leaving node iv.
      ARP(IV)=-1
c is node iv the root of a block.
   60 IF (LOWL(IV).LT.NUMB(IV)) GO TO 90
c
c order nodes in a block.
      NUM=NUM+1
      IST1=N+1-IST
      LCNT=ICNT+1
c peel block off the top of the stack starting at the top and
c     working down to the root of the block.
      DO 70 STP=IST1,N
      IW=IB(STP)
      LOWL(IW)=N+1
      ICNT=ICNT+1
      NUMB(IW)=ICNT
      IF (IW.EQ.IV) GO TO 80
   70 CONTINUE
   80 IST=N-STP
      IB(NUM)=LCNT
c are there any nodes left on the stack.
      IF (IST.NE.0) GO TO 90
c have all the nodes been ordered.
      IF (ICNT.LT.N) GO TO 120
      GO TO 130
c
c backtrack to previous node on path.
   90 IW=IV
      IV=PREV(IV)
c update value of lowl(iv) if necessary.
      LOWL(IV)=MIN0(LOWL(IV),LOWL(IW))
      GO TO 110
c
c put new node on the stack.
 100  ARP(IV)=I2-II-1
      PREV(IW)=IV
      IV=IW
      IST=IST+1
      LOWL(IV)=IST
      NUMB(IV)=IST
      K=N+1-IST
      IB(K)=IV
  110 CONTINUE
c
  120 CONTINUE
c
c
c put permutation in the required form.
  130 DO 140 I=1,N
      II=NUMB(I)
 140  ARP(II)=I
      RETURN
      END
c######date   01 jan 1984     copyright ukaea, harwell.
c######alias mc20ad mc20bd
      SUBROUTINE MC20AD(NC,MAXA,A,INUM,JPTR,JNUM,JDISP)
c
      INTEGER   INUM(MAXA),JNUM(MAXA)
      DOUBLE PRECISION A(MAXA),ACE,ACEP
      DIMENSION JPTR(NC)
c
c     ******************************************************************
c
      NULL=-JDISP
c**      clear jptr
      DO 60 J=1,NC
   60 JPTR(J)=0
c**      count the number of elements in each column.
      DO 120 K=1,MAXA
      J=JNUM(K)+JDISP
      JPTR(J)=JPTR(J)+1
  120 CONTINUE
c**      set the jptr array
      K=1
      DO 150 J=1,NC
      KR=K+JPTR(J)
      JPTR(J)=K
  150 K=KR
c
c**      reorder the elements into column order.  the algorithm is an
c        in-place sort and is of order maxa.
      DO 230 I=1,MAXA
c        establish the current entry.
      JCE=JNUM(I)+JDISP
      IF(JCE.EQ.0) GO TO 230
      ACE=A(I)
      ICE=INUM(I)
c        clear the location vacated.
      JNUM(I)=NULL
c        chain from current entry to store items.
      DO 200 J=1,MAXA
c        current entry not in correct position.  determine correct
c        position to store entry.
      LOC=JPTR(JCE)
      JPTR(JCE)=JPTR(JCE)+1
c        save contents of that location.
      ACEP=A(LOC)
      ICEP=INUM(LOC)
      JCEP=JNUM(LOC)
c        store current entry.
      A(LOC)=ACE
      INUM(LOC)=ICE
      JNUM(LOC)=NULL
c        check if next current entry needs to be processed.
      IF(JCEP.EQ.NULL) GO TO 230
c        it does.  copy into current entry.
      ACE=ACEP
      ICE=ICEP
  200 JCE=JCEP+JDISP
c
  230 CONTINUE
c
c**      reset jptr vector.
      JA=1
      DO 250 J=1,NC
      JB=JPTR(J)
      JPTR(J)=JA
  250 JA=JB
      RETURN
      END
      SUBROUTINE MC20BD(NC,MAXA,A,INUM,JPTR)
C     REAL A(MAXA)                                                   JS/
      DOUBLE PRECISION A(MAXA),ACE
      INTEGER   INUM(MAXA)
C     INTEGER*2 INUM(MAXA)                                           DS/
      DIMENSION JPTR(NC)
C
C     ******************************************************************
C
      KMAX=MAXA
      DO 30 JJ=1,NC
      J=NC+1-JJ
      KLO=JPTR(J)+1
      IF(KLO.GT.KMAX)GO TO 30
      KOR=KMAX
      DO 25 KDUMMY=KLO,KMAX
C ITEMS KOR, KOR+1, .... ,KMAX ARE IN ORDER
      ACE=A(KOR-1)
      ICE=INUM(KOR-1)
      DO 10 K=KOR,KMAX
      IK=INUM(K)
      IF(IABS(ICE).LE.IABS(IK))GO TO 20
      INUM(K-1)=IK
10    A(K-1)=A(K)
      K=KMAX+1
20    INUM(K-1)=ICE
      A(K-1)=ACE
25    KOR=KOR-1
C        NEXT COLUMN
30    KMAX=KLO-2
      RETURN
      END
c######date   01 jan 1984     copyright ukaea, harwell.
c######alias mc21a
      SUBROUTINE MC21A(N,ICN,LICN,IP,LENR,IPERM,NUMNZ,IW)
      INTEGER IP(N)
      INTEGER ICN(LICN),LENR(N),IPERM(N),IW(N,4)
      CALL MC21B(N,ICN,LICN,IP,LENR,IPERM,NUMNZ,IW(1,1),IW(1,2),IW(1,3),
     1IW(1,4))
      RETURN
      END
      SUBROUTINE MC21B(N,ICN,LICN,IP,LENR,IPERM,NUMNZ,PR,ARP,CV,OUT)
      INTEGER IP(N)
c   pr(i) is the previous row to i in the depth first search.
c it is used as a work array in the sorting algorithm.
c   elements (iperm(i),i) i=1, ... n  are non-zero at the end of the
c algorithm unless n assignments have not been made.  in which case
c (iperm(i),i) will be zero for n-numnz entries.
c   cv(i) is the most recent row extension at which column i
c was visited.
c   arp(i) is one less than the number of non-zeros in row i
c which have not been scanned when looking for a cheap assignment.
c   out(i) is one less than the number of non-zeros in row i
c which have not been scanned during one pass through the main loop.
      INTEGER ICN(LICN),LENR(N),IPERM(N),PR(N),CV(N),
     1ARP(N),OUT(N)
c
c   initialization of arrays.
      DO 10 I=1,N
      ARP(I)=LENR(I)-1
      CV(I)=0
   10 IPERM(I)=0
      NUMNZ=0
c
c
c   main loop.
c   each pass round this loop either results in a new assignment
c or gives a row with no assignment.
      DO 130 JORD=1,N
      J=JORD
      PR(J)=-1
      DO 100 K=1,JORD
c look for a cheap assignment
      IN1=ARP(J)
      IF (IN1.LT.0) GO TO 60
      IN2=IP(J)+LENR(J)-1
      IN1=IN2-IN1
      DO 50 II=IN1,IN2
      I=ICN(II)
      IF (IPERM(I).EQ.0) GO TO 110
   50 CONTINUE
c   no cheap assignment in row.
      ARP(J)=-1
c   begin looking for assignment chain starting with row j.
   60 OUT(J)=LENR(J)-1
c inner loop.  extends chain by one or backtracks.
      DO 90 KK=1,JORD
      IN1=OUT(J)
      IF (IN1.LT.0) GO TO 80
      IN2=IP(J)+LENR(J)-1
      IN1=IN2-IN1
c forward scan.
      DO 70 II=IN1,IN2
      I=ICN(II)
      IF (CV(I).EQ.JORD) GO TO 70
c   column i has not yet been accessed during this pass.
      J1=J
      J=IPERM(I)
      CV(I)=JORD
      PR(J)=J1
      OUT(J1)=IN2-II-1
      GO TO 100
   70 CONTINUE
c
c   backtracking step.
   80 J=PR(J)
      IF (J.EQ.-1) GO TO 130
   90 CONTINUE
c
  100 CONTINUE
c
c   new assignment is made.
  110 IPERM(I)=J
      ARP(J)=IN2-II-1
      NUMNZ=NUMNZ+1
      DO 120 K=1,JORD
      J=PR(J)
      IF (J.EQ.-1) GO TO 130
      II=IP(J)+LENR(J)-OUT(J)-2
      I=ICN(II)
      IPERM(I)=J
  120 CONTINUE
c
  130 CONTINUE
c
c   if matrix is structurally singular, we now complete the
c permutation iperm.
      IF (NUMNZ.EQ.N) RETURN
      DO 140 I=1,N
  140 ARP(I)=0
      K=0
      DO 160 I=1,N
      IF (IPERM(I).NE.0) GO TO 150
      K=K+1
      OUT(K)=I
      GO TO 160
  150 J=IPERM(I)
      ARP(J)=I
  160 CONTINUE
      K=0
      DO 170 I=1,N
      IF (ARP(I).NE.0) GO TO 170
      K=K+1
      IOUTK=OUT(K)
      IPERM(IOUTK)=I
  170 CONTINUE
      RETURN
      END
c######date   01 jan 1984     copyright ukaea, harwell.
c######alias mc22ad
      SUBROUTINE MC22AD(N,ICN,A,NZ,LENROW,IP,IQ,IW,IW1)
      DOUBLE PRECISION A(NZ),AVAL
      INTEGER IW(N,2)
      INTEGER   ICN(NZ),LENROW(N),IP(N),IQ(N),IW1(NZ)
      IF (NZ.LE.0) GO TO 1000
      IF (N.LE.0) GO TO 1000
c set start of row i in iw(i,1) and lenrow(i) in iw(i,2)
      IW(1,1)=1
      IW(1,2)=LENROW(1)
      DO 10 I=2,N
      IW(I,1)=IW(I-1,1)+LENROW(I-1)
 10   IW(I,2)=LENROW(I)
c permute lenrow according to ip.  set off-sets for new position
c     of row iold in iw(iold,1) and put old row indices in iw1 in
c     positions corresponding to the new position of this row in a/icn.
      JJ=1
      DO 20 I=1,N
      IOLD=IP(I)
      IOLD=IABS(IOLD)
      LENGTH=IW(IOLD,2)
      LENROW(I)=LENGTH
      IF (LENGTH.EQ.0) GO TO 20
      IW(IOLD,1)=IW(IOLD,1)-JJ
      J2=JJ+LENGTH-1
      DO 15 J=JJ,J2
 15   IW1(J)=IOLD
      JJ=J2+1
 20   CONTINUE
c set inverse permutation to iq in iw(.,2).
      DO 30 I=1,N
      IOLD=IQ(I)
      IOLD=IABS(IOLD)
 30   IW(IOLD,2)=I
c permute a and icn in place, changing to new column numbers.
c
c ***   main loop   ***
c each pass through this loop places a closed chain of column indices
c     in their new (and final) positions ... this is recorded by
c     setting the iw1 entry to zero so that any which are subsequently
c     encountered during this major scan can be bypassed.
      DO 200 I=1,NZ
      IOLD=IW1(I)
      IF (IOLD.EQ.0) GO TO 200
      IPOS=I
      JVAL=ICN(I)
c if row iold is in same positions after permutation go to 150.
      IF (IW(IOLD,1).EQ.0) GO TO 150
      AVAL=A(I)
c **  chain loop  **
c each pass through this loop places one (permuted) column index
c     in its final position  .. viz. ipos.
      DO 100 ICHAIN=1,NZ
c newpos is the original position in a/icn of the element to be placed
c in position ipos.  it is also the position of the next element in
c     the chain.
      NEWPOS=IPOS+IW(IOLD,1)
c is chain complete ?
      IF (NEWPOS.EQ.I) GO TO 130
      A(IPOS)=A(NEWPOS)
      JNUM=ICN(NEWPOS)
      ICN(IPOS)=IW(JNUM,2)
      IPOS=NEWPOS
      IOLD=IW1(IPOS)
      IW1(IPOS)=0
c **  end of chain loop  **
 100  CONTINUE
 130  A(IPOS)=AVAL
 150  ICN(IPOS)=IW(JVAL,2)
c ***   end of main loop   ***
 200  CONTINUE
c
 1000 RETURN
      END
c######date   01 jan 1984     copyright ukaea, harwell.
c######alias mc23ad
c###### calls   mc13    mc21
      SUBROUTINE MC23AD(N,ICN,A,LICN,LENR,IDISP,IP,IQ,LENOFF,IW,IW1)
      DOUBLE PRECISION A(LICN)
      INTEGER IDISP(2),IW1(N,2)
      LOGICAL ABORT
      INTEGER   ICN(LICN),LENR(N),IP(N),IQ(N),LENOFF(N),IW(N,5)
      COMMON /MC23BD/ LP,NUMNZ,NUM,LARGE,ABORT
c input ... n,icn .. a,icn,lenr ....
c
c set up pointers iw(.,1) to the beginning of the rows and set lenoff
c     equal to lenr.
      IW1(1,1)=1
      LENOFF(1)=LENR(1)
      IF (N.EQ.1) GO TO 20
      DO 10 I=2,N
      LENOFF(I)=LENR(I)
   10 IW1(I,1)=IW1(I-1,1)+LENR(I-1)
c idisp(1) points to the first position in a/icn after the
c     off-diagonal blocks and untreated rows.
   20 IDISP(1)=IW1(N,1)+LENR(N)
c
c find row permutation ip to make diagonal zero-free.
      CALL MC21A(N,ICN,LICN,IW1,LENR,IP,NUMNZ,IW)
c
c possible error return for structurally singular matrices.
      IF (NUMNZ.NE.N.AND.ABORT) GO TO 170
c
c iw1(.,2) and lenr are permutations of iw1(.,1) and lenr/lenoff
c     suitable for entry
c     to mc13d since matrix with these row pointer and length arrays
c     has maximum number of non-zeros on the diagonal.
      DO 30 II=1,N
      I=IP(II)
      IW1(II,2)=IW1(I,1)
   30 LENR(II)=LENOFF(I)
c
c find symmetric permutation iq to block lower triangular form.
      CALL MC13D(N,ICN,LICN,IW1(1,2),LENR,IQ,IW(1,4),NUM,IW)
c
      IF (NUM.NE.1) GO TO 60
c
c action taken if matrix is irreducible.
c whole matrix is just moved to the end of the storage.
      DO 40 I=1,N
      LENR(I)=LENOFF(I)
      IP(I)=I
   40 IQ(I)=I
      LENOFF(1)=-1
c idisp(1) is the first position after the last element in the
c     off-diagonal blocks and untreated rows.
      NZ=IDISP(1)-1
      IDISP(1)=1
c idisp(2) is the position in a/icn of the first element in the
c     diagonal blocks.
      IDISP(2)=LICN-NZ+1
      LARGE=N
      IF (NZ.EQ.LICN) GO TO 230
      DO 50 K=1,NZ
      J=NZ-K+1
      JJ=LICN-K+1
      A(JJ)=A(J)
   50 ICN(JJ)=ICN(J)
c 230 = return
      GO TO 230
c
c data structure reordered.
c
c form composite row permutation ... ip(i) = ip(iq(i)).
   60 DO 70 II=1,N
      I=IQ(II)
   70 IW(II,1)=IP(I)
      DO 80 I=1,N
   80 IP(I)=IW(I,1)
c
c run through blocks in reverse order separating diagonal blocks
c     which are moved to the end of the storage.  elements in
c     off-diagonal blocks are left in place unless a compress is
c     necessary.
c
c ibeg indicates the lowest value of j for which icn(j) has been
c     set to zero when element in position j was moved to the
c     diagonal block part of storage.
      IBEG=LICN+1
c iend is the position of the first element of those treated rows
c     which are in diagonal blocks.
      IEND=LICN+1
c large is the dimension of the largest block encountered so far.
      LARGE=0
c
c num is the number of diagonal blocks.
      DO 150 K=1,NUM
      IBLOCK=NUM-K+1
c i1 is first row (in permuted form) of block iblock.
c i2 is last row (in permuted form) of block iblock.
      I1=IW(IBLOCK,4)
      I2=N
      IF (K.NE.1) I2=IW(IBLOCK+1,4)-1
      LARGE=MAX0(LARGE,I2-I1+1)
c go through the rows of block iblock in the reverse order.
      DO 140 II=I1,I2
      INEW=I2-II+I1
c we now deal with row inew in permuted form (row iold in original
c     matrix).
      IOLD=IP(INEW)
c if there is space to move up diagonal block portion of row go to 110
      IF (IEND-IDISP(1).GE.LENOFF(IOLD)) GO TO 110
c
c in-line compress.
c moves separated off-diagonal elements and untreated rows to
c     front of storage.
      JNPOS=IBEG
      ILEND=IDISP(1)-1
      IF (ILEND.LT.IBEG) GO TO 190
      DO 90 J=IBEG,ILEND
      IF (ICN(J).EQ.0) GO TO 90
      ICN(JNPOS)=ICN(J)
      A(JNPOS)=A(J)
      JNPOS=JNPOS+1
   90 CONTINUE
      IDISP(1)=JNPOS
      IF (IEND-JNPOS.LT.LENOFF(IOLD)) GO TO 190
      IBEG=LICN+1
c reset pointers to the beginning of the rows.
      DO 100 I=2,N
  100 IW1(I,1)=IW1(I-1,1)+LENOFF(I-1)
c
c row iold is now split into diag. and off-diag. parts.
  110 IROWB=IW1(IOLD,1)
      LENI=0
      IROWE=IROWB+LENOFF(IOLD)-1
c backward scan of whole of row iold (in original matrix).
      IF (IROWE.LT.IROWB) GO TO 130
      DO 120 JJ=IROWB,IROWE
      J=IROWE-JJ+IROWB
      JOLD=ICN(J)
c iw(.,2) holds the inverse permutation to iq.
c     ..... it was set to this in mc13d.
      JNEW=IW(JOLD,2)
c if (jnew.lt.i1) then ....
c element is in off-diagonal block and so is left in situ.
      IF (JNEW.LT.I1) GO TO 120
c element is in diagonal block and is moved to the end of the storage.
      IEND=IEND-1
      A(IEND)=A(J)
      ICN(IEND)=JNEW
      IBEG=MIN0(IBEG,J)
      ICN(J)=0
      LENI=LENI+1
  120 CONTINUE
c
      LENOFF(IOLD)=LENOFF(IOLD)-LENI
  130 LENR(INEW)=LENI
  140 CONTINUE
c
      IP(I2)=-IP(I2)
  150 CONTINUE
c resets ip(n) to positive value.
      IP(N)=-IP(N)
c idisp(2) is position of first element in diagonal blocks.
      IDISP(2)=IEND
c
c this compress is used to move all off-diagonal elements to the
c     front of the storage.
      IF (IBEG.GT.LICN) GO TO 230
      JNPOS=IBEG
      ILEND=IDISP(1)-1
      DO 160 J=IBEG,ILEND
      IF (ICN(J).EQ.0) GO TO 160
      ICN(JNPOS)=ICN(J)
      A(JNPOS)=A(J)
      JNPOS=JNPOS+1
  160 CONTINUE
c idisp(1) is first position after last element of off-diagonal blocks.
      IDISP(1)=JNPOS
      GO TO 230
c
c
c error return
  170 IF (LP.NE.0) WRITE(LP,180) NUMNZ
  180 FORMAT(33X,41H MATRIX IS STRUCTURALLY SINGULAR, RANK = ,I6)
      IDISP(1)=-1
      GO TO 210
  190 IF (LP.NE.0) WRITE(LP,200) N
  200 FORMAT(33X,33H LICN NOT BIG ENOUGH INCREASE BY ,I6)
      IDISP(1)=-2
  210 IF (LP.NE.0) WRITE(LP,220)
  220 FORMAT(33H+ERROR RETURN FROM MC23AD BECAUSE)
c
  230 RETURN
      END
c######date   01 jan 1984     copyright ukaea, harwell.
c######alias mc24ad
      SUBROUTINE MC24AD(N,ICN,A,LICN,LENR,LENRL,W)
      DOUBLE PRECISION A(LICN),W(N),AMAXL,WROWL,AMAXU,ZERO
      INTEGER   ICN(LICN),LENR(N),LENRL(N)
      DATA ZERO/0.0D0/
      AMAXL=ZERO
      DO 10 I=1,N
 10   W(I)=ZERO
      J0=1
      DO 100 I=1,N
      IF (LENR(I).EQ.0) GO TO 100
      J2=J0+LENR(I)-1
      IF (LENRL(I).EQ.0) GO TO 50
c calculation of 1-norm of l.
      J1=J0+LENRL(I)-1
      WROWL=ZERO
      DO 30 JJ=J0,J1
 30   WROWL=WROWL+DABS(A(JJ))
c amaxl is the maximum norm of columns of l so far found.
      AMAXL=DMAX1(AMAXL,WROWL)
      J0=J1+1
c calculation of norms of columns of u (max-norms).
 50   J0=J0+1
      IF (J0.GT.J2) GO TO 90
      DO 80 JJ=J0,J2
      J=ICN(JJ)
 80   W(J)=DMAX1(DABS(A(JJ)),W(J))
 90   J0=J2+1
 100  CONTINUE
c amaxu is set to maximum max-norm of columns of u.
      AMAXU=ZERO
      DO 200 I=1,N
 200  AMAXU=DMAX1(AMAXU,W(I))
c grofac is max u max-norm times max l 1-norm.
      W(1)=AMAXL*AMAXU
      RETURN
      END
      BLOCK DATA MABLD1
c
c comments on all the common block variables are given here even
c     though some are not initialized by block data.
c lp,mp are used by the subroutine as the unit numbers for its warning
c     and diagnostic messages. default value for both is 6 (for line
c     printer output). the user can either reset them to a different
c     stream number or suppress the output by setting them to zero.
c     while lp directs the output of error diagnostics from the
c     principal subroutines and internally called subroutines, mp
c     controls only the output of a message which warns the user that he
c     has input two or more non-zeros a(i), . . ,a(k) with the same row
c     and column indices.  the action taken in this case is to proceed
c     using a numerical value of a(i)+...+a(k). in the absence of other
c     errors, iflag will equal -14 on exit.
c lblock is a logical variable which controls an option of first
c     preordering the matrix to block lower triangular form (using
c     harwell subroutine mc23a). the preordering is performed if lblock
c     is equal to its default value of .true. if lblock is set to
c     .false. , the option is not invoked and the space allocated to
c     ikeep can be reduced to 4*n+1.
c grow is a logical variable. if it is left at its default value of
c     .true. , then on return from ma28a/ad or ma28b/bd, w(1) will give
c     an estimate (an upper bound) of the increase in size of elements
c     encountered during the decomposition. if the matrix is well
c     scaled, then a high value for w(1), relative to the largest entry
c     in the input matrix, indicates that the lu decomposition may be
c     inaccurate and the user should be wary of his results and perhaps
c     increase u for subsequent runs.  we would like to emphasise that
c     this value only relates to the accuracy of our lu decomposition
c     and gives no indication as to the singularity of the matrix or the
c     accuracy of the solution.  this upper bound can be a significant
c     overestimate particularly if the matrix is badly scaled. if an
c     accurate value for the growth is required, lbig (q.v.) should be
c     set to .true.
c eps,rmin are real variables. if, on entry to ma28b/bd, eps is less
c     than one, then rmin will give the smallest ratio of the pivot to
c     the largest element in the corresponding row of the upper
c     triangular factor thus monitoring the stability of successive
c     factorizations. if rmin becomes very large and w(1) from
c     ma28b/bd is also very large, it may be advisable to perform a
c     new decomposition using ma28a/ad.
c resid is a real variable which on exit from ma28c/cd gives the value
c     of the maximum residual over all the equations unsatisfied because
c     of dependency (zero pivots).
c irncp,icncp are integer variables which monitor the adequacy of "elbow
c     room" in irn and a/icn respectively. if either is quite large (say
c     greater than n/10), it will probably pay to increase the size of
c     the corresponding array for subsequent runs. if either is very low
c     or zero then one can perhaps save storage by reducing the size of
c     the corresponding array.
c minirn,minicn are integer variables which, in the event of a
c     successful return (iflag ge 0 or iflag=-14) give the minimum size
c     of irn and a/icn respectively which would enable a successful run
c     on an identical matrix. on an exit with iflag equal to -5, minicn
c     gives the minimum value of icn for success on subsequent runs on
c     an identical matrix. in the event of failure with iflag= -6, -4,
c     -3, -2, or -1, then minicn and minirn give the minimum value of
c     licn and lirn respectively which would be required for a
c     successful decomposition up to the point at which the failure
c     occurred.
c irank is an integer variable which gives an upper bound on the rank of
c     the matrix.
c abort1 is a logical variable with default value .true.  if abort1 is
c     set to .false.  then ma28a/ad will decompose structurally singular
c     matrices (including rectangular ones).
c abort2 is a logical variable with default value .true.  if abort2 is
c     set to .false. then ma28a/ad will decompose numerically singular
c     matrices.
c idisp is an integer array of length 2. on output from ma28a/ad, the
c     indices of the diagonal blocks of the factors lie in positions
c     idisp(1) to idisp(2) of a/icn. this array must be preserved
c     between a call to ma28a/ad and subsequent calls to ma28b/bd,
c     ma28c/cd or ma28i/id.
c tol is a real variable.  if it is set to a positive value, then any
c     non-zero whose modulus is less than tol will be dropped from the
c     factorization.  the factorization will then require less storage
c     but will be inaccurate.  after a run of ma28a/ad with tol positive
c     it is not possible to use ma28b/bd and the user is recommended to
c     use ma28i/id to obtain the solution.  the default value for tol is
c     0.0.
c themax is a real variable.  on exit from ma28a/ad, it will hold the
c     largest entry of the original matrix.
c big is a real variable. if lbig has been set to .true., big will hold
c     the largest entry encountered during the factorization by ma28a/ad
c     or ma28b/bd.
c dxmax is a real variable. on exit from ma28i/id, dxmax will be set to
c     the largest component of the solution.
c errmax is a real variable.  on exit from ma28i/id, if maxit is
c     positive, errmax will be set to the largest component in the
c     estimate of the error.
c dres is a real variable.  on exit from ma28i/id, if maxit is positive,
c     dres will be set to the largest component of the residual.
c cgce is a real variable. it is used by ma28i/id to check the
c     convergence rate.  if the ratio of successive corrections is
c     not less than cgce then we terminate since the convergence
c     rate is adjudged too slow.
c ndrop is an integer variable. if tol has been set positive, on exit
c     from ma28a/ad, ndrop will hold the number of entries dropped from
c     the data structure.
c maxit is an integer variable. it is the maximum number of iterations
c     performed by ma28i/id. it has a default value of 16.
c noiter is an integer variable. it is set by ma28i/id to the number of
c     iterative refinement iterations actually used.
c nsrch is an integer variable. if nsrch is set to a value less than n,
c     then a different pivot option will be employed by ma28a/ad.  this
c     may result in different fill-in and execution time for ma28a/ad.
c     if nsrch is less than or equal to n, the workspace array iw can be
c     reduced in length.  the default value for nsrch is 32768.
c istart is an integer variable. if istart is set to a value other than
c     zero, then the user must supply an estimate of the solution to
c     ma28i/id.  the default value for istart is zero.
c lbig is a logical variable. if lbig is set to .true., the value of the
c     largest element encountered in the factorization by ma28a/ad or
c     ma28b/bd is returned in big.  setting lbig to .true.  will
c     increase the time for ma28a/ad marginally and that for ma28b/bd
c     by about 20%.  the default value for lbig is .false.
c
      DOUBLE PRECISION EPS, RMIN, RESID, TOL, THEMAX, BIG, DXMAX,
     * ERRMAX, DRES, CGCE
      LOGICAL LBLOCK, GROW, ABORT1, ABORT2, LBIG
      COMMON /MA28ED/ LP, MP, LBLOCK, GROW
      COMMON /MA28FD/ EPS, RMIN, RESID, IRNCP, ICNCP, MINIRN, MINICN,
     * IRANK, ABORT1, ABORT2
c     common /ma28gd/ idisp(2)
      COMMON /MA28HD/ TOL, THEMAX, BIG, DXMAX, ERRMAX, DRES, CGCE,
     * NDROP, MAXIT, NOITER, NSRCH, ISTART, LBIG
      DATA EPS /1.0D-4/, TOL /0.0D0/, CGCE /0.5D0/
      DATA MAXIT /16/
      DATA LP /6/, MP /6/, NSRCH /32768/, ISTART /0/
      DATA LBLOCK /.TRUE./, GROW /.TRUE./, LBIG /.FALSE./
      DATA ABORT1 /.TRUE./, ABORT2 /.TRUE./
      END
      BLOCK DATA MABLD2
c although all common block variables do not have default values,
c     we comment on all the common block variables here.
c
c common block ma30e/ed holds control parameters ....
c     common /ma30ed/ lp, abort1, abort2, abort3
c the integer lp is the unit number to which the error messages are
c     sent. lp has a default value of 6.  this default value can be
c     reset by the user, if desired.  a value of 0 suppresses all
c     messages.
c the logical variables abort1,abort2,abort3 are used to control the
c     conditions under which the subroutine will terminate.
c if abort1 is .true. then the subroutine will exit  immediately on
c     detecting structural singularity.
c if abort2 is .true. then the subroutine will exit immediately on
c     detecting numerical singularity.
c if abort3 is .true. then the subroutine will exit immediately when
c     the available space in a/icn is filled up by the previously
c     decomposed, active, and undecomposed parts of the matrix.
c the default values for abort1,abort2,abort3 are set to .true.,.true.
c     and .false. respectively.
c
c the variables in the common block ma30f/fd are used to provide the
c     user with information on the decomposition.
c     common /ma30fd/ irncp, icncp, irank, minirn, minicn
c irncp and icncp are integer variables used to monitor the adequacy
c     of the allocated space in arrays irn and a/icn respectively, by
c     taking account of the number of data management compresses
c     required on these arrays. if irncp or icncp is fairly large (say
c     greater than n/10), it may be advantageous to increase the size
c     of the corresponding array(s).  irncp and icncp are initialized
c     to zero on entry to ma30a/ad and are incremented each time the
c     compressing routine ma30d/dd is entered.
c icncp is the number of compresses on a/icn.
c irncp is the number of compresses on irn.
c irank is an integer variable which gives an estimate (actually an
c     upper bound) of the rank of the matrix. on an exit with iflag
c     equal to 0, this will be equal to n.
c minirn is an integer variable which, after a successful call to
c     ma30a/ad, indicates the minimum length to which irn can be
c     reduced while still permitting a successful decomposition of the
c     same matrix. if, however, the user were to decrease the length
c     of irn to that size, the number of compresses (irncp) may be
c     very high and quite costly. if lirn is not large enough to begin
c     the decomposition on a diagonal block, minirn will be equal to
c     the value required to continue the decomposition and iflag will
c     be set to -3 or -6. a value of lirn slightly greater than this
c     (say about n/2) will usually provide enough space to complete
c     the decomposition on that block. in the event of any other
c     failure minirn gives the minimum size of irn required for a
c     successful decomposition up to that point.
c minicn is an integer variable which after a successful call to
c     ma30a/ad, indicates the minimum size of licn required to enable
c     a successful decomposition. in the event of failure with iflag=
c     -5, minicn will, if abort3 is left set to .false., indicate the
c     minimum length that would be sufficient to prevent this error in
c     a subsequent run on an identical matrix. again the user may
c     prefer to use a value of icn slightly greater than minicn for
c     subsequent runs to avoid too many conpresses (icncp). in the
c     event of failure with iflag equal to any negative value except
c     -4, minicn will give the minimum length to which licn could be
c     reduced to enable a successful decomposition to the point at
c     which failure occurred.  notice that, on a successful entry
c     idisp(2) gives the amount of space in a/icn required for the
c     decomposition while minicn will usually be slightly greater
c     because of the need for "elbow room".  if the user is very
c     unsure how large to make licn, the variable minicn can be used
c     to provide that information. a preliminary run should be
c     performed with abort3 left set to .false. and licn about 3/2
c     times as big as the number of non-zeros in the original matrix.
c     unless the initial problem is very sparse (when the run will be
c     successful) or fills in extremely badly (giving an error return
c     with iflag equal to -4), an error return with iflag equal to -5
c     should result and minicn will give the amount of space required
c     for a successful decomposition.
c
c common block ma30g/gd is used by the ma30b/bd entry only.
c     common /ma30gd/ eps, rmin
c eps is a real/double precision variable. it is used to test for
c     small pivots. its default value is 1.0e-4 (1.0d-4 in d version).
c     if the user sets eps to any value greater than 1.0, then no
c     check is made on the size of the pivots. although the absence of
c     such a check would fail to warn the user of bad instability, its
c     absence will enable ma30b/bd to run slightly faster. an  a
c     posteriori  check on the stability of the factorization can be
c     obtained from mc24a/ad.
c rmin is a real/double precision variable which gives the user some
c     information about the stability of the decomposition.  at each
c     stage of the lu decomposition the magnitude of the pivot apiv
c     is compared with the largest off-diagonal entry currently in its
c     row (row of u), rowmax say. if the ratio
c                       min (apiv/rowmax)
c     where the minimum is taken over all the rows, is less than eps
c     then rmin is set to this minimum value and iflag is returned
c     with the value +i where i is the row in which this minimum
c     occurs.  if the user sets eps greater than one, then this test
c     is not performed. in this case, and when there are no small
c     pivots rmin will be set equal to eps.
c
c common block ma30h/hd is used by ma30c/cd only.
c     common /ma30hd/ resid
c resid is a real/double precision variable. in the case of singular
c     or rectangular matrices its final value will be equal to the
c     maximum residual for the unsatisfied equations; otherwise its
c     value will be set to zero.
c
c common  block ma30i/id controls the use of drop tolerances, the
c     modified pivot option and the the calculation of the largest
c     entry in the factorization process. this common block was added
c     to the ma30 package in february, 1983.
c     common /ma30id/ tol, big, ndrop, nsrch, lbig
c tol is a real/double precision variable.  if it is set to a positive
c     value, then ma30a/ad will drop from the factors any non-zero
c     whose modulus is less than tol.  the factorization will then
c     require less storage but will be inaccurate.  after a run of
c     ma30a/ad where entries have been dropped, ma30b/bd  should not
c     be called.  the default value for tol is 0.0.
c big is a real/double precision variable.  if lbig has been set to
c     .true., big will be set to the largest entry encountered during
c     the factorization.
c ndrop is an integer variable. if tol has been set positive, on exit
c     from ma30a/ad, ndrop will hold the number of entries dropped
c     from the data structure.
c nsrch is an integer variable. if nsrch is set to a value less than
c     or equal to n, then a different pivot option will be employed by
c     ma30a/ad.  this may result in different fill-in and execution
c     time for ma30a/ad. if nsrch is less than or equal to n, the
c     workspace arrays lastc and nextc are not referenced by ma30a/ad.
c     the default value for nsrch is 32768.
c lbig is a logical variable. if lbig is set to .true., the value of
c     the largest entry encountered in the factorization by ma30a/ad
c     is returned in big.  setting lbig to .true.  will marginally
c     increase the factorization time for ma30a/ad and will increase
c     that for ma30b/bd by about 20%.  the default value for lbig is
c     .false.
c
      DOUBLE PRECISION EPS, RMIN, TOL, BIG
      LOGICAL ABORT1, ABORT2, ABORT3, LBIG
      COMMON /MA30ED/ LP, ABORT1, ABORT2, ABORT3
      COMMON /MA30GD/ EPS, RMIN
      COMMON /MA30ID/ TOL, BIG, NDROP, NSRCH, LBIG
      DATA EPS /1.0D-4/, TOL /0.0D0/, BIG /0.0D0/
      DATA LP /6/, NSRCH /32768/
      DATA LBIG /.FALSE./
      DATA ABORT1 /.TRUE./, ABORT2 /.TRUE./, ABORT3 /.FALSE./
      END
      BLOCK DATA MABLD3
      LOGICAL ABORT
      COMMON /MC23BD/ LP,NUMNZ,NUM,LARGE,ABORT
      DATA LP/6/,ABORT/.FALSE./
      END
