LAPACK dgetrs vs dgesv

Question

In LAPACK documentation,

dgesv Solves a general system of linear equations AX=B.

dgetrs Solves a general system of linear equations AX=B, AT X=B or AH X=B, using the LU factorization computed by DGETRF.

What is the difference ? Also, I made following C/C++ program and both give different result. Why is it so ?

#include <iostream>
#include <iomanip>
#include <vector>
#include <math.h>
#include <time.h>
#include "stdafx.h"

using namespace std;

extern "C" {
// 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_(&N, &N, A, &N, IPIV, &INFO);
dgetrs_(&C, &N, &NRHS, A, &N, IPIV, B, &N, &INFO);

//dgetri_(&N, A, &N, IPIV, WORK, &LWORK, &INFO);

delete IPIV;
delete WORK;
}

double comparematrices(double* A, double* B, int N)
{
double sum = 0.;
for (int i = 0; i < N; ++i)
    sum += fabs(A[i] - B[i]);
return sum;
}

int main() {
int dim;
std::cout << "Enter Dimensions of Matrix : \n";
std::cin >> dim;
clock_t c1, c2;
c1 = clock();

vector<double> a(dim*dim);
vector<double> b(dim);
vector<double> c(dim);
vector<int> ipiv(dim);
srand(1);              // seed the random # generator with a known value
double maxr = (double)RAND_MAX;
for (int r = 0; r < dim; r++) {  // set a to a random matrix, i to the identity
    for (int c = 0; c < dim; c++) {
        a[r + c*dim] = rand() / maxr;
    }
    b[r] = rand() / maxr;
    c[r] = b[r];
}
c2 = clock();
cout << endl << dim << " x " << dim << " matrix generated in : " << double(c2 - c1) / CLK_TCK << " seconds " << endl;
// C is identical matrix to B because we are solving by two methods.

c1 = clock();
solvelineq(&*a.begin(), &*c.begin(), dim);
c2 = clock();
cout << endl << " dgetrf_ and dgetrs_ completed in : " << double(c2 - c1) / CLK_TCK << " seconds " << endl;
vector<double> a1(a);
vector<double> b1(b);
int info;
int one = 1;
c1 = clock();
dgesv_(&dim, &one, &*a.begin(), &dim, &*ipiv.begin(), &*b.begin(), &dim, &info);
c2 = clock();
cout << endl << " dgesv_ completed in : " << double(c2 - c1) / CLK_TCK << " seconds " << endl;
cout << "info is " << info << endl;
double eps = 0.;
c1 = clock();
for (int i = 0; i < dim; ++i)
{
    double sum = 0.;
    for (int j = 0; j < dim; ++j)
        sum += a1[i + j*dim] * b[j];
    eps += fabs(b1[i] - sum);
}
c2 = clock();
cout << endl << " dgesv_ check completed in : " << double(c2 - c1) / CLK_TCK << " seconds " << endl;
cout << "check is " << eps << endl;
cout << "\nFinal Matrix 1 : " << endl;
for (int i = 0; i < dim; ++i)
    cout << b[i] << endl;

cout << "\nFinal Matrix 2 : " << endl;
for (int i = 0; i < dim; ++i)
    cout << c[i] << endl;
double diff;
c1 = clock();
diff = comparematrices(&*b.begin(), &*c.begin(), dim);
c2 = clock();
cout << endl << " Both functions compared in : " << double(c2 - c1) / CLK_TCK << " seconds " << endl;
cout << "Sum of difference is : " << diff << endl;
return 0;
}

My Result : Enter Dimensions of Matrix : 5

5 x 5 matrix generated in : 0 seconds

dgetrf_ and dgetrs_ completed in : 0.001 seconds

dgesv_ completed in : 0 seconds info is 0

dgesv_ check completed in : 0 seconds check is 2.02009e-15

Final Matrix 1 : -2.97629 4.83948 2.00769 -1.98887 -5.15561

Final Matrix 2 : -1.40622 2.29029 0.480829 -1.63597 0.71962

Both functions compared in : 0 seconds

Sum of difference is : 11.8743

score 2 · Answer 1 · answered Mar 17 '16 at 15:53

Calling dgesv is equivalent of calling dgetrf and dgetrs. The source code of dgesv is quite simple. It mainly contains the calls to dgetrf and dgetrs functions.

The results are different in the example because dgetrf changes the A matrix. In lapack's reference it states:

A is DOUBLE PRECISION array, dimension (LDA,N)
On entry, the M-by-N matrix to be factored.
On exit, the factors L and U from the factorization
A = P*L*U; the unit diagonal elements of L are not stored.

In your example the second time is used a different A matrix. Since you are planning to use lapack, I suggest you also try to use blas routines also (daxpy,dnrm2).

I created a simple example comparing dgesv, dgetrf and dgetrs

#include <iostream>
#include <vector>
#include <stdlib.h>

using namespace std;

extern "C" {
  void daxpy_(int* n,double* alpha,double* dx,int* incx,double* dy,int* incy);
  double dnrm2_(int* n,double* x, int* incx);

  void dgetrf_(int* M, int *N, double* A, int* lda, int* IPIV, 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 print(const char* head, int m, int n, double* a){
  cout << endl << head << endl << "---------" << endl;
  for(int i=0; i<m; ++i){
    for(int j=0; j<n; ++j){
      cout << a[i+j*m]<< ' ';
    }
    cout << endl;
  }
  cout << "---------" << endl;
}

int main() {
  int dim;
  std::cout << "Enter Dimensions of Matrix : \n";
  std::cin >> dim;

  vector<double> a(dim*dim);
  vector<double> b(dim);
  srand(1);              // seed the random # generator with a known value
  double maxr = (double)RAND_MAX;
  for (int r = 0; r < dim; r++) {  // set a to a random matrix, i to the identity
    for (int i = 0; i < dim; i++) {
      a[r + i*dim] = rand() / maxr;
    }
    b[r] = rand() / maxr;
  }

  int info;
  int one = 1;
  char N = 'N';
  vector<int> ipiv(dim);

  vector<double> a1(a);
  vector<double> b1(b);
  dgesv_(&dim, &one, a1.data(), &dim, ipiv.data(), b1.data(), &dim, &info);
  print("B1",5,1,b1.data());

  vector<double> a2(a);
  vector<double> b2(b);
  dgetrf_(&dim, &dim, a2.data(), &dim, ipiv.data(), &info);
  dgetrs_(&N, &dim, &one, a2.data(), &dim, ipiv.data(), b2.data(), &dim, &info);
  print("B2",5,1,b2.data());

  double d_m_one = -1e0;
  daxpy_(&dim,&d_m_one,b2.data(),&one,b1.data(),&one);
  cout << endl << "diff = " << dnrm2_(&dim,b1.data(),&one) << endl;

  return 0;
}

In my PC running the dim=5 case it gave:

B1
---------
0.782902 
3.35317 
4.40269 
-5.10512 
-1.26615 
---------

B2
---------
0.782902 
3.35317 
4.40269 
-5.10512 
-1.26615 
---------

diff = 0

LAPACK dgetrs vs dgesv

1 Answers1