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(adsbygoogle = window.adsbygoogle || []).push({});

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 cant see where the bug is. Could anyone please have a look through my code to see if there are any obvious issues:

Code (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

**Physics Forums | Science Articles, Homework Help, Discussion**

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

# Fortran I can't find a small bug in this Fortran code

**Physics Forums | Science Articles, Homework Help, Discussion**