|
Invert3dMatrix, Библиотека "JINRLIB" F499
Lin2dsysccmSolver,
Lin3dsysccmSolver,
LinsysccmSolver,
PseudsysccmSolver,
Init_Const Вы
Авторы: Г.А.Емельяненко, Э.Б.Душанов, М.Г.Емельяненко, 
Т.Т.Рахмонов, А.П.Сапожников посетитель
Язык: Фортран
JINRLINPACK - МАШИННО-НЕЗАВИСИМЫЙ ПАКЕТ ПРОГРАММ ДЛЯ ОБРАЩЕНИЯ
ДВУХДИАГОНАЛЬНЫХ, ТРЕХДИАГОНАЛЬНЫХ И ПРЯМОУГОЛЬНЫХ (ЗАПОЛНЕННЫХ) МАТРИЦ
И РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ УРАВНЕНИЙ
С СООТВЕТСТВУЮЩИМИ МАТРИЦАМИ КОЭФФИЦИЕНТОВ
Init_Const - программа получения констант вещественной арифметики ЭВМ
(подробное описание).
Invert3dMatrix - программа обращения трехдиагональных матриц общего вида,
включая плохо обусловленные (подробное описание).
Lin2dsysccmSolver - программа обращения двухдиагональных матриц,
вычисления определителя и решения системы линейных алгебраических
уравнений с двухдиагональной матрицей коэффициентов (подробное описание).
Lin3dsysccmSolver - программа обращения трехдиагональных матриц, решения
систем линейных алгебраических уравнений с трехдиагональной матрицей
коэффициентов, а также вычисления определителя матрицы (подробное описание).
LinsysccmSolver - программа нахождения матрицы, обратной прямоугольной,
вычисления определителя и решения системы линейных алгебраических
уравнений с заполненными матрицами коэффициентов (подробное описание).
PseudsysccmSolver - программа решения систем линейных уравнений AX=Y
со слабовырожденной матрицей A (у A имеется лишь одно нулевое собственное
значение) (подробное описание).
Init_Const - ОПРЕДЕЛЕНИЕ КОНСТАНТ ВЕЩЕСТВЕННОЙ АРИФМЕТИКИ ЭВМ
Программа Init_Const определяет константы вещественной арифметики ЭВМ:
beta - основание системы счисления,
(l<0,u>0) - нижнюю и верхнюю границы порядков представления
нормализованного числа в ЭВМ,
t - число beta-ичных разрядов, отведенных для хранения
мантиссы числа в ЭВМ,
e_8 - максимальное число, представимое в данной ЭВМ,
e_0 - минимальное число, представимое в данной ЭВМ,
e_1 - относительную погрешность вычислений с нормализованными
числами данной ЭВМ,
а также получает параметры настройки функций GTEXP,CHEXP [1].
Структура:
----------
Тип: SUBROUTINE
Имена входа для пользователя: INIT_CONST,INITP,GTEXP,CHEXP
Общие блоки: COMMON /F499_RCONST/ OVFLO,UNDFLO,EPS,EPSMIN
COMMON /F499_ICONST/ RADIX,MAXEXP,MINEXP,MANTSZ
COMMON /F499_ICHGT/ BINV,INO,ILT,MNE1,MNE2,INB,LWT,KB
Обpащение:
----------
CALL INIT_CONST, если вывод вычисленных констант не требуется;
CALL INITP, если константы нужно вывести.
Выходные данные:
----------------
вещественные константы e_8,e_0,e_1,e_2=e_1/beta
и целые константы beta,u,l,t. Выдаются в блоках F499_RCONST
и F499_ICONST в указанных последовательностях соответственно.
В блоке F499_ICHGT содержатся следующие выделенные параметры
настройки функции GTEXP, CHEXP:
BINV - (real*8) вещественное число, равное 1/beta;
KB - (integer) номер половины формата, в котором хранится
младший разряд порядка нормализованного числа;
LWT - (integer) представление целого порядка e=0 в разрядах
формата, отведенных для порядка;
INO - (integer) количество beta-ичных разрядов, отведенных для
порядка с учетом его знака (в формате для хранения
вещественного числа);
INB - (integer) количество beta-ичных разрядов, отведенных для
мантиссы с учетом ее знака (в формате для хранения
вещественного числа);
MNE1 - (integer) признак кода (e>0)-положительного порядка.
При этом MNE1=[e]^k_(код)-e;
MNE2 - (integer) признак кода (e<0)-отрицательного порядка.
При этом MNE2=[e]^k_(код)-e.
Подпрограмма-функция GTEXP - выделение мантиссы и порядка
вещественного числа.
Общие блоки: COMMON /F499_ICHGT/ BINV,INO,ILT,MNE1,MNE2,INB,LWT,KB
Обpащение: RM=GTEXP(R,E)
Входные данные: R - (real*8) вещественное число, порядок и мантисса
которого выделяется. Блок F499_ICHGT содержит настроечные
параметры, определенные в Init_Const (или заранее
известные);
Выходные данные: RM - (real*8) содержит выделенную мантиссу числа R;
E - (integer) содержит выделенный порядок числа R
Подпрограмма-функция CHEXP - восстановление вещественного числа
по его мантиссе и порядку.
Общие блоки: COMMON /F499_ICONST/ RADIX,MAXEXP,MINEXP,MANTSZ
COMMON /F499_ICHGT/ BINV,INO,ILT,MNE1,MNE2,INB,LWT,KB
Обpащение: R=CHEXP(RM,E)
Входные данные: RM - (real*8) мантисса числа R; E - (integer) порядок
числа R. Блоки F499_ICONST и F499_ICHGT содержат элементы,
определенные в Init_Const (или заранее известные);
Выходные данные: R - (real*8) содержит вещественное число beta^E RM=R.
Метод: алгоритм автоматического получения констант вещественной
арифметики ЭВМ [1].
Пример:
-------
Пусть, например, ЭВМ есть IBM PC. Известно, что машинные константы
компьютеров такого типа есть: beta=2,t=53,l=-1021,u=1024.
Построим программу TEST, которая, используя программы INIT_CONST,GTEXP и
CHEXP, определит эти константы, а также константы e_1=1/beta*beta^{2-t} -
относительная погрешность машинной арифметики, e_0=1/beta*beta^l -
"машинный" нуль и e_8=(1-1/beta^t)*beta^u - "машинная" бесконечность,
и при этом выделит мантиссы и порядки вещественных констант:
1,e_1,e_0,e_8, а также обратно восстановит эти
константы по выделенным мантиссам и порядкам.
Ниже приводится текст программы этой задачи, а также таблицы 1, 2 и 3 -
результаты работы программы. При этом в программе TEST использована
подпрограмма SUB, которая восстанавливает одновременно группу указанных
чисел по выделенной мантиссе и порядку. Программа TEST печатает
результаты своей работы - таблицы 1, 2, 3.
PROGRAM TEST
IMPLICIT REAL*8 (A-H, O-Z)
INTEGER RADIX
COMMON /F499_RCONST/ OVFLO,UNDFLO,EPS,EPSMIN
COMMON /F499_ICONST/ RADIX,MAXEXP,MINEXP,MANTSZ
CALL INITP
D1=GTEXP(1D0,I)
DE=GTEXP(EPS,J)
DO=GTEXP(OVFLO,L)
DU=GTEXP(UNDFLO,K)
D=1D0
PRINT1,D,D1,I,EPS,DE,J,UNDFLO,DU,K,OVFLO,DO,L,
1 FORMAT(10X,'ЧИСЛО',10X,':',6X,'МАНТИССА',6X,': ПОРЯДОК',/,
* 25('-'),
* '+',20('-'),'+',8('-'),4(/,E24.15E3,' :',F19.16,' : ',I5),/)
CALL SUB(D1,I,DE,J,DO,L,DU,K)
STOP
END
C
SUBROUTINE SUB(D,I,DE,J,DO,L,DU,K)
IMPLICIT REAL*8 (A-H, O-Z)
D1=CHEXP(D,I)
EPS=CHEXP(DE,J)
OVFL=CHEXP(DO,L)
UNDF=CHEXP(DU,K)
PRINT1,D,I,D1,DE,J,EPS,DE,K,UNDF,DO,L,OVFL
1 FORMAT(6X,'МАНТИССА',6X,': ПОРЯДОК:',10X,'ЧИСЛО',/,20('-'),
* '+',8('-'),'+',25('-'),4(/,F19.16,' : ',I5,' :',E24.15E3))
RETURN
END
Таблица 1 констант вещественной арифметики ЭВМ:
The constants of mashin`s arithmetic ...
MAXEXP= 1024 OVFLO= 1.797693134862316E+308
MINEXP= -1021 UNDFLO= 2.225073858507201E-308
RADIX= 2 EPS= 2.220446049250313E-016
MANTSZ= 53 EPSMIN= 1.110223024625157E-016
Таблица 2 выделенных мантисс и порядков вещественных чисел:
1=1/beta*beta^1,e_1=1/beta*beta^{2-t},e_0=1/beta*beta^l и
e_8=(1-1/beta^t)*beta^u
ЧИСЛО : МАНТИССА : ПОРЯДОК
------------------------+--------------------+--------
.100000000000000E+001 : .5000000000000000 : 1
.222044604925031E-015 : .5000000000000000 : -51
.222507385850720E-307 : .5000000000000000 : -1021
.179769313486232E+309 : .9999999999999999 : 1024
Таблица 3 восстановленных по мантиссам и порядкам вещественных чисел:
1/beta*beta^1=1,1/beta*beta^{2-t}=e_1,1/beta*beta^l=e_0 и
(1-1/beta^t)*beta^u=e_8
МАНТИССА : ПОРЯДОК: ЧИСЛО
-------------------+--------+-------------------------
.5000000000000000 : 1 : .100000000000000E+001
.5000000000000000 : -51 : .222044604925031E-015
.5000000000000000 : -1021 : .222507385850720E-307
.9999999999999999 : 1024 : .179769313486232E+309
Invert3dMatrix - ОБРАЩЕНИЕ ТРЕХДИАГОНАЛЬНЫХ МАТРИЦ ОБЩЕГО ВИДА
(ВКЛЮЧАЯ ПЛОХООБУСЛОВЛЕННЫЕ)
Программа Invert3dMatrix - обращение трехдиагональных матриц C_3.
Структура:
----------
Тип: SUBROUTINE
Имена входа для пользователя: INVERT3DMATRIX
Общие блоки: COMMON /F499_RCONST/ OVFLO,UNDFLO,EPS,EPSMIN
Обpащение:
----------
CALL Invert3dMatrix(M,A,INF,R1,R2,R,IZ)
Входные данные: блок F499_RCONST содержит вещественные константы
e_8,e_0,e_1,e_2=e_1/beta, определенные в Init_Const
(или заранее известные);
M - (integer) поpядок квадратной матpицы
_ _
| q_1 r_2 0|
| p_2 q_2 \ |
C_3=| \ \ r_M |; (1)
| 0 p_M q_M |
- -
A - (real*8) массив размерности [M,M], который на входе
содержит заданную матрицу, а на выходе - инвертированную
матрицу;
INF - (integer) выходной параметр:
INF= 0 при нормальном завершении работы программы;
INF=-1, если исходная матрица вырождена.
R1,R2,R - (real*8) одномерные рабочие массивы размерности M;
IZ - (integer) одномерный рабочий массив размерности M.
Используется метод критических компонент [2-5].
Lin2dsysccmSolver - ОБРАЩЕНИЕ ДВУХДИАГОНАЛЬНЫХ МАТРИЦ
И РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ
С ДВУХДИАГОНАЛЬНОЙ МАТРИЦЕЙ КОЭФФИЦИЕНТОВ
Программа Lin2dsysccmSolver - обращение двухдиагональных матриц C_2
и решение систем линейных уравнений C_2X=Y (метод критических
компонент [2-5]).
Структура:
----------
Тип: SUBROUTINE
Имена входа для пользователя: LIN2DSYSCCMSOLVER
Общие блоки: COMMON /F499_CONST/ OVFLO,UNDFLO,EPS,EPSMIN
Обpащение:
----------
CALL Lin2dsysccmSolver(M,R,N,INF,IM,JM,B,DET)
Входные данные:
блок F499_RCONST содержит вещественные константы
e_8,e_0,e_1,e_2=e_1/beta, определенные в Init_Const
(или заранее известные);
M - (integer) поpядок квадратной матpицы
_ _
| q_1 r_2 0|
| q_2 \ |
C_2=| \ r_M |;
| 0 q_M |
- -
N - (integer) первая размерность массива B (N=M);
R - (real*8) одномерный (или двумеpный) массив
pазмеpности 2M. На входе в массиве R размещается
исходная матрица C_2 в виде:
R=[*,q_1,r_2,q_2,...,r_{M-1},q_{M-1},r_M,q_M]. (*)
Здесь {r_i}^m_{i=2} - внедиагональные элементы матрицы
C_2, символом * отмечены ячейки, которые при вводе не
заполняются;
INF - (integer): INF=1, если задача решается с
верхнедвухдиагональной матрицей; INF=-1, если задача
решается с нижнедвухдиагональной матpицей;
IM,JM - (integer). Если JM=0, то решается система C_2X=Y. При
этом IM - число правых частей. Если же JM=k (k<=IM),
то вычисляются элементы обратной к C_2 матрицы,
находящиеся в столбцах с номерами k,k+1,...,IM. При
этом k - номер последнего вычисляемого столбца
обратной матрицы;
B - (real*8) двумеpный массив. Если отведенная память
компьютера не достаточна для хранения массива B
размерности [N,IM], то обращение к подпрограмме
Lin2dsysccmSolver организуется самим пользователем в
цикле с соответствующим выделением из B блоков
нужных размеров. В частности, B может быть одномерным,
если IM=1). Если решается система C_2X=Y, т.е. JM=0,
то массив B имеет размерность [N,IM]. При этом в B
хранятся векторы правых частей Y_j в виде
Y_1 Y_2 ... Y_{IM}
_ _
| b_{11} b_{12} ... b_{1IM} |
| b_{21} b_{22} ... b_{2IM} |
B=| : : : : |. (**)
| b_{N1} b_{N2} ... b_{NIM} |
- -
Если вычисляются элементы обратной к C_2 матрицы,
т.е. JM!=0, то массив B имеет размерность [N,IM-JM+1];
DET - (real*8), если на входе DET=1, то на выходе
DET=det(C_2) - определитель матрицы C_2, если на входе
DET=0, то вычисления определителя не производится.
Выходные данные:
INF - (integer):
INF=0 - ноpмальное завеpшение pаботы подпpогpамм;
INF=-1 - исходная матpица выpождена;
B - (real*8) массив содеpжит матpицу pешений и элементов
обратной матрицы в виде (**);
DET - (real*8) содержит вычисленный определитель матрицы C_2,
если на входе DET=1.
Используется метод критических компонент [2-5].
Пример:
-------
Вычислить решения системы [1,4]:
_ _ _ _
| 7/5 11/3 | | y_1 |
| 7/5 11/3 | | y_2 |
| \ \ | X= | : |, (1)
| 7/5 11/3 | | : |
| 7/5 | | y_M |
- - - -
y_i=(152i+118)/(15(2i+1)(2i+3)),y_M=7/(5(2M+1)),i=1,2,...,M-1.
Вычисление обратной матрицы и определителя матрицы системы (1).
Ниже приводится текст программы решения этой задачи, а также
таблица - результаты работы программы.
PROGRAM TEST
IMPLICIT REAL*8 (A-H, O-Z)
INTEGER RADIX
DIMENSION A(10,10),XB(10,5),R(15),X(10)
COMMON /F499_RCONST/ OVFLO,UNDFLO,EPS,EPSMIN
COMMON /F499_ICONST/ RADIX,MAXEXP,MINEXP,MANTSZ
CALL INIT_CONST
M=3
N=10
DET=1D0
JM=-1
IM=1
2 JM=JM+1
INF=1
DO I=1,3
R(2*I)=1D0
R(2*I-1)=2D0
X(I)=XB(I,1)
ENDDO
XB(1,1)=2D0
XB(2,1)=7D0/6D0
XB(3,1)=1D0/3D0
CALL Lin2dsysccmSolver(M,R,N,INF,IM,JM,XB,DET)
IM=M
IF(JM.EQ.0) GOTO 2
PRINT 7
7 FORMAT('Решение системы (1) и обращение (двухдиагональной)'
*' матрицы этой системы'/' (Программа Lin2dsysccmSolver)'/
*' Решение Обратная матрица')
write(6,5) (X(I),(XB(I,J),J=1,IM),I=1,M)
write(6,6) INF,DET
5 FORMAT(D10.3,' : ',3d10.3)
6 FORMAT(' INF=',I3,', DET=',d10.4)
STOP
END
Таблица решения системы (1) и обратной матрицы к матрице системы.
Решение Обратная матрица
.100D+01 : .100D+01 -.200D+01 .400D+01
.500D+00 : .000D+00 .100D+01 -.200D+01
.333D+00 : .000D+00 .000D+00 .100D+01
INF= 0, DET= .1000D+01
Lin3dsysccmSolver - ОБРАЩЕНИЕ ТРЕХДИАГОНАЛЬНЫХ МАТРИЦ
И РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ
С ТРЕХДИАГОНАЛЬНОЙ МАТРИЦЕЙ КОЭФФИЦИЕНТОВ
Программа Lin3dsysccmSolver - обращение трехдиагональных матриц C_3
и решение систем линейных уравнений C_3X=Y (метод критических
компонент [2-5]).
Структура:
----------
Тип: SUBROUTINE
Имена входа для пользователя: LIN3DSYSCCMSOLVER
Общие блоки: COMMON /F499_RCONST/ OVFLO,UNDFLO,EPS,EPSMIN
Обpащение:
----------
CALL Lin3dsysccmSolver(M,R,N,INF,IM,JM,B,A,DET)
Входные данные: блок F499_RCONST содержит вещественные константы
e_8,e_0,e_1,e_2=e_1/beta, определенные в Init_Const
(или заранее известные);
M - (integer) поpядок квадратной матpицы
_ _
| q_1 r_2 0|
| p_2 q_2 \ |
C_3=| \ \ r_M |;
| 0 p_M q_M |
- -
N - (integer) первая размерность массива B (N=M);
R - (real*8) одномерный (или двумеpный) массив
pазмеpности 3M. На входе в массиве R размещается
исходная матрица C_3, если она несимметричная, в виде:
R=[*,q_1,p_2,r_2,q_2,...,p_M,r_M,q_M].
Если она симметричная, то в виде
R=[*,q_1,r_2,q_2,...,r_M,q_M].
Здесь символом * отмечены ячейки, которые при вводе не
заполняются;
INF - (integer): INF=1, если задача решается с симметричной
трехдиагональной матрицей; INF=-1, если задача решается
с несимметричной трехдиагональной матpицей;
IM,JM - (integer). Если JM=0, то решается система C_3X=Y. При
этом IM - число правых частей. Если же JM=k (k<=IM),
то вычисляются элементы обратной к C_3 матрицы,
находящиеся в столбцах с номерами k,k+1,...,IM. При
этом k - номер последнего вычисляемого столбца
обратной матрицы;
B - (real*8) двумеpный массив. Если отведенная память
компьютера не достаточна для хранения массива B
размерности [N,IM], то обращение к подпрограмме
Lin3dsysccmSolver организуется самим пользователем в
цикле с соответствующим выделением из B блоков
нужных размеров. В частности, B может быть одномерным,
если IM=1). Если решается система C_3X=Y, т.е. JM=0,
то массив B имеет размерность [N,IM]. При этом в B
хранятся векторы правых частей Y_j в виде
Y_1 Y_2 ... Y_{IM}
_ _
| b_{11} b_{12} ... b_{1IM} |
| b_{21} b_{22} ... b_{2IM} |
B=| : : : : |. (*)
| b_{N1} b_{N2} ... b_{NIM} |
- -
Если вычисляются элементы обратной к C_3 матрицы,
т.е. JM!=0, то массив B имеет размерность [N,IM-JM+1];
A - (real*8) рабочий массив размерности [3,M];
DET - (real*8) если на входе DET=1, то на выходе
DET=det(C_3) - определитель матрицы C_3, если на входе
DET=0, то вычисления определителя не производится.
Выходные данные:
INF - (integer):
INF=0 - ноpмальное завеpшение pаботы подпpогpамм;
INF=-1 - исходная матpица выpождена;
B - (real*8) массив содеpжит матpицу pешений и элементов
обратной матрицы в виде (*);
DET - (real*8) содержит вычисленный определитель матрицы C_3,
если на входе DET=1.
Используется метод критических компонент [2-5].
Пример:
-------
Вычислить решения системы [1,4]:
_ _ _ _
| 6 3 | | 9 |
| 4 6 \ | | 13 |
| \ \ 3 | X= | : |, (1)
| 4 6 2 | | 13 |
| 4 6 | | 10 |
- - - -
Вычислить обратную матрицу и определитель матрицы системы (1).
Ниже приводится текст программы решения этой задачи, а также
таблица - результаты работы программы.
PROGRAM TEST
IMPLICIT REAL*8 (A-H, O-Z)
INTEGER RADIX
DIMENSION A(10,10),XB(10,5),R(15),X(10)
COMMON /F499_RCONST/ OVFLO,UNDFLO,EPS,EPSMIN
COMMON /F499_ICONST/ RADIX,MAXEXP,MINEXP,MANTSZ
CALL INIT_CONST
M=3
N=10
det=1.
JM=-1
IM=1
3 JM=JM+1
inf=-1
DO I=1,3
R(3*I)=3D0
R(3*I-1)=6D0
R(3*I-2)=4D0
X(I)=XB(I,1)
ENDDO
XB(1,1)=10D0
XB(2,1)=13D0
XB(3,1)=9D0
CALL Lin3dsysccmSolver(M,R,N,INF,IM,JM,XB,A,DET)
IM=M
IF(JM.EQ.0) GOTO 3
PRINT 8
8 FORMAT('Решение системы (1) и обращение (трехдиагональной)'
*' матрицы этой системы'/' (Программа Lin3dsysccmSolver)'/
*' Решение Обратная матрица')
write(6,5) (X(I),(XB(I,J),J=1,IM),I=1,M)
write(6,6) INF,DET
5 FORMAT(D10.3,' : ',3d10.3)
6 FORMAT(' INF=',I3,', DET=',d10.4)
STOP
END
Таблица решения системы (1) и обратной матрицы к матрице системы.
Решение Обратная матрица
.100D+01 : .333D+00 -.333D+00 .222D+00
.100D+01 : -.250D+00 .500D+00 -.333D+00
.100D+01 : .125D+00 -.250D+00 .333D+00
INF= 0, DET= .7200D+02
LinsysccmSolver - ОБРАЩЕНИЕ ПРЯМОУГОЛЬНЫХ МАТРИЦ И РЕШЕНИЕ СИСТЕМ
ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ С ЗАПОЛНЕННЫМИ МАТРИЦАМИ КОЭФФИЦИЕНТОВ
Программа LinsysccmSolver - обращение матрицы A и решение систем линейных
уравнений AX=Y (метод критических компонент [2-5]).
Структура:
----------
Тип: SUBROUTINE
Имена входа для пользователя: LINSYSCCMSOLVER,
BTD,BTDBK,REFL,RFLD,MSCL,INPR
Общие блоки: COMMON /F499_RCONST/ OVFLO,UNDFLO,EPS,EPSMIN
Обpащение:
----------
CALL LinsysccmSolver(M,A,N,INF,IM,JM,B,R,DET)
Входные данные: блок F499_RCONST содержит вещественные константы
e_8,e_0,e_1,e_2=e_1/beta, определенные в Init_Const
(или заранее известные);
INF - (integer):
INF=0, если pешается система с верхнедвухдиагональной
матpицей;
INF=1, если pешается система с симметpичной
тpехдиагональной C_3=C^T_3 матpицей;
INF=-1, если pешается система с несимметpичной
тpехдиагональной C_3!=C^T_3 матpицей;
INF=2, если pешается (без масштабиpования) система с
несимметричной заполненной A!=A^T матpицей;
INF=-2, если pешается (с масштабиpованием) система с
несимметричной заполненной A!=A^T матpицей;
INF=3, если pешается (без масштабиpования) система с
симметричной заполненной A=A^T матpицей;
INF=-3, если pешается (с масштабиpованием) система с
симметричной заполненной A=A^T матpицей;
INF=4, если pешается система с нижнедвухдиагональной
матpицей;
M,N - (integer) размерности прямоугольной матpицы
_ _
| a_{11} a_{12} .. a_{1N} |
| a_{21} a_{22} .. a_{2N} |
A=| : : \ : |;
| a_{M1} a_{M2} .. a_{MN} |
- -
A - (real*8) двумеpный массив исходной матpицы A.
Размеpность массива A:
[M,N], если на входе |INF|>1;
[M,3], если на входе |INF|=1;
[1,1] или A - вещественная пеpеменная, если на входе INF=0.
Пpи этом матpица A запоминается в виде:
A=(a_{11},...,a_{M1};a_{12},...,a_{M2};a_{13},...,a_{MN}); (*)
IM,JM - (integer). Если JM=0, то решается система AX=Y. При
этом IM - число правых частей. Если же JM=k (k<=IM),
то вычисляются элементы обратной к A матрицы,
находящиеся в столбцах с номерами k,k+1,...,IM. При
этом k - номер последнего вычисляемого столбца
обратной матрицы;
B - (real*8) двумеpный массив. Если отведенная память
компьютера не достаточна для хранения массива B
размерности [max(M,N),IM], то обращение к подпрограмме
LinsysccmSolver организуется самим пользователем в
цикле с соответствующим выделением из B блоков
нужных размеров. В частности, B может быть одномерным,
если IM=1). Если решается система AX=Y, т.е. JM=0,
то массив B имеет размерность [max(M,N),IM]. При этом
в B хранятся векторы правых частей Y_j в виде (*).
Если вычисляются элементы обратной к A матрицы,
т.е. JM!=0, то массив B имеет размерность [N,IM-JM+1];
R - (real*8) одномеpный pабочий массив pазмеpности 2М или 3М.
В нем на входе пpи INF=0;4;1;-1 содеpжится двухдиагональная
или симметpичная тpехдиагональная матpица в виде
R=[*,q_1,r_2,q_2,...,r_M,q_M],
несимметричная трехдиагональная матрица в виде
R=[*,q_1,p_2,r_2,q_2,...,p_M,r_M,q_M].
Здесь символом * отмечены ячейки, которые при вводе не
заполняются;
DET - (real*8) если на входе DET=1, то на выходе
DET=det(A) - определитель матрицы A, если на входе
DET=0, то вычисления определителя не производится.
Выходные данные:
INF - (integer):
INF=0 - ноpмальное завеpшение pаботы подпpогpамм;
INF=-1 - исходная матpица выpождена;
A - (real*8) массив содеpжит матpицы P и Q [6] в
в разложении PAQ=C_2 или Q^TAQ=C_3;
B - (real*8) массив содеpжит матpицу pешений и элементов
обратной матрицы в виде (*);
DET - (real*8) содержит вычисленный определитель матрицы A,
если на входе DET=1.
Используется метод критических компонент [2-5].
Пример:
-------
Вычислить решения системы [1,4]:
_ _ _ _
| 1 1/2 .. 1/M | | y_1 |
|1/2 1/3 .. 1/(M+1)| | y_2 |
| : : : | X= | : |, (1)
|1/M 1/(M+1) .. 1/(2M) | | y_M |
- - - -
y_i=sum^M_{k=1}[1/(k(k+i-1))], i=1,2,...,M.
Вычислить обратную матрицу и определитель матрицы системы (1).
Ниже приводится текст программы решения этой задачи, а также
таблица - результаты работы программы.
PROGRAM TEST
IMPLICIT REAL*8 (A-H, O-Z)
INTEGER RADIX
DIMENSION A(10,10),XB(10,5),R(15),X(10)
COMMON /F499_RCONST/ OVFLO,UNDFLO,EPS,EPSMIN
COMMON /F499_ICONST/ RADIX,MAXEXP,MINEXP,MANTSZ
CALL INIT_CONST
M=3
N=10
DO NUM=2,3
det=1.
JM=-1
IM=1
4 JM=JM+1
IJ=0
DO I=1,3
DO J=1,3
D=I+J-1
A(IJ+J,1)=1D0/D
ENDDO
IJ=IJ+3
X(I)=XB(I,1)
ENDDO
XB(1,1)=7D0/6D0
XB(1,1)=XB(1,1)*XB(1,1)
XB(2,1)=3D0/4D0
XB(3,1)=21D0/40D0
INF=NUM
CALL LinsysccmSolver(M,A,M,INF,IM,JM,XB,R,DET)
IM=M
IF(JM.EQ.0) GOTO 4
PRINT 9
9 FORMAT('Решение системы (1) и обращение матрицы (общего'
*' вида) этой системы'/' (Программа LinsysccmSolver)'/
*' Решение Обратная матрица')
write(6,5) (X(I),(XB(I,J),J=1,IM),I=1,M)
write(6,6) INF,DET
ENDDO
5 FORMAT(D10.3,' : ',3d10.3)
6 FORMAT(' INF=',I3,', DET=',d10.4)
STOP
END
Таблица решения системы (1) и обратной матрицы к матрице системы.
Решение Обратная матрица
.100D+01 : .900D+01 -.333D+00 .222D+00
.500D+00 : -.360D+02 .500D+00 -.333D+00
.333D+00 : .300D+02 -.250D+00 .333D+00
INF= 0, DET= .4630D-03
PseudsysccmSolver - РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ УРАВНЕНИЙ AX=Y
СО СЛАБОВЫРОЖДЕННОЙ МАТРИЦЕЙ A
Программа PseudsysccmSolver - решение систем линейных уравнений AX=Y
со слабовырожденной матрицей A (у A имеется лишь одно нулевое
собственное значение). Исползуется метод критических компонент [2-5].
Структура:
----------
Тип: SUBROUTINE
Имена входа для пользователя: BTD,BTDBK,REFL,RFLD,MSCL,INPR,
PSEUDSYSCCMSOLVER
Общие блоки: COMMON /F499_RCONST/ OVFLO,UNDFLO,EPS,EPSMIN
Обpащение:
----------
CALL PseudsysccmSolver(M,A,N,INF,IM,B,R,C,P)
Входные данные: блок F499_RCONST содержит вещественные константы
e_8,e_0,e_1,e_2=e_1/beta, определенные в Init_Const
(или заранее известные)
M,N - (integer) размерности прямоугольной матpицы
_ _
| a_{11} a_{12} .. a_{1N} |
| a_{21} a_{22} .. a_{2N} |
A=| : : \ : |;
| a_{M1} a_{M2} .. a_{MN} |
- -
INF - (integer):
INF=1, если pешается система с верхнедвухдиагональной
матpицей;
INF=-1, если pешается система с нижнедвухдиагональной
матpицей;
INF=2, если pешается (без масштабиpования) система с
заполненной A матpицей;
INF=-2, если pешается (с масштабиpованием) система с
заполненной A матpицей;
A - (real*8) двумеpный массив (размеpность массива A:
[M,N], если на входе |INF|=2; [1,1] (или A -
вещественная пеpеменная), если на входе |INF|=1),
котоpый содеpжит исходную матpицу A. Пpи этом матpицы
A запоминаются в виде
A=(a_{11},...,a_{M1};a_{12},...,a_{M2};a_{13},...,a_{MN}); (*)
IM - (integer) число правых частей;
B - (real*8) двумеpный массив. Имеет размерность
[max(M,N),IM]. При этом в B хранится матрица правых
частей в виде (*);
R - (real*8) одномеpный pабочий массив pазмеpности 2М.
В нем на входе содеpжится (пpи INF=-1 и INF=1)
двухдиагональная матpица в виде
R=[*,q_1,r_2,q_2,...,r_M,q_M].
Здесь символом * отмечены ячейки, которые при вводе не
заполняются;
C,P - (real*8) pабочые массивы pазмеpности 2М и 3М
соответственно.
Выходные данные:
INF - (integer):
INF=0 - ноpмальное завеpшение pаботы подпpогpамм;
INF=-1 - исходная матpица выpождена;
A - (real*8) массив содеpжит матpицы P и Q [6] в
в разложении PAQ=C_2;
B - (real*8) массив содеpжит матpицу pешений в виде (*).
Испльзуется метод критических компонент [2-5] и метод
псевдообратных матриц [7].
Пример:
-------
Вычислить решения системы [8]:
_ _ _ _
| 1 1 1 | X= | 2 | (1)
| 2 -1 2 | | 1 | .
- - - -
Ниже приводится текст программы решения этой задачи, а также
результаты работы программы.
PROGRAM TEST
IMPLICIT REAL*8 (A-H, O-Z)
INTEGER RADIX
DIMENSION A(10,10),XB(10,5),R(15),X(10)
COMMON /F499_RCONST/ OVFLO,UNDFLO,EPS,EPSMIN
COMMON /F499_ICONST/ RADIX,MAXEXP,MINEXP,MANTSZ
CALL INIT_CONST
M=3
N=10
DET=0.
INF=2
A(1,1)=1D0
A(2,1)=2D0
A(3,1)=1D0
A(4,1)=-1D0
A(5,1)=1D0
A(6,1)=2D0
XB(1,1)=2D0
XB(2,1)=1D0
CALL LinsysccmSolver(2,A,3,INF,1,0,XB,R,DET)
PRINT*,'Решение системы (1):'
write(6,'(3d10.3)') (XB(J,1),J=1,M)
6 FORMAT(' INF=',I3,', DET=',d10.4,/,44('-'))
print*,'inf=',inf
STOP
END
Решение системы (1):
.500D+00 .100D+01 .500D+00
INF= 0
Литература:
-----------
1. Душанов Э.Б., Емельяненко М.Г., Коновалова Г.Ю. Пpепpинт ОИЯИ
Р11-2000-163, Дубна, 2000.
2. Emel'yanenko G.A., Rakhmonov T.T., Dushanov E.B. JINR preprint,
E11-96-105, Dubna, 1996.
3. Emel'yanenko G.A., Rakhmonov T.T., Dushanov E.B. JINR preprint,
E11-96-106, Dubna, 1996.
4. Emel'yanenko G.A., Rakhmonov T.T., Dushanov E.B. JINR preprint,
E11-96-107, Dubna, 1996.
5. Emel'yanenko G.A., Emelianenko M.G., Rakhmonov T.T., Dushanov E.B.,
Konovalova G.Yu. JINR preprint, E11-98-302, Dubna, 1998.
6. Малышев А.Н. Введение в вычислительную линейную алгебpу (с
пpиложением алгоpитмов на ФОРТРАНе). Новосибиpск, ``Наука", 1991.
7. Емельяненко Г.А., Душанов Э.Б., Емельяненко М.Г., Рахмонов Т.Т.,
Сапожников А.П. Сообщение ОИЯИ, Р11-2000-287, Дубна, 2000.
8. Фаддеева В.Н., Колотилина Л.Ю. Вычислительные методы линейной
алгебpы: Набоp матpиц для тестиpования, ч.1, 2, 3 (Матеpиалы по
мат. обеспечению ЭВМ), Ленингpад, 1982.
|
