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