ADMB Documentation  11.5.3254
 All Classes Files Functions Variables Typedefs Friends Defines
dmat28.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id$
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include "fvar.hpp"
00012 
00013 //void get_eigenv(const dvector& _d,const dvector& _e,const dmatrix& _z);
00014 
00015 #if !defined(OPT_LIB)
00016 
00021 dvector banded_symmetric_dmatrix::operator () (int i)
00022 {
00023   return d(i);
00024 }
00025 
00030 const dvector banded_symmetric_dmatrix::operator()(int i) const
00031 {
00032   return d(i);
00033 }
00034 
00039 const double& banded_symmetric_dmatrix::operator()(int i, int j) const
00040 {
00041   return d(i-j,i);
00042 }
00043 
00048 double& banded_symmetric_dmatrix::operator () (int i,int j)
00049 {
00050   return d(i-j,i);
00051 }
00052 
00057 dvector banded_lower_triangular_dmatrix::operator () (int i)
00058 {
00059   return d(i);
00060 }
00061 
00066 const dvector banded_lower_triangular_dmatrix::operator()(int i) const
00067 {
00068   return d(i);
00069 }
00070 
00075 double& banded_lower_triangular_dmatrix::operator () (int i,int j)
00076 {
00077   return d(i-j,i);
00078 }
00079 
00084 const double& banded_lower_triangular_dmatrix::operator()(int i, int j) const
00085 {
00086   return d(i-j,i);
00087 }
00088 
00089 #endif
00090 
00095 banded_symmetric_dmatrix::banded_symmetric_dmatrix(
00096   const banded_symmetric_dmatrix& BB,int _lb,int _ub) : bw(BB.bw), d(0,BB.bw-1)
00097 {
00098   banded_symmetric_dmatrix& B= (banded_symmetric_dmatrix&) BB;
00099   if (_lb<B.indexmin() || _ub>B.indexmax())
00100   {
00101     cerr << "bounds error" << endl;
00102     ad_exit(1);
00103   }
00104   d.allocate(0,bw-1);
00105   for (int i=0;i<=bw-1;i++)
00106   {
00107     d(i)=B.d(i)(i+_lb,_ub);
00108   }
00109 }
00110 
00115 void banded_symmetric_dmatrix::shift(int j)
00116 {
00117   for (int i=0;i<bw;i++)
00118     d(i).shift(j+i);
00119 }
00120 
00125 void banded_lower_triangular_dmatrix::shift(int j)
00126 {
00127   for (int i=0;i<bw;i++)
00128     d(i).shift(j+i);
00129 }
00130 
00135  banded_symmetric_dmatrix::banded_symmetric_dmatrix
00136   (int _min,int _max,int _bw)
00137 {
00138   bw=_bw;
00139   ivector  lb(0,bw-1);
00140   lb.fill_seqadd(_min,1);
00141   d.allocate(0,bw-1,lb,_max);
00142 }
00143 
00148 banded_symmetric_dmatrix::banded_symmetric_dmatrix(
00149   const dvar_matrix_position& pos)
00150 {
00151   int nrl=pos.row_min;
00152   int nrh=pos.row_max;
00153   int cmin=pos.lb(nrl);
00154   int cmax=pos.ub(nrl);
00155   bw=nrh-nrl+1;
00156   ivector lb(nrl,nrh);
00157   lb.fill_seqadd(cmin,1);
00158   d.allocate(nrl,nrh,lb,cmax);
00159 }
00160 
00165 banded_lower_triangular_dmatrix::banded_lower_triangular_dmatrix(
00166   const dvar_matrix_position& pos)
00167 {
00168   int nrl=pos.row_min;
00169   int nrh=pos.row_max;
00170   int cmin=pos.lb(nrl);
00171   int cmax=pos.ub(nrl);
00172   bw=nrh-nrl+1;
00173   ivector  lb(nrl,nrh);
00174   lb.fill_seqadd(cmin,1);
00175   d.allocate(nrl,nrh,lb,cmax);
00176 }
00177 
00182 dmatrix::dmatrix(const banded_lower_triangular_dmatrix& S)
00183 {
00184   int imin=S.indexmin();
00185   int imax=S.indexmax();
00186   int bw=S.bandwidth();
00187   allocate(imin,imax,imin,imax);
00188   for (int i=imin;i<=imax;i++)
00189   {
00190     for (int j=imin;j<=imax;j++)
00191     {
00192       if (j<=i)
00193       {
00194         if ( (i-j) < bw)
00195           (*this)(i,j)=S(i,j);
00196         else
00197           (*this)(i,j)=0.0;
00198       }
00199       else
00200       {
00201         (*this)(i,j)=0.0;
00202       }
00203     }
00204   }
00205 }
00206 
00211 dmatrix::dmatrix(const banded_symmetric_dmatrix& S)
00212 {
00213   int imin=S.indexmin();
00214   int imax=S.indexmax();
00215   int bw=S.bandwidth();
00216   allocate(imin,imax,imin,imax);
00217   int i1;
00218   int j1;
00219   for (int i=imin;i<=imax;i++)
00220   {
00221     for (int j=imin;j<=imax;j++)
00222     {
00223       if (j<=i)
00224       {
00225         j1=j;
00226         i1=i;
00227       }
00228       else
00229       {
00230         j1=i;
00231         i1=j;
00232       }
00233       if ( (i1-j1) < bw)
00234         (*this)(i,j)=S(i1,j1);
00235       else
00236         (*this)(i,j)=0.0;
00237     }
00238   }
00239 }
00240 
00245 banded_lower_triangular_dmatrix::banded_lower_triangular_dmatrix
00246   (int _min,int _max,int _bw)
00247 {
00248   bw=_bw;
00249   ivector  lb(0,_bw-1);
00250   lb.fill_seqadd(_min,1);
00251   d.allocate(0,_bw-1,lb,_max);
00252 }
00253 
00258 ostream& operator<<(const ostream& _ofs,
00259   const banded_lower_triangular_dmatrix& S1)
00260 {
00261   ostream & ofs = (ostream&) _ofs;
00262   banded_lower_triangular_dmatrix& S=(banded_lower_triangular_dmatrix&)(S1);
00263   int imin=S.indexmin();
00264   int imax=S.indexmax();
00265   int bw=S.bandwidth();
00266   for (int i=imin;i<=imax;i++)
00267   {
00268     for (int j=imin;j<=imax;j++)
00269     {
00270       if (j<=i)
00271       {
00272         if ( (i-j) < bw)
00273           ofs << S(i,j) << " ";
00274         else
00275           ofs << 0.0 << " ";
00276       }
00277       else
00278       {
00279         ofs << 0.0 << " ";
00280       }
00281     }
00282     if (i<imax) ofs << endl;
00283   }
00284   return ofs;
00285 }
00286 
00287 banded_lower_triangular_dmatrix choleski_decomp(
00288   const banded_symmetric_dmatrix& MM)
00289 {
00290   int ierr = 0;
00291   return choleski_decomp(MM, ierr);
00292 }
00293 banded_lower_triangular_dmatrix choleski_decomp(
00294   const banded_symmetric_dmatrix& _M, int& ierr)
00295 {
00296   ADUNCONST(banded_symmetric_dmatrix,M)
00297   int minsave=M.indexmin();
00298   M.shift(1);
00299   int n=M.indexmax();
00300 
00301   int bw=M.bandwidth();
00302   banded_lower_triangular_dmatrix L(1,n,bw);
00303 #ifndef SAFE_INITIALIZE
00304     L.initialize();
00305 #endif
00306 
00307   int i,j,k;
00308   double tmp;
00309     if (M(1,1)<=0)
00310     {
00311       if (ierr==0)
00312         cerr << "Error matrix not positive definite in choleski_decomp"
00313           <<endl;
00314       ierr=1;
00315       return L;
00316     }
00317   L(1,1)=sqrt(M(1,1));
00318   for (i=2;i<=bw;i++)
00319   {
00320     L(i,1)=M(i,1)/L(1,1);
00321   }
00322 
00323   for (i=2;i<=n;i++)
00324   {
00325     for (j=i-bw+1;j<=i-1;j++)
00326     {
00327       if (j>1)
00328       {
00329         tmp=M(i,j);
00330         for (k=i-bw+1;k<=j-1;k++)
00331         {
00332           if (k>0 && k>j-bw)
00333             tmp-=L(i,k)*L(j,k);
00334         }
00335         L(i,j)=tmp/L(j,j);
00336       }
00337     }
00338     tmp=M(i,i);
00339     for (k=i-bw+1;k<=i-1;k++)
00340     {
00341       if (k>0)
00342         tmp-=L(i,k)*L(i,k);
00343     }
00344     if (tmp<=0)
00345     {
00346       if (ierr==0)
00347         cerr << "Error matrix not positive definite in choleski_decomp"
00348           <<endl;
00349       ierr=1;
00350       return L;
00351     }
00352     L(i,i)=sqrt(tmp);
00353   }
00354   M.shift(minsave);
00355   L.shift(minsave);
00356 
00357   return L;
00358 }
00359 
00364 banded_symmetric_dmatrix& banded_symmetric_dmatrix::operator =
00365   (const banded_symmetric_dmatrix& M)
00366 {
00367   int _bw=M.bandwidth();
00368   int mmin=M.indexmin();
00369   int mmax=M.indexmax();
00370   if (!allocated(d))
00371   {
00372     bw=_bw;
00373     ivector  lb(0,_bw-1);
00374     lb.fill_seqadd(mmin,1);
00375     d.allocate(0,_bw-1,lb,mmax);
00376   }
00377   if (bw<bandwidth()) initialize();
00378 
00379   if (M.indexmin() != indexmin() || M.indexmax() != indexmax()
00380     || M.bandwidth() > bandwidth() )
00381   {
00382     cerr << "incompatible shape in symmetric_dmatrix::operator = "
00383          << endl;
00384     ad_exit(1);
00385   }
00386 
00387   for (int i=mmin;i<=mmax;i++)
00388   {
00389     int jmin=admax(mmin,i-bw+1);
00390     for (int j=jmin;j<=i;j++)
00391     {
00392       (*this)(i,j)=M(i,j);
00393     }
00394   }
00395   return *this;
00396 }
00397 
00402   banded_symmetric_dmatrix banded_symmetric_dmatrix::sub(int l,int u)
00403   {
00404     return banded_symmetric_dmatrix(*this,l,u);
00405   }
00406 
00411   dvector eigenvalues(const banded_symmetric_dmatrix& _SS)
00412   {
00413     banded_symmetric_dmatrix& S = (banded_symmetric_dmatrix&) _SS;
00414     if (S.bandwidth() !=2)
00415     {
00416       cerr << "error bandwidth not equal 2" << endl;
00417       ad_exit(1);
00418     }
00419 
00420     int lb=S.indexmin();
00421     int ub=S.indexmax();
00422     int bw=S.bandwidth();
00423     dmatrix M(lb,ub,lb,ub);
00424     M.initialize();
00425 
00426     for(int i=lb;i<=ub;i++)
00427     {
00428       for(int j=i;j<=min(bw+i-1,ub);j++)
00429       {
00430         M(j,i) = S(j,i);
00431         if(i!=j) M(i,j)=M(j,i);
00432       }
00433     }
00434 
00435     return eigenvalues(M);
00436   }
00437 
00442   dmatrix eigenvectors(const banded_symmetric_dmatrix& _SS,const dvector& _e)
00443   {
00444     banded_symmetric_dmatrix& S = (banded_symmetric_dmatrix&) _SS;
00445     if (S.bandwidth() !=2)
00446     {
00447       cerr << "error bandwidth not equal 2" << endl;
00448       ad_exit(1);
00449     }
00450 
00451     int lb=S.indexmin();
00452     int ub=S.indexmax();
00453     int bw=S.bandwidth();
00454     dmatrix M(lb,ub,lb,ub);
00455     M.initialize();
00456 
00457     for(int i=lb;i<=ub;i++)
00458     {
00459       for(int j=i;j<=min(bw+i-1,ub);j++)
00460       {
00461         M(j,i) = S(j,i);
00462         if(i!=j) M(i,j)=M(j,i);
00463       }
00464     }
00465 
00466     return eigenvectors(M);
00467   }
00468 /*
00469   dvector eigenvalues(const banded_symmetric_dmatrix& _SS)
00470   {
00471     banded_symmetric_dmatrix& S = (banded_symmetric_dmatrix&) _SS;
00472     if (S.bandwidth() !=2)
00473     {
00474       cerr << "error bandwidth not equal 2" << endl;
00475       ad_exit(1);
00476     }
00477     int mmin=S.indexmin();
00478     int mmax=S.indexmax();
00479     dvector diag(mmin,mmax);
00480     dvector offdiag(mmin,mmax);
00481     diag=S(0);
00482     offdiag(mmin+1,mmax)=S(1);
00483     offdiag(mmin)=0.0;
00484 
00485     return get_eigen_values(diag.shift(1),offdiag.shift(1));
00486   }
00487 
00488   dmatrix eigenvectors(const banded_symmetric_dmatrix& _SS,const dvector& _e)
00489   {
00490     banded_symmetric_dmatrix& S = (banded_symmetric_dmatrix&) _SS;
00491     ADUNCONST(dvector,e);
00492     if (S.bandwidth() !=2)
00493     {
00494       cerr << "error bandwidth not equal 2" << endl;
00495       ad_exit(1);
00496     }
00497     int minsave=S.indexmin();
00498     S.shift(1);
00499     e.shift(1);
00500     int mmin=1;
00501     int mmax=S.indexmax();
00502     dvector diag(mmin,mmax);
00503     dvector offdiag(mmin,mmax);
00504     diag=S(0);
00505     offdiag(mmin+1,mmax)=S(1);
00506     offdiag(mmin)=0.0;
00507     dmatrix z(mmin,mmax,mmin,mmax);
00508     z.initialize();
00509     for (int i=mmin;i<=mmax;i++)
00510     {
00511       z(i,i)=1.0;
00512     }
00513 
00514     get_eigenv(diag,offdiag,z);
00515     e=diag;
00516     e.shift(minsave);
00517     S.shift(minsave);
00518     z.colshift(minsave);
00519     z.rowshift(minsave);
00520     return z;
00521   }
00522 */