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