LU-разложение матрицы с частичным выбором ведущего элемента - Pascal - Ответ 10773928

{ФИО: Харламенко Артём Валерьевич
Вариант:9}
Program lr6_2;
const n = 4;
type rmatr = array[1..n, 1..n] of real;
     ivect = array[1..n] of integer;
var i,j:integer;
    A:rmatr;
    ip:ivect;
    eps, det: real;
    ifsolve: boolean;
 
{
Назначение:
Выполняет LU разложение матрицы с частичным выбором
ведущего элемента
ФОРМАЛЬНЫЕ ПАРАМЕТРЫ
 n – размер матрицы А(входной) 
 A – при входе – элементы матрицы А, 
     при выходе – искомые элементы матриц L и U 
 eps – при входе - относительный, 
       при выходе- абсолютный (eps=eps*||A||)
       уровень оценки ведущего элемента 
       как равного нулю 
 ip – массив перестановок строк матрицы А (выходной)
 det – определитель матрицы А (выходной)
 ifsolve – при выходе из подпрограммы равен:  true
           если разложение выполнено успешно,
           false – если хотя бы один ведущий 
           элемент по модулю меньше eps и 
           разложение невозможно
}
Procedure DECOMP_2(n:integer; var A:rmatr; var eps: real;
                    var ip:ivect;
                    var det:real;
                    var ifsolve:boolean);
  label 10;
  var i,j,k,km,ipm: integer;
      s,norm,am: real;
  begin
    ifsolve:=true;
    det:=1.0;
    norm:=0.0;
  
    for i:=1 to n do begin 
      ip[i]:= i;
      s:=0.0;
      for j:=1 to n do s:=s+abs(A[i,j]);
      if norm<s then norm:=s;
    end;
    
    eps := eps*norm;
    
    for k:=1 to n-1 do begin
      am:=abs(A[k,k]);
      km:=k;
      
      for i:=k+1 to n do begin
        if abs(A[i,k])>am then begin
          am:=abs(A[i,k]);
          km:=i;
        end;
        if km<>k then begin
          det := (-det);
          ipm:=ip[km];
          ip[km]:=ip[k];
          ip[k]:=ipm;
          for j:=1 to n do begin
            am := A[km, j];
            A[km,j]:=A[k,j];
            A[k,j]:=am;
          end;
        end;
      end;
    
     if abs(A[k,k])<eps then begin 
       ifsolve:=false; 
       goto 10; 
     end;
     
     for i:=1 to n do begin 
      am:=A[i,k]/A[k,k];
      A[i,k]:=am;
      for j:=1 to n do A[i,j]:=A[i,j]-am*A[k,j];
      end;
     end;
     
  10: if ifsolve and (abs(A[n,n])>= eps) then
        for k:=1 to n do det:=det*A[k,k] else det:=0.0;
  
  end;
  
begin // основной блок программы
for i:=1 to n do
  for j:=1 to n do begin
  //write('A[',i,',',j,']= ');
  //readln(A[i,j]);
  A[i,j]:=4;
  end;
det:=1.0;
//write('Введите точность вычисления eps= ');
//readln(eps);
eps:=0.000001;
DECOMP_2(n, A,eps,ip,det,ifsolve);
writeln('Определитель= ',det); // выводим определитель
writeln('U',a); // печатаем ступенчатый вид матрицы
writeln('Массив перестановок', ip); 
end.