ADMB Documentation  11.1x.2716
 All Classes Files Functions Variables Typedefs Friends Defines
dmat43.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: dmat43.cpp 2591 2014-11-07 22:20:29Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include "fvar.hpp"
00012 
00017 banded_lower_triangular_dmatrix::
00018   banded_lower_triangular_dmatrix(const banded_lower_triangular_dmatrix& mm) :
00019   bw(mm.bw), d(mm.d)
00020 {}
00021 
00026 banded_lower_triangular_dmatrix & banded_lower_triangular_dmatrix::operator =
00027   (const banded_lower_triangular_dmatrix& mm)
00028 {
00029   if (mm.bw!=bw)
00030   {
00031     cerr << "shape error" << endl;
00032     ad_exit(1);
00033   }
00034   else
00035   {
00036     for (int i=0;i<=bw-1;i++)
00037     {
00038       d(i)=mm.d(i);
00039     }
00040   }
00041   return *this;
00042 }
00043 
00048 banded_lower_triangular_dmatrix choleski_decomp_trust_bound(
00049   const banded_symmetric_dmatrix& _M,const int& _ierr)
00050 {
00051   int & ierr = (int &) _ierr;
00052   ADUNCONST(banded_symmetric_dmatrix,M)
00053   int minsave=M.indexmin();
00054   M.shift(1);
00055   int n=M.indexmax();
00056 
00057   double delta=0.0;
00058   int bw=M.bandwidth();
00059   banded_lower_triangular_dmatrix L(1,n,bw);
00060 #ifndef SAFE_INITIALIZE
00061     L.initialize();
00062 #endif
00063 
00064   int i;
00065   double tmp;
00066     if (M(1,1)<=0)
00067     {
00068       if (ierr==0)
00069         cerr << "Error matrix not positive definite in choleski_decomp"
00070           <<endl;
00071       ierr=1;
00072       return L;
00073     }
00074   L(1,1)=sqrt(M(1,1));
00075   for (i=2;i<=bw;i++)
00076   {
00077     L(i,1)=M(i,1)/L(1,1);
00078   }
00079 
00080   for (i=2;i<=n;i++)
00081   {
00082     for (int j=i-bw+1;j<=i-1;j++)
00083     {
00084       if (j>1)
00085       {
00086         tmp=M(i,j);
00087         for (int k=i-bw+1;k<=j-1;k++)
00088         {
00089           if (k>0 && k>j-bw)
00090             tmp-=L(i,k)*L(j,k);
00091         }
00092         L(i,j)=tmp/L(j,j);
00093       }
00094     }
00095     tmp=M(i,i);
00096     for (int k=i-bw+1;k<=i-1;k++)
00097     {
00098       if (k>0)
00099         tmp-=L(i,k)*L(i,k);
00100     }
00101     if (tmp<=0)
00102     {
00103       delta=-tmp;
00104       ierr=1;
00105       break;
00106     }
00107     L(i,i)=sqrt(tmp);
00108     if (i==n)
00109     {
00110       ierr=0;
00111     }
00112   }
00113   dvector v(1,n);
00114   if (ierr==1)
00115   {
00116     int k=i;
00117     v.initialize();
00118     v(k)=1.0;
00119     for (int i=k-1;i>=1;i--)
00120     {
00121       double ssum=0.0;
00122       int jmax=admin(n,i+bw-1);
00123       for (int j=i+1;j<=jmax;j++)
00124       {
00125         ssum+=L(j,i)*v(j);
00126       }
00127       v(i)=-ssum/L(i,i);
00128     }
00129     L(1,1)=delta/norm2(v);
00130   }
00131 
00132   M.shift(minsave);
00133   L.shift(minsave);
00134 
00135   return L;
00136 }
00137 //***********************************************************
00138 //***********************************************************
00139 /*
00140   int i,j;
00141 
00142   for (i=mmax;i>=mmin;i--)
00143   {
00144     double sum=0.0;
00145     int jmax=admin(mmax,i+bw-1);
00146     for (j=i+1;j<=jmax;j++)
00147     {
00148       sum+=M(j,i)*x(j);
00149     }
00150     x(i)=(y(i)-sum)/M(i,i);
00151   }
00152 
00153   return x;
00154 }
00155 */
00156 //***********************************************************
00157 //***********************************************************
00158 //***********************************************************
00159