subroutine lu6sol( mode, m, n, v, w, &
lena, luparm, parmlu, &
a, indc, indr, p, q, &
lenc, lenr, locc, locr, &
inform )
integer(ip), intent(in) :: mode, m, n, lena
integer(ip), intent(in) :: indc(lena), indr(lena), p(m), q(n), &
lenc(n), lenr(m), locc(n), locr(m)
real(rp), intent(in) :: a(lena)
integer(ip), intent(inout) :: luparm(30)
real(rp), intent(inout) :: parmlu(30), v(m), w(n)
integer(ip), intent(out) :: inform
!-----------------------------------------------------------------------
! lu6sol uses the factorization A = L U as follows:
!
! mode
! 1 v solves L v = v(input). w is not touched.
! 2 v solves L'v = v(input). w is not touched.
! 3 w solves U w = v. v is not altered.
! 4 v solves U'v = w. w is destroyed.
! 5 w solves A w = v. v is altered as in 1.
! 6 v solves A'v = w. w is destroyed.
!
! If mode = 3,4,5,6, v and w must not be the same arrays.
!
! If lu1fac has just been used to factorize a symmetric matrix A
! (which must be definite or quasi-definite), the factors A = L U
! may be regarded as A = LDL', where D = diag(U). In such cases,
!
! mode
! 7 v solves A v = L D L'v = v(input). w is not touched.
! 8 v solves L |D| L'v = v(input). w is not touched.
!
! p(*), q(*) hold row and column numbers in pivotal order.
! lenc(k) is the length of the k-th column of initial L.
! lenr(i) is the length of the i-th row of U.
! locc(*) is not used.
! locr(i) is the start of the i-th row of U.
!
! U is assumed to be in upper-trapezoidal form (nrank by n).
! The first entry for each row is the diagonal element
! (according to the permutations p, q). It is stored at
! location locr(i) in a(*), indr(*).
!
! On exit, inform = 0 except as follows.
! If mode = 3,4,5,6 and if U (and hence A) is singular, then
! inform = 1 if there is a nonzero residual in solving the system
! involving U. parmlu(20) returns the norm of the residual.
!
! July 1987: Early version.
! 09 May 1988: f77 version.
! 27 Apr 2000: Abolished the dreaded "computed go to".
! But hard to change other "go to"s to "if then else".
! 15 Dec 2002: lu6L, lu6Lt, lu6U, lu6Ut added to modularize lu6sol.
! 13 Dec 2011: First f90 version.
!--------------------------------------------------------------------
if (mode == 1) then ! Solve L v(new) = v.
call lu6L ( inform, m, n, v, &
lena, luparm, parmlu, a, indc, indr, lenc )
else if (mode == 2) then ! Solve L'v(new) = v.
call lu6Lt ( inform, m, n, v, &
lena, luparm, parmlu, a, indc, indr, lenc )
else if (mode == 3) then ! Solve U w = v.
call lu6U ( inform, m, n, v, w, &
lena, luparm, parmlu, a, indr, p, q, lenr, locr )
else if (mode == 4) then ! Solve U'v = w.
call lu6Ut ( inform, m, n, v, w, &
lena, luparm, parmlu, a, indr, p, q, lenr, locr )
else if (mode == 5) then ! Solve A w = v
! via L v(new) = v
! and U w = v(new).
call lu6L ( inform, m, n, v, &
lena, luparm, parmlu, a, indc, indr, lenc )
call lu6U ( inform, m, n, v, w, &
lena, luparm, parmlu, a, indr, p, q, lenr, locr )
else if (mode == 6) then ! Solve A'v = w
! via U'v = w
! and L'v(new) = v.
call lu6Ut ( inform, m, n, v, w, &
lena, luparm, parmlu, a, indr, p, q, lenr, locr )
call lu6Lt ( inform, m, n, v, &
lena, luparm, parmlu, a, indc, indr, lenc )
else if (mode == 7) then ! Solve LDv(bar) = v
! and L'v(new) = v(bar).
call lu6LD ( inform, i1, m, n, v, &
lena, luparm, parmlu, a, indc, indr, lenc, locr )
call lu6Lt ( inform, m, n, v, &
lena, luparm, parmlu, a, indc, indr, lenc )
else if (mode == 8) then ! Solve L|D|v(bar) = v
! and L'v(new) = v(bar).
call lu6LD ( inform, i2, m, n, v, &
lena, luparm, parmlu, a, indc, indr, lenc, locr )
call lu6Lt ( inform, m, n, v, &
lena, luparm, parmlu, a, indc, indr, lenc )
end if
end subroutine lu6sol