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
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.