00001 /*************************************************************************** 00002 eigenvalues.h - description 00003 ------------------- 00004 begin : Tue May 16 2000 00005 copyright : (C) 2000 by Michael Peeters 00006 email : Michael.Peeters@vub.ac.be 00007 ***************************************************************************/ 00008 00009 /*************************************************************************** 00010 * * 00011 * This program is free software; you can redistribute it and/or modify * 00012 * it under the terms of the GNU General Public License as published by * 00013 * the Free Software Foundation; either version 2 of the License, or * 00014 * (at your option) any later version. * 00015 * * 00016 ***************************************************************************/ 00017 00018 #ifndef EIGENVALUES_H 00019 #define EIGENVALUES_H 00020 00021 #include "numerictypes.h" 00022 #include "numerictraits.h" 00023 #include "invariant.h" 00024 #include "utility.h" 00025 #include <stdexcept> 00026 00027 namespace MODEL 00028 { 00029 00033 template<integer dims, typename nelem=number, class NT = NumericTraits<nelem,dims> > 00034 class Eigenvalues 00035 { 00036 public: 00037 00038 typedef typename NT::number numT; 00039 typedef typename NT::vect vect; 00040 typedef typename NT::matrix matrix; 00041 typedef typename NT::vf vf; 00042 00043 typedef typename NumericTraits<complex,dims>::vect EVect; 00044 00046 Eigenvalues(const matrix& ma) : m(ma) 00047 {calculate();} 00048 00050 const vect& real(void) {return wr;} 00051 00053 const vect& imag(void) {return wi;} 00054 00055 private: 00059 const matrix& balance(void); 00060 00063 const matrix& hessenberg(void); 00064 00066 void calculate(void); 00067 00069 NO_COPY(Eigenvalues); 00070 00071 private: 00072 matrix m; 00073 00074 vect wr; // real part 00075 vect wi; // imaginary part 00076 00077 static const number radix=2.; 00078 }; 00079 00082 template <integer dims, typename nelem, class NT > 00083 void 00084 Eigenvalues<dims,nelem,NT>::calculate(void) 00085 { 00086 // Make sure we have balanced it and reduced to upper hessenberg 00087 balance(); 00088 hessenberg(); 00089 00090 numT anorm(abs(m[0][0])); 00091 00092 for (integer i=0;i<dims;i++) // compute matrix norm 00093 for (integer j=max(i-1,0);j<dims;j++) // Because we start from upper Hessenberg ! 00094 anorm += abs(m[i][j]); 00095 00096 integer nn(dims-1); 00097 numT t(0.); // Gets changed only by an exceptional shift. 00098 00099 while (nn >= 0) // Search next eigenvalue 00100 { 00101 integer its(0); // iterations 00102 integer l(0); 00103 do 00104 { 00105 for (l=nn;l>=1;l--) // look for single small subdiagonal 00106 { 00107 numT s=abs(m[l-1][l-1])+abs(m[l][l]); 00108 if (s == 0.0) s=anorm; 00109 if ( ( abs(m[l][l-1]) + s ) == s) break; 00110 } 00111 numT x(m[nn][nn]),y(0.),z(0.); 00112 if (l == nn) // One root found 00113 { 00114 wr[nn]=x+t; // Compensate for shift 00115 wi[nn]=0.0; 00116 nn--; 00117 } 00118 else 00119 { 00120 y = m[nn-1][nn-1]; 00121 numT w(m[nn][nn-1]*m[nn-1][nn]); 00122 if (l == (nn-1)) 00123 { // Two roots found 00124 numT p = 0.5*(y-x); 00125 numT q=p*p+w; 00126 z=sqrt(abs(q)); 00127 x += t; // Compensate for shift 00128 if (q >= 0.0) 00129 { // a real pair 00130 z= p + ((p>=0.)?abs(z):-abs(z)); 00131 wr[nn-1]=wr[nn]=x+z; 00132 if (z) wr[nn]=x-w/z; 00133 wi[nn-1]=wi[nn]=0.0; 00134 } 00135 else 00136 { // a complex pair 00137 wr[nn-1]=wr[nn]=x+p; 00138 wi[nn-1]= -(wi[nn]=z); 00139 } 00140 nn -= 2; 00141 } 00142 else 00143 { // no roots found, continue 00144 numT p(0.),q(0.),r(0.),s(0.); 00145 00146 if (its == 30) throw std::logic_error("Too many iterations in eigenvalue loop"); 00147 if (its == 10 || its == 20) 00148 { // do a exceptional shift 00149 t += x; 00150 for (integer i=0;i<=nn;i++) m[i][i] -= x; 00151 s = abs(m[nn][nn-1]) + abs(m[nn-1][nn-2]); 00152 y=x=0.75*s; 00153 w = -0.4375*s*s; 00154 } 00155 00156 ++its; 00157 00158 integer mm(0); 00159 00160 for (mm=(nn-2);mm>=l;mm--) 00161 { // do shift and look for 2 consequtive small elements 00162 z=m[mm][mm]; 00163 r=x-z; 00164 s=y-z; 00165 00166 p=(r*s-w)/m[mm+1][mm]+m[mm][mm+1]; 00167 q=m[mm+1][mm+1]-z-r-s; 00168 00169 r=m[mm+2][mm+1]; 00170 s=abs(p)+abs(q)+abs(r); 00171 00172 p /= s; 00173 q /= s; 00174 r /= s; 00175 00176 if (mm == l) break; 00177 numT u= abs(m[mm][mm-1])*(abs(q)+abs(r)); 00178 numT v= abs(p)*(abs(m[mm-1][mm-1])+abs(z)+abs(m[mm+1][mm+1])); 00179 if ((u+v) == v) break; 00180 } 00181 00182 for (integer i=mm+2;i<=nn;i++) 00183 { 00184 m[i][i-2]=0.0; 00185 if (i != (mm+2)) m[i][i-3]=0.0; 00186 } 00187 00188 for (integer k=mm;k<=nn-1;k++) 00189 { // Double QR step 00190 if (k != mm) 00191 { 00192 p=m[k][k-1]; // Householder vector setup 00193 q=m[k+1][k-1]; 00194 r=0.; 00195 if (k != (nn-1)) r=m[k+2][k-1]; 00196 x=abs(p)+abs(q)+abs(r); 00197 if (x != 0.0) 00198 { 00199 p /= x; // Scale to prevent overflow 00200 q /= x; 00201 r /= x; 00202 } 00203 } 00204 s=sqrt(p*p+q*q+r*r); 00205 s=((p>=0.)?s:-s); 00206 if (s != 0.0) 00207 { 00208 if (k == mm) 00209 { 00210 if (l != mm) m[k][k-1] = -m[k][k-1]; 00211 } 00212 else m[k][k-1] = -s*x; 00213 p += s; 00214 00215 x= p/s; 00216 y= q/s; 00217 z= r/s; 00218 00219 q /= p; 00220 r /= p; 00221 00222 for (integer j=k;j<=nn;j++) 00223 { 00224 p=m[k][j]+q*m[k+1][j]; 00225 if (k != (nn-1)) 00226 { 00227 p += r*m[k+2][j]; 00228 m[k+2][j] -= p*z; 00229 } 00230 m[k+1][j] -= p*y; 00231 m[k][j] -= p*x; 00232 } 00233 integer mmin = nn<k+3 ? nn : k+3; 00234 for (integer i=l;i<=mmin;i++) 00235 { 00236 p=x*m[i][k]+y*m[i][k+1]; 00237 if (k != (nn-1)) 00238 { 00239 p += z*m[i][k+2]; 00240 m[i][k+2] -= p*r; 00241 } 00242 m[i][k+1] -= p*q; 00243 m[i][k] -= p; 00244 } 00245 } 00246 } 00247 } 00248 } 00249 } while (l < nn); 00250 } 00251 } 00252 00253 template <integer dims, typename nelem, class NT > 00254 const Eigenvalues<dims,nelem,NT>::matrix& 00255 Eigenvalues<dims,nelem,NT>::balance(void) 00256 { 00257 integer last(0); 00258 numT sqrdx(radix*radix); 00259 00260 while (last == 0) 00261 { 00262 last=1; 00263 for (integer i=0;i<dims;i++) 00264 { 00265 // Calculate row and column norms. 00266 numT r(0.),c(0.),g(0.); 00267 for (integer j=0;j<dims;j++) 00268 if (j != i) 00269 { 00270 c += abs(m[j][i]); 00271 r += abs(m[i][j]); 00272 } 00273 00274 if ((c!=0.) && (r!=0.)) 00275 { 00276 // If both are nonzero, 00277 numT f(1.); 00278 numT s(c+r); 00279 00280 g=r/radix; 00281 while (c<g) 00282 { 00283 // nd the integer power of the machine radix that comes closest to balancing the matrix. 00284 f *= radix; 00285 c *= sqrdx; 00286 } 00287 00288 g=r*radix; 00289 while (c>g) 00290 { 00291 f /= radix; 00292 c /= sqrdx; 00293 } 00294 if (((c+r)/f) < (0.95*s)) 00295 { 00296 last=0; 00297 g=1.0/f; 00298 // numT a=m[0][0]; 00299 for (integer k=0;k<dims;k++) m[i][k] *= g; // Apply similarity transforma- tion. 00300 // numT b=m[0][0]; 00301 for (integer k=0;k<dims;k++) m[k][i] *= f; 00302 } 00303 } 00304 } 00305 } 00306 return m; 00307 } 00308 00310 template<typename T> 00311 void swap(T& a,T& b) 00312 { 00313 T temp(a); 00314 a=b; 00315 b=temp; 00316 } 00317 00318 template <integer dims, typename nelem, class NT > 00319 const Eigenvalues<dims,nelem,NT>::matrix& 00320 Eigenvalues<dims,nelem,NT>::hessenberg(void) 00321 { 00322 integer i(0); 00323 for (integer rp1=1;rp1<(dims-1);rp1++) 00324 { //rp1 is r + 1in NRC. 00325 numT x(0.); 00326 i=rp1; 00327 00328 for (integer j=rp1;j<dims;j++) 00329 { // Find the pivot. 00330 if (abs(m[j][rp1-1]) > abs(x)) 00331 { 00332 x=m[j][rp1-1]; 00333 i=j; 00334 } 00335 } 00336 if (i != rp1) 00337 { // Interchange rows and columns. 00338 for (integer j=rp1-1;j<dims;j++) swap(m[i][j],m[rp1][j]); 00339 for (integer j=0;j<dims;j++) swap(m[j][i],m[j][rp1]); 00340 } 00341 if (x!=0.) 00342 { 00343 //Carry out the elimination. 00344 for (integer k=rp1+1;k<dims;k++) 00345 { 00346 numT y(0.); 00347 if ((y=m[k][rp1-1]) != 0.0) 00348 { 00349 y /= x; 00350 m[k][rp1-1]=y; 00351 for (integer j=rp1;j<dims;j++) m[k][j] -= y*m[rp1][j]; 00352 for (integer j=0;j<dims;j++) m[j][rp1] += y*m[j][k]; 00353 } 00354 } 00355 } 00356 } 00357 return m; 00358 } 00359 00360 } // end namespace 00361 #endif 00362 00363 /********************************************************************* 00364 $Id: eigenvalues.h,v 1.1 2001/05/22 10:54:55 mpeeters Exp $ 00365 ********************************************************************** 00366 00367 $Log: eigenvalues.h,v $ 00368 Revision 1.1 2001/05/22 10:54:55 mpeeters 00369 Moved sources and headers for libModel to model/ subdirectory, in an attempt to rationalize the source tree. This should make things "netter". 00370 00371 Revision 1.2 2000/09/15 10:26:31 mpeeters 00372 Cleaned out header and added CVS tails to files, removed superfuous 00373 @author comments, inserted dates 00374 00375 *********************************************************************/
More Info? Michael Peeters. Also, check our research website: www.alna.vub.ac.be
Last update: June 2002.
Looking for Open Source? Check out
!