Obtención de Base ortogonal mediante Gram-Schmidt

El siguiente programa está dedicado al cálculo de una base ortogonal a partir de un conjunto de vectores linealmente independiente, tomándose como producto escalar el producto escalar usual, mediante el método de Gram. Este método se puede deducir de manera sencilla, y es una de las maneras más simples de calcular una base ortogonal.
El mecanismo del programa es el siguiente:

1.   Se introducen las componentes de los vectores en una matriz, que se triangulariza superior. En el caso de que el último pivote sea 0, quiere decir que hay por lo menos un vector que es combinación lineal de los demás, con lo cual no forman una base, si acaso un sistema generador de un espacio vectorial.
2. Se aplica el método de Gram (que en realidad se atribuye a Laplace) escogiendo los vectores en el orden en que han sido introducidos (dependiendo del orden de los mismos la base será distinta, pero todas generarán el mismo espacio).
3. Se dividen los vectores de la base ortognal por su norma, de manera que sean vectores unitarios que forman una base ortonormal.

{programa que apartir de unos vectores que forman una base de un espacio vectorial
(deben ser linealmente independientes) devuelve una base ortonormal respecto al
producto escalar usual. Basado en el metodo de Gram-Schmidt}

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

type
  tpDimension=1..tanMax;
  tpContador=1..MAXINT;
  tpVector=array[1..tanMax] of real;
  tpMatriz=array[1..tanMax] of tpVector;

var
  mat,matorto:tpMatriz;
  vectores,componentes:tpDimension;


procedure EscribirMatriz(m:tpMatriz; filas,columnas:tpDimension);
  var
    i,j:tpContador;

  begin
    for i:=1 to filas do begin
      write('Vector:',i,':    (');
      for j:=1 to columnas do
        write(m[i,j]:7:2,' ');
      writeln(')');
    end
  end;

procedure CargaVectores(VAR v,c:tpDimension; VAR m:tpMatriz);
  var
    i,j:tpContador;
  begin
    write('Introduzca el numero de vectores: ');readln(v);
    write('Introduzca el numero de componentes: ');readln(c);
    for i:=1 to v do begin
      writeln('Vector numero ',i,':');
        for j:=1 to c do begin
          write('Componenetes numero ',j,': ');readln(m[i,j]) end;
      writeln
    end
  end;

function lineal(matriz:tpMatriz; columnas,filas:tpDimension):boolean;
  Const
    precision=0.000001;

  Var
    min,pasoF,pasoC:tpDimension;
    cont:tpContador;

  Procedure SeleccionarFila(f,c,pF,pC:tpDimension; VAR m:tpMatriz);
    Var
      i,j:tpContador;
      aux:real;

    Begin
      i:=pF;
      while (m[i,pC]<precision) and (i<f) and (m[i,pC]>-precision) do
        i:=i+1;
      If i<>f then
        begin
          for j:=pC to c do begin
            aux:=m[pF,j];
            m[pF,j]:=m[i,j];
            m[i,j]:=aux
          end
        end
    end;


  Procedure restarFilas(f,c:tpDimension; VAR pF:tpDimension; VAR pC:tpDimension; VAR m:tpMatriz);
    var
      i,j:tpContador;
      aux:real;

    begin
      if m[pF,pC]<>0 then begin
        for i:=(pF+1) to f do begin
          aux:=m[i,Pc];
          for j:=pC to c do
            m[i,j]:=m[i,j]-aux*(m[pF,j]/m[pF,pC]);
        end;
        pF:=pF+1
      end;
      pC:=pC+1
    end;


  procedure MinDim(f,c:tpDimension; VAR mini:tpDimension);
    begin
      if c<f then mini:=c else mini:=f
    end;

  begin
    MinDim(filas,columnas,min);
    pasoF:=1; pasoC:=1;
    for cont:=1 to min do
      begin
        SeleccionarFila(filas,columnas,pasoF,pasoC,matriz);
        restarFilas(filas,columnas,pasoF,pasoC,matriz)
      end;
    if (matriz[filas,columnas]<precision) and (matriz[filas,columnas]>-precision) then
      lineal:=false
    else
      lineal:=true
  end;
{final de la function lineal}

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

procedure ortogonaliza(VAR morto:tpMatriz; m:tpMatriz; comp:tpContador);
  var
    i,j,n:tpContador;
    a:real;

  begin
  for n:=1 to vectores do begin
    for j:=1 to comp do morto[n,j]:=m[n,j];
    for i:=1 to (n-1) do begin
      a:=aux(morto,m,n,i,comp);
      for j:=1 to comp do
        morto[n,j]:=morto[n,j]-a*morto[i,j]
      end
    end;
  end;

procedure normaliza(VAR morto:tpMatriz);
  var
    i,j:tpContador;
    a:real;
  begin
    for i:=1 to vectores do begin
      a:=0;
      for j:=1 to componentes do
        a:=a+morto[i,j]*morto[i,j];
      a:=sqrt(a);
      for j:=1 to componentes do
        morto[i,j]:=morto[i,j]/a
    end
  end;

{comienzo programa}
begin
  writeln('OBTIENE BASE ORTOGONAL Y OROTONORMAL MEDIANTE GRAM-SCHMIDT');
  writeln;
  CargaVectores(vectores,componentes,mat);
  writeln;
  writeln('Vectores Introducidos:');
  EscribirMatriz(mat,vectores,componentes);writeln;writeln;
  If lineal(mat,componentes,vectores) then
    begin
      ortogonaliza(matorto,mat,componentes);
      writeln('Base ortogonal del Espacio engendrado por los vectores introducidos: ');
      EscribirMatriz(matorto,vectores,componentes);writeln;writeln;
      normaliza(matorto);
      writeln('Base ortonormal del Espacio engendrado por los vectores introducidos: ');
      EscribirMatriz(matorto,vectores,componentes);writeln;writeln;
    end
  else
    writeln('El conjunto de vectores no es base puesto que no es linealmente independiente');
  writeln;
  readln
end.