problem using c++ version of dgesv

[geowolf]

Dear all,

please find below two times the same sample code.
Purpose is to to solve a system of linear equations of matrix A and matrix B. However, the resulting Matrix B differs.
The Fortran version gives the correct, the C++ a wrong result.
Can anyone tell me what is wrong with the C++ source?

This is the (incorrect) C++ version:

[cpp]#include "mkl.h"
#include 
#include 
using namespace std;

#define NI 4
#define NJ 3


int main(int argc, char* argv[])
{
	double a[NI][NI];
	double b[NI][NJ];

	a[0][0] = 1;
	a[0][1] = 1;
	a[0][2] = 1;
	a[0][3] = 1;
	a[1][0] = 7.4895;
	a[1][1] = 7.5292;
	a[1][2] = 7.3905;
	a[1][3] = 7.4141;
	a[2][0] = 1.1759;
	a[2][1] = 1.2956;
	a[2][2] = 1.2342;
	a[2][3] = 1.3149;
	a[3][0] = 0.1496;
	a[3][1] = 0.3214;
	a[3][2] = 0.2887;
	a[3][3] = 0.1598;

	b[0][0] = 0;
	b[0][1] = 0;
	b[0][2] = 0;
	b[1][0] = 1;
	b[1][1] = 0;
	b[1][2] = 0;
	b[2][0] = 0;
	b[2][1] = 1;
	b[2][2] = 0;
	b[3][0] = 0;
	b[3][1] = 0;
	b[3][2] = 1;

	cout << "A1 = \n";
	for(int i = 0; i < NI; i++)
	{
		for(int j = 0; j < NI; j++)
		{
			cout << a << "\t";
		}
		cout << "\n";
	}
	cout << "\n";

	int ipiv[NI];
	for (int i = 0; i < NI; i++){
		ipiv = 0;
	}

	cout << "B1 = \n";
	for(int i = 0; i < NI; i++)
	{
		for(int j = 0; j < NJ; j++)
		{
			cout << b << "\t";
		}
		cout << "\n";
	}
	cout << "\n";
	int n = NI;
	int nrhs = NJ;
	int lda = NI;
	int ldb = NI;
	int info = -1;
	dgesv(&n, &nrhs, *a, &lda, ipiv, *b, &ldb, &info);

	cout << "INFO = " << info << "\n";

	cout << "G = \n";
	cout.setf(std::ios::scientific);
	for(int i = 0; i < NI; i++)
	{
		for(int j = 0; j < NJ; j++)
		{
			cout << b << "\t";
		}
		cout << "\n";
	}
	cout << "\n";

	return 0;
}

[/cpp]

And here is the (correct) Fortran version:

[fxfortran]      PROGRAMM TEST
      IMPLICIT NONE    

      DOUBLE PRECISION A(4,4), B(4,3)
      INTEGER N,NRHS,LDA,INFO,LDB,I,J
      DOUBLE PRECISION IPIV(4)
      
      N = 4
      NRHS = 3
      LDA = 4
      LDB = 4
      
 
      A(1,1) = 1
      A(1,2) = 1
      A(1,3) = 1
      A(1,4) = 1
      A(2,1) = 7.4895  
      A(2,2) = 7.5292  
      A(2,3) = 7.3905  
      A(2,4) = 7.4141
      A(3,1) = 1.1759  
      A(3,2) = 1.2956  
      A(3,3) = 1.2342  
      A(3,4) = 1.3149
      A(4,1) = 0.1496  
      A(4,2) = 0.3214  
      A(4,3) = 0.2887  
      A(4,4) = 0.1598

      B(1,1) = 0
      B(1,2) = 0
      B(1,3) = 0
      B(2,1) = 1
      B(2,2) = 0
      B(2,3) = 0
      B(3,1) = 0
      B(3,2) = 1
      B(3,3) = 0
      B(4,1) = 0
      B(4,2) = 0
      B(4,3) = 1

      DO I = 1, 4
      IPIV(I) = 0.0
      WRITE(*,'(4(F10.4))') (A(I,J), J = 1, 4)
      ENDDO
      
      DO I = 1, 4
      WRITE(*,'(4(F10.4))') (B(I,J), J = 1, 3)
      ENDDO
      
      CALL DGESV(N,NRHS,A,LDA,IPIV,B,LDB,INFO)
      
      WRITE(*,*) 'INFO = ', INFO
      
      DO I = 1, 4
      WRITE(*,'(4(F10.4))') (B(I,J), J = 1, 3)
      ENDDO
      
      END
