ADMB Documentation  11.2.2828
 All Classes Files Functions Variables Typedefs Friends Defines
fvar_m41.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: fvar_m41.cpp 1919 2014-04-22 22:02:01Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include "fvar.hpp"
00012 void dfbltsolve(void);
00013 
00018 dvar_vector solve(const banded_lower_triangular_dvar_matrix& m,
00019   const dvar_vector &v)
00020 {
00021   int bw=m.bandwidth();
00022   int imin=m.indexmin();
00023   int imax=m.indexmax();
00024   dvar_vector x(imin,imax);
00025   x.elem_value(imin)=v.elem_value(imin)/m.elem_value(imin,imin);
00026   for (int i=imin+1;i<=imax;i++)
00027   {
00028     int jmin=admax(imin,i-bw+1);
00029     double ssum=0.0;
00030     for (int j=jmin;j<=i-1;j++)
00031     {
00032       ssum+=m.elem_value(i,j)*x.elem_value(j);
00033     }
00034     x.elem_value(i)=(v.elem_value(i)-ssum)/m.elem_value(i,i);
00035   }
00036   save_identifier_string("rt");
00037   m.save_dvar_matrix_value();
00038   m.save_dvar_matrix_position();
00039   v.save_dvar_vector_value();
00040   v.save_dvar_vector_position();
00041   x.save_dvar_vector_position();
00042   save_identifier_string("ww");
00043   gradient_structure::GRAD_STACK1->
00044       set_gradient_stack(dfbltsolve);
00045   return x;
00046 }
00047 
00048 /*
00049 dvar_vector solve(const banded_lower_triangular_dvar_matrix& m,
00050   const dvar_vector &v, dvariable& lndet)
00051 {
00052   int bw=m.bandwidth();
00053   int imin=m.indexmin();
00054   int imax=m.indexmax();
00055   dvar_vector x(imin,imax);
00056   value(x(imin))=value(v(1))/value(m(1,1));
00057   for (int i=2;i<=imax;i++)
00058   {
00059     int jmin=admax(1,i-bw+1);
00060     double ssum=0.0;
00061     for (int j=jmin;j<=i-1;j++)
00062     {
00063       ssum+=value(m(i,j))*value(x(j));
00064     }
00065     value(x(i))=(value(v(i))-ssum)/value(m(i,i));
00066   }
00067   save_identifier_string("rt");
00068   m.save_dvar_matrix_value();
00069   m.save_dvar_matrix_position();
00070   v.save_dvar_vector_value();
00071   v.save_dvar_vector_position();
00072   x.save_dvar_vector_position();
00073   save_identifier_string("ww");
00074   gradient_structure::GRAD_STACK1->
00075       set_gradient_stack(dfbltsolve);
00076   return x;
00077 }
00078 */
00079 
00084 void dfbltsolve(void)
00085 {
00086   verify_identifier_string("ww");
00087   dvar_vector_position xpos=restore_dvar_vector_position();
00088   dvar_vector_position vpos=restore_dvar_vector_position();
00089   dvector v=restore_dvar_vector_value(vpos);
00090   dvar_matrix_position mpos=restore_dvar_matrix_position();
00091   banded_lower_triangular_dmatrix m=
00092     restore_banded_lower_triangular_dvar_matrix_value(mpos);
00093   verify_identifier_string("rt");
00094   dvector dfx=
00095     restore_dvar_vector_derivatives(xpos);
00096 
00097   int bw=m.bandwidth();
00098   int imin=m.indexmin();
00099   int imax=m.indexmax();
00100 
00101   banded_lower_triangular_dmatrix dfm(imin,imax,bw);
00102   dvector x(imin,imax);
00103   dvector dfv(imin,imax);
00104   dfm.initialize();
00105   dfv.initialize();
00106   x(imin)=v(1)/m(1,1);
00107   dvector dfsum(imin,imax);
00108   dfsum.initialize();
00109   dvector ssum(imin,imax);
00110   ssum.initialize();
00111   int i;
00112   for (i=2;i<=imax;i++)
00113   {
00114     int jmin=admax(1,i-bw+1);
00115     for (int j=jmin;j<=i-1;j++)
00116     {
00117       ssum(i)+=m(i,j)*x(j);
00118     }
00119     x(i)=(v(i)-ssum(i))/m(i,i);
00120   }
00121 
00122   for (i=imax;i>=2;i--)
00123   {
00124     int jmin=admax(1,i-bw+1);
00125     //x(i)=(v(i)-ssum(i))/m(i,i);
00126     dfv(i)+=dfx(i)/m(i,i);
00127     dfsum(i)-=dfx(i)/m(i,i);
00128     dfm(i,i)-=dfx(i)*(v(i)-ssum(i))/(m(i,i)*m(i,i));
00129     dfx(i)=0.0;
00130     for (int j=i-1;j>=jmin;j--)
00131     {
00132       //ssum(i)+=m(i,j)*x(j);
00133       dfm(i,j)+=dfsum(i)*x(j);
00134       dfx(j)+=dfsum(i)*m(i,j);
00135     }
00136   }
00137   //x(imin)=v(1)/m(1,1);
00138   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00139   //dfm(imax,imax)=0.0;
00140   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00141   //x.elem_value(imin)=v.elem_value(imin)/m.elem_value(imin,imin);
00142   dfv(imin)+=dfx(imin)/m(imin,imin);
00143   dfm(imin,imin)-=dfx(imin)*v(imin)/(m(imin,imin)*m(imin,imin));
00144   dfx(imin)=0.0;
00145 
00146   dfm.save_dmatrix_derivatives(mpos);
00147   dfv.save_dvector_derivatives(vpos);
00148 }