ОБЪЕДИНЕННЫЙ   ИНСТИТУТ   ЯДЕРНЫХ   ИССЛЕДОВАНИЙ
lit
БИБЛИОТЕКА   ПРОГРАММ   JINRLIB

JINRLINPACK - машинно-независимый пакет программ для обращения двухдиагональных, трехдиагональных и прямоугольных (заполненных) матриц и решение систем линейных уравнений с соответствующими матрицами коэффициентов

F499

 Авторы: Г.А.Емельяненко, Э.Б.Душанов, М.Г.Емельяненко,
Т.Т.Рахмонов, А.П.Сапожников		 Язык: Фортран

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

Обращение:

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) содержит выделенный порядок числа.

Подпрограмма-функция 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

Обращение:

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

Обращение:

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_CONST/ OVFLO,UNDFLO,EPS,EPSMIN

Обращение:

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_CONST/ OVFLO,UNDFLO,EPS,EPSMIN

Обращение:

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_CONST/ OVFLO,UNDFLO,EPS,EPSMIN

Обращение:

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.

home up e-mail