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

Problems with mkl_dcootrsv

Sören_Plönnigs
Beginner
527 Views
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
0 Kudos
1 Solution
mecej4
Honored Contributor III
527 Views
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.

View solution in original post

0 Kudos
9 Replies
Ying_H_Intel
Employee
527 Views
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
0 Kudos
barragan_villanueva_
Valued Contributor I
527 Views
Hi,

I can see that you got approximatly the same result. The most difference is about 0.006.
Do you expect more accuracy?
0 Kudos
Sören_Plönnigs
Beginner
527 Views
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
0 Kudos
Sören_Plönnigs
Beginner
527 Views
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
0 Kudos
Sören_Plönnigs
Beginner
527 Views
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
0 Kudos
Sören_Plönnigs
Beginner
527 Views
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
0 Kudos
mecej4
Honored Contributor III
528 Views
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.
0 Kudos
Ying_H_Intel
Employee
527 Views
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

0 Kudos
Sören_Plönnigs
Beginner
527 Views
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.
0 Kudos
Reply