title: QR decomposition to solve linear systems in CUDA
body: <p>I'm writing an image restoration algorithm on GPU, details in </p>
  
  <p><a href="https://stackoverflow.com/questions/22284533/cuda-least-square-solving-poor-in-speed">Cuda: least square solving , poor in speed</a></p>
  
  <p>The QR decomposition method to solve the linear system</p>
  
  <pre><code>Ax=b  
  </code></pre>
  
  <p>works as follows</p>
  
  <pre><code>min||Ax-b|| ---&gt; ||QRx-b||  ---&gt; ||(Q^T)QRx-(Q^T)b|| ---&gt; ||Rx-(Q^T)b||
  </code></pre>
  
  <p>where <code>R</code> is the upper triangular matrix. The resulting upper triangular linear system is easy to solve.</p>
  
  <p>I want to use CULA tools to implement this method. The CULA routine <code>GEQRF</code> computes a QR factorization. The manual says:</p>
  
  <blockquote>
    <p>On exit, the elements on and above the diagonal of the array contain
    the <code>min(M,N)-by-N</code> upper trapezoidal matrix <code>R</code> (<code>R</code> is upper
    triangular if <code>m &gt;= n)</code>; the elements below the diagonal, with the
    array <code>TAU</code>, represent the orthogonal/unitary matrix <code>Q</code> as a product
    of <code>min(m,n)</code> elementary reflectors.</p>
  </blockquote>
  
  <p>I cannot figure out where <code>Q</code> is stored, and the algorithm seems too complex for me. Could you give any advice?</p>
  
  <p>Thanks!</p>