Problems with mkl_dcootrsv

[Sören_Plönnigs]

Hi,

I used the MKL to solve the following problem:

Matrix:

row column value

0 0 1

1 0 -1

1 1 1.73205

2 1 -0.57735

2 2 1.29099

3 1 -0.57735

3 3 1.29099

4 1 -0.57735

4 4 1.91485

5 2 -0.774597

5 4 -0.522233

5 5 1.76841

6 3 -0.774597

6 4 -0.522233

6 6 1.76841

7 7 1

8 8 1

9 4 -0.522233

9 9 1.93061

10 6 -0.56548

10 9 -0.51797

10 10 1.18825

11 5 -0.56548

11 9 -0.51797

11 11 1.18825

12 9 -0.51797

12 12 0.855399

rhs:

0 0 0 0 0 -2 0 -2 0 0 0 0 0

why do I get

3.35965e-322 1.4822e-322 1.4822e-322 9.88131e-323 9.38725e-323 -1.13096 0 -2 2.47033e-323 -0.538215 5.92879e-323 1.97626e-323 0

as the result instead of

0 0 0 0 0 -1.1364 0 -2 0 0 0 0 0

which would be the exact result?

Sren

[Ying_H_Intel]

Hello Soren,

could you please specify what version ofMKL do you use, what is the target cpu (32bit or X64), what is your OS (windows, linux or MacOS), simple test case would help to investigate that problem.

Regards,
Ying

[barragan_villanueva_]

Hi,

I can see that you got approximatly the same result. The most difference is about 0.006.
Do you expect more accuracy?

[Sören_Plönnigs]

Hello Ying,

I'm using MKL10.2.5.035. The Problem occurs on ia32 as well as intel64. my OS is linux.

I wrote this little testprogram for the above problem:

matrixtest.cpp:

#include "mkl_spblas.h"

#include "mkl_types.h"

#include

using namespace std;

