Zheev特征值和向量不正确 [英] Inaccurate zheev eigen values and vectors

问题描述

我在Fortran科学代码中使用LAPACK zheev例程来计算不太大(可能永远不会超过1000大小)的矩阵的特征值和向量.

I use the LAPACK zheev routine in a Fortran scientific code to compute both eigenvalues and vectors of a matrix not too big (likely to never exceed size 1000).

由于此步骤发生在计算的开始，所以我必须获得很高的精度，以避免重要的错误传播.问题在于，在我的测试用例中(仅使用12x12矩阵)，计算精度仅为1e-9左右，这根本不够.

Since this step happens at the beginning of the computation I have to achieve a great accuracy to avoid important error propagation. The problem is that the accuracy of the computation is only about 1e-9 in my test case (with a 12x12 matrix only) and that is not enough at all.

我与numpy.linalg.eigh进行了比较，后者给出了可笑的更好结果，我想知道如何使用zheev达到这样的精度.

I compared with numpy.linalg.eigh which give ridiculously better results and I wonder how I can achieve such a precision with zheev.

为说明这一点，它是一个测试Python代码，其中有一些Fortran代码的输出作为数据:

To illustrate here is a test Python code with some output from the Fortran code as datas:

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文件重现数据: main.F90

Then follow the two fortran files to reproduce the data : main.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

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存储库(图集)

I compiled with gfortran -ffree-form -llapack main.F90 diag_occ.F90 -o main and the LAPCK lib should come from CentOS repository (atlas)

我使用以下脚本将主要可执行文件的输出转换为numpy数组: fortran2num.py

I use the following script to convert outpout from the main executable into numpy array : 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:

I used these commands to create F, d and 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文件:

EDIT : I was in a rush last time I edited this and I forgot some data. Here come the fortran.111 file that is the input of the main executable :

 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

然后重现我的情况，程序是:

Then to reproduce my situations the procedure is :

• 将fortran代码编译为主要可执行文件
• 运行$ ./main > f以对角化矩阵

• 运行

• Compile the fortran code into main executable
• run $ ./main > f to diagonalise matrix

• run

 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脚本中

and copy-paste the output arrays into the test_mat.py script

• 运行$ python3 ./test_mat.py 9 15

这应该输出类似的内容(不完全是因为它来自稍作修改的fortran代码，但确实如此):

This should output something like that (not exactly since this come from a slightly modified fortran code but it is realy alike) :

$ ./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.

If so then lapack is far less precise than python. I have no clue what the problem is since I tried both zheev and zheevd from both MKL and ATLAS.

推荐答案

我真是个白痴，我没有注意到我将fortran输出转换为numpy数组的方式会截断原始数据中的许多小数！确实，我在numpy.array上使用过repr来生成可复制复制的numpy数组，并用它来进行测试，但我忘了提高numpy数组的打印精度.这就是明显不正确的原因.

I am such an idiot that I didn't noticed that the way I convert fortran output to numpy array I truncate a lot of decimal from the original data ! Indeed I used repr on numpy.array to produce copy-pastable numpy arrays and use it to lead my test but I forgot to increase the printed precision of numpy arrays. That was the cause of the apparent inaccuracy.

解决方案是将np.set_printoptions(precision=17)放入我的转换脚本中.

The solution was just to put np.set_printoptions(precision=17) in my conversion script.

我真的为这种困惑感到羞耻，看来真正的结果和numpy的结果一样好.

I am truly ashamed of this confusion, it appear that the true results are as good as the ones from numpy.

这篇关于Zheev特征值和向量不正确的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！