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

for(i=0;i

for(zsym=0;zsym<2;zsym){

//computations

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

I then solve a linear system with:

int info;

vector

vector

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.

Link Copied

5 Replies

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

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?

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

then passed by reference to:

void myFun(vector

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.

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

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

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

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".

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

*> 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.

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