ADMB Documentation  11.1.2249
 All Classes Files Functions Variables Typedefs Friends Defines
dveigenv.cpp
Go to the documentation of this file.
00001 
00010 #define EIGEN_VECTORS
00011 
00012 #include <fvar.hpp>
00013 #ifdef EIGEN_VECTORS
00014   void tri_dagv(const dvar_matrix& ,const dvar_vector& , const dvar_vector& );
00015 #else
00016   void tri_dag(const dvar_matrix& m,const dvar_vector& d,const dvar_vector& e);
00017 #endif
00018 #ifdef EIGEN_VECTORS
00019   void get_eigenv(const dvar_vector& d, const dvar_vector& e,
00020     const dvar_matrix& z);
00021 #else
00022   void get_eigen(const dvar_vector& d, const dvar_vector& e,
00023     const dvar_matrix& z);
00024 #endif
00025 
00032 dvar_matrix eigenvectors(const dvar_matrix& m)
00033 {
00034   if (m.rowsize()!=m.colsize())
00035   {
00036     cerr << "Error -- non square matrix passed to "
00037     "dvector eigen(const dmatrix& m)\n";
00038     ad_exit(1);
00039   }
00040 
00041   dvar_matrix m1=symmetrize(m);
00042   int n=m1.rowsize();
00043   m1.colshift(1);     // set minimum column and row indices to 1
00044   m1.rowshift(1);
00045   dvar_vector diag(1,n);
00046   dvar_vector off_diag(1,n);
00047 
00048   #ifdef EIGEN_VECTORS
00049     tri_dagv(m1,diag,off_diag);
00050   #else
00051     tri_dag(m1,diag,off_diag);
00052   #endif
00053 
00054   #ifdef EIGEN_VECTORS
00055     get_eigenv(diag,off_diag,m1); // eigenvalues are returned in diag
00056   #else
00057     get_eigen(diag,off_diag,m1); // eigenvalues are returned in diag
00058   #endif
00059            // eigenvalues are returned in columns of z
00060   return m1;
00061 }
00062 
00074 #ifdef EIGEN_VECTORS
00075 void tri_dagv(const dvar_matrix& _m,const dvar_vector& _d,
00076   const dvar_vector& _e)
00077 #else
00078 void tri_dagv(const dvar_matrix& _m,const dvar_vector& _d,
00079   const dvar_vector& _e)
00080 #endif
00081 {
00082   ADUNCONST(dvar_vector,d)
00083   ADUNCONST(dvar_vector,e)
00084   dvar_matrix& m=(dvar_matrix&) _m;
00085   if (m.rowsize() != m.colsize())
00086   {
00087     cerr << "Error -- non square matrix passed to "
00088     "void tridag(const dvar_matrix& m)\n";
00089     ad_exit(1);
00090   }
00091   if (m.rowsize() != d.size() || m.rowsize() != e.size()
00092     || d.indexmin() != 1 || e.indexmin() !=1 )
00093   {
00094     cerr <<"Error -- incorrect vector size passed to "
00095     "void tridag(const dmatrix& m)\n";
00096     ad_exit(1);
00097   }
00098   int n=m.rowsize();
00099   int l,k,j,i;
00100   dvariable scale,hh,h,g,f;
00101 
00102   for (i=n;i>=2;i--)
00103   {
00104     l=i-1;
00105     h=scale=0.0;
00106     if (l > 1)
00107     {
00108       for (k=1;k<=l;k++)
00109         scale += fabs(m[i][k]);
00110       if (scale == 0.0)
00111         e[i]=m[i][l];
00112       else
00113       {
00114         for (k=1;k<=l;k++)
00115         {
00116           m[i][k] /= scale;
00117           h += m[i][k]*m[i][k];
00118         }
00119         f=m[i][l];
00120         g = f>0. ? -sqrt(h) : sqrt(h);
00121         e[i]=scale*g;
00122         h -= f*g;
00123         m[i][l]=f-g;
00124         f=0.0;
00125         for (j=1;j<=l;j++)
00126         {
00127         #ifdef EIGEN_VECTORS
00128         /* Next statement can be omitted if eigenvectors not wanted */
00129           m[j][i]=m[i][j]/h;
00130         #endif
00131           g=0.0;
00132           for (k=1;k<=j;k++)
00133             g += m[j][k]*m[i][k];
00134           for (k=j+1;k<=l;k++)
00135             g += m[k][j]*m[i][k];
00136           e[j]=g/h;
00137           f += e[j]*m[i][j];
00138         }
00139         hh=f/(h+h);
00140         for (j=1;j<=l;j++)
00141         {
00142           f=m[i][j];
00143           e[j]=g=e[j]-hh*f;
00144           for (k=1;k<=j;k++)
00145             m[j][k] -= (f*e[k]+g*m[i][k]);
00146         }
00147       }
00148     }
00149     else
00150     {
00151       e[i]=m[i][l];
00152     }
00153     d[i]=h;
00154   }
00155   /* Next statement can be omitted if eigenvectors not wanted */
00156   d[1]=0.0;
00157   e[1]=0.0;
00158   /* Contents of this loop can be omitted if eigenvectors not
00159       wanted except for statement d[i]=a[i][i]; */
00160   #ifdef EIGEN_VECTORS
00161   for (i=1;i<=n;i++)
00162   {
00163     l=i-1;
00164     if (value( d[i]) )
00165     {
00166       for (j=1;j<=l;j++)
00167       {
00168         g=0.0;
00169         for (k=1;k<=l;k++)
00170           g += m[i][k]*m[k][j];
00171         for (k=1;k<=l;k++)
00172           m[k][j] -= g*m[k][i];
00173       }
00174     }
00175     d[i]=m[i][i];
00176     m[i][i]=1.0;
00177     for (j=1;j<=l;j++) m[j][i]=m[i][j]=0.0;
00178   }
00179   #endif
00180 }
00181 
00187   dvariable SIGNV(const prevariable& x, const prevariable& y)
00188 {
00189   if (y<0.)
00190   {
00191     return -fabs(x);
00192   }
00193   else
00194   {
00195     return fabs(x);
00196   }
00197 }
00198 
00209 #ifdef EIGEN_VECTORS
00210   void get_eigenv(const dvar_vector& _d, const dvar_vector& _e,
00211     const dvar_matrix& _z)
00212 #else
00213   void get_eigen(const dvar_vector& _d, const dvar_vector& _e,
00214     const dvar_matrix& _z)
00215 #endif
00216 {
00217   ADUNCONST(dvar_matrix,z)
00218   dvar_vector& e=(dvar_vector&) _e;
00219   dvar_vector& d=(dvar_vector&) _d;
00220   int n=d.size();
00221   int m,l,iter,i,k;
00222   dvariable s,r,p,g,f,dd,c,b;
00223 
00224   for (i=2;i<=n;i++) e[i-1]=e[i];
00225   e[n]=0.0;
00226   for (l=1;l<=n;l++) {
00227     iter=0;
00228     do {
00229       for (m=l;m<=n-1;m++) {
00230         dd=fabs(d[m])+fabs(d[m+1]);
00231         if (fabs(e[m])+dd == dd) break;
00232       }
00233       if (m != l)
00234       {
00235         if (iter++ == 30)
00236         {
00237           cerr << "Maximum number of iterations exceeded in"
00238           " dvector eigen(const dmatrix& m)\n";
00239           ad_exit(1);
00240         }
00241         g=(d[l+1]-d[l])/(2.0*e[l]);
00242         r=sqrt((g*g)+1.0);
00243         g=d[m]-d[l]+e[l]/(g+SIGNV(r,g));
00244         s=c=1.0;
00245         p=0.0;
00246         for (i=m-1;i>=l;i--) {
00247           f=s*e[i];
00248           b=c*e[i];
00249           if (fabs(f) >= fabs(g)) {
00250             c=g/f;
00251             r=sqrt((c*c)+1.0);
00252             e[i+1]=f*r;
00253             c *= (s=1.0/r);
00254           } else {
00255             s=f/g;
00256             r=sqrt((s*s)+1.0);
00257             e[i+1]=g*r;
00258             s *= (c=1.0/r);
00259           }
00260           g=d[i+1]-p;
00261           r=(d[i]-g)*s+2.0*c*b;
00262           p=s*r;
00263           d[i+1]=g+p;
00264           g=c*r-b;
00265           /* Next loop can be omitted if eigenvectors not wanted */
00266           #ifdef EIGEN_VECTORS
00267             for (k=1;k<=n;k++)
00268             {
00269               f=z[k][i+1];
00270               z[k][i+1]=s*z[k][i]+c*f;
00271               z[k][i]=c*z[k][i]-s*f;
00272             }
00273           #endif
00274         }
00275         d[l]=d[l]-p;
00276         e[l]=g;
00277         e[m]=0.0;
00278       }
00279     } while (m != l);
00280   }
00281 }
00282 #undef EIGEN_VECTORS