Текст подпрограммы и версий
ammhr_p.zip , ammhe_p.zip , ammhc_p.zip
Тексты тестовых примеров
tammhr_p.zip , tammhe_p.zip , tammhc_p.zip

Подпрограмма:  AMMHR (модуль AMMHR_p)

Назначение

Умножение матрицы, обратной к заданной в компактной форме вещественной ленточной матрице A, на прямоугольную матрицу и оценка числа обусловленности матрицы A.

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

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

      RCOND  =  1 / ( || A ||1 - || A-1 ||1 )  ,   где
                                           N
        || A ||1 =      max     {   ∑   | ai j |   }
                         j=1,...N       i=1 

и затем для заданной вещественной прямоугольной матрицы B размера N на NN вычисляется прямоугольная матрица С = А- 1 В размера N на NN путем решения NN систем линейных алгебраических уравнений А*С (J) = В(J), где С(J) и В(J) суть J - е столбцы матриц C и B,  J = 1,...,NN.

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

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

procedure AMMHR(var A :Array of Real; var MA :Integer; var N :Integer;
                var ML :Integer; var MU :Integer;
                var NLEAD :Array of Integer; var B :Array of Real;
                MB :Integer; NN :Integer; var RCOND :Real;
                var Z :Array of Real; var IERR :Integer); 

Параметры

A - вещественный двумерный массив размера МА*N, в первых МL+МU+1 столбцах которого задается в компактном виде ленточная матрица А порядка N; на выходе в первых МL столбцах массива находятся нижние кодиагонали ленточной матрицы L1 *...* LN - 1 (Li, i = 1,...,N - 1, суть элементарные матрицы исключения метода Гаусса), в следующих МL+МU+1 столбцах содержится в компактном виде матрица U;
MA - первая размерность массива А в вызывающей программе (тип: целый);
N - порядок матрицы А и число строк матриц B и C (тип: целый);
ML - число нижних кодиагоналей матрицы А (тип: целый);
MU - число верхних кодиагоналей матрицы А (тип: целый);
NLEAE - целый вектор длины N, содержащий на выходе информацию о выполненных в ходе факторизации перестановках (см. замечания по использованию);
B - вещественный двумерный массив размера МВ*NN, в котором задается матрица B; на выходе на соответствующих местах находятся элементы матрицы С = A- 1 B;
MB - первая размерность массива B в вызывающей программе (тип: целый);
NN - число столбцов матрицы B и C (тип: целый);
RCONE - вещественная переменная, содержащая на выходе вычисленное значение величины 1 / (|| A ||1 - || A- 1 ||1) (см. замечания по использованию);
Z - вещественный рабочий вектор длины N;
IERR - целая переменная, содержащая на выходе информацию о прохождении счета, при этом:
IЕRR=65 - если хотя бы одна из переменных MA, N, MB, NN имеет значение, меньшее единицы;
IЕRR=66 - если в процессе счета произошло переполнение (это говорит о том, что либо || A ||1, либо некоторые элементы матрицы U или матрицы C превосходят по абсолютной величине максимальное представимое на данной машине число);
IЕRR= -К - если в результате факторизации в К - ой строке матрицы U диагональный элемент равен нулю (это свидетельствует о вырождености матрицы А). Если таких строк у матрицы U несколько, то значение K полагается равным номеру последней из них;
IЕRR=67 - если для некоторого J, 1 ≤ J ≤ NN, система А*С (J) = В(J) несовместна.

Версии

AMMHE - умножение матрицы, обратной к заданной с расширенной (Extended) точностью в компактной форме вещественной ленточной матрице A, на вещественную прямоугольную матрицу, заданную с расширенной (Extended) точностью, и оценка числа обусловленности матрицы А;
AMMHC - умножение матрицы, обратной к заданной в компактной форме комплексной ленточной матрице A, на комплексную прямоугольную матрицу и оценка числа обусловленности матрицы А.

Вызываемые подпрограммы

AFB2R - подпрограмма треугольной факторизации и оценки числа обусловленности ленточной матрицы А.
UTAFSI - подпрограмма выдачи диагностических сообщений.

Замечания по использованию

  1. 

В подпрограмме АММНE массивы А, B, Z и переменная RСОNE имеют тип Extended, для треугольного разложения и оценки числа обусловленности ленточной матрицы А вызывается подпрограмма АFВ2E.

  2. 

В подпрограмме AMMHC массивы А, B и Z имеют тип Complex, для треугольного разложения и оценки числа обусловленности ленточной матрицы А вызывается подпрограмма АFВ2С.

  3. 

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

  4. 

