title: Solving a linear system with Lapack's dgeqrf_
body: <p>I am trying to factorize a matrix with the  <a href="http://en.wikipedia.org/wiki/QR_decomposition" rel="nofollow">QR factorization</a> in C++, using Lapack's functions in order to solve a system of linear equations (Ax=b)</p>
  
  <p>As far as I understood, <em>dgeqrf</em> computes the QR factorization and overwrites the input matrix. The output clearly contains values for L (upper triangular), but how do I obtain Q?</p>
  
  <p>I tried <em>dormqr</em>, which is said to calculate Q from <em>dgeqrf</em>'s output, but the result is the same matrix as in the previous call.</p>
  
  <p>Here's my complete code:</p>
  
  <pre><code>boost::numeric::ublas::matrix&lt;double&gt; in_A(4, 3);
  in_A(0, 0) = 1.0;
  in_A(0, 1) = 2.0;
  in_A(0, 2) = 3.0;
  
  in_A(1, 1) = -3.0;
  in_A(1, 2) = 2.0;
  in_A(1, 3) = 1.0;
  
  in_A(2, 1) = 2.0;
  in_A(2, 2) = 0.0;
  in_A(2, 3) = -1.0;
  
  in_A(3, 1) = 3.0;
  in_A(3, 2) = -1.0;
  in_A(3, 3) = 2.0;
  
  boost::numeric::ublas::vector&lt;double&gt; in_b(4);
  in_b(0) = 2;
  in_b(1) = 4;
  in_b(2) = 6;
  in_b(3) = 8;
  
  int rows = in_A.size1();
  int cols = in_A.size2();
  double *A = (double *)malloc(rows*cols*sizeof(double));
  double *b = (double *)malloc(in_b.size()*sizeof(double));
  
  //Lapack has column-major order
  for(size_t col=0; col&lt;in_A.size2(); ++col)
  {
      for(size_t row = 0; row&lt;in_A.size1(); ++row)
  {
      int D1_idx = col*in_A.size1() + row;
      A[D1_idx] = in_A(row, col);
  }
  b[col] = in_b(col);
  }
  
  integer m = rows;
  integer n = cols;
  
  integer info = 0;
  integer k = n;          /* k = min(m,n);       */
  integer lda = m;        /* lda = max(m,1);     */
  integer lwork = n;      /* lwork = max(n,1);   */
  int max = lwork;    /* max = max(lwork,1); */
  
  double *work;
  double *tau;
  
  char *side = "L";
  char *TR    = "T";
  integer one = 1;
  int i;
  
  double *vec;
  
  work = (double *) malloc( max * sizeof( double ) );
  tau  = (double *) malloc( k * sizeof( double ) );
  vec  = (double *) malloc( m * sizeof( double ) );
  
  memset(work, 0, max * sizeof(double));
  memset(tau, 0, k * sizeof(double));
  std::cout &lt;&lt; std::endl;
  for(size_t row = 0; row &lt; rows; ++row)
  {
  for(size_t col = 0; col &lt; cols; ++col)
  {
  size_t idx = col*rows + row;
  std::cout &lt;&lt; A[idx] &lt;&lt; " ";
  }
  std::cout &lt;&lt; std::endl;
  }
  dgeqrf_(&amp;m, &amp;n, A, &amp;lda, tau, work, &amp;lwork, &amp;info);
  //printf("tau[0] = %f tau[1] = %f\n",tau[0],tau[1]);
  
  std::cout &lt;&lt; std::endl;
  for(size_t row = 0; row &lt; rows; ++row)
  {
    for(size_t col = 0; col &lt; cols; ++col)
    {
    size_t idx = col*rows + row;
    std::cout &lt;&lt; A[idx] &lt;&lt; " ";
    }
  std::cout &lt;&lt; std::endl;
  }
  
  memset(vec, 0, m * sizeof(double));
  vec[2] = 1.0;
  
  dormqr_(side, TR, &amp;m, &amp;one, &amp;k, A, &amp;lda, tau, vec, &amp;lda, work, &amp;lwork, &amp;info);
  
  free(vec);
  free(tau);
  free(work);
  </code></pre>
  
  <p>What's wrong with my code?</p>
  
  <p>How can I factorize a matrix and solve a corresponding system of linear equations?</p>