! a(ldiagU + nrowu) = abest ! This was in pivot order. !!! DEBUG
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer(kind=ip), | intent(in) | :: | m | |||
| integer(kind=ip), | intent(in) | :: | n | |||
| integer(kind=ip), | intent(in) | :: | nelem | |||
| integer(kind=ip), | intent(in) | :: | lena | |||
| integer(kind=ip), | intent(inout) | :: | luparm(30) | |||
| real(kind=rp), | intent(inout) | :: | parmlu(30) | |||
| real(kind=rp), | intent(inout) | :: | a(lena) | |||
| integer(kind=ip), | intent(inout) | :: | indc(lena) | |||
| integer(kind=ip), | intent(inout) | :: | indr(lena) | |||
| integer(kind=ip), | intent(inout) | :: | p(m) | |||
| integer(kind=ip), | intent(inout) | :: | q(n) | |||
| integer(kind=ip), | intent(inout) | :: | lenc(n) | |||
| integer(kind=ip), | intent(inout) | :: | lenr(m) | |||
| integer(kind=ip), | intent(inout) | :: | locc(n) | |||
| integer(kind=ip), | intent(inout) | :: | locr(m) | |||
| integer(kind=ip), | intent(inout) | :: | iploc(n) | |||
| integer(kind=ip), | intent(inout) | :: | iqloc(m) | |||
| integer(kind=ip), | intent(inout) | :: | ipinv(m) | |||
| integer(kind=ip), | intent(inout) | :: | iqinv(n) | |||
| real(kind=rp), | intent(inout) | :: | w(n) | |||
| integer(kind=ip), | intent(in) | :: | lenH | |||
| real(kind=rp), | intent(inout) | :: | Ha(lenH) | |||
| integer(kind=ip), | intent(inout) | :: | Hj(lenH) | |||
| integer(kind=ip), | intent(inout) | :: | Hk(lenH) | |||
| real(kind=rp), | intent(inout) | :: | Amaxr(m) | |||
| integer(kind=ip), | intent(out) | :: | inform | |||
| integer(kind=ip), | intent(out) | :: | lenL | |||
| integer(kind=ip), | intent(out) | :: | lenU | |||
| integer(kind=ip), | intent(out) | :: | minlen | |||
| integer(kind=ip), | intent(out) | :: | mersum | |||
| integer(kind=ip), | intent(out) | :: | nUtri | |||
| integer(kind=ip), | intent(out) | :: | nLtri | |||
| integer(kind=ip), | intent(out) | :: | ndens1 | |||
| integer(kind=ip), | intent(out) | :: | ndens2 | |||
| integer(kind=ip), | intent(out) | :: | nrank | |||
| integer(kind=ip), | intent(in) | :: | nslack | |||
| real(kind=rp), | intent(out) | :: | Lmax | |||
| real(kind=rp), | intent(out) | :: | Umax | |||
| real(kind=rp), | intent(out) | :: | DUmax | |||
| real(kind=rp), | intent(out) | :: | DUmin | |||
| real(kind=rp), | intent(out) | :: | Akmax |
subroutine lu1fad( m , n , nelem , lena , luparm, parmlu, & a , indc , indr , p , q , & lenc , lenr , locc , locr , & iploc , iqloc , ipinv , iqinv , w , & lenH , Ha , Hj , Hk , Amaxr , & inform, lenL , lenU , minlen, mersum, & nUtri , nLtri , ndens1, ndens2, nrank , nslack, & Lmax , Umax , DUmax , DUmin , Akmax ) integer(ip), intent(in) :: m, n, nelem, lena, lenH, nslack integer(ip), intent(inout) :: luparm(30) real(rp), intent(inout) :: parmlu(30), a(lena), Amaxr(m), & w(n), Ha(lenH) integer(ip), intent(inout) :: indc(lena), indr(lena), & p(m) , q(n) , & lenc(n) , lenr(m) , & locc(n) , locr(m) , & iploc(n) , iqloc(m) , & ipinv(m), iqinv(n), & Hj(lenH) , Hk(lenH) integer(ip), intent(out) :: inform, lenL , lenU , & minlen, mersum, & nUtri , nLtri , ndens1, ndens2, nrank real(rp), intent(out) :: Lmax, Umax, DUmax, DUmin, Akmax !------------------------------------------------------------------ ! lu1fad is a driver for the numerical phase of lu1fac. ! At each stage it computes a column of L and a row of U, ! using a Markowitz criterion to select the pivot element, ! subject to a stability criterion that bounds the elements of L. ! ! 00 Jan 1986 Version documented in LUSOL paper: ! Gill, Murray, Saunders and Wright (1987), ! Maintaining LU factors of a general sparse matrix, ! Linear algebra and its applications 88/89, 239-270. ! ! 02 Feb 1989 Following Suhl and Aittoniemi (1987), the largest ! element in each column is now kept at the start of ! the column, i.e. in position locc(j) of a and indc. ! This should speed up the Markowitz searches. ! To save time on highly triangular matrices, we wait ! until there are no further columns of length 1 ! before setting and maintaining that property. ! ! 12 Apr 1989 ipinv and iqinv added (inverses of p and q) ! to save searching p and q for rows and columns ! altered in each elimination step. (Used in lu1pq2) ! ! 19 Apr 1989 Code segmented to reduce its size. ! lu1gau does most of the Gaussian elimination work. ! lu1mar does just the Markowitz search. ! lu1mxc moves biggest elements to top of columns. ! lu1pen deals with pending fill-in in the row list. ! lu1pq2 updates the row and column permutations. ! ! 26 Apr 1989 maxtie replaced by maxcol, maxrow in the Markowitz ! search. maxcol, maxrow change as density increases. ! ! 25 Oct 1993 keepLU implemented. ! ! 07 Feb 1994 Exit main loop early to finish off with a dense LU. ! densLU tells lu1fad whether to do it. ! 21 Dec 1994 Bug fixed. nrank was wrong after the call to lu1ful. ! 12 Nov 1999 A parallel version of dcopy gave trouble in lu1ful ! during left-shift of dense matrix D within a(*). ! Fixed this unexpected problem here in lu1fad ! by making sure the first and second D don't overlap. ! ! 13 Sep 2000 TCP (Threshold Complete Pivoting) implemented. ! lu2max added ! (finds aijmax from biggest elems in each col). ! Utri, Ltri and Spars1 phases apply. ! No switch to Dense CP yet. (Only TPP switches.) ! 14 Sep 2000 imax needed to remember row containing aijmax. ! 22 Sep 2000 For simplicity, lu1mxc always fixes ! all modified cols. ! (TPP spars2 used to fix just the first maxcol cols.) ! 08 Nov 2000: Speed up search for aijmax. ! Don't need to search all columns if the elimination ! didn't alter the col containing the current aijmax. ! 21 Nov 2000: lu1slk implemented for Utri phase with TCP ! to guard against deceptive triangular matrices. ! (Utri used to have aijtol >= 0.9999 to include ! slacks, but this allows other 1s to be accepted.) ! Utri now accepts slacks, but applies normal aijtol ! test to other pivots. ! 28 Nov 2000: TCP with empty cols must call lu1mxc and lu2max ! with ( lq1, n, ... ), not just ( 1, n, ... ). ! 23 Mar 2001: lu1fad bug with TCP. ! A col of length 1 might not be accepted as a pivot. ! Later it appears in a pivot row and temporarily ! has length 0 (when pivot row is removed ! but before the column is filled in). If it is the ! last column in storage, the preceding col also thinks ! it is "last". Trouble arises when the preceding col ! needs fill-in -- it overlaps the real "last" column. ! (Very rarely, same trouble might have happened if ! the drop tolerance caused columns to have length 0.) ! ! Introduced ilast to record the last row in row file, ! jlast to record the last col in col file. ! lu1rec returns ilast = indr(lrow + 1) ! or jlast = indc(lcol + 1). ! *** (Should be an output parameter, but didn't want to ! alter lu1rec's parameter list.) ! lu1rec also treats empty rows or cols safely. ! (Doesn't eliminate them!) ! *** 20 Dec 2015: Made ilast an output as it should be. ! ! 26 Apr 2002: Heap routines added for TCP. ! lu2max no longer needed. ! imax, jmax used only for printing. ! 01 May 2002: lu1DCP implemented (dense complete pivoting). ! Both TPP and TCP now switch to dense LU ! when density exceeds dens2. ! 06 May 2002: In dense mode, store diag(U) in natural order. ! 09 May 2002: lu1mCP implemented (Markowitz TCP via heap). ! 11 Jun 2002: lu1mRP implemented (Markowitz TRP). ! 28 Jun 2002: Fixed call to lu1mxr. ! 14 Dec 2002: lu1mSP implemented (Markowitz TSP). ! 15 Dec 2002: Both TPP and TSP can grab cols of length 1 ! during Utri. ! 19 Dec 2004: Hdelete(...) has new input argument Hlenin. ! 26 Mar 2006: lu1fad returns nrank = min( mrank, nrank ) ! and ignores nsing from lu1ful ! ! 10 Jan 2010: First f90 version. ! 03 Apr 2013: lu1mxr recoded to improve efficiency of TRP. ! 12 Dec 2015: nslack is now an input. ! 20 Dec 2015: lu1rec returns ilast as output parameter. !------------------------------------------------------------------ logical :: Utri, Ltri, spars1, spars2, dense, & densLU, keepLU, TCP, TPP, TRP, TSP real(rp) :: abest, aijmax, aijtol, amax, & dens1, dens2, diag, & Lij, Ltol, small, Uspace !------------------------------------------------------------------ ! Local variables !--------------- ! ! lcol is the length of the column file. It points to the last ! nonzero in the column list. ! lrow is the analogous quantity for the row file. ! lfile is the file length (lcol or lrow) after the most recent ! compression of the column list or row list. ! nrowd and ncold are the number of rows and columns in the ! matrix defined by the pivot column and row. They are the ! dimensions of the submatrix D being altered at this stage. ! melim and nelim are the number of rows and columns in the ! same matrix D, excluding the pivot column and row. ! mleft and nleft are the number of rows and columns ! still left to be factored. ! nzchng is the increase in nonzeros in the matrix that remains ! to be factored after the current elimination ! (usually negative). ! nzleft is the number of nonzeros still left to be factored. ! nspare is the space we leave at the end of the last row or ! column whenever a row or column is being moved to the end ! of its file. nspare = 1 or 2 might help reduce the ! number of file compressions when storage is tight. ! ! The row and column ordering permutes A into the form ! ! ------------------------ ! \ | ! \ U1 | ! \ | ! -------------------- ! |\ ! | \ ! | \ ! P A Q = | \ ! | \ ! | -------------- ! | | | ! | | | ! | L1 | A2 | ! | | | ! | | | ! -------------------- ! ! where the block A2 is factored as A2 = L2 U2. ! The phases of the factorization are as follows. ! ! Utri is true when U1 is being determined. ! Any column of length 1 is accepted immediately (if TPP). ! ! Ltri is true when L1 is being determined. ! lu1mar exits as soon as an acceptable pivot is found ! in a row of length 1. ! ! spars1 is true while the density of the (modified) A2 is less ! than the parameter dens1 = parmlu(7) = 0.3 say. ! lu1mar searches maxcol columns and maxrow rows, ! where maxcol = luparm(3), maxrow = maxcol - 1. ! lu1mxc is used to keep the biggest element at the top ! of all remaining columns. ! ! spars2 is true while the density of the modified A2 is less ! than the parameter dens2 = parmlu(8) = 0.6 say. ! lu1mar searches maxcol columns and no rows. ! lu1mxc could fix up only the first maxcol cols (with TPP). ! 22 Sep 2000: For simplicity, lu1mxc fixes all modified cols. ! ! dense is true once the density of A2 reaches dens2. ! lu1mar searches only 1 column (the shortest). ! lu1mxc could fix up only the first column (with TPP). ! 22 Sep 2000: For simplicity, lu1mxc fixes all modified cols. !------------------------------------------------------------------ integer(ip) :: Hlen, Hlenin, hops, h, & i, ibest, ilast, imax, & j, jbest, jlast, jmax, lPiv, & k, kbest, kk, kslack, & l, last, lc, lc1, lcol, & lD, ldiagU, lenD, leni, lenj, & lfile, lfirst, lfree, limit, & ll, ll1, lpivc, lpivc1, lpivc2, & lpivr, lpivr1, lpivr2, lprint, & lq, lq1, lq2, lr, lr1, & lrow, ls, lsave, lu, lu1, & mark, maxcol, maxmn, maxrow, mbest, & melim, minfre, minmn, mleft, & mrank, ncold, nelim, nfill, & nfree, nleft, nout, nrowd, nrowu, & nsing, nspare, nzchng, nzleft integer(ip) :: markc(n), markr(m) real(rp) :: v nout = luparm(1) lprint = luparm(2) maxcol = luparm(3) lPiv = luparm(6) keepLU = luparm(8) /= 0 TPP = lPiv == 0 ! Threshold Partial Pivoting (normal). TRP = lPiv == 1 ! Threshold Rook Pivoting TCP = lPiv == 2 ! Threshold Complete Pivoting. TSP = lPiv == 3 ! Threshold Symmetric Pivoting. densLU = .false. maxrow = maxcol - 1 ilast = m ! Assume row m is last in the row file. jlast = n ! Assume col n is last in the col file. lfile = nelem lrow = nelem lcol = nelem minmn = min( m, n ) maxmn = max( m, n ) nzleft = nelem nspare = 1 ldiagU = 0 ! Keep -Wmaybe-uninitialized happy. if ( keepLU ) then lu1 = lena + 1 else ! Store only the diagonals of U in the top of memory. ldiagU = lena - n lu1 = ldiagU + 1 end if Ltol = parmlu(1) small = parmlu(3) Uspace = parmlu(6) dens1 = parmlu(7) dens2 = parmlu(8) Utri = .true. Ltri = .false. spars1 = .false. spars2 = .false. dense = .false. kslack = 0 ! 12 Dec 2015: Count slacks accepted during Utri. ! Check parameters. Ltol = max( Ltol, 1.0001_rp ) dens1 = min( dens1, dens2 ) ! Initialize output parameters. ! lenL, lenU, minlen, mersum, nUtri, nLtri, ndens1, ndens2, nrank, ! nslack, are already initialized by lu1fac. lenL = 0 lenU = 0 minlen = 0 mersum = 0 nLtri = 0 nUtri = 0 ndens1 = 0 ndens2 = 0 nrank = 0 Lmax = zero Umax = zero DUmax = zero DUmin = 1.0e+20_rp if (nelem == 0) Dumin = zero Akmax = zero hops = 0 ! More initialization. if (TPP .or. TSP) then ! Don't worry yet about lu1mxc. aijmax = zero aijtol = zero Hlen = 1 else ! TRP or TCP ! Move biggest element to top of each column. ! Set w(*) to mark slack columns (unit vectors). ! 12 Dec 2015: lu1fac (lu1slk) sets w(*) before lu1fad. ! 13 Dec 2015: lu1mxc fixed (empty cols caused trouble). call lu1mxc( i1, n, q, a, indc, lenc, locc ) ! call lu1slk( m, n, lena, q, iqloc, a, locc, w ) end if if (TRP) then ! Find biggest element in each row. mark = 0 call lu1mxr( mark, i1, m, m, n, lena, inform, & a, indc, lenc, locc, indr, lenr, locr, & p, markc, markr, Amaxr ) if (inform > 0) go to 981 end if if (TCP) then ! Set Ha(1:Hlen) = biggest element in each column, ! Hj(1:Hlen) = corresponding column indices. ! 17 Dec 2015: Allow for empty columns. Hlen = 0 do kk = 1, n Hlen = Hlen + 1 j = q(kk) if (lenc(j) > 0) then lc = locc(j) amax = abs( a(lc) ) else amax = zero end if Ha(Hlen) = amax Hj(Hlen) = j Hk(j) = Hlen end do ! Build the heap, creating new Ha, Hj and setting Hk(1:Hlen). call Hbuild( Ha, Hj, Hk, Hlen, Hlen, hops ) end if !------------------------------------------------------------------ ! Start of main loop. !------------------------------------------------------------------ mleft = m + 1 nleft = n + 1 do 800 nrowu = 1, minmn ! mktime = (nrowu / ntime) + 4 ! eltime = (nrowu / ntime) + 9 mleft = mleft - 1 nleft = nleft - 1 ! Bail out if there are no nonzero rows left. if (iploc(1) > m) go to 900 ! For TCP, the largest Aij is at the top of the heap. if ( TCP ) then aijmax = Ha(1) ! Marvelously easy ! Akmax = max( Akmax, aijmax ) aijtol = aijmax / Ltol end if !=============================================================== ! Find a suitable pivot element. !=============================================================== if ( Utri ) then !------------------------------------------------------------ ! So far all columns have had length 1. ! We are still looking for the (backward) triangular part of A ! that forms the first rows and columns of U. ! 12 Dec 2015: Use nslack and kslack to choose slacks first. !------------------------------------------------------------ lq1 = iqloc(1) lq2 = n if (m > 1) lq2 = iqloc(2) - 1 if (kslack < nslack) then do lq = lq1, lq2 j = q(lq) if (w(j) > zero) then ! Accept a slack kslack = kslack + 1 jbest = j lc = locc(jbest) ibest = indc(lc) abest = a(lc) mbest = 0 go to 300 end if end do ! DEBUG ERROR ! write(*,*) 'slack not found' ! write(*,*) 'kslack, nslack =', kslack, nslack ! stop else if (kslack == nslack) then ! Maybe print msg if (lprint >= 50) then write(nout,*) 'Slacks ended. nslack =', nslack end if kslack = nslack + 1 ! So print happens once end if ! All slacks will be grabbed before we get here. if (lq1 <= lq2) then ! There are more cols of length 1. if (TPP .or. TSP) then jbest = q(lq1) ! Grab the first one. else ! TRP or TCP ! Scan all columns of length 1. jbest = 0 do lq = lq1, lq2 j = q(lq) ! 12 Dec 2015: Slacks grabbed earlier. ! if (w(j) > zero) then ! Accept a slack ! jbest = j ! go to 250 ! end if lc = locc(j) amax = abs( a(lc) ) if (TRP) then i = indc(lc) aijtol = Amaxr(i) / Ltol end if if (amax >= aijtol) then jbest = j go to 250 end if end do end if 250 if (jbest > 0) then lc = locc(jbest) ibest = indc(lc) mbest = 0 go to 300 end if end if ! This is the end of the U triangle. ! We will not return to this part of the code. ! TPP and TSP call lu1mxc for the first time ! (to move biggest element to top of each column). if (lprint >= 50) then write(nout, 1100) 'Utri ended. spars1 = true' end if Utri = .false. Ltri = .true. spars1 = .true. nUtri = nrowu - 1 if (TPP .or. TSP) then call lu1mxc( lq1, n, q, a, indc, lenc, locc ) end if end if if ( spars1 ) then !------------------------------------------------------------ ! Perform a Markowitz search. ! Search cols of length 1, then rows of length 1, ! then cols of length 2, then rows of length 2, etc. !------------------------------------------------------------ ! if (TPP) then ! 12 Jun 2002: Next line disables lu1mCP below if (TPP .or. TCP) then call lu1mar( m , n , lena , maxmn, & TCP , aijtol, Ltol , maxcol, maxrow, & ibest, jbest , mbest , & a , indc , indr , p , q, & lenc , lenr , locc , locr , & iploc, iqloc ) else if (TRP) then call lu1mRP( m , n , lena , maxmn, & Ltol , maxcol, maxrow, & ibest, jbest , mbest , & a , indc , indr , p , q, & lenc , lenr , locc , locr , & iploc, iqloc , Amaxr ) ! else if (TCP) then ! Disabled by test above ! call lu1mCP( m , n , lena , aijtol, & ! ibest, jbest , mbest , & ! a , indc , indr , & ! lenc , lenr , locc , & ! Hlen , Ha , Hj ) else if (TSP) then call lu1mSP( m , n , lena , maxmn, & Ltol , maxcol, & ibest, jbest , mbest , & a , indc , q , locc , iqloc ) if (ibest == 0) go to 990 end if if ( Ltri ) then ! So far all rows have had length 1. ! We are still looking for the (forward) triangle of A ! that forms the first rows and columns of L. if (mbest > 0) then Ltri = .false. nLtri = nrowu - 1 - nUtri if (lprint >= 50) then write(nout, 1100) 'Ltri ended.' end if end if else ! See if what's left is as dense as dens1. if (nzleft >= (dens1 * mleft) * nleft) then spars1 = .false. spars2 = .true. ndens1 = nleft maxrow = 0 if (lprint >= 50) then write(nout, 1100) 'spars1 ended. spars2 = true' end if end if end if else if ( spars2 .or. dense ) then !------------------------------------------------------------ ! Perform a restricted Markowitz search, ! looking at only the first maxcol columns. (maxrow = 0.) !------------------------------------------------------------ ! if (TPP) then ! 12 Jun 2002: Next line disables lu1mCP below if (TPP .or. TCP) then call lu1mar( m , n , lena , maxmn, & TCP , aijtol, Ltol , maxcol, maxrow, & ibest, jbest , mbest , & a , indc , indr , p , q, & lenc , lenr , locc , locr , & iploc, iqloc ) else if (TRP) then call lu1mRP( m , n , lena , maxmn, & Ltol , maxcol, maxrow, & ibest, jbest , mbest , & a , indc , indr , p , q, & lenc , lenr , locc , locr , & iploc, iqloc , Amaxr ) ! else if (TCP) then ! Disabled by test above ! call lu1mCP( m , n , lena , aijtol, & ! ibest, jbest , mbest , & ! a , indc , indr , & ! lenc , lenr , locc , & ! Hlen , Ha , Hj ) else if (TSP) then call lu1mSP( m , n , lena , maxmn, & Ltol , maxcol, & ibest, jbest , mbest , & a , indc , q , locc , iqloc ) if (ibest == 0) go to 985 end if ! See if what's left is as dense as dens2. if ( spars2 ) then if (nzleft >= (dens2 * mleft) * nleft) then spars2 = .false. dense = .true. ndens2 = nleft maxcol = 1 if (lprint >= 50) then write(nout, 1100) 'spars2 ended. dense = true' end if end if end if end if !--------------------------------------------------------------- ! See if we can finish quickly. !--------------------------------------------------------------- if ( dense ) then lenD = mleft * nleft nfree = lu1 - 1 ! 28 Sep 2015: Change 2 to 3 for safety. if (nfree >= 3 * lenD) then ! There is room to treat the remaining matrix as ! a dense matrix D. ! We may have to compress the column file first. ! 12 Nov 1999: D used to be put at the ! beginning of free storage (lD = lcol + 1). ! Now put it at the end (lD = lu1 - lenD) ! so the left-shift in lu1ful will not ! involve overlapping storage ! (fatal with parallel dcopy). densLU = .true. ndens2 = nleft lD = lu1 - lenD if (lcol >= lD) then call lu1rec( n, .true., luparm, lcol, jlast, & lena, a, indc, lenc, locc ) lfile = lcol end if go to 900 end if end if !=============================================================== ! The best aij has been found. ! The pivot row ibest and the pivot column jbest ! define a dense matrix D of size nrowd x ncold. !=============================================================== 300 ncold = lenr(ibest) nrowd = lenc(jbest) melim = nrowd - 1 nelim = ncold - 1 mersum = mersum + mbest lenL = lenL + melim lenU = lenU + ncold if (lprint >= 50) then if (nrowu == 1) then write(nout, 1100) 'lu1fad debug:' end if if ( TPP .or. TRP .or. TSP ) then write(nout, 1200) nrowu, ibest, jbest, nrowd, ncold else ! TCP jmax = Hj(1) imax = indc(locc(jmax)) write(nout, 1200) nrowu, ibest, jbest, nrowd, ncold, & imax , jmax , aijmax end if end if !=============================================================== ! Allocate storage for the next column of L and next row of U. ! Initially the top of a, indc, indr are used as follows: ! ! ncold melim ncold melim ! ! a |...........|...........|ujbest..ujn|li1......lim| ! ! indc |...........| lenr(i) | lenc(j) | markl(i) | ! ! indr |...........| iqloc(i) | jfill(j) | ifill(i) | ! ! ^ ^ ^ ^ ^ ! lfree lsave lu1 ll1 oldlu1 ! ! Later the correct indices are inserted: ! ! indc | | | |i1........im| ! ! indr | | |jbest....jn|ibest..ibest| ! !=============================================================== if ( keepLU ) then ! relax else ! Always point to the top spot. ! Only the current column of L and row of U will ! take up space, overwriting the previous ones. lu1 = ldiagU + 1 end if ll1 = lu1 - melim lu1 = ll1 - ncold lsave = lu1 - nrowd lfree = lsave - ncold ! Make sure the column file has room. ! Also force a compression if its length exceeds a certain limit. limit = int(Uspace*real(lfile)) + m + n + 1000 minfre = ncold + melim nfree = lfree - lcol if (nfree < minfre .or. lcol > limit) then call lu1rec( n, .true., luparm, lcol, jlast, & lena, a, indc, lenc, locc ) lfile = lcol nfree = lfree - lcol if (nfree < minfre) go to 970 end if ! Make sure the row file has room. minfre = melim + ncold nfree = lfree - lrow if (nfree < minfre .or. lrow > limit) then call lu1rec( m, .false., luparm, lrow, ilast, & lena, a, indr, lenr, locr ) lfile = lrow nfree = lfree - lrow if (nfree < minfre) go to 970 end if !=============================================================== ! Move the pivot element to the front of its row ! and to the top of its column. !=============================================================== lpivr = locr(ibest) lpivr1 = lpivr + 1 lpivr2 = lpivr + nelim do l = lpivr, lpivr2 if (indr(l) == jbest) exit end do indr(l) = indr(lpivr) indr(lpivr) = jbest lpivc = locc(jbest) lpivc1 = lpivc + 1 lpivc2 = lpivc + melim do l = lpivc, lpivc2 if (indc(l) == ibest) exit end do indc(l) = indc(lpivc) indc(lpivc) = ibest abest = a(l) a(l) = a(lpivc) a(lpivc) = abest if ( keepLU ) then ! relax else ! Store just the diagonal of U, in natural order. !!! a(ldiagU + nrowu) = abest ! This was in pivot order. a(ldiagU + jbest) = abest end if !============================================================== ! Delete pivot col from heap. ! Hk tells us where it is in the heap. !============================================================== if ( TCP ) then kbest = Hk(jbest) Hlenin = Hlen call Hdelete( Ha, Hj, Hk, Hlenin, Hlen, n, kbest, h ) hops = hops + h end if !=============================================================== ! Delete the pivot row from the column file ! and store it as the next row of U. ! Set indr(lu) = 0 to initialize jfill ptrs on columns of D, ! indc(lu) = lenj to save the original column lengths. !=============================================================== a(lu1) = abest indr(lu1) = jbest indc(lu1) = nrowd lu = lu1 diag = abs( abest ) Umax = max( Umax, diag ) DUmax = max( DUmax, diag ) DUmin = min( DUmin, diag ) do lr = lpivr1, lpivr2 lu = lu + 1 j = indr(lr) lenj = lenc(j) lenc(j) = lenj - 1 lc1 = locc(j) last = lc1 + lenc(j) do l = lc1, last if (indc(l) == ibest) exit end do a(lu) = a(l) indr(lu) = 0 indc(lu) = lenj Umax = max( Umax, abs( a(lu) ) ) a(l) = a(last) indc(l) = indc(last) indc(last) = 0 ! Free entry if (j == jlast) lcol = lcol - 1 end do !=============================================================== ! Delete the pivot column from the row file ! and store the nonzeros of the next column of L. ! Set indc(ll) = 0 to initialize markl(*) markers, ! indr(ll) = 0 to initialize ifill(*) row fill-in cntrs, ! indc(ls) = leni to save the original row lengths, ! indr(ls) = iqloc(i) to save parts of iqloc(*), ! iqloc(i) = lsave - ls to point to the nonzeros of L ! = -1, -2, -3, ... in mark(*). !=============================================================== indc(lsave) = ncold if (melim == 0) go to 700 ll = ll1 - 1 ls = lsave abest = one / abest do lc = lpivc1, lpivc2 ll = ll + 1 ls = ls + 1 i = indc(lc) leni = lenr(i) lenr(i) = leni - 1 lr1 = locr(i) last = lr1 + lenr(i) do l = lr1, last if (indr(l) == jbest) exit end do indr(l) = indr(last) indr(last) = 0 ! Free entry if (i == ilast) lrow = lrow - 1 a(ll) = - a(lc) * abest Lij = abs( a(ll) ) Lmax = max( Lmax, Lij ) !!!!! DEBUG ! if (Lij > Ltol) then ! write( * ,*) ' Big Lij!!!', nrowu ! write(nout,*) ' Big Lij!!!', nrowu ! end if indc(ll) = 0 indr(ll) = 0 indc(ls) = leni indr(ls) = iqloc(i) iqloc(i) = lsave - ls end do !=============================================================== ! Do the Gaussian elimination. ! This involves adding a multiple of the pivot column ! to all other columns in the pivot row. ! ! Sometimes more than one call to lu1gau is needed to allow ! compression of the column file. ! lfirst says which column the elimination should start with. ! minfre is a bound on the storage needed for any one column. ! lu points to off-diagonals of u. ! nfill keeps track of pending fill-in in the row file. !=============================================================== if (nelim == 0) go to 700 lfirst = lpivr1 minfre = mleft + nspare lu = 1 nfill = 0 400 call lu1gau( m , melim , ncold , nspare, small , & lpivc1, lpivc2, lfirst, lpivr2, lfree , minfre, & ilast , jlast , lrow , lcol , lu , nfill , & a , indc , indr , & lenc , lenr , locc , locr , & iqloc , a(ll1), indc(ll1), & a(lu1), indr(ll1), indr(lu1) ) if (lfirst > 0) then ! The elimination was interrupted. ! Compress the column file and try again. ! lfirst, lu and nfill have appropriate new values. call lu1rec( n, .true., luparm, lcol, jlast, & lena, a, indc, lenc, locc ) lfile = lcol lpivc = locc(jbest) lpivc1 = lpivc + 1 lpivc2 = lpivc + melim nfree = lfree - lcol if (nfree < minfre) go to 970 go to 400 end if !=============================================================== ! The column file has been fully updated. ! Deal with any pending fill-in in the row file. !=============================================================== if (nfill > 0) then ! Compress the row file if necessary. ! lu1gau has set nfill to be the number of pending fill-ins ! plus the current length of any rows that need to be moved. minfre = nfill nfree = lfree - lrow if (nfree < minfre) then call lu1rec( m, .false., luparm, lrow, ilast, & lena, a, indr, lenr, locr ) lfile = lrow lpivr = locr(ibest) lpivr1 = lpivr + 1 lpivr2 = lpivr + nelim nfree = lfree - lrow if (nfree < minfre) go to 970 end if ! Move rows that have pending fill-in to end of the row file. ! Then insert the fill-in. call lu1pen( m , melim , ncold , nspare, ilast, & lpivc1, lpivc2, lpivr1, lpivr2, lrow , & lenc , lenr , locc , locr , & indc , indr , indr(ll1), indr(lu1) ) end if !=============================================================== ! Restore the saved values of iqloc. ! Insert the correct indices for the col of L and the row of U. !=============================================================== 700 lenr(ibest) = 0 lenc(jbest) = 0 ll = ll1 - 1 ls = lsave do lc = lpivc1, lpivc2 ll = ll + 1 ls = ls + 1 i = indc(lc) iqloc(i) = indr(ls) indc(ll) = i indr(ll) = ibest end do lu = lu1 - 1 do lr = lpivr, lpivr2 lu = lu + 1 indr(lu) = indr(lr) end do !=============================================================== ! Free the space occupied by the pivot row ! and update the column permutation. ! Then free the space occupied by the pivot column ! and update the row permutation. ! ! nzchng is found in both calls to lu1pq2, but we use it only ! after the second. !=============================================================== call lu1pq2( ncold, nzchng, & indr(lpivr), indc( lu1 ), lenc, iqloc, q, iqinv ) call lu1pq2( nrowd, nzchng, & indc(lpivc), indc(lsave), lenr, iploc, p, ipinv ) nzleft = nzleft + nzchng !=============================================================== ! lu1mxr resets Amaxr(i) in each modified row i. ! lu1mxc moves the largest aij to the top of each modified col j. ! 28 Jun 2002: Note that cols of L have an implicit diag of 1.0, ! so lu1mxr is called with ll1, not ll1+1, whereas ! lu1mxc is called with lu1+1. !=============================================================== if (Utri .and. TPP) then ! Relax -- we're not keeping big elements at the top yet. else if (TRP .and. melim > 0) then ! Beware: The parts of p that we need are in indc(ll1:ll) ! 28 Sep 2015: inform is now an output. mark = mark + 1 call lu1mxr( mark, ll1, ll, m, n, lena, inform, & a, indc, lenc, locc, indr, lenr, locr, & indc, markc, markr, Amaxr ) ! ^^^^ Here are the p(k1:k2) needed by lu1mxr. if (inform > 0) go to 981 end if if (nelim > 0) then call lu1mxc( lu1+1, lu, indr, a, indc, lenc, locc ) if (TCP) then ! Update modified columns in heap ! 20 Dec 2015: Allow for empty columns. do kk = lu1+1, lu j = indr(kk) k = Hk(j) if (lenc(j) > 0) then v = abs( a(locc(j)) ) ! Biggest aij in column j else v = zero end if call Hchange( Ha, Hj, Hk, Hlen, n, k, v, j, h ) hops = hops + h end do end if end if end if !=============================================================== ! Negate lengths of pivot row and column so they will be ! eliminated during compressions. !=============================================================== lenr(ibest) = - ncold lenc(jbest) = - nrowd ! Test for fatal bug: row or column lists overwriting L and U. if (lrow > lsave) go to 980 if (lcol > lsave) go to 980 ! Reset the file lengths if pivot row or col was at the end. if (ibest == ilast) then lrow = locr(ibest) end if if (jbest == jlast) then lcol = locc(jbest) end if 800 end do !------------------------------------------------------------------ ! End of main loop. !------------------------------------------------------------------ !------------------------------------------------------------------ ! Normal exit. ! Move empty rows and cols to the end of p, q. ! Then finish with a dense LU if necessary. !------------------------------------------------------------------ 900 inform = 0 call lu1pq3( m, lenr, p, ipinv, mrank ) call lu1pq3( n, lenc, q, iqinv, nrank ) nrank = min( mrank, nrank ) if ( densLU ) then call lu1ful( m , n , lena , lenD , lu1 , TPP, & mleft , nleft, nrank, nrowu, & lenL , lenU , nsing, & keepLU, small, & a , a(lD), indc , indr , p , q, & lenc , lenr , locc , ipinv, locr ) !*** 21 Dec 1994: Bug in next line. !*** nrank = nrank - nsing. Changed to next line: !*** nrank = minmn - nsing !*** 26 Mar 2006: Previous line caused bug with m<n and nsing>0. ! Don't mess with nrank any more. Let end of lu1fac handle it. end if minlen = lenL + lenU + 2*(m + n) go to 990 ! Not enough space free after a compress. ! Set minlen to an estimate of the necessary value of lena. 970 inform = 7 minlen = lena + lfile + 2*(m + n) go to 990 ! Fatal error. This will never happen! ! (Famous last words.) 980 inform = 8 go to 990 ! Fatal error in lu1mxr. This will never happen! 981 inform = 10 go to 990 ! Fatal error with TSP. Diagonal pivot not found. 985 inform = 9 ! Exit. 990 return 1100 format(/ 1x, a) 1200 format(' nrowu', i7, ' i,jbest', 2i7, ' nrowd,ncold', 2i6, & ' i,jmax', 2i7, ' aijmax', es10.2) end subroutine lu1fad