Подпрограммы умножения, транспонирования и обращения матриц

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

Подпрограммы представлены на двух языках: Fortran и Delphi.

 1. Fortran

================ Подпрограмма умножения матриц =================================
A(M,N)- 1-я матрица, размерностью (M,N)

B(N,M)- 2-я матрица, размерностью (N,M)

С(K,K)- Матрица – результат умножения, размерностью (K,K)

SUBROUTINE MAMU(A,M,N,B,K,C)
DIMENSION A(M,N),B(N,M),C(K,K)
DO 2 I=1,M
DO 2 J=1,K
C(I,J)=0.
DO 1 L=1,N
1 C(I,J)=C(I,J)+A(I,L)*B(L,J)
2 CONTINUE
RETURN
END

================ Подпрограмма транспонирования матриц ==========================   

X(N,K)- 1-я матрица, размерностью (N,K)

XTR(K,N)- Транспонированная матрица, размерностью (K,N)

SUBROUTINE MTRAN(X,N,K,XTR)
DIMENSION X(N,K),XTR(K,N)
DO 1 I=1,N
DO 1 J=1,K
1 XTR(J,I)=X(I,J)
OPEN(3,FILE='XTR_FOR.TXT')
WRITE(3,541)X
541 FORMAT(7(F8.3,3X))
CLOSE(3)
RETURN
END</pre>

================ Подпрограмма обращения матриц =================================

A(N,N)- Входная матрица, размерностью (N,N)

E(N,N)- Обращенная матрица, размерностью (N,N)

SUBROUTINE MATDP(A,E,N)
DIMENSION A(N,N),E(N,N)
25 FORMAT(5X,'Сбой при обращении матрицы')
DO 11 I=1,N
DO 11 J=1,N
IF(I.EQ.J)GO TO 10
E(I,J)=0.
GO TO 11
10 E(I,J)=1.
11 CONTINUE
DO 27 I=1,N
J=I
L=I
LD=0
12 IF(ABS(A(L,J)))16,13,16
13 LD=1
IF(L-N)15,14,14
14 WRITE(3,25)
GOTO 27
15 L=L+1
GO TO 12
16 IF(LD-1)26,17,26
17 DO 18 K=1,N
T=A(L,K)
A(L,K)=A(I,K)
A(I,K)=T
T=E(L,K)
E(L,K)=E(I,K)
18 E(I,K)=T
26 BM=1./A(I,I)
DO 19 K=1,N
19 E(I,K)=E(I,K)*BM
I1=I+1
DO 20 KI=I1,N
K=N+I1-KI
20 A(I,K)=A(I,K)*BM
DO 23 L=1,N
IF(L.EQ.I)GOTO 23
BM=A(L,I)
DO 21 K=1,N
21 E(L,K)=E(L,K)-E(I,K)*BM
I1=I+1
DO 22 KI=I1,N
K=N+I1-KI
22 A(L,K)=A(L,K)-A(I,K)*BM
23 CONTINUE
27 CONTINUE
RETURN
END

2. Delphi

==================== MAMU – УМНОЖЕНИЕ МАТРИЦ ========================

A(M10,N10)- 1-я матрица, размерностью (M10,N10)

B(N10,M10)- 2-я матрица, размерностью (N10,M10)

С(K10,K10)- Матрица – результат умножения, размерностью (K10,K10)

mas=array[1..M,1..N] of real;

Procedure MAMU(A:MAS;M10,N10:INTEGER;B:MAS;K10:INTEGER; VAR C:MAS);
Var
I,J,L:integer;
BEGIN
  For I:=1 To M10 do
  Begin
    For J:=1 To K10 do
    Begin
      C[I,J]:=0.0;
      For L:=1 To N10 do C[I,J]:=C[I,J]+A[I,L]*B[L,J];
    End;
  End;
END;

//========= Транспонирование матриц MTRAN =================================

mas=array[1..M,1..N] of real;

X(N,K)- 1-я матрица, размерностью (N,K)

XTR(K,N)- Транспонированная матрица, размерностью (K,N)

procedure MTRAN(X:mas;N,K:integer; var XTR:mas);
Var
  I,J:integer;
BEGIN
  for I := 1 to N do
  for j := 1 to K do
  XTR[J,I]:=X[I,J];
END;

================== MATDP – ОБРАЩЕНИЕ МАТРИЦ =========================

A(N,N)- Входная матрица, размерностью (N,N)

E(N,N)- Обращенная матрица, размерностью (N,N)

mas=array[1..M,1..N] of real;

Procedure MATDP(A:mas; var E:mas;N:integer);
Var
  L,LD,I1,KI,K : integer;
  T,BM : real;
BEGIN
   for I := 1 to N do
   for j := 1 to N do IF (I<>J) then E[i,j]:=0 else E[i,j]:=1;
   for I:=1 to N do
   Begin
     J:=I;
     L:=I;
     LD:=0;
     IF (ABS(A[L,J])=0) then
     begin
       LD:=1;
       IF (L-N=0) then
       Begin
         ShowMessage('Сбой при обращении матрицы! При вычислении выражения
         обнаружена сингулярность.'); Exit;
       End
       else L:=L+1;
     end;
     IF((LD-1)=0) then
     Begin
       For K:=1 To N do
       begin
         T:=A[L,K];
         A[L,K]:=A[I,K];
         A[I,K]:=T;
         T:=E[L,K];
         E[L,K]:=E[I,K];
         E[I,K]:=T;
       end;
     End;
     if A[I,I]<>0 then BM:=1./A[I,I] else BM:=0.0001;
     For K:=1 to N do E[I,K]:=E[I,K]*BM;
     I1:=I+1;
     for KI:=I1 to N do
     begin
       K:=N+I1-KI;
       A[I,K]:=A[I,K]*BM;
     end;
     For L:=1 to N do
     Begin
       if (L<>I) then
       begin
         BM:=A[L,I];
         for K:=1 to N do E[L,K]:=E[L,K]-E[I,K]*BM;
         I1:=I+1;
         for KI:=I1 to N do
         begin
           K:=N+I1-KI;
           A[L,K]:=A[L,K]-A[I,K]*BM;
         end;
       end;
     End;
   End;
END;