title: How to control lapack solver presion
body: <p>I'm trying to use Lapack <code>ZHESV</code> subroutine in Fortran to solve a linear system, but the precision seems not good.</p>
  
  <p>Here is the code:</p>
  
  <pre><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),&amp;
            &amp;(0.,0.),(1.,0.0352595),(0.,-4.51893E-11),(0.,-4.51893E-11),&amp;
            &amp;(0.,-6.94908E-13),(0.,-4.51893E-11),(1.,0.0376938),(0.,0.),&amp;
            &amp;(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
  </code></pre>
  
  <p>The outputs are</p>
  
  <pre><code> --------------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)
  </code></pre>
  
  <p>The error is quite large compare to some element of B.</p>
  
  <p>In contrary, I tried to solve this in Mathematica, and the error is quite small.</p>
  
  <pre><code>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}
  </code></pre>
  
  <p>So how can I control the precision of the Lapack solver to achieve the same precision of Mathematica's <code>LinearSolver</code>?</p>