QR descomposition (Mediante Gram)

In linear algebra, a QR decomposition (also called a QR factorization) of a matrix is a decomposition of a matrix A into a product A=QR of an orthogonal matrix Q and an upper triangular matrix R. QR decomposition is often used to solve the linear least squares problem,

We can use QR decomposition to find the absolute value of the determinant of a square matrix;

det(A)=det(Q)*det(R)=1*det(R) 

and determiant of a upper triangular matrix is very easy:

program QR(input,output);
const
  tanMax=20;

type
  tpDimension=1..tanMax;
  tpContador=1..MAXINT;
  tpVector=array[1..tanMax] of real;
  tpMatriz=record
    m:array[1..tanMax] of tpVector;
    f,c:tpDimension
  end;

var
  A,Q,R,B:tpMatriz;{x,y,z,b}

procedure CargaDatos(VAR x:tpMatriz);
  var
    i,j:tpContador;
  begin
    writeln('QR DESCOMPOSITION:'); writeln;
    write('Number of rows: ');readln(x.f);
    write('Number of Columns: ');readln(x.c);
    writeln;
    for i:=1 to x.f do begin
        for j:=1 to x.c do begin
          write('Position (',i,',',j,'): ');readln(x.m[i,j]) end;
      writeln
    end
  end;

procedure EscribirMatriz(x:tpMatriz);
var
  i,j:tpContador;

begin
  writeln;
  for i:=1 to x.f do begin
    for j:=1 to x.c do
      write(x.m[i,j]:7:3,' ');
    writeln
  end;
  writeln; writeln
end;

function aux(y,x:tpMatriz; n,i:tpContador):real;
  var
    j:tpContador;
    aux1,aux2:real;
  begin
    aux1:=0;
    aux2:=0;
    for j:=1 to x.f do begin
      aux1:=aux1+x.m[j,n]*y.m[j,i];
      aux2:=aux2+y.m[j,i]*y.m[j,i]
    end;
    aux:=aux1/aux2
  end;

procedure ortogonaliza(VAR y:tpMatriz; x:tpMatriz);
  var
    i,j,n:tpContador;
    a:real;

  begin
  y.f:=x.f; y.c:=x.c;
  for n:=1 to x.f do begin
    for j:=1 to x.f do y.m[j,n]:=x.m[j,n];
    for i:=1 to (n-1) do begin
      a:=aux(y,x,n,i);
      for j:=1 to x.f do
        y.m[j,n]:=y.m[j,n]-a*y.m[j,i]
      end
    end;
  end;

procedure normaliza(VAR y:tpMatriz);
  var
    i,j:tpContador;
    a:real;
  begin
    for i:=1 to y.c do begin
      a:=0;
      for j:=1 to y.f do
        a:=a+y.m[j,i]*y.m[j,i];
      a:=sqrt(a);
      for j:=1 to y.f do
        y.m[j,i]:=y.m[j,i]/a
    end
  end;

procedure transponer(VAR b:tpMatriz; x:tpMatriz);
  var
    i,j:tpContador;
  begin
    b.f:=x.c; b.c:=x.f;
    for i:=1 to b.f do
      for j:=i to b.c do begin
        b.m[i,j]:=x.m[j,i];
        b.m[j,i]:=x.m[i,j]
      end
  end;

procedure multiplica(VAR x:tpMatriz; y,z:tpMatriz);
  {multiplica Y*Z en este orden y lo guarda en X (Y.c=Z.f)}
  var
    i,j,k:tpDimension;
  begin
    x.f:=y.f; x.c:=z.c;
    for i:=1 to y.f do
      for j:=1 to z.c do begin
        x.m[i,j]:=0;
        for k:=1 to y.c do x.m[i,j]:=x.m[i,j]+y.m[i,k]*z.m[k,j]
      end
  end;


{comienzo del programa}
begin
  CargaDatos(A);
  EscribirMatriz(A);
  ortogonaliza(Q,A);
  normaliza(Q);
  transponer(B,Q);
  multiplica(R,B,A);
  writeln('Matrix Q:'); EscribirMatriz(Q);
  writeln('Matrix R:'); EscribirMatriz(R);
  multiplica(B,Q,R);
  writeln('QR is:'); EscribirMatriz(B);
  writeln('It should to be like A');
  readln
end.

Salud.