ADMB Documentation  11.1.1927
 All Classes Files Functions Variables Typedefs Friends Defines
fvar_m50.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: fvar_m50.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 
00013 
00014 void dfcholeski_decomp_banded_positive(void);
00015 
00020 banded_lower_triangular_dvar_matrix choleski_decomp_positive(
00021   const banded_symmetric_dvar_matrix& MM, double eps,
00022   dvariable& _fpen)
00023 {
00024   // kludge to deal with constantness
00025   banded_symmetric_dvar_matrix& M = (banded_symmetric_dvar_matrix&) MM;
00026   int n=M.indexmax();
00027 
00028   int bw=M.bandwidth();
00029   banded_lower_triangular_dvar_matrix L(1,n,bw);
00030 #ifndef SAFE_INITIALIZE
00031     L.initialize();
00032 #endif
00033 
00034   int i,j,k;
00035   double tmp;
00036   double ptmp;
00037   double fpen=0.0;
00038   ptmp=posfun(value(M(1,1)),eps,fpen);
00039   L.elem_value(1,1)=sqrt(ptmp);
00040   for (i=2;i<=bw;i++)
00041   {
00042     L.elem_value(i,1)=value(M(i,1))/L.elem_value(1,1);
00043   }
00044 
00045   for (i=2;i<=n;i++)
00046   {
00047     int jmin=admax(2,i-bw+1);
00048     for (j=jmin;j<=i-1;j++)
00049     {
00050       tmp=value(M(i,j));
00051       int kmin=max(1,j-bw+1,i-bw+1);
00052       for (k=kmin;k<=j-1;k++)
00053       {
00054         tmp-=L.elem_value(i,k)*L.elem_value(j,k);
00055       }
00056       L.elem_value(i,j)=tmp/L.elem_value(j,j);
00057     }
00058     tmp=value(M(i,i));
00059     int kmin=admax(i-bw+1,1);
00060     for (k=kmin;k<=i-1;k++)
00061     {
00062       tmp-=L.elem_value(i,k)*L.elem_value(i,k);
00063     }
00064     ptmp=posfun(tmp,eps,fpen);
00065     L.elem_value(i,i)=sqrt(ptmp);
00066   }
00067 
00068 
00069   if (fpen>0)
00070     cout << "fpen = " << fpen << endl;
00071 
00072    value(_fpen)=fpen;
00073   //banded_lower_triangular_dvar_matrix vc=nograd_assign(L);
00074   save_identifier_string("qs");
00075   _fpen.save_prevariable_position();
00076   save_double_value(eps);
00077   save_identifier_string("rs");
00078   L.save_dvar_matrix_position();
00079   save_identifier_string("rt");
00080   MM.save_dvar_matrix_value();
00081   save_identifier_string("rl");
00082   MM.save_dvar_matrix_position();
00083   save_identifier_string("ro");
00084   gradient_structure::GRAD_STACK1->
00085       set_gradient_stack(dfcholeski_decomp_banded_positive);
00086 
00087   return L;
00088 }
00089 
00094 void dfcholeski_decomp_banded_positive(void)
00095 {
00096   verify_identifier_string("ro");
00097   dvar_matrix_position MMpos=restore_dvar_matrix_position();
00098   verify_identifier_string("rl");
00099   banded_symmetric_dmatrix M=
00100     restore_banded_symmetric_dvar_matrix_value(MMpos);
00101   verify_identifier_string("rt");
00102   dvar_matrix_position vcpos=restore_dvar_matrix_position();
00103   verify_identifier_string("rs");
00104   banded_lower_triangular_dmatrix dfL=
00105     restore_banded_lower_triangular_dvar_matrix_derivatives(vcpos);
00106   double eps=restore_double_value();
00107   prevariable_position fpenpos=restore_prevariable_position();
00108   verify_identifier_string("qs");
00109   double dfpen=restore_prevariable_derivative(fpenpos);
00110 
00111   int n=M.indexmax();
00112   int bw=M.bandwidth();
00113 
00114   double fpen=0.0;
00115   banded_lower_triangular_dmatrix L(1,n,bw);
00116   banded_lower_triangular_dmatrix tmp1(1,n,bw);
00117   banded_lower_triangular_dmatrix dftmp1(1,n,bw);
00118   banded_symmetric_dmatrix dfM(1,n,bw);
00119   dvector tmp(1,n);
00120   dvector ptmp(1,n);
00121   dvector dftmp(1,n);
00122   dvector dfptmp(1,n);
00123   tmp.initialize();
00124   ptmp.initialize();
00125   tmp1.initialize();
00126   dftmp.initialize();
00127   dfptmp.initialize();
00128   dftmp1.initialize();
00129   dfM.initialize();
00130 #ifndef SAFE_INITIALIZE
00131     L.initialize();
00132 #endif
00133 
00134   int i,j,k;
00135   ptmp(1)=posfun(M(1,1),eps,fpen);
00136   L(1,1)=sqrt(ptmp(1));
00137   L(1,1)=sqrt(M(1,1));
00138   for (i=2;i<=bw;i++)
00139   {
00140     L(i,1)=M(i,1)/L(1,1);
00141   }
00142 
00143   for (i=2;i<=n;i++)
00144   {
00145     int jmin=admax(2,i-bw+1);
00146     for (j=jmin;j<=i-1;j++)
00147     {
00148       tmp1(i,j)=M(i,j);
00149       int kmin=max(1,j-bw+1,i-bw+1);
00150       for (k=kmin;k<=j-1;k++)
00151       {
00152         tmp1(i,j)-=L(i,k)*L(j,k);
00153       }
00154       L(i,j)=tmp1(i,j)/L(j,j);
00155     }
00156     tmp(i)=M(i,i);
00157     int kmin=admax(i-bw+1,1);
00158     for (k=kmin;k<=i-1;k++)
00159     {
00160       tmp(i)-=L(i,k)*L(i,k);
00161     }
00162     double pen;
00163     ptmp(i)=posfun(tmp(i),eps,pen);
00164     L(i,i)=sqrt(ptmp(i));
00165   }
00166 
00167   dfptmp.initialize();
00168 
00169   for (i=n;i>=2;i--)
00170   {
00171     //L(i,i)=sqrt(ptmp(i));
00172     dfptmp(i)+=dfL(i,i)/(2.0*L(i,i));
00173     dfL(i,i)=0.0;
00174     // ptmp(i)=posfun(tmp(i),eps);
00175     dftmp(i)=dfptmp(i)*dfposfun(tmp(i),eps);
00176     dftmp(i)+=dfpen*dfposfun1(tmp(i),eps);
00177     dfptmp(i)=0.0;
00178 
00179     int kmin=admax(i-bw+1,1);
00180     for (k=i-1;k>=kmin;k--)
00181     {
00182       //tmp(i)-=L(i,k)*L(i,k);
00183       dfL(i,k)-=2.*dftmp(i)*L(i,k);
00184     }
00185     //tmp(i)=M(i,i);
00186     dfM(i,i)+=dftmp(i);
00187     dftmp(i)=0.0;
00188     int jmin=admax(2,i-bw+1);
00189     for (j=i-1;j>=jmin;j--)
00190     {
00191       //L(i,j)=tmp1(i,j)/L(j,j);
00192       double linv=1./L(j,j);
00193       dftmp1(i,j)+=dfL(i,j)*linv;
00194       dfL(j,j)-=dfL(i,j)*tmp1(i,j)*linv*linv;
00195       dfL(i,j)=0.0;
00196       int kmin=max(1,j-bw+1,i-bw+1);
00197       for (k=j-1;k>=kmin;k--)
00198       {
00199         //tmp(i,j)-=L(i,k)*L(j,k);
00200         dfL(i,k)-=dftmp1(i,j)*L(j,k);
00201         dfL(j,k)-=dftmp1(i,j)*L(i,k);
00202       }
00203       //tmp(i,j)=M(i,j);
00204       dfM(i,j)+=dftmp1(i,j);
00205       dftmp1(i,j)=0.0;
00206     }
00207   }
00208   double linv=1./L(1,1);
00209   for (i=bw;i>=2;i--)
00210   {
00211     //L(i,1)=M(i,1)/L(1,1);
00212     dfM(i,1)+=dfL(i,1)*linv;
00213     dfL(1,1)-=dfL(i,1)*M(i,1)*linv*linv;
00214     dfL(i,1)=0.0;
00215   }
00216   //L(1,1)=sqrt(M(1,1));
00217   //dfM(1,1)+=dfL(1,1)/(2.*L(1,1));
00218 
00219   //L(1,1)=sqrt(ptmp(1));
00220   dfptmp(1)+=dfL(1,1)/(2.*L(1,1));
00221   dfL(1,1)=0.0;
00222   //ptmp(1)=posfun(M(1,1),eps,fpen);
00223   dfM(1,1)+=dfptmp(1)*dfposfun(M(1,1),eps);
00224   dfM(1,1)+=dfpen*dfposfun1(M(1,1),eps);;
00225   dfptmp(1)=0.0;
00226 
00227   dfM.save_dmatrix_derivatives(MMpos);
00228 }