title: LAPACK dgetrs vs dgesv
body: <p>In LAPACK documentation,</p>
  <p><strong>dgesv</strong> Solves a general system of linear equations AX=B.</p>
  <p><strong>dgetrs</strong> Solves a general system of linear equations AX=B, A<strong>T X=B or A</strong>H X=B, using the LU factorization computed by DGETRF.</p>
  <p>What is the difference ? Also, I made following C/C++ program and both give different result. Why is it so ?</p>
  <pre><code>#include &lt;iostream&gt;
  #include &lt;iomanip&gt;
  #include &lt;vector&gt;
  #include &lt;math.h&gt;
  #include &lt;time.h&gt;
  #include &quot;stdafx.h&quot;
  
  using namespace std;
  
  extern &quot;C&quot; {
  // LU decomoposition of a general matrix
  void dgetrf_(int* M, int *N, double* A, int* lda, int* IPIV, int* INFO);
  //// generate inverse of a matrix given its LU decomposition
  //void dgetri_(int* N, double* A, int* lda, int* IPIV, double* WORK, int*    lwork, int* INFO);
  void dgetrs_(char* C, int* N, int* NRHS, double* A, int* LDA, int* IPIV, double* B, int* LDB, int* INFO);
  void dgesv_(int *n, int *nrhs, double *a, int *lda, int *ipiv, double *b, int *ldb, int *info);
  }
  
  void solvelineq(double* A, double* B, int N)
  {
  int *IPIV = new int[N + 1];
  int LWORK = N*N;
  double *WORK = new double[LWORK];
  int INFO;
  char C = 'N';
  int NRHS = 1;
  dgetrf_(&amp;N, &amp;N, A, &amp;N, IPIV, &amp;INFO);
  dgetrs_(&amp;C, &amp;N, &amp;NRHS, A, &amp;N, IPIV, B, &amp;N, &amp;INFO);
  
  //dgetri_(&amp;N, A, &amp;N, IPIV, WORK, &amp;LWORK, &amp;INFO);
  
  delete IPIV;
  delete WORK;
  }
  
  double comparematrices(double* A, double* B, int N)
  {
  double sum = 0.;
  for (int i = 0; i &lt; N; ++i)
      sum += fabs(A[i] - B[i]);
  return sum;
  }
  
  int main() {
  int dim;
  std::cout &lt;&lt; &quot;Enter Dimensions of Matrix : \n&quot;;
  std::cin &gt;&gt; dim;
  clock_t c1, c2;
  c1 = clock();
  
  vector&lt;double&gt; a(dim*dim);
  vector&lt;double&gt; b(dim);
  vector&lt;double&gt; c(dim);
  vector&lt;int&gt; ipiv(dim);
  srand(1);              // seed the random # generator with a known value
  double maxr = (double)RAND_MAX;
  for (int r = 0; r &lt; dim; r++) {  // set a to a random matrix, i to the identity
      for (int c = 0; c &lt; dim; c++) {
          a[r + c*dim] = rand() / maxr;
      }
      b[r] = rand() / maxr;
      c[r] = b[r];
  }
  c2 = clock();
  cout &lt;&lt; endl &lt;&lt; dim &lt;&lt; &quot; x &quot; &lt;&lt; dim &lt;&lt; &quot; matrix generated in : &quot; &lt;&lt; double(c2 - c1) / CLK_TCK &lt;&lt; &quot; seconds &quot; &lt;&lt; endl;
  // C is identical matrix to B because we are solving by two methods.
  
  c1 = clock();
  solvelineq(&amp;*a.begin(), &amp;*c.begin(), dim);
  c2 = clock();
  cout &lt;&lt; endl &lt;&lt; &quot; dgetrf_ and dgetrs_ completed in : &quot; &lt;&lt; double(c2 - c1) / CLK_TCK &lt;&lt; &quot; seconds &quot; &lt;&lt; endl;
  vector&lt;double&gt; a1(a);
  vector&lt;double&gt; b1(b);
  int info;
  int one = 1;
  c1 = clock();
  dgesv_(&amp;dim, &amp;one, &amp;*a.begin(), &amp;dim, &amp;*ipiv.begin(), &amp;*b.begin(), &amp;dim, &amp;info);
  c2 = clock();
  cout &lt;&lt; endl &lt;&lt; &quot; dgesv_ completed in : &quot; &lt;&lt; double(c2 - c1) / CLK_TCK &lt;&lt; &quot; seconds &quot; &lt;&lt; endl;
  cout &lt;&lt; &quot;info is &quot; &lt;&lt; info &lt;&lt; endl;
  double eps = 0.;
  c1 = clock();
  for (int i = 0; i &lt; dim; ++i)
  {
      double sum = 0.;
      for (int j = 0; j &lt; dim; ++j)
          sum += a1[i + j*dim] * b[j];
      eps += fabs(b1[i] - sum);
  }
  c2 = clock();
  cout &lt;&lt; endl &lt;&lt; &quot; dgesv_ check completed in : &quot; &lt;&lt; double(c2 - c1) / CLK_TCK &lt;&lt; &quot; seconds &quot; &lt;&lt; endl;
  cout &lt;&lt; &quot;check is &quot; &lt;&lt; eps &lt;&lt; endl;
  cout &lt;&lt; &quot;\nFinal Matrix 1 : &quot; &lt;&lt; endl;
  for (int i = 0; i &lt; dim; ++i)
      cout &lt;&lt; b[i] &lt;&lt; endl;
  
  cout &lt;&lt; &quot;\nFinal Matrix 2 : &quot; &lt;&lt; endl;
  for (int i = 0; i &lt; dim; ++i)
      cout &lt;&lt; c[i] &lt;&lt; endl;
  double diff;
  c1 = clock();
  diff = comparematrices(&amp;*b.begin(), &amp;*c.begin(), dim);
  c2 = clock();
  cout &lt;&lt; endl &lt;&lt; &quot; Both functions compared in : &quot; &lt;&lt; double(c2 - c1) / CLK_TCK &lt;&lt; &quot; seconds &quot; &lt;&lt; endl;
  cout &lt;&lt; &quot;Sum of difference is : &quot; &lt;&lt; diff &lt;&lt; endl;
  return 0;
  }
  </code></pre>
  <p><em><strong>My Result :</strong></em>
  Enter Dimensions of Matrix :
  5</p>
  <p>5 x 5 matrix generated in : 0 seconds</p>
  <p>dgetrf_ and dgetrs_ completed in : 0.001 seconds</p>
  <p>dgesv_ completed in : 0 seconds
  info is 0</p>
  <p>dgesv_ check completed in : 0 seconds
  check is 2.02009e-15</p>
  <p>Final Matrix 1 :
  -2.97629
  4.83948
  2.00769
  -1.98887
  -5.15561</p>
  <p>Final Matrix 2 :
  -1.40622
  2.29029
  0.480829
  -1.63597
  0.71962</p>
  <p>Both functions compared in : 0 seconds</p>
  <p><strong>Sum of difference is : 11.8743</strong></p>