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