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