LAPACK dgetrs vs dgesv

Open discussion regarding features, bugs, issues, vendors, etc.

LAPACK dgetrs vs dgesv

Postby mailmaverick » Thu Mar 17, 2016 11:47 am

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 ?

Code: Select all
#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 :
Code: Select all
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
mailmaverick
 
Posts: 2
Joined: Thu Mar 17, 2016 11:37 am

Re: LAPACK dgetrs vs dgesv

Postby CyLith » Mon Mar 21, 2016 2:45 am

If you look at the code for DGESV, it does nothing more than call DGETRF and then DGETRS:

http://www.netlib.org/lapack/explore-3.1.1-html/dgesv.f.html

The reason you get different results is because DGETRF replaces the matrix with the LU factorization in-place, so if you call DGESV afterwards, it will think the LU factorization is your original matrix, and give you something completely different. You need to save the original matrix before calling
Code: Select all
solvelineq
in order to do a comparison. Also, pet peeve: you need to use "delete []" instead of just "delete" on arrays.
CyLith
 
Posts: 40
Joined: Sun Feb 08, 2009 7:23 am
Location: Stanford, CA


Return to User Discussion

Who is online

Users browsing this forum: Bing [Bot] and 3 guests