XLLS_Solver         Библиотека "JINRLIB"               F498      

    Авторы: П.Г.Акишин,А.П.Сапожников
    Язык: Фортран

            РЕШЕНИЕ СВЕРХБОЛЬШИХ СИСТЕМ ЛИНЕЙНЫХ УРАВНЕНИЙ
        С ИСПОЛЬЗОВАНИЕМ  ВНЕШНЕЙ  РАБОЧЕЙ  ПАМЯТИ  НА ДИСКАХ

    Как это ни странно, проблема использования дисковой памяти при
    решении больших систем линейных уравнений все еще остается актуальной.
    В самом деле, уже для матрицы 2000x2000 требуется 32 мегабайта памяти,
    что доступно далеко не на каждой "персоналке". На больших же компьютерах
    бездумное сканирование матрицы "по столбцам и по строкам" порождает
    огромное количество замещений страниц виртуальной памяти, что резко
    замедляет работу.
    Предлагаемая программа, используя классическую схему Гаусса с выбором
    главного элемента по столбцу, хранит матрицу на диске и обрабатывает ее
    порциями в буфере, заданном пользователем.  В частном случае, когда
    буфер настолько велик, что вмещает всю матрицу, работа происходит без
    обращения к диску.  Способ размещения матрицы и доступа к ней полностью
    скрыт от пользователя, от него требуются только две вещи:

    1. заготовить подпрограмму-функцию, вычисляющую элементы
       исходной матрицы и правую часть системы;
    2. заготовить массив для буферизации матрицы, чем больше он будет,
       тем меньше будет обращений к диску.

    Итак, решается система    A*X = C   размерностью  N,
    где A - матрица [ N : N ],   X,C - вектора длиной N.
    На практике очень часто решают целое семейство из  NC  систем

       A * X(k) = C(k)          k=1,2,...,NC

    с одинаковой матрицей A. Вначале матрица факторизуется, превращаясь
    в две треугольные матрицы  L и R, такие, что  A = ~L * R .
    После этого для любой правой части  C(k) системы  решение X(k)
    получается тривиальным перемножением:  X(k) = ~R * ( L * C(k) ).

    Работа пользователя происходит в 3 стадии:

    1. CALL XLLS_Solver(buf,lbuf,n,elm)  - прием матрицы и ее факторизация.
       n    - размерность системы;
       buf  - одномерный массив, описанный в вызывающей программе как REAL*8
              длиной не менее lbuf слов;
       lbuf - число слов массива buf, разрешенное к использованию в качестве
              рабочей памяти Solver-a.  lbuf не должно быть меньше чем 3*n.
              С другой стороны, для работы без участия дисков программе
              требуется не более чем 2*n*(n+nc+1)+n слов.
              При работе с дисками используются логические номера 78 и 79;
       elm(i,j) - (REAL*8) функция, вычисляющая элемент А(i,j) исходной матрицы.
              Solver запрашивает исходную матрицу А, обращаясь к функции elm
              строго "по столбцам", т.е. вот в такой последовательности:
              (1,1),(2,1),...,(n,1), (1,2),(2,2),...,(n,2),...,(n-1,n),(n,n).

    2. CALL XLLS_Solver(buf,0,nc,elc)  - решение системы с заданным комплектом
                                         ее правых частей.
       buf  - тот же самый рабочий массив, что и в стадии 1;
       nc   - количество правых частей системы;
       elc(i,j) - (REAL*8) функция, вычисляющая элемент C(i,j) матрицы C
              правых частей системы (n строк, nc столбцов).
              Обращения к elc также производятся сугубо "по столбцам".
              Полученное решение системы остается на дисках.

    3. CALL XLLS_Extractor(buf,x,k)  - извлечение решения X(k) системы,
                                       соответствующего правой части C(k).
       buf  - тот же самый рабочий массив, что и в стадии 1;
       k    - 1,2,...,nc - номер правой части системы;
       x    - (REAL*8) массив пользователя длиной n. 
              Сюда будет записано k-e решение.

    Если требуется не решать систему, а получить матрицу ~А, обратную
    исходной, стадии 1 и 2 совмещаются в одно обращение, когда размерность
    системы задается со знаком "-" :
       CALL XLLS_Solver(buf,lbuf,-n,elm)
    При этом Solver, произведя факторизацию исходной матрицы А, сразу же
    начнет решать систему  A * X = E  с единичной матрицей правых частей.

    Перед первым обращением к Solver-у можно подкорректировать 
    некоторые его глобальные параметры, обратившись к 
       subroutine XLLS_Params(np), где:
       np = 1 - установка режима экономного расходования дисковой памяти.
                При этом на диск будут писаться только старшие, наиболее
                значимые, биты всех матричных элементов. Это позволит
                сократить требуемый обьем диска c 16*n*(n+nc) вдвое,
                хотя несколько ухудшит точность вычислений.
       np = 2 - восстановление стандартного режима использования диска.
       np = 7 - печать счетчика дисковых обменов.
       np = 8 - на некоторых машинах, в частности на VAX, фортранная
                подсистема ввода-вывода требует задавать размер записи файла
                в словах, а не в байтах!