2

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
Nobody
  • 79
  • 13

1 Answers1

5

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.

Raul Laasner
  • 1,475
  • 1
  • 17
  • 30
  • Thanks for helping :) What did you find R(transformation matrix)? – Nobody May 22 '18 at 22:13
  • This is what it showed on my machine: `262144.906 -3.0 -6.49038268E+32 0.624999940`. – Raul Laasner May 23 '18 at 01:13
  • Isn't the kind of input arrays (real) and the subroutine dgeqrf (double precision) inconsistent? http://www.netlib.org/lapack/lapack-3.1.1/html/sgeqrf.f.html http://www.netlib.org/lapack/explore-3.1.1-html/dgeqrf.f.html – roygvib May 23 '18 at 04:57
  • if A=[4 1;3 1],the result should be R=[ -5.00 -1.40 ;0.00 0.20]. But the code result=[-5 0.33;-1.4 0.2] something is wrong because A(2,1) should be zero for orthogonal transformation – Nobody May 23 '18 at 07:17
  • Actually you are right, because I am using real but when I try sqeqrf, the result is so different – Nobody May 23 '18 at 07:18
  • @Nobody I think the 0.33 in the R(2,1) element represents part of the "v" vectors (see the manual pages above). That is, the returned A matrix contains the values of R in the upper triangular part (including the diagonals), while the lower triangular part of A gives info about the Q matrix (in terms of the v vectors). The correct (2,1) element of the R should be 0 (probably). – roygvib May 23 '18 at 07:40
  • Also please note that the 2-dim array "A" in Fortran is column-major, so if you print "A" (or "R") as "print *, A", then the values are printed in this order: A(1,1), A(2,1), A(1,2), A(2,2). [As for sgeqrf(), I got similar results when used with "real" arrays, so something seems strange...(note that types need to match with routines.)] – roygvib May 23 '18 at 07:42
  • @roygvib Yes you are right R(2,1) represents the "v" vector, How can I obtain only R=[-5 -14.4;0 0.2] because I don't need the diagonal element for QR factorization or do I need to use other subroutine for Qr factorization? – Nobody May 23 '18 at 08:06
  • @Nobody I think you can get R from the upper triangular part of the returned A matrix (by ignoring the lower triangular part), i.e., R(:,:) = 0 then R(i,j) = A(i,j) for all 1 <= i <= j <= N. – roygvib May 23 '18 at 08:21
  • @roygvib You are right that it should be double precision. I can't believe I missed that! – Raul Laasner May 23 '18 at 14:04
  • Thanks I changed it :) – Nobody May 25 '18 at 07:48