program hqr_test(input,output,ein,aus);

use mv_ari;

var m,n : integer;
    ein, aus : text;


procedure H_QR(A : rmatrix; var Q,R : rmatrix; var err : boolean);

var u,v,pQ : rvector[1..ub(A,1)];
    pA : rvector[1..ub(A,2)];
    dtd,d,msd,grN : real;
    m,n,k,l,j,i : integer;

begin

  m := ub(A,1); n := ub(A,2);

  {Initialisierungen}
  err := false;
  Q := id(m);  {Q = (m x m)-Einheitsmatrix}
  k := 1;

  while k <= n do
    begin

      dtd := #*(for i := k to m sum(A[i,k]*A[i,k])); {ev. Unter-o.Ueberlaufgefahr!}

      if dtd = 0.0 then begin err := true; k := n+1; end
      else {dtd > 0.0}
        begin

          d := sqrt(dtd);

          {Ausloeschung verhindernde Festlegung von sigma}
          if A[k,k] > 0.0 then msd := d
          else msd := -d;

          {im k-ten Schritt festzulegende Elemente von R}
          for i := 1 to k-1 do
            R[i,k] := A[i,k];
          R[k,k] := -msd;

          {Komponenten u[k],...,u[m] von u bez. des aktuellen Orthonormalsystems}
          u[k] := A[k,k] + msd;
          for i := k+1 to m do
          u[i] := A[i,k];

          {Berechnung von grN und v[i] := u[i]/grN, i=k,...,m}
          grN := #*(dtd + msd*A[k,k]); { = dtd + d*abs(A[k,k]) }
          for i := k to m do
            v[i] := u[i]/grN;

          {Berechnung der neuen Koeffizienten von A}
          {1. Berechnung von pA[k+1],...,pA[n]}
          for l := k+1 to n do
            pA[l] := #*(for j := k to m sum(A[j,l]*u[j]));
          {2. Berechnung von A[i,l], l=k+1,...,n,i=k,...,m}
          for l := k+1 to n do
            for i := k to m do
              A[i,l] := #*(A[i,l] - v[i]*pA[l]);

          {Berechnung der neuen Koeffizienten von Q}
          {1. Berechnung von pQ[1],...,pQ[m]}
          for l := 1 to m do
            pQ[l] := #*(for j := k to m sum(Q[l,j]*u[j]));
          {2. Berechnung von Q[l,i], l = 1,...,m,i=k,...,m}
          for l := 1 to m do
            for i := k to m do
              Q[l,i] := #*(Q[l,i] - v[i]*pQ[l]);

          k := k+1;

        end;{else}
    end;{while}
end;{HQ_R}


procedure main(m,n : integer);

var i,j : integer;
    A : rmatrix[1..m,1..n];
    Q : rmatrix[1..m,1..m];
    R : rmatrix[1..n,1..n];
    err : boolean;

begin

  {Zeilenweises Einlesen von A aus der Datei ein}
  for i := 1 to m do
    for j := 1 to n do
      read(ein,A[i,j]);

  H_QR(A,Q,R,err);

  {Ausgabe in die Datei aus}
  rewrite(aus);
  if err = true then {irregulaerer Fall}
    begin
      writeln(aus);
      writeln(aus,'Spalten von A (num.) l.a. !');
      writeln(aus,'Es wird keine QR-Zerlegung berechnet!');
    end
  else {regulaerer Fall}
    begin
      writeln(aus);
      writeln(aus,'Berechnete QR-Zerlegung der eingelesenen Matrix A');
      writeln(aus);
      for j := 1 to n do
        begin
          for i := 1 to m do
            writeln(aus,'Q[',i,',',j,'] = ',Q[i,j]);
          writeln(aus);
        end;
      writeln(aus);
      for j := 1 to n do
        begin
          for i := 1 to j do
            writeln(aus,'R[',i,',',j,'] = ',R[i,j]);
          writeln(aus);
        end;
    end;{else}

end;{main}


begin

  reset(ein);
  read(ein,m,n); {Daten in ein : m n A[1,1]...A[1,n] A[2,1]......A[m,n]}
  main(m,n);

end.
  



7 4
1 3 4 2
5 1 2 6
4 2 1 3
2 6 5 1
4 2 1 3 
5 1 2 6
1 4 3 2





Berechnete QR-Zerlegung der eingelesenen Matrix A

Q[1,1] = -1.066003581778053E-001
Q[2,1] = -5.330017908890261E-001
Q[3,1] = -4.264014327112209E-001
Q[4,1] = -2.132007163556104E-001
Q[5,1] = -4.264014327112209E-001
Q[6,1] = -5.330017908890261E-001
Q[7,1] = -1.066003581778052E-001

Q[1,2] =  3.592462455458031E-001
Q[2,2] = -2.247339526930366E-001
Q[3,2] = -6.561575261110509E-003
Q[4,2] =  7.184924910916060E-001
Q[5,2] = -6.561575261110509E-003
Q[6,2] = -2.247339526930364E-001
Q[7,2] =  5.036009012902350E-001

Q[1,3] =  6.363806738297116E-001
Q[2,3] =  3.071969994754338E-001
Q[3,3] = -4.375759266775494E-001
Q[4,3] = -2.887035346240708E-002
Q[5,3] = -4.375759266775496E-001
Q[6,3] =  3.071969994754340E-001
Q[7,3] = -1.500025482388407E-001

Q[1,4] = -1.147026382966364E-001
Q[2,4] =  1.213566587503087E-001
Q[3,4] = -1.276938210871395E-001
Q[4,4] = -4.607117018875944E-001
Q[5,4] = -1.276938210871393E-001
Q[6,4] =  1.213566587503085E-001
Q[7,4] =  8.441100232658539E-001


R[1,1] = -9.380831519646859E+000

R[1,2] = -4.797016118001236E+000
R[2,2] =  6.927383081917469E+000

R[1,3] = -4.797016118001236E+000
R[2,3] =  5.628191180217581E+000
R[3,3] =  2.304799427836926E+000

R[1,4] = -9.594032236002471E+000
R[2,4] = -2.919900991194196E-001
R[3,4] =  2.004794331359245E+000
R[4,4] =  1.688220046531708E+000