FORTRAN неправильно работает с матрицами!

[simplermethod]

Провел разными способами расчет обращения матрицы,
Результаты - не совпадают!
Чему верить?

Входная матрица
|1.0|7.|11.|
|7.0|3.|9.0|
|11.|9.|5.0|

1. Расчет на сайте http://matrixcalc.org/

Определитель = 712.000

Обратная матрица (-1)
|-0.093|00.090|00.042|
|00.090|-0.163|00.096|
|00.042|00.096|-0.065|

2. Расчет на Fortran п/п MINV
(см. код ниже)

Определитель = 316.364

Обратная матрица (-1)
|000.043|-0.153|00.644|
|-00.095|00.237|-0.392|
|-10.545|-0.009|00.056|

---------------------------

Используется синтаксис Fortran

! Обращение матрицы (общего вида)
! A - входная матрица, заменяется обратной 
! N - порядок матрицы A
! D - вычисляемый определитель
! L,M - рабочие векторы размерности (N)
      SUBROUTINE MINV(A,N,D,L,M)
      DIMENSION A(1),L(1),M(1)
      D=1.0
      NK=-N
      DO 80 K=1,N
      NK=NK+N
      L(K)=K
      M(K)=K
      KK=NK+K
      BIGA=A(KK)
      DO 20 J=K,N
      IZ=N*(J-1)
      DO 20 I=K,N
      IJ=IZ+I
   10 IF( ABS(BIGA)- ABS(A(IJ))) 15,20,20
   15 BIGA=A(IJ)
      L(K)=I
      M(K)=J
   20 CONTINUE
      J=L(K)
      IF(J-K) 35,35,25
   25 KI=K-N
      DO 30 I=1,N
      KI=KI+N
      HOLD=-A(KI)
      JI=KI-K+J
      A(KI)=A(JI)
   30 A(JI) =HOLD
   35 I=M(K)
      IF(I-K) 45,45,38
   38 JP=N*(I-1)
      DO 40 J=1,N
      JK=NK+J
      JI=JP+J
      HOLD=-A(JK)
      A(JK)=A(JI)
   40 A(JI) =HOLD
   45 IF(BIGA) 48,46,48
   46 D=0.0
      RETURN
   48 DO 55 I=1,N
      IF(I-K) 50,55,50
   50 IK=NK+I
      A(IK)=A(IK)/(-BIGA)
   55 CONTINUE
      DO 65 I=1,N
      IK=NK+I
      HOLD=A(IK)
      IJ=I-N
      DO 65 J=1,N
      IJ=IJ+N
      IF(I-K) 60,65,60
   60 IF(J-K) 62,65,62
   62 KJ=IJ-I+K
      A(IJ)=HOLD*A(KJ)+A(IJ)
   65 CONTINUE
      KJ=K-N
      DO 75 J=1,N
      KJ=KJ+N
      IF(J-K) 70,75,70
   70 A(KJ)=A(KJ)/BIGA
   75 CONTINUE
      D=D*BIGA
      A(KK)=1.0/BIGA
   80 CONTINUE
      K=N
  100 K=(K-1)
      IF(K) 150,150,105
  105 I=L(K)
      IF(I-K) 120,120,108
  108 JQ=N*(K-1)
      JR=N*(I-1)
      DO 110 J=1,N
      JK=JQ+J
      HOLD=A(JK)
      JI=JR+J
      A(JK)=-A(JI)
  110 A(JI)=HOLD
  120 J=M(K)
      IF(J-K) 100,100,125
  125 KI=K-N
      DO 130 I=1,N
      KI=KI+N
      HOLD=A(KI)
      JI=KI-K+J
      A(KI)=-A(JI)
  130 A(JI)=HOLD
      GO TO 100
  150 RETURN
      END
!==============================
   program test_f
   dimension OMI(9),OM(3,3)
! Cчитывание из файла "input.txt" линейного массива OMI(9)
   open(11,file='C:\test_f\input.txt',status='old')
   read(11,*)OMI
   close(11)
! Преобращование линейного массива OMI(9) в массив OM(3,3)
      II=0
   DO 21 I1=1,3
      DO 21 J1=1,3
      II=II+1
   21 OM(I1,J1)=OMI(II)
  N=3 
  L=N
  M=N
! Вывод в файл входных данных
  open(12, file='C:\test_f\output.txt')
print 53,N,L,M
write (12,53)N,L,M
   DO 22 I1=1,3
! Печать входной матрицы OM(3,3)
write (12,55)(OM(I1,J1),J1=1,3)
22 print 55,(OM(I1,J1),J1=1,3)
! Вызов п/п обращения матрицы OM(3,3)
   CALL MINV(OM,N,DD,L,M)
! Обращение матрицы (общего вида)
! SUBROUTINE MINV(A,N,D,L,M)
! A - входная матрица, заменяется обратной 
! N - порядок матрицы A
! D - вычисляемый определитель
! L,M - рабочие векторы размерности (N)
!---
print 54,N,L,M,DD
write (12,54)N,L,M,DD
! Печать обратной матрицы OM(3,3)
   DO 23 I1=1,3
print 55,(OM(I1,J1),J1=1,3)
23 write (12,55)(OM(I1,J1),J1=1,3)
   close(12)
53 format ('N,L,M=',3i3)
54 format ('N,L,M,D=',3i3,f9.3)
55 format ('OM:',99f8.3)
   end program test_f
 

[Munin]

А руками? Матрица маленькая, числа целые.

[simplermethod]

Фортран неправильно считает!

Исходники MINV взял из 2-х мест (см. ниже) - считают одинаково.
1) http://www.cyberforum.ru/fortran/thread168070.html
2) Сборник научных программ на Фортране http://mirknig.com/knigi/programming/11 ... trane.html
3) "разноцветный" код можно посмотреть здесь:
topic84682.html

"input.txt"
1.,7.,11.,7.,3.,9.,11.,9.,5.

"output.txt"
Вот входная матрица и результаты расчета:
N,L,M= 3 3 3
OM: 1.000 7.000 11.000
OM: 7.000 3.000 9.000
OM: 11.000 9.000 5.000

Результаты расчета:
N,L,M,D= 3 2 1 316.364
OM: 0.043 -0.153 0.644
OM: -0.095 0.237 -0.392
OM: -10.545 -0.009 0.056

[Munin]

Первым делом:
1) всё ли вы правильно поняли в документации MINV, и как её использовать?
2) проверьте, всё ли набрали без ошибок.

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

[simplermethod]

L и M - это массивы!!!

-- 26.05.2014, 16:43 --

Проблема решена!

L,M - массивы!
Dimension L(N),M(N)

Матрица A подается в п/п в виде линейного массива

[mihailm]

simplermethod, так фортран правильно считает или нет?

[Munin]

Теперь, видимо, да :-)