#=================================================#
# Dies ist der Maple source zur Q-R-Zerlegung     #
# mit Givensrotationen                            #
#                                                 #
# Erwartet wird, dass die Dimmension m, n, die    # 
# (m x n)-Matrix A schon deklariert sind.         #
#=================================================#

Q:=evalm(array(identity,1..m,1..m));
R:=evalm(A);

for j from 1 to n do 
  for id from 1 to 1 do
    a:=R[j,j];
    aa:=abs(a);
  od;
  for i from j+1 to m do
    b:=R[i,j]; 
    ab:=abs(b);
    if (ab > 0) then # hier muss eps hin!!
      a:=R[j,j];
      if (ab > abs(a)) then
        tau:=a/b; s:=1/sqrt(1+tau^2); c:=s*tau;
printf(`Spalte %d, Zeile %d: tau=%11.5f, s=%11.5f, c=%11.5f\n`, j,i, tau, s,c);
      else
        tau:=b/a; c:=1/sqrt(1+tau^2); s:=c*tau;
printf(`Spalte %d, Zeile %d: tau=%11.5f, c=%11.5f, s=%11.5f\n`, j,i, tau, c,s);
      fi;
      # update R
      for k from j to n do
        h:=R[j,k];
	R[j,k]:= c*h+s*R[i,k];
	R[i,k]:=-s*h+c*R[i,k];
      od;
printf(`R: %d. und %d. Zeile\n`, j,i);
for k from j to n do printf(`%11.5f`, R[j,k]); od; printf(`\n`):
for k from j to n do printf(`%11.5f`, R[i,k]); od; printf(`\n`):
      # now update Q'
      for k from 1 to m do
        h:=Q[j,k];
	Q[j,k]:= c*h+s*Q[i,k];
	Q[i,k]:=-s*h+c*Q[i,k];
      od;
printf(`Q': %d. und %d. Zeile\n`, j,i);
for k from 1 to m do printf(`%11.5f`, Q[j,k]); od; printf(`\n`):
for k from 1 to m do printf(`%11.5f`, Q[i,k]); od; printf(`\n`):
    fi;  

  od;
od;

R:=evalm(R);
Q:=evalm(transpose(Q));
