Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

## problem using c++ version of dgesv Beginner
251 Views
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 3int main(int argc, char* argv[]){	double a[NI][NI];	double b[NI][NJ];	a = 1;	a = 1;	a = 1;	a = 1;	a = 7.4895;	a = 7.5292;	a = 7.3905;	a = 7.4141;	a = 1.1759;	a = 1.2956;	a = 1.2342;	a = 1.3149;	a = 0.1496;	a = 0.3214;	a = 0.2887;	a = 0.1598;	b = 0;	b = 0;	b = 0;	b = 1;	b = 0;	b = 0;	b = 0;	b = 1;	b = 0;	b = 0;	b = 0;	b = 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
3 Replies Black Belt
251 Views
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. Moderator
251 Views
please see the DGESV C Example here. Beginner
251 Views
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#...
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_dg...

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 = 1;
a = 1;
a = 1;
a = 1;
a = 7.4895;
a = 7.5292;
a = 7.3905;
a = 7.4141;
a = 1.1759;
a = 1.2956;
a = 1.2342;
a = 1.3149;
a = 0.1496;
a = 0.3214;
a = 0.2887;
a = 0.1598;

b = 0;
b = 0;
b = 0;
b = 1;
b = 0;
b = 0;
b = 0;
b = 1;
b = 0;
b = 0;
b = 0;
b = 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 