Автор cайта: Владимир Потемкин fortran-90@yandex.ru

Обращение матриц

LA50. Обращение матрицы A, преобразованной в произведение двух треугольных матриц L и U при помощи программы LA10.

LA52. Обращение симметричной матрицы А с ненулевыми главными минорами на основе L^TDL - разложения.

LA50

Программа LA50 вычисляет матрицу C, обратную к матрице A (C=A^-1). Матрица A должна быть предварительно преобразована в произведение двух треугольных матриц L и U при помощи программы LA10. Программа LA10 должна закончиться без ошибок (Error=0).
Матрица C, обратная к матрице A, находится из уравнения
CLU = E,
где L и U - треугольные матрицы, E - единичная матрица.
Вычисление матрицы C проводится в два этапа. На первом этапе вычисляется матрица D исходя из уравнения DU = E. Расчетные формулы для вычисления матрицы D имеют следующий вид

d_ii = 1 / u_ii
j-1
d_ij = -( ∑ d_iku_kj) / d_jj    (j = i+1, i+2, ..., N)
^k=i
для i от 1 до N.
   Матрица C вычисляется на втором этапе исходя из уравнения CL = D по формулам
_N
c_ij = d_ij - ∑ c_ikl_kj    (i = 1, 2, ..., j)
^k=j+1
_N
w_i = - ∑ c_ikl_kj
^k=j+1
c_ij = w_i    (i = j+1, j+2, ..., N)
для j от N-1 до 1.
После завершения вычислений в матрице C переставляются столбцы согласно массиву T, который содержит информацию о перестановках строк при работе программы LA10.

Вызов программы

CALL LA50(AF, T)

Параметры программы

AF - входной и выходной параметр;
T - входной параметр.

Real AF(1:N, 1:N) - массив, содержащий элементы треугольных матриц L и U. На выходе массив содержит элементы матрицы C;
Integer T(1:N) - массив, содержащий информацию о перестановках строк матрицы A.

Пример

!Обращение матрицы
program TestLA50
use ConvertCyr
use NML
implicit none
integer, parameter:: n=4
real:: A(n,n), AF(n,n), D, Z(n,n)
integer:: T(n), Er
data A/1.0, 1.1, 1.2, 1.4, 1.1, 1.1, 1.2, 1.3, 1.2, 1.2, 1.2, 1.3, 1.4, 1.3, 1.3, 1.3/
integer i, j
!begin
  call LA10(A, AF, T, D, Er)
  if (Er>0) then
    print *, WinToDos('Матрица особенная. Программа остановлена.')
    stop
  end if
  call LA50(AF, T)
  Z=matmul(A,AF)
  print 10
  print 11, ((A(i,j), j=1,n), i=1,n)
  print 12
  print 13, ((Z(i,j), j=1,n), i=1,n)
  print 14, D
10 format(/'  The initial matrix:'/)
11 format(8x,<n>F5.1)
12 format(/'  The resulting identity matrix:'/)
13 format(8x,<n>F11.7)
14 format(/'  Determinant =',E15.7)
end

The initial matrix:

        1.0  1.1  1.2  1.4
        1.1  1.1  1.2  1.3
        1.2  1.2  1.2  1.3
        1.4  1.3  1.3  1.3

The resulting identity matrix:

        0.9999983  0.0000078  0.0000061 -0.0000154
        0.0000004  1.0000042  0.0000045 -0.0000080
        0.0000004  0.0000042  1.0000045 -0.0000233
       -0.0000149  0.0000042  0.0000197  0.9999768

Determinant = -0.9999890E-04

Вернуться к оглавлению    Скачать LA50

LA52

Программа LA52 предназначена для обращения симметричной матрицы с ненулевыми главными минорами с помщью изменённого метода квадратных корней. Преимущество этого метода состоит в том, что требуется только n рабочих ячеек, не нужно никакой единичной матрицы, не извлекаются никакие квадратные корни, выполняется только n делений. Алгоритм основан на решении уравнения
CL^TDL = E,
где C - искомая матрица, E - единичная матрица, L^TDL = A - разложение матрицы A на верхнюю треугольную L, транспонированную к ней нижнюю треугольную L^T и диагональную D. Вычисление матрицы C производится одновременно с разложением матрицы A по циклическому алгоритму
a_jN = a_1,j+1 / a₁₁    (j = 1, 2, ..., N-1),
a_ij = a_i+1,j+1 - a_i+1,1a_1,j+1 / a₁₁    (j = i, ..., N-1)
для всех i = 1, ..., N-1,
a_NN = -1 / a₁₁
Обратная матрица C получается после изменения знака у вычисленной матрицы. Элемент a₁₁ и все главные миноры исходной матрицы должны быть отличны от нуля.

Алгоритм 66б из [Б9]

Вызов программы

CALL LA52(A, Error)

Параметры программы

A - входной и выходной параметр;
   Error - выходной параметр.

   Real A(1:N, 1:N) - массив, содержащий элементы исходной матрицы. На выходе выше главной диагонали и на главной диагонали находятся элементы обратной матрицы. Элементы ниже главной диагонали не используются и остаются без изменений;
   Integer Error - индикатор ошибки. Error=0, если ошибок нет, и равно 65, если программа завершилась аврийно.

Пример

!Обращение симметричной матрицы
program TestLA52
use ConvertCyr
use NML
implicit none
integer, parameter:: n=4
real:: A(n,n), AA(n,n), Z(n,n)
integer:: T(n), Er
data A/1.0, 1.1, 1.2, 1.4, 1.1, 1.1, 1.2, 1.3, 1.2, 1.2, 1.2, 1.3, 1.4, 1.3, 1.3, 1.3/
integer i, j
!begin
  AA=A
  call LA52(AA, Er)
  if (Er>0) then
    print *, WinToDos('Матрица особенная. Программа остановлена.')
    stop
  end if
  forall(i=2:n, j=1:n-1, i>j) AA(i,j)=AA(j,i)
  Z=matmul(A,AA)
  print 10
  print 11, ((A(i,j), j=1,n), i=1,n)
  print 12
  print 13, ((Z(i,j), j=1,n), i=1,n)
10 format(/'  The initial matrix:'/)
11 format(8x,<n>F5.1)
12 format(/'  The resulting identity matrix:'/)
13 format(8x,<n>F11.7)
end

The initial matrix:

        1.0  1.1  1.2  1.4
        1.1  1.1  1.2  1.3
        1.2  1.2  1.2  1.3
        1.4  1.3  1.3  1.3

The resulting identity matrix:

        0.9999983 -0.0000074 -0.0000062 -0.0000078
        0.0000004  0.9999889 -0.0000047 -0.0000042
        0.0000004 -0.0000111  0.9999954 -0.0000042
        0.0000004 -0.0000111 -0.0000047  0.9999958

Вернуться к оглавлению    Скачать LA52