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