- #1
O S
- 1
- 0
So I'm writing a short subroutine to compute a QR factorization of an (mxn) matrix using the Householder method. To test this, I'm computing the factorization of the matrix
12 -51 4
A = 6 167 -68 .
-4 24 -41
My householder method should return the following upper triangular matrix:
14 21 -14
R = 0 175 -70.
0 0 35
However, what I actually get is
-14 -21 14
R = 0 -175 70.
-0 0 35
So as you can see, it looks almost correct, but with some strange sign changes. I've been staring at the code for hours trying to see where these sign changes are being introduced, but I just can't see where the bug is. Could anyone please have a look through my code to see if there are any obvious issues:
12 -51 4
A = 6 167 -68 .
-4 24 -41
My householder method should return the following upper triangular matrix:
14 21 -14
R = 0 175 -70.
0 0 35
However, what I actually get is
-14 -21 14
R = 0 -175 70.
-0 0 35
So as you can see, it looks almost correct, but with some strange sign changes. I've been staring at the code for hours trying to see where these sign changes are being introduced, but I just can't see where the bug is. Could anyone please have a look through my code to see if there are any obvious issues:
Fortran:
subroutine full_QR(A, Q, m, n)
!==========================================================
! Subroutine for full QR factorization using the
! Householder Method.
!==========================================================
integer, parameter :: dp = selected_real_kind(15)
integer, intent(in) :: m, n
real(dp), dimension(m,n), intent(inout) :: A
real(dp), dimension(m,m), intent(out) :: Q
integer :: k
real(dp) :: two_norm
real(dp), dimension(m) :: x, e1
real(dp), dimension(m,n) :: v
real(dp), dimension(m,m) :: outprod_vv
v = 0.0_dp
do k=1,m
Q(k,k) = 1.0_dp
end do
!Householder triangularization
do k=1,n
e1(k) = 1.0_dp
x(k:m) = A(k:m,k)
v(k:m,k) = sign( sqrt(dot_product(x(k:m),x(k:m))), x(k) )* &
e1(k:m) + x(k:m)
v(k:m,k) = v(k:m,k)/(sqrt(dot_product(v(k:m,k),v(k:m,k))))
call outer_product(v(k:m,k), m-k+1, outprod_vv(k:m,k:m))
A(k:m,k:n) = A(k:m,k:n) - &
2.0_dp*matmul(outprod_vv(k:m,k:m), A(k:m,k:n))
!Form Q implicitly
Q(k:m,k:m) = Q(k:m,k:m) - 2.0_dp* &
matmul(outprod_vv(k:m,k:m), Q(k:m,k:m))
end do
Q = transpose(Q)
end subroutine full_QR
Last edited: