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 имеют следующий вид dii = 1 / uii
j-1
dij = -( ∑ dikukj) / djj (j = i+1, i+2, ..., N)
k=i
для i от 1 до N.
Матрица C вычисляется на втором этапе исходя из уравнения CL = D по формулам
N
cij = dij - ∑ ciklkj (i = 1, 2, ..., j)
k=j+1
N
wi = - ∑ ciklkj
k=j+1
cij = wi (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
LA52
Программа LA52 предназначена для обращения симметричной матрицы с ненулевыми главными минорами с помщью изменённого метода квадратных корней. Преимущество этого метода состоит в том, что требуется только n рабочих ячеек, не нужно никакой единичной матрицы, не извлекаются никакие квадратные корни, выполняется только n делений. Алгоритм основан на решении уравнения
CLTDL = E,
где C - искомая матрица, E - единичная матрица, LTDL = A - разложение матрицы A на верхнюю треугольную L, транспонированную к ней нижнюю треугольную LT и диагональную D. Вычисление матрицы C производится одновременно с разложением матрицы A по циклическому алгоритму
ajN = a1,j+1 / a11 (j = 1, 2, ..., N-1), aij = ai+1,j+1 - ai+1,1a1,j+1 / a11 (j = i, ..., N-1)
для всех i = 1, ..., N-1, aNN = -1 / a11
Обратная матрица C получается после изменения знака у вычисленной матрицы. Элемент a11 и все главные миноры исходной матрицы должны быть отличны от нуля.
Алгоритм 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
AIG2R
Программа AIG2R обращает матрицу A методом Жордана с выбором ведущего элемента по всей матрице. Матрица A должна быть неособенная, иначе программа завершится аварийно.
Исходная матрица A порядка NxN приводится к единичной матрице E при помощи последовательности элементарных преобразований Жордана L1, L2, ... LN и матрицами перестановок Р и Q так, что
LNLN-1 ... L1PAQ = E
при этом А-1 = QLNLN-1 ... L1Р. Матрицы P и Q осуществляют соответственно перестановку строк и столбцов матрицы A и обеспечивают стратегию выбора ведущего элемента по всей матрице.
http://num-anal.srcc.msu.ru/lib_na/cat/ai/aig2r.htm
Вызов программы
call AIG2R(A, N)
Параметры программы
N - входной параметр;
A - входной и выходной параметр;
Real A(1:N, 1:N) - массив, содержащий исходную матрицу, на выходе из программы содержит обратную матрицу;
Integer N - размерность матрицы A.
Пример
!Обращение матрицы методом Жордана с выбором ведущего элемента по всей матрице
program Jordan
use NML
implicit none
integer, parameter:: n=4
real A(n,n), AS(n,n), AZ(n,n)
data A /7.9, 8.5, 4.3, 3.2, 5.6, -4.8, 4.2, -1.4, &
5.7, .8, -3.2, -8.9, -7.2, 3.5, 9.3, 3.3/
integer i, j
!
AS=A
call AIG2R (A, n)
AZ=matmul(A,AS)
!
print 4
print 10, ((AS(i,j), j=1,n), i=1,n)
print 5
print 10, ((AZ(i,j), j=1,n), i=1,n)
!
4 format(/4X,'The initial matrix')
5 format(/4X,'The result of multiplying the initial matrix with the inverse')
10 format(/(4X,<n>F10.5))
end program Jordan
The initial matrix
7.90000 5.60000 5.70000 -7.20000
8.50000 -4.80000 0.80000 3.50000
4.30000 4.20000 -3.20000 9.30000
3.20000 -1.40000 -8.90000 3.30000
The result of multiplying the initial matrix with the inverse
1.00000 0.00000 0.00000 -0.00000
0.00000 1.00000 0.00000 -0.00000
0.00000 0.00000 1.00000 0.00000
0.00000 -0.00000 0.00000 1.00000