ADMB Documentation  11.1.2528
 All Classes Files Functions Variables Typedefs Friends Defines
eigen.cpp
Go to the documentation of this file.
00001 
00008 //#define eigen_vectors
00009 
00010 #include <fvar.hpp>
00011 
00012 
00013 void tri_dag(const dmatrix& m, const dvector& d,const dvector& e);
00014 void get_eigen(const dvector& d, const dvector& e, const dmatrix& z);
00015 
00020 dvector eigenvalues(const dmatrix& m)
00021 {
00022   if (m.rowsize()!=m.colsize())
00023   {
00024     cerr << "error -- non square matrix passed to "
00025     "dvector eigen(const dmatrix& m)\n";
00026     ad_exit(1);
00027   }
00028   dmatrix m1=symmetrize(m);
00029   m1.colshift(1);     // set minimum column and row indices to 1
00030   m1.rowshift(1);
00031   int n=m1.rowsize();
00032   dvector diag(1,n);
00033   dvector off_diag(1,n);
00034 
00035   tri_dag(m1,diag,off_diag);
00036 
00037   get_eigen(diag,off_diag,m1); // eigenvalues are returned in diag
00038            // eigenvalues are returned in columns of z
00039   return diag;
00040 }
00041 
00052 void tri_dag(const dmatrix& _m, const dvector& _d, const dvector& _e)
00053 {
00054   dvector& d = (dvector&) _d;
00055   dvector& e = (dvector&) _e;
00056   dmatrix& m = (dmatrix&) _m;
00057   if (m.rowsize() != m.colsize())
00058   {
00059     cerr << "Error -- non square matrix passed to "
00060     "void tridag(const dmatrix& m)\n";
00061     ad_exit(1);
00062   }
00063   if (m.rowsize() != d.size() || m.rowsize() != e.size()
00064     || d.indexmin() != 1 || e.indexmin() !=1 )
00065   {
00066     cerr <<"Error -- incorrect vector size passed to "
00067     "void tridag(const dmatrix& m)\n";
00068     ad_exit(1);
00069   }
00070   int n=m.rowsize();
00071   int l,j,i;
00072   double scale,hh,h,g,f;
00073 
00074   for (i=n;i>=2;i--)
00075   {
00076     l=i-1;
00077     h=scale=0.0;
00078     if (l > 1)
00079     {
00080       for (int k=1;k<=l;k++)
00081         scale += fabs(m[i][k]);
00082       if (scale == 0.0)
00083         e[i]=m[i][l];
00084       else
00085       {
00086         for (int k=1;k<=l;k++)
00087         {
00088           m[i][k] /= scale;
00089           h += m[i][k]*m[i][k];
00090         }
00091         f=m[i][l];
00092         g = f>0 ? -sqrt(h) : sqrt(h);
00093         e[i]=scale*g;
00094         h -= f*g;
00095         m[i][l]=f-g;
00096         f=0.0;
00097         for (j=1;j<=l;j++)
00098         {
00099         #ifdef EIGEN_VECTORS
00100         /* Next statement can be omitted if eigenvectors not wanted */
00101           m[j][i]=m[i][j]/h;
00102         #endif
00103           g=0.0;
00104           for (int k=1;k<=j;k++)
00105             g += m[j][k]*m[i][k];
00106           for (int k=j+1;k<=l;k++)
00107             g += m[k][j]*m[i][k];
00108           e[j]=g/h;
00109           f += e[j]*m[i][j];
00110         }
00111         hh=f/(h+h);
00112         for (j=1;j<=l;j++)
00113         {
00114           f=m[i][j];
00115           e[j]=g=e[j]-hh*f;
00116           for (int k=1;k<=j;k++)
00117             m[j][k] -= (f*e[k]+g*m[i][k]);
00118         }
00119       }
00120     }
00121     else
00122     {
00123       e[i]=m[i][l];
00124     }
00125     d[i]=h;
00126   }
00127   /* Next statement can be omitted if eigenvectors not wanted */
00128   d[1]=0.0;
00129   e[1]=0.0;
00130   /* Contents of this loop can be omitted if eigenvectors not
00131       wanted except for statement d[i]=a[i][i]; */
00132   #ifdef EIGEN_VECTORS
00133   for (i=1;i<=n;i++)
00134   {
00135     l=i-1;
00136     if (d[i])
00137     {
00138       for (j=1;j<=l;j++)
00139       {
00140         g=0.0;
00141         for (int k=1;k<=l;k++)
00142           g += m[i][k]*m[k][j];
00143         for (int k=1;k<=l;k++)
00144           m[k][j] -= g*m[k][i];
00145       }
00146     }
00147     d[i]=m[i][i];
00148     m[i][i]=1.0;
00149     for (j=1;j<=l;j++) m[j][i]=m[i][j]=0.0;
00150   }
00151   #else
00152   for (i=1;i<=n;i++)
00153   {
00154     d[i]=m[i][i];
00155   }
00156   #endif
00157 }
00158 
00159 double SIGN(const double x, double y)
00160 {
00161   if (y<0)
00162   {
00163     return -fabs(x);
00164   }
00165   else
00166   {
00167     return fabs(x);
00168   }
00169 }
00170 //#define SIGN(a,b) ((b)<0 ? -fabs(a) : fabs(a))
00171 
00182 void get_eigen(const dvector& _d, const dvector& _e, const dmatrix& _z)
00183 {
00184   dvector& d = (dvector&) _d;
00185   dvector& e = (dvector&) _e;
00186 #ifdef EIGEN_VECTORS
00187   dmatrix& z = (dmatrix&) _z;
00188 #endif
00189   int max_iterations=30;
00190   int n=d.size();
00191   max_iterations+=10*(n/100);
00192   int m,l,iter,i;
00193   double s,r,p,g,f,dd,c,b;
00194 
00195   for (i=2;i<=n;i++) e[i-1]=e[i];
00196   e[n]=0.0;
00197   for (l=1;l<=n;l++) {
00198     iter=0;
00199     do {
00200       for (m=l;m<=n-1;m++) {
00201         dd=fabs(d[m])+fabs(d[m+1]);
00202         if (fabs(e[m])+dd == dd) break;
00203       }
00204       if (m != l)
00205       {
00206         if (iter++ == max_iterations)
00207         {
00208           cerr << "Maximum number of iterations exceeded in"
00209           " dvector eigen(const dmatrix& m)\n";
00210           ad_exit(1);
00211         }
00212         g=(d[l+1]-d[l])/(2.0*e[l]);
00213         r=sqrt((g*g)+1.0);
00214         g=d[m]-d[l]+e[l]/(g+SIGN(r,g));
00215         s=c=1.0;
00216         p=0.0;
00217         for (i=m-1;i>=l;i--) {
00218           f=s*e[i];
00219           b=c*e[i];
00220           if (fabs(f) >= fabs(g)) {
00221             c=g/f;
00222             r=sqrt((c*c)+1.0);
00223             e[i+1]=f*r;
00224             c *= (s=1.0/r);
00225           } else {
00226             s=f/g;
00227             r=sqrt((s*s)+1.0);
00228             e[i+1]=g*r;
00229             s *= (c=1.0/r);
00230           }
00231           g=d[i+1]-p;
00232           r=(d[i]-g)*s+2.0*c*b;
00233           p=s*r;
00234           d[i+1]=g+p;
00235           g=c*r-b;
00236           /* Next loop can be omitted if eigenvectors not wanted */
00237           #ifdef EIGEN_VECTORS
00238             for (int k=1;k<=n;k++)
00239             {
00240               f=z[k][i+1];
00241               z[k][i+1]=s*z[k][i]+c*f;
00242               z[k][i]=c*z[k][i]-s*f;
00243             }
00244           #endif
00245         }
00246         d[l]=d[l]-p;
00247         e[l]=g;
00248         e[m]=0.0;
00249       }
00250     } while (m != l);
00251   }
00252 }
00253 
00264 dvector get_eigen_values(const dvector& _d,const dvector& _e)
00265 {
00266   dvector& d = (dvector&) _d;
00267   dvector& e = (dvector&) _e;
00268 
00269   int max_iterations=30;
00270   int n=d.size();
00271   max_iterations+=10*(n/100);
00272   int m,l,iter,i;
00273   double s,r,p,g,f,dd,c,b;
00274 
00275   for (i=2;i<=n;i++) e[i-1]=e[i];
00276   e[n]=0.0;
00277   for (l=1;l<=n;l++) {
00278     iter=0;
00279     do {
00280       for (m=l;m<=n-1;m++) {
00281         dd=fabs(d[m])+fabs(d[m+1]);
00282         if (fabs(e[m])+dd == dd) break;
00283       }
00284       if (m != l)
00285       {
00286         if (iter++ == max_iterations)
00287         {
00288           cerr << "Maximum number of iterations exceeded in"
00289           " dvector eigen(const dmatrix& m)\n";
00290           ad_exit(1);
00291         }
00292         g=(d[l+1]-d[l])/(2.0*e[l]);
00293         r=sqrt((g*g)+1.0);
00294         g=d[m]-d[l]+e[l]/(g+SIGN(r,g));
00295         s=c=1.0;
00296         p=0.0;
00297         for (i=m-1;i>=l;i--) {
00298           f=s*e[i];
00299           b=c*e[i];
00300           if (fabs(f) >= fabs(g)) {
00301             c=g/f;
00302             r=sqrt((c*c)+1.0);
00303             e[i+1]=f*r;
00304             c *= (s=1.0/r);
00305           } else {
00306             s=f/g;
00307             r=sqrt((s*s)+1.0);
00308             e[i+1]=g*r;
00309             s *= (c=1.0/r);
00310           }
00311           g=d[i+1]-p;
00312           r=(d[i]-g)*s+2.0*c*b;
00313           p=s*r;
00314           d[i+1]=g+p;
00315           g=c*r-b;
00316         }
00317         d[l]=d[l]-p;
00318         e[l]=g;
00319         e[m]=0.0;
00320       }
00321     } while (m != l);
00322   }
00323   return d;
00324 }
00325 
00337 dvector get_eigen_values(const dvector& _d,const dvector& _e,
00338   const dmatrix& _z)
00339   //eigenvectors are returned in z
00340 {
00341   dvector& xd = (dvector&) _d;
00342   dvector& xe = (dvector&) _e;
00343   dmatrix& z = (dmatrix&) _z;
00344 
00345   dvector d(xd.indexmin(),xd.indexmax());
00346   dvector e(xe.indexmin(),xe.indexmax());
00347 
00348   d=xd;
00349   e=xe;
00350 
00351   int max_iterations=30;
00352   int n=d.size();
00353   max_iterations+=10*(n/100);
00354   int m,l,iter,i;
00355   double s,r,p,g,f,dd,c,b;
00356 
00357   for (i=2;i<=n;i++) e[i-1]=e[i];
00358   e[n]=0.0;
00359   for (l=1;l<=n;l++) {
00360     iter=0;
00361     do {
00362       for (m=l;m<=n-1;m++) {
00363         dd=fabs(d[m])+fabs(d[m+1]);
00364         if (fabs(e[m])+dd == dd) break;
00365       }
00366       if (m != l)
00367       {
00368         if (iter++ == max_iterations)
00369         {
00370           cerr << "Maximum number of iterations exceeded in"
00371           " dvector eigen(const dmatrix& m)\n";
00372           ad_exit(1);
00373         }
00374         g=(d[l+1]-d[l])/(2.0*e[l]);
00375         r=sqrt((g*g)+1.0);
00376         g=d[m]-d[l]+e[l]/(g+SIGN(r,g));
00377         s=c=1.0;
00378         p=0.0;
00379         for (i=m-1;i>=l;i--) {
00380           f=s*e[i];
00381           b=c*e[i];
00382           if (fabs(f) >= fabs(g)) {
00383             c=g/f;
00384             r=sqrt((c*c)+1.0);
00385             e[i+1]=f*r;
00386             c *= (s=1.0/r);
00387           } else {
00388             s=f/g;
00389             r=sqrt((s*s)+1.0);
00390             e[i+1]=g*r;
00391             s *= (c=1.0/r);
00392           }
00393           g=d[i+1]-p;
00394           r=(d[i]-g)*s+2.0*c*b;
00395           p=s*r;
00396           d[i+1]=g+p;
00397           g=c*r-b;
00398           /* Next loop can be omitted if eigenvectors not wanted */
00399           for (int k=1;k<=n;k++)
00400           {
00401             f=z[k][i+1];
00402             z[k][i+1]=s*z[k][i]+c*f;
00403             z[k][i]=c*z[k][i]-s*f;
00404           }
00405         }
00406         d[l]=d[l]-p;
00407         e[l]=g;
00408         e[m]=0.0;
00409       }
00410     } while (m != l);
00411   }
00412   return d;
00413 }
00414 #undef EIGEN_VECTORS