[Developers] bug in file fvar_m41.cpp
dave fournier
davef at otter-rsch.com
Fri Oct 2 12:46:05 PDT 2009
dave fournier wrote:
> at line 11
> replace
> x.elem_value(imin)=v.elem_value(1)/m.elem_value(imin,imin);
> by
> x.elem_value(imin)=v.elem_value(imin)/m.elem_value(imin,imin);
>
>
But the adjoint code needs to be fixed as well.
void dfbltsolve(void)
{
verify_identifier_string("ww");
dvar_vector_position xpos=restore_dvar_vector_position();
dvar_vector_position vpos=restore_dvar_vector_position();
dvector v=restore_dvar_vector_value(vpos);
dvar_matrix_position mpos=restore_dvar_matrix_position();
banded_lower_triangular_dmatrix m=
restore_banded_lower_triangular_dvar_matrix_value(mpos);
verify_identifier_string("rt");
dvector dfx=
restore_dvar_vector_derivatives(xpos);
int bw=m.bandwidth();
int imin=m.indexmin();
int imax=m.indexmax();
banded_lower_triangular_dmatrix dfm(imin,imax,bw);
dvector x(imin,imax);
dvector dfv(imin,imax);
dfm.initialize();
dfv.initialize();
x(imin)=v(imin)/m(imin,imin);
dvector dfsum(imin,imax);
dfsum.initialize();
dvector ssum(imin,imax);
ssum.initialize();
int i;
for (i=imin+1;i<=imax;i++)
{
int jmin=admax(imin,i-bw+1);
for (int j=jmin;j<=i-1;j++)
{
ssum(i)+=m(i,j)*x(j);
}
x(i)=(v(i)-ssum(i))/m(i,i);
}
for (i=imax;i>=imin+1;i--)
{
int jmin=admax(imin,i-bw+1);
//x(i)=(v(i)-ssum(i))/m(i,i);
dfv(i)+=dfx(i)/m(i,i);
dfsum(i)-=dfx(i)/m(i,i);
dfm(i,i)-=dfx(i)*(v(i)-ssum(i))/(m(i,i)*m(i,i));
dfx(i)=0.0;
for (int j=i-1;j>=jmin;j--)
{
//ssum(i)+=m(i,j)*x(j);
dfm(i,j)+=dfsum(i)*x(j);
dfx(j)+=dfsum(i)*m(i,j);
}
}
//x(imin)=v(1)/m(imin,imin);
// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
//dfm(imax,imax)=0.0;
// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
dfv(imin)+=dfx(imin)/m(imin,imin);
dfm(imin,imin)-=dfx(imin)*v(imin)/(m(imin,imin)*m(imin,imin));
dfx(imin)=0.0;
dfm.save_dmatrix_derivatives(mpos);
dfv.save_dvector_derivatives(vpos);
}
--
David A. Fournier
P.O. Box 2040,
Sidney, B.C. V8l 3S3
Canada
Phone/FAX 250-655-3364
http://otter-rsch.com
More information about the Developers
mailing list