unit UMatrix;

interface
type TMatrix = array of array of real;
     TVector = array of real;
     TInt = array of integer;

// ShowErrorMatrix - надо ли показывать окна с сообщениями
// об ошибках
// Если ShowErrorMatrix = false, то ErrorMatrix показывает,
// были ли ошибки, а MessageMatrix содержит
// последнее сообщение об ошибке
var ShowErrorMatrix: boolean = false;
    ErrorMatrix: boolean = false;
    MessageMatrix: string;

// ******************************************
// Операции с векторами
// ******************************************

// Создание вектора размерностью N
function CreateVector(N: integer): TVector; overload;

// Создание вектора из массива Arr
function CreateVector(const Arr: array of real): TVector; overload;

// Возвращает размерность вектора
function DimVector(const V: TVector): word;

// Возвращает последний индекс вектора
function HighVector(const V: TVector): integer;

// Модуль вектора
function ModVector(const V: TVector): real;

// Эквивалентность двух векторов
function EqVector(const V1, V2: TVector): boolean;

// Сложение двух векторов
// Размерность векторов должна быть одинаковой
function SumVector(const V1, V2: TVector): TVector;

// Разность двух векторов
// Размерность векторов должна быть одинаковой
function SubVector(const V1, V2: TVector): TVector;

// Произведение скаляра (числа) на вектор
function CMultVector(C: real; const V: TVector): TVector;

// Скалярное умножение двух векторов
// Размерность векторов должна быть одинаковой
function MultSVector(const V1, V2: TVector): real;

// Косинус угла между двумя векторами
// Размерность векторов должна быть одинаковой
function CosVector(const V1, V2: TVector): real;

// Угол в радианах между двумя векторами
// Размерность векторов должна быть одинаковой
function FVector(const V1, V2: TVector): real;

// Векторное умножение двух трехмерных векторов
function MultVVector(const V1, V2:  array of real): TVector;

// ******************************************
// Операции с матрицами
// ******************************************

// Создание матрицы M x N: M - число строк, N - число столбцов
function CreateMatrix(M, N: integer): TMatrix; overload;

// Возвращает число столбцов N
function ColMatrix(const A: TMatrix): integer;

// Возвращает число строк M
function RowMatrix(const A: TMatrix): word;

// Возвращает последний индекс столбца
function HighColMatrix(const A: TMatrix): integer;

// Возвращает последний индекс строки
function HighRowMatrix(const A: TMatrix): integer;

// Транспонирование матрицы
function TransMatrix(const A: TMatrix): TMatrix;

// Эквивалентность двух матриц
function EqMatrix(const A, B: TMatrix): boolean;

// Сумма двух матриц
// Размерность матриц должна быть одинаковой
function SumMatrix(const A, B: TMatrix): TMatrix;

// Разность двух матриц
// Размерность матриц должна быть одинаковой
function SubMatrix(const A, B: TMatrix): TMatrix;

// Произведение скаляра (числа) на матрицу
function CMultMatrix(C: real; const A: TMatrix): TMatrix;

// Произведение двух матриц
// Размерности матриц должны быть M x K и K x N
function MultMatrix(const A, B: TMatrix): TMatrix;

// Произведение матрицы на вектор
// Размерности матрица должны быть M x N, а вектора N
function MultMatrixVector(const A: TMatrix; const V: TVector): TVector;

// Вычисление определителя квадратной матрицы
function DetMatrix(const A: TMatrix): real;

// Прямое вычисление обратной матрицы
// Матрица A должна быть квадратной
function ReversMatrix(const A: TMatrix): TMatrix;
function ReversMatrix2(const A: TMatrix): TMatrix;

// Решение системы линейных уравнений
// Матрица A должна быть квадратной
// Вектор VRight той же размерности, что и A
function LSystem(const A: TMatrix; const VRight: TVector): TVector;

// LU разложение матрицы
// Matr - заданная матрица, заменяемая в результате разложенной
// V - массив индексов, определяющих последовательность 
// уравнений с учетом перестановок в процессе LU-разложения
procedure LU(var Matr: TMatrix; var V: TInt);

// LU разложение матрицы с выбором главного элнмента
// Matr - заданная матрица, заменяемая в результате разложенной
// V - массив индексов, определяющих последовательность 
// уравнений с учетом перестановок в процессе LU-разложения
procedure LUMax(var Matr: TMatrix; var V: TInt);

