[ADOL-C] sparse_jac return unreasonable NaN's
何旭
mailhexu at gmail.com
Sun Aug 28 04:27:15 EDT 2011
Hello,
I made some simple test on sparse_jac, and I found that if my function is
with condassign() and if any branch of condassign() returns a nan, the
sparse_jac returns a nan,too. I think it there might be a bug in condassign
or sparse_jac.
My function is very simple. y=x/(exp(x)-1) . However if x is too large
(x>100 or so) ,exp(x) results in nan, if fabs(x) is too small
(fabs(x)<1.0e-3), exp(x)-1 ->0. So I wrote some conditional statement to
avoid nan's. But that seems to be useless with sparse_jac. In my test,
x=2000, which means the function return 0 and the jacobian be 0. But that
does not happen. That's result of test1.
I tried to substitute y=x/(exp(x)-1) with y=x, which never return nan. And
the result is right. (test 3)
So I guess there might be some logical error in sparse_jac .
I put the test results here. test2 seems to be of no use. Just ignore it.
Result of test1:
x= 20000:
0.000000e+00
dense:0.000000e+00
number of sparse matrix elements: 1
sparse: -nan
Result of test2:
x= 20000:
0.000000e+00
dense:0.000000e+00
number of sparse matrix elements: 0
sparse: -3.848143e-42
Result of test3:
x= 20000:
0.000000e+00
dense:0.000000e+00
number of sparse matrix elements: 1
sparse: 0.000000e+00
And here's my testing code.
---------------------------------------------------------------------------------------
#include <iostream>
#include <adolc/adolc.h>
#include <adolc/adolc_sparse.h>
using namespace std;
adouble myfunc(adouble x)
{
adouble y;
condassign(y,1.0e-03-fabs(x),1.0);
condassign(y,x-100.0,0.0);
condassign(y,(1.0e-3-fabs(x))*(x-100.0),x/(exp(x)-1.0)); // TEST 1
// NOTHING HERE
//TEST 2
//condassign(y,(1.0e-3-fabs(x))*(x-100.0),x); //TEST
3
return y;
}
int main()
{
adouble ax,ay;
double x,y;
x=200.0;
trace_on(1);
ax<<=x;
ay=myfunc(ax);
ay>>=y;
trace_off();
double **jac;
jac=(double **)malloc(sizeof(double*));
jac[0]=(double *)malloc(sizeof(double));
x=20000.0;
cout<<"x= "<<x<<":"<<endl;
function(1,1,1,&x,&y);
cout<<scientific<<y<<endl;
// jacobian
jacobian(1,1,1,&x,jac);
cout<<"dense:"<<jac[0][0]<<endl;
free(jac[0]);
free(jac);
//sparse_jac
int nnz;
unsigned int *rindex=NULL;
unsigned int *cindex=NULL;
double *sparse_Jacobi_value=NULL;
int options[4]={0,0,0,1};
sparse_jac(1,1,1,0,&x,&nnz,&rindex,&cindex,&sparse_Jacobi_value,options);
cout<<"number of sparse matrix elements: "<<nnz<<'\n'<<"sparse:
"<<sparse_Jacobi_value[0]<<endl;
free(rindex);
free(cindex);
free(sparse_Jacobi_value);
return 0;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://list.coin-or.org/pipermail/adol-c/attachments/20110828/5fdbec71/attachment.html
More information about the ADOL-C
mailing list