При обращении к подпрограмме необходимо выполнение условия МВ ≥ N. Кроме того, так как в результате выполненных в ходе факторизации перестановок число верхних кодиагоналей матрицы U равно МU+МL, а также в силу некоторых конструктивных особенностей подпрограммы, для правильной ее работы необходимо выполнение условия МА ≥ N > МU + 2 МL + 1. Если же МU + 2 МL + 1 ≥ N, то более целесообразно, задав матрицу А не в компактной, а в полной форме, обратиться к подпрограмме АММGR.

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

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

Unit TAMMHR_p;
interface
uses
SysUtils, Math, { Delphi }
LStruct, Lfunc, UtRes_p, AMMHR_p;

function TAMMHR: String; 

implementation

function TAMMHR: String;
var
MA,N,ML,MU,MB,NN,I,J,J0,J1,K,IERR :Integer;
RCONE :Real;
A :Array [0..24] of Real;
B :Array [0..49] of Real;
Z :Array [0..4] of Real;
NLEAE :Array [0..4] of Integer;
label
_1,_2,_3,_4;
begin
Result := '';
МА := 5;
N := 5;
ML := 1;
MU := 1;
МВ := 5;
NN := 10;
for I:=1 to МА do
 begin
  for J:=1 to N do
   begin
    A[(I-1)+(J-1)*5] := 0.0;
    B[(I-1)+(J-1)*5] := 0.0;
    B[(I-1)+(J+5-1)*5] := 0.0;
_1:
   end;
_2:
 end;
for I:=1 to МА do
 begin
  J0 := Max0(1,I-ML);
  J1 := Min0(N,I+MU);
  for J:=J0 to J1 do
   begin
    K := J-I+ML+1;
    A[(I-1)+(K-1)*5] := (I*10+J);
    B[(I-1)+(J-1)*5] := A[(I-1)+(K-1)*5];
    B[(I-1)+(J+5-1)*5] := A[(I-1)+(K-1)*5];
_3:
   end;
_4:
 end;
Result := Result + Format('%s',['   A=']);
Result := Result + #$0D#$0A;
for J:=1 to N do
 begin
for I:=1 to МА do
 begin
  Result := Result + Format('%20.16f ',[A[(I-1)+(J-1)*5]]) + #$0D#$0A;
 end;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['   B=']);
Result := Result + #$0D#$0A;
for J:=1 to NN do
 begin
for I:=1 to МВ do
 begin
  Result := Result + Format('%20.16f ',[B[(I-1)+(J-1)*5]]) + #$0D#$0A;
 end;
 end;
Result := Result + #$0D#$0A;
AMMHR(A,MA,N,ML,MU,NLEAE,B,MB,NN,RCONE,Z,IERR);
Result := Result + Format('%s',['   NLEAE=']);
Result := Result + #$0D#$0A;
for J:=1 to N do
 begin
  Result := Result + Format('%3d ',[NLEAD[J-1]]) + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['   A=']);
Result := Result + #$0D#$0A;
for J:=1 to N do
 begin
for I:=1 to МА do
 begin
  Result := Result + Format('%20.16f ',[A[(I-1)+(J-1)*5]]) + #$0D#$0A;
 end;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['   IERR=']);
Result := Result + Format('%3d',[IERR]) + #$0D#$0A;
Result := Result + Format('%s',['     RCONE=']);
Result := Result + Format('%20.16f ',[RCOND]) + #$0D#$0A;
Result := Result + Format('%s',['   B=']);
Result := Result + #$0D#$0A;
for J:=1 to NN do
 begin
for I:=1 to МВ do
 begin
  Result := Result + Format('%20.16f ',[B[(I-1)+(J-1)*5]]) + #$0D#$0A;
 end;
 end;
Result := Result + #$0D#$0A;
UtRes('TAMMHR',Result);  { вывод результатов в файл TAMMHR.res }
exit;
end;

end.

Результаты:

            |   0         21.0     22.0   23.0    0 |
            | -0.524   32.0     33.0   34.0    0 |
   A  =  | -0.015   43.0     44.0   45.0    0 |
            |  0.292   54.0     55.0     0       0 |
            | -0.228    0.569    0        0       0 |

   NLEAD  =  (2, 3, 4, 5, 5) ,   RCOND  =  1.47362E-03 ,

            | 1.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0 |
            | 0.0  1.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0 |
   B  =  | 0.0  0.0  1.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0 |
            | 0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  1.0  0.0 |
            | 1.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  1.0 |