I then solve a linear system with:
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.
You have not shown any initialization of the elements of A. Do you perform the initialization somewhere else? Otherwise
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?
then passed by reference to:
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.
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 ).
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".
> 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.