given an n by n upper triangular matrix r, this subroutine computes the qr decomposition of the matrix formed when a row is added to r. if the row is specified by the vector w, then rwupdt determines an orthogonal matrix q such that when the n+1 by n matrix composed of r augmented by w is premultiplied by (q transpose), the resulting matrix is upper trapezoidal. the matrix (q transpose) is the product of n transformations
g(n)*g(n-1)* ... *g(1)
where g(i) is a givens rotation in the (i,n+1) plane which eliminates elements in the (n+1)-st plane. rwupdt also computes the product (q transpose)*c where c is the (n+1)-vector (b,alpha). q itself is not accumulated, rather the information to recover the g rotations is supplied.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | n | a positive integer input variable set to the order of r. |
||
| real(kind=wp), | intent(inout) | :: | r(Ldr,n) | an n by n array. on input the upper triangular part of r must contain the matrix to be updated. on output r contains the updated triangular matrix. |
||
| integer, | intent(in) | :: | Ldr | a positive integer input variable not less than n which specifies the leading dimension of the array r. |
||
| real(kind=wp), | intent(in) | :: | w(n) | an input array of length n which must contain the row vector to be added to r. |
||
| real(kind=wp), | intent(inout) | :: | b(n) | an array of length n. on input b must contain the first n elements of the vector c. on output b contains the first n elements of the vector (q transpose)*c. |
||
| real(kind=wp), | intent(inout) | :: | Alpha | a variable. on input alpha must contain the (n+1)-st element of the vector c. on output alpha contains the (n+1)-st element of the vector (q transpose)*c. |
||
| real(kind=wp), | intent(out) | :: | Cos(n) | an output array of length n which contains the cosines of the transforming givens rotations. |
||
| real(kind=wp), | intent(out) | :: | Sin(n) | an output array of length n which contains the sines of the transforming givens rotations. |
subroutine rwupdt(n, r, Ldr, w, b, Alpha, Cos, Sin)
implicit none
integer, intent(in) :: n !! a positive integer input variable set to the order of r.
integer, intent(in) :: Ldr !! a positive integer input variable not less than n
!! which specifies the leading dimension of the array r.
real(wp), intent(inout) :: Alpha !! a variable. on input alpha must contain the
!! (n+1)-st element of the vector c. on output alpha contains
!! the (n+1)-st element of the vector (q transpose)*c.
real(wp), intent(inout) :: r(Ldr, n) !! an n by n array. on input the upper triangular part of
!! r must contain the matrix to be updated. on output r
!! contains the updated triangular matrix.
real(wp), intent(in) :: w(n) !! an input array of length n which must contain the row
!! vector to be added to r.
real(wp), intent(inout) :: b(n) !! an array of length n. on input b must contain the
!! first n elements of the vector c. on output b contains
!! the first n elements of the vector (q transpose)*c.
real(wp), intent(out) :: Cos(n) !! an output array of length n which contains the
!! cosines of the transforming givens rotations.
real(wp), intent(out) :: Sin(n) !! an output array of length n which contains the
!! sines of the transforming givens rotations.
integer :: i, j, jm1
real(wp) :: cotan, rowj, tan, temp
real(wp), parameter :: p5 = 5.0e-1_wp
real(wp), parameter :: p25 = 2.5e-1_wp
do j = 1, n
rowj = w(j)
jm1 = j - 1
! apply the previous transformations to
! r(i,j), i=1,2,...,j-1, and to w(j).
if (jm1 >= 1) then
do i = 1, jm1
temp = Cos(i)*r(i, j) + Sin(i)*rowj
rowj = -Sin(i)*r(i, j) + Cos(i)*rowj
r(i, j) = temp
end do
end if
! determine a givens rotation which eliminates w(j).
Cos(j) = one
Sin(j) = zero
if (rowj /= zero) then
if (abs(r(j, j)) >= abs(rowj)) then
tan = rowj/r(j, j)
Cos(j) = p5/sqrt(p25 + p25*tan**2)
Sin(j) = Cos(j)*tan
else
cotan = r(j, j)/rowj
Sin(j) = p5/sqrt(p25 + p25*cotan**2)
Cos(j) = Sin(j)*cotan
end if
! apply the current transformation to r(j,j), b(j), and alpha.
r(j, j) = Cos(j)*r(j, j) + Sin(j)*rowj
temp = Cos(j)*b(j) + Sin(j)*Alpha
Alpha = -Sin(j)*b(j) + Cos(j)*Alpha
b(j) = temp
end if
end do
end subroutine rwupdt