Hello everybody,<br><br>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. <br><br>I tried ADOL-C 2.1.12 ,2.2.0 and 2.2.1 and all of them produced the same result.<br>
ColPack version : 1.0.3.<br>gcc 4.5.2<br><br>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.<br><br>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.<br>
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 .<br>Sparse:<br>col row value<br>0 0 1.000000e+00<br>
1 1 1.000000e+00<br>2 2 nan<br>3 0 3.644718e-07<br>3 3 -1.483424e-05<br>3 4 -8.319726e-62<br>3 5 nan<br>3 6 3.644718e-07<br>3 8 nan<br>4 1 -6.626760e-08<br>4 2 nan<br>
4 3 1.602189e-10<br>4 4 1.325352e-07<br>4 5 nan<br><br>Dense<br>1.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 <br>0.000000e+00 1.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 <br>
0.000000e+00 0.000000e+00 1.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 <br>3.644718e-07 0.000000e+00 0.000000e+00 -1.483424e-05 -8.319726e-62 5.376117e-32 <br>0.000000e+00 -6.626760e-08 -1.882395e+19 1.602189e-10 1.325352e-07 3.764790e+19 <br>
0.000000e+00 0.000000e+00 1.000000e+00 -4.522029e-27 4.522029e-27 -2.000000e+00<br><br>Thanks in advance.<br><br><br><br>And the codes related are put here.<br><br>adouble norm_abern(adouble x)
<br>{
<br> adouble nx,y;
<br> nx=x*qe/(kb*T);
<br> condassign(y,1.0e-5-fabs(nx),1.0,nx/(exp(nx)-1.0));
<br> condassign(y,nx-100.0,0.0);
<br> return y;<br>}<br><br>void block_1d::self_consistent_coupled_mktape()<br>{<br> //dependent vars<br> adouble *aFn=new adouble[N];<br> adouble *aFp=new adouble[N];<br> adouble *aFV=new adouble[N];<br><br> double *Fn=new double[N];<br>
double *Fp=new double[N];<br> double *FV=new double[N];<br> <br> // independent vars<br> adouble *an=new adouble[N];<br> adouble *ap=new adouble[N];<br> adouble *aV=new adouble[N];<br><br><br><br> trace_on(1);<br>
//dependent vars : an0, ap0, aV0, an1, ap1, ......<br> for(int i=0;i<N;i++)<br> {<br> an[i]<<=n[i];<br> ap[i]<<=p[i];<br> aV[i]<<=V[i];<br> }<br><br> adouble *aR=new adouble[N];<br> for(int i=0;i<N;i++)<br>
{<br> aR[i]=(an[i]*ap[i]-sqr(ni[i]))/(tp[i]*(an[i]+ni[i])+tn[i]*(ap[i]+ni[i]));<br> }<br> <br><br> for(int i=0;i<1;i++)<br> {<br> aFn[i]=an[i]-n[i];<br> aFp[i]=ap[i]-p[i];<br> aFV[i]=aV[i]-V[i];<br> }<br>
<br> for(int i=1;i<N-1;i++)<br> {<br> //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];<br>
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];<br><br><br> 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];<br>
<br> //condassign();<br> int k;<br> k=i-1;<br> <br> aFV[i]=(aV[i-1]+aV[i+1]-2.0*aV[i])-sqr(hm[i])*qe*(an[i]-ap[i]-dfile[i])/Epsilon[i];<br> }<br><br> for(int i=N-1;i<N;i++)<br> {<br> aFn[i]=an[i]-n[i];<br>
aFp[i]=ap[i]-p[i];<br> aFV[i]=aV[i]-V[i];<br> }<br><br><br> for(int i=0;i<N;i++)<br> {<br> aFn[i]>>=Fn[i];<br> aFp[i]>>=Fp[i];<br> aFV[i]>>=FV[i];<br> }<br> <br> trace_off();<br>
<br><br> delete []aFn;<br> delete []aFp;<br> delete []aFV;<br> delete []Fn;<br> delete []Fp;<br> delete []FV;<br> delete []an;<br> delete []ap;<br> delete []aV;<br>}<br><br><br>void block_1d::self_consistent_solve_Jacobi()<br>
{<br> double *delta_npV=new double[3*N];<br> double *array_npV=new double[3*N];<br> double *array_F=new double[3*N];<br><br> double **jac;<br> jac=(double **)malloc(3*N*sizeof(double*));<br> for (int i=0;i<3*N;i++)<br>
{<br> jac[i]=(double *)malloc(3*N*sizeof(double));<br> }<br> //cout<<"allocated"<<endl;<br> <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> for (int i=0;i<N;i++)<br> {<br> array_npV[3*i]=n[i];<br> array_npV[3*i+1]=p[i];<br> array_npV[3*i+2]=V[i];<br> }<br> <br> function(1,3*N,3*N,array_npV,array_F);<br>
for (int i=0;i<3*N;i++)<br> {<br> array_F[i]=-array_F[i];<br> }<br><br> <br> //array2file(array_F,N*3,"F.dat");<br> //array2file(array_F,3*N-6,"./array_f.dat");<br> <br> <br> jacobian(1,3*N,3*N,array_npV,jac);<br>
cout<<"Jac done"<<endl;<br> ofstream jaco("jac3.dat");<br> for (int i=0;i<6;i++)<br> {<br> for (int j=0;j<6;j++)<br> {<br> jaco<<scientific<<jac[i][j]<<'\t';<br>
}<br> jaco<<'\n';<br> }<br> jaco.close();<br> <br> <br> sparse_jac(1,3*N,3*N,0,array_npV,&nnz,&rindex,&cindex,&sparse_Jacobi_value,options);<br><br>......<br><br>}<br>
<br><br>Xu He<br><br>