0
我在驗證Lapack的SVD分解時遇到了奇怪的結果。這些例程通常很健壯，所以我相信這個錯誤在我身邊。任何幫助將不勝感激。我的矩陣是五對，大小爲n×n個，並代碼如下：SVD分解後的矩陣重構
! 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))
編輯：添加變量聲明：
! 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
由於