- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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"And here is the (correct) Fortran version:

#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]

[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

Link Copied

3 Replies

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

Topic Options

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page