*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ * LUMOD 2.0 (29 Jun 1999) * The following routines are members of a package for maintaining * L*C = U, an LU factorization of a dense square matrix C: * * LUmod Lprod Usolve * elm elmgen LUback LUforw * Also needed: * daxpy dcopy ddot (from the BLAS) * * The first three routines may be used as follows: * * call LUmod ( 1, ... ) to build up the matrix C, * call LUmod ( 1,2,3 or 4, ...) to modify the matrix C, * * call Lprod ( 1, ... b, x ) * call Usolve( 1, ... x ) to solve C x = b, * * call Usolve( 2, ... b ) * call Lprod ( 2, ... b, z ) to solve C'x = b. * * BEWARE: LUmod uses a machine-dependent constant eps. * * Michael Saunders * Systems Optimization Laboratory, * Department of EESOR, Stanford University. * * 05 Apr 1981: Original version used 2-D arrays and plane rotations. * 09 Oct 1986: Plane rotations optionally replaced by eliminations. * 17 Apr 1990: Mode 2 altered to put new column at end. * 26 Apr 1990: L and U stored by rows in 1-D arrays. * 21 Aug 1991: LUmod created from QRmod3. * 29 Jun 1999: Minor clean-up. * Comments added for each mode to show what updates * should be made to two lists, say listr(*) and listc(*), * that the user will probably need to maintain. * They contain the row and column indices of some bigger * matrix that are now the rows and columns of C. * 29 Jun 1999: Mode 2 changed back to original: put new column * in same place as the column being replaced. * This generally needs more arithmetic, but it * simplifies updating listr(*) and listc(*), * and it keeps C symmetric where possible. * 30 Jun 1999: Mode 4 now replaces a row and column by * swapping them with the last row and column. * ------------------------------------------------------------------ subroutine LUmod ( mode, maxmod, n, krow, kcol, L, U, y, z, w ) implicit double precision (a-h,o-z) integer mode, maxmod, n, krow, kcol double precision L(*), U(*), y(n), z(n), w(n) * ------------------------------------------------------------------ * LUmod modifies the matrix factorization L*C = U where * L and C are square, L is a product of 2x2 stabilized elimination * matrices, and U is upper triangular. * * The modifications depend upon 'mode' as follows. * n is not altered by any of the entries to LUmod. * * ================================================================== * mode=1. Expand L and U by adding a new row and column to C. * The dimension expands from n-1 to n (not from n to n+1). * ================================================================== * On entry: * y(*), z(*) contain the new row and column. * y(n) contains the new diagonal element of C. * On exit: * y(*) is altered. * Not used: * krow, kcol, w(*). * Updates (before call to LUmod): * n = n + 1 * listr(n) = new row index * listc(n) = new column index * * ================================================================== * mode=2. Replace the kcol-th column of C. * ================================================================== * On entry: * kcol says which column of C is being replaced. * z(*) contains the new column. * On exit: * y(*), w(*) are altered. * z(*) is not altered. * Not used: * krow * Updates: * listc(kcol) = new column index * * ================================================================== * mode=3. Replace the krow-th row of C. * ================================================================== * On entry: * krow says which row of C is being replaced. * y(*) contains the new row. * On exit: * y(*), z(*), w(*) are altered. * Not used: * kcol * Updates: * listr(krow) = new row index * * ================================================================== * mode=4. Shrink L and U by deleting the * krow-th row and the kcol-th column of C * (swapping them with the n-th row and column). * n should be decreased to n-1 by the user after the call. * ================================================================== * On entry: * krow, kcol say which row and col of C are being replaced. * On exit: * y(*), z(*), w(*) are altered. * Updates (after call to LUmod): * listr(krow) = listr(n) * listc(kcol) = listc(n) * n = n - 1 * * * ================================================================== * The LU factors are well defined even if C is singular * (in which case U is also singular). At some stage the factors * will be used to solve systems of linear equations of the form * C x = b or C' x = b. * The diagonals of U should then be reasonably different from zero. * ================================================================== * * * Other input parameters: * * maxmod The maximum dimension of C. * * L(*) An array of length at least maxmod*maxmod. * When L has maximum dimension maxmod, its rows are stored * contiguously in L(*). * For lower dimensions, each row of L starts in the same * place but fills only the front of its allowable space. * Row i of L starts in position (i - 1)*maxmod + 1. * * U(*) An array of length at least maxmod*(maxmod + 1)/2. * When U has maximum dimension maxmod, the upper-triangular * part of its rows are stored contiguously in U(*). * For lower dimensions, each row of U starts in the same * place but fills only the front of its allowable space. * Row i of U starts in position (i-1)*maxmod + (3-i)*i/2. * * w(n) A work vector. *----------------------------------------------------------------------- integer first double precision zero , one parameter ( zero = 0.0d+0, one = 1.0d+0 ) * MACHINE-DEPENDENT CONSTANT: * eps is the machine precision -- typically 2.22d-16. * A value that is slightly too large will be OK. double precision eps data eps / 2.22d-16 / n1 = n - 1 if (mode .eq. 1) then * --------------------------------------------------------------- * mode = 1. Add a row y and a column z. * The LU factors will expand in dimension from n-1 to n. * The new diagonal element of C is in y(n). * --------------------------------------------------------------- lastu = n1*maxmod + (3 - n)*n/2 lastl = n1*maxmod + n ls = n1*maxmod + 1 L(lastl) = one if (n .eq. 1) then U(lastu) = y(n) return end if * Compute L*z and temporarily store it * in the new bottom row of L. call Lprod ( 1, maxmod, n1, L, z, L(ls) ) * Copy L*z into the new last column of U. * Border L with zeros. ll = ls lu = n incu = maxmod - 1 do 110 j = 1, n1 U(lu) = L(ll) L(ll) = zero ll = ll + 1 lu = lu + incu incu = incu - 1 110 continue ll = n do 120 i = 1, n1 L(ll) = zero ll = ll + maxmod 120 continue * Add row y to the factorization * using a forward sweep of eliminations. last = n call LUforw( 1, last, n, n, maxmod, eps, L, U, y ) else if (mode .eq. 2) then * --------------------------------------------------------------- * mode=2. Replace the kcol-th column of C by the vector z. * --------------------------------------------------------------- * Compute w = L*z. call Lprod ( 1, maxmod, n, L, z, w ) * Copy the top of w into column kcol of U. lu = kcol incu = maxmod - 1 do i = 1, kcol U(lu) = w(i) lu = lu + incu incu = incu - 1 end do if (kcol .lt. n) then * Find w(last), the last nonzero in the bottom part of w. * Eliminate elements last-1, last-2, ... kcol+1 of w(*) * using a partial backward sweep of eliminations. first = kcol + 1 last = n call LUback( first, last, n, n, maxmod, eps, L, U, y, w ) y(kcol) = w(last) * Eliminate elements kcol, kcol+1, ... last-1 of y(*) * using a partial forward sweep of eliminations. call LUforw( kcol, last, n, n, maxmod, eps, L, U, y ) end if else if (mode .eq. 3) then * --------------------------------------------------------------- * mode=3. Replace the krow-th row of C by the vector y. * --------------------------------------------------------------- if (n .eq. 1) then L(1) = one U(1) = y(1) return end if * Copy the krow-th column of L into w, and zero the column. ll = krow do 330 i = 1, n w(i) = L(ll) L(ll) = zero ll = ll + maxmod 330 continue * Reduce the krow-th column of L to the unit vector e(last). * where 'last' is determined by LUback. * This is done by eliminating elements last-1, last-2, ..., 1 * using a backward sweep of eliminations. * On exit, row 'last' of U is a spike stored in z, whose first * nonzero entry is in z(first). However, z will be discarded. first = 1 last = n call LUback( first, last, n, n, maxmod, eps, L, U, z, w ) * Replace the 'last' row of L by the krow-th unit vector. ll = (last - 1)*maxmod do 340 j = 1, n L(ll + j) = zero 340 continue L(ll + krow) = one * Eliminate the elements of the new row y, * using a forward sweep of eliminations. call LUforw( 1, last, n, n, maxmod, eps, L, U, y ) else if (mode .eq. 4) then * --------------------------------------------------------------- * mode=4. Delete the krow-th row and the kcol-th column of C. * Replace them by the last row and column respectively. * --------------------------------------------------------------- * First, move the last column into position kcol. if (kcol .lt. n) then * Set w = last column of U. lu = n incu = maxmod - 1 do i = 1, n w(i) = U(lu) lu = lu + incu incu = incu - 1 end do * Copy the top of w into column kcol of U. lu = kcol incu = maxmod - 1 do i = 1, kcol U(lu) = w(i) lu = lu + incu incu = incu - 1 end do * U now has only n-1 columns. * Find w(last), the last nonzero in the bottom part of w. * Eliminate elements last-1, last-2, ... kcol+1 of w(*) * using a partial backward sweep of eliminations. first = kcol + 1 last = n call LUback( first, last, n, n1, maxmod, eps, L, U, y, w ) y(kcol) = w(last) * Eliminate elements kcol, kcol+1, ... last-1 of y(*) * using a partial forward sweep of eliminations. call LUforw( kcol, last, n, n1, maxmod, eps, L, U, y ) end if * Now, move the last row into position krow. * Swap columns krow and n of L, using w = krow-th column of L. call dcopy ( n, L(krow), maxmod, w, 1 ) if (krow .lt. n) then call dcopy ( n, L(n), maxmod, L(krow), maxmod ) end if * Reduce the last column of L (in w) to the unit vector e(n). * This is done by eliminating elements n-1, n-2, ..., 1 * using a backward sweep of eliminations. last = - n call LUback( 1, last, n, n1, maxmod, eps, L, U, z, w ) end if * End of LUmod end *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine Lprod ( mode, maxmod, n, L, y, z ) implicit double precision (a-h,o-z) double precision L(*), y(*), z(*) * ------------------------------------------------------------------ * If mode = 1, Lprod computes z = L*y. * If mode = 2, Lprod computes z = L(transpose)*y. * L is stored by rows in L(*). It is equivalent to storing * L(transpose) by columns in a 2-D array L(maxmod,n). * y is not altered. * ------------------------------------------------------------------ external daxpy, ddot ll = 1 if (mode .eq. 1) then do 100 i = 1, n z(i) = ddot ( n, L(ll), 1, y, 1 ) ll = ll + maxmod 100 continue else * call dzero ( n, z, 1 ) do 120 i = 1, n z(i) = 0.0d+0 120 continue do 200 i = 1, n call daxpy ( n, y(i), L(ll), 1, z, 1 ) ll = ll + maxmod 200 continue end if * End of Lprod end *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine Usolve( mode, maxmod, n, U, y ) implicit double precision (a-h,o-z) double precision U(*), y(n) * ------------------------------------------------------------------ * If mode = 1, Usolve solves U * y(new) = y(old). * If mode = 2, Usolve solves U(transpose) * y(new) = y(old). * U is upper triangular, stored by rows. * y is overwritten by the solution. * ------------------------------------------------------------------ if (mode .eq. 1) then lu = (n - 1)*maxmod + (3 - n)*n/2 y(n) = y(n) / U(lu) incu = maxmod + 1 - n nu = 0 do 100 i = n-1, 1, -1 nu = nu + 1 incu = incu + 1 lu = lu - incu sum = y(i) - ddot ( nu, U(lu+1), 1, y(i+1), 1 ) y(i) = sum / U(lu) 100 continue else lu = 1 incu = maxmod nu = n - 1 do 200 i = 1, n - 1 y(i) = y(i) / U(lu) call daxpy ( nu, (- y(i)), U(lu+1), 1, y(i+1), 1 ) lu = lu + incu incu = incu - 1 nu = nu - 1 200 continue y(n) = y(n) / U(lu) end if * End of Usolve end *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine elm ( n, x, y, cs, sn ) implicit double precision (a-h,o-z) double precision x(n), y(n) * ------------------------------------------------------------------ * elm computes the transformation (x y)*E and returns the result * in (x y), where the 2 by 2 matrix E is defined by cs and sn * as follows: * * E = ( 1 sn ) if cs .ge. zero, E = ( 1 ) otherwise. * ( 1 ) ( 1 sn ) * ------------------------------------------------------------------ parameter ( zero = 0.0d+0 ) if (cs .ge. zero) then if (sn .ne. zero) then *>>> call daxpy ( n, sn, x, 1, y, 1 ) do 10 i = 1, n y(i) = x(i)*sn + y(i) 10 continue end if else if (sn .eq. zero) then * Treat an interchange specially. *>>> call dswap ( n, x, 1, y, 1 ) do 20 i = 1, n xi = x(i) x(i) = y(i) y(i) = xi 20 continue else do 30 i = 1, n xi = x(i) x(i) = y(i) y(i) = xi + y(i)*sn 30 continue end if end if * End of elm end *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine elmgen( x, y, eps, cs, sn ) implicit double precision (a-h,o-z) * ------------------------------------------------------------------ * elmgen generates an elimination transformation E such that * (x y)*E = (x 0) or (y 0), depending on the relative * sizes of x and y. * eps is an input parameter -- the machine precision. * * CAUTION: It is assumed that the data generating x and y * are in general well-scaled, so that if both x and y * are smaller than 'tiny', they are both changed to zero. * This is an attempt to save work and reduce underflow. * ------------------------------------------------------------------ intrinsic abs parameter ( zero = 0.0d+0, one = 1.0d+0 ) tiny = eps * 1.0d-4 cs = zero if (abs(x) .ge. abs(y)) then if (abs(x) .le. tiny) then sn = zero x = zero else sn = - y/x end if else if (abs(y) .le. tiny) then sn = zero x = zero else cs = - one sn = - x/y x = y end if end if y = zero * End of elmgen end *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine LUback( first, last, n, nu, $ maxmod, eps, L, U, y, z ) implicit double precision (a-h,o-z) integer first, last, n, maxmod double precision L(*), U(*), y(n), z(n), eps * ------------------------------------------------------------------ * LUback updates the factors LC = U by performing a backward sweep * to eliminate all but the 'last' nonzero in the column vector z(*), * stopping at z(first). * * If 'last' is positive, LUback searches backwards for a nonzero * element in z and possibly alters 'last' accordingly. * Otherwise, 'last' will be reset to abs(last) and so used. * * L is n by n. * U is n by nu. * y(*) will eventually contain a row spike in row 'last' of U. * The 'spike' row of L begins at L(ls). * * 18 Mar 1990: First version with L and U stored row-wise. * 29 Jun 1999: Save w(last) = zlast at end so column replace * can do correct forward sweep. * ------------------------------------------------------------------ parameter ( zero = 0.0d+0 ) if (last .gt. 0) then * Find the last significant element in z(*). do 100 i = last, first, -1 if (abs(z(i)) .gt. eps) go to 120 100 continue i = first 120 last = i else last = abs(last) end if * Load the 'last' row of U into the end of y * and do the backward sweep. zlast = z(last) lu = (last - 1)*maxmod + (3 - last)*last/2 ll = (last - 1)*maxmod + 1 ls = ll incu = maxmod + 1 - last numu = nu + 1 - last if (numu .gt. 0) call dcopy ( numu, U(lu), 1, y(last), 1 ) do 700 lz = last - 1, first, -1 ly = lz ll = ll - maxmod incu = incu + 1 lu = lu - incu y(ly) = zero * See if this element of z is worth eliminating. * We compare z(lz) with the current last nonzero, zlast. if ( abs(z(lz)) .le. eps * abs(zlast) ) go to 700 * Generate a 2x2 elimination and apply it to U and L. numu = nu + 1 - ly call elmgen( zlast, z(lz), eps , cs, sn ) call elm ( numu , y(ly), U(lu), cs, sn ) call elm ( n , L(ls), L(ll), cs, sn ) 700 continue z(last) = zlast * end of LUback end *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine LUforw( first, last, n, nu, maxmod, eps, L, U, y ) implicit double precision (a-h,o-z) integer first, last, n, maxmod double precision L(*), U(*), y(n), eps * ------------------------------------------------------------------ * LUforw updates the factors LC = U by performing a forward sweep * to eliminate subdiagonals in a row spike in row 'last' of U. * * L is n by n. * U is n by nu. * y(*) contains the row spike. The first nonzero to be eliminated * is in y(first) or later. * The 'spike' row of L begins at L(ls). * * 18 Mar 1990: First version with L and U stored row-wise. * ------------------------------------------------------------------ lu = (first - 1)*maxmod + (3 - first)*first/2 ll = (first - 1)*maxmod + 1 ls = (last - 1)*maxmod + 1 incu = maxmod + 1 - first do 700 ly = first, last - 1 * See if this element of y is worth eliminating. * We compare y(ly) with the corresponding diagonal of U. if ( abs(y(ly)) .le. eps * abs(U(lu)) ) go to 680 * Generate a 2x2 elimination and apply it to U and L. numu = nu - ly call elmgen( U(lu), y(ly) , eps , cs, sn ) if (numu .gt. 0) $ call elm ( numu , U(lu+1), y(ly+1), cs, sn ) call elm ( n , L(ll ), L(ls ), cs, sn ) 680 ll = ll + maxmod lu = lu + incu incu = incu - 1 700 continue * Copy the remaining part of y into U. numu = nu - last + 1 if (numu .gt. 0) call dcopy ( numu, y(last), 1, U(lu), 1 ) * end of LUforw end