Dismiss Notice
Join Physics Forums Today!
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

  1. May 3, 2015 #1

    O S

    User Avatar

    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 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
    Last edited: May 3, 2015
  2. jcsd
  3. May 3, 2015 #2


    User Avatar
    Gold Member

    "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.
    Last edited by a moderator: May 4, 2015
  4. May 4, 2015 #3
    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! :-)
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook