Текст подпрограммы и версий
afg4r_p.zip , afg4e_p.zip , afg4c_p.zip 		Тексты тестовых примеров
tafg4r_p.zip , tafg4e_p.zip , tafg4c_p.zip

Подпрограмма:  AFG4R (модуль AFG4R_p)

Назначение

Треугольное разложение вещественной матрицы общего вида методом Гаусса с выбором ведущего элемента по столбцу и оценка числа обусловленности матрицы.

Математическое описание

Для заданной квадратной вещественной матрицы А порядка N выполняется треугольная факторизация А = LU, где U - верхняя треугольная матрица, а L - нижняя треугольная матрица, и вычисляется величина RCOND, обратная числу обусловленности матрицы А: 

     RCOND = 1 / || A ||1 * || A-1 ||1

Для вычисления || A- 1 ||1 решается система АZ = Y и полагается 

     || A-1 ||1 ≈ || Z ||1 / || Y ||1 ,

где Y - решение системы АTY = Е; компоненты вектора Е выбираются таким образом, чтобы максимизировать норму || W ||1, где W - решение системы UTW = Е.

Вычисления проводятся в следующей последовательности: 

1) вычисляется  || A ||1 = max || aj ||1 ,      для   j = 1, ..., N ,
        где аj  -  j-й  вектор-столбец матрицы  A, а
       || аj ||1 = | a1 j | + | a2 j | +...+ | aN j | ;

2) выполняется треугольная факторизация А = LU ;

3) решается система UTW = E
    одновременно с выбором компонент вектора E;

4) решается система LTY = W;

5) решается система LV = Y;

6) решается система UZ = V и полагается 

         RCOND = || Y ||1 / ( || A ||1 * || Z ||1 )

Дж. Форсайт, М. Малькольм, К. Моулер. Машинные методы математических вычислений. М.: Мир, 1980.

Использование 

procedure AFG4R(var A :Array of Real; M :Integer; N :Integer;
                var NLEAD :Array of Integer; var RCOND :Real;
                var Z :Array of Real; var IERR :Integer);


   Параметры
A - вещественный двумерный массив размера М*N, в котором задается исходная квадратная матрица порядка N, на выходе на соответствующих местах находятся элементы матрицы U и поддиагональные элементы элементарных матриц исключения Li,  i = 1, ..., N - 1;
M - первая размерность массива А в вызывающей программе (тип: целый);
N - порядок матрицы А (тип: целый);
NLEAD - целый вектор длины N, содержащий на выходе информацию о выполненых в процессе факторизации перестановках (см. замечания по использованию);
RCOND - вещественная переменная, содержащая на выходе вычисленное значение величины, обратной числу обусловленности матрицы А (см. замечания по использованию);
Z - вещественный рабочий вектор длины N (см. замечания по использованию);
IERR - целая переменная, содержащая на выходе информацию о прохождении счета, при этом:
IЕRR=65 - если значение М ≤ 0 или N ≤ 0;
IЕRR=66 - если в процессе работы подпрограммы произошло переполнение (это говорит о том, что норма || A ||1 либо некоторые элементы матрицы U превосходят по абсолютной величине максимальное представимое на данной машине число);
IЕRR=-К - если в результате факторизации диагональный элемент в К - й строке матрицы U равен нулю (это свидетельствует о вырожденности матрицы А). Если таких строк у матрицы несколько, то значение К полагается равным номеру последней из них. 
   Версии
AFG4E - треугольное разложение вещественной матрицы, заданной с расширенной (Extended) точностью, методом Гаусса с выбором ведущего элемента по столбцу и оценка числа обусловленности матрицы.
AFG4C - треугольное разложение комплексной матрицы методом Гаусса с выбором ведущего элемента по столбцу и оценка числа обусловленности матрицы. 
   Вызываемые подпрограммы
UTAFSI - подпрограмма выдачи диагностических сообщений. 
   Замечания по использованию
  1. 

В подпрограмме АFG4С массивы А и Z имеют тип Complex.

  2. 

В подпрограмме АFG4E массивы А, Z и переменная RСОND имеют тип Extended.

  3. 

На выходе К - й элемент вектора NLЕАE равен номеру строки, переставленной на К - м шаге факторизации с К - й строкой матрицы А. Поскольку факторизация Гаусса требует N - 1 шагов, то NLЕАD (N) = N.

  4. 

Если переменной IЕRR присвоено значение, отличное от нуля, то величина RСОNE полагается равной нулю, выдается соответствующее диагностическое сообщение и происходит выход из подпрограммы.

  5.  Если IЕRR = 0, то на выходе вектор Z удовлетворяет условию || AZ ||1 = RСОND * || A ||1 * || Z ||1 .

Пример использования 

Unit TAFG4R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, AFG4R_p;

function TAFG4R: String;

implementation

function TAFG4R: String;
var
M,N,J,I,JI,JJ,II,IERR :Integer;
RCONE :Real;
NLEAE :Array [0..3] of Integer;
Z :Array [0..3] of Real;
const
A :Array [0..15] of Real = ( 1.0,0.42,0.54,0.66,0.42,1.0,0.32,0.44,0.54,0.32,
1.0,0.22,0.66,0.44,0.22,1.0 );
begin
Result := '';  { результат функции }
M := 4;
N := M;
Result := Result + #$0D#$0A;
for I:=1 to M do
 begin
  for J:=1 to N do
   begin
    Result := Result + Format('%20.16f ',[A[(I-1)+(J-1)*4]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
AFG4R(A,M,N,NLEAE,RCONE,Z,IERR);
Result := Result + #$0D#$0A;
for JI:=1 to N do
 begin
  Result := Result + Format('   NLEAE= + #$0D#$0A   %3d ',
 [NLEAD[JI-1]]) + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
for II:=1 to M do
 begin
  for JJ:=1 to N do
   begin
    Result := Result + Format('%20.16f ',
 [A[(II-1)+(JJ-1)*4]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['   IERR=']);
Result := Result + Format('%3d ',[IERR]);
Result := Result + Format('%s',['     RCONE=']);
Result := Result + Format('%20.16f ',[RCOND]) + #$0D#$0A;
UtRes('TAFG4R',Result);  { вывод результатов в файл TAFG4R.res }
exit;
end;

end.

 Для данной матрицы известны собственные значения:
  
      λ1  =  2.32274
      λ2  =  0.24226
      λ3  =  0.63828
      λ4  =  0.79670

 Спектральное  число  обусловленности   λ1 / λ2 ≈ 9.588
  

Результаты:

               |   1.0       0.42          0.54           0.66      |
               | -0.42     0.82360     0.09320     0.16280 |
      A  =  | -0.54    -0.11316     0.69785   -0.15482 |
               | -0.66    -0.19767     0.22186     0.49787 |

      NLEAD  =   (1, 2, 3, 4)
      RCOND  =  0.09880