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

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

Итак, решается система A*X = C размерностью N,
где A - матрица [ N : N ], X,C - вектора длиной N.
На практике очень часто решают целое семейство из 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 совмещаются в одно обращение, когда размерность системы задается со знаком "-" :

При этом Solver, произведя факторизацию исходной матрицы А, сразу же начнет решать систему A * X = E с единичной матрицей правых частей.

Перед первым обращением к Solver-у можно подкорректировать некоторые его глобальные параметры, обратившись к

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

---