title: QR Decomposition Algorithm Using Givens Rotations
body: <p>I am coding a QR decomposition algorithm in MATLAB, just to make sure I have the mechanics correct.  Here is the code for the main function:</p>
  
  <pre><code>    function [Q,R] = QRgivens(A)
          n = length(A(:,1));
          Q = eye(n);
          R = A;
  
          for j = 1:(n-1)
              for i = n:(-1):(j+1)
                  G = eye(n);
                  [c,s] = GivensRotation( A(i-1,j),A(i,j) );
                  G(i-1,(i-1):i) = [c s];
                  G(i,(i-1):i)   = [-s c];
                  Q = Q*G';
                  R = G*R;
              end
          end
      end
  </code></pre>
  
  <p>The sub function GivensRotation is given below:</p>
  
  <pre><code>    function [c,s] = GivensRotation(a,b)
          if b == 0
              c = 1;
              s = 0;
          else
              if abs(b) &gt; abs(a)
                  r = -a / b;
                  s = 1 / sqrt(1 + r^2);
                  c = s*r;
              else
                  r = -b / a;
                  c = 1 / sqrt(1 + r^2);
                  s = c*r;
              end
          end
      end
  </code></pre>
  
  <p>I've done research and I'm pretty sure this is one of the most straightforward ways to implement this decomposition, in MATLAB especially.  But when I test it on a matrix A, the R produced is not right triangular as it should be.  The Q is orthogonal, and Q*R = A, so the algorithm is doing some things right, but it is not producing exactly the correct factorization.  Perhaps I've just been staring at the problem too long, but any insight as to what I've overlooked would be appreciated.</p>