int main() {

int m=13;

int nz=27;

char uplo = 'l';

char transa = 'n';

char diag = 'n';

// The Matrix

int rowind[] = {0,1,1,2,2,3,3,4,4,5,5,5,6,6,6,7,8,9,9,10,10,10,11,11,11,12,12};

int colind[] = {0,0,1,1,2,1,3,1,4,2,4,5,3,4,6,7,8,4,9, 6, 9,10, 5, 9,11, 9,12};

double value[] = {1,-1,1.73205,-0.57735,1.29099,-0.57735,1.29099,-0.57735,1.91485,-0.774597,-0.522233,1.76841,

-0.774597,-0.522233,1.76841,1,1,-0.522233,1.93061,-0.56548,-0.51797,1.18825,-0.56548,-0.51797,1.18825,-0.51797,0.855399};

// The Vectors

double x[] = { 0,0,0,0,0,2,0,2,0,0,0,0,0 };

double *y = new double();

mkl_dcootrsv(&uplo,&transa,&diag,&m,value,rowind,colind,&nz,x,y);

for (int i=0;i

cout << y << " ";

cout << endl;

return 0;

}

With this program I found out that the inaccuracy results from the debugging flag.

output with debugging flag:

icpc -L${MKLPATH} ${MKLPATH}/libmkl_solver.a -Wl,--start-group -lmkl_intel -lmkl_intel_thread -lmkl_core -Wl,--end-group -openmp -lpthread matrixtest.cpp -g -o matrixtest

./matrixtest

1.3844e-312 6.19126e-313 6.19126e-313 4.17414e-313 3.94456e-313 1.13096 0 2 1.12911e-313 0.538216 2.36938e-313 6.83711e-314 0

output without debugging flag:

icpc -L${MKLPATH} ${MKLPATH}/libmkl_solver.a -Wl,--start-group -lmkl_intel -lmkl_intel_thread -lmkl_core -Wl,--end-group -openmp -lpthread matrixtest.cpp -o matrixtest

./matrixtest

0 0 0 0 0 1.13096 0 2 0 0.538216 0 0 0

What can I do to avoid this problem?

Regards,

Sren

[Sören_Plönnigs]

Hi Victor,

I need much more accuracy. I need this result for further calculations, and the error will rise with each iteration (a few thousand iterations). The result of my program has to reach an accuracy of e-12, which is not possible with this inaccuracy.

Please read my above posting to Ying to see the result of my further investigations.

Sren

[Sören_Plönnigs]

I could resolv the error with the debug flag. The problem was, that I didn't read the mkl_dcootrsvmanual careful enough, so I didn't realized that one-based indexing is needed for it to work.

Now I got the next problem. I want to use the the solver for solving thq Equation M*M'*y=x. The first step, solving M*z=x, works now. But the second step, solving M'*y=z, does not work. The Matrix is the same as described above. The vector x is the result vector from the above test program (0 0 0 0 0 1.13096 0 2 0 0.538216 0 0 0). Here is my little test program:

#include "mkl_spblas.h"

#include "mkl_types.h"

#include

using namespace std;

int main() {

int m=13;

int nz=27;

char uplo = 'l';

char transa = 'n';

char diag = 'n';

// The Matrix

int rowind[] = {0,1,1,2,2,3,3,4,4,5,5,5,6,6,6,7,8,9,9,10,10,10,11,11,11,12,12};

int colind[] = {0,0,1,1,2,1,3,1,4,2,4,5,3,4,6,7,8,4,9, 6, 9,10, 5, 9,11, 9,12};

double value[] = {1,-1,1.73205,-0.57735,1.29099,-0.57735,1.29099,-0.57735,1.91485,-0.774597,-0.522233,1.76841,

-0.774597,-0.522233,1.76841,1,1,-0.522233,1.93061,-0.56548,-0.51797,1.18825,-0.56548,-0.51797,1.18825,-0.51797,0.855399};

// Shift the index for 1-based indexing

for (int i=0;i

rowind += 1;

colind += 1;

}

// The Vectors

double input[] = { 0,0,0,0,0,2,0,2,0,0,0,0,0 };

double *x = new double();

double *y = new double();

mkl_dcootrsv(&uplo,&transa,&diag,&m,value,rowind,colind,&nz,input,x);

cout << "First result (vector x):\n";

for (int i=0;i

cout << x << " ";

cout << endl;

transa = 't';

mkl_dcootrsv(&uplo,&transa,&diag,&m,value,rowind,colind,&nz,x,y);

cout << "Solution of the equation (vector y):\n";

for (int i=0;i

cout << y << " ";

cout << endl;

return 0;

}

As the result for the equation I get the vector

0.350333 0.350333 0.470625 0.366454 0.21392 0.784373 0 2 0.121523 0.452948 0 0 0

which is definitely wrong, as the last non-zero element shows. The last three zeros are correct, as the rhs is zero in this points. The equation for the last non-zero element is1.93061*y[10]=0.538216, so y[10] should be 0.27878, not0.452948.

If I go one step back with0.452948, I get1.18825 as the matrix element, which is exactly the diagonal element in row 11 of the matrix. This is really strange. Are there still problems with indexing?

Sren

[Sören_Plönnigs]

Does nobody have an idea? This problem is really important to me, as I need to implement this for my bachelor thesis. I don't know where else I could get help on this lib...

Regards,

Sren

[mecej4 · Accepted solution]

In line 17 of your program, you should have nz instead of m as the loop count.

Until you gain confidence in the correctness of your code, it would help to print out the matrix in canonical form, for that would have shown you that 14 of the 27 non-zero entries in your matrix were misplaced.

[Ying_H_Intel]

Thank youMecej4.

Yes, I have run the fix with and without debugging. It can run fine. The right result should be
First result (vector x):
0 0 0 0 0 1.13096 0 2 0 0 0 0.538216 0
Solution of the equation (vector y):
0.23923 0.23923 0.470625 0 0.247063 0.784373 0 2 0 0.121523 0 0.452948 0
Press any key to continue . . .

Regards,
Ying

[Sören_Plönnigs]

Thank you for correcting my error! I now wrote a function for transforming my coordinate matrices to canonical form. I hope in future I won't ask such stupid questions anymore.

I still have some problems, but they have nothing to do with this thread. I will mark it as solved.