Householder QR-factorization using Lapack and Fortran

Question

I have a problem with QR-factorization using dgeqrf with my fortran code. I need to create Q and R matrices, so that A = Q * R and then compute err=||A - Q*R|| / ||A||, where || || is a Frobenius norm (same to Euclidian for two-dimensional arrays). I've read about dgeqrf here https://www.netlib.org/lapack/explore-html/dd/d9a/group__double_g_ecomputational_ga3766ea903391b5cf9008132f7440ec7b.html. The mistake is definitely somewhere in my usage of the function, because I have a similar code with my own QR factrorization, where everything works. Please, help me.I think that there is a mistake in the program because err should be somewhere around 10^(-15) (here it is somewhere around 10 ^(0)), that is because A should be equal to Q * R with a small error.

program terra
    implicit none
    integer(4) :: n, diag = 1, h, wdth, u
    real(8), dimension(:, :), allocatable :: A, Q, R, M_aux
    real(8), dimension(:), allocatable :: arr, aux
    real(8) :: err
    n = 1024
    allocate (A(n,n), Q(n,n), R(n, n), M_aux(n, n), arr(n), aux(n))

    !filling matrices
    do h = 1, n
        do wdth = 1, n
            A(h, wdth) = 1 / real(h + wdth, 8)
            R(h,wdth) = A(h,wdth)
            Q(h, wdth) = 0
        end do
        Q(diag, diag) = 1
        diag = diag + 1
    end do

    !doing the factrorisation
    call dgeqrf(n, n, R, n , arr, aux, n, u)

    !getting Q matrix
    do h = 2, n
        diag = h - 1
        do wdth = 1, diag
            Q(h, wdth) = R(h, wdth)
            R(h, wdth) = 0
            !write(*, *) h, wdth, R(h,wdth), Q(h, wdth)
        end do
    end do

    !Transposing it
    do wdth = 1, n
        do h = 1, n
            M_aux(wdth, h) = Q(wdth, h)
        end do
    end do
    do wdth = 1, n
        do h = 1, n
            Q(wdth, h) = M_aux(h, wdth)
        end do
    end do

    !computing err
    do h = 1, n
        do wdth = 1, n
            M_aux(h, wdth) = 0
            do diag = 1, n
                M_aux(h, wdth) = M_aux(h, wdth) + Q(h, diag) * R(diag, wdth)
            end do
        end do
    end do
    M_aux = M_aux - A
    err = norm2(M_aux) / norm2(A)
    write(*,*) err
end program

Double check your formatting of your code, especially the first few lines. Either preface every line of code with 4 spaces, or wrap the entire thing in triple backticks. See formatting guide — blackbrandt
Please tell us why you think you code is wrong. What should it do for it to be considered right? — Ian Bush
err should be somewhere around 10^(-15), whereas it is somewhere around 10^0 — Layvis
RSI playing up so brief. DGEQRF doesn't return Q, you need to call DORG2R as well: netlib.org/lapack/explore-html/da/dba/… — Ian Bush

Layvis Layvis · Accepted Answer · 2020-12-02T19:27:13

program terra
    implicit none
    integer(4) :: n, diag = 1, h, wdth, u
    real(8), dimension(:, :), allocatable :: A, Q, R, M_aux
    real(8), dimension(:), allocatable :: arr, aux
    real(8) :: err
    n = 1024
    allocate (A(n,n), Q(n,n), R(n, n), M_aux(n, n), arr(n), aux(n), work(n))

    !filling matrices
    do h = 1, n
        do wdth = 1, n
            A(h, wdth) = 1 / real(h + wdth, 8)
            R(h,wdth) = A(h,wdth)
            Q(h, wdth) = 0
        end do
    end do

    !doing the factrorization
    call dgeqrf(n, n, R, n , arr, aux, n, u)
    Q = R

    !getting R matrix
    do h = 2, n
        diag = h - 1
        do wdth = 1, diag
            R(h, wdth) = 0
        end do
    end do

    !getting Q matrix
    call DORG2R(n, n, n, Q, n, arr, aux, u)

    !computing err
    do h = 1, n
        do wdth = 1, n
            M_aux(h, wdth) = 0
            do diag = 1, n
                M_aux(h, wdth) = M_aux(h, wdth) + Q(h, diag) * R(diag, wdth)
            end do
        end do
    end do
    M_aux = M_aux - A
    err = norm2(M_aux) / norm2(A)
    write(*,*) err
end program

Now this works. err is around 10^(-16). Ian Bush, my sincere thanks to you! Can I put a like or something like that on your comment? I don't understand how do it.

Householder QR-factorization using Lapack and Fortran

2 Answers