Search code examples
visual-studiofortrantransformationqr-decompositionplato

QR decomposition Fortran error


I have a problem with QR decomposition method. I use dgeqrf subroutine for decomposition but there is no error in the compiler but it gives a problem after that. I haven't found where is the mistake. Another question is, A=Q*R=> if A matrix has zero, Can decomposition be zero or lose the rank.

program decomposition

!CONTAINS
!subroutine Qrdecomposition(A_mat, R)
real,dimension(2,2)   :: A_mat    !real,dimension(2,2),intent(inout)   
:: A_mat
real,dimension(2,2)   :: R        !real,dimension(2,2),intent(out)     
:: R
real,dimension(2,2)                  :: A
integer                              :: M,N,LDA,LWORK,INFO
real,allocatable, dimension(:,:)     :: TAU
real,allocatable, dimension(:,:)     :: WORK
external   dgeqrf
M=2
N=2
LDA=2
LWORK=2
INFO=0
A_mat(1,1)=4
A_mat(1,2)=1
A_mat(2,1)=3
A_mat(2,2)=1
A=A_mat

call dgeqrf(M,N,A,TAU,WORK,LWORK,INFO)
R=A
print *,R,WORK,LWORK

!end subroutine Qrdecomposition
end program decomposition

Solution

  • I see three mistakes in your code:

    1) You forgot the LDA argument to dgeqrf,

    2) TAU and WORK must be explicitly allocated,

    3) All arrays should be declared with double precision to be consistent with dgeqrf interface:

    program decomposition
    
    !CONTAINS
    !subroutine Qrdecomposition(A_mat, R)
    ! Note: using '8' for the kind parameter is not the best style but I'm doing it here for brevity.
    real(8),dimension(2,2)   :: A_mat    !real,dimension(2,2),intent(inout)
    real(8),dimension(2,2)   :: R        !real,dimension(2,2),intent(out)
    real(8),dimension(2,2)                  :: A
    integer                              :: M,N,LDA,LWORK,INFO
    real(8),allocatable, dimension(:,:)     :: TAU
    real(8),allocatable, dimension(:,:)     :: WORK
    external   dgeqrf
    M=2
    N=2
    LDA=2
    LWORK=2
    INFO=0
    A_mat(1,1)=4
    A_mat(1,2)=1
    A_mat(2,1)=3
    A_mat(2,2)=1
    A=A_mat
    
    allocate(tau(M,N), work(M,N))
    call dgeqrf(M,N,A,LDA,TAU,WORK,LWORK,INFO)
    R=A
    print *,R,WORK,LWORK
    
    !end subroutine Qrdecomposition
    end program decomposition
    

    In certain situations, Fortran does perform automatic allocation of arrays, but it should generally not be counted on and it is not the case here.

    EDIT Point 3 was pointed out by roygvib, see below.