0

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

1 Answer 1

2

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.