Beim Givens-Verfahren wird eine QR-Zerlegung

  A = Q R

einer Matrix A zur Lösung eines linearen Gleichungssystems verwendet. Q ist dabei eine orthogonale Matrix (das Produkt von Q mit dem transponierten Q ergibt die Einheitsmatrix) und R eine obere Dreiecksmatrix. Das Gleichungssystem

  A x = b

ist dann äquivalent zu

  R x = transponiert(Q) b.

Als Pseudo Code sieht die Realisierung des Verfahrens wie folgt aus:

(*
Löst Ax=y nach x und verändert A,x und y
wobei A eine n mal m Matrix,
x ein n dimensionaler
und y ein m dimensionaler Vektor ist.
*)
procedure givensInverse(m,n:integer;
                        var A:Matrix;
                        var x,y:Vector;
                        var ok:boolean);
(* QR-Zerlegung nach Givens *)
var i,j,k:integer;
    c,s,
    aki,aii,gamma,aij,akj,
    yi,yk:extended;
    summe,aik:extended;
begin
  (* QR-Zerlegung nach Givens *)
  for i:=1 to min(n, pred(m)) do begin
    for k:=succ(i) to m do begin
      aki:=A[k,i];
      if aki <> 0.0 then begin
        (* A[k,i] ungleich 0 -> zu 0 machen *)
        aii:=A[i,i];
        rot(aii, aki, c, s, gamma);
        A[i,i]:=gamma;
        A[k,i]:=0.0;
        (* i-te und k-te Zeile werden verändert *)
        for j:=succ(i) to n do begin
          aij:=A[i,j];
          akj:=A[k,j];
          A[i,j]:=c * aij + s * akj;
          A[k,j]:=c * akj - s * aij;
        end;
        yi:=y[i];
        yk:=y[k];
        y[i]:=c * yi + s * yk;
        y[k]:=c * yk - s * yi;
      end; (* if *)
    end; (* for *)
  end; (* for *)
  (* nun liegt in A eine obere Dreiecksmatrix R = QT A vor
     und y wurde mit QT multipliziert *)
  (* Rückwärtseinsetzen *)
  ok:=true;
  for i:=n downto 1 do begin
    if ok then begin
      (* summiere die bekannten Produkte der Zeile auf *)
      summe:=0.0;
      for k:=succ(i) to n do begin
        aik:=A[i,k];
        summe:=summe + x[k] * aik;
      end;
      aii:=A[i,i];
      (* wenn a[i,i] = 0 -> Matrix ist singulär *)
      ok:=(aii <> 0.0);
      if ok then begin
        (* bestimme xi *)
        x[i]:=(y[i] - summe) / aii;
      end;
    end;
  end;
end;

(* Berechnet die Parameter für die Givens-Rotation *)
procedure rot(const alpha,beta:extended;
              var c,s,gamma:extended);
var t,u:extended;
begin
  if beta = 0.0 then begin
    c:=1.0;
    s:=0.0;
    gamma:=alpha;
  end else if abs(beta) >= abs(alpha) then begin
    t:=alpha / beta;
    u:=sqrt(1.0 + sqr(t));
    s:=1.0 / u;
    c:=s * t;
    gamma:=beta * u;
  end else begin
    t:=beta / alpha;
    u:=sqrt(1.0 + sqr(t));
    c:=1.0 / u;
    s:=c * t;
    gamma:=alpha * u;
  end;
end;

(* Minimum zweier Ganzzahlen *)
function min(a,b:integer):integer;
begin
  if a < b then min:=a else min:=b;
end;