C The codes MA28**, MA30**, MC*** are part of the Harwell Subroutine Library,
C which is developed
C and maintained by staff at the Harwell Laboratory of the United Kingdom
C Atomic Energy Authority. The researcher may use this code SOLELY for
C research purposes subject to the conditions below.
C
C
C On receipt of this code the requester must return the following
C form (signed) to Mr. Sid Marlow, Computer Science and Systems Division,
C B7.12, AERE Harwell, Didcot, Oxon, England OX11 0RA.
C
C
C The conditions as specified below are accepted and shall be complied with.
C
C
C Signature ........................... Date ............................
C
C Name (BLOCK CAPITALS please) ............................................
C
C on behalf of (BLOCK CAPITALS please) ....................................
C
C Address..................................................................
C
C ........................................................................
C
C Position ................................................................
C
C
C The conditions of use of the Harwell Subroutine Library or any part
C thereof by any external user are as follows:
C
C (i) The subroutines may only be used for research purposes by the person or
C organization to whom they are supplied ("the User"), and must not be
C copied by the User for use by any other persons or organizations.
C All information on the subroutines is provided by the Atomic Energy
C Authority ("the Authority") to the User on the understanding that
C the details thereof are confidential.
C
C (ii) All publications issued by the User which include results obtained
C with the help of one or more of the subroutines supplied shall
C acknowledge the use of the subroutines.
C
C (iii) The subroutines supplied may be modified by or on behalf of the User
C for such use in research applications but at no time shall such
C subroutines or modifications thereof become the property of the User.
C The Authority shall not be liable for any direct or consequential
C loss or damage whatsoever arising out of the use of subroutines by
C the User.
C
C (iv) Any use of the subroutines supplied in any commercial application
C shall be subject to prior written agreement between the Authority
C and the User on suitable terms and conditions (including financial
C conditions).
C
c ---------------------------------------------------------
C
C########################################################################
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
c
c==============================================================
c
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
c
c============================================================
c
SUBROUTINE MA28CD(N, A, LICN, ICN, IKEEP, RHS, W, MTYPE)
c
c this subroutine uses the factors from ma28a/ad or ma28b/bd to
c solve a system of equations without iterative refinement.
c the parameters are ...
c n integer order of matrix not altered by subroutine.
c a real/double precision array length licn. the same array as
c was used in the most recent call to ma28a/ad or ma28b/bd.
c licn integer length of arrays a and icn. not altered by
c subroutine.
c icn integer array of length licn. same array as output from
c ma28a/ad. unchanged by ma28c/cd.
c ikeep integer array of length 5*n. same array as output from
c ma28a/ad. unchanged by ma28c/cd.
c rhs real/double precision array length n. on entry, it holds the
c right hand side. on exit, the solution vector.
c w real/double precision array length n. used as workspace by
c ma30c/cd.
c mtype integer used to tell ma30c/cd to solve the direct equation
c (mtype=1) or its transpose (mtype.ne.1).
c
DOUBLE PRECISION A(LICN), RHS(N), W(N), RESID, MRESID, EPS, RMIN
INTEGER IDISP(2)
INTEGER ICN(LICN), IKEEP(N,5)
LOGICAL ABORT1, ABORT2
c common block variables.
c unless otherwise stated common block variables are as in ma28a/ad.
c those variables referenced by ma28c/cd are mentioned below.
c resid real/double precision variable returns maximum residual of
c equations where pivot was zero.
c mresid real/double precision variable used by ma28c/cd to
c communicate between ma28f/fd and ma30h/hd.
c idisp integer array length 2 the same as that used by ma28a/ad.
c it is unchanged by ma28b/bd.
c
c further information on common block variables can be found in block
c data or ma28a/ad.
COMMON /MA28FD/ EPS, RMIN, RESID, IRNCP, ICNCP, MINIRN, MINICN,
* IRANK, ABORT1, ABORT2
COMMON /MA28GD/ IDISP
COMMON /MA30HD/ MRESID
c
c this call performs the solution of the set of equations.
CALL MA30CD(N, ICN, A, LICN, IKEEP, IKEEP(1,4), IKEEP(1,5),
* IDISP, IKEEP(1,2), IKEEP(1,3), RHS, W, MTYPE)
c transfer common block information.
RESID = MRESID
RETURN
END
SUBROUTINE MA28DD(N, A, LICN, IVECT, JVECT, NZ, ICN, LENR, LENRL,
* LENOFF, IP, IQ, IW1, IW, W1, IFLAG)
c this subroutine need never be called by the user directly.
c it sorts the user's matrix into the structure of the decomposed
c form and checks for the presence of duplicate entries or
c non-zeros lying outside the sparsity pattern of the decomposition
c it also calculates the largest element in the input matrix.
DOUBLE PRECISION A(LICN), ZERO, W1, AA
INTEGER IW(N,2), IDISP(2)
INTEGER ICN(LICN), IVECT(NZ), JVECT(NZ), IP(N), IQ(N),
* LENR(N), IW1(N,3), LENRL(N), LENOFF(N)
LOGICAL LBLOCK, GROW, BLOCKL
COMMON /MA28ED/ LP, MP, LBLOCK, GROW
COMMON /MA28GD/ IDISP
DATA ZERO /0.0D0/
BLOCKL = LENOFF(1).GE.0
c iw1(i,3) is set to the block in which row i lies and the
c inverse permutations to ip and iq are set in iw1(.,1) and
c iw1(.,2) resp.
c pointers to beginning of the part of row i in diagonal and
c off-diagonal blocks are set in iw(i,2) and iw(i,1) resp.
IBLOCK = 1
IW(1,1) = 1
IW(1,2) = IDISP(1)
DO 10 I=1,N
IW1(I,3) = IBLOCK
IF (IP(I).LT.0) IBLOCK = IBLOCK + 1
II = IABS(IP(I)+0)
IW1(II,1) = I
JJ = IQ(I)
JJ = IABS(JJ)
IW1(JJ,2) = I
IF (I.EQ.1) GO TO 10
IF (BLOCKL) IW(I,1) = IW(I-1,1) + LENOFF(I-1)
IW(I,2) = IW(I-1,2) + LENR(I-1)
10 CONTINUE
c place each non-zero in turn into its correct location
c in the a/icn array.
IDISP2 = IDISP(2)
DO 170 I=1,NZ
c necessary to avoid reference to unassigned element of icn.
IF (I.GT.IDISP2) GO TO 20
IF (ICN(I).LT.0) GO TO 170
20 IOLD = IVECT(I)
JOLD = JVECT(I)
AA = A(I)
c this is a dummy loop for following a chain of interchanges.
c it will be executed nz times in total.
DO 140 IDUMMY=1,NZ
c perform some validity checks on iold and jold.
IF (IOLD.LE.N .AND. IOLD.GT.0 .AND. JOLD.LE.N .AND.
* JOLD.GT.0) GO TO 30
IF (LP.NE.0) WRITE (LP,99999) I, A(I), IOLD, JOLD
IFLAG = -12
GO TO 180
30 INEW = IW1(IOLD,1)
JNEW = IW1(JOLD,2)
c are we in a valid block and is it diagonal or off-diagonal?
IF (IW1(INEW,3)-IW1(JNEW,3)) 40, 60, 50
40 IFLAG = -13
IF (LP.NE.0) WRITE (LP,99998) IOLD, JOLD
GO TO 180
50 J1 = IW(INEW,1)
J2 = J1 + LENOFF(INEW) - 1
GO TO 110
c element is in diagonal block.
60 J1 = IW(INEW,2)
IF (INEW.GT.JNEW) GO TO 70
J2 = J1 + LENR(INEW) - 1
J1 = J1 + LENRL(INEW)
GO TO 110
70 J2 = J1 + LENRL(INEW)
c binary search of ordered list .. element in l part of row.
DO 100 JDUMMY=1,N
MIDPT = (J1+J2)/2
JCOMP = IABS(ICN(MIDPT)+0)
IF (JNEW-JCOMP) 80, 130, 90
80 J2 = MIDPT
GO TO 100
90 J1 = MIDPT
100 CONTINUE
IFLAG = -13
IF (LP.NE.0) WRITE (LP,99997) IOLD, JOLD
GO TO 180
c linear search ... element in l part of row or off-diagonal blocks.
110 DO 120 MIDPT=J1,J2
IF (IABS(ICN(MIDPT)+0).EQ.JNEW) GO TO 130
120 CONTINUE
IFLAG = -13
IF (LP.NE.0) WRITE (LP,99997) IOLD, JOLD
GO TO 180
c equivalent element of icn is in position midpt.
130 IF (ICN(MIDPT).LT.0) GO TO 160
IF (MIDPT.GT.NZ .OR. MIDPT.LE.I) GO TO 150
W1 = A(MIDPT)
A(MIDPT) = AA
AA = W1
IOLD = IVECT(MIDPT)
JOLD = JVECT(MIDPT)
ICN(MIDPT) = -ICN(MIDPT)
140 CONTINUE
150 A(MIDPT) = AA
ICN(MIDPT) = -ICN(MIDPT)
GO TO 170
160 A(MIDPT) = A(MIDPT) + AA
c set flag for duplicate elements.
IFLAG = N + 1
170 CONTINUE
c reset icn array and zero elements in l/u but not in a.
c also calculate maximum element of a.
180 W1 = ZERO
DO 200 I=1,IDISP2
IF (ICN(I).LT.0) GO TO 190
A(I) = ZERO
GO TO 200
190 ICN(I) = -ICN(I)
W1 = DMAX1(W1,DABS(A(I)))
200 CONTINUE
RETURN
99999 FORMAT (9H ELEMENT , I6, 12H WITH VALUE , 1PD22.14, 10H HAS INDIC,
* 3HES , I8, 2H ,, I8/36X, 20HINDICES OUT OF RANGE)
99998 FORMAT (36X, 8HNON-ZERO, I7, 2H ,, I6, 23H IN ZERO OFF-DIAGONAL B,
* 4HLOCK)
99997 FORMAT (36X, 8H ELEMENT, I6, 2H ,, I6, 23H WAS NOT IN L/U PATTERN)
END
c######date 01 jan 1984 copyright ukaea, harwell.
c######alias ma30ad
c
c===================================================================
c
SUBROUTINE MA30AD(NN, ICN, A, LICN, LENR, LENRL, IDISP, IP, IQ,
* IRN, LIRN, LENC, IFIRST, LASTR, NEXTR, LASTC, NEXTC, IPTR, IPC,
* U, IFLAG)
c if the user requires a more convenient data interface then the ma28
c package should be used. the ma28 subroutines call the ma30
c subroutines after checking the user's input data and optionally
c using mc23a/ad to permute the matrix to block triangular form.
c this package of subroutines (ma30a/ad, ma30b/bd, ma30c/cd and
c ma30d/dd) performs operations pertinent to the solution of a
c general sparse n by n system of linear equations (i.e. solve
c ax=b). structually singular matrices are permitted including
c those with row or columns consisting entirely of zeros (i.e.
c including rectangular matrices). it is assumed that the
c non-zeros of the matrix a do not differ widely in size. if
c necessary a prior call of the scaling subroutine mc19a/ad may be
c made.
c a discussion of the design of these subroutines is given by duff and
c reid (acm trans math software 5 pp 18-35,1979 (css 48)) while
c fuller details of the implementation are given in duff (harwell
c report aere-r 8730,1977). the additional pivoting option in
c ma30a/ad and the use of drop tolerances (see common block
c ma30i/id) were added to the package after joint work with reid,
c schaumburg, wasniewski and zlatev (duff, reid, schaumburg,
c wasniewski and zlatev, harwell report css 135, 1983).
c
c ma30a/ad performs the lu decomposition of the diagonal blocks of the
c permutation paq of a sparse matrix a, where input permutations
c p1 and q1 are used to define the diagonal blocks. there may be
c non-zeros in the off-diagonal blocks but they are unaffected by
c ma30a/ad. p and p1 differ only within blocks as do q and q1. the
c permutations p1 and q1 may be found by calling mc23a/ad or the
c matrix may be treated as a single block by using p1=q1=i. the
c matrix non-zeros should be held compactly by rows, although it
c should be noted that the user can supply the matrix by columns
c to get the lu decomposition of a transpose.
c
c the parameters are...
c this description should also be consulted for further information on
c most of the parameters of ma30b/bd and ma30c/cd.
c
c n is an integer variable which must be set by the user to the order
c of the matrix. it is not altered by ma30a/ad.
c icn is an integer array of length licn. positions idisp(2) to
c licn must be set by the user to contain the column indices of
c the non-zeros in the diagonal blocks of p1*a*q1. those belonging
c to a single row must be contiguous but the ordering of column
c indices with each row is unimportant. the non-zeros of row i
c precede those of row i+1,i=1,...,n-1 and no wasted space is
c allowed between the rows. on output the column indices of the
c lu decomposition of paq are held in positions idisp(1) to
c idisp(2), the rows are in pivotal order, and the column indices
c of the l part of each row are in pivotal order and precede those
c of u. again there is no wasted space either within a row or
c between the rows. icn(1) to icn(idisp(1)-1), are neither
c required nor altered. if mc23a/ad been called, these will hold
c information about the off-diagonal blocks.
c a is a real/double precision array of length licn whose entries
c idisp(2) to licn must be set by the user to the values of the
c non-zero entries of the matrix in the order indicated by icn.
c on output a will hold the lu factors of the matrix where again
c the position in the matrix is determined by the corresponding
c values in icn. a(1) to a(idisp(1)-1) are neither required nor
c altered.
c licn is an integer variable which must be set by the user to the
c length of arrays icn and a. it must be big enough for a and icn
c to hold all the non-zeros of l and u and leave some "elbow
c room". it is possible to calculate a minimum value for licn by
c a preliminary run of ma30a/ad. the adequacy of the elbow room
c can be judged by the size of the common block variable icncp. it
c is not altered by ma30a/ad.
c lenr is an integer array of length n. on input, lenr(i) should
c equal the number of non-zeros in row i, i=1,...,n of the
c diagonal blocks of p1*a*q1. on output, lenr(i) will equal the
c total number of non-zeros in row i of l and row i of u.
c lenrl is an integer array of length n. on output from ma30a/ad,
c lenrl(i) will hold the number of non-zeros in row i of l.
c idisp is an integer array of length 2. the user should set idisp(1)
c to be the first available position in a/icn for the lu
c decomposition while idisp(2) is set to the position in a/icn of
c the first non-zero in the diagonal blocks of p1*a*q1. on output,
c idisp(1) will be unaltered while idisp(2) will be set to the
c position in a/icn of the last non-zero of the lu decomposition.
c ip is an integer array of length n which holds a permutation of
c the integers 1 to n. on input to ma30a/ad, the absolute value of
c ip(i) must be set to the row of a which is row i of p1*a*q1. a
c negative value for ip(i) indicates that row i is at the end of a
c diagonal block. on output from ma30a/ad, ip(i) indicates the row
c of a which is the i th row in paq. ip(i) will still be negative
c for the last row of each block (except the last).
c iq is an integer array of length n which again holds a
c permutation of the integers 1 to n. on input to ma30a/ad, iq(j)
c must be set to the column of a which is column j of p1*a*q1. on
c output from ma30a/ad, the absolute value of iq(j) indicates the
c column of a which is the j th in paq. for rows, i say, in which
c structural or numerical singularity is detected iq(i) is
c negated.
c irn is an integer array of length lirn used as workspace by
c ma30a/ad.
c lirn is an integer variable. it should be greater than the
c largest number of non-zeros in a diagonal block of p1*a*q1 but
c need not be as large as licn. it is the length of array irn and
c should be large enough to hold the active part of any block,
c plus some "elbow room", the a posteriori adequacy of which can
c be estimated by examining the size of common block variable
c irncp.
c lenc,ifirst,lastr,nextr,lastc,nextc are all integer arrays of
c length n which are used as workspace by ma30a/ad. if nsrch is
c set to a value less than or equal to n, then arrays lastc and
c nextc are not referenced by ma30a/ad and so can be dummied in
c the call to ma30a/ad.
c iptr,ipc are integer arrays of length n which are used as workspace
c by ma30a/ad.
c u is a real/double precision variable which should be set by the
c user to a value between 0. and 1.0. if less than zero it is
c reset to zero and if its value is 1.0 or greater it is reset to
c 0.9999 (0.999999999 in d version). it determines the balance
c between pivoting for sparsity and for stability, values near
c zero emphasizing sparsity and values near one emphasizing
c stability. we recommend u=0.1 as a posible first trial value.
c the stability can be judged by a later call to mc24a/ad or by
c setting lbig to .true.
c iflag is an integer variable. it will have a non-negative value if
c ma30a/ad is successful. negative values indicate error
c conditions while positive values indicate that the matrix has
c been successfully decomposed but is singular. for each non-zero
c value, an appropriate message is output on unit lp. possible
c non-zero values for iflag are ...
c
c -1 the matrix is structually singular with rank given by irank in
c common block ma30f/fd.
c +1 if, however, the user wants the lu decomposition of a
c structurally singular matrix and sets the common block variable
c abort1 to .false., then, in the event of singularity and a
c successful decomposition, iflag is returned with the value +1
c and no message is output.
c -2 the matrix is numerically singular (it may also be structually
c singular) with estimated rank given by irank in common block
c ma30f/fd.
c +2 the user can choose to continue the decomposition even when a
c zero pivot is encountered by setting common block variable
c abort2 to .false. if a singularity is encountered, iflag will
c then return with a value of +2, and no message is output if the
c decomposition has been completed successfully.
c -3 lirn has not been large enough to continue with the
c decomposition. if the stage was zero then common block variable
c minirn gives the length sufficient to start the decomposition on
c this block. for a successful decomposition on this block the
c user should make lirn slightly (say about n/2) greater than this
c value.
c -4 licn not large enough to continue with the decomposition.
c -5 the decomposition has been completed but some of the lu factors
c have been discarded to create enough room in a/icn to continue
c the decomposition. the variable minicn in common block ma30f/fd
c then gives the size that licn should be to enable the
c factorization to be successful. if the user sets common block
c variable abort3 to .true., then the subroutine will exit
c immediately instead of destroying any factors and continuing.
c -6 both licn and lirn are too small. termination has been caused by
c lack of space in irn (see error iflag= -3), but already some of
c the lu factors in a/icn have been lost (see error iflag= -5).
c minicn gives the minimum amount of space required in a/icn for
c decomposition up to this point.
c
DOUBLE PRECISION A(LICN), U, AU, UMAX, AMAX, ZERO, PIVRAT, PIVR,
* TOL, BIG, ANEW, AANEW, SCALE
INTEGER IPTR(NN), PIVOT, PIVEND, DISPC, OLDPIV, OLDEND, PIVROW,
* ROWI, IPC(NN), IDISP(2), COLUPD
INTEGER ICN(LICN), LENR(NN), LENRL(NN), IP(NN), IQ(NN),
* LENC(NN), IRN(LIRN), IFIRST(NN), LASTR(NN), NEXTR(NN),
* LASTC(NN), NEXTC(NN)
LOGICAL ABORT1, ABORT2, ABORT3, LBIG
c for comments of common block variables see block data subprogram.
COMMON /MA30ED/ LP, ABORT1, ABORT2, ABORT3
COMMON /MA30FD/ IRNCP, ICNCP, IRANK, MINIRN, MINICN
COMMON /MA30ID/ TOL, BIG, NDROP, NSRCH, LBIG
COMMON /LPIVOT/ LPIV(10),LNPIV(10),MAPIV,MANPIV,IAVPIV,
* IANPIV,KOUNTL
c
DATA UMAX/.999999999D0/
DATA ZERO /0.0D0/
MSRCH = NSRCH
NDROP = 0
DO 1272 KK=1,10
LNPIV(KK)=0
LPIV(KK)=0
1272 CONTINUE
MAPIV = 0
MANPIV = 0
IAVPIV = 0
IANPIV = 0
KOUNTL = 0
MINIRN = 0
MINICN = IDISP(1) - 1
MOREI = 0
IRANK = NN
IRNCP = 0
ICNCP = 0
IFLAG = 0
c reset u if necessary.
U = DMIN1(U,UMAX)
c ibeg is the position of the next pivot row after elimination step
c using it.
U = DMAX1(U,ZERO)
IBEG = IDISP(1)
c iactiv is the position of the first entry in the active part of a/icn.
IACTIV = IDISP(2)
c nzrow is current number of non-zeros in active and unprocessed part
c of row file icn.
NZROW = LICN - IACTIV + 1
MINICN = NZROW + MINICN
c
c count the number of diagonal blocks and set up pointers to the
c beginnings of the rows.
c num is the number of diagonal blocks.
NUM = 1
IPTR(1) = IACTIV
IF (NN.EQ.1) GO TO 20
NNM1 = NN - 1
DO 10 I=1,NNM1
IF (IP(I).LT.0) NUM = NUM + 1
IPTR(I+1) = IPTR(I) + LENR(I)
10 CONTINUE
c ilast is the last row in the previous block.
20 ILAST = 0
c
c ***********************************************
c **** lu decomposition of block nblock ****
c ***********************************************
c
c each pass through this loop performs lu decomposition on one
c of the diagonal blocks.
DO 1000 NBLOCK=1,NUM
ISTART = ILAST + 1
DO 30 IROWS=ISTART,NN
IF (IP(IROWS).LT.0) GO TO 40
30 CONTINUE
IROWS = NN
40 ILAST = IROWS
c n is the number of rows in the current block.
c istart is the index of the first row in the current block.
c ilast is the index of the last row in the current block.
c iactiv is the position of the first entry in the block.
c itop is the position of the last entry in the block.
N = ILAST - ISTART + 1
IF (N.NE.1) GO TO 90
c
c code for dealing with 1x1 block.
LENRL(ILAST) = 0
ISING = ISTART
IF (LENR(ILAST).NE.0) GO TO 50
c block is structurally singular.
IRANK = IRANK - 1
ISING = -ISING
IF (IFLAG.NE.2 .AND. IFLAG.NE.-5) IFLAG = 1
IF (.NOT.ABORT1) GO TO 80
IDISP(2) = IACTIV
IFLAG = -1
IF (LP.NE.0) WRITE (LP,99999)
c return
GO TO 1120
50 SCALE = DABS(A(IACTIV))
IF (SCALE.EQ.ZERO) GO TO 60
IF (LBIG) BIG = DMAX1(BIG,SCALE)
GO TO 70
60 ISING = -ISING
IRANK = IRANK - 1
IPTR(ILAST) = 0
IF (IFLAG.NE.-5) IFLAG = 2
IF (.NOT.ABORT2) GO TO 70
IDISP(2) = IACTIV
IFLAG = -2
IF (LP.NE.0) WRITE (LP,99998)
GO TO 1120
70 A(IBEG) = A(IACTIV)
ICN(IBEG) = ICN(IACTIV)
IACTIV = IACTIV + 1
IPTR(ISTART) = 0
IBEG = IBEG + 1
NZROW = NZROW - 1
80 LASTR(ISTART) = ISTART
IPC(ISTART) = -ISING
GO TO 1000
c
c non-trivial block.
90 ITOP = LICN
IF (ILAST.NE.NN) ITOP = IPTR(ILAST+1) - 1
c
c set up column oriented storage.
DO 100 I=ISTART,ILAST
LENRL(I) = 0
LENC(I) = 0
100 CONTINUE
IF (ITOP-IACTIV.LT.LIRN) GO TO 110
MINIRN = ITOP - IACTIV + 1
PIVOT = ISTART - 1
GO TO 1100
c
c calculate column counts.
110 DO 120 II=IACTIV,ITOP
I = ICN(II)
LENC(I) = LENC(I) + 1
120 CONTINUE
c set up column pointers so that ipc(j) points to position after end
c of column j in column file.
IPC(ILAST) = LIRN + 1
J1 = ISTART + 1
DO 130 JJ=J1,ILAST
J = ILAST - JJ + J1 - 1
IPC(J) = IPC(J+1) - LENC(J+1)
130 CONTINUE
DO 150 INDROW=ISTART,ILAST
J1 = IPTR(INDROW)
J2 = J1 + LENR(INDROW) - 1
IF (J1.GT.J2) GO TO 150
DO 140 JJ=J1,J2
J = ICN(JJ)
IPOS = IPC(J) - 1
IRN(IPOS) = INDROW
IPC(J) = IPOS
140 CONTINUE
150 CONTINUE
c dispc is the lowest indexed active location in the column file.
DISPC = IPC(ISTART)
NZCOL = LIRN - DISPC + 1
MINIRN = MAX0(NZCOL,MINIRN)
NZMIN = 1
c
c initialize array ifirst. ifirst(i) = +/- k indicates that row/col
c k has i non-zeros. if ifirst(i) = 0, there is no row or column
c with i non zeros.
DO 160 I=1,N
IFIRST(I) = 0
160 CONTINUE
c
c compute ordering of row and column counts.
c first run through columns (from column n to column 1).
DO 180 JJ=ISTART,ILAST
J = ILAST - JJ + ISTART
NZ = LENC(J)
IF (NZ.NE.0) GO TO 170
IPC(J) = 0
GO TO 180
170 IF (NSRCH.LE.NN) GO TO 180
ISW = IFIRST(NZ)
IFIRST(NZ) = -J
LASTC(J) = 0
NEXTC(J) = -ISW
ISW1 = IABS(ISW)
IF (ISW.NE.0) LASTC(ISW1) = J
180 CONTINUE
c now run through rows (again from n to 1).
DO 210 II=ISTART,ILAST
I = ILAST - II + ISTART
NZ = LENR(I)
IF (NZ.NE.0) GO TO 190
IPTR(I) = 0
LASTR(I) = 0
GO TO 210
190 ISW = IFIRST(NZ)
IFIRST(NZ) = I
IF (ISW.GT.0) GO TO 200
NEXTR(I) = 0
LASTR(I) = ISW
GO TO 210
200 NEXTR(I) = ISW
LASTR(I) = LASTR(ISW)
LASTR(ISW) = I
210 CONTINUE
c
c
c **********************************************
c **** start of main elimination loop ****
c **********************************************
DO 980 PIVOT=ISTART,ILAST
c
c first find the pivot using markowitz criterion with stability
c control.
c jcost is the markowitz cost of the best pivot so far,.. this
c pivot is in row ipiv and column jpiv.
NZ2 = NZMIN
JCOST = N*N
c
c examine rows/columns in order of ascending count.
DO 340 L=1,2
PIVRAT = ZERO
ISRCH = 1
LL = L
c a pass with l equal to 2 is only performed in the case of singularity.
DO 330 NZ=NZ2,N
IF (JCOST.LE.(NZ-1)**2) GO TO 420
IJFIR = IFIRST(NZ)
IF (IJFIR) 230, 220, 240
220 IF (LL.EQ.1) NZMIN = NZ + 1
GO TO 330
230 LL = 2
IJFIR = -IJFIR
GO TO 290
240 LL = 2
c scan rows with nz non-zeros.
DO 270 IDUMMY=1,N
IF (JCOST.LE.(NZ-1)**2) GO TO 420
IF (ISRCH.GT.MSRCH) GO TO 420
IF (IJFIR.EQ.0) GO TO 280
c row ijfir is now examined.
I = IJFIR
IJFIR = NEXTR(I)
c first calculate multiplier threshold level.
AMAX = ZERO
J1 = IPTR(I) + LENRL(I)
J2 = IPTR(I) + LENR(I) - 1
DO 250 JJ=J1,J2
AMAX = DMAX1(AMAX,DABS(A(JJ)))
250 CONTINUE
AU = AMAX*U
ISRCH = ISRCH + 1
c scan row for possible pivots
DO 260 JJ=J1,J2
IF (DABS(A(JJ)).LE.AU .AND. L.EQ.1) GO TO 260
J = ICN(JJ)
KCOST = (NZ-1)*(LENC(J)-1)
IF (KCOST.GT.JCOST) GO TO 260
PIVR = ZERO
IF (AMAX.NE.ZERO) PIVR = DABS(A(JJ))/AMAX
IF (KCOST.EQ.JCOST .AND. (PIVR.LE.PIVRAT .OR.
* NSRCH.GT.NN+1)) GO TO 260
c best pivot so far is found.
JCOST = KCOST
IJPOS = JJ
IPIV = I
JPIV = J
IF (MSRCH.GT.NN+1 .AND. JCOST.LE.(NZ-1)**2) GO TO 420
PIVRAT = PIVR
260 CONTINUE
270 CONTINUE
c
c columns with nz non-zeros now examined.
280 IJFIR = IFIRST(NZ)
IJFIR = -LASTR(IJFIR)
290 IF (JCOST.LE.NZ*(NZ-1)) GO TO 420
IF (MSRCH.LE.NN) GO TO 330
DO 320 IDUMMY=1,N
IF (IJFIR.EQ.0) GO TO 330
J = IJFIR
IJFIR = NEXTC(IJFIR)
I1 = IPC(J)
I2 = I1 + NZ - 1
c scan column j.
DO 310 II=I1,I2
I = IRN(II)
KCOST = (NZ-1)*(LENR(I)-LENRL(I)-1)
IF (KCOST.GE.JCOST) GO TO 310
c pivot has best markowitz count so far ... now check its
c suitability on numeric grounds by examining the other non-zeros
c in its row.
J1 = IPTR(I) + LENRL(I)
J2 = IPTR(I) + LENR(I) - 1
c we need a stability check on singleton columns because of possible
c problems with underdetermined systems.
AMAX = ZERO
DO 300 JJ=J1,J2
AMAX = DMAX1(AMAX,DABS(A(JJ)))
IF (ICN(JJ).EQ.J) JPOS = JJ
300 CONTINUE
IF (DABS(A(JPOS)).LE.AMAX*U .AND. L.EQ.1) GO TO 310
JCOST = KCOST
IPIV = I
JPIV = J
IJPOS = JPOS
IF (AMAX.NE.ZERO) PIVRAT = DABS(A(JPOS))/AMAX
IF (JCOST.LE.NZ*(NZ-1)) GO TO 420
310 CONTINUE
c
320 CONTINUE
c
330 CONTINUE
c in the event of singularity, we must make sure all rows and columns
c are tested.
MSRCH = N
c
c matrix is numerically or structurally singular ... which it is will
c be diagnosed later.
IRANK = IRANK - 1
340 CONTINUE
c assign rest of rows and columns to ordering array.
c matrix is structurally singular.
IF (IFLAG.NE.2 .AND. IFLAG.NE.-5) IFLAG = 1
IRANK = IRANK - ILAST + PIVOT + 1
IF (.NOT.ABORT1) GO TO 350
IDISP(2) = IACTIV
IFLAG = -1
IF (LP.NE.0) WRITE (LP,99999)
GO TO 1120
350 K = PIVOT - 1
DO 390 I=ISTART,ILAST
IF (LASTR(I).NE.0) GO TO 390
K = K + 1
LASTR(I) = K
IF (LENRL(I).EQ.0) GO TO 380
MINICN = MAX0(MINICN,NZROW+IBEG-1+MOREI+LENRL(I))
IF (IACTIV-IBEG.GE.LENRL(I)) GO TO 360
CALL MA30DD(A, ICN, IPTR(ISTART), N, IACTIV, ITOP, .TRUE.)
c check now to see if ma30d/dd has created enough available space.
IF (IACTIV-IBEG.GE.LENRL(I)) GO TO 360
c create more space by destroying previously created lu factors.
MOREI = MOREI + IBEG - IDISP(1)
IBEG = IDISP(1)
IF (LP.NE.0) WRITE (LP,99997)
IFLAG = -5
IF (ABORT3) GO TO 1090
360 J1 = IPTR(I)
J2 = J1 + LENRL(I) - 1
IPTR(I) = 0
DO 370 JJ=J1,J2
A(IBEG) = A(JJ)
ICN(IBEG) = ICN(JJ)
ICN(JJ) = 0
IBEG = IBEG + 1
370 CONTINUE
NZROW = NZROW - LENRL(I)
380 IF (K.EQ.ILAST) GO TO 400
390 CONTINUE
400 K = PIVOT - 1
DO 410 I=ISTART,ILAST
IF (IPC(I).NE.0) GO TO 410
K = K + 1
IPC(I) = K
IF (K.EQ.ILAST) GO TO 990
410 CONTINUE
c
c the pivot has now been found in position (ipiv,jpiv) in location
c ijpos in row file.
c update column and row ordering arrays to correspond with removal
c of the active part of the matrix.
420 ISING = PIVOT
IF (A(IJPOS).NE.ZERO) GO TO 430
c numerical singularity is recorded here.
ISING = -ISING
IF (IFLAG.NE.-5) IFLAG = 2
IF (.NOT.ABORT2) GO TO 430
IDISP(2) = IACTIV
IFLAG = -2
IF (LP.NE.0) WRITE (LP,99998)
GO TO 1120
430 OLDPIV = IPTR(IPIV) + LENRL(IPIV)
OLDEND = IPTR(IPIV) + LENR(IPIV) - 1
c changes to column ordering.
IF (NSRCH.LE.NN) GO TO 460
COLUPD = NN + 1
LENPP = OLDEND-OLDPIV+1
IF (LENPP.LT.4) LPIV(1) = LPIV(1) + 1
IF (LENPP.GE.4 .AND. LENPP.LE.6) LPIV(2) = LPIV(2) + 1
IF (LENPP.GE.7 .AND. LENPP.LE.10) LPIV(3) = LPIV(3) + 1
IF (LENPP.GE.11 .AND. LENPP.LE.15) LPIV(4) = LPIV(4) + 1
IF (LENPP.GE.16 .AND. LENPP.LE.20) LPIV(5) = LPIV(5) + 1
IF (LENPP.GE.21 .AND. LENPP.LE.30) LPIV(6) = LPIV(6) + 1
IF (LENPP.GE.31 .AND. LENPP.LE.50) LPIV(7) = LPIV(7) + 1
IF (LENPP.GE.51 .AND. LENPP.LE.70) LPIV(8) = LPIV(8) + 1
IF (LENPP.GE.71 .AND. LENPP.LE.100) LPIV(9) = LPIV(9) + 1
IF (LENPP.GE.101) LPIV(10) = LPIV(10) + 1
MAPIV = MAX0(MAPIV,LENPP)
IAVPIV = IAVPIV + LENPP
DO 450 JJ=OLDPIV,OLDEND
J = ICN(JJ)
LC = LASTC(J)
NC = NEXTC(J)
NEXTC(J) = -COLUPD
IF (JJ.NE.IJPOS) COLUPD = J
IF (NC.NE.0) LASTC(NC) = LC
IF (LC.EQ.0) GO TO 440
NEXTC(LC) = NC
GO TO 450
440 NZ = LENC(J)
ISW = IFIRST(NZ)
IF (ISW.GT.0) LASTR(ISW) = -NC
IF (ISW.LT.0) IFIRST(NZ) = -NC
450 CONTINUE
c changes to row ordering.
460 I1 = IPC(JPIV)
I2 = I1 + LENC(JPIV) - 1
DO 480 II=I1,I2
I = IRN(II)
LR = LASTR(I)
NR = NEXTR(I)
IF (NR.NE.0) LASTR(NR) = LR
IF (LR.LE.0) GO TO 470
NEXTR(LR) = NR
GO TO 480
470 NZ = LENR(I) - LENRL(I)
IF (NR.NE.0) IFIRST(NZ) = NR
IF (NR.EQ.0) IFIRST(NZ) = LR
480 CONTINUE
c
c move pivot to position lenrl+1 in pivot row and move pivot row
c to the beginning of the available storage.
c the l part and the pivot in the old copy of the pivot row is
c nullified while, in the strictly upper triangular part, the
c column indices, j say, are overwritten by the corresponding
c entry of iq (iq(j)) and iq(j) is set to the negative of the
c displacement of the column index from the pivot entry.
IF (OLDPIV.EQ.IJPOS) GO TO 490
AU = A(OLDPIV)
A(OLDPIV) = A(IJPOS)
A(IJPOS) = AU
ICN(IJPOS) = ICN(OLDPIV)
ICN(OLDPIV) = JPIV
c check to see if there is space immediately available in a/icn to
c hold new copy of pivot row.
490 MINICN = MAX0(MINICN,NZROW+IBEG-1+MOREI+LENR(IPIV))
IF (IACTIV-IBEG.GE.LENR(IPIV)) GO TO 500
CALL MA30DD(A, ICN, IPTR(ISTART), N, IACTIV, ITOP, .TRUE.)
OLDPIV = IPTR(IPIV) + LENRL(IPIV)
OLDEND = IPTR(IPIV) + LENR(IPIV) - 1
c check now to see if ma30d/dd has created enough available space.
IF (IACTIV-IBEG.GE.LENR(IPIV)) GO TO 500
c create more space by destroying previously created lu factors.
MOREI = MOREI + IBEG - IDISP(1)
IBEG = IDISP(1)
IF (LP.NE.0) WRITE (LP,99997)
IFLAG = -5
IF (ABORT3) GO TO 1090
IF (IACTIV-IBEG.GE.LENR(IPIV)) GO TO 500
c there is still not enough room in a/icn.
IFLAG = -4
GO TO 1090
c copy pivot row and set up iq array.
500 IJPOS = 0
J1 = IPTR(IPIV)
c
DO 530 JJ=J1,OLDEND
A(IBEG) = A(JJ)
ICN(IBEG) = ICN(JJ)
IF (IJPOS.NE.0) GO TO 510
IF (ICN(JJ).EQ.JPIV) IJPOS = IBEG
ICN(JJ) = 0
GO TO 520
510 K = IBEG - IJPOS
J = ICN(JJ)
ICN(JJ) = IQ(J)
IQ(J) = -K
520 IBEG = IBEG + 1
530 CONTINUE
c
IJP1 = IJPOS + 1
PIVEND = IBEG - 1
LENPIV = PIVEND - IJPOS
NZROW = NZROW - LENRL(IPIV) - 1
IPTR(IPIV) = OLDPIV + 1
IF (LENPIV.EQ.0) IPTR(IPIV) = 0
c
c remove pivot row (including pivot) from column oriented file.
DO 560 JJ=IJPOS,PIVEND
J = ICN(JJ)
I1 = IPC(J)
LENC(J) = LENC(J) - 1
c i2 is last position in new column.
I2 = IPC(J) + LENC(J) - 1
IF (I2.LT.I1) GO TO 550
DO 540 II=I1,I2
IF (IRN(II).NE.IPIV) GO TO 540
IRN(II) = IRN(I2+1)
GO TO 550
540 CONTINUE
550 IRN(I2+1) = 0
560 CONTINUE
NZCOL = NZCOL - LENPIV - 1
c
c go down the pivot column and for each row with a non-zero add
c the appropriate multiple of the pivot row to it.
c we loop on the number of non-zeros in the pivot column since
c ma30d/dd may change its actual position.
c
NZPC = LENC(JPIV)
IF (NZPC.EQ.0) GO TO 900
DO 840 III=1,NZPC
II = IPC(JPIV) + III - 1
I = IRN(II)
c search row i for non-zero to be eliminated, calculate multiplier,
c and place it in position lenrl+1 in its row.
c idrop is the number of non-zero entries dropped from row i
c because these fall beneath tolerance level.
c
IDROP = 0
J1 = IPTR(I) + LENRL(I)
IEND = IPTR(I) + LENR(I) - 1
DO 570 JJ=J1,IEND
IF (ICN(JJ).NE.JPIV) GO TO 570
c if pivot is zero, rest of column is and so multiplier is zero.
AU = ZERO
IF (A(IJPOS).NE.ZERO) AU = -A(JJ)/A(IJPOS)
IF (LBIG) BIG = DMAX1(BIG,DABS(AU))
A(JJ) = A(J1)
A(J1) = AU
ICN(JJ) = ICN(J1)
ICN(J1) = JPIV
LENRL(I) = LENRL(I) + 1
GO TO 580
570 CONTINUE
c jump if pivot row is a singleton.
580 IF (LENPIV.EQ.0) GO TO 840
c now perform necessary operations on rest of non-pivot row i.
ROWI = J1 + 1
IOP = 0
c jump if all the pivot row causes fill-in.
IF (ROWI.GT.IEND) GO TO 650
c perform operations on current non-zeros in row i.
c innermost loop.
LENPP = IEND-ROWI+1
IF (LENPP.LT.4) LNPIV(1) = LNPIV(1) + 1
IF (LENPP.GE.4 .AND. LENPP.LE.6) LNPIV(2) = LNPIV(2) + 1
IF (LENPP.GE.7 .AND. LENPP.LE.10) LNPIV(3) = LNPIV(3) + 1
IF (LENPP.GE.11 .AND. LENPP.LE.15) LNPIV(4) = LNPIV(4) + 1
IF (LENPP.GE.16 .AND. LENPP.LE.20) LNPIV(5) = LNPIV(5) + 1
IF (LENPP.GE.21 .AND. LENPP.LE.30) LNPIV(6) = LNPIV(6) + 1
IF (LENPP.GE.31 .AND. LENPP.LE.50) LNPIV(7) = LNPIV(7) + 1
IF (LENPP.GE.51 .AND. LENPP.LE.70) LNPIV(8) = LNPIV(8) + 1
IF (LENPP.GE.71 .AND. LENPP.LE.100) LNPIV(9) = LNPIV(9) + 1
IF (LENPP.GE.101) LNPIV(10) = LNPIV(10) + 1
MANPIV = MAX0(MANPIV,LENPP)
IANPIV = IANPIV + LENPP
KOUNTL = KOUNTL + 1
DO 590 JJ=ROWI,IEND
J = ICN(JJ)
IF (IQ(J).GT.0) GO TO 590
IOP = IOP + 1
PIVROW = IJPOS - IQ(J)
A(JJ) = A(JJ) + AU*A(PIVROW)
IF (LBIG) BIG = DMAX1(DABS(A(JJ)),BIG)
ICN(PIVROW) = -ICN(PIVROW)
IF (DABS(A(JJ)).LT.TOL) IDROP = IDROP + 1
590 CONTINUE
c
c jump if no non-zeros in non-pivot row have been removed
c because these are beneath the drop-tolerance tol.
c
IF (IDROP.EQ.0) GO TO 650
c
c run through non-pivot row compressing row so that only
c non-zeros greater than tol are stored. all non-zeros
c less than tol are also removed from the column structure.
c
JNEW = ROWI
DO 630 JJ=ROWI,IEND
IF (DABS(A(JJ)).LT.TOL) GO TO 600
A(JNEW) = A(JJ)
ICN(JNEW) = ICN(JJ)
JNEW = JNEW + 1
GO TO 630
c
c remove non-zero entry from column structure.
c
600 J = ICN(JJ)
I1 = IPC(J)
I2 = I1 + LENC(J) - 1
DO 610 II=I1,I2
IF (IRN(II).EQ.I) GO TO 620
610 CONTINUE
620 IRN(II) = IRN(I2)
IRN(I2) = 0
LENC(J) = LENC(J) - 1
IF (NSRCH.LE.NN) GO TO 630
c remove column from column chain and place in update chain.
IF (NEXTC(J).LT.0) GO TO 630
c jump if column already in update chain.
LC = LASTC(J)
NC = NEXTC(J)
NEXTC(J) = -COLUPD
COLUPD = J
IF (NC.NE.0) LASTC(NC) = LC
IF (LC.EQ.0) GO TO 622
NEXTC(LC) = NC
GO TO 630
622 NZ = LENC(J) + 1
ISW = IFIRST(NZ)
IF (ISW.GT.0) LASTR(ISW) = -NC
IF (ISW.LT.0) IFIRST(NZ) = -NC
630 CONTINUE
DO 640 JJ=JNEW,IEND
ICN(JJ) = 0
640 CONTINUE
c the value of idrop might be different from that calculated earlier
c because, we may now have dropped some non-zeros which were not
c modified by the pivot row.
IDROP = IEND + 1 - JNEW
IEND = JNEW - 1
LENR(I) = LENR(I) - IDROP
NZROW = NZROW - IDROP
NZCOL = NZCOL - IDROP
NDROP = NDROP + IDROP
650 IFILL = LENPIV - IOP
c jump is if there is no fill-in.
IF (IFILL.EQ.0) GO TO 750
c now for the fill-in.
MINICN = MAX0(MINICN,MOREI+IBEG-1+NZROW+IFILL+LENR(I))
c see if there is room for fill-in.
c get maximum space for row i in situ.
DO 660 JDIFF=1,IFILL
JNPOS = IEND + JDIFF
IF (JNPOS.GT.LICN) GO TO 670
IF (ICN(JNPOS).NE.0) GO TO 670
660 CONTINUE
c there is room for all the fill-in after the end of the row so it
c can be left in situ.
c next available space for fill-in.
IEND = IEND + 1
GO TO 750
c jmore spaces for fill-in are required in front of row.
670 JMORE = IFILL - JDIFF + 1
I1 = IPTR(I)
c we now look in front of the row to see if there is space for
c the rest of the fill-in.
DO 680 JDIFF=1,JMORE
JNPOS = I1 - JDIFF
IF (JNPOS.LT.IACTIV) GO TO 690
IF (ICN(JNPOS).NE.0) GO TO 700
680 CONTINUE
690 JNPOS = I1 - JMORE
GO TO 710
c whole row must be moved to the beginning of available storage.
700 JNPOS = IACTIV - LENR(I) - IFILL
c jump if there is space immediately available for the shifted row.
710 IF (JNPOS.GE.IBEG) GO TO 730
CALL MA30DD(A, ICN, IPTR(ISTART), N, IACTIV, ITOP, .TRUE.)
I1 = IPTR(I)
IEND = I1 + LENR(I) - 1
JNPOS = IACTIV - LENR(I) - IFILL
IF (JNPOS.GE.IBEG) GO TO 730
c no space available so try to create some by throwing away previous
c lu decomposition.
MOREI = MOREI + IBEG - IDISP(1) - LENPIV - 1
IF (LP.NE.0) WRITE (LP,99997)
IFLAG = -5
IF (ABORT3) GO TO 1090
c keep record of current pivot row.
IBEG = IDISP(1)
ICN(IBEG) = JPIV
A(IBEG) = A(IJPOS)
IJPOS = IBEG
DO 720 JJ=IJP1,PIVEND
IBEG = IBEG + 1
A(IBEG) = A(JJ)
ICN(IBEG) = ICN(JJ)
720 CONTINUE
IJP1 = IJPOS + 1
PIVEND = IBEG
IBEG = IBEG + 1
IF (JNPOS.GE.IBEG) GO TO 730
c this still does not give enough room.
IFLAG = -4
GO TO 1090
730 IACTIV = MIN0(IACTIV,JNPOS)
c move non-pivot row i.
IPTR(I) = JNPOS
DO 740 JJ=I1,IEND
A(JNPOS) = A(JJ)
ICN(JNPOS) = ICN(JJ)
JNPOS = JNPOS + 1
ICN(JJ) = 0
740 CONTINUE
c first new available space.
IEND = JNPOS
750 NZROW = NZROW + IFILL
c innermost fill-in loop which also resets icn.
IDROP = 0
DO 830 JJ=IJP1,PIVEND
J = ICN(JJ)
IF (J.LT.0) GO TO 820
ANEW = AU*A(JJ)
AANEW = DABS(ANEW)
IF (AANEW.GE.TOL) GO TO 760
IDROP = IDROP + 1
NDROP = NDROP + 1
NZROW = NZROW - 1
MINICN = MINICN - 1
IFILL = IFILL - 1
GO TO 830
760 IF (LBIG) BIG = DMAX1(AANEW,BIG)
A(IEND) = ANEW
ICN(IEND) = J
IEND = IEND + 1
c
c put new entry in column file.
MINIRN = MAX0(MINIRN,NZCOL+LENC(J)+1)
JEND = IPC(J) + LENC(J)
JROOM = NZPC - III + 1 + LENC(J)
IF (JEND.GT.LIRN) GO TO 770
IF (IRN(JEND).EQ.0) GO TO 810
770 IF (JROOM.LT.DISPC) GO TO 780
c compress column file to obtain space for new copy of column.
CALL MA30DD(A, IRN, IPC(ISTART), N, DISPC, LIRN, .FALSE.)
IF (JROOM.LT.DISPC) GO TO 780
JROOM = DISPC - 1
IF (JROOM.GE.LENC(J)+1) GO TO 780
c column file is not large enough.
GO TO 1100
c copy column to beginning of file.
780 JBEG = IPC(J)
JEND = IPC(J) + LENC(J) - 1
JZERO = DISPC - 1
DISPC = DISPC - JROOM
IDISPC = DISPC
DO 790 II=JBEG,JEND
IRN(IDISPC) = IRN(II)
IRN(II) = 0
IDISPC = IDISPC + 1
790 CONTINUE
IPC(J) = DISPC
JEND = IDISPC
DO 800 II=JEND,JZERO
IRN(II) = 0
800 CONTINUE
810 IRN(JEND) = I
NZCOL = NZCOL + 1
LENC(J) = LENC(J) + 1
c end of adjustment to column file.
GO TO 830
c
820 ICN(JJ) = -J
830 CONTINUE
IF (IDROP.EQ.0) GO TO 834
DO 832 KDROP=1,IDROP
ICN(IEND) = 0
IEND = IEND + 1
832 CONTINUE
834 LENR(I) = LENR(I) + IFILL
c end of scan of pivot column.
840 CONTINUE
c
c
c remove pivot column from column oriented storage and update row
c ordering arrays.
I1 = IPC(JPIV)
I2 = IPC(JPIV) + LENC(JPIV) - 1
NZCOL = NZCOL - LENC(JPIV)
DO 890 II=I1,I2
I = IRN(II)
IRN(II) = 0
NZ = LENR(I) - LENRL(I)
IF (NZ.NE.0) GO TO 850
LASTR(I) = 0
GO TO 890
850 IFIR = IFIRST(NZ)
IFIRST(NZ) = I
IF (IFIR) 860, 880, 870
860 LASTR(I) = IFIR
NEXTR(I) = 0
GO TO 890
870 LASTR(I) = LASTR(IFIR)
NEXTR(I) = IFIR
LASTR(IFIR) = I
GO TO 890
880 LASTR(I) = 0
NEXTR(I) = 0
NZMIN = MIN0(NZMIN,NZ)
890 CONTINUE
c restore iq and nullify u part of old pivot row.
c record the column permutation in lastc(jpiv) and the row
c permutation in lastr(ipiv).
900 IPC(JPIV) = -ISING
LASTR(IPIV) = PIVOT
IF (LENPIV.EQ.0) GO TO 980
NZROW = NZROW - LENPIV
JVAL = IJP1
JZER = IPTR(IPIV)
IPTR(IPIV) = 0
DO 910 JCOUNT=1,LENPIV
J = ICN(JVAL)
IQ(J) = ICN(JZER)
ICN(JZER) = 0
JVAL = JVAL + 1
JZER = JZER + 1
910 CONTINUE
c adjust column ordering arrays.
IF (NSRCH.GT.NN) GO TO 920
DO 916 JJ=IJP1,PIVEND
J = ICN(JJ)
NZ = LENC(J)
IF (NZ.NE.0) GO TO 914
IPC(J) = 0
GO TO 916
914 NZMIN = MIN0(NZMIN,NZ)
916 CONTINUE
GO TO 980
920 JJ = COLUPD
DO 970 JDUMMY=1,NN
J = JJ
IF (J.EQ.NN+1) GO TO 980
JJ = -NEXTC(J)
NZ = LENC(J)
IF (NZ.NE.0) GO TO 924
IPC(J) = 0
GO TO 970
924 IFIR = IFIRST(NZ)
LASTC(J) = 0
IF (IFIR) 930, 940, 950
930 IFIRST(NZ) = -J
IFIR = -IFIR
LASTC(IFIR) = J
NEXTC(J) = IFIR
GO TO 970
940 IFIRST(NZ) = -J
NEXTC(J) = 0
GO TO 960
950 LC = -LASTR(IFIR)
LASTR(IFIR) = -J
NEXTC(J) = LC
IF (LC.NE.0) LASTC(LC) = J
960 NZMIN = MIN0(NZMIN,NZ)
970 CONTINUE
980 CONTINUE
c ********************************************
c **** end of main elimination loop ****
c ********************************************
c
c reset iactiv to point to the beginning of the next block.
990 IF (ILAST.NE.NN) IACTIV = IPTR(ILAST+1)
1000 CONTINUE
c
c ********************************************
c **** end of deomposition of block ****
c ********************************************
c
c record singularity (if any) in iq array.
IF (IRANK.EQ.NN) GO TO 1020
DO 1010 I=1,NN
IF (IPC(I).LT.0) GO TO 1010
ISING = IPC(I)
IQ(ISING) = -IQ(ISING)
IPC(I) = -ISING
1010 CONTINUE
c
c run through lu decomposition changing column indices to that of new
c order and permuting lenr and lenrl arrays according to pivot
c permutations.
1020 ISTART = IDISP(1)
IEND = IBEG - 1
IF (IEND.LT.ISTART) GO TO 1040
DO 1030 JJ=ISTART,IEND
JOLD = ICN(JJ)
ICN(JJ) = -IPC(JOLD)
1030 CONTINUE
1040 DO 1050 II=1,NN
I = LASTR(II)
NEXTR(I) = LENR(II)
IPTR(I) = LENRL(II)
1050 CONTINUE
DO 1060 I=1,NN
LENRL(I) = IPTR(I)
LENR(I) = NEXTR(I)
1060 CONTINUE
c
c update permutation arrays ip and iq.
DO 1070 II=1,NN
I = LASTR(II)
J = -IPC(II)
NEXTR(I) = IABS(IP(II)+0)
IPTR(J) = IABS(IQ(II)+0)
1070 CONTINUE
DO 1080 I=1,NN
IF (IP(I).LT.0) NEXTR(I) = -NEXTR(I)
IP(I) = NEXTR(I)
IF (IQ(I).LT.0) IPTR(I) = -IPTR(I)
IQ(I) = IPTR(I)
1080 CONTINUE
IP(NN) = IABS(IP(NN)+0)
IDISP(2) = IEND
GO TO 1120
c
c *** error returns ***
1090 IDISP(2) = IACTIV
IF (LP.EQ.0) GO TO 1120
WRITE (LP,99996)
GO TO 1110
1100 IF (IFLAG.EQ.-5) IFLAG = -6
IF (IFLAG.NE.-6) IFLAG = -3
IDISP(2) = IACTIV
IF (LP.EQ.0) GO TO 1120
IF (IFLAG.EQ.-3) WRITE (LP,99995)
IF (IFLAG.EQ.-6) WRITE (LP,99994)
1110 PIVOT = PIVOT - ISTART + 1
WRITE (LP,99993) PIVOT, NBLOCK, ISTART, ILAST
IF (PIVOT.EQ.0) WRITE (LP,99992) MINIRN
c
c
1120 RETURN
99999 FORMAT (54H ERROR RETURN FROM MA30A/AD BECAUSE MATRIX IS STRUCTUR,
* 13HALLY SINGULAR)
99998 FORMAT (54H ERROR RETURN FROM MA30A/AD BECAUSE MATRIX IS NUMERICA,
* 12HLLY SINGULAR)
99997 FORMAT (48H LU DECOMPOSITION DESTROYED TO CREATE MORE SPACE)
99996 FORMAT (54H ERROR RETURN FROM MA30A/AD BECAUSE LICN NOT BIG ENOUG,
* 1HH)
99995 FORMAT (54H ERROR RETURN FROM MA30A/AD BECAUSE LIRN NOT BIG ENOUG,
* 1HH)
99994 FORMAT (51H ERROR RETURN FROM MA30A/AD LIRN AND LICN TOO SMALL)
99993 FORMAT (10H AT STAGE , I5, 10H IN BLOCK , I5, 16H WITH FIRST ROW ,
* I5, 14H AND LAST ROW , I5)
99992 FORMAT (34H TO CONTINUE SET LIRN TO AT LEAST , I8)
END
c
c=====================================================================
c
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
c
c======================================================================
c
SUBROUTINE MA30CD(N, ICN, A, LICN, LENR, LENRL, LENOFF, IDISP, IP,
* IQ, X, W, MTYPE)
c ma30c/cd uses the factors produced by ma30a/ad or ma30b/bd to solve
c ax=b or a transpose x=b when the matrix p1*a*q1 (paq) is block
c lower triangular (including the case of only one diagonal
c block).
c
c we now describe the argument list for ma30c/cd.
c n is an integer variable set to the order of the matrix. it is not
c altered by the subroutine.
c icn is an integer array of length licn. entries idisp(1) to
c idisp(2) should be unchanged since the last call to ma30a/ad. if
c the matrix has more than one diagonal block, then column indices
c corresponding to non-zeros in sub-diagonal blocks of paq must
c appear in positions 1 to idisp(1)-1. for the same row those
c entries must be contiguous, with those in row i preceding those
c in row i+1 (i=1,...,n-1) and no wasted space between rows.
c entries may be in any order within each row. it is not altered
c by ma30c/cd.
c a is a real/double precision array of length licn. entries
c idisp(1) to idisp(2) should be unchanged since the last call to
c ma30a/ad or ma30b/bd. if the matrix has more than one diagonal
c block, then the values of the non-zeros in sub-diagonal blocks
c must be in positions 1 to idisp(1)-1 in the order given by icn.
c it is not altered by ma30c/cd.
c licn is an integer variable set to the size of arrays icn and a.
c it is not altered by ma30c/cd.
c lenr,lenrl are integer arrays of length n which should be
c unchanged since the last call to ma30a/ad. they are not altered
c by ma30c/cd.
c lenoff is an integer array of length n. if the matrix paq (or
c p1*a*q1) has more than one diagonal block, then lenoff(i),
c i=1,...,n should be set to the number of non-zeros in row i of
c the matrix paq which are in sub-diagonal blocks. if there is
c only one diagonal block then lenoff(1) may be set to -1, in
c which case the other entries of lenoff are never accessed. it is
c not altered by ma30c/cd.
c idisp is an integer array of length 2 which should be unchanged
c since the last call to ma30a/ad. it is not altered by ma30c/cd.
c ip,iq are integer arrays of length n which should be unchanged
c since the last call to ma30a/ad. they are not altered by
c ma30c/cd.
c x is a real/double precision array of length n. it must be set by
c the user to the values of the right hand side vector b for the
c equations being solved. on exit from ma30c/cd it will be equal
c to the solution x required.
c w is a real/double precision array of length n which is used as
c workspace by ma30c/cd.
c mtype is an integer variable which must be set by the user. if
c mtype=1, then the solution to the system ax=b is returned; any
c other value for mtype will return the solution to the system a
c transpose x=b. it is not altered by ma30c/cd.
c
DOUBLE PRECISION A(LICN), X(N), W(N), WII, WI, RESID, ZERO
LOGICAL NEG, NOBLOC
INTEGER IDISP(2)
INTEGER ICN(LICN), LENR(N), LENRL(N), LENOFF(N), IP(N), IQ(N)
c see block data for comments on variables in common.
COMMON /MA30HD/ RESID
DATA ZERO /0.0D0/
c
c the final value of resid is the maximum residual for an inconsistent
c set of equations.
RESID = ZERO
c nobloc is .true. if subroutine block has been used previously and
c is .false. otherwise. the value .false. means that lenoff
c will not be subsequently accessed.
NOBLOC = LENOFF(1).LT.0
IF (MTYPE.NE.1) GO TO 140
c
c we now solve a * x = b.
c neg is used to indicate when the last row in a block has been
c reached. it is then set to true whereafter backsubstitution is
c performed on the block.
NEG = .FALSE.
c ip(n) is negated so that the last row of the last block can be
c recognised. it is reset to its positive value on exit.
IP(N) = -IP(N)
c preorder vector ... w(i) = x(ip(i))
DO 10 II=1,N
I = IP(II)
I = IABS(I)
W(II) = X(I)
10 CONTINUE
c lt holds the position of the first non-zero in the current row of the
c off-diagonal blocks.
LT = 1
c ifirst holds the index of the first row in the current block.
IFIRST = 1
c iblock holds the position of the first non-zero in the current row
c of the lu decomposition of the diagonal blocks.
IBLOCK = IDISP(1)
c if i is not the last row of a block, then a pass through this loop
c adds the inner product of row i of the off-diagonal blocks and w
c to w and performs forward elimination using row i of the lu
c decomposition. if i is the last row of a block then, after
c performing these aforementioned operations, backsubstitution is
c performed using the rows of the block.
DO 120 I=1,N
WI = W(I)
IF (NOBLOC) GO TO 30
IF (LENOFF(I).EQ.0) GO TO 30
c operations using lower triangular blocks.
c ltend is the end of row i in the off-diagonal blocks.
LTEND = LT + LENOFF(I) - 1
DO 20 JJ=LT,LTEND
J = ICN(JJ)
WI = WI - A(JJ)*W(J)
20 CONTINUE
c lt is set the beginning of the next off-diagonal row.
LT = LTEND + 1
c set neg to .true. if we are on the last row of the block.
30 IF (IP(I).LT.0) NEG = .TRUE.
IF (LENRL(I).EQ.0) GO TO 50
c forward elimination phase.
c iend is the end of the l part of row i in the lu decomposition.
IEND = IBLOCK + LENRL(I) - 1
DO 40 JJ=IBLOCK,IEND
J = ICN(JJ)
WI = WI + A(JJ)*W(J)
40 CONTINUE
c iblock is adjusted to point to the start of the next row.
50 IBLOCK = IBLOCK + LENR(I)
W(I) = WI
IF (.NOT.NEG) GO TO 120
c back substitution phase.
c j1 is position in a/icn after end of block beginning in row ifirst
c and ending in row i.
J1 = IBLOCK
c are there any singularities in this block? if not, continue with
c the backsubstitution.
IB = I
IF (IQ(I).GT.0) GO TO 70
DO 60 III=IFIRST,I
IB = I - III + IFIRST
IF (IQ(IB).GT.0) GO TO 70
J1 = J1 - LENR(IB)
RESID = DMAX1(RESID,DABS(W(IB)))
W(IB) = ZERO
60 CONTINUE
c entire block is singular.
GO TO 110
c each pass through this loop performs the back-substitution
c operations for a single row, starting at the end of the block and
c working through it in reverse order.
70 DO 100 III=IFIRST,IB
II = IB - III + IFIRST
c j2 is end of row ii.
J2 = J1 - 1
c j1 is beginning of row ii.
J1 = J1 - LENR(II)
c jpiv is the position of the pivot in row ii.
JPIV = J1 + LENRL(II)
JPIVP1 = JPIV + 1
c jump if row ii of u has no non-zeros.
IF (J2.LT.JPIVP1) GO TO 90
WII = W(II)
DO 80 JJ=JPIVP1,J2
J = ICN(JJ)
WII = WII - A(JJ)*W(J)
80 CONTINUE
W(II) = WII
90 W(II) = W(II)/A(JPIV)
100 CONTINUE
110 IFIRST = I + 1
NEG = .FALSE.
120 CONTINUE
c
c reorder solution vector ... x(i) = w(iqinverse(i))
DO 130 II=1,N
I = IQ(II)
I = IABS(I)
X(I) = W(II)
130 CONTINUE
IP(N) = -IP(N)
GO TO 320
c
c
c we now solve atranspose * x = b.
c preorder vector ... w(i)=x(iq(i))
140 DO 150 II=1,N
I = IQ(II)
I = IABS(I)
W(II) = X(I)
150 CONTINUE
c lj1 points to the beginning the current row in the off-diagonal
c blocks.
LJ1 = IDISP(1)
c iblock is initialized to point to the beginning of the block after
c the last one ]
IBLOCK = IDISP(2) + 1
c ilast is the last row in the current block.
ILAST = N
c iblend points to the position after the last non-zero in the
c current block.
IBLEND = IBLOCK
c each pass through this loop operates with one diagonal block and
c the off-diagonal part of the matrix corresponding to the rows
c of this block. the blocks are taken in reverse order and the
c number of times the loop is entered is min(n,no. blocks+1).
DO 290 NUMBLK=1,N
IF (ILAST.EQ.0) GO TO 300
IBLOCK = IBLOCK - LENR(ILAST)
c this loop finds the index of the first row in the current block..
c it is first and iblock is set to the position of the beginning
c of this first row.
DO 160 K=1,N
II = ILAST - K
IF (II.EQ.0) GO TO 170
IF (IP(II).LT.0) GO TO 170
IBLOCK = IBLOCK - LENR(II)
160 CONTINUE
170 IFIRST = II + 1
c j1 points to the position of the beginning of row i (lt part) or pivot
J1 = IBLOCK
c forward elimination.
c each pass through this loop performs the operations for one row of the
c block. if the corresponding entry of w is zero then the
c operations can be avoided.
DO 210 I=IFIRST,ILAST
IF (W(I).EQ.ZERO) GO TO 200
c jump if row i singular.
IF (IQ(I).LT.0) GO TO 220
c j2 first points to the pivot in row i and then is made to point to the
c first non-zero in the u transpose part of the row.
J2 = J1 + LENRL(I)
WI = W(I)/A(J2)
IF (LENR(I)-LENRL(I).EQ.1) GO TO 190
J2 = J2 + 1
c j3 points to the end of row i.
J3 = J1 + LENR(I) - 1
DO 180 JJ=J2,J3
J = ICN(JJ)
W(J) = W(J) - A(JJ)*WI
180 CONTINUE
190 W(I) = WI
200 J1 = J1 + LENR(I)
210 CONTINUE
GO TO 240
c deals with rest of block which is singular.
220 DO 230 II=I,ILAST
RESID = DMAX1(RESID,DABS(W(II)))
W(II) = ZERO
230 CONTINUE
c back substitution.
c this loop does the back substitution on the rows of the block in
c the reverse order doing it simultaneously on the l transpose part
c of the diagonal blocks and the off-diagonal blocks.
240 J1 = IBLEND
DO 280 IBACK=IFIRST,ILAST
I = ILAST - IBACK + IFIRST
c j1 points to the beginning of row i.
J1 = J1 - LENR(I)
IF (LENRL(I).EQ.0) GO TO 260
c j2 points to the end of the l transpose part of row i.
J2 = J1 + LENRL(I) - 1
DO 250 JJ=J1,J2
J = ICN(JJ)
W(J) = W(J) + A(JJ)*W(I)
250 CONTINUE
260 IF (NOBLOC) GO TO 280
c operations using lower triangular blocks.
IF (LENOFF(I).EQ.0) GO TO 280
c lj2 points to the end of row i of the off-diagonal blocks.
LJ2 = LJ1 - 1
c lj1 points to the beginning of row i of the off-diagonal blocks.
LJ1 = LJ1 - LENOFF(I)
DO 270 JJ=LJ1,LJ2
J = ICN(JJ)
W(J) = W(J) - A(JJ)*W(I)
270 CONTINUE
280 CONTINUE
IBLEND = J1
ILAST = IFIRST - 1
290 CONTINUE
c reorder solution vector ... x(i)=w(ipinverse(i))
300 DO 310 II=1,N
I = IP(II)
I = IABS(I)
X(I) = W(II)
310 CONTINUE
c
320 RETURN
END
SUBROUTINE MA30DD(A, ICN, IPTR, N, IACTIV, ITOP, REALS)
c this subroutine performs garbage collection operations on the
c arrays a, icn and irn.
c iactiv is the first position in arrays a/icn from which the compress
c starts. on exit, iactiv equals the position of the first entry
c in the compressed part of a/icn
c
DOUBLE PRECISION A(ITOP)
LOGICAL REALS
INTEGER IPTR(N)
INTEGER ICN(ITOP)
c see block data for comments on variables in common.
COMMON /MA30FD/ IRNCP, ICNCP, IRANK, MINIRN, MINICN
c
IF (REALS) ICNCP = ICNCP + 1
IF (.NOT.REALS) IRNCP = IRNCP + 1
c set the first non-zero entry in each row to the negative of the
c row/col number and hold this row/col index in the row/col
c pointer. this is so that the beginning of each row/col can
c be recognized in the subsequent scan.
DO 10 J=1,N
K = IPTR(J)
IF (K.LT.IACTIV) GO TO 10
IPTR(J) = ICN(K)
ICN(K) = -J
10 CONTINUE
KN = ITOP + 1
KL = ITOP - IACTIV + 1
c go through arrays in reverse order compressing to the back so
c that there are no zeros held in positions iactiv to itop in icn.
c reset first entry of each row/col and pointer array iptr.
DO 30 K=1,KL
JPOS = ITOP - K + 1
IF (ICN(JPOS).EQ.0) GO TO 30
KN = KN - 1
IF (REALS) A(KN) = A(JPOS)
IF (ICN(JPOS).GE.0) GO TO 20
c first non-zero of row/col has been located
J = -ICN(JPOS)
ICN(JPOS) = IPTR(J)
IPTR(J) = KN
20 ICN(KN) = ICN(JPOS)
30 CONTINUE
IACTIV = KN
RETURN
END
c######date 01 jan 1984 copyright ukaea, harwell.
c######alias mc13d
SUBROUTINE MC13D(N,ICN,LICN,IP,LENR,IOR,IB,NUM,IW)
INTEGER IP(N)
INTEGER ICN(LICN),LENR(N),IOR(N),IB(N),IW(N,3)
CALL MC13E(N,ICN,LICN,IP,LENR,IOR,IB,NUM,IW(1,1),IW(1,2),IW(1,3))
RETURN
END
SUBROUTINE MC13E(N,ICN,LICN,IP,LENR,ARP,IB,NUM,LOWL,NUMB,PREV)
INTEGER STP,DUMMY
INTEGER IP(N)
c
c arp(i) is one less than the number of unsearched edges leaving
c node i. at the end of the algorithm it is set to a
c permutation which puts the matrix in block lower
c triangular form.
c ib(i) is the position in the ordering of the start of the ith
c block. ib(n+1-i) holds the node number of the ith node
c on the stack.
c lowl(i) is the smallest stack position of any node to which a path
c from node i has been found. it is set to n+1 when node i
c is removed from the stack.
c numb(i) is the position of node i in the stack if it is on
c it, is the permuted order of node i for those nodes
c whose final position has been found and is otherwise zero.
c prev(i) is the node at the end of the path when node i was
c placed on the stack.
INTEGER ICN(LICN),LENR(N),ARP(N),IB(N),LOWL(N),NUMB(N),
1PREV(N)
c
c
c icnt is the number of nodes whose positions in final ordering have
c been found.
ICNT=0
c num is the number of blocks that have been found.
NUM=0
NNM1=N+N-1
c
c initialization of arrays.
DO 20 J=1,N
NUMB(J)=0
ARP(J)=LENR(J)-1
20 CONTINUE
c
c
DO 120 ISN=1,N
c look for a starting node
IF (NUMB(ISN).NE.0) GO TO 120
IV=ISN
c ist is the number of nodes on the stack ... it is the stack pointer.
IST=1
c put node iv at beginning of stack.
LOWL(IV)=1
NUMB(IV)=1
IB(N)=IV
c
c the body of this loop puts a new node on the stack or backtracks.
DO 110 DUMMY=1,NNM1
I1=ARP(IV)
c have all edges leaving node iv been searched.
IF (I1.LT.0) GO TO 60
I2=IP(IV)+LENR(IV)-1
I1=I2-I1
c
c look at edges leaving node iv until one enters a new node or
c all edges are exhausted.
DO 50 II=I1,I2
IW=ICN(II)
c has node iw been on stack already.
IF (NUMB(IW).EQ.0) GO TO 100
c update value of lowl(iv) if necessary.
50 LOWL(IV)=MIN0(LOWL(IV),LOWL(IW))
c
c there are no more edges leaving node iv.
ARP(IV)=-1
c is node iv the root of a block.
60 IF (LOWL(IV).LT.NUMB(IV)) GO TO 90
c
c order nodes in a block.
NUM=NUM+1
IST1=N+1-IST
LCNT=ICNT+1
c peel block off the top of the stack starting at the top and
c working down to the root of the block.
DO 70 STP=IST1,N
IW=IB(STP)
LOWL(IW)=N+1
ICNT=ICNT+1
NUMB(IW)=ICNT
IF (IW.EQ.IV) GO TO 80
70 CONTINUE
80 IST=N-STP
IB(NUM)=LCNT
c are there any nodes left on the stack.
IF (IST.NE.0) GO TO 90
c have all the nodes been ordered.
IF (ICNT.LT.N) GO TO 120
GO TO 130
c
c backtrack to previous node on path.
90 IW=IV
IV=PREV(IV)
c update value of lowl(iv) if necessary.
LOWL(IV)=MIN0(LOWL(IV),LOWL(IW))
GO TO 110
c
c put new node on the stack.
100 ARP(IV)=I2-II-1
PREV(IW)=IV
IV=IW
IST=IST+1
LOWL(IV)=IST
NUMB(IV)=IST
K=N+1-IST
IB(K)=IV
110 CONTINUE
c
120 CONTINUE
c
c
c put permutation in the required form.
130 DO 140 I=1,N
II=NUMB(I)
140 ARP(II)=I
RETURN
END
c######date 01 jan 1984 copyright ukaea, harwell.
c######alias mc20ad mc20bd
SUBROUTINE MC20AD(NC,MAXA,A,INUM,JPTR,JNUM,JDISP)
c
INTEGER INUM(MAXA),JNUM(MAXA)
DOUBLE PRECISION A(MAXA),ACE,ACEP
DIMENSION JPTR(NC)
c
c ******************************************************************
c
NULL=-JDISP
c** clear jptr
DO 60 J=1,NC
60 JPTR(J)=0
c** count the number of elements in each column.
DO 120 K=1,MAXA
J=JNUM(K)+JDISP
JPTR(J)=JPTR(J)+1
120 CONTINUE
c** set the jptr array
K=1
DO 150 J=1,NC
KR=K+JPTR(J)
JPTR(J)=K
150 K=KR
c
c** reorder the elements into column order. the algorithm is an
c in-place sort and is of order maxa.
DO 230 I=1,MAXA
c establish the current entry.
JCE=JNUM(I)+JDISP
IF(JCE.EQ.0) GO TO 230
ACE=A(I)
ICE=INUM(I)
c clear the location vacated.
JNUM(I)=NULL
c chain from current entry to store items.
DO 200 J=1,MAXA
c current entry not in correct position. determine correct
c position to store entry.
LOC=JPTR(JCE)
JPTR(JCE)=JPTR(JCE)+1
c save contents of that location.
ACEP=A(LOC)
ICEP=INUM(LOC)
JCEP=JNUM(LOC)
c store current entry.
A(LOC)=ACE
INUM(LOC)=ICE
JNUM(LOC)=NULL
c check if next current entry needs to be processed.
IF(JCEP.EQ.NULL) GO TO 230
c it does. copy into current entry.
ACE=ACEP
ICE=ICEP
200 JCE=JCEP+JDISP
c
230 CONTINUE
c
c** reset jptr vector.
JA=1
DO 250 J=1,NC
JB=JPTR(J)
JPTR(J)=JA
250 JA=JB
RETURN
END
c######date 01 jan 1984 copyright ukaea, harwell.
c######alias mc21a
SUBROUTINE MC21A(N,ICN,LICN,IP,LENR,IPERM,NUMNZ,IW)
INTEGER IP(N)
INTEGER ICN(LICN),LENR(N),IPERM(N),IW(N,4)
CALL MC21B(N,ICN,LICN,IP,LENR,IPERM,NUMNZ,IW(1,1),IW(1,2),IW(1,3),
1IW(1,4))
RETURN
END
SUBROUTINE MC21B(N,ICN,LICN,IP,LENR,IPERM,NUMNZ,PR,ARP,CV,OUT)
INTEGER IP(N)
c pr(i) is the previous row to i in the depth first search.
c it is used as a work array in the sorting algorithm.
c elements (iperm(i),i) i=1, ... n are non-zero at the end of the
c algorithm unless n assignments have not been made. in which case
c (iperm(i),i) will be zero for n-numnz entries.
c cv(i) is the most recent row extension at which column i
c was visited.
c arp(i) is one less than the number of non-zeros in row i
c which have not been scanned when looking for a cheap assignment.
c out(i) is one less than the number of non-zeros in row i
c which have not been scanned during one pass through the main loop.
INTEGER ICN(LICN),LENR(N),IPERM(N),PR(N),CV(N),
1ARP(N),OUT(N)
c
c initialization of arrays.
DO 10 I=1,N
ARP(I)=LENR(I)-1
CV(I)=0
10 IPERM(I)=0
NUMNZ=0
c
c
c main loop.
c each pass round this loop either results in a new assignment
c or gives a row with no assignment.
DO 130 JORD=1,N
J=JORD
PR(J)=-1
DO 100 K=1,JORD
c look for a cheap assignment
IN1=ARP(J)
IF (IN1.LT.0) GO TO 60
IN2=IP(J)+LENR(J)-1
IN1=IN2-IN1
DO 50 II=IN1,IN2
I=ICN(II)
IF (IPERM(I).EQ.0) GO TO 110
50 CONTINUE
c no cheap assignment in row.
ARP(J)=-1
c begin looking for assignment chain starting with row j.
60 OUT(J)=LENR(J)-1
c inner loop. extends chain by one or backtracks.
DO 90 KK=1,JORD
IN1=OUT(J)
IF (IN1.LT.0) GO TO 80
IN2=IP(J)+LENR(J)-1
IN1=IN2-IN1
c forward scan.
DO 70 II=IN1,IN2
I=ICN(II)
IF (CV(I).EQ.JORD) GO TO 70
c column i has not yet been accessed during this pass.
J1=J
J=IPERM(I)
CV(I)=JORD
PR(J)=J1
OUT(J1)=IN2-II-1
GO TO 100
70 CONTINUE
c
c backtracking step.
80 J=PR(J)
IF (J.EQ.-1) GO TO 130
90 CONTINUE
c
100 CONTINUE
c
c new assignment is made.
110 IPERM(I)=J
ARP(J)=IN2-II-1
NUMNZ=NUMNZ+1
DO 120 K=1,JORD
J=PR(J)
IF (J.EQ.-1) GO TO 130
II=IP(J)+LENR(J)-OUT(J)-2
I=ICN(II)
IPERM(I)=J
120 CONTINUE
c
130 CONTINUE
c
c if matrix is structurally singular, we now complete the
c permutation iperm.
IF (NUMNZ.EQ.N) RETURN
DO 140 I=1,N
140 ARP(I)=0
K=0
DO 160 I=1,N
IF (IPERM(I).NE.0) GO TO 150
K=K+1
OUT(K)=I
GO TO 160
150 J=IPERM(I)
ARP(J)=I
160 CONTINUE
K=0
DO 170 I=1,N
IF (ARP(I).NE.0) GO TO 170
K=K+1
IOUTK=OUT(K)
IPERM(IOUTK)=I
170 CONTINUE
RETURN
END
c######date 01 jan 1984 copyright ukaea, harwell.
c######alias mc22ad
SUBROUTINE MC22AD(N,ICN,A,NZ,LENROW,IP,IQ,IW,IW1)
DOUBLE PRECISION A(NZ),AVAL
INTEGER IW(N,2)
INTEGER ICN(NZ),LENROW(N),IP(N),IQ(N),IW1(NZ)
IF (NZ.LE.0) GO TO 1000
IF (N.LE.0) GO TO 1000
c set start of row i in iw(i,1) and lenrow(i) in iw(i,2)
IW(1,1)=1
IW(1,2)=LENROW(1)
DO 10 I=2,N
IW(I,1)=IW(I-1,1)+LENROW(I-1)
10 IW(I,2)=LENROW(I)
c permute lenrow according to ip. set off-sets for new position
c of row iold in iw(iold,1) and put old row indices in iw1 in
c positions corresponding to the new position of this row in a/icn.
JJ=1
DO 20 I=1,N
IOLD=IP(I)
IOLD=IABS(IOLD)
LENGTH=IW(IOLD,2)
LENROW(I)=LENGTH
IF (LENGTH.EQ.0) GO TO 20
IW(IOLD,1)=IW(IOLD,1)-JJ
J2=JJ+LENGTH-1
DO 15 J=JJ,J2
15 IW1(J)=IOLD
JJ=J2+1
20 CONTINUE
c set inverse permutation to iq in iw(.,2).
DO 30 I=1,N
IOLD=IQ(I)
IOLD=IABS(IOLD)
30 IW(IOLD,2)=I
c permute a and icn in place, changing to new column numbers.
c
c *** main loop ***
c each pass through this loop places a closed chain of column indices
c in their new (and final) positions ... this is recorded by
c setting the iw1 entry to zero so that any which are subsequently
c encountered during this major scan can be bypassed.
DO 200 I=1,NZ
IOLD=IW1(I)
IF (IOLD.EQ.0) GO TO 200
IPOS=I
JVAL=ICN(I)
c if row iold is in same positions after permutation go to 150.
IF (IW(IOLD,1).EQ.0) GO TO 150
AVAL=A(I)
c ** chain loop **
c each pass through this loop places one (permuted) column index
c in its final position .. viz. ipos.
DO 100 ICHAIN=1,NZ
c newpos is the original position in a/icn of the element to be placed
c in position ipos. it is also the position of the next element in
c the chain.
NEWPOS=IPOS+IW(IOLD,1)
c is chain complete ?
IF (NEWPOS.EQ.I) GO TO 130
A(IPOS)=A(NEWPOS)
JNUM=ICN(NEWPOS)
ICN(IPOS)=IW(JNUM,2)
IPOS=NEWPOS
IOLD=IW1(IPOS)
IW1(IPOS)=0
c ** end of chain loop **
100 CONTINUE
130 A(IPOS)=AVAL
150 ICN(IPOS)=IW(JVAL,2)
c *** end of main loop ***
200 CONTINUE
c
1000 RETURN
END
c######date 01 jan 1984 copyright ukaea, harwell.
c######alias mc23ad
c###### calls mc13 mc21
SUBROUTINE MC23AD(N,ICN,A,LICN,LENR,IDISP,IP,IQ,LENOFF,IW,IW1)
DOUBLE PRECISION A(LICN)
INTEGER IDISP(2),IW1(N,2)
LOGICAL ABORT
INTEGER ICN(LICN),LENR(N),IP(N),IQ(N),LENOFF(N),IW(N,5)
COMMON /MC23BD/ LP,NUMNZ,NUM,LARGE,ABORT
c input ... n,icn .. a,icn,lenr ....
c
c set up pointers iw(.,1) to the beginning of the rows and set lenoff
c equal to lenr.
IW1(1,1)=1
LENOFF(1)=LENR(1)
IF (N.EQ.1) GO TO 20
DO 10 I=2,N
LENOFF(I)=LENR(I)
10 IW1(I,1)=IW1(I-1,1)+LENR(I-1)
c idisp(1) points to the first position in a/icn after the
c off-diagonal blocks and untreated rows.
20 IDISP(1)=IW1(N,1)+LENR(N)
c
c find row permutation ip to make diagonal zero-free.
CALL MC21A(N,ICN,LICN,IW1,LENR,IP,NUMNZ,IW)
c
c possible error return for structurally singular matrices.
IF (NUMNZ.NE.N.AND.ABORT) GO TO 170
c
c iw1(.,2) and lenr are permutations of iw1(.,1) and lenr/lenoff
c suitable for entry
c to mc13d since matrix with these row pointer and length arrays
c has maximum number of non-zeros on the diagonal.
DO 30 II=1,N
I=IP(II)
IW1(II,2)=IW1(I,1)
30 LENR(II)=LENOFF(I)
c
c find symmetric permutation iq to block lower triangular form.
CALL MC13D(N,ICN,LICN,IW1(1,2),LENR,IQ,IW(1,4),NUM,IW)
c
IF (NUM.NE.1) GO TO 60
c
c action taken if matrix is irreducible.
c whole matrix is just moved to the end of the storage.
DO 40 I=1,N
LENR(I)=LENOFF(I)
IP(I)=I
40 IQ(I)=I
LENOFF(1)=-1
c idisp(1) is the first position after the last element in the
c off-diagonal blocks and untreated rows.
NZ=IDISP(1)-1
IDISP(1)=1
c idisp(2) is the position in a/icn of the first element in the
c diagonal blocks.
IDISP(2)=LICN-NZ+1
LARGE=N
IF (NZ.EQ.LICN) GO TO 230
DO 50 K=1,NZ
J=NZ-K+1
JJ=LICN-K+1
A(JJ)=A(J)
50 ICN(JJ)=ICN(J)
c 230 = return
GO TO 230
c
c data structure reordered.
c
c form composite row permutation ... ip(i) = ip(iq(i)).
60 DO 70 II=1,N
I=IQ(II)
70 IW(II,1)=IP(I)
DO 80 I=1,N
80 IP(I)=IW(I,1)
c
c run through blocks in reverse order separating diagonal blocks
c which are moved to the end of the storage. elements in
c off-diagonal blocks are left in place unless a compress is
c necessary.
c
c ibeg indicates the lowest value of j for which icn(j) has been
c set to zero when element in position j was moved to the
c diagonal block part of storage.
IBEG=LICN+1
c iend is the position of the first element of those treated rows
c which are in diagonal blocks.
IEND=LICN+1
c large is the dimension of the largest block encountered so far.
LARGE=0
c
c num is the number of diagonal blocks.
DO 150 K=1,NUM
IBLOCK=NUM-K+1
c i1 is first row (in permuted form) of block iblock.
c i2 is last row (in permuted form) of block iblock.
I1=IW(IBLOCK,4)
I2=N
IF (K.NE.1) I2=IW(IBLOCK+1,4)-1
LARGE=MAX0(LARGE,I2-I1+1)
c go through the rows of block iblock in the reverse order.
DO 140 II=I1,I2
INEW=I2-II+I1
c we now deal with row inew in permuted form (row iold in original
c matrix).
IOLD=IP(INEW)
c if there is space to move up diagonal block portion of row go to 110
IF (IEND-IDISP(1).GE.LENOFF(IOLD)) GO TO 110
c
c in-line compress.
c moves separated off-diagonal elements and untreated rows to
c front of storage.
JNPOS=IBEG
ILEND=IDISP(1)-1
IF (ILEND.LT.IBEG) GO TO 190
DO 90 J=IBEG,ILEND
IF (ICN(J).EQ.0) GO TO 90
ICN(JNPOS)=ICN(J)
A(JNPOS)=A(J)
JNPOS=JNPOS+1
90 CONTINUE
IDISP(1)=JNPOS
IF (IEND-JNPOS.LT.LENOFF(IOLD)) GO TO 190
IBEG=LICN+1
c reset pointers to the beginning of the rows.
DO 100 I=2,N
100 IW1(I,1)=IW1(I-1,1)+LENOFF(I-1)
c
c row iold is now split into diag. and off-diag. parts.
110 IROWB=IW1(IOLD,1)
LENI=0
IROWE=IROWB+LENOFF(IOLD)-1
c backward scan of whole of row iold (in original matrix).
IF (IROWE.LT.IROWB) GO TO 130
DO 120 JJ=IROWB,IROWE
J=IROWE-JJ+IROWB
JOLD=ICN(J)
c iw(.,2) holds the inverse permutation to iq.
c ..... it was set to this in mc13d.
JNEW=IW(JOLD,2)
c if (jnew.lt.i1) then ....
c element is in off-diagonal block and so is left in situ.
IF (JNEW.LT.I1) GO TO 120
c element is in diagonal block and is moved to the end of the storage.
IEND=IEND-1
A(IEND)=A(J)
ICN(IEND)=JNEW
IBEG=MIN0(IBEG,J)
ICN(J)=0
LENI=LENI+1
120 CONTINUE
c
LENOFF(IOLD)=LENOFF(IOLD)-LENI
130 LENR(INEW)=LENI
140 CONTINUE
c
IP(I2)=-IP(I2)
150 CONTINUE
c resets ip(n) to positive value.
IP(N)=-IP(N)
c idisp(2) is position of first element in diagonal blocks.
IDISP(2)=IEND
c
c this compress is used to move all off-diagonal elements to the
c front of the storage.
IF (IBEG.GT.LICN) GO TO 230
JNPOS=IBEG
ILEND=IDISP(1)-1
DO 160 J=IBEG,ILEND
IF (ICN(J).EQ.0) GO TO 160
ICN(JNPOS)=ICN(J)
A(JNPOS)=A(J)
JNPOS=JNPOS+1
160 CONTINUE
c idisp(1) is first position after last element of off-diagonal blocks.
IDISP(1)=JNPOS
GO TO 230
c
c
c error return
170 IF (LP.NE.0) WRITE(LP,180) NUMNZ
180 FORMAT(33X,41H MATRIX IS STRUCTURALLY SINGULAR, RANK = ,I6)
IDISP(1)=-1
GO TO 210
190 IF (LP.NE.0) WRITE(LP,200) N
200 FORMAT(33X,33H LICN NOT BIG ENOUGH INCREASE BY ,I6)
IDISP(1)=-2
210 IF (LP.NE.0) WRITE(LP,220)
220 FORMAT(33H+ERROR RETURN FROM MC23AD BECAUSE)
c
230 RETURN
END
c######date 01 jan 1984 copyright ukaea, harwell.
c######alias mc24ad
SUBROUTINE MC24AD(N,ICN,A,LICN,LENR,LENRL,W)
DOUBLE PRECISION A(LICN),W(N),AMAXL,WROWL,AMAXU,ZERO
INTEGER ICN(LICN),LENR(N),LENRL(N)
DATA ZERO/0.0D0/
AMAXL=ZERO
DO 10 I=1,N
10 W(I)=ZERO
J0=1
DO 100 I=1,N
IF (LENR(I).EQ.0) GO TO 100
J2=J0+LENR(I)-1
IF (LENRL(I).EQ.0) GO TO 50
c calculation of 1-norm of l.
J1=J0+LENRL(I)-1
WROWL=ZERO
DO 30 JJ=J0,J1
30 WROWL=WROWL+DABS(A(JJ))
c amaxl is the maximum norm of columns of l so far found.
AMAXL=DMAX1(AMAXL,WROWL)
J0=J1+1
c calculation of norms of columns of u (max-norms).
50 J0=J0+1
IF (J0.GT.J2) GO TO 90
DO 80 JJ=J0,J2
J=ICN(JJ)
80 W(J)=DMAX1(DABS(A(JJ)),W(J))
90 J0=J2+1
100 CONTINUE
c amaxu is set to maximum max-norm of columns of u.
AMAXU=ZERO
DO 200 I=1,N
200 AMAXU=DMAX1(AMAXU,W(I))
c grofac is max u max-norm times max l 1-norm.
W(1)=AMAXL*AMAXU
RETURN
END
BLOCK DATA MABLD1
c
c comments on all the common block variables are given here even
c though some are not initialized by block data.
c lp,mp are used by the subroutine as the unit numbers for its warning
c and diagnostic messages. default value for both is 6 (for line
c printer output). the user can either reset them to a different
c stream number or suppress the output by setting them to zero.
c while lp directs the output of error diagnostics from the
c principal subroutines and internally called subroutines, mp
c controls only the output of a message which warns the user that he
c has input two or more non-zeros a(i), . . ,a(k) with the same row
c and column indices. the action taken in this case is to proceed
c using a numerical value of a(i)+...+a(k). in the absence of other
c errors, iflag will equal -14 on exit.
c lblock is a logical variable which controls an option of first
c preordering the matrix to block lower triangular form (using
c harwell subroutine mc23a). the preordering is performed if lblock
c is equal to its default value of .true. if lblock is set to
c .false. , the option is not invoked and the space allocated to
c ikeep can be reduced to 4*n+1.
c grow is a logical variable. if it is left at its default value of
c .true. , then on return from ma28a/ad or ma28b/bd, w(1) will give
c an estimate (an upper bound) of the increase in size of elements
c encountered during the decomposition. if the matrix is well
c scaled, then a high value for w(1), relative to the largest entry
c in the input matrix, indicates that the lu decomposition may be
c inaccurate and the user should be wary of his results and perhaps
c increase u for subsequent runs. we would like to emphasise that
c this value only relates to the accuracy of our lu decomposition
c and gives no indication as to the singularity of the matrix or the
c accuracy of the solution. this upper bound can be a significant
c overestimate particularly if the matrix is badly scaled. if an
c accurate value for the growth is required, lbig (q.v.) should be
c set to .true.
c eps,rmin are real variables. if, on entry to ma28b/bd, eps is less
c than one, then rmin will give the smallest ratio of the pivot to
c the largest element in the corresponding row of the upper
c triangular factor thus monitoring the stability of successive
c factorizations. if rmin becomes very large and w(1) from
c ma28b/bd is also very large, it may be advisable to perform a
c new decomposition using ma28a/ad.
c resid is a real variable which on exit from ma28c/cd gives the value
c of the maximum residual over all the equations unsatisfied because
c of dependency (zero pivots).
c irncp,icncp are integer variables which monitor the adequacy of "elbow
c room" in irn and a/icn respectively. if either is quite large (say
c greater than n/10), it will probably pay to increase the size of
c the corresponding array for subsequent runs. if either is very low
c or zero then one can perhaps save storage by reducing the size of
c the corresponding array.
c minirn,minicn are integer variables which, in the event of a
c successful return (iflag ge 0 or iflag=-14) give the minimum size
c of irn and a/icn respectively which would enable a successful run
c on an identical matrix. on an exit with iflag equal to -5, minicn
c gives the minimum value of icn for success on subsequent runs on
c an identical matrix. in the event of failure with iflag= -6, -4,
c -3, -2, or -1, then minicn and minirn give the minimum value of
c licn and lirn respectively which would be required for a
c successful decomposition up to the point at which the failure
c occurred.
c irank is an integer variable which gives an upper bound on the rank of
c the matrix.
c abort1 is a logical variable with default value .true. if abort1 is
c set to .false. then ma28a/ad will decompose structurally singular
c matrices (including rectangular ones).
c abort2 is a logical variable with default value .true. if abort2 is
c set to .false. then ma28a/ad will decompose numerically singular
c matrices.
c idisp is an integer array of length 2. on output from ma28a/ad, the
c indices of the diagonal blocks of the factors lie in positions
c idisp(1) to idisp(2) of a/icn. this array must be preserved
c between a call to ma28a/ad and subsequent calls to ma28b/bd,
c ma28c/cd or ma28i/id.
c tol is a real variable. if it is set to a positive value, then any
c non-zero whose modulus is less than tol will be dropped from the
c factorization. the factorization will then require less storage
c but will be inaccurate. after a run of ma28a/ad with tol positive
c it is not possible to use ma28b/bd and the user is recommended to
c use ma28i/id to obtain the solution. the default value for tol is
c 0.0.
c themax is a real variable. on exit from ma28a/ad, it will hold the
c largest entry of the original matrix.
c big is a real variable. if lbig has been set to .true., big will hold
c the largest entry encountered during the factorization by ma28a/ad
c or ma28b/bd.
c dxmax is a real variable. on exit from ma28i/id, dxmax will be set to
c the largest component of the solution.
c errmax is a real variable. on exit from ma28i/id, if maxit is
c positive, errmax will be set to the largest component in the
c estimate of the error.
c dres is a real variable. on exit from ma28i/id, if maxit is positive,
c dres will be set to the largest component of the residual.
c cgce is a real variable. it is used by ma28i/id to check the
c convergence rate. if the ratio of successive corrections is
c not less than cgce then we terminate since the convergence
c rate is adjudged too slow.
c ndrop is an integer variable. if tol has been set positive, on exit
c from ma28a/ad, ndrop will hold the number of entries dropped from
c the data structure.
c maxit is an integer variable. it is the maximum number of iterations
c performed by ma28i/id. it has a default value of 16.
c noiter is an integer variable. it is set by ma28i/id to the number of
c iterative refinement iterations actually used.
c nsrch is an integer variable. if nsrch is set to a value less than n,
c then a different pivot option will be employed by ma28a/ad. this
c may result in different fill-in and execution time for ma28a/ad.
c if nsrch is less than or equal to n, the workspace array iw can be
c reduced in length. the default value for nsrch is 32768.
c istart is an integer variable. if istart is set to a value other than
c zero, then the user must supply an estimate of the solution to
c ma28i/id. the default value for istart is zero.
c lbig is a logical variable. if lbig is set to .true., the value of the
c largest element encountered in the factorization by ma28a/ad or
c ma28b/bd is returned in big. setting lbig to .true. will
c increase the time for ma28a/ad marginally and that for ma28b/bd
c by about 20%. the default value for lbig is .false.
c
DOUBLE PRECISION EPS, RMIN, RESID, TOL, THEMAX, BIG, DXMAX,
* ERRMAX, DRES, CGCE
LOGICAL LBLOCK, GROW, ABORT1, ABORT2, LBIG
COMMON /MA28ED/ LP, MP, LBLOCK, GROW
COMMON /MA28FD/ EPS, RMIN, RESID, IRNCP, ICNCP, MINIRN, MINICN,
* IRANK, ABORT1, ABORT2
c common /ma28gd/ idisp(2)
COMMON /MA28HD/ TOL, THEMAX, BIG, DXMAX, ERRMAX, DRES, CGCE,
* NDROP, MAXIT, NOITER, NSRCH, ISTART, LBIG
DATA EPS /1.0D-4/, TOL /0.0D0/, CGCE /0.5D0/
DATA MAXIT /16/
DATA LP /6/, MP /6/, NSRCH /32768/, ISTART /0/
DATA LBLOCK /.TRUE./, GROW /.TRUE./, LBIG /.FALSE./
DATA ABORT1 /.TRUE./, ABORT2 /.TRUE./
END
BLOCK DATA MABLD2
c although all common block variables do not have default values,
c we comment on all the common block variables here.
c
c common block ma30e/ed holds control parameters ....
c common /ma30ed/ lp, abort1, abort2, abort3
c the integer lp is the unit number to which the error messages are
c sent. lp has a default value of 6. this default value can be
c reset by the user, if desired. a value of 0 suppresses all
c messages.
c the logical variables abort1,abort2,abort3 are used to control the
c conditions under which the subroutine will terminate.
c if abort1 is .true. then the subroutine will exit immediately on
c detecting structural singularity.
c if abort2 is .true. then the subroutine will exit immediately on
c detecting numerical singularity.
c if abort3 is .true. then the subroutine will exit immediately when
c the available space in a/icn is filled up by the previously
c decomposed, active, and undecomposed parts of the matrix.
c the default values for abort1,abort2,abort3 are set to .true.,.true.
c and .false. respectively.
c
c the variables in the common block ma30f/fd are used to provide the
c user with information on the decomposition.
c common /ma30fd/ irncp, icncp, irank, minirn, minicn
c irncp and icncp are integer variables used to monitor the adequacy
c of the allocated space in arrays irn and a/icn respectively, by
c taking account of the number of data management compresses
c required on these arrays. if irncp or icncp is fairly large (say
c greater than n/10), it may be advantageous to increase the size
c of the corresponding array(s). irncp and icncp are initialized
c to zero on entry to ma30a/ad and are incremented each time the
c compressing routine ma30d/dd is entered.
c icncp is the number of compresses on a/icn.
c irncp is the number of compresses on irn.
c irank is an integer variable which gives an estimate (actually an
c upper bound) of the rank of the matrix. on an exit with iflag
c equal to 0, this will be equal to n.
c minirn is an integer variable which, after a successful call to
c ma30a/ad, indicates the minimum length to which irn can be
c reduced while still permitting a successful decomposition of the
c same matrix. if, however, the user were to decrease the length
c of irn to that size, the number of compresses (irncp) may be
c very high and quite costly. if lirn is not large enough to begin
c the decomposition on a diagonal block, minirn will be equal to
c the value required to continue the decomposition and iflag will
c be set to -3 or -6. a value of lirn slightly greater than this
c (say about n/2) will usually provide enough space to complete
c the decomposition on that block. in the event of any other
c failure minirn gives the minimum size of irn required for a
c successful decomposition up to that point.
c minicn is an integer variable which after a successful call to
c ma30a/ad, indicates the minimum size of licn required to enable
c a successful decomposition. in the event of failure with iflag=
c -5, minicn will, if abort3 is left set to .false., indicate the
c minimum length that would be sufficient to prevent this error in
c a subsequent run on an identical matrix. again the user may
c prefer to use a value of icn slightly greater than minicn for
c subsequent runs to avoid too many conpresses (icncp). in the
c event of failure with iflag equal to any negative value except
c -4, minicn will give the minimum length to which licn could be
c reduced to enable a successful decomposition to the point at
c which failure occurred. notice that, on a successful entry
c idisp(2) gives the amount of space in a/icn required for the
c decomposition while minicn will usually be slightly greater
c because of the need for "elbow room". if the user is very
c unsure how large to make licn, the variable minicn can be used
c to provide that information. a preliminary run should be
c performed with abort3 left set to .false. and licn about 3/2
c times as big as the number of non-zeros in the original matrix.
c unless the initial problem is very sparse (when the run will be
c successful) or fills in extremely badly (giving an error return
c with iflag equal to -4), an error return with iflag equal to -5
c should result and minicn will give the amount of space required
c for a successful decomposition.
c
c common block ma30g/gd is used by the ma30b/bd entry only.
c common /ma30gd/ eps, rmin
c eps is a real/double precision variable. it is used to test for
c small pivots. its default value is 1.0e-4 (1.0d-4 in d version).
c if the user sets eps to any value greater than 1.0, then no
c check is made on the size of the pivots. although the absence of
c such a check would fail to warn the user of bad instability, its
c absence will enable ma30b/bd to run slightly faster. an a
c posteriori check on the stability of the factorization can be
c obtained from mc24a/ad.
c rmin is a real/double precision variable which gives the user some
c information about the stability of the decomposition. at each
c stage of the lu decomposition the magnitude of the pivot apiv
c is compared with the largest off-diagonal entry currently in its
c row (row of u), rowmax say. if the ratio
c min (apiv/rowmax)
c where the minimum is taken over all the rows, is less than eps
c then rmin is set to this minimum value and iflag is returned
c with the value +i where i is the row in which this minimum
c occurs. if the user sets eps greater than one, then this test
c is not performed. in this case, and when there are no small
c pivots rmin will be set equal to eps.
c
c common block ma30h/hd is used by ma30c/cd only.
c common /ma30hd/ resid
c resid is a real/double precision variable. in the case of singular
c or rectangular matrices its final value will be equal to the
c maximum residual for the unsatisfied equations; otherwise its
c value will be set to zero.
c
c common block ma30i/id controls the use of drop tolerances, the
c modified pivot option and the the calculation of the largest
c entry in the factorization process. this common block was added
c to the ma30 package in february, 1983.
c common /ma30id/ tol, big, ndrop, nsrch, lbig
c tol is a real/double precision variable. if it is set to a positive
c value, then ma30a/ad will drop from the factors any non-zero
c whose modulus is less than tol. the factorization will then
c require less storage but will be inaccurate. after a run of
c ma30a/ad where entries have been dropped, ma30b/bd should not
c be called. the default value for tol is 0.0.
c big is a real/double precision variable. if lbig has been set to
c .true., big will be set to the largest entry encountered during
c the factorization.
c ndrop is an integer variable. if tol has been set positive, on exit
c from ma30a/ad, ndrop will hold the number of entries dropped
c from the data structure.
c nsrch is an integer variable. if nsrch is set to a value less than
c or equal to n, then a different pivot option will be employed by
c ma30a/ad. this may result in different fill-in and execution
c time for ma30a/ad. if nsrch is less than or equal to n, the
c workspace arrays lastc and nextc are not referenced by ma30a/ad.
c the default value for nsrch is 32768.
c lbig is a logical variable. if lbig is set to .true., the value of
c the largest entry encountered in the factorization by ma30a/ad
c is returned in big. setting lbig to .true. will marginally
c increase the factorization time for ma30a/ad and will increase
c that for ma30b/bd by about 20%. the default value for lbig is
c .false.
c
DOUBLE PRECISION EPS, RMIN, TOL, BIG
LOGICAL ABORT1, ABORT2, ABORT3, LBIG
COMMON /MA30ED/ LP, ABORT1, ABORT2, ABORT3
COMMON /MA30GD/ EPS, RMIN
COMMON /MA30ID/ TOL, BIG, NDROP, NSRCH, LBIG
DATA EPS /1.0D-4/, TOL /0.0D0/, BIG /0.0D0/
DATA LP /6/, NSRCH /32768/
DATA LBIG /.FALSE./
DATA ABORT1 /.TRUE./, ABORT2 /.TRUE./, ABORT3 /.FALSE./
END
BLOCK DATA MABLD3
LOGICAL ABORT
COMMON /MC23BD/ LP,NUMNZ,NUM,LARGE,ABORT
DATA LP/6/,ABORT/.FALSE./
END