Вычисление обратной матрицы Fortran с использованием SGETR (F, I) работает только с одинарной точностью.

Я пытаюсь вычислить обратную квадратную матрицу, но это не работает. Я проверил предыдущие сообщения, но логика такая же, но я до сих пор не нашел, в чем проблема. Я также делюсь результатом Matlab, например, матрицей.

program test

Implicit none

real,allocatable,dimension(:,:)         :: A       
real,allocatable,dimension(:)           :: WORK
integer ,allocatable,dimension(:)       :: ipiv
integer                                 :: n,info,M
external     SGETRF,SGETRI
M=8
allocate(A(M,M),WORK(M),IPIV(M))

A(1,:)=(/3.74E-4, 0.0, 0.0, 4.98E-5, 0.0, 0.0, 0.0, 0.0/)
A(2,:)=(/0.0 , 3.74E-4, 0.0, 0.0, 4.98E-5 ,0.0 ,0.0 ,0.0 /)
A(3,:)=(/0.0 , 0.0 ,3.74E-4, 0.0 ,0.0, 4.98E-5, 0.0 ,0.0/)
A(4,:)=(/4.98E-5 ,0.0 ,0.0 ,6.64e-6, 0.0 ,0.0, 0.0, 0.0 /)
A(5,:)=(/0.0 , 4.98E-5, 0.0, 0.0 ,6.64E-6 ,0.0 ,0.0 ,0.0 /)
A(6,:)=(/0.0, 0.0, 4.98E-5, 0.0 ,0.0, 6.64E-6, 0.0 ,0.0 /)
A(7,:)=(/0.0, 0.0 ,0.0, 0.0 ,0.0 ,0.0 ,1.49E-11, 0.0 /)
A(8,:)=(/0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0, 0.0 ,1.49E-11 /)


call SGETRF(M,M,A,M,IPIV,info)
if(info .eq. 0) then
   Print *,'succeded'
else
   Print *,'failed'
end if

call SGETRI(M,A,M,IPIV,WORK,M,info)
if(info .eq. 0) then
  Print *,'succeded'
else
  Print *,'failed'
end if
Print *,A

deallocate(A,IPIV,WORK)

end 

!!!!! Matlab Result
!1.0e+10 *
! 0.0002     0       0   -0.0015       0      0        0   0
!     0      0.0002  0       0       -0.0015  0        0   0
!     0      0    0.0002     0         0     -0.0015   0   0
! -0.0015    0       0     0.0113      0      0        0   0
!     0     -0.0015  0       0       0.0113   0        0   0
!     0      0   -0.0015     0         0    0.0113     0   0
!     0      0       0       0         0      0     6.7114 0
!     0      0       0       0         0      0        0   6.7114

fortran matrix-inverse lapack fortran90 fortran95
person Nobody    schedule 12.06.2018    source источник
comment
Пожалуйста, оставьте вопрос как есть, потому что теперь на него есть принятый ответ. Не переписывайте код с новым вопросом, потому что теперь у вас есть новый вопрос stackoverflow.com/questions/50833567/ для этого. Сохраните это в новом вопросе.   -  person Vladimir F    schedule 13.06.2018

Ответы (1)


arrow_upward
4
arrow_downward

Ваши reals имеют только одинарную точность. Префикс lapack D подразумевает двойную точность. Два исправления:

  1. Измените свои DGs на SGs
  2. Сохраняйте DG и используйте двойную точность
person ptb    schedule 12.06.2018
comment
Я изменил для двойной точности, но результат неверен, но для одинарной точности работает! - person Nobody; 13.06.2018
comment
Основываясь на новом заголовке вопроса, кажется, что вы пытаетесь вызвать подпрограммы lapack с одинарной точностью с данными двойной точности. Это правильно? - person ptb; 13.06.2018
comment
Корпус в плохом состоянии для матрицы. Он также работает с двойной точностью, но результат немного странный. После обратной матрицы я разложу матрицу на разложение Холецкого. Но если я это сделаю, матрица будет сингулярной, а не несингулярной. - person Nobody; 13.06.2018