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