How to control lapack solver presion

Question

I'm trying to use Lapack ZHESV subroutine in Fortran to solve a linear system, but the precision seems not good.

Here is the code:

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

The outputs are

 --------------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)

The error is quite large compare to some element of B.

In contrary, I tried to solve this in Mathematica, and the error is quite small.

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}

So how can I control the precision of the Lapack solver to achieve the same precision of Mathematica's LinearSolver?

Answer 1

ZHESV computes the solution to a complex system of linear equations A * X = B, where A is an N-by-N Hermitian matrix and X and B are N-by-NRHS matrices.

But, your matrix A is not a Hermitian matrix.

Answer 2

Recall that, after the subroutine ZHESV call, X is no longer the initial value (define initially by B) but the solution to A*X=B. When you compute your error there, you are actually computing matmul(A, solution) - initial when you really want matmul(A, initial) - solution. Fixing this (i.e., swapping X and B in that line), I get

 -------------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)

Pretty good error, if you ask me.