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