ADMB Documentation  11.1x.2711
 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 (value(y) < 0)
00159   {
00160     return -fabs(x);
00161   }
00162   else
00163   {
00164     return fabs(x);
00165   }
00166 }
00167 //#define SIGN(a,b) ((b)<0 ? -fabs(a) : fabs(a))
00168 
00179 void get_eigen(const dvar_vector& _d,const dvar_vector& _e,
00180   const dvar_matrix& z)
00181 {
00182   ADUNCONST(dvar_vector,d)
00183   ADUNCONST(dvar_vector,e)
00184   int n=d.size();
00185   int m,l,iter,i;
00186   dvariable s,r,p,g,f,dd,c,b;
00187 
00188   for (i=2;i<=n;i++) e[i-1]=e[i];
00189   e[n]=0.0;
00190   for (l=1;l<=n;l++) {
00191     iter=0;
00192     do {
00193       for (m=l;m<=n-1;m++) {
00194         dd=fabs(d[m])+fabs(d[m+1]);
00195         if (fabs(e[m])+dd == dd) break;
00196       }
00197       if (m != l)
00198       {
00199         if (iter++ == 30)
00200         {
00201           cerr << "Maximum number of iterations exceeded in"
00202           " dvector eigen(const dmatrix& m)\n";
00203           ad_exit(1);
00204         }
00205         g=(d[l+1]-d[l])/(2.0*e[l]);
00206         r=sqrt((g*g)+1.0);
00207         g=d[m]-d[l]+e[l]/(g+SIGN(r,g));
00208         s=c=1.0;
00209         p=0.0;
00210         for (i=m-1;i>=l;i--) {
00211           f=s*e[i];
00212           b=c*e[i];
00213           if (fabs(f) >= fabs(g)) {
00214             c=g/f;
00215             r=sqrt((c*c)+1.0);
00216             e[i+1]=f*r;
00217             c *= (s=1.0/r);
00218           } else {
00219             s=f/g;
00220             r=sqrt((s*s)+1.0);
00221             e[i+1]=g*r;
00222             s *= (c=1.0/r);
00223           }
00224           g=d[i+1]-p;
00225           r=(d[i]-g)*s+2.0*c*b;
00226           p=s*r;
00227           d[i+1]=g+p;
00228           g=c*r-b;
00229           /* Next loop can be omitted if eigenvectors not wanted */
00230           #ifdef EIGEN_VECTORS
00231             for (int k=1;k<=n;k++)
00232             {
00233               f=z[k][i+1];
00234               z[k][i+1]=s*z[k][i]+c*f;
00235               z[k][i]=c*z[k][i]-s*f;
00236             }
00237           #endif
00238         }
00239         d[l]=d[l]-p;
00240         e[l]=g;
00241         e[m]=0.0;
00242       }
00243     } while (m != l);
00244   }
00245 }
00246 
00258 dvar_vector get_eigen_values(const dvar_vector& _ddd,const dvar_vector& _eee)
00259 {
00260   ADUNCONST(dvar_vector,ddd)
00261   ADUNCONST(dvar_vector,eee)
00262 
00263   dvar_vector d(ddd.indexmin(),ddd.indexmax());
00264   dvar_vector e(eee.indexmin(),eee.indexmax());
00265 
00266   d=ddd;
00267   e=eee;
00268 
00269   int n=d.size();
00270   int m,l,iter,i;
00271   dvariable s,r,p,g,f,dd,c,b;
00272 
00273   for (i=2;i<=n;i++) e[i-1]=e[i];
00274   e[n]=0.0;
00275   for (l=1;l<=n;l++) {
00276     iter=0;
00277     do {
00278       for (m=l;m<=n-1;m++) {
00279         dd=fabs(d[m])+fabs(d[m+1]);
00280         if (fabs(e[m])+dd == dd) break;
00281       }
00282       if (m != l)
00283       {
00284         if (iter++ == 30)
00285         {
00286           cerr << "Maximum number of iterations exceeded in"
00287           " dvector eigen(const dmatrix& m)\n";
00288           ad_exit(1);
00289         }
00290         g=(d[l+1]-d[l])/(2.0*e[l]);
00291         r=sqrt((g*g)+1.0);
00292         g=d[m]-d[l]+e[l]/(g+SIGN(r,g));
00293         s=c=1.0;
00294         p=0.0;
00295         for (i=m-1;i>=l;i--) {
00296           f=s*e[i];
00297           b=c*e[i];
00298           if (fabs(f) >= fabs(g)) {
00299             c=g/f;
00300             r=sqrt((c*c)+1.0);
00301             e[i+1]=f*r;
00302             c *= (s=1.0/r);
00303           } else {
00304             s=f/g;
00305             r=sqrt((s*s)+1.0);
00306             e[i+1]=g*r;
00307             s *= (c=1.0/r);
00308           }
00309           g=d[i+1]-p;
00310           r=(d[i]-g)*s+2.0*c*b;
00311           p=s*r;
00312           d[i+1]=g+p;
00313           g=c*r-b;
00314           /* Next loop can be omitted if eigenvectors not wanted */
00315         }
00316         d[l]=d[l]-p;
00317         e[l]=g;
00318         e[m]=0.0;
00319       }
00320     } while (m != l);
00321   }
00322   return d;
00323 }
00324 #undef EIGEN_VECTORS