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