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.