Выходные параметры: A_LÎR^n,m1- массив, содержащий нижнюю треугольную матрицу, полученную при разложении матрицы А; intÎZⁿ – целочисленный массив, содержащий информацию о перестановке строк, выполненных при приведении матрицы А к треугольной форме.

Код возврата: retcode=0 при успешном решении и retcode = 1 в противном случае.

Описание: Алгоритм реализует разложение матрицы с частичным выбором ведущего элемента с учетом ленточности матрицы. Для общности будем предполагать, что рассматриваемая ленточная матрица несимметрична, т.е. ненулевыми являются только элемента, для которых j>i+m2 и i>j+m1. Специальный алгоритм для обработки таких матриц состоит из N-1 основных шагов. На текущем k-ом шаге алгоритма исключаются поддиагональные элементы k-столбца. На текущем шаге алгоритма, кроме m1 последних шагов, исключается m1 ненулевых поддиагональных элементов. Текущий k-ый шаг алгоритма включает в себя следующие операции:

1. Выбор главного элемента по столбцу среди диагональных  и поддиагональных элементов. Если несколько поддиагональных элементов оказались равными, то у главного элемента номер строки наименьший. Номер этой строки записан в k-м элементе массива Int, так что главный элемент k-го столбца удовлетворяет условию

;

2. перестановку строк с номерами k и int[k],
3. вычисление величин

M[i,k]=A[i,k]/A[k,k],|M[i,k]|<=1,

после чего из каждой строки с номером  , где 1=min(n,k+mi), вычитается k-я строка, умноженная на коэффициент M[i,k].

Необходимая экономия памяти при размещении ленточных  матриц достигается с помощью  следующего приема.

Память: В алгоритме  матрица А представлена в виде массива А размера n x (m1+m2+1). При такой форме записи не удается использовать свойства симметрии, если матрица А ими обладает. Для удобства столбцы нового массива нумеруются от -m1 до m2, причем диагональные элементы исходной матрицы А расположены в столбце с номером 0. Расположение элементов в массиве А для случая n=7, m1=2, m2=1 показано ниже.

                  

Заметим, что  в первую, вторую и последнюю строки массива А требуется ввести нулевые  элементы.

В процессе разложения матрицы массив А перестраивается  таким образом, чтобы элементы, участвующие в текущем преобразовании, образовывали столбец. Соответствующее расположение элементов в массиве А непосредственно перед исключением переменных х1 и х2 показано ниже.

    

Строки, участвующие  в текущем преобразовании, выделены прямоугольником. Изменение массива  А производится на всех шагах, за исключением  первого.

Алгоритм:

1. retcode<=0
2. CALL BANDECOMP (n,m1,m2,macheps,A,A_L,int,retcode)

{Алг. А1.5.1}

3. IF retcode=0 THEN {разложение успешно}

CALL BANSLU(n,m1,m2,A,A_L,int,b) {Алг. А1.5.2}

{RETURN для алг. A1.5)

{END для Алг. A1.5)

 

1.2.2. Алгоритм А1.5.1 (BANDECOMP). Триангуляция ленточной матрицы

Назначение: Находит разложение ленточной матрицы в виде произведения нижней L и верхней U треугольных матриц.

Входные параметры: nÎZ, m1ÎZ, m2ÎZ, machepsÎR.

bÎRⁿ , L и U в массиве АÎR^n,n; intÎZⁿ, A0ÎRⁿ.

Входно-выходные параметры: АÎR^n,m1+m2+1 – массив, состоящий из элементов ленточной матрицы (после разложения содержит нижний треугольный сомножитель L).

Выходные параметры: АÎR^n,m1 – массив, содержащий нижнюю треугольную матрицу L, полученную при разложении матрицы А; intÎZⁿ – целочисленный массив, содержащий информацию о перестановке строк, выполненных при приведении матрицы А к треугольной форме.

Код возврата: : retcode=0 при успешном решении и retcode = 1 – если матрица вырождена.

Память: Исходная ленточная матрица А порядка n записана в виде массива A[1..n,-m1..m2] размера n x (m1+m2+1). Полученная из нее матрица разлагается на произведение нижней и верхней треугольных матриц с использованием перестановки строк. Ненулевые элементы нижней треугольной матрицы хранятся в массиве A_L[1..n,1..m1], а ненулевые элементы верхней треугольной матрицы записываются в массив А. Массив Int[1..n] содержит информацию о перестановках строк при разложении матрицы А.

Алгоритм:

{подготовка к первому шагу}

1.l <=m+1

2. FOR i=1 TO m1 DO

      2.1 FOR  j=l-1 TO m2 DO

      A[i,j-1]<=A[i,j]

      2.2. l <= l-1

      2.3. FOR j=m2-l TO m2 DO

      A[i,j]<=0

{Гауссово исключение с частичным выбором ведущего элемента}

3. l<=m1

FOR k=1 TO n DO

4.1 x<=A[k,-m1]

4.2 i<=k

4.3 IF l<k THEN

l<=l+1

4.4 FOR j=k+1 TO l DO

     4.4.1 IF |A[j,-m1]|>|x| THEN

           4.4.1.1 x<=|A[j,-m1]|

           4.4.1.2 i<=j

4.5 int[k]<=i

{проверка на вырожденность}

4.6 IF x<macheps THEN

      4.6.1 retcode<=1

      4.6.2 RETURN

