SDS-разложение матрицы (метод Холецкого или квадратных корней)

Type
  TMatrix=array[1..10,1..10] of real;
  BVector=array[1..10] of real;
Var
  Mat:TMatrix; b:BVector; f,n,i,k,j:integer;
procedure SDS;
var i,j,k:integer;s1:real;s,d:TMatrix;f:boolean;
begin
f:=true;
for i:=1  to n do
  for j:=1 to n do
    if Mat[i,j]<>Mat[j,i] then f:=false;
if f=true then
begin
//d[1,1]:=sgn(Mat[1,1]);
s[1,1]:=sqrt(Mat[1,1]);
for i:=1 to n do
  begin
    s1:=0;
    for k:=1 to (i-1) do
    s1:=s1+s[k,i]*s[k,i];
    s[i,i]:=sqrt(Mat[i,i]-s1);
    for j:=(i+1) to n do
      begin
      s1:=0;
      for k:=1 to (i-1) do
      s1:=s1+s[k,i]*s[k,j];
      //d[i,i]:=sign(Mat[i,i]-s1*s1*d[k,k]);
      s[i,j]:=(Mat[i,j]-s1)/s[i,i];
      end;
  end;
writeln('Matrica S*:');
for i:=1 to n do
  begin
    for j:=1 to n do
    write(s[j,i]:2:0,' ');writeln;
  end;
writeln('Matrica S:');
for i:=1 to n do
  begin
    for j:=1 to n do
    write(s[i,j]:2:0,' ');writeln;
  end;
end
  else writeln('невозможно разложить S*DS-методом, т.к. матрица не эрмитовая!');
end;
begin
Write(' n=');
    Readln(n);
    Writeln('M(',n,'x',n,'):');
    For i:=1 to n do
      For j:=1 to n do
        Begin
          Write('m[',i,',',j,']=');
          Readln(Mat[i,j]);
        End;
    Writeln;
    Writeln(' b(',n,'x','1):');
    For i:=1 to n do
      Begin
        Write('b[',i,',','1]=');
        Readln(b[i]);
      End;
      SDS;
end.