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

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

Входная матрица
|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|


!==============================
! Обращение матрицы (общего вида)
! 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)
! Обращение матрицы (общего вида)
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 [/more] [/more]

Проверить очень просто - перемножив прямую и обратную матрицы должна получиться единичная матрица.
Путем расчета на калькуляторе у меня получился определитель 712, обратную матрицу на калькуляторе трудно )
Смотрите свою MINV и правильность входных данных, может, у вас строки со столбцами перепутаны.

Исходники MINV взял из 2-х мест (см. ниже) - считают одинаково.
1) http://www.cyberforum.ru/fortran/thread168070.html
2) Сборник научных программ на Фортране

"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


Добавлено:
2) Сборник научных программ на Фортране
http://mirknig.com/knigi/programming/1181545162-sbornik-nauchnyh-programm-na-fortrane.html

Добавлено:
"разноцветный" код можно посмотреть здесь:
http://dxdy.ru/topic84682.html

вот как работает твоя программа без каких-либо изменений...

выдача на консоль:
Код: 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 3 1 712.000
OM: -0.093 0.090 0.042
OM: 0.090 -0.163 0.096
OM: 0.042 0.096 -0.065

Я не увидел в тексте ни одного явного указания используемого типа данных, например, REAL*8.

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

вообще-то такие вопосы
http://forum.ru-board.com/topic.cgi?forum=33&topic=7860#1
вот здесь обсуждаются...

а вообще за такой стиль написания надо руки вырывать...
в основной программе L -- единичное целое, а в п/п шлётся и обрабатывается там да ещё возвращается как массив...

тaк что прежде чем кричать
Цитата:

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

посмотри на себя и на свой код...

akaGM

Ниче не понимаю!
кстати, в Subroutine MINV матрица похоже идет в виде одномерного массива
(в фортране по умолчанию сначала меняется строка, потом столбец).

Gold Member    Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
вот как работает твоя программа без каких-либо изменений...
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 3 1 712.000
OM: -0.093 0.090 0.042
OM: 0.090 -0.163 0.096
OM: 0.042 0.096 -0.065

Добавлено:
akaGM
в основной программе L -- единичное целое, а в п/п шлётся и обрабатывается там да ещё возвращается как массив...

Спасиб!
Да! L и M - это массивы!!!

program test_f
...
L=N
первое появление L в основной программе...
какого она типа и что из себя представляет?

CALL MINV(OM,N,DD,L,M)

...

SUBROUTINE MINV(A,N,D,L,M)
DIMENSION A(1),L(1),M(1)
а здесь?

DO 80 K=1,N
NK=NK+N
L(K)=K
да ещё вот так...

сидел всю ночь - видать мозги стали отключаться ((

Добавлено:
Проблема решена!

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

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

возьми себе за правило описывать _все_ переменные _до_ их использования...