0

Fortran でLapackZHESVサブルーチンを使用して線形システムを解こうとしていますが、精度が良くないようです。

コードは次のとおりです。

program main
implicit none

integer,parameter::N=4
integer::LDA=N,IPIV(N),LDB=N,LWORK=N*N,info,i
complex*16::A(N,N),B(N),work(N,N),X(N)


A=reshape( (/(1.,0.),(0.,0.),(0.,-6.94908E-13),(0.,-6.94908E-13),&
          &(0.,0.),(1.,0.0352595),(0.,-4.51893E-11),(0.,-4.51893E-11),&
          &(0.,-6.94908E-13),(0.,-4.51893E-11),(1.,0.0376938),(0.,0.),&
          &(0.,-6.94908E-13),(0.,-4.51893E-11),(0.,0.),(1.,0.0378932)/),shape(A))

A=TRANSPOSE(A)

B=(/(1.,0.),(0.,0.),(0.,6.94908E-13),(0.,6.94908E-13)/)
X=B

write(*,*) "--------------B----------------"
write(*,99999) B


CALL ZHESV('Upper',N,1,A,LDA,IPIV,X,N,WORK,LWORK,INFO)

write(*,*) "--------------x----------------"
write(*,99999) X
write(*,*) "-------------INFO--------------"
write(*,*) INFO

write(*,*) "-------------error-------------"
write(*,99999) matmul(A,X)-B

99999 FORMAT ((3X,4(' (',E15.8,',',E15.8,')',:)))
end program main

出力は

 --------------B----------------
    ( 0.10000000E+01, 0.00000000E+00) ( 0.00000000E+00, 0.00000000E+00) ( 0.00000000E+00, 0.69490800E-12) ( 0.00000000E+00, 0.69490800E-12)
 --------------x----------------
    ( 0.10000000E+01, 0.00000000E+00) ( 0.00000000E+00, 0.00000000E+00) ( 0.00000000E+00, 0.00000000E+00) ( 0.00000000E+00, 0.00000000E+00)
 -------------INFO--------------
           0
 -------------error-------------
    ( 0.00000000E+00, 0.00000000E+00) ( 0.00000000E+00, 0.00000000E+00) ( 0.00000000E+00,-0.13898160E-11) ( 0.00000000E+00,-0.13898160E-11)

この誤差は、B の一部の要素に比べてかなり大きいです。

逆にこれを Mathematica で解いてみましたが、誤差はかなり小さいです。

A.LinearSolve[A, B] - B

{0. + 0. I, 0. + 3.67342*10^-40 I, 4.13609*10^-34 + 0. I, 4.13609*10^-34 + 0. I}

では、Lapack ソルバーの精度を制御して Mathematica と同じ精度を達成するにはどうすればよいLinearSolverでしょうか?

4

2 に答える 2

0

サブルーチンZHESV呼び出しの後、Xはもはや初期値 (最初に で定義B) ではなく、 の解であることを思い出してくださいA*X=B。そこで誤差を計算すると、matmul(A, solution) - initial本当に必要なときに実際に計算されますmatmul(A, initial) - solution。これを修正すると (つまり、スワップXしB、その行で)、次のようになります。

 -------------error-------------
(    0.00000E+00,    0.00000E+00)(    0.62805E-22,    0.00000E+00)(    0.00000E+00,    0.00000E+00)(    0.00000E+00,    0.00000E+00)

あなたが私に尋ねると、かなり良いエラーです。

于 2013-07-24T02:36:53.473 に答える