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

LAPACKE_dgetrf _dgetrf input error

marcoxp
Beginner
685 Views
I have a vectorA(n*n), initialized inside a loop like:

for(i=0;ifor(j=0;jfor(ysym=0;ysym<2;ysym){
for(zsym=0;zsym<2;zsym){

//computations

A[i*n+j]+=//values


I then solve a linear system with:

int info;
vectorb(n) //initialized to values
vectorpivots(n);

info=LAPACKE_dgetrf(LAPACK_ROW_MAJOR,n,n,&(A[0]),n,&pivots[0]);

info2=LAPACKE_dgetrs(LAPACK_ROW_MAJOR,'N',n,1,&A[0],n,&pivots[0],&b[0],1);

Results: info=-4;info2=-5; that means problems on A.

If I change the A initialization to:
A[i*n+j]=//values (plus is removed)

then the program runs, info=0, info2=0. But the solution b differs from the solution of the same system obtained in another solver.

How to get out of this? Thanks.
0 Kudos
5 Replies
mecej4
Honored Contributor III
685 Views
Show the actual code, not a paraphrase.

You have not shown any initialization of the elements of A. Do you perform the initialization somewhere else? Otherwise

A[i*n+j]+=

adds a defined quantity to an uninitialized array element, leaving the result undefined.

Why are you surprised when you solve two different problems with two different solvers and obtain different results?
0 Kudos
marcoxp
Beginner
685 Views
A is declared and initialized in the main function as: vectorA(n*n,0);

then passed by reference to:

void myFun(vector &A)

where the A[i*n+j]+= happens.

On exit of myFun all elements of A are finite.

All values of A and b are then also exported manually to matlab and sytem solved there too. The matrix is not ill conditioned so I don't see how the two problems are different.
0 Kudos
Alexander_K_Intel3
685 Views

Hello,

The INFO=-4 for LAPACKE_DGETRFand INFO=-5 for LAPACKE_DGETRSare set by NaN checker. This means that your initial matrix contains NaNs. This correlates with your replacing += with =, because NaN+=value is still NaN. And with = operator you put your correct values substituting all NaNs.

And how large is the difference between solutions you mentioned? If you compare your solutions with just equality (==) then this is expected, few last digits of floating point numbers could differ due to different order of simple operations in algorithms, this results other sequence of rounding operations on intermediate data. Correct way to compare two matrices is to check residual norm(A1 - A2) / ( N * norm(A1) * EPS ).

W.B.R.,
Alexander

0 Kudos
marcoxp
Beginner
685 Views
My apologies, indeed A contains NaNs. When importing into matlab to test it, those NaNs magically go. The differences are not few decimals, they're quite big (when using = instead of +=) but at this point I don't trust the whole import method, so I'll just fix A.
Thanks for pointing to the NaN check, do you have a link with the whole list of parameters? The reference I use says only: " -i means incorrect i-th input paramater".
0 Kudos
mecej4
Honored Contributor III
685 Views
NaN arithmetic is so special, I feel, that it makes no sense to say that when they were used in a calculation the results were found to differ "much" from some other results. You have to be able to certify that no NaNs were used before you can discuss any property of the results.

> do you have a link with the whole list of parameters?

When you call Lapack routines through Lapack-E or Lapack-95, the return code from the F77 Lapack routines tend to be of little use, unless you are willing to step through the interface code. It is best to check that the arguments in your call meet the specifications.
0 Kudos
Reply