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