ADMB Documentation  11.1.2495
 All Classes Files Functions Variables Typedefs Friends Defines
dmat28.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: dmat28.cpp 2435 2014-09-30 23:43: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 //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 
00298 banded_lower_triangular_dmatrix choleski_decomp(
00299   const banded_symmetric_dmatrix& _M,const int& _ierr)
00300 {
00301   int & ierr = (int &) _ierr;
00302   ADUNCONST(banded_symmetric_dmatrix,M)
00303   int minsave=M.indexmin();
00304   M.shift(1);
00305   int n=M.indexmax();
00306 
00307   int bw=M.bandwidth();
00308   banded_lower_triangular_dmatrix L(1,n,bw);
00309 #ifndef SAFE_INITIALIZE
00310     L.initialize();
00311 #endif
00312 
00313   int i,j,k;
00314   double tmp;
00315     if (M(1,1)<=0)
00316     {
00317       if (ierr==0)
00318         cerr << "Error matrix not positive definite in choleski_decomp"
00319           <<endl;
00320       ierr=1;
00321       return L;
00322     }
00323   L(1,1)=sqrt(M(1,1));
00324   for (i=2;i<=bw;i++)
00325   {
00326     L(i,1)=M(i,1)/L(1,1);
00327   }
00328 
00329   for (i=2;i<=n;i++)
00330   {
00331     for (j=i-bw+1;j<=i-1;j++)
00332     {
00333       if (j>1)
00334       {
00335         tmp=M(i,j);
00336         for (k=i-bw+1;k<=j-1;k++)
00337         {
00338           if (k>0 && k>j-bw)
00339             tmp-=L(i,k)*L(j,k);
00340         }
00341         L(i,j)=tmp/L(j,j);
00342       }
00343     }
00344     tmp=M(i,i);
00345     for (k=i-bw+1;k<=i-1;k++)
00346     {
00347       if (k>0)
00348         tmp-=L(i,k)*L(i,k);
00349     }
00350     if (tmp<=0)
00351     {
00352       if (ierr==0)
00353         cerr << "Error matrix not positive definite in choleski_decomp"
00354           <<endl;
00355       ierr=1;
00356       return L;
00357     }
00358     L(i,i)=sqrt(tmp);
00359   }
00360   M.shift(minsave);
00361   L.shift(minsave);
00362 
00363   return L;
00364 }
00365 
00370 banded_symmetric_dmatrix& banded_symmetric_dmatrix::operator =
00371   (const banded_symmetric_dmatrix& M)
00372 {
00373   int _bw=M.bandwidth();
00374   int mmin=M.indexmin();
00375   int mmax=M.indexmax();
00376   if (!allocated(d))
00377   {
00378     bw=_bw;
00379     ivector  lb(0,_bw-1);
00380     lb.fill_seqadd(mmin,1);
00381     d.allocate(0,_bw-1,lb,mmax);
00382   }
00383   if (bw<bandwidth()) initialize();
00384 
00385   if (M.indexmin() != indexmin() || M.indexmax() != indexmax()
00386     || M.bandwidth() > bandwidth() )
00387   {
00388     cerr << "incompatible shape in symmetric_dmatrix::operator = "
00389          << endl;
00390     ad_exit(1);
00391   }
00392 
00393   for (int i=mmin;i<=mmax;i++)
00394   {
00395     int jmin=admax(mmin,i-bw+1);
00396     for (int j=jmin;j<=i;j++)
00397     {
00398       (*this)(i,j)=M(i,j);
00399     }
00400   }
00401   return *this;
00402 }
00403 
00408   banded_symmetric_dmatrix banded_symmetric_dmatrix::sub(int l,int u)
00409   {
00410     return banded_symmetric_dmatrix(*this,l,u);
00411   }
00412 
00417   dvector eigenvalues(const banded_symmetric_dmatrix& _SS)
00418   {
00419     banded_symmetric_dmatrix& S = (banded_symmetric_dmatrix&) _SS;
00420     if (S.bandwidth() !=2)
00421     {
00422       cerr << "error bandwidth not equal 2" << endl;
00423       ad_exit(1);
00424     }
00425 
00426     int lb=S.indexmin();
00427     int ub=S.indexmax();
00428     int bw=S.bandwidth();
00429     dmatrix M(lb,ub,lb,ub);
00430     M.initialize();
00431 
00432     for(int i=lb;i<=ub;i++)
00433     {
00434       for(int j=i;j<=min(bw+i-1,ub);j++)
00435       {
00436         M(j,i) = S(j,i);
00437         if(i!=j) M(i,j)=M(j,i);
00438       }
00439     }
00440 
00441     return eigenvalues(M);
00442   }
00443 
00448   dmatrix eigenvectors(const banded_symmetric_dmatrix& _SS,const dvector& _e)
00449   {
00450     banded_symmetric_dmatrix& S = (banded_symmetric_dmatrix&) _SS;
00451     if (S.bandwidth() !=2)
00452     {
00453       cerr << "error bandwidth not equal 2" << endl;
00454       ad_exit(1);
00455     }
00456 
00457     int lb=S.indexmin();
00458     int ub=S.indexmax();
00459     int bw=S.bandwidth();
00460     dmatrix M(lb,ub,lb,ub);
00461     M.initialize();
00462 
00463     for(int i=lb;i<=ub;i++)
00464     {
00465       for(int j=i;j<=min(bw+i-1,ub);j++)
00466       {
00467         M(j,i) = S(j,i);
00468         if(i!=j) M(i,j)=M(j,i);
00469       }
00470     }
00471 
00472     return eigenvectors(M);
00473   }
00474 /*
00475   dvector eigenvalues(const banded_symmetric_dmatrix& _SS)
00476   {
00477     banded_symmetric_dmatrix& S = (banded_symmetric_dmatrix&) _SS;
00478     if (S.bandwidth() !=2)
00479     {
00480       cerr << "error bandwidth not equal 2" << endl;
00481       ad_exit(1);
00482     }
00483     int mmin=S.indexmin();
00484     int mmax=S.indexmax();
00485     dvector diag(mmin,mmax);
00486     dvector offdiag(mmin,mmax);
00487     diag=S(0);
00488     offdiag(mmin+1,mmax)=S(1);
00489     offdiag(mmin)=0.0;
00490 
00491     return get_eigen_values(diag.shift(1),offdiag.shift(1));
00492   }
00493 
00494   dmatrix eigenvectors(const banded_symmetric_dmatrix& _SS,const dvector& _e)
00495   {
00496     banded_symmetric_dmatrix& S = (banded_symmetric_dmatrix&) _SS;
00497     ADUNCONST(dvector,e);
00498     if (S.bandwidth() !=2)
00499     {
00500       cerr << "error bandwidth not equal 2" << endl;
00501       ad_exit(1);
00502     }
00503     int minsave=S.indexmin();
00504     S.shift(1);
00505     e.shift(1);
00506     int mmin=1;
00507     int mmax=S.indexmax();
00508     dvector diag(mmin,mmax);
00509     dvector offdiag(mmin,mmax);
00510     diag=S(0);
00511     offdiag(mmin+1,mmax)=S(1);
00512     offdiag(mmin)=0.0;
00513     dmatrix z(mmin,mmax,mmin,mmax);
00514     z.initialize();
00515     for (int i=mmin;i<=mmax;i++)
00516     {
00517       z(i,i)=1.0;
00518     }
00519 
00520     get_eigenv(diag,offdiag,z);
00521     e=diag;
00522     e.shift(minsave);
00523     S.shift(minsave);
00524     z.colshift(minsave);
00525     z.rowshift(minsave);
00526     return z;
00527   }
00528 */