// Решение системы линейных уравнений
// с уже проведенным LU-разложением (реультаты в Matr и V)
// и с изменяющимся вектором правых частей B
function LUVector(const B: TVector;
                  const Matr: TMatrix; const V: TInt): TVector;

// Полное решение системы линейных уравнений
// с помощью LU-разложения без выбора главного элемента
function LUSystem(const A: TMatrix; const VRight: TVector):
                                                       TVector;

// Полное решение системы линейных уравнений
// с помощью LU-разложения с выбором главного элемента
function LUMaxSystem(const A: TMatrix; const VRight: TVector):
                                                       TVector;

// ****************************************************
implementation

uses Forms,   // Для переменной Application
     Windows, // Для обозначений флагов кнопок Application.MessageBox
     SysUtils, Math;// Для функции Abort

procedure Messages(S:string);
begin
 if ShowErrorMatrix then begin
   Application.MessageBox(PChar(S),'Ошибка',MB_OK+MB_ICONSTOP);
   Abort;
 end
 else begin
  ErrorMatrix := true;
  MessageMatrix := S;
 end;
end;

function CreateVector(N: integer): TVector;
// Создание вектора размерностью N
begin
 SetLength(Result, N);
end;

function CreateVector(const Arr: array of real): TVector;
// Создание вектора из массива Arr
var i, N:word;
begin
ErrorMatrix := false;
 if (Length(Arr) = 0)
 then begin
   Messages('Функция CreateVector :' + #13'Задан нулевой параметр');
   exit;
 end;
 N := High(Arr);
 SetLength(Result, N+1);
 for i:=0 to N do
   Result[i] := Arr[i];
end;

function DimVector(const V: TVector): word;
// Возвращает размерность вектора
begin
 Result := Length(V);
end;

function HighVector(const V: TVector): integer;
// Возвращает последний индекс вектора
// Если размерность вектора 0, возвращается -1
begin
 Result := Length(V) - 1;
end;

function ModVector(const V: TVector): real;
// Модуль вектора
var i, N: word;
begin
 Result := 0;
 if(Length(V) = 0)
  // Если длина вектора 0, возвращается 0
  then exit;
 N := High(V);
 for i:=0 to N do
   Result := Result + Sqr(V[i]);
 Result := Sqrt(Result);
end;

function EqVector(const V1, V2: TVector): boolean;
// Эквивалентность двух векторов
var i, N: word;
begin
 Result := false;
 if(Length(V1) = 0) or (Length(V2) = 0)
  then exit;
 N := High(V1);
 if (N <> High(V2))
   then begin
    exit;
   end;
 Result := true;
 for i:=0 to N do
   if( V1[i] <> V2[i])
    then begin
     Result := false;
     exit;
   end;
end;

function SumVector(const V1, V2: TVector): TVector;
// Сложение двух векторов
var i, N: word;
    Vector: TVector;
begin
 ErrorMatrix := false;
 if(Length(V1) = 0) or (Length(V2) = 0)
 then begin
   Messages('Функция SumVector :'#13 +
       'Не задан вектор');
   exit;
 end;
 N := High(V1);
// Размерность векторов должна быть одинаковой
 if(N <> High(V2))
 then begin
   Messages('Функция SumVector :'#13 +
       'Заданы векторы разных размерностей :'#13 +
       IntToStr(N+1) + ' и ' + IntToStr(Length(V2)));
   exit;
 end;
 SetLength(Vector, N+1);
 for i:=0 to N do
     Vector[i] := V1[i] + V2[i];
 Result := Vector;
end;

function SubVector(const V1, V2: TVector): TVector;
// Разность двух векторов
var i, N: word;
    Vector: TVector;
begin
 ErrorMatrix := false;
 if(Length(V1) = 0) or (Length(V2) = 0)
 then begin
   Messages('Функция SubVector :'#13 +
       'Не задан вектор');
   exit;
 end;
 N := High(V1);
// Размерность векторов должна быть одинаковой
 if(N <> High(V2))
 then begin
   Messages('Функция SumVector :'#13 +
       'Заданы векторы разных размерностей :'#13 +
       IntToStr(N+1) + ' и ' + IntToStr(Length(V2)));
   exit;
 end;
 SetLength(Vector, N+1);
 for i:=0 to N do
     Vector[i] := V1[i] - V2[i];
 Result := Vector;
end;

function CMultVector(C: real; const V: TVector): TVector;
// Произведение скаляра (числа) на вектор
var i, N: word;
    Vector: TVector;
begin
 ErrorMatrix := false;
 if(Length(V) = 0)
 then begin
   Messages('Функция CMultVector :'#13 +
       'Не задан вектор');
   exit;
 end;
 N := High(V);
 SetLength(Vector, N+1);
 for i:=0 to N do
   Vector[i] := C * V[i];
 Result := Vector;
end;

function MultSVector(const V1, V2: TVector): real;
// Скалярное умножение двух векторов
var i, N: word;
begin
 ErrorMatrix := false;
 Result := 0;
 if(Length(V1) = 0) or (Length(V2) = 0)
 then begin
   Messages('Функция MultSVector :'#13 +
       'Не задан вектор');
   exit;
 end;
 N := High(V1);
 if (N <> High(V2))
 then begin
   Messages('Функция MultSVector :'#13 +
       'Заданы векторы разных размерностей :'#13 +
       IntToStr(N+1) + ' и ' + IntToStr(Length(V2)));
   exit;
 end;
 for i:=0 to N do
  Result := Result + V1[i] * V2[i];
end;

function CosVector(const V1, V2: TVector): real;
// Косинус угла между двумя векторами
// Реализация далеко не самая быстрая, но компактная
begin
 Result := MultSVector(V1, V2) / (ModVector(V1) * ModVector(V2));
end;

function FVector(const V1, V2: TVector): real;
// Угол в радианах между двумя векторами
begin
 Result := ArcCos(CosVector(V1,V2));
end;

function MultVVector(const V1, V2:  array of real): TVector;
// Векторное умножение двух трехмерных векторов
begin
 ErrorMatrix := false;
 if(Length(V1) <> 3) or (Length(V2) <> 3)
 then begin
   Messages('Функция MultVVector :'#13 +
       'Заданы не трехмерные векторы');
   exit;
 end;
  SetLength(Result, 3);
  Result[0] := V1[1] * V2[2] - V1[2] * V2[1];
  Result[1] := V1[2] * V2[0] - V1[0] * V2[2];
  Result[2] := V1[0] * V2[1] - V1[1] * V2[0];
end;

// **************************************************

function CreateMatrix(M, N: integer): TMatrix;
// Создание матрицы M x N: M - число строк, N - число столбцов
begin
 SetLength(Result, M, N);
end;

function ColMatrix(const A: TMatrix): integer;
// Возвращает число столбцов N
// Если размерность матрицы 0, возвращается -1
begin
 if (Length(A) = 0)
   then Result := -1
   else Result := Length(A[0]);
end;

function RowMatrix(const A: TMatrix): word;
// Возвращает число строк M
begin
 Result := Length(A);
end;

function HighColMatrix(const A: TMatrix): integer;
// Возвращает последний индекс столбца
// Если размерность матрицы 0, возвращается -1
begin
 if (Length(A) = 0)
   then Result := -1
 else Result := High(A[0]);
end;

function HighRowMatrix(const A: TMatrix): integer;
// Возвращает последний индекс строки
// Если размерность матрицы 0, возвращается -1
begin
 if (Length(A) = 0)
   then Result := -1
   else Result := High(A);
end;

function TransMatrix(const A: TMatrix): TMatrix;
// Транспонирование матрицы
var i, j, M, N:word;
    Matr:TMatrix;
begin
 if (Length(A) = 0)
   then begin
    Result := A;
    exit;
   end;
 M := High(A);
 N := High(A[0]);
 SetLength(Matr, N+1, M+1);
 for i:=0 to M do
   for j:=0 to N do
     Matr[j,i] := A[i,j];
 Result := Matr;
end;

function EqMatrix(const A, B: TMatrix): boolean;
// Эквивалентность двух матриц
var i, j, N, M: word;
begin
 Result := false;
 if(Length(A) = 0) or (Length(B) = 0)
  then exit;
 N := High(A);
 M := High(A[0]);
 if (N <> High(B)) or (M <> High(B[0]))
   then exit;
// Цикл сравнения элементов
 Result := true;
 for i:=0 to N do
   for j:=0 to M do
     if( A[i,j] <> B[i,j])
     then begin
       Result := false;
       exit;
     end;
end;

function SumMatrix(const A, B: TMatrix): TMatrix;
// Сумма двух матриц
var i, j, M, N:word;
    Matr:TMatrix;
begin
 ErrorMatrix := false;
 if (Length(A) = 0) or (Length(B) = 0)
 then begin
   Messages('Функция SumMatrix :'#13 +
       'Не задана матрица');
   exit;
 end;
 M := High(A);
 N := High(A[0]);
// Размерность матриц должна быть одинаковой
 if(N <> High(B[0])) or
   (M <> High(B))
 then begin
   Messages('Функция SumMatrix :'#13 +
       'Заданы матрицы разных размерностей :'#13 +
       IntToStr(M+1) + ' x ' + IntToStr(N+1) +
       ' и ' + IntToStr(Length(B)) + ' x ' + IntToStr(Length(B[0])));
   exit;
 end;
 SetLength(Matr, M+1, N+1);
 for i:=0 to M do
   for j:=0 to N do
     Matr[i, j] := A[i, j] + B[i, j];
 Result := Matr;
end;

function SubMatrix(const A, B: TMatrix): TMatrix;
// Разность двух матриц
var i, j, M, N:word;
    Matr:TMatrix;
begin
 ErrorMatrix := false;
 M := High(A);
 if (Length(A) = 0) or (Length(B) = 0)
 then begin
   Messages('Функция SubMatrix :'#13 +
       'Не задана матрица');
   exit;
 end;
 N := High(A[0]);
// Размерность матриц должна быть одинаковой
 if(N <> High(B[0])) or
   (M <> High(B))
 then begin
   Messages('Функция SumMatrix :'#13 +
       'Заданы матрицы разных размерностей :'#13 +
       IntToStr(M+1) + ' x ' + IntToStr(N+1) +
       ' и ' + IntToStr(Length(B)) + ' x ' + IntToStr(Length(B[0])));
   exit;
 end;
 SetLength(Matr, M+1, N+1);
 for i:=0 to M do
   for j:=0 to N do
     Matr[i, j] := A[i, j] - B[i, j];
 Result := Matr;
end;

function CMultMatrix(C: real; const A: TMatrix): TMatrix;
// Произведение скаляра (числа) на матрицу
var i, j, M, N:word;
    Matr:TMatrix;
begin
 ErrorMatrix := false;
 if (Length(A) = 0)
 then begin
   Messages('Функция CMultMatrix :'#13 +
       'Не задана матрица');
   exit;
 end;
 M := High(A);
 N := High(A[0]);
 SetLength(Matr, M+1, N+1);
 for i:=0 to M do
   for j:=0 to N do
     Matr[i, j] := C * A[i, j];
 Result := Matr;
end;

function MultMatrix(const A, B: TMatrix): TMatrix;
// Произведение двух матриц
var i, j, k, M1, N, N2:word;
    Matr:TMatrix;
begin
 ErrorMatrix := false;
 if (Length(A) = 0) or (Length(B) = 0)
 then begin
   Messages('Функция MultMatrix :'#13 +
       'Не задана матрица');
   exit;
 end;
 M1 := High(A);
 N := High(A[0]);
 N2 := High(B[0]);
// Размерности матриц должны быть M x K и K x N
 if(N <> High(B))
 then begin
   Messages('Функция MultMatrix :'#13 +
       'Число столбцов первой матрицы (' + IntToStr(N+1) +
       ')'#13'не совпадает с числом строк второй матрицы (' +
        IntToStr(Length(B)) + ')');
   exit;
 end;
 SetLength(Matr, M1+1, N2+1);
 for i:=0 to M1 do
   for j:=0 to N2 do begin
     Matr[i, j] := 0;
     for k:=0 to N do
       Matr[i, j] := Matr[i, j] + A[i, k] * B[k, j];
   end;
 Result := Matr;
end;

function MultMatrixVector(const A: TMatrix; const V: TVector): TVector;
// Произведение матрицы на вектор
// Размерности матрица должны быть M x N, а вектора N
var i, j, M, N:word;
    Vect:TVector;
begin
 ErrorMatrix := false;
 if (Length(A) = 0)
 then begin
   Messages('Функция MultMatrixVector :'#13 +
       'Не задана матрица');
   exit;
 end;
 M := High(A);
 N := High(A[0]);
// Размерности матрица должны быть M x N, а вектора N
 if(N <> High(V))
 then begin
   Messages('Функция MultMatrixVector :'#13 +
       'Число столбцов  матрицы (' + IntToStr(N+1) +
       ')'#13'не совпадает с размерностью векора (' +
        IntToStr(Length(V)) + ')');
   exit;
 end;
 SetLength(Vect, M+1);
 for i:=0 to M do begin
     Vect[i] := 0;
     for j:=0 to N do
       Vect[i] := Vect[i] + A[i, j] * V[j];
   end;
 Result := Vect;

end;

function DetMatrix(const A: TMatrix): real;
// Вычисление определителя квадратной матрицы
var i, j, N: word;
    Matr: TMatrix;
    R: real;
procedure Iterat(var R:real; const Mat: TMatrix; N: word);
var i, j: word;
    Matr: TMatrix;
    Rtmp, a11:real;
begin
 if(N = 0)
 then begin
   R := R * Mat[0,0];
   exit;
 end;
 SetLength(Matr, N, N);
 a11 := Mat[0,0];
 if(a11 = 0)
 then begin
  // Поиск строки с ненулевым элементом
  for i:=1 to N do
    if (Mat[i,0] <> 0)
     then begin
      for j:= 0 to N do begin
       Rtmp := Mat[0,j];
       Mat[0,j] := Mat[i,j];
       Mat[i,j] := Rtmp;
      end;
      break;
     end;
 a11 := Mat[0,0];
 if(a11 = 0) then begin
   R := 0;
   exit;
 end
 else R := -R;
 end;
 R := R * a11;
 a11 := 1. / a11;
 for i:=0 to N-1 do
   for j:=0 to N-1 do begin
     Matr[i,j] := Mat[i+1,j+1] - Mat[i+1,0] * Mat[0,j+1] * a11;
     if(Abs(Matr[i,j]) < Abs(Mat[i+1,j+1]) * 1E-10)
      then Matr[i,j] := 0;
   end;
 Iterat(R, Matr, N-1);
end;

begin
 ErrorMatrix := false;
 Result := 0;
 if (Length(A) = 0)
 then begin
   Messages('Функция DetMatrix :'#13 +
       'Не задана матрица');
   exit;
 end;
 N:= High(A);
 if(N <> High(A[0]))
 then begin
   Messages('Функция DetMatrix '#13'матрица не квадратная,'#13+
                 'ее размерность ' + IntToStr(N+1) + ' x ' +
                 IntToStr(High(A[0])+1));
   exit;
 end;
 SetLength(Matr, N+1, N+1);
 for i:= 0 to N do
   for j:=0 to N do
     Matr[i,j] := A[i,j];
 R := 1;
 Iterat(R, Matr, N);
 Result := R;
end;

function ReversMatrix(const A: TMatrix): TMatrix;
// Прямое вычисление обратной матрицы
// Матрица A должна быть квадратной
var i, j, i1, j1, k1, k2, N: word;
    Matr, Dop: TMatrix;
    Rtmp, Det: real;
begin
 ErrorMatrix := false;
 if (Length(A) = 0)
 then begin
   Messages('Функция ReversMatrix :'#13 +
       'Не задана матрица');
   exit;
 end;
 N:= High(A);
 if(N <> High(A[0]))
 then begin
   Messages('Функция DetMatrix :'#13'матрица не квадратная,'#13+
                 'ее размерность ' + IntToStr(N+1) + ' x ' +
                 IntToStr(High(A[0])+1));
   exit;
 end;
// Вариант 1 экономнее
// (не вычисляется определитель всей матрицы),
// но при разложении возникает > ошибок округления
// и часто вместо Det = 0 получается малое число
 { Det := 0; // Вариант 1}
 Det := DetMatrix(A); // Вариант 2
 if(Det = 0) // Вариант 2
 then begin
   Messages('Функция DetMatrix '#13'матрица особенная');
   exit;
 end;
 SetLength(Matr, N+1, N+1);
 SetLength(Dop, N, N);
 for i1:= 0 to N do
    for j1:=0 to N do begin
     k1 := 0;
     for i:= 0 to N do begin
      if(i = i1) then Continue;
      k2 := 0;
      for j:=0 to N do begin
       if (j = j1) then Continue;
       Dop[k1,k2] := A[i,j];
       Inc(k2);
      end;
      Inc(k1);
     end;
     Rtmp := DetMatrix(Dop);
     if (Odd(i1 + j1))
       then Rtmp := -Rtmp;
{     Matr[j1,i1] := Rtmp; // Вариант 1}
     Matr[j1,i1] := Rtmp / Det; // Вариант 2
{     if(i1 = 0) then Det := Det + A[0,j1] * Rtmp; // Вариант 1}
    end;
{ if(Det = 0) // Вариант 1
 then begin
   Messages('Функция DetMatrix '#13'матрица особенная');
   exit;
 end;
 for i:= 0 to N do
   for j:=0 to N do
     Matr[i,j] := Matr[i,j] / Det;}
 Result := Matr;
end;

function LSystem(const A: TMatrix; const VRight: TVector): TVector;
// Решение системы линейных уравнений
// Матрица A должна быть квадратной: N x N
// Вектор VRight - той же размерности, что и A
var Matr: TMatrix;
begin
 ErrorMatrix := false;
 if (Length(A) = 0)
 then begin
   Messages('Функция LSystem :'#13 +
       'Не задана матрица');
   exit;
 end;
 Matr := ReversMatrix(A);
 Result := MultMatrixVector(Matr,VRight);
end;

function ReversMatrix2(const A: TMatrix): TMatrix;
// Прямое вычисление обратной матрицы
// с вычислением определителя через алгебраические дополнения
// Матрица A должна быть квадратной
var i, j, i1, j1, k1, k2, N: word;
    Matr, Dop: TMatrix;
    Rtmp, Det: real;
begin
 ErrorMatrix := false;
 if (Length(A) = 0)
 then begin
   Messages('Функция ReversMatrix2 :'#13 +
            'Не задана матрица');
   exit;
 end;
 N:= High(A);
 if(N <> High(A[0]))
 then begin
   Messages('Функция ReversMatrix2 :' +
            #13'матрица не квадратная,'#13 +
            'ее размерность ' + IntToStr(N+1) + ' x ' +
            IntToStr(High(A[0])+1));
   exit;
 end;
 Det := 0;
 SetLength(Matr, N+1, N+1);
 SetLength(Dop, N, N);
 for i1:= 0 to N do
    for j1:=0 to N do begin
     k1 := 0;
     for i:= 0 to N do begin
      if(i = i1) then Continue;
      k2 := 0;
      for j:=0 to N do begin
       if (j = j1) then Continue;
       Dop[k1,k2] := A[i,j];
       Inc(k2);
      end;
      Inc(k1);
     end;
     Rtmp := DetMatrix(Dop);
     if (Odd(i1 + j1))
       then Rtmp := -Rtmp;
     Matr[j1,i1] := Rtmp;
     if(i1 = 0) then begin
      Det := Det + A[0,j1] * Rtmp;
     end;
    end;
 if(Det = 0)
 then begin
   Messages('Функция ReversMatrix2 '#13'матрица особенная');
   exit;
 end;
 for i:= 0 to N do
   for j:=0 to N do
     Matr[i,j] := Matr[i,j] / Det;
 Result := Matr;
end;

procedure LU(var Matr: TMatrix; var V: TInt);
// LU разложение матрицы
// Matr - заданная матрица, заменяемая в результате разложенной
// V - массив индексов, определяющих последовательность 
// уравнений с учетом перестановок в процессе LU-разложения
var i, j, k, N: word;
    Rtmp, Rtmp2: real;
begin
 ErrorMatrix := false;
 if (Length(Matr) = 0)
 then begin
   Messages('Функция LU :'#13 + 'Не задана матрица');
   exit;
 end;
 N:= High(Matr);
 if(N <> High(Matr[0]))
 then begin
   Messages('Функция LU :'#13'матрица не квадратная,'#13+
            'ее размерность ' + IntToStr(N+1) + ' x ' +
            IntToStr(High(Matr[0])+1));
   exit;
 end;
// Заполнение V начальными значениями
 SetLength(V, N+1);
 for i:=0 to N do
   V[i] := i;

 for i:=0 to N do begin
  if(Matr[i, i] = 0)
   then begin
    // Поиск первого ненулевого элемента в строке i
    for j:=i+1 to N do
     if(Matr[i, j] <> 0)
      then begin
    // Перестановка столбцов i и j
       for k:=0 to N do begin
        Rtmp := Matr[k, i];
        Matr[k, i] := Matr[k, j];
        Matr[k, j] := Rtmp;
       end;
       // Перестановка последовательности индексов столбцов
       k := V[i];
       V[i] := V[j];
       V[j] := k;
       break;
      end;
    if(Matr[i, i] = 0)
    // Ненулевой элемент не найден
     then begin
      Messages('Функция LU :'#13'матрица особенная');
      exit;
     end;
   end;
  Matr[i, i] :=  1. / Matr[i, i];

  // Формирование столбца L-матрицы
  for j:=i+1 to N do begin
   Rtmp := 0;
   if(i > 0) then
    for k:=0 to i-1 do begin
     Rtmp2 := Rtmp;
     Rtmp := Rtmp + Matr[j, k] * Matr[k,i];
     if (Abs(Rtmp) < Abs(Rtmp2) * 1E-10) then Rtmp := 0;
    end;
   Matr[j, i] := (Matr[j, i] - Rtmp) * Matr[i, i];
   if (Abs(Matr[j, i]) < Abs(Rtmp) * 1E-10)
    then Matr[j, i] := 0;
  end;

  // Формирование (i+1)-ой строки U-матрицы
  for j:=i+1 to N do
   for k:=0 to i do begin
    Rtmp := Abs(Matr[i+1, j]);
    Matr[i+1, j] := Matr[i+1, j] - Matr[i+1,k] * Matr[k,j];
    if(Abs(Matr[i+1, j]) < Rtmp*1e-10) then Matr[i+1, j] := 0;
   end;
 end;
end;

procedure LUMax(var Matr: TMatrix; var V: TInt);
// LU разложение матрицы с выбором главного элемента
// Matr - заданная матрица, заменяемая в результате разложенной
// V - массив индексов, определяющих последовательность 
// уравнений с учетом перестановок в процессе LU-разложения
var i, j, k, im, N: word;
    Rtmp, Rtmp2: real;
begin
 ErrorMatrix := false;
 if (Length(Matr) = 0)
 then begin
   Messages('Функция LU :'#13 + 'Не задана матрица');
   exit;
 end;
 N:= High(Matr);
 if(N <> High(Matr[0]))
 then begin
   Messages('Функция LU :'#13'матрица не квадратная,'#13 +
            'ее размерность ' + IntToStr(N+1) + ' x ' +
            IntToStr(High(Matr[0])+1));
   exit;
 end;
// Заполнение V начальными значениями
 SetLength(V, N+1);
 for i:=0 to N do
   V[i] := i;

 for i:=0 to N do begin
  Rtmp := Abs(Matr[i, i]);
  im := i;             // Индекс максимального элемента
  // Поиск индекса максимального элемента в строке i
  for j:=i+1 to N do
   if(Abs(Matr[i, j]) > Rtmp)
    then begin
     Rtmp := Abs(Matr[i, j]);
     im := j;
    end;
  if (im > i)
   then begin
    // Перестановка столбцов i и im
       for k:=0 to N do begin
        Rtmp := Matr[k, i];
        Matr[k, i] := Matr[k, im];
        Matr[k, im] := Rtmp;
       end;
       // Перестановка последовательности индексов столбцов
       k := V[i];
       V[i] := V[im];
       V[im] := k;
   end;
  if(Matr[i, i] = 0)
    // Ненулевой элемент не найден
   then begin
    Messages('Функция LU :'#13'матрица особенная');
    exit;
   end;
  Matr[i, i] :=  1. / Matr[i, i];

  // Формирование столбца L-матрицы
  for j:=i+1 to N do begin
   Rtmp := 0;
   if(i > 0) then
    for k:=0 to i-1 do begin
     Rtmp2 := Rtmp;
     Rtmp := Rtmp + Matr[j, k] * Matr[k,i];
     if (Abs(Rtmp) < Abs(Rtmp2) * 1E-10) then Rtmp := 0;
    end;
   Matr[j, i] := (Matr[j, i] - Rtmp) * Matr[i, i];
   if (Abs(Matr[j, i]) < Abs(Rtmp) * 1E-10)
    then Matr[j, i] := 0;
  end;

  // Формирование (i+1)-ой строки U-матрицы
  for j:=i+1 to N do
   for k:=0 to i do begin
    Rtmp := Abs(Matr[i+1, j]);
    Matr[i+1, j] := Matr[i+1, j] - Matr[i+1,k] * Matr[k,j];
    if(Abs(Matr[i+1, j]) < Rtmp*1e-10) then Matr[i+1, j] := 0;
   end;
 end;
end;

function LUVector(const B: TVector;
                  const Matr: TMatrix; const V: TInt): TVector;
// Решение системы линейных уравнений
// с уже проведенным LU-разложением (реультаты в Matr и V)
// и с заданным вектором правых частей B
var i, j, N: word;
    Rtmp: real;
    Vz: TVector;
begin
 N:= High(Matr);
 if(N <> High(Matr[0]))
 then begin
   Messages('Функция LUVector :'#13'матрица не квадратная,'#13+
            'ее размерность ' + IntToStr(N+1) + ' x ' +
            IntToStr(High(Matr[0])+1));
   exit;
 end;
 if(N <> High(V))
 then begin
  Messages('Функция LUVector :' +
         #13'размер вектора не соответствует размеру матрицы');
  exit;
 end;
 SetLength(Vz, N+1);
 SetLength(Result, N+1);
 // Прямой ход
 Vz[0] := B[0];
 for i:=1 to N do begin
   Rtmp := 0;
   for j:=0 to i-1 do
    Rtmp := Rtmp + Matr[i, j] * Vz[j];
   Vz[i] := B[i] - Rtmp;
 end;
 // Обратный ход
 Vz[N] := Vz[N] * Matr[N, N];
 for i:=N-1 downto 0 do begin
   Rtmp := 0;
   for j:=i+1 to N do
    Rtmp := Rtmp + Matr[i, j] * Vz[j];
   Vz[i] := (Vz[i] - Rtmp) * Matr[i, i];
 end;
 // Восстановление начальной последовательности переменных
 for i:=0 to N do
  Result[V[i]] := Vz[i];
end;

function LUSystem(const A: TMatrix;
                       const VRight: TVector): TVector;
// Решение системы линейных уравнений
// с помощью LU-разложения без выбора главного элемента
// с сохранением исходной матрицы
var Matr: TMatrix;
    V: TInt;      // массив индексов с учетом перестановок
    Vz: TVector;
    i, j, N: word;
    Rtmp: real;
begin
 N:= High(A);
 if(N <> High(A[0]))
 then begin
   Messages('Функция LUSystem :' +
            #13'матрица не квадратная,'#13+
            'ее размерность ' + IntToStr(N+1) + ' x ' +
            IntToStr(High(A[0])+1));
   exit;
 end;
 SetLength(Matr, N+1, N+1);
 for i:=0 to N do
   for j:=0 to N do
     Matr[i, j] := A[i, j];
 SetLength(Vz, N+1);
 SetLength(Result, N+1);
 LU(Matr, V);
 if(ErrorMatrix) then exit;
 // Прямой ход
 Vz[0] := VRight[0];
 for i:=1 to N do begin
   Rtmp := 0;
   for j:=0 to i-1 do
    Rtmp := Rtmp + Matr[i, j] * Vz[j];
   Vz[i] := VRight[i] - Rtmp;
 end;
 // Обратный ход
 Vz[N] := Vz[N] * Matr[N, N];
 for i:=N-1 downto 0 do begin
   Rtmp := 0;
   for j:=i+1 to N do
    Rtmp := Rtmp + Matr[i, j] * Vz[j];
   Vz[i] := (Vz[i] - Rtmp) * Matr[i, i];
 end;
 // Восстановление начальной последовательности переменных
 for i:=0 to N do
  Result[V[i]] := Vz[i];
end;

function LUMaxSystem(const A: TMatrix; const VRight: TVector):
                                                       TVector;
// Решение системы линейных уравнений
// с помощью LU-разложения с выбором главного элемента
// с сохранением исходной матрицы
var Matr: TMatrix;
    V: TInt;      // массив индексов с учетом перестановок
    Vz: TVector;
    i, j, N: word;
    Rtmp: real;
begin
 N:= High(A);
 if(N <> High(A[0]))
 then begin
   Messages('Функция LUSystem :' +
            #13'матрица не квадратная,'#13 +
            'ее размерность ' + IntToStr(N+1) + ' x ' +
            IntToStr(High(A[0])+1));
   exit;
 end;
 SetLength(Matr, N+1, N+1);
 for i:=0 to N do
   for j:=0 to N do
     Matr[i, j] := A[i, j];
 SetLength(Vz, N+1);
 SetLength(Result, N+1);
 LUMax(Matr, V);
 if(ErrorMatrix) then exit;
 // Прямой ход
 Vz[0] := VRight[0];
 for i:=1 to N do begin
   Rtmp := 0;
   for j:=0 to i-1 do
    Rtmp := Rtmp + Matr[i, j] * Vz[j];
   Vz[i] := VRight[i] - Rtmp;
 end;
 // Обратный ход
 Vz[N] := Vz[N] * Matr[N, N];
 for i:=N-1 downto 0 do begin
   Rtmp := 0;
   for j:=i+1 to N do
    Rtmp := Rtmp + Matr[i, j] * Vz[j];
   Vz[i] := (Vz[i] - Rtmp) * Matr[i, i];
 end;
 // Восстановление начальной последовательности переменных
 for i:=0 to N do
  Result[V[i]] := Vz[i];
end;

end.