4.7 IF i<>k THEN

{перестановка строк i и k}

      4.7.1 FOR j=-m1 TO m2 DO

      переставить местами A[k,j] и A[i,j]

4.8 FOR i=k+1 TO l DO

      4.8.1 x<=A[i,-m1]/A[k,-m1]

      4.8.2 A_L[k,i-k]<=x

      4.8.3 FOR j=1-m1 TO m2 DO

      A[i,j-1]<=A[i,j]-x*A[k,j]

      4.8.4 A[i,m1]<=0

{RETURN для Алг. A1.5.1}

{END для Алг. A1.5.1}

 

1.2.3 Алгоритм А1.5.2 (BANSLU).

Решение СЛУ  Ax=b относительно х с ранее триангулированной ленточной матрицы A=LU

Назначение: Ищется решение СЛУ Ax=b, где матрица предварительно разложена на треугольные с помощью алгоритма А1.5.1.

Входные параметры: nÎZ, m1ÎZ, m2ÎZ, АÎR^n,m1+m2+1 – массив, содержащий верхнюю треугольную матрицу U разложения для исходной матрицы А; A_L ÎR^n,m1 – массив, содержащий нижнюю треугольную матрицу L, полученную при разложении матрицы А; intÎZⁿ.

Входно-выходные параметры: вектор правых частей bÎRⁿ (после прямой и обратной подстановок содержит решение исходной системы).

Алгоритм:

1. i<=m1

{прямая подстановка}

2. FOR k=1 TO n DO

      2.1 i<=int[k]

      2.2 IF i<>k THEN

      l<=l+1

      2.4 FOR i=k+1 TO l DO

      b[i]<=b[i]-A_L[k,i-k]*b[k]

{обратная подстановка}

3. l<=-m1

4. FOR i=n DOWNTO 1 TO

      4.1

      4.2 IF l<m2 THEN

      l<=l+1

{RETURN для Алг. A1.5.2}

{END для Алг. A1.5.2}

 

2. Программная реализация

2.1. Тексты программ

program Kraut;

uses crt,dos;

const nn=40;

mm1=20;

mm2=20;

type

   matr=array[1..nn,-mm1..mm2] of real;

   matr0=array[1..nn,1..nn] of real;

   vect=array[1..nn] of real;

   vect0=array[1..nn] of integer;

var

   a,a_l:matr;

   a0: matr0;

   x,y,b:vect;

   int: vect0;

   retcode,i,j,k,N,m,m1,m2:integer;

   macheps,max,delta:real;

   hour,minute,second,sot:word;

   hour0,minute0,second0,sot0:word; 

{пpоцедуpа генеpации матpицы}

procedure GenMatr(var a0:matr0;var b:vect;m1,m2,n: integer);

var

   p: real;

   kod: integer;

function Urand (kod: integer):real;

label 10,20,99;

const itwo=2;

      ia:integer=0;

      ic:integer=0;

      mic:integer=0;

      iy:integer=0;

      m2:integer=0;

      s:real=0;

var

      m:integer;

      halfm:real;

begin

   if kod<0 then begin

   m2:=0;

   goto 99;

   end;

   if m2<>0 then goto 20;

   m:=1;

   iy:=kod;

   10:m2:=m;

   m:=itwo*m2;

   if m>m2 then goto 10;

   halfm:=m2;

   ia:=trunc (halfm*arctan(1.0)/8.0)+5;

   ic:=trunc (halfm*(0.5-sqrt(3.0)/6.0)+1);

   mic:=(m2-ic)+m2;

   s:=0.5/halfm;

   20:iy:=iy*ia;

   if iy>mic then iy:=(iy-m2)-m2;

   iy:=iy+ic;

   if iy/2>m2 then iy:=(iy-m2)-m2;

   if iy=0 then iy:=(iy+m2)+m2;

   99: urand:=iy*s;

end;

begin

  write ('Код генерации матрицы - ');

  readln (kod);

  p:=urand(-1); 

{фоpмиpование матpицы  A}

       for i:=1 to n do begin

         for j:=1 to n do begin

             if (j>i+m2) or (i>j+m1) then a0[i,j]:=0

             else a0[i,j]:=urand(kod);

             write (a0[i,j]:6:2);

         end;

         writeln;

         end;

{фоpмиpование вектоpа b}

  for i:=1 to n do begin

    b[i]:=0;

    for j:=1 to n do

      b[i]:=b[i]+a0[i,j];

  end; 

  write ('>...');

  readkey;

end; 

{Пpоцедуpа генеpации матpицы 4}

procedure Matrix4(var a0: matr0;var b:vect;m1,m2,n:integer);

var alfa: real;

begin

  clrscr;

  alfa:=5;

{Данный алгоритм  производит формирование матрицы}

  for i:=1 to n do begin

    for j:=1 to n do begin

             if (j>i+m2) or (i>j+m1) then a0[i,j]:=0

             else  if i=j then a0[i,j]:=alfa

             else a0[i,j]:=1; 

             write (a0[i,j]:6:2);

    end;

    writeln;

    end;

{Формирование  вектора b как суммы элементов  по строкам матрицы}

  for i:=1 to n do

   begin

    b[i]:=0;

    for j:=1 to n do

        b[i]:=b[i]+a0[i,j];

   end;

  write ('>...');

  readkey;