ADMB Documentation  11.2.2853
 All Classes Files Functions Variables Typedefs Friends Defines
fvar_m39.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: fvar_m39.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 #ifdef __TURBOC__
00014 #pragma hdrstop
00015 #include <iostream.h>
00016 #endif
00017 
00018 #ifdef __ZTC__
00019 #include <iostream.hpp>
00020 #endif
00021 
00022 #ifdef __SUN__
00023 #include <iostream.h>
00024 #endif
00025 #ifdef __NDPX__
00026 #include <iostream.h>
00027 #endif
00028 void dfcholeski_decomp(void);
00029 
00034 dvar_matrix choleski_decomp(const dvar_matrix& MM)
00035 {
00036   // kludge to deal with constantness
00037   if (MM.colsize() != MM.rowsize())
00038   {
00039     cerr << "Error in chol_decomp. Matrix not square" << endl;
00040     ad_exit(1);
00041   }
00042   if (MM.colsize()==1)
00043   {
00044     int mmin=MM.indexmin();
00045     int mmax=MM.indexmax();
00046     if (MM(mmin,mmin)<=0)
00047     {
00048       cerr << "Error matrix not positive definite in choleski_decomp"
00049         <<endl;
00050       ad_exit(1);
00051     }
00052     dvar_matrix  vc(mmin,mmax,mmin,mmax);
00053     vc(mmin,mmin)=sqrt(MM(mmin,mmin));
00054     return vc;
00055   }
00056   dmatrix M=value(MM);
00057   int rowsave=M.rowmin();
00058   int colsave=M.colmin();
00059   M.rowshift(1);
00060   M.colshift(1);
00061   int n=M.rowmax();
00062 
00063   dmatrix L(1,n,1,n);
00064   //dmatrix C(1,n,1,n);
00065   //imatrix B(1,n,1,n);
00066   //B.initialize();
00067 #ifndef SAFE_INITIALIZE
00068     L.initialize();
00069 #endif
00070 
00071   int i,j,k;
00072   double tmp;
00073     if (M(1,1)<=0)
00074     {
00075       cerr << "Error matrix not positive definite in choleski_decomp"
00076         <<endl;
00077       ad_exit(1);
00078     }
00079   L(1,1)=sqrt(M(1,1));
00080   for (i=2;i<=n;i++)
00081   {
00082     /*
00083     if (!B(1,1))
00084     {
00085       C(1,1)=L(1,1);
00086       B(1,1)=1;
00087     }
00088     */
00089     L(i,1)=M(i,1)/L(1,1);
00090   }
00091 
00092   for (i=2;i<=n;i++)
00093   {
00094     for (j=2;j<=i-1;j++)
00095     {
00096       tmp=M(i,j);
00097       for (k=1;k<=j-1;k++)
00098       {
00099   /*
00100   if (!B(i,k))
00101   {
00102     C(i,k)=L(i,k);
00103     B(i,k)=1;
00104         }
00105   if (!B(j,k))
00106   {
00107     C(j,k)=L(j,k);
00108     B(j,k)=1;
00109         }
00110   */
00111         tmp-=L(i,k)*L(j,k);
00112       }
00113       /*
00114       if (!B(j,j))
00115       {
00116         C(j,j)=L(j,j);
00117         B(j,j)=1;
00118       }
00119       */
00120       L(i,j)=tmp/L(j,j);
00121     }
00122     tmp=M(i,i);
00123     for (k=1;k<=i-1;k++)
00124     {
00125       /*
00126       if (!B(i,k))
00127       {
00128         C(i,k)=L(i,k);
00129         B(i,k)=1;
00130       }
00131       */
00132       tmp-=L(i,k)*L(i,k);
00133     }
00134     if (tmp<=0)
00135     {
00136       cerr << "Error matrix not positive definite in choleski_decomp"
00137         <<endl;
00138       ad_exit(1);
00139     }
00140     L(i,i)=sqrt(tmp);
00141   }
00142   //C(n,n)=L(n,n);
00143   //B(n,n)=1;
00144   L.rowshift(rowsave);
00145   L.colshift(colsave);
00146   //cout << norm2(C-L) << endl;
00147   //cout << B << endl;
00148   //cout << L << endl;
00149 
00150   dvar_matrix vc=nograd_assign(L);
00151   save_identifier_string("rs");
00152   vc.save_dvar_matrix_position();
00153   save_identifier_string("rt");
00154   MM.save_dvar_matrix_value();
00155   save_identifier_string("rl");
00156   MM.save_dvar_matrix_position();
00157   save_identifier_string("ro");
00158   gradient_structure::GRAD_STACK1->
00159       set_gradient_stack(dfcholeski_decomp);
00160   return vc;
00161 }
00162 
00167 void dfcholeski_decomp(void)
00168 {
00169   verify_identifier_string("ro");
00170   dvar_matrix_position MMpos=restore_dvar_matrix_position();
00171   verify_identifier_string("rl");
00172   dmatrix M=restore_dvar_matrix_value(MMpos);
00173   verify_identifier_string("rt");
00174   dvar_matrix_position vcpos=restore_dvar_matrix_position();
00175   verify_identifier_string("rs");
00176   dmatrix dfL=restore_dvar_matrix_derivatives(vcpos);
00177 
00178   if (M.colsize() != M.rowsize())
00179   {
00180     cerr << "Error in chol_decomp. Matrix not square" << endl;
00181     ad_exit(1);
00182   }
00183   int rowsave=M.rowmin();
00184   int colsave=M.colmin();
00185   dfL.rowshift(1);
00186   dfL.colshift(1);
00187   M.rowshift(1);
00188   M.colshift(1);
00189   int n=M.rowmax();
00190 
00191   dmatrix L(1,n,1,n);
00192   dvector tmp(1,n);
00193   dmatrix tmp1(1,n,1,n);
00194   dmatrix dftmp1(1,n,1,n);
00195   dmatrix dfM(1,n,1,n);
00196   dvector dftmp(1,n);
00197   tmp.initialize();
00198   tmp1.initialize();
00199   dftmp.initialize();
00200   dftmp1.initialize();
00201   dfM.initialize();
00202 #ifndef SAFE_INITIALIZE
00203     L.initialize();
00204 #endif
00205 
00206   int i,j,k;
00207   if (M(1,1)<=0)
00208   {
00209     cerr << "Error matrix not positive definite in choleski_decomp"
00210       <<endl;
00211     ad_exit(1);
00212   }
00213   L(1,1)=sqrt(M(1,1));
00214   for (i=2;i<=n;i++)
00215   {
00216     L(i,1)=M(i,1)/L(1,1);
00217   }
00218 
00219   for (i=2;i<=n;i++)
00220   {
00221     for (j=2;j<=i-1;j++)
00222     {
00223       tmp1(i,j)=M(i,j);
00224       for (k=1;k<=j-1;k++)
00225       {
00226         tmp1(i,j)-=L(i,k)*L(j,k);
00227       }
00228       L(i,j)=tmp1(i,j)/L(j,j);
00229     }
00230     tmp(i)=M(i,i);
00231     for (k=1;k<=i-1;k++)
00232     {
00233       tmp(i)-=L(i,k)*L(i,k);
00234     }
00235     if (tmp(i)<=0)
00236     {
00237       cerr << "Error matrix not positive definite in choleski_decomp"
00238         <<endl;
00239       ad_exit(1);
00240     }
00241     L(i,i)=sqrt(tmp(i));
00242   }
00243  //*******************************************************************8
00244  //*******************************************************************8
00245  //*******************************************************************8
00246   for (i=n;i>=2;i--)
00247   {
00248     //L(i,i)=sqrt(tmp(i));
00249     dftmp(i)+=dfL(i,i)/(2.0*L(i,i));
00250     dfL(i,i)=0.0;
00251     for (k=i-1;k>=1;k--)
00252     {
00253       //tmp(i)-=L(i,k)*L(i,k);
00254       dfL(i,k)-=2.*dftmp(i)*L(i,k);
00255     }
00256     //tmp(i)=M(i,i);
00257     dfM(i,i)+=dftmp(i);
00258     dftmp(i)=0.0;
00259     for (j=i-1;j>=2;j--)
00260     {
00261       //L(i,j)=tmp1(i,j)/L(j,j);
00262       double linv=1./L(j,j);
00263       dftmp1(i,j)+=dfL(i,j)*linv;
00264       dfL(j,j)-=dfL(i,j)*tmp1(i,j)*linv*linv;
00265       dfL(i,j)=0.0;
00266       for (k=j-1;k>=1;k--)
00267       {
00268         //tmp(i,j)-=L(i,k)*L(j,k);
00269         dfL(i,k)-=dftmp1(i,j)*L(j,k);
00270         dfL(j,k)-=dftmp1(i,j)*L(i,k);
00271       }
00272       //tmp(i,j)=M(i,j);
00273       dfM(i,j)+=dftmp1(i,j);
00274       dftmp1(i,j)=0.0;
00275     }
00276   }
00277   double linv=1./L(1,1);
00278   for (i=n;i>=2;i--)
00279   {
00280     //L(i,1)=M(i,1)/L(1,1);
00281     dfM(i,1)+=dfL(i,1)*linv;
00282     dfL(1,1)-=dfL(i,1)*M(i,1)*linv*linv;
00283     dfL(i,1)=0.0;
00284   }
00285   //L(1,1)=sqrt(M(1,1));
00286   dfM(1,1)+=dfL(1,1)/(2.*L(1,1));
00287 
00288 
00289  //*******************************************************************8
00290  //*******************************************************************8
00291  //*******************************************************************8
00292   dfM.rowshift(rowsave);
00293   dfM.colshift(colsave);
00294 
00295   dfM.save_dmatrix_derivatives(MMpos);
00296 }