0

我正在尝试使用ZHESVFortran 中的 Lapack 子程序来求解线性系统，但精度似乎并不好。

这是代码：

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

4

ZHESV 计算复杂线性方程组 的解A * X = B，其中A是N×N Hermitian 矩阵，X和B是N×NRHS 矩阵。

但是，您的矩阵A不是 Hermitian 矩阵。

于 2013-07-24T16:52:17.587 回答
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 回答