不准确的 zheev 特征值和向量
Inaccurate zheev eigen values and vectors
我在 Fortran 科学代码中使用 LAPACK zheev 例程来计算不太大（可能永远不会超过大小 1000）的矩阵的特征值和向量。
由于这一步发生在计算的开始，我必须达到很高的准确性以避免重要的错误传播。问题是在我的测试用例中（仅使用 12x12 矩阵）计算的精度仅为 1e-9 左右，这根本不够。
我与 numpy.linalg.eigh 进行了比较，结果好得离谱，我想知道如何使用 zheev 实现如此精确。
这里是一个测试 Python 代码，其中一些来自 Fortran 代码的输出作为数据：
from numpy import array, linalg, diagflat
from sys import argv
prec = 10
if len(argv) > 1:
prec = int(argv[1])
tolZHEEV = 0.1**prec
tolEIGH = tolZHEEV
if len(argv) > 2:
prec2 = int(argv[2])
tolEIGH = 0.1**prec2
F = array([[ 1.99964605e+00+0.00000000e+00j, -6.77332886e-04-1.58965282e-04j,
2.50377192e-16-7.43492040e-16j, -1.03545717e-03-4.75346808e-04j,
-1.33562677e-16-7.96956118e-17j, 2.69297453e-04+6.01705394e-05j,
-2.42109713e-04-5.83896237e-05j, 1.03264964e-16-4.49556750e-17j,
-9.63072787e-04-1.58486337e-04j, 1.32803698e-16-6.77095449e-17j,
3.28130435e-15-1.74181973e-15j, 1.28210970e-02+2.86865646e-03j],
[-6.77332886e-04+1.58965282e-04j, 1.99814145e+00+0.00000000e+00j,
3.11074872e-16-2.03650443e-15j, -2.86327278e-03-5.79965447e-04j,
4.87808786e-19-4.06383000e-17j, 7.07826164e-04-7.57149613e-06j,
-6.39134657e-04-3.91798299e-06j, -3.05291091e-17-5.12700538e-17j,
-2.49967233e-03+1.68782808e-04j, 8.54644546e-16-1.36907397e-15j,
4.06124564e-15-6.47961846e-15j, 3.32552948e-02-3.45914007e-04j],
[ 2.50377192e-16+7.43492040e-16j, 3.11074872e-16+2.03650443e-15j,
1.99801950e+00+0.00000000e+00j, -1.38347852e-16+1.96637141e-15j,
-1.50312544e-03-2.16900444e-05j, -7.20316942e-17-8.68605005e-16j,
3.19121715e-17+4.14687698e-16j, 5.27702831e-03-4.38145588e-05j,
5.69734505e-16+2.88527199e-15j, -3.28447952e-02-6.09751924e-04j,
2.01063189e-02+2.33609810e-04j, -8.18479961e-15-4.06359927e-14j],
[-1.03545717e-03+4.75346808e-04j, -2.86327278e-03+5.79965447e-04j,
-1.38347852e-16-1.96637141e-15j, 1.99514423e+00+0.00000000e+00j,
8.58581114e-17+1.24995999e-15j, 1.13414603e-03-2.42381989e-04j,
-1.02783476e-03+2.01640217e-04j, -1.40117884e-16-3.89960963e-15j,
-3.96054944e-03+1.08447880e-03j, 3.89988223e-16+2.05430996e-14j,
3.71606023e-15-2.76569268e-14j, 5.28364719e-02-1.12755531e-02j],
[-1.33562677e-16+7.96956118e-17j, 4.87808786e-19+4.06383000e-17j,
-1.50312544e-03+2.16900444e-05j, 8.58581114e-17-1.24995999e-15j,
1.98788997e+00+0.00000000e+00j, 2.43606700e-16-2.72581191e-16j,
2.00791472e-17+2.16983168e-16j, 1.39389509e-03-3.16910018e-05j,
-3.80959351e-17+8.54418747e-17j, 3.31334855e-02+1.36959698e-04j,
7.93487223e-02-2.23030670e-04j, -1.08947296e-14-1.69482780e-14j],
[ 2.69297453e-04-6.01705394e-05j, 7.07826164e-04+7.57149613e-06j,
-7.20316942e-17+8.68605005e-16j, 1.13414603e-03+2.42381989e-04j,
2.43606700e-16+2.72581191e-16j, 1.99966113e+00+0.00000000e+00j,
2.54460391e-04+4.28207738e-06j, -1.04566138e-16-9.67432153e-17j,
9.96762100e-04-5.66002680e-05j, -3.99532806e-16+4.90873168e-16j,
-3.30484846e-15+6.18171590e-16j, -1.25976974e-02-3.71663247e-06j],
[-2.42109713e-04+5.83896237e-05j, -6.39134657e-04+3.91798299e-06j,
3.19121715e-17-4.14687698e-16j, -1.02783476e-03-2.01640217e-04j,
2.00791472e-17-2.16983168e-16j, 2.54460391e-04-4.28207738e-06j,
1.99972038e+00+0.00000000e+00j, 2.08257020e-16-1.01061133e-15j,
-9.01577811e-04+6.64306049e-05j, -1.55427085e-15+7.39572778e-15j,
1.64292158e-15-2.72307945e-15j, 1.09937571e-02-1.81759386e-04j],
[ 1.03264964e-16+4.49556750e-17j, -3.05291091e-17+5.12700538e-17j,
5.27702831e-03+4.38145588e-05j, -1.40117884e-16+3.89960963e-15j,
1.39389509e-03+3.16910018e-05j, -1.04566138e-16+9.67432153e-17j,
2.08257020e-16+1.01061133e-15j, 1.98462892e+00+0.00000000e+00j,
1.29258085e-16-9.03812998e-17j, 1.06287458e-01+2.85612157e-03j,
-3.58021992e-02-7.13306294e-04j, 3.48395880e-15+3.39493033e-15j],
[-9.63072787e-04+1.58486337e-04j, -2.49967233e-03-1.68782808e-04j,
5.69734505e-16-2.88527199e-15j, -3.96054944e-03-1.08447880e-03j,
-3.80959351e-17-8.54418747e-17j, 9.96762100e-04+5.66002680e-05j,
-9.01577811e-04-6.64306049e-05j, 1.29258085e-16+9.03812998e-17j,
1.99640522e+00+0.00000000e+00j, 3.20322529e-16-2.99500218e-15j,
6.23226621e-15-7.41164282e-15j, 4.21794574e-02+2.40760808e-03j],
[ 1.32803698e-16+6.77095449e-17j, 8.54644546e-16+1.36907397e-15j,
-3.28447952e-02+6.09751924e-04j, 3.89988223e-16-2.05430996e-14j,
3.31334855e-02-1.36959698e-04j, -3.99532806e-16-4.90873168e-16j,
-1.55427085e-15-7.39572778e-15j, 1.06287458e-01-2.85612157e-03j,
3.20322529e-16+2.99500218e-15j, 5.44503041e-01+0.00000000e+00j,
3.71105458e-04-2.57711214e-06j, 4.58329172e-16-1.52262708e-16j],
[ 3.28130435e-15+1.74181973e-15j, 4.06124564e-15+6.47961846e-15j,
2.01063189e-02-2.33609810e-04j, 3.71606023e-15+2.76569268e-14j,
7.93487223e-02+2.23030670e-04j, -3.30484846e-15-6.18171590e-16j,
1.64292158e-15+2.72307945e-15j, -3.58021992e-02+7.13306294e-04j,
6.23226621e-15+7.41164282e-15j, 3.71105458e-04+2.57711214e-06j,
3.21340892e-01+0.00000000e+00j, -1.54683439e-15+3.05485859e-17j],
[ 1.28210970e-02-2.86865646e-03j, 3.32552948e-02+3.45914007e-04j,
-8.18479961e-15+4.06359927e-14j, 5.28364719e-02+1.12755531e-02j,
-1.08947296e-14+1.69482780e-14j, -1.25976974e-02+3.71663247e-06j,
1.09937571e-02+1.81759386e-04j, 3.48395880e-15-3.39493033e-15j,
4.21794574e-02-2.40760808e-03j, 4.58329172e-16+1.52262708e-16j,
-1.54683439e-15-3.05485859e-17j, 3.20142168e-01+0.00000000e+00j]])
d =array([[1.99995093],
[1.9999434 ],
[1.99993708],
[1.99992886],
[1.99992603],
[1.99991767],
[1.99477538],
[1.99278704],
[1.98991992],
[0.53521096],
[0.3165472 ],
[0.31639848]])
D = diagflat(d)
P = array([[ 8.50795081e-04+1.90361152e-04j, 6.41403660e-03+1.43510867e-03j,
2.93305520e-02+6.56256458e-03j, -2.73319405e-12-4.62616407e-12j,
2.36860608e-02+5.29963787e-03j, 9.63512391e-01+2.15581088e-01j,
-1.18829015e-14+6.38499316e-15j, 1.49823188e-01+3.35221904e-02j,
4.63016072e-15-4.09382827e-15j, 5.35804136e-17-2.23033645e-17j,
-7.22943096e-14-1.78616099e-14j, -7.65412277e-03-1.71257176e-03j],
[ 5.96040942e-03-6.19988216e-05j, 4.13411608e-02-4.30021346e-04j,
1.97753784e-01-2.05698985e-03j, 1.61937250e-11+7.38548825e-12j,
8.90791562e-01-9.26581108e-03j, -8.99773456e-02+9.35923868e-04j,
-6.78075986e-14+2.27240115e-14j, 3.96283860e-01-4.12205452e-03j,
3.62263322e-14+2.70989836e-14j, 2.80882500e-17-1.14706336e-15j,
-1.89947123e-13-2.17176505e-15j, -1.98547508e-02+2.06524599e-04j],
[ 1.43774119e-12-2.04841175e-13j, -3.84634090e-12+1.02008935e-12j,
2.14225297e-12-1.58626154e-12j, 5.76965748e-01+7.48380143e-01j,
-1.64985984e-11-8.05494677e-12j, 5.49834259e-12+8.08747081e-13j,
1.64047208e-01+6.01119828e-02j, -8.41828091e-14-4.52343007e-13j,
1.15974306e-01+2.49814055e-01j, -1.89718719e-02+1.21355877e-02j,
1.21027166e-02-1.40334323e-03j, -1.12121788e-13+3.78501430e-14j],
[ 1.49140349e-02-3.18272560e-03j, 1.07045322e-01-2.28439783e-02j,
5.83924798e-01-1.24612315e-01j, -7.68362748e-12-3.20504318e-12j,
-4.29872936e-01+9.17369184e-02j, -1.07323514e-01+2.29033458e-02j,
2.52998000e-14-3.37978544e-14j, 6.37564522e-01-1.36059286e-01j,
2.28120526e-13-8.82579171e-14j, 8.15806584e-15+1.14028011e-14j,
-3.05058814e-13+4.83321942e-14j, -3.15470869e-02+6.73229760e-03j],
[-1.50646374e-13+4.30474354e-14j, 2.88653601e-13-7.82348111e-14j,
-1.56416337e-13+6.81259624e-14j, -4.56052348e-02-5.74215381e-02j,
1.24973393e-12+6.35459037e-13j, -4.16040415e-13-4.47866383e-14j,
-6.68174568e-01-2.33961052e-01j, -7.77941534e-14-1.12876492e-13j,
3.04089801e-01+6.31021660e-01j, 1.92382317e-02-1.27008073e-02j,
4.69043183e-02-6.12576169e-03j, -4.46806257e-13+6.93146032e-14j],
[-2.00864022e-02-5.92598703e-06j, -9.12820255e-01-2.69304562e-04j,
3.72566080e-01+1.09916213e-04j, -1.06450784e-12-4.32276375e-12j,
3.36150018e-02+9.91725789e-06j, 1.90596074e-02+5.62305535e-06j,
3.49708145e-14+4.34006214e-15j, -1.61246105e-01-4.75715912e-05j,
8.38875458e-15-4.53541474e-15j, -2.47012312e-16+4.74760834e-16j,
7.10584649e-14+6.19010235e-16j, 7.52355819e-03+2.21963585e-06j],
[ 9.62420518e-01-1.59116634e-02j, -1.24950681e-01+2.06580506e-03j,
-1.86710046e-01+3.08687039e-03j, -1.30253808e-12-1.41868926e-12j,
-2.69408788e-02+4.45412570e-04j, -1.67620621e-02+2.77126563e-04j,
-1.12422533e-16-5.40641968e-14j, 1.48156214e-01-2.44946130e-03j,
5.33921595e-14+3.20699063e-14j, 2.10253369e-15+4.74023210e-15j,
-6.27766220e-14-7.29612349e-16j, -6.56715101e-03+1.08574469e-04j],
[ 5.05414777e-13+9.40307821e-15j, -1.32687279e-12+3.49123357e-13j,
7.54773974e-13-3.92282159e-13j, 1.92570500e-01+2.54118137e-01j,
-5.58410396e-12-2.86684025e-12j, 1.83884514e-12+2.64297906e-13j,
-6.36994666e-01-2.39431942e-01j, 1.93041775e-14-4.12911291e-14j,
-2.70995935e-01-5.96659039e-01j, 6.18161118e-02-3.88219877e-02j,
-2.18180305e-02+2.34644894e-03j, 2.08901179e-13-2.47813917e-14j],
[-2.69415450e-01-1.53782636e-02j, -3.70083377e-01-2.11244000e-02j,
-6.54530261e-01-3.73606595e-02j, -1.63145928e-12-8.78784663e-14j,
-1.03526300e-01-5.90929264e-03j, -6.61072940e-02-3.77341164e-03j,
-4.04464076e-14+5.89878883e-14j, 5.85724599e-01+3.34332248e-02j,
3.15506629e-14+2.84647065e-14j, -1.09481211e-15-1.77104489e-15j,
-2.40425891e-13-1.87628386e-14j, -2.51994911e-02-1.43838973e-03j],
[ 2.54752203e-16+2.12657752e-16j, -3.21054014e-15+3.13509411e-16j,
1.50876469e-15-4.65388714e-16j, 1.84339353e-04+2.30141854e-04j,
-7.60919954e-15-2.39290177e-15j, 1.07829490e-15+1.72556420e-16j,
-6.61876184e-02-2.28688765e-02j, 9.75386453e-17-3.60186046e-15j,
-1.68109710e-02-3.45191353e-02j, -8.29576318e-01+5.52611651e-01j,
6.68078399e-03-9.00620322e-04j, -6.62742104e-14+9.74645368e-15j],
[-4.83726841e-16-3.76055196e-17j, -3.28302552e-15+1.41917600e-15j,
2.80242692e-15-7.23389671e-16j, 6.51697903e-04+8.25305197e-04j,
-1.99649588e-14-9.17555906e-15j, 7.51565389e-15+9.10141876e-16j,
-1.59564626e-02-5.63754717e-03j, 2.04075321e-15+5.12238396e-15j,
2.18750663e-02+4.57214672e-02j, -6.28420329e-03+4.12342677e-03j,
-9.90493516e-01+1.26529179e-01j, 9.57441797e-12-1.22296432e-12j],
[ 2.79384410e-04+0.00000000e+00j, 1.09500742e-03+0.00000000e+00j,
2.84867177e-03+0.00000000e+00j, 9.54279829e-15+0.00000000e+00j,
6.53433434e-04+0.00000000e+00j, 4.93180792e-04+0.00000000e+00j,
1.04504620e-15+0.00000000e+00j, 4.71495396e-02+0.00000000e+00j,
-7.15402040e-15+0.00000000e+00j, -1.97732525e-15+0.00000000e+00j,
9.65544745e-12+0.00000000e+00j, 9.98882805e-01+0.00000000e+00j]])
fp = F.dot(P) # FP = PD
dp = P.dot(D)
dd, PP = linalg.eigh(-F, 'U')
dd = -dd
DD = diagflat(dd)
ffp = F.dot(PP) # FP = PD
ddp = PP.dot(DD)
zer = fp - dp
zernp = ffp - ddp
print("\n== Innacurate (module > " + str(tolZHEEV) + ") commponants for zheev")
for n in range(12):
for p in range(12):
if abs(zer[n,p]) > tolZHEEV:
print(n, p, zer[n,p])
print("\n== Innacurate (module > " + str(tolEIGH) + ") commponants for eigh")
for n in range(12):
for p in range(12):
if abs(zernp[n,p]) > tolEIGH:
print(n, p, zernp[n,p])
print("\n== Error on eigen values")
for i in range(dd.shape[0]):
print(dd[i], dd[i] - d[i])
然后按照两个fortran文件重现数据：
主.F90
program main
implicit none
integer, parameter :: nband = 12
complex*16 :: occ_nd_cpx(nband, nband)
real*8 :: occ_diag(nband), a, b
integer :: i, j
do i=1,nband
do j=1,nband
read (111,'(2D25.16)') a, b
occ_nd_cpx(i,j) = cmplx(a, b, kind=8)
end do
end do
call diag_occ(occ_nd_cpx, nband, occ_diag)
end program main
diag_occ.F90
subroutine diag_occ(occ_nd_cpx, nband, occ_diag)
implicit none
integer, parameter :: dp=kind(1.0d0), dpc = dp
integer,intent(in) :: nband
complex(kind=dpc), intent(inout) :: occ_nd_cpx(nband, nband)
real(kind=dp), intent(out) :: occ_diag(nband)
integer :: info, lwork, n
character(len=500) :: message
real(kind=dp) :: rwork(3*nband-1)
complex(kind=dpc), allocatable :: work(:)
do n=1,nband
write (*,*) '!!!1', occ_nd_cpx(n, :)
end do
!! Initialisation
rwork = 0.0_dp
info = 0
!! use the opposite to have zheev orders the eigenvalues descending (once the
!! opposite have been taken again)
occ_nd_cpx = -occ_nd_cpx
!! Get diagonal occupations and associeted base
! Compute the optimal working array size
allocate(work(1))
work = (0.0_dp, 0.0_dp)
call zheev('V', 'U', nband, occ_nd_cpx, nband, occ_diag, work, -1, rwork, info)
lwork = work(1)
write (*,*) 'lwork, work(1)', lwork, work(1)
deallocate(work)
! Compute the eigenvalues (occ_diag) and vectors
allocate(work(lwork))
work = (0.0_dp, 0.0_dp)
call zheev('V', 'U', nband, occ_nd_cpx, nband, occ_diag, work, lwork, rwork, info)
!! obtain the true eigen values of occupation matrix in descending order
occ_diag = -occ_diag
do n=1,nband
write (*,*) '!!!2', occ_diag(n)
end do
do n=1,nband
write (*,*) '!!!3', occ_nd_cpx(n, :)
end do
info = 0
deallocate(work)
end subroutine diag_occ
我用gfortran -ffree-form -llapack main.F90 diag_occ.F90 -o main编译
LAPCK 库应该来自 CentOS 存储库 (atlas)
我使用以下脚本将主可执行文件的输出转换为 numpy 数组：
fortran2num.py
import sys
import numpy as np
## Format
FSIZE = 25
CPX = True
N = 12
IGNORE = 0
if len(sys.argv) > 1 and sys.argv[1] == '-h':
print(sys.argv[0] + ' [FSIZE [N [CMPX]]]')
exit()
if len(sys.argv) > 1:
FSIZE = int(sys.argv[1])
if len(sys.argv) > 2:
N = int(sys.argv[2])
if len(sys.argv) > 3:
CPX = sys.argv[3] in ['true', 'True', '1', 'T', 't']
if len(sys.argv) > 4:
IGNORE = int(sys.argv[4])
src = sys.stdin.read()
lst = [[]]
c = IGNORE
k = 0
while c < len(src):
k += 1
if CPX:
step = 2*FSIZE+3
a,b = src[c+1:c+step-1].split(',')
lst[-1].append(float(a)+1j*float(b))
else:
step = FSIZE
a = src[c:c+step]
lst[-1].append(float(a))
c += step + 1
if k >= N:
k = 0
c += IGNORE
lst.append([])
if len(lst[-1]) == 0:
lst.pop()
print(repr(np.array(lst)))
我使用这些命令创建了 F、d 和 P：
cat f | grep '!!!1' | ./fortran2num.py 25 12 T 6
cat f | grep '!!!2' | ./fortran2num.py 25 1 F 6
cat f | grep '!!!3' | ./fortran2num.py 25 12 T 6
编辑：
上次编辑这个我很匆忙，忘记了一些数据。这是主要可执行文件的输入 fortran.111 文件：
1.9996460536245420 0.0000000000000000
-6.7733288642473765E-004 -1.5896528158261059E-004
2.5037719177319647E-016 -7.4349203998305547E-016
-1.0354571719663188E-003 -4.7534680782504999E-004
-1.3356267732320494E-016 -7.9695611788720466E-017
2.6929745274964923E-004 6.0170539417274264E-005
-2.4210971257827060E-004 -5.8389623659711954E-005
1.0326496350503192E-016 -4.4955675025152728E-017
-9.6307278676681905E-004 -1.5848633659367058E-004
1.3280369753342969E-016 -6.7709544933739064E-017
3.2813043510687810E-015 -1.7418197254477417E-015
1.2821097042473820E-002 2.8686564596901464E-003
-6.7733288642473765E-004 1.5896528158261059E-004
1.9981414542927907 0.0000000000000000
3.1107487236320250E-016 -2.0365044300513661E-015
-2.8632727795318471E-003 -5.7996544650846344E-004
4.8780878630335091E-019 -4.0638299964232169E-017
7.0782616404454677E-004 -7.5714961311373116E-006
-6.3913465705948242E-004 -3.9179829909342934E-006
-3.0529109137315778E-017 -5.1270053829163420E-017
-2.4996723334157079E-003 1.6878280752610168E-004
8.5464454629174907E-016 -1.3690739743679114E-015
4.0612456352043776E-015 -6.4796184564793087E-015
3.3255294771216873E-002 -3.4591400749335414E-004
2.5037719177319647E-016 7.4349203998305547E-016
3.1107487236320250E-016 2.0365044300513661E-015
1.9980194987385831 0.0000000000000000
-1.3834785212208944E-016 1.9663714097332382E-015
-1.5031254367981852E-003 -2.1690044364606380E-005
-7.2031694214561962E-017 -8.6860500453088960E-016
3.1912171522844715E-017 4.1468769820119327E-016
5.2770283128296177E-003 -4.3814558771467556E-005
5.6973450505311861E-016 2.8852719875277487E-015
-3.2844795244561291E-002 -6.0975192388691354E-004
2.0106318907807723E-002 2.3360980952464097E-004
-8.1847996124173807E-015 -4.0635992692925845E-014
-1.0354571719663188E-003 4.7534680782504999E-004
-2.8632727795318471E-003 5.7996544650846344E-004
-1.3834785212208944E-016 -1.9663714097332382E-015
1.9951442266000110 0.0000000000000000
8.5858111362132558E-017 1.2499599851488651E-015
1.1341460255180225E-003 -2.4238198932919101E-004
-1.0278347636459521E-003 2.0164021732396069E-004
-1.4011788447525168E-016 -3.8996096313994688E-015
-3.9605494357122413E-003 1.0844788031769959E-003
3.8998822316673752E-016 2.0543099601558386E-014
3.7160602329464569E-015 -2.7656926837772642E-014
5.2836471927730132E-002 -1.1275553097512523E-002
-1.3356267732320494E-016 7.9695611788720466E-017
4.8780878630335091E-019 4.0638299964232169E-017
-1.5031254367981852E-003 2.1690044364606380E-005
8.5858111362132558E-017 -1.2499599851488651E-015
1.9878899727642418 0.0000000000000000
2.4360669979466240E-016 -2.7258119127922638E-016
2.0079147210768701E-017 2.1698316754064389E-016
1.3938950857567451E-003 -3.1691001834563675E-005
-3.8095935113214172E-017 8.5441874687637690E-017
3.3133485502343657E-002 1.3695969769169161E-004
7.9348722336522293E-002 -2.2303067036207741E-004
-1.0894729641781911E-014 -1.6948278015743026E-014
2.6929745274964923E-004 -6.0170539417274264E-005
7.0782616404454677E-004 7.5714961311373116E-006
-7.2031694214561962E-017 8.6860500453088960E-016
1.1341460255180225E-003 2.4238198932919101E-004
2.4360669979466240E-016 2.7258119127922638E-016
1.9996611328056231 0.0000000000000000
2.5446039125414028E-004 4.2820773824964989E-006
-1.0456613821891448E-016 -9.6743215325285459E-017
9.9676209991647049E-004 -5.6600268038821031E-005
-3.9953280624214323E-016 4.9087316827147771E-016
-3.3048484636195375E-015 6.1817158977131618E-016
-1.2597697359199645E-002 -3.7166324717129856E-006
-2.4210971257827060E-004 5.8389623659711954E-005
-6.3913465705948242E-004 3.9179829909342934E-006
3.1912171522844715E-017 -4.1468769820119327E-016
-1.0278347636459521E-003 -2.0164021732396069E-004
2.0079147210768701E-017 -2.1698316754064389E-016
2.5446039125414028E-004 -4.2820773824964989E-006
1.9997203826784054 0.0000000000000000
2.0825701982571277E-016 -1.0106113319512224E-015
-9.0157781139720526E-004 6.6430604897648815E-005
-1.5542708508403283E-015 7.3957277816789575E-015
1.6429215849485919E-015 -2.7230794497799261E-015
1.0993757087111443E-002 -1.8175938563548891E-004
1.0326496350503192E-016 4.4955675025152728E-017
-3.0529109137315778E-017 5.1270053829163420E-017
5.2770283128296177E-003 4.3814558771467556E-005
-1.4011788447525168E-016 3.8996096313994688E-015
1.3938950857567451E-003 3.1691001834563675E-005
-1.0456613821891448E-016 9.6743215325285459E-017
2.0825701982571277E-016 1.0106113319512224E-015
1.9846289158358426 0.0000000000000000
1.2925808521756067E-016 -9.0381299816088434E-017
0.10628745779824479 2.8561215655848321E-003
-3.5802199235378813E-002 -7.1330629364437653E-004
3.4839587994840887E-015 3.3949303282829035E-015
-9.6307278676681905E-004 1.5848633659367058E-004
-2.4996723334157079E-003 -1.6878280752610168E-004
5.6973450505311861E-016 -2.8852719875277487E-015
-3.9605494357122413E-003 -1.0844788031769959E-003
-3.8095935113214172E-017 -8.5441874687637690E-017
9.9676209991647049E-004 5.6600268038821031E-005
-9.0157781139720526E-004 -6.6430604897648815E-005
1.2925808521756067E-016 9.0381299816088434E-017
1.9964052195342945 0.0000000000000000
3.2032252891242674E-016 -2.9950021822287419E-015
6.2322662101272650E-015 -7.4116428229773892E-015
4.2179457397055385E-002 2.4076080835713082E-003
1.3280369753342969E-016 6.7709544933739064E-017
8.5464454629174907E-016 1.3690739743679114E-015
-3.2844795244561291E-002 6.0975192388691354E-004
3.8998822316673752E-016 -2.0543099601558386E-014
3.3133485502343657E-002 -1.3695969769169161E-004
-3.9953280624214323E-016 -4.9087316827147771E-016
-1.5542708508403283E-015 -7.3957277816789575E-015
0.10628745779824479 -2.8561215655848321E-003
3.2032252891242674E-016 2.9950021822287419E-015
0.54450304079461742 0.0000000000000000
3.7110545796761472E-004 -2.5771121374435629E-006
4.5832917242384341E-016 -1.5226270793952235E-016
3.2813043510687810E-015 1.7418197254477417E-015
4.0612456352043776E-015 6.4796184564793087E-015
2.0106318907807723E-002 -2.3360980952464097E-004
3.7160602329464569E-015 2.7656926837772642E-014
7.9348722336522293E-002 2.2303067036207741E-004
-3.3048484636195375E-015 -6.1817158977131618E-016
1.6429215849485919E-015 2.7230794497799261E-015
-3.5802199235378813E-002 7.1330629364437653E-004
6.2322662101272650E-015 7.4116428229773892E-015
3.7110545796761472E-004 2.5771121374435629E-006
0.32134089194039417 0.0000000000000000
-1.5468343862666688E-015 3.0548585852135864E-017
1.2821097042473820E-002 -2.8686564596901464E-003
3.3255294771216873E-002 3.4591400749335414E-004
-8.1847996124173807E-015 4.0635992692925845E-014
5.2836471927730132E-002 1.1275553097512523E-002
-1.0894729641781911E-014 1.6948278015743026E-014
-1.2597697359199645E-002 3.7166324717129856E-006
1.0993757087111443E-002 1.8175938563548891E-004
3.4839587994840887E-015 -3.3949303282829035E-015
4.2179457397055385E-002 -2.4076080835713082E-003
4.5832917242384341E-016 1.5226270793952235E-016
-1.5468343862666688E-015 -3.0548585852135864E-017
0.32014216847529692 0.0000000000000000
然后重现我的情况的程序是：
- 将 fortran 代码编译成主要可执行文件
- 运行
$ ./main > f 对角化矩阵
运行
cat f | grep '!!!1' | python3 ./fortran2num.py 25 12 T 6
cat f | grep '!!!2' | python3 ./fortran2num.py 25 1 F 6
cat f | grep '!!!3' | python3 ./fortran2num.py 25 12 T 6
并将输出数组复制粘贴到 test_mat.py 脚本中
- 运行
$ python3 ./test_mat.py 9 15
这应该输出类似的东西（不完全是因为它来自稍微修改过的 fortran 代码，但它真的很相似）：
$ ./test_mat.py 9 15
== Innacurate (module > 1.0000000000000005e-09) commponants for zheev
0 5 (-1.8814689806134766e-09-4.2103970310236605e-10j)
1 4 (-3.309330809386779e-09+3.460515207720505e-11j)
3 2 (1.6517163192730777e-09-3.476553767089996e-10j)
3 4 (-1.7113597206019904e-09+3.661825842549149e-10j)
3 7 (4.525875452188188e-09-9.69861357891233e-10j)
4 8 (-6.080326242496881e-10-1.2655512193759932e-09j)
5 1 (3.867461240147918e-09+8.605002561196118e-13j)
5 2 (-1.259992332691695e-09-1.5106162848388394e-12j)
6 0 (-2.4973547674278507e-09+4.124683233852622e-11j)
7 6 (-4.174611545337825e-09-1.5728970348405369e-09j)
7 8 (-1.3306069579499535e-09-2.9363069753429727e-09j)
8 7 (2.440851298857183e-09+1.4367365630540974e-10j)
9 9 (2.9735734430325067e-09-1.9784610061357455e-09j)
10 10 (-2.305021784554384e-09+2.927970713106909e-10j)
11 11 (3.5043191126682416e-09-3.879731280166665e-13j)
== Innacurate (module > 1.0000000000000009e-15) commponants for eigh
9 9 (-1.0547118733938987e-15+7.216449660063518e-16j)
10 10 (-1.1102230246251565e-15+1.1102230246251565e-16j)
== Error on eigen values
1.9999509276358423 [-2.36415776e-09]
1.99994339628587 [-3.71413011e-09]
1.9999370801794079 [1.79407822e-10]
1.9999288599859764 [-1.40236711e-11]
1.9999260278258282 [-2.17417173e-09]
1.9999176681452422 [-1.8547579e-09]
1.994775382976095 [2.97609493e-09]
1.992787044413837 [4.41383685e-09]
1.9899199212878338 [1.28783384e-09]
0.5352109564222124 [-3.57778762e-09]
0.31654720232788414 [2.32788416e-09]
0.31639848351397176 [3.51397178e-09]
如果是这样，那么 lapack 远不如 python 精确。
我不知道问题出在哪里，因为我尝试了 MKL 和 ATLAS 的 zheev 和 zheevd。
我真是个白痴，我没有注意到我将 fortran 输出转换为 numpy 数组的方式我从原始数据中截断了很多小数点！事实上，我在 numpy.array 上使用 repr 来生成可复制粘贴的 numpy 数组并用它来引导我的测试，但我忘记提高 numpy 数组的打印精度。这就是明显不准确的原因。
解决方案就是将 np.set_printoptions(precision=17) 放入我的转换脚本中。
真是惭愧啊，看来真实的结果和numpy的一样好。
我在 Fortran 科学代码中使用 LAPACK zheev 例程来计算不太大（可能永远不会超过大小 1000）的矩阵的特征值和向量。
由于这一步发生在计算的开始，我必须达到很高的准确性以避免重要的错误传播。问题是在我的测试用例中（仅使用 12x12 矩阵）计算的精度仅为 1e-9 左右，这根本不够。
我与 numpy.linalg.eigh 进行了比较，结果好得离谱，我想知道如何使用 zheev 实现如此精确。
这里是一个测试 Python 代码，其中一些来自 Fortran 代码的输出作为数据：
from numpy import array, linalg, diagflat
from sys import argv
prec = 10
if len(argv) > 1:
prec = int(argv[1])
tolZHEEV = 0.1**prec
tolEIGH = tolZHEEV
if len(argv) > 2:
prec2 = int(argv[2])
tolEIGH = 0.1**prec2
F = array([[ 1.99964605e+00+0.00000000e+00j, -6.77332886e-04-1.58965282e-04j,
2.50377192e-16-7.43492040e-16j, -1.03545717e-03-4.75346808e-04j,
-1.33562677e-16-7.96956118e-17j, 2.69297453e-04+6.01705394e-05j,
-2.42109713e-04-5.83896237e-05j, 1.03264964e-16-4.49556750e-17j,
-9.63072787e-04-1.58486337e-04j, 1.32803698e-16-6.77095449e-17j,
3.28130435e-15-1.74181973e-15j, 1.28210970e-02+2.86865646e-03j],
[-6.77332886e-04+1.58965282e-04j, 1.99814145e+00+0.00000000e+00j,
3.11074872e-16-2.03650443e-15j, -2.86327278e-03-5.79965447e-04j,
4.87808786e-19-4.06383000e-17j, 7.07826164e-04-7.57149613e-06j,
-6.39134657e-04-3.91798299e-06j, -3.05291091e-17-5.12700538e-17j,
-2.49967233e-03+1.68782808e-04j, 8.54644546e-16-1.36907397e-15j,
4.06124564e-15-6.47961846e-15j, 3.32552948e-02-3.45914007e-04j],
[ 2.50377192e-16+7.43492040e-16j, 3.11074872e-16+2.03650443e-15j,
1.99801950e+00+0.00000000e+00j, -1.38347852e-16+1.96637141e-15j,
-1.50312544e-03-2.16900444e-05j, -7.20316942e-17-8.68605005e-16j,
3.19121715e-17+4.14687698e-16j, 5.27702831e-03-4.38145588e-05j,
5.69734505e-16+2.88527199e-15j, -3.28447952e-02-6.09751924e-04j,
2.01063189e-02+2.33609810e-04j, -8.18479961e-15-4.06359927e-14j],
[-1.03545717e-03+4.75346808e-04j, -2.86327278e-03+5.79965447e-04j,
-1.38347852e-16-1.96637141e-15j, 1.99514423e+00+0.00000000e+00j,
8.58581114e-17+1.24995999e-15j, 1.13414603e-03-2.42381989e-04j,
-1.02783476e-03+2.01640217e-04j, -1.40117884e-16-3.89960963e-15j,
-3.96054944e-03+1.08447880e-03j, 3.89988223e-16+2.05430996e-14j,
3.71606023e-15-2.76569268e-14j, 5.28364719e-02-1.12755531e-02j],
[-1.33562677e-16+7.96956118e-17j, 4.87808786e-19+4.06383000e-17j,
-1.50312544e-03+2.16900444e-05j, 8.58581114e-17-1.24995999e-15j,
1.98788997e+00+0.00000000e+00j, 2.43606700e-16-2.72581191e-16j,
2.00791472e-17+2.16983168e-16j, 1.39389509e-03-3.16910018e-05j,
-3.80959351e-17+8.54418747e-17j, 3.31334855e-02+1.36959698e-04j,
7.93487223e-02-2.23030670e-04j, -1.08947296e-14-1.69482780e-14j],
[ 2.69297453e-04-6.01705394e-05j, 7.07826164e-04+7.57149613e-06j,
-7.20316942e-17+8.68605005e-16j, 1.13414603e-03+2.42381989e-04j,
2.43606700e-16+2.72581191e-16j, 1.99966113e+00+0.00000000e+00j,
2.54460391e-04+4.28207738e-06j, -1.04566138e-16-9.67432153e-17j,
9.96762100e-04-5.66002680e-05j, -3.99532806e-16+4.90873168e-16j,
-3.30484846e-15+6.18171590e-16j, -1.25976974e-02-3.71663247e-06j],
[-2.42109713e-04+5.83896237e-05j, -6.39134657e-04+3.91798299e-06j,
3.19121715e-17-4.14687698e-16j, -1.02783476e-03-2.01640217e-04j,
2.00791472e-17-2.16983168e-16j, 2.54460391e-04-4.28207738e-06j,
1.99972038e+00+0.00000000e+00j, 2.08257020e-16-1.01061133e-15j,
-9.01577811e-04+6.64306049e-05j, -1.55427085e-15+7.39572778e-15j,
1.64292158e-15-2.72307945e-15j, 1.09937571e-02-1.81759386e-04j],
[ 1.03264964e-16+4.49556750e-17j, -3.05291091e-17+5.12700538e-17j,
5.27702831e-03+4.38145588e-05j, -1.40117884e-16+3.89960963e-15j,
1.39389509e-03+3.16910018e-05j, -1.04566138e-16+9.67432153e-17j,
2.08257020e-16+1.01061133e-15j, 1.98462892e+00+0.00000000e+00j,
1.29258085e-16-9.03812998e-17j, 1.06287458e-01+2.85612157e-03j,
-3.58021992e-02-7.13306294e-04j, 3.48395880e-15+3.39493033e-15j],
[-9.63072787e-04+1.58486337e-04j, -2.49967233e-03-1.68782808e-04j,
5.69734505e-16-2.88527199e-15j, -3.96054944e-03-1.08447880e-03j,
-3.80959351e-17-8.54418747e-17j, 9.96762100e-04+5.66002680e-05j,
-9.01577811e-04-6.64306049e-05j, 1.29258085e-16+9.03812998e-17j,
1.99640522e+00+0.00000000e+00j, 3.20322529e-16-2.99500218e-15j,
6.23226621e-15-7.41164282e-15j, 4.21794574e-02+2.40760808e-03j],
[ 1.32803698e-16+6.77095449e-17j, 8.54644546e-16+1.36907397e-15j,
-3.28447952e-02+6.09751924e-04j, 3.89988223e-16-2.05430996e-14j,
3.31334855e-02-1.36959698e-04j, -3.99532806e-16-4.90873168e-16j,
-1.55427085e-15-7.39572778e-15j, 1.06287458e-01-2.85612157e-03j,
3.20322529e-16+2.99500218e-15j, 5.44503041e-01+0.00000000e+00j,
3.71105458e-04-2.57711214e-06j, 4.58329172e-16-1.52262708e-16j],
[ 3.28130435e-15+1.74181973e-15j, 4.06124564e-15+6.47961846e-15j,
2.01063189e-02-2.33609810e-04j, 3.71606023e-15+2.76569268e-14j,
7.93487223e-02+2.23030670e-04j, -3.30484846e-15-6.18171590e-16j,
1.64292158e-15+2.72307945e-15j, -3.58021992e-02+7.13306294e-04j,
6.23226621e-15+7.41164282e-15j, 3.71105458e-04+2.57711214e-06j,
3.21340892e-01+0.00000000e+00j, -1.54683439e-15+3.05485859e-17j],
[ 1.28210970e-02-2.86865646e-03j, 3.32552948e-02+3.45914007e-04j,
-8.18479961e-15+4.06359927e-14j, 5.28364719e-02+1.12755531e-02j,
-1.08947296e-14+1.69482780e-14j, -1.25976974e-02+3.71663247e-06j,
1.09937571e-02+1.81759386e-04j, 3.48395880e-15-3.39493033e-15j,
4.21794574e-02-2.40760808e-03j, 4.58329172e-16+1.52262708e-16j,
-1.54683439e-15-3.05485859e-17j, 3.20142168e-01+0.00000000e+00j]])
d =array([[1.99995093],
[1.9999434 ],
[1.99993708],
[1.99992886],
[1.99992603],
[1.99991767],
[1.99477538],
[1.99278704],
[1.98991992],
[0.53521096],
[0.3165472 ],
[0.31639848]])
D = diagflat(d)
P = array([[ 8.50795081e-04+1.90361152e-04j, 6.41403660e-03+1.43510867e-03j,
2.93305520e-02+6.56256458e-03j, -2.73319405e-12-4.62616407e-12j,
2.36860608e-02+5.29963787e-03j, 9.63512391e-01+2.15581088e-01j,
-1.18829015e-14+6.38499316e-15j, 1.49823188e-01+3.35221904e-02j,
4.63016072e-15-4.09382827e-15j, 5.35804136e-17-2.23033645e-17j,
-7.22943096e-14-1.78616099e-14j, -7.65412277e-03-1.71257176e-03j],
[ 5.96040942e-03-6.19988216e-05j, 4.13411608e-02-4.30021346e-04j,
1.97753784e-01-2.05698985e-03j, 1.61937250e-11+7.38548825e-12j,
8.90791562e-01-9.26581108e-03j, -8.99773456e-02+9.35923868e-04j,
-6.78075986e-14+2.27240115e-14j, 3.96283860e-01-4.12205452e-03j,
3.62263322e-14+2.70989836e-14j, 2.80882500e-17-1.14706336e-15j,
-1.89947123e-13-2.17176505e-15j, -1.98547508e-02+2.06524599e-04j],
[ 1.43774119e-12-2.04841175e-13j, -3.84634090e-12+1.02008935e-12j,
2.14225297e-12-1.58626154e-12j, 5.76965748e-01+7.48380143e-01j,
-1.64985984e-11-8.05494677e-12j, 5.49834259e-12+8.08747081e-13j,
1.64047208e-01+6.01119828e-02j, -8.41828091e-14-4.52343007e-13j,
1.15974306e-01+2.49814055e-01j, -1.89718719e-02+1.21355877e-02j,
1.21027166e-02-1.40334323e-03j, -1.12121788e-13+3.78501430e-14j],
[ 1.49140349e-02-3.18272560e-03j, 1.07045322e-01-2.28439783e-02j,
5.83924798e-01-1.24612315e-01j, -7.68362748e-12-3.20504318e-12j,
-4.29872936e-01+9.17369184e-02j, -1.07323514e-01+2.29033458e-02j,
2.52998000e-14-3.37978544e-14j, 6.37564522e-01-1.36059286e-01j,
2.28120526e-13-8.82579171e-14j, 8.15806584e-15+1.14028011e-14j,
-3.05058814e-13+4.83321942e-14j, -3.15470869e-02+6.73229760e-03j],
[-1.50646374e-13+4.30474354e-14j, 2.88653601e-13-7.82348111e-14j,
-1.56416337e-13+6.81259624e-14j, -4.56052348e-02-5.74215381e-02j,
1.24973393e-12+6.35459037e-13j, -4.16040415e-13-4.47866383e-14j,
-6.68174568e-01-2.33961052e-01j, -7.77941534e-14-1.12876492e-13j,
3.04089801e-01+6.31021660e-01j, 1.92382317e-02-1.27008073e-02j,
4.69043183e-02-6.12576169e-03j, -4.46806257e-13+6.93146032e-14j],
[-2.00864022e-02-5.92598703e-06j, -9.12820255e-01-2.69304562e-04j,
3.72566080e-01+1.09916213e-04j, -1.06450784e-12-4.32276375e-12j,
3.36150018e-02+9.91725789e-06j, 1.90596074e-02+5.62305535e-06j,
3.49708145e-14+4.34006214e-15j, -1.61246105e-01-4.75715912e-05j,
8.38875458e-15-4.53541474e-15j, -2.47012312e-16+4.74760834e-16j,
7.10584649e-14+6.19010235e-16j, 7.52355819e-03+2.21963585e-06j],
[ 9.62420518e-01-1.59116634e-02j, -1.24950681e-01+2.06580506e-03j,
-1.86710046e-01+3.08687039e-03j, -1.30253808e-12-1.41868926e-12j,
-2.69408788e-02+4.45412570e-04j, -1.67620621e-02+2.77126563e-04j,
-1.12422533e-16-5.40641968e-14j, 1.48156214e-01-2.44946130e-03j,
5.33921595e-14+3.20699063e-14j, 2.10253369e-15+4.74023210e-15j,
-6.27766220e-14-7.29612349e-16j, -6.56715101e-03+1.08574469e-04j],
[ 5.05414777e-13+9.40307821e-15j, -1.32687279e-12+3.49123357e-13j,
7.54773974e-13-3.92282159e-13j, 1.92570500e-01+2.54118137e-01j,
-5.58410396e-12-2.86684025e-12j, 1.83884514e-12+2.64297906e-13j,
-6.36994666e-01-2.39431942e-01j, 1.93041775e-14-4.12911291e-14j,
-2.70995935e-01-5.96659039e-01j, 6.18161118e-02-3.88219877e-02j,
-2.18180305e-02+2.34644894e-03j, 2.08901179e-13-2.47813917e-14j],
[-2.69415450e-01-1.53782636e-02j, -3.70083377e-01-2.11244000e-02j,
-6.54530261e-01-3.73606595e-02j, -1.63145928e-12-8.78784663e-14j,
-1.03526300e-01-5.90929264e-03j, -6.61072940e-02-3.77341164e-03j,
-4.04464076e-14+5.89878883e-14j, 5.85724599e-01+3.34332248e-02j,
3.15506629e-14+2.84647065e-14j, -1.09481211e-15-1.77104489e-15j,
-2.40425891e-13-1.87628386e-14j, -2.51994911e-02-1.43838973e-03j],
[ 2.54752203e-16+2.12657752e-16j, -3.21054014e-15+3.13509411e-16j,
1.50876469e-15-4.65388714e-16j, 1.84339353e-04+2.30141854e-04j,
-7.60919954e-15-2.39290177e-15j, 1.07829490e-15+1.72556420e-16j,
-6.61876184e-02-2.28688765e-02j, 9.75386453e-17-3.60186046e-15j,
-1.68109710e-02-3.45191353e-02j, -8.29576318e-01+5.52611651e-01j,
6.68078399e-03-9.00620322e-04j, -6.62742104e-14+9.74645368e-15j],
[-4.83726841e-16-3.76055196e-17j, -3.28302552e-15+1.41917600e-15j,
2.80242692e-15-7.23389671e-16j, 6.51697903e-04+8.25305197e-04j,
-1.99649588e-14-9.17555906e-15j, 7.51565389e-15+9.10141876e-16j,
-1.59564626e-02-5.63754717e-03j, 2.04075321e-15+5.12238396e-15j,
2.18750663e-02+4.57214672e-02j, -6.28420329e-03+4.12342677e-03j,
-9.90493516e-01+1.26529179e-01j, 9.57441797e-12-1.22296432e-12j],
[ 2.79384410e-04+0.00000000e+00j, 1.09500742e-03+0.00000000e+00j,
2.84867177e-03+0.00000000e+00j, 9.54279829e-15+0.00000000e+00j,
6.53433434e-04+0.00000000e+00j, 4.93180792e-04+0.00000000e+00j,
1.04504620e-15+0.00000000e+00j, 4.71495396e-02+0.00000000e+00j,
-7.15402040e-15+0.00000000e+00j, -1.97732525e-15+0.00000000e+00j,
9.65544745e-12+0.00000000e+00j, 9.98882805e-01+0.00000000e+00j]])
fp = F.dot(P) # FP = PD
dp = P.dot(D)
dd, PP = linalg.eigh(-F, 'U')
dd = -dd
DD = diagflat(dd)
ffp = F.dot(PP) # FP = PD
ddp = PP.dot(DD)
zer = fp - dp
zernp = ffp - ddp
print("\n== Innacurate (module > " + str(tolZHEEV) + ") commponants for zheev")
for n in range(12):
for p in range(12):
if abs(zer[n,p]) > tolZHEEV:
print(n, p, zer[n,p])
print("\n== Innacurate (module > " + str(tolEIGH) + ") commponants for eigh")
for n in range(12):
for p in range(12):
if abs(zernp[n,p]) > tolEIGH:
print(n, p, zernp[n,p])
print("\n== Error on eigen values")
for i in range(dd.shape[0]):
print(dd[i], dd[i] - d[i])
然后按照两个fortran文件重现数据： 主.F90
program main
implicit none
integer, parameter :: nband = 12
complex*16 :: occ_nd_cpx(nband, nband)
real*8 :: occ_diag(nband), a, b
integer :: i, j
do i=1,nband
do j=1,nband
read (111,'(2D25.16)') a, b
occ_nd_cpx(i,j) = cmplx(a, b, kind=8)
end do
end do
call diag_occ(occ_nd_cpx, nband, occ_diag)
end program main
diag_occ.F90
subroutine diag_occ(occ_nd_cpx, nband, occ_diag)
implicit none
integer, parameter :: dp=kind(1.0d0), dpc = dp
integer,intent(in) :: nband
complex(kind=dpc), intent(inout) :: occ_nd_cpx(nband, nband)
real(kind=dp), intent(out) :: occ_diag(nband)
integer :: info, lwork, n
character(len=500) :: message
real(kind=dp) :: rwork(3*nband-1)
complex(kind=dpc), allocatable :: work(:)
do n=1,nband
write (*,*) '!!!1', occ_nd_cpx(n, :)
end do
!! Initialisation
rwork = 0.0_dp
info = 0
!! use the opposite to have zheev orders the eigenvalues descending (once the
!! opposite have been taken again)
occ_nd_cpx = -occ_nd_cpx
!! Get diagonal occupations and associeted base
! Compute the optimal working array size
allocate(work(1))
work = (0.0_dp, 0.0_dp)
call zheev('V', 'U', nband, occ_nd_cpx, nband, occ_diag, work, -1, rwork, info)
lwork = work(1)
write (*,*) 'lwork, work(1)', lwork, work(1)
deallocate(work)
! Compute the eigenvalues (occ_diag) and vectors
allocate(work(lwork))
work = (0.0_dp, 0.0_dp)
call zheev('V', 'U', nband, occ_nd_cpx, nband, occ_diag, work, lwork, rwork, info)
!! obtain the true eigen values of occupation matrix in descending order
occ_diag = -occ_diag
do n=1,nband
write (*,*) '!!!2', occ_diag(n)
end do
do n=1,nband
write (*,*) '!!!3', occ_nd_cpx(n, :)
end do
info = 0
deallocate(work)
end subroutine diag_occ
我用gfortran -ffree-form -llapack main.F90 diag_occ.F90 -o main编译
LAPCK 库应该来自 CentOS 存储库 (atlas)
我使用以下脚本将主可执行文件的输出转换为 numpy 数组： fortran2num.py
import sys
import numpy as np
## Format
FSIZE = 25
CPX = True
N = 12
IGNORE = 0
if len(sys.argv) > 1 and sys.argv[1] == '-h':
print(sys.argv[0] + ' [FSIZE [N [CMPX]]]')
exit()
if len(sys.argv) > 1:
FSIZE = int(sys.argv[1])
if len(sys.argv) > 2:
N = int(sys.argv[2])
if len(sys.argv) > 3:
CPX = sys.argv[3] in ['true', 'True', '1', 'T', 't']
if len(sys.argv) > 4:
IGNORE = int(sys.argv[4])
src = sys.stdin.read()
lst = [[]]
c = IGNORE
k = 0
while c < len(src):
k += 1
if CPX:
step = 2*FSIZE+3
a,b = src[c+1:c+step-1].split(',')
lst[-1].append(float(a)+1j*float(b))
else:
step = FSIZE
a = src[c:c+step]
lst[-1].append(float(a))
c += step + 1
if k >= N:
k = 0
c += IGNORE
lst.append([])
if len(lst[-1]) == 0:
lst.pop()
print(repr(np.array(lst)))
我使用这些命令创建了 F、d 和 P：
cat f | grep '!!!1' | ./fortran2num.py 25 12 T 6
cat f | grep '!!!2' | ./fortran2num.py 25 1 F 6
cat f | grep '!!!3' | ./fortran2num.py 25 12 T 6
编辑： 上次编辑这个我很匆忙，忘记了一些数据。这是主要可执行文件的输入 fortran.111 文件：
1.9996460536245420 0.0000000000000000
-6.7733288642473765E-004 -1.5896528158261059E-004
2.5037719177319647E-016 -7.4349203998305547E-016
-1.0354571719663188E-003 -4.7534680782504999E-004
-1.3356267732320494E-016 -7.9695611788720466E-017
2.6929745274964923E-004 6.0170539417274264E-005
-2.4210971257827060E-004 -5.8389623659711954E-005
1.0326496350503192E-016 -4.4955675025152728E-017
-9.6307278676681905E-004 -1.5848633659367058E-004
1.3280369753342969E-016 -6.7709544933739064E-017
3.2813043510687810E-015 -1.7418197254477417E-015
1.2821097042473820E-002 2.8686564596901464E-003
-6.7733288642473765E-004 1.5896528158261059E-004
1.9981414542927907 0.0000000000000000
3.1107487236320250E-016 -2.0365044300513661E-015
-2.8632727795318471E-003 -5.7996544650846344E-004
4.8780878630335091E-019 -4.0638299964232169E-017
7.0782616404454677E-004 -7.5714961311373116E-006
-6.3913465705948242E-004 -3.9179829909342934E-006
-3.0529109137315778E-017 -5.1270053829163420E-017
-2.4996723334157079E-003 1.6878280752610168E-004
8.5464454629174907E-016 -1.3690739743679114E-015
4.0612456352043776E-015 -6.4796184564793087E-015
3.3255294771216873E-002 -3.4591400749335414E-004
2.5037719177319647E-016 7.4349203998305547E-016
3.1107487236320250E-016 2.0365044300513661E-015
1.9980194987385831 0.0000000000000000
-1.3834785212208944E-016 1.9663714097332382E-015
-1.5031254367981852E-003 -2.1690044364606380E-005
-7.2031694214561962E-017 -8.6860500453088960E-016
3.1912171522844715E-017 4.1468769820119327E-016
5.2770283128296177E-003 -4.3814558771467556E-005
5.6973450505311861E-016 2.8852719875277487E-015
-3.2844795244561291E-002 -6.0975192388691354E-004
2.0106318907807723E-002 2.3360980952464097E-004
-8.1847996124173807E-015 -4.0635992692925845E-014
-1.0354571719663188E-003 4.7534680782504999E-004
-2.8632727795318471E-003 5.7996544650846344E-004
-1.3834785212208944E-016 -1.9663714097332382E-015
1.9951442266000110 0.0000000000000000
8.5858111362132558E-017 1.2499599851488651E-015
1.1341460255180225E-003 -2.4238198932919101E-004
-1.0278347636459521E-003 2.0164021732396069E-004
-1.4011788447525168E-016 -3.8996096313994688E-015
-3.9605494357122413E-003 1.0844788031769959E-003
3.8998822316673752E-016 2.0543099601558386E-014
3.7160602329464569E-015 -2.7656926837772642E-014
5.2836471927730132E-002 -1.1275553097512523E-002
-1.3356267732320494E-016 7.9695611788720466E-017
4.8780878630335091E-019 4.0638299964232169E-017
-1.5031254367981852E-003 2.1690044364606380E-005
8.5858111362132558E-017 -1.2499599851488651E-015
1.9878899727642418 0.0000000000000000
2.4360669979466240E-016 -2.7258119127922638E-016
2.0079147210768701E-017 2.1698316754064389E-016
1.3938950857567451E-003 -3.1691001834563675E-005
-3.8095935113214172E-017 8.5441874687637690E-017
3.3133485502343657E-002 1.3695969769169161E-004
7.9348722336522293E-002 -2.2303067036207741E-004
-1.0894729641781911E-014 -1.6948278015743026E-014
2.6929745274964923E-004 -6.0170539417274264E-005
7.0782616404454677E-004 7.5714961311373116E-006
-7.2031694214561962E-017 8.6860500453088960E-016
1.1341460255180225E-003 2.4238198932919101E-004
2.4360669979466240E-016 2.7258119127922638E-016
1.9996611328056231 0.0000000000000000
2.5446039125414028E-004 4.2820773824964989E-006
-1.0456613821891448E-016 -9.6743215325285459E-017
9.9676209991647049E-004 -5.6600268038821031E-005
-3.9953280624214323E-016 4.9087316827147771E-016
-3.3048484636195375E-015 6.1817158977131618E-016
-1.2597697359199645E-002 -3.7166324717129856E-006
-2.4210971257827060E-004 5.8389623659711954E-005
-6.3913465705948242E-004 3.9179829909342934E-006
3.1912171522844715E-017 -4.1468769820119327E-016
-1.0278347636459521E-003 -2.0164021732396069E-004
2.0079147210768701E-017 -2.1698316754064389E-016
2.5446039125414028E-004 -4.2820773824964989E-006
1.9997203826784054 0.0000000000000000
2.0825701982571277E-016 -1.0106113319512224E-015
-9.0157781139720526E-004 6.6430604897648815E-005
-1.5542708508403283E-015 7.3957277816789575E-015
1.6429215849485919E-015 -2.7230794497799261E-015
1.0993757087111443E-002 -1.8175938563548891E-004
1.0326496350503192E-016 4.4955675025152728E-017
-3.0529109137315778E-017 5.1270053829163420E-017
5.2770283128296177E-003 4.3814558771467556E-005
-1.4011788447525168E-016 3.8996096313994688E-015
1.3938950857567451E-003 3.1691001834563675E-005
-1.0456613821891448E-016 9.6743215325285459E-017
2.0825701982571277E-016 1.0106113319512224E-015
1.9846289158358426 0.0000000000000000
1.2925808521756067E-016 -9.0381299816088434E-017
0.10628745779824479 2.8561215655848321E-003
-3.5802199235378813E-002 -7.1330629364437653E-004
3.4839587994840887E-015 3.3949303282829035E-015
-9.6307278676681905E-004 1.5848633659367058E-004
-2.4996723334157079E-003 -1.6878280752610168E-004
5.6973450505311861E-016 -2.8852719875277487E-015
-3.9605494357122413E-003 -1.0844788031769959E-003
-3.8095935113214172E-017 -8.5441874687637690E-017
9.9676209991647049E-004 5.6600268038821031E-005
-9.0157781139720526E-004 -6.6430604897648815E-005
1.2925808521756067E-016 9.0381299816088434E-017
1.9964052195342945 0.0000000000000000
3.2032252891242674E-016 -2.9950021822287419E-015
6.2322662101272650E-015 -7.4116428229773892E-015
4.2179457397055385E-002 2.4076080835713082E-003
1.3280369753342969E-016 6.7709544933739064E-017
8.5464454629174907E-016 1.3690739743679114E-015
-3.2844795244561291E-002 6.0975192388691354E-004
3.8998822316673752E-016 -2.0543099601558386E-014
3.3133485502343657E-002 -1.3695969769169161E-004
-3.9953280624214323E-016 -4.9087316827147771E-016
-1.5542708508403283E-015 -7.3957277816789575E-015
0.10628745779824479 -2.8561215655848321E-003
3.2032252891242674E-016 2.9950021822287419E-015
0.54450304079461742 0.0000000000000000
3.7110545796761472E-004 -2.5771121374435629E-006
4.5832917242384341E-016 -1.5226270793952235E-016
3.2813043510687810E-015 1.7418197254477417E-015
4.0612456352043776E-015 6.4796184564793087E-015
2.0106318907807723E-002 -2.3360980952464097E-004
3.7160602329464569E-015 2.7656926837772642E-014
7.9348722336522293E-002 2.2303067036207741E-004
-3.3048484636195375E-015 -6.1817158977131618E-016
1.6429215849485919E-015 2.7230794497799261E-015
-3.5802199235378813E-002 7.1330629364437653E-004
6.2322662101272650E-015 7.4116428229773892E-015
3.7110545796761472E-004 2.5771121374435629E-006
0.32134089194039417 0.0000000000000000
-1.5468343862666688E-015 3.0548585852135864E-017
1.2821097042473820E-002 -2.8686564596901464E-003
3.3255294771216873E-002 3.4591400749335414E-004
-8.1847996124173807E-015 4.0635992692925845E-014
5.2836471927730132E-002 1.1275553097512523E-002
-1.0894729641781911E-014 1.6948278015743026E-014
-1.2597697359199645E-002 3.7166324717129856E-006
1.0993757087111443E-002 1.8175938563548891E-004
3.4839587994840887E-015 -3.3949303282829035E-015
4.2179457397055385E-002 -2.4076080835713082E-003
4.5832917242384341E-016 1.5226270793952235E-016
-1.5468343862666688E-015 -3.0548585852135864E-017
0.32014216847529692 0.0000000000000000
然后重现我的情况的程序是：
- 将 fortran 代码编译成主要可执行文件
- 运行
$ ./main > f对角化矩阵 运行
cat f | grep '!!!1' | python3 ./fortran2num.py 25 12 T 6 cat f | grep '!!!2' | python3 ./fortran2num.py 25 1 F 6 cat f | grep '!!!3' | python3 ./fortran2num.py 25 12 T 6
并将输出数组复制粘贴到 test_mat.py 脚本中
- 运行
$ python3 ./test_mat.py 9 15
这应该输出类似的东西（不完全是因为它来自稍微修改过的 fortran 代码，但它真的很相似）：
$ ./test_mat.py 9 15
== Innacurate (module > 1.0000000000000005e-09) commponants for zheev
0 5 (-1.8814689806134766e-09-4.2103970310236605e-10j)
1 4 (-3.309330809386779e-09+3.460515207720505e-11j)
3 2 (1.6517163192730777e-09-3.476553767089996e-10j)
3 4 (-1.7113597206019904e-09+3.661825842549149e-10j)
3 7 (4.525875452188188e-09-9.69861357891233e-10j)
4 8 (-6.080326242496881e-10-1.2655512193759932e-09j)
5 1 (3.867461240147918e-09+8.605002561196118e-13j)
5 2 (-1.259992332691695e-09-1.5106162848388394e-12j)
6 0 (-2.4973547674278507e-09+4.124683233852622e-11j)
7 6 (-4.174611545337825e-09-1.5728970348405369e-09j)
7 8 (-1.3306069579499535e-09-2.9363069753429727e-09j)
8 7 (2.440851298857183e-09+1.4367365630540974e-10j)
9 9 (2.9735734430325067e-09-1.9784610061357455e-09j)
10 10 (-2.305021784554384e-09+2.927970713106909e-10j)
11 11 (3.5043191126682416e-09-3.879731280166665e-13j)
== Innacurate (module > 1.0000000000000009e-15) commponants for eigh
9 9 (-1.0547118733938987e-15+7.216449660063518e-16j)
10 10 (-1.1102230246251565e-15+1.1102230246251565e-16j)
== Error on eigen values
1.9999509276358423 [-2.36415776e-09]
1.99994339628587 [-3.71413011e-09]
1.9999370801794079 [1.79407822e-10]
1.9999288599859764 [-1.40236711e-11]
1.9999260278258282 [-2.17417173e-09]
1.9999176681452422 [-1.8547579e-09]
1.994775382976095 [2.97609493e-09]
1.992787044413837 [4.41383685e-09]
1.9899199212878338 [1.28783384e-09]
0.5352109564222124 [-3.57778762e-09]
0.31654720232788414 [2.32788416e-09]
0.31639848351397176 [3.51397178e-09]
如果是这样，那么 lapack 远不如 python 精确。 我不知道问题出在哪里，因为我尝试了 MKL 和 ATLAS 的 zheev 和 zheevd。
我真是个白痴，我没有注意到我将 fortran 输出转换为 numpy 数组的方式我从原始数据中截断了很多小数点！事实上，我在 numpy.array 上使用 repr 来生成可复制粘贴的 numpy 数组并用它来引导我的测试，但我忘记提高 numpy 数组的打印精度。这就是明显不准确的原因。
解决方案就是将 np.set_printoptions(precision=17) 放入我的转换脚本中。
真是惭愧啊，看来真实的结果和numpy的一样好。