разложить матрицу методом qr отражений

Procedure QR(A:TMatrix);
Var
  R,Q,H,T,E:TMatrix;
  i,j,m,k:integer;
  s,b,p,y:real;
Begin
  s:=0;
  For i:=1 to n do
    For j:=1 to n do
      Begin
        T[i,j]:=0;
      End;
  E:=T;
  For i:=1 to n do E[i,i]:=1;
  R:=T; Q:=E; H:=E;  {1.0}
  For i:=1 to n-1 do  {2.0}
    Begin
      For k:=i to n do s:=s+sqr(A[k,i]);  {2.1}
      If s=0 then continue;  {2.2}
      If A[i,i]<0 then b:=sqrt(s) else b:=(-1)*sqrt(s); {2.3}
      y:=1/(b*(b-A[i,i]));  {2.4}
      A[i,i]:=A[i,i]-b;  {2.5}
      H:=E;  {2.6}
      For m:=i to n do  {2.7}
        For j:=i to n do
          Begin
            If (j<>m) then H[m,j]:=(-1)*y*A[m,i]*A[j,i]
              Else H[m,j]:=1-y*sqr(A[m,i]);
            H[j,m]:=H[m,j];
          End;
      For m:=i to n do  {2.8}
        For j:=i+1 to n do
          Begin
            p:=0;
            For k:=i to n do p:=p+A[k,i]*A[k,j];
            T[m,j]:=A[m,j]-y*A[m,i]*p;
          End;
      For m:=i to n do  {2.9}
        For j:=i+1 to n do A[m,j]:=T[m,j];
      A[i,i]:=b;  {2.10}
      For m:=i+1 to n do A[m,i]:=0;  {2.11}
      For m:=1 to n do  {2.12}
        For j:=1 to n do
          Begin
            p:=0;
            for k:=i to n do p:=p+Q[m,k]*H[k,j];
            T[m,j]:=p;
          End;
      For m:=1 to n do  {2.13}
        For j:=1 to n do Q[m,j]:=T[m,j];
    End;
  For m:=1 to n do  {3.0}
    For j:=1 to n do R[m,j]:=A[m,j];
 
  Writeln('Вывод матрицы R (правая треугольная):'); Writeln;
  For i:=1 to n do
    Begin
      For j:=1 to n do Write(R[i,j]:7:2);
      Writeln; Writeln;
    End;
  Writeln('Вывод матрицы Q (ортогональная):'); Writeln;
  For i:=1 to n do
    Begin
      For j:=1 to n do Write(Q[i,j]:7:2);
      Writeln; Writeln;
    End;
  Writeln('Нажмите Enter, чтобы выполнить проверку...');
  Readln;  Clrscr;
  Writeln('Проверка A=QR:');
  Writeln('Вывод матрицы A:'); Writeln;
  For k:=1 to n do
    Begin
      For j:=1 to n do
        Begin
          A[k,j]:=0;
          For i:=1 to n do
            A[k,j]:=A[k,j]+Q[k,i]*R[i,j];
          Write(A[k,j]:7:2);
        End;
      Writeln; Writeln;
    End;
  Writeln('Исходная матрица M:'); Writeln;
  For i:=1 to n do
    Begin
      For j:=1 to n do Write(Mat[i,j]:7:2);
      Writeln; Writeln;
    End;
End;