0

在验证来自 Lapack 的 SVD 分解时，我得到了奇怪的结果。这些例程通常很健壮，所以我相信这个错误就在我身边。任何帮助将不胜感激。我的矩阵是五边形的，大小为 n*n，代码如下：

! 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

谢谢

4

1 回答 1

2

例程 ZGBBRD 修改输入数组 AB。在调用例程之前，它应该保存在另一个数组中。使用这种预防措施看起来效果很好。谢谢。

于 2013-10-07T14:35:32.203 回答