[ADOL-C] sparse_jac return unreasonable NaN's

Fri Aug 26 10:56:59 EDT 2011

Hello  everybody,

I used sparse_jac to generate a Jacobian and got a lot of NaN elements which
seems to be incorrect. So I tried jacobian and I found  the dense matrix
correct.

I tried ADOL-C 2.1.12 ,2.2.0 and 2.2.1 and all of them produced the same
result.
ColPack version : 1.0.3.
gcc 4.5.2

I'd like to ask  the reason why the two function return different results
and how I can make the sparse matrix same with the dense one.

As my code is too complicated to be put here ,I'll try to make a neat
version which may have similar probelem if in need. I put part of the two
matrix here.
It is notable that the third element (col 2, row 2) in the sparse matrix
should be 1.0 apparently, because  y_3=x_3- some_const. And I can not find
any other operation on y_3 .
Sparse:
col  row value
0    0    1.000000e+00
1    1    1.000000e+00
2    2    nan
3    0    3.644718e-07
3    3    -1.483424e-05
3    4    -8.319726e-62
3    5    nan
3    6    3.644718e-07
3    8    nan
4    1    -6.626760e-08
4    2    nan
4    3    1.602189e-10
4    4    1.325352e-07
4    5    nan

Dense
1.000000e+00    0.000000e+00    0.000000e+00    0.000000e+00
0.000000e+00    0.000000e+00
0.000000e+00    1.000000e+00    0.000000e+00    0.000000e+00
0.000000e+00    0.000000e+00
0.000000e+00    0.000000e+00    1.000000e+00    0.000000e+00
0.000000e+00    0.000000e+00
3.644718e-07    0.000000e+00    0.000000e+00    -1.483424e-05
-8.319726e-62    5.376117e-32
0.000000e+00    -6.626760e-08    -1.882395e+19    1.602189e-10
1.325352e-07    3.764790e+19
0.000000e+00    0.000000e+00    1.000000e+00    -4.522029e-27
4.522029e-27    -2.000000e+00

Thanks in advance.

And the codes related are put here.

adouble norm_abern(adouble x)
{
  adouble nx,y;
  nx=x*qe/(kb*T);
  condassign(y,1.0e-5-fabs(nx),1.0,nx/(exp(nx)-1.0));
  condassign(y,nx-100.0,0.0);
  return y;
}

void block_1d::self_consistent_coupled_mktape()
{
  //dependent vars
  adouble *aFn=new adouble[N];
  adouble *aFp=new adouble[N];
  adouble *aFV=new adouble[N];

  double *Fn=new double[N];
  double *Fp=new double[N];
  double *FV=new double[N];

  // independent vars
  adouble *an=new adouble[N];
  adouble *ap=new adouble[N];
  adouble *aV=new adouble[N];

    trace_on(1);
//dependent vars :  an0, ap0, aV0, an1, ap1, ......
  for(int i=0;i<N;i++)
  {
    an[i]<<=n[i];
    ap[i]<<=p[i];
    aV[i]<<=V[i];
  }

  adouble *aR=new adouble[N];
  for(int i=0;i<N;i++)
  {

aR[i]=(an[i]*ap[i]-sqr(ni[i]))/(tp[i]*(an[i]+ni[i])+tn[i]*(ap[i]+ni[i]));
  }

  for(int i=0;i<1;i++)
  {
    aFn[i]=an[i]-n[i];
    aFp[i]=ap[i]-p[i];
    aFV[i]=aV[i]-V[i];
  }

  for(int i=1;i<N-1;i++)
  {

//aFn[i]=un[i]*kb*T*(an[i-1]-2.0*an[i]+an[i+1])/sqr(hm[i])-un[i]*qe*(an[i+1]-an[i-1])*(aV[i+1]-aV[i-1])/sqr(2.0*hm[i])-un[i]*qe*an[i]*(aV[i-1]-2.0*aV[i]+aV[i+1])/sqr(hm[i])-qe*aR[i];

aFp[i]=-up[i]*kb*T*(ap[i-1]-2.0*ap[i]+ap[i+1])/sqr(hm[i])-up[i]*qe*(ap[i+1]-ap[i-1])*(aV[i+1]-aV[i-1])/sqr(2.0*hm[i])-up[i]*qe*ap[i]*(aV[i-1]-2.0*aV[i]+aV[i+1])/sqr(hm[i])+qe*aR[i];

    aFn[i]=(
un[i+1]*kb*T*(an[i+1]-an[i])/hm[i]*norm_abern(V[i]-V[i+1])+un[i+1]*qe*an[i+1]*(aV[i]-aV[i+1])/hm[i]
-  un[i]*kb*T*(an[i]-an[i-1])/hm[i]*norm_abern(aV[i-1]-aV[i])
-un[i]*qe*an[i]*(norm_abern(aV[i-1]-aV[i]))/hm[i]  )/hm[i]-qe*aR[i];

    //condassign();
    int k;
    k=i-1;

aFV[i]=(aV[i-1]+aV[i+1]-2.0*aV[i])-sqr(hm[i])*qe*(an[i]-ap[i]-dfile[i])/Epsilon[i];
  }

  for(int i=N-1;i<N;i++)
  {
    aFn[i]=an[i]-n[i];
    aFp[i]=ap[i]-p[i];
    aFV[i]=aV[i]-V[i];
  }

  for(int i=0;i<N;i++)
  {
    aFn[i]>>=Fn[i];
    aFp[i]>>=Fp[i];
    aFV[i]>>=FV[i];
  }

  trace_off();

  delete []aFn;
  delete []aFp;
  delete []aFV;
  delete []Fn;
  delete []Fp;
  delete []FV;
  delete []an;
  delete []ap;
  delete []aV;
}

void block_1d::self_consistent_solve_Jacobi()
{
  double *delta_npV=new double[3*N];
  double *array_npV=new double[3*N];
  double *array_F=new double[3*N];

  double **jac;
  jac=(double **)malloc(3*N*sizeof(double*));
  for (int i=0;i<3*N;i++)
  {
  jac[i]=(double *)malloc(3*N*sizeof(double));
  }
  //cout<<"allocated"<<endl;

  int nnz;
  unsigned int *rindex=NULL;
  unsigned int *cindex=NULL;
  double *sparse_Jacobi_value=NULL;

  int options[4]={0,0,0,1};
  for (int i=0;i<N;i++)
  {
    array_npV[3*i]=n[i];
    array_npV[3*i+1]=p[i];
    array_npV[3*i+2]=V[i];
  }

  function(1,3*N,3*N,array_npV,array_F);
  for (int i=0;i<3*N;i++)
  {
    array_F[i]=-array_F[i];
  }

  //array2file(array_F,N*3,"F.dat");
  //array2file(array_F,3*N-6,"./array_f.dat");

    jacobian(1,3*N,3*N,array_npV,jac);
    cout<<"Jac done"<<endl;
    ofstream jaco("jac3.dat");
    for (int i=0;i<6;i++)
    {
      for (int j=0;j<6;j++)
      {
    jaco<<scientific<<jac[i][j]<<'\t';
      }
      jaco<<'\n';
    }
    jaco.close();

sparse_jac(1,3*N,3*N,0,array_npV,&nnz,&rindex,&cindex,&sparse_Jacobi_value,options);

......

}

Xu He
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://list.coin-or.org/pipermail/adol-c/attachments/20110826/71fa02c5/attachment.html 