[/fxfortran]

Kind regards
Wolfram

[mecej4]

C stores a matrix with row-1 first, then row-2,...

Fortran stores a matrix with col-1 first, then col-2,...

Unless your matrix is symmetric, you need to take these facts into consideration. Alternatively, use the C-Lapack interface which gives you wrappers around the Fortran routines with the ability to specify the storage order of the matrices, etc.

[Gennady_F_Intel]

please see the DGESV C Example here.

[geowolf]

Dear all,

the respective discrete links to the mentioned examples are:
http://software.intel.com/sites/products/documentation/hpc/mkl/lapack/mkl_lapack_examples/dgesv.htm#top
and for the version which is the one which applies here (row major) it is
http://software.intel.com/sites/products/documentation/hpc/mkl/lapack/mkl_lapack_examples/lapacke_dgesv_row.c.htm

Please find below a corrected version of the C++ source which gives a result identical to the Fortran version.
The main differences are the usage of:
info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, n, nrhs, *a, lda, ipiv, *b, ldb);
and a fix concerning the size of ldb which has to be identical to nrhs (i.e. NJ). (The size of a is [NI][NI] while the size of b is [NI][NJ]!)

[cpp]#include "mkl.h"
#include "mkl_lapacke.h"
#include 
#include 
using namespace std;
#define NI 4 
#define NJ 3 

int main(int argc, char* argv[]) 
{ 
	double a[NI][NI]; 
	double b[NI][NJ]; 

	a[0][0] = 1; 
	a[0][1] = 1; 
	a[0][2] = 1; 
	a[0][3] = 1; 
	a[1][0] = 7.4895; 
	a[1][1] = 7.5292; 
	a[1][2] = 7.3905; 
	a[1][3] = 7.4141; 
	a[2][0] = 1.1759; 
	a[2][1] = 1.2956; 
	a[2][2] = 1.2342; 
	a[2][3] = 1.3149; 
	a[3][0] = 0.1496; 
	a[3][1] = 0.3214; 
	a[3][2] = 0.2887; 
	a[3][3] = 0.1598; 

	b[0][0] = 0; 
	b[0][1] = 0; 
	b[0][2] = 0; 
	b[1][0] = 1; 
	b[1][1] = 0; 
	b[1][2] = 0; 
	b[2][0] = 0; 
	b[2][1] = 1; 
	b[2][2] = 0; 
	b[3][0] = 0; 
	b[3][1] = 0; 
	b[3][2] = 1; 

	cout << "A1 = n"; 
	for(int i = 0; i < NI; i++) { 
		for(int j = 0; j < NI; j++) { 
			cout << a << "t"; 
		} 
		cout << "n"; 
	} 
	cout << "n"; 

	MKL_INT ipiv[NI]; 
	for (int i = 0; i < NI; i++){ 
		ipiv = 0; 
	} 

	cout << "B1 = n"; 
	for(int i = 0; i < NI; i++) { 
		for(int j = 0; j < NJ; j++) { 
			cout << b << "t"; 
		} 
		cout << "n"; 
	} 
	cout << "n"; 
	MKL_INT n = NI;
	MKL_INT nrhs = NJ;
	MKL_INT lda = NI;
	MKL_INT ldb = NJ;
	MKL_INT info = -1;
	info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, n, nrhs, *a, lda, ipiv, *b, ldb);

	cout << "INFO = " << info << "n"; 

	cout << "G = n"; 
	cout.setf(std::ios::scientific); 
	for(int i = 0; i < NI; i++) { 
		for(int j = 0; j < NJ; j++) { 
			cout << b << "t"; 
		} 
		cout << "n"; 
	} 
	cout << "n"; 

	return 0; 
} 
[/cpp]

Kind regards
Wolfram