(*$B- Aufgabe 5.10B *)
PROGRAM qr ;
CONST nn=20;
TYPE vektor = ARRAY[1..nn] OF real;
     matrix = ARRAY[1..nn] OF vektor;
VAR a,aa: matrix; b,x:vektor; i,j,n:integer;

PROCEDURE eingabe(VAR a:matrix);
VAR i,j : integer;
BEGIN
  writeln('Matrix zeilenweise eingeben');
  FOR i:=1 TO n DO
  BEGIN
    writeln('Zeile Nr.',i:3);
    FOR j:=1 TO n DO read(a[i,j]);
  END
END;

PROCEDURE singulaer;
BEGIN
  writeln('Matrix ist singulaer');
END;

PROCEDURE givenselimination;
VAR i,j,k: integer; cot,co,si,h: real;
BEGIN
 FOR i:=1 TO n DO
 BEGIN
  FOR k := i+1 TO n DO
   IF a[k,i]<>0 THEN
   BEGIN
     cot:=a[i,i]/a[k,i];
     si:=1/sqrt(1+cot*cot); co:=si*cot;
     a[i,i]:=a[i,i]*co+a[k,i]*si;
     FOR j:=i+1 TO n DO
     BEGIN
       h:=a[i,j]*co+a[k,j]*si;
       a[k,j]:=-a[i,j]*si+a[k,j]*co;
       a[i,j]:=h
     END;
     h:=b[i]*co+b[k]*si;
     b[k]:=-b[i]*si+b[k]*co;
     b[i]:=h
   END;
   IF a[i,i]=0 THEN singulaer
  END
END;

BEGIN
 writeln('Ordnung der Matrix');read(n);
 eingabe(aa);
 writeln('Matrix Q =');
 FOR i := 1 TO n DO
 BEGIN
   a := aa;
   (* b = i-ter Einheitsvektor  *)
   FOR j := 1 TO n DO b[j]:=0; b[i]:=1;
   givenselimination;
   (* b enthaelt jetzt die i-te Kolonnen von Q^T
      d.h. die i-te Zeile von Q  *)
   FOR j := 1 TO n DO write(b[j]:10:5); writeln
 END;
 writeln('Matrix R');
 FOR i :=1 TO n DO
 BEGIN
   FOR j := 1 TO i-1 DO write(' ':10);
   FOR j := i TO n DO write(a[i,j]:10:5);
   writeln
 END
END.