ADMB Documentation  11.1x.2738
 All Classes Files Functions Variables Typedefs Friends Defines
dmat35.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: dmat35.cpp 1112 2013-07-12 21:41:41Z 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 
00033 dmatrix choleski_decomp_positive(const dmatrix& MM,double bound)
00034 {
00035   // kludge to deal with constantness
00036   dmatrix & M= * (dmatrix *) &MM;
00037   if (M.colsize() != M.rowsize())
00038   {
00039     cerr << "Error in chol_decomp. Matrix not square" << endl;
00040     ad_exit(1);
00041   }
00042   int rowsave=M.rowmin();
00043   int colsave=M.colmin();
00044   M.rowshift(1);
00045   M.colshift(1);
00046   int n=M.rowmax();
00047 
00048   dmatrix L(1,n,1,n);
00049 #ifndef SAFE_INITIALIZE
00050     L.initialize();
00051 #endif
00052 
00053   int i,j,k;
00054   double tmp;
00055     if (M(1,1)<=bound)
00056     {
00057       M(1,1)=bound;
00058     }
00059   L(1,1)=sqrt(M(1,1));
00060   for (i=2;i<=n;i++)
00061   {
00062     L(i,1)=M(i,1)/L(1,1);
00063   }
00064 
00065   for (i=2;i<=n;i++)
00066   {
00067     for (j=2;j<=i-1;j++)
00068     {
00069       tmp=M(i,j);
00070       for (k=1;k<=j-1;k++)
00071       {
00072         tmp-=L(i,k)*L(j,k);
00073       }
00074       L(i,j)=tmp/L(j,j);
00075     }
00076     tmp=M(i,i);
00077     for (k=1;k<=i-1;k++)
00078     {
00079       tmp-=L(i,k)*L(i,k);
00080       if (tmp<=bound) break;
00081     }
00082     if (tmp<=bound)
00083     {
00084       tmp=bound;
00085     }
00086     L(i,i)=sqrt(tmp);
00087   }
00088   L.rowshift(rowsave);
00089   L.colshift(colsave);
00090 
00091   return L;
00092 }
00093 
00098 dmatrix choleski_decomp_positive(const dmatrix& MM,const int& _ierr)
00099 {
00100   // kludge to deal with constantness
00101   int & ierr = (int &)(_ierr);
00102   ierr=0;
00103   dmatrix & M= * (dmatrix *) &MM;
00104   if (M.colsize() != M.rowsize())
00105   {
00106     cerr << "Error in chol_decomp. Matrix not square" << endl;
00107     ad_exit(1);
00108   }
00109   int rowsave=M.rowmin();
00110   int colsave=M.colmin();
00111   M.rowshift(1);
00112   M.colshift(1);
00113   int n=M.rowmax();
00114 
00115   dmatrix L(1,n,1,n);
00116 #ifndef SAFE_INITIALIZE
00117     L.initialize();
00118 #endif
00119 
00120   int i,j,k;
00121   double tmp;
00122     if (M(1,1)<=0.0)
00123     {
00124       ierr =1;
00125       L.rowshift(rowsave);
00126       L.colshift(colsave);
00127       return L;
00128     }
00129   L(1,1)=sqrt(M(1,1));
00130   for (i=2;i<=n;i++)
00131   {
00132     L(i,1)=M(i,1)/L(1,1);
00133   }
00134 
00135   for (i=2;i<=n;i++)
00136   {
00137     for (j=2;j<=i-1;j++)
00138     {
00139       tmp=M(i,j);
00140       for (k=1;k<=j-1;k++)
00141       {
00142         tmp-=L(i,k)*L(j,k);
00143       }
00144       L(i,j)=tmp/L(j,j);
00145     }
00146     tmp=M(i,i);
00147     for (k=1;k<=i-1;k++)
00148     {
00149       tmp-=L(i,k)*L(i,k);
00150     }
00151     if (tmp<=0.0)
00152     {
00153       ierr =1;
00154       L.rowshift(rowsave);
00155       L.colshift(colsave);
00156       return L;
00157     }
00158     L(i,i)=sqrt(tmp);
00159   }
00160   L.rowshift(rowsave);
00161   L.colshift(colsave);
00162 
00163   return L;
00164 }
00165 
00170 lower_triangular_dmatrix  lower_triangular_choleski_decomp_positive
00171   (const dmatrix& MM,const int& _ierr)
00172 {
00173   // kludge to deal with constantness
00174   int & ierr = (int &)(_ierr);
00175   ierr=0;
00176   dmatrix & M= * (dmatrix *) &MM;
00177   if (M.colsize() != M.rowsize())
00178   {
00179     cerr << "Error in chol_decomp. Matrix not square" << endl;
00180     ad_exit(1);
00181   }
00182   int rowsave=M.rowmin();
00183   int colsave=M.colmin();
00184   M.rowshift(1);
00185   M.colshift(1);
00186   int n=M.rowmax();
00187 
00188   lower_triangular_dmatrix L(1,n);
00189 #ifndef SAFE_INITIALIZE
00190     L.initialize();
00191 #endif
00192 
00193   int i,j,k;
00194   double tmp;
00195     if (M(1,1)<=0.0)
00196     {
00197       ierr =1;
00198       L.rowshift(rowsave);
00199       L.colshift(colsave);
00200       return L;
00201     }
00202   L(1,1)=sqrt(M(1,1));
00203   for (i=2;i<=n;i++)
00204   {
00205     L(i,1)=M(i,1)/L(1,1);
00206   }
00207 
00208   for (i=2;i<=n;i++)
00209   {
00210     for (j=2;j<=i-1;j++)
00211     {
00212       tmp=M(i,j);
00213       for (k=1;k<=j-1;k++)
00214       {
00215         tmp-=L(i,k)*L(j,k);
00216       }
00217       L(i,j)=tmp/L(j,j);
00218     }
00219     tmp=M(i,i);
00220     for (k=1;k<=i-1;k++)
00221     {
00222       tmp-=L(i,k)*L(i,k);
00223     }
00224     if (tmp<=0.0)
00225     {
00226       ierr =1;
00227       L.rowshift(rowsave);
00228       L.colshift(colsave);
00229       return L;
00230     }
00231     L(i,i)=sqrt(tmp);
00232   }
00233   return L;
00234 }