- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I have a vectorA(n*n), initialized inside a loop like:
for(i=0;i for(j=0;j for(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.
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
- Report Inappropriate Content
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?
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
- Report Inappropriate Content
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.
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
- 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
- Report Inappropriate Content
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".
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
- Report Inappropriate Content
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.
> 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.

Reply
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