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

R:=evalm(A):


for jGRj from 1 to n-1 do 
  a:=R[jGRj,jGRj];
  aa:=abs(a);
  for iGRi from jGRj+1 to m do
    b:=R[iGRi,jGRj]; 
    ab:=abs(b);
    if (ab > 0) then # hier muss eps hin!!
      a:=R[jGRj,jGRj];
      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`, jGRj,iGRi, 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`, jGRj,iGRi, tau, c,s);
      fi;
      for k from jGRj to n do
        h:=R[jGRj,k];
	R[jGRj,k]:= c*h+s*R[iGRi,k];
	R[iGRi,k]:=-s*h+c*R[iGRi,k];
      od;
# printf(`R: %d. und %d. Zeile\n`, jGRj,iGRi);
# for k from jGRj to n do printf(`%11.5f`, R[jGRj,k]); od; printf(`\n`):
# for k from jGRj to n do printf(`%11.5f`, R[iGRi,k]); od; printf(`\n`):
    fi;  
  od:
od:

R:=evalm(R);
