0
Lapack에서 SVD 분해를 검증 할 때 이상한 결과가 나타납니다. 이러한 루틴은 대개 견고하기 때문에 버그가 내 측면에 있다고 믿습니다. 어떤 도움을 주시면 감사하겠습니다. 내 행렬은 크기 n의 * n을 pentadiagonal이며 코드는 다음과 같습니다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))
편집 : 입력 배열 AB를 수정
! 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
감사