Как контролировать точность решателя lapack

Я пытаюсь использовать подпрограмму Lapack ZHESV в Fortran для решения линейной системы, но точность кажется не очень хорошей.

Вот код:

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}

Итак, как я могу контролировать точность решателя Лапака, чтобы достичь той же точности, что и у Mathematica LinearSolver?


precision fortran lapack
person xslittlegrass    schedule 24.07.2013    source источник

Ответы (2)


arrow_upward
4
arrow_downward

ZHESV вычисляет решение сложной системы линейных уравнений A * X = B, где A — это эрмитова матрица размера N на N, а X и B — это матрицы размера N на N.

Но ваша матрица A не является эрмитовой матрицей.

person Yong Yang    schedule 24.07.2013
comment
Как я скучал по этому! Спасибо большое!! - person xslittlegrass; 25.07.2013

arrow_upward
0
arrow_downward

Напомним, что после вызова подпрограммы 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)

Довольно хорошая ошибка, если вы спросите меня.

person Kyle Kanos    schedule 24.07.2013
comment
Спасибо за ответ. Если я правильно понимаю, что initial=B, X является решением AX=B, то AX-B должно быть маленьким, верно? Почему мы должны использовать matmul(A,initial)-solution для проверки A*B-X? - person xslittlegrass; 24.07.2013
comment
Хм, наверное, мне следует перестать пытаться помочь, когда я лишен сна (новые двойняшки). Я собирался решить A'A*X = A'B--›X = A'B--›0=A'B-X и пропустил часть, где нужно транспонировать (снова) A. Это дает мне 0 по всем направлениям. - person Kyle Kanos; 24.07.2013