ADMB Documentation  11.2.2828
 All Classes Files Functions Variables Typedefs Friends Defines
eigenv.cpp
Go to the documentation of this file.
00001 
00011 #define EIGEN_VECTORS
00012 
00013 #include <fvar.hpp>
00014 
00015 #ifdef ISZERO
00016   #undef ISZERO
00017 #endif
00018 #define ISZERO(d) ((d)==0.0)
00019 
00020 #if !defined(EIGEN_VECTORS)
00021 #  define EIGEN_VECTORS
00022 #endif
00023 
00024 #ifdef EIGEN_VECTORS
00025   void tri_dagv(const dmatrix& m,const dvector& d,const dvector& e);
00026 #else
00027   void tri_dag(const dmatrix& m,const dvector& d,const dvector& e);
00028 #endif
00029 #ifdef EIGEN_VECTORS
00030   void get_eigenv(const dvector& d,const dvector& e,const dmatrix& z);
00031 #else
00032   void get_eigen(const dvector& d,const dvector& e,const dmatrix& z);
00033 #endif
00034 
00040 dmatrix eigenvectors(const dmatrix& m)  //,_CONST dvector& diag)
00041 {
00042   if (m.rowsize()!=m.colsize())
00043   {
00044     cerr << "Error -- non square matrix passed to dvector"
00045            " eigen(const dmatrix& m)" << endl;
00046     ad_exit(1);
00047   }
00048 
00049   dmatrix m1=symmetrize(m);
00050   int n=m1.rowsize();
00051   m1.colshift(1);     // set minimum column and row indices to 1
00052   m1.rowshift(1);
00053   dvector diag(1,n);
00054   dvector off_diag(1,n);
00055 
00056   #ifdef EIGEN_VECTORS
00057     tri_dagv(m1,diag,off_diag);
00058   #else
00059     tri_dag(m1,diag,off_diag);
00060   #endif
00061 
00062   #ifdef EIGEN_VECTORS
00063     get_eigenv(diag,off_diag,m1); // eigenvalues are returned in diag
00064   #else
00065     get_eigen(diag,off_diag,m1); // eigenvalues are returned in diag
00066   #endif
00067   return m1;
00068 }
00069 
00076 dmatrix eigenvectors(const dmatrix& m,const dvector& _diag)
00077   //,_CONST dvector& diag)
00078 {
00079   ADUNCONST(dvector,diag)
00080   if (m.rowsize()!=m.colsize())
00081   {
00082     cerr << "Error -- non square matrix passed to dvector"
00083            " eigen(_CONST dmatrix& m)" << endl;
00084     ad_exit(1);
00085   }
00086 
00087   dmatrix m1=symmetrize(m);
00088   int n=m1.rowsize();
00089   m1.colshift(1);     // set minimum column and row indices to 1
00090   m1.rowshift(1);
00091   diag.shift(1);
00092   if (diag.size()!=m.rowsize())
00093   {
00094     cerr << "incompatible sizes in vector and matrix passed to"
00095           " dmatrix eigenvector routine" << endl;
00096   }
00097   dvector off_diag(1,n);
00098 
00099   #ifdef EIGEN_VECTORS
00100     tri_dagv(m1,diag,off_diag);
00101   #else
00102     tri_dag(m1,diag,off_diag);
00103   #endif
00104 
00105   #ifdef EIGEN_VECTORS
00106     get_eigenv(diag,off_diag,m1); // eigenvalues are returned in diag
00107   #else
00108     get_eigen(diag,off_diag,m1); // eigenvalues are returned in diag
00109   #endif
00110   return m1;
00111 }
00112 
00113 
00124 #ifdef EIGEN_VECTORS
00125   void tri_dagv(const dmatrix& _m,const dvector& _d,const dvector& _e)
00126 #else
00127   void tri_dagv(const dmatrix& _m,const dvector& _d,const dvector& _e)
00128 #endif
00129 {
00130   dvector& d = (dvector&) _d;
00131   dvector& e = (dvector&) _e;
00132   dmatrix& m = (dmatrix&) _m;
00133   if (m.rowsize() != m.colsize())
00134   {
00135     cerr << "Error -- non square matrix passed to "
00136     "void tridag(const dmatrix& m)\n";
00137     ad_exit(1);
00138   }
00139   if (m.rowsize() != d.size() || m.rowsize() != e.size()
00140     || d.indexmin() != 1 || e.indexmin() !=1 )
00141   {
00142     cerr <<"Error -- incorrect vector size passed to "
00143     "void tridag(const dmatrix& m)\n";
00144     ad_exit(1);
00145   }
00146   int n=m.rowsize();
00147   int l,k,j,i;
00148   double scale,hh,h,g,f;
00149 
00150   for (i=n;i>=2;i--)
00151   {
00152     l=i-1;
00153     h=scale=0.0;
00154     if (l > 1)
00155     {
00156       for (k=1;k<=l;k++)
00157         scale += fabs(m[i][k]);
00158       if (scale == 0.0)
00159         e[i]=m[i][l];
00160       else
00161       {
00162         for (k=1;k<=l;k++)
00163         {
00164           m[i][k] /= scale;
00165           h += m[i][k]*m[i][k];
00166         }
00167         f=m[i][l];
00168         g = f>0 ? -sqrt(h) : sqrt(h);
00169         e[i]=scale*g;
00170         h -= f*g;
00171         m[i][l]=f-g;
00172         f=0.0;
00173         for (j=1;j<=l;j++)
00174         {
00175         #ifdef EIGEN_VECTORS
00176         /* Next statement can be omitted if eigenvectors not wanted */
00177           m[j][i]=m[i][j]/h;
00178         #endif
00179           g=0.0;
00180           for (k=1;k<=j;k++)
00181             g += m[j][k]*m[i][k];
00182           for (k=j+1;k<=l;k++)
00183             g += m[k][j]*m[i][k];
00184           e[j]=g/h;
00185           f += e[j]*m[i][j];
00186         }
00187         hh=f/(h+h);
00188         for (j=1;j<=l;j++)
00189         {
00190           f=m[i][j];
00191           e[j]=g=e[j]-hh*f;
00192           for (k=1;k<=j;k++)
00193             m[j][k] -= (f*e[k]+g*m[i][k]);
00194         }
00195       }
00196     }
00197     else
00198     {
00199       e[i]=m[i][l];
00200     }
00201     d[i]=h;
00202   }
00203   /* Next statement can be omitted if eigenvectors not wanted */
00204   d[1]=0.0;
00205   e[1]=0.0;
00206   /* Contents of this loop can be omitted if eigenvectors not
00207       wanted except for statement d[i]=a[i][i]; */
00208   #ifdef EIGEN_VECTORS
00209   for (i=1;i<=n;i++)
00210   {
00211     l=i-1;
00212     if (!ISZERO(d[i]))
00213     {
00214       for (j=1;j<=l;j++)
00215       {
00216         g=0.0;
00217         for (k=1;k<=l;k++)
00218           g += m[i][k]*m[k][j];
00219         for (k=1;k<=l;k++)
00220           m[k][j] -= g*m[k][i];
00221       }
00222     }
00223     d[i]=m[i][i];
00224     m[i][i]=1.0;
00225     for (j=1;j<=l;j++) m[j][i]=m[i][j]=0.0;
00226   }
00227   #endif
00228 }
00229 
00235 double SIGNV(const double x, double y)
00236 {
00237   if (y<0)
00238   {
00239     return -fabs(x);
00240   }
00241   else
00242   {
00243     return fabs(x);
00244   }
00245 }
00246 
00256 #ifdef EIGEN_VECTORS
00257   void get_eigenv(const dvector& _d,const dvector& _e,const dmatrix& _z)
00258 #else
00259   void get_eigen(const dvector& _d,const dvector& _e,const dmatrix& _z)
00260 #endif
00261 {
00262   dvector& d = (dvector&) _d;
00263   dvector& e = (dvector&) _e;
00264   dmatrix& z = (dmatrix&) _z;
00265   int n=d.size();
00266   int m,l,iter,i;
00267   double s,r,p,g,f,dd,c,b;
00268 
00269   for (i=2;i<=n;i++) e[i-1]=e[i];
00270   e[n]=0.0;
00271   for (l=1;l<=n;l++) {
00272     iter=0;
00273     do {
00274       for (m=l;m<=n-1;m++) {
00275         dd=fabs(d[m])+fabs(d[m+1]);
00276         if (fabs(e[m])+dd == dd) break;
00277       }
00278       if (m != l)
00279       {
00280         if (iter++ == 30)
00281         {
00282           cerr << "Maximum number of iterations exceeded in"
00283           " dvector eigen(const dmatrix& m)\n";
00284           ad_exit(1);
00285         }
00286         g=(d[l+1]-d[l])/(2.0*e[l]);
00287         r=sqrt((g*g)+1.0);
00288         g=d[m]-d[l]+e[l]/(g+SIGNV(r,g));
00289         s=c=1.0;
00290         p=0.0;
00291         for (i=m-1;i>=l;i--) {
00292           f=s*e[i];
00293           b=c*e[i];
00294           if (fabs(f) >= fabs(g)) {
00295             c=g/f;
00296             r=sqrt((c*c)+1.0);
00297             e[i+1]=f*r;
00298             c *= (s=1.0/r);
00299           } else {
00300             s=f/g;
00301             r=sqrt((s*s)+1.0);
00302             e[i+1]=g*r;
00303             s *= (c=1.0/r);
00304           }
00305           g=d[i+1]-p;
00306           r=(d[i]-g)*s+2.0*c*b;
00307           p=s*r;
00308           d[i+1]=g+p;
00309           g=c*r-b;
00310           /* Next loop can be omitted if eigenvectors not wanted */
00311           #ifdef EIGEN_VECTORS
00312             for (int k=1;k<=n;k++)
00313             {
00314               f=z[k][i+1];
00315               z[k][i+1]=s*z[k][i]+c*f;
00316               z[k][i]=c*z[k][i]-s*f;
00317             }
00318           #endif
00319         }
00320         d[l]=d[l]-p;
00321         e[l]=g;
00322         e[m]=0.0;
00323       }
00324     } while (m != l);
00325   }
00326 }
00327 #undef EIGEN_VECTORS