The routine ZGBBRD modify the input array AB. It should be saved in another array before calling the routine. Looks like it works perfectly using this precaution. Thanks.
Matrix reconstruction after SVD decomposition
Question
I am having strange results when verifying the SVD decomposition from Lapack. Those routines are usually robust so I believe the bug is on my side. Any help will be highly appreciated. My matrix is pentadiagonal, size n*n and the code looks like :
! Compute real bi-diag from complex pentadiag
call ZGBBRD('B', n, n, 0, 2, 2, ab, 5, &
d, e, q, n, pt, n, dummy_argc, 1, work, rwork, info)
if (info.ne.0) then
print *,'Call to *GBBRD failed, info : ',info
call exit(0)
endif
! Compute diagonal matrix from bi-diagonal one
call dbdsdc('U', 'I', n, d, e, u, n, vt, n, &
dummy_argr, 1, work_big, iwork, info)
if (info.ne.0) then
print *,'Call to *BDSDC failed, info : ',info
call exit(0)
endif
print *,'singular values min/max : ',minval(d),maxval(d)
do ii=1,n
do jj=1,n
tmpqu(ii,jj)=0.
do kk=1,n
tmpqu(ii,jj)=tmpqu(ii,jj)+q(ii,kk)*u(kk,jj) ! Q*U
enddo
enddo
enddo
do ii=1,n
do jj=1,n
tmpqu(ii,jj)=tmpqu(ii,jj)*d(jj) ! Q*U*sigma
enddo
enddo
do ii=1,n
do jj=1,n
tmptot(ii,jj)=0.
do kk=1,n
tmptot(ii,jj) = tmptot(ii,jj) + & ! Q*U*sigma*VT
tmpqu(ii,kk)*vt(kk,jj)
enddo
enddo
enddo
tmpqu=tmptot
do ii=1,n
do jj=1,n
tmptot(ii,jj)=0.
do kk=1,n
tmptot(ii,jj) = tmptot(ii,jj) + & ! Q*U*sigma*VT*P
tmpqu(ii,kk)*pt(kk,jj)
enddo
enddo
enddo
tmpa=0.
do ii=1,n
tmpa(ii,ii)=ab(3,ii) ! diag
enddo
do ii=2,n
tmpa(ii-1,ii)=ab(2,ii) ! diag+1
enddo
do ii=3,n
tmpa(ii-2,ii)=ab(1,ii) ! diag+2
enddo
do ii=1,n-1
tmpa(ii+1,ii)=ab(4,ii) ! diag-1
enddo
do ii=1,n-2
tmpa(ii+2,ii)=ab(5,ii) ! diag-2
enddo
print *, 'maxabs delta',maxval(abs(tmptot-tmpa)), maxloc(abs(tmptot-tmpa))
EDIT : add variable declaration :
! Local variables
integer :: i,j,k
integer :: info, info2, code
! From pentadiagonale to bi-diagonale
complex(mytype), dimension(5,n) :: ab ! matrice pentadiagonale
real(mytype), dimension(n) :: d ! diagonale matrice bidiagonale
real(mytype), dimension(n-1) :: e ! diag+1 matrice bidiagonale
complex(mytype), dimension(n,n) :: q ! unitary matrix Q
complex(mytype), dimension(n,n) :: pt ! Unitary matrix P'
complex(mytype) :: dummy_argc
complex(mytype), dimension(n) :: work
real(mytype), dimension(n) :: rwork
! From bi-diagonale to SVD
real(mytype), dimension(n,n) :: u ! Left singular vectors
real(mytype), dimension(n,n) :: vt ! Right singular vectors
real(mytype) :: dummy_argr
real(mytype), dimension(3*n*n+4*n) :: work_big
integer, dimension(8*n) :: iwork
! Temp verif sigma
integer :: ii,jj,kk
complex(mytype), dimension(n,n) :: tmpa, tmpqu, tmptot
Thanks
Solution
Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow