c######date 01 jan 1984 copyright ukaea, harwell. c######alias ma30ad 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