Hello,<br><br>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.<br>
<br>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.<br>
<br>I tried to substitute y=x/(exp(x)-1) with y=x, which never return nan. And the result is right. (test 3)<br><br>So I guess there might be some logical error in sparse_jac . <br><br>I put the test results here. test2 seems to be of no use. Just ignore it.<br>
<br>Result of test1:<br>x= 20000:<br>0.000000e+00<br>dense:0.000000e+00<br>number of sparse matrix elements: 1<br>sparse: -nan<br><br>Result of test2:<br>x= 20000:<br>0.000000e+00<br>dense:0.000000e+00<br>number of sparse matrix elements: 0<br>
sparse: -3.848143e-42<br><br>Result of test3:<br>x= 20000:<br>0.000000e+00<br>dense:0.000000e+00<br>number of sparse matrix elements: 1<br>sparse: 0.000000e+00<br><br><br><br>And here's my testing code.<br>---------------------------------------------------------------------------------------<br>
#include <iostream><br>#include <adolc/adolc.h><br>#include <adolc/adolc_sparse.h><br>using namespace std;<br><br>adouble myfunc(adouble x)<br>{<br> adouble y;<br> condassign(y,1.0e-03-fabs(x),1.0);<br> condassign(y,x-100.0,0.0);<br>
condassign(y,(1.0e-3-fabs(x))*(x-100.0),x/(exp(x)-1.0)); // TEST 1<br>// NOTHING HERE //TEST 2 <br>//condassign(y,(1.0e-3-fabs(x))*(x-100.0),x); //TEST 3<br>
return y;<br>}<br><br><br>int main()<br>{<br> adouble ax,ay;<br> double x,y;<br> x=200.0;<br> trace_on(1);<br> ax<<=x;<br> ay=myfunc(ax);<br> ay>>=y;<br> trace_off();<br><br><br> double **jac;<br> jac=(double **)malloc(sizeof(double*));<br>
jac[0]=(double *)malloc(sizeof(double));<br><br> x=20000.0;<br> cout<<"x= "<<x<<":"<<endl;<br> function(1,1,1,&x,&y);<br> cout<<scientific<<y<<endl;<br>
<br><br>// jacobian<br> jacobian(1,1,1,&x,jac);<br> cout<<"dense:"<<jac[0][0]<<endl;<br> <br> free(jac[0]);<br> free(jac);<br> <br><br>//sparse_jac <br> int nnz;<br> unsigned int *rindex=NULL;<br>
unsigned int *cindex=NULL;<br> double *sparse_Jacobi_value=NULL;<br> <br> int options[4]={0,0,0,1}; <br> sparse_jac(1,1,1,0,&x,&nnz,&rindex,&cindex,&sparse_Jacobi_value,options);<br><br> cout<<"number of sparse matrix elements: "<<nnz<<'\n'<<"sparse: "<<sparse_Jacobi_value[0]<<endl;<br>
<br> free(rindex);<br> free(cindex);<br> free(sparse_Jacobi_value);<br> return 0;<br> <br>}<br><br>