[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