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