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

[O S]

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:

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

 

[phinds]

dfg

"dfg" is your whole set of code? Doesn't seem like it's going to do much

Mod note (Mark44): I have restored the next-to-last version of the post from O S.

 

[Carno Raar]

Staring at it for hours? Would be quicker to change a sign, recompile and repeat until it worked.

Sometimes programming isn't as clever as they make it out to be in the adverts. We have a whole discipline (TDD) devoted to building effective tests, just so that we could debug in this manner! :-)