|
realroot_doc 0.1.1
|
|
realroot_doc 0.1.1
|
00001 /******************************************************************** 00002 * This file is part of the source code of the realroot library. 00003 * Author(s): J.P. Pavone, GALAAD, INRIA 00004 * $Id: householder.hpp,v 1.1 2005/07/11 10:03:55 jppavone Exp $ 00005 ********************************************************************/ 00006 #ifndef realroot_SOLVE_SBDSLV_LINEAR_HOUSEHOLDER_HPP 00007 #define realroot_SOLVE_SBDSLV_LINEAR_HOUSEHOLDER_HPP 00008 //-------------------------------------------------------------------- 00009 #include <math.h> 00010 #include <realroot/system_descartes1d.hpp> 00011 //-------------------------------------------------------------------- 00012 namespace mmx { 00013 //-------------------------------------------------------------------- 00014 namespace linear 00015 { 00016 /* 00017 * M: matrice ring_sparseetrique a reduire, a la fin contient la matrice de transformation 00018 * d: elements diagonaux 00019 * e: au dessus de la diagonale (e[0] = 0 ) 00020 */ 00021 template< typename real_t > 00022 void householder( real_t * M, real_t * d, real_t * e, int n ) 00023 { 00024 int i,j,k,l; 00025 real_t scale, hh, h, g, f; 00026 real_t * mi = M + (n-1)*n; 00027 for ( i = n-1; i >= 1; i--, mi -= n ) 00028 { 00029 l = i-1; 00030 h = scale = 0.0; 00031 if ( l > 0 ) 00032 { 00033 for ( k = 0; k <= l; scale += fabs(mi[k]), k ++ ) ; 00034 if ( scale == 0.0 ) e[i] = mi[l]; 00035 else 00036 { 00037 for ( k = 0; k <= l; mi[k] /= scale, h += mi[k]*mi[k], k ++ ) ; 00038 f = mi[l]; 00039 g = (f >= 0.0 ? -std::sqrt(h) : std::sqrt(h)); 00040 e[i] = scale*g; 00041 h -= f*g; 00042 mi[l] = f - g; 00043 f = 0.0; 00044 real_t * mj = M; 00045 for ( j = 0; j <= l; j ++, mj += n ) { 00046 mj[i] = mi[j]/h; 00047 g = 0.0; 00048 for ( k = 0; k <= j; g += mj[k]*mi[k], k ++ ) ; 00049 for ( k = j+1; k<= l; g += M[k*n+j]*mi[k], k ++ ) ; 00050 e[j] = g/h; 00051 f += e[j]*mi[j]; 00052 }; 00053 hh = f /(h+h); 00054 for ( j = 0; j <= l; j ++ ) 00055 { 00056 f = mi[j]; 00057 e[j] = g = e[j] - hh*f; 00058 for ( k = 0; k <= j; M[j*n+k] -= (f*e[k]+g*mi[k]), k ++ ) ; 00059 }; 00060 } 00061 } 00062 else 00063 e[i] = mi[l]; 00064 d[i] = h; 00065 }; 00066 00067 d[0] = 0; 00068 e[0] = 0; 00069 mi = M; 00070 for ( i = 0; i < n; i ++, mi += n ) 00071 { 00072 l = i - 1; 00073 if ( d[i] ) 00074 for ( j = 0; j <= l; j ++ ) 00075 { 00076 g = 0.0; 00077 for ( k = 0; k <= l; g += mi[k]*M[k*n+j], k ++ ) ; 00078 for ( k = 0; k <= l; M[k*n+j] -= g*M[k*n+i], k ++ ) ; 00079 }; 00080 d[i] = mi[i]; 00081 mi[i] = 1.0; 00082 for ( j = 0; j <= l; j ++ ) M[j*n+i] = mi[j] = 0.0; 00083 }; 00084 }; 00085 00086 template< typename real_t > 00087 real_t pythag( const real_t& a, const real_t& b) 00088 { 00089 real_t fa = std::abs(a); 00090 real_t fb = std::abs(b); 00091 real_t c; 00092 if ( fa > fb ) 00093 { 00094 c = fb/fa; 00095 return fa * std::sqrt(1.0+c*c); 00096 } 00097 else 00098 { 00099 c = fa/fb; 00100 return fb * std::sqrt(1.0+c*c); 00101 }; 00102 }; 00103 00104 00105 template< typename real_t > 00106 void tqli( real_t * d, real_t * e, int n, real_t * z ) 00107 { 00108 int m,l,iter,i,k; 00109 real_t s,r,p,g,f,dd,c,b; 00110 for ( i = 1; i < n; i ++ ) e[i-1] = e[i]; 00111 e[n-1] = 0.0; 00112 for ( l = 0; l < n; l++ )//for ( l = 1; l <= n; l ++ ) 00113 { 00114 00115 iter = 0; 00116 do 00117 { 00118 00119 for ( m = l; m < n-1; m++ )// for ( m = 1; m <=n-1; m++ ) 00120 { 00121 dd = std::abs(d[m])+std::abs(d[m+1]); 00122 if ( (real_t)(std::abs(e[m])+dd) == dd ) 00123 { break; }; 00124 }; 00125 if ( m != l ) 00126 { 00127 iter++; 00128 assert(iter<30); 00129 g = (d[l+1]-d[l])/(2.0*e[l]); 00130 r = pythag(g,(real_t)1.0); 00131 if ( g >= 0.0 ) g = d[m]-d[l]+e[l]/(g+std::abs(r)); 00132 else g = d[m]-d[l]+e[l]/(g-std::abs(r)); 00133 s = c = 1.0; 00134 p = 0.0; 00135 for ( i = m-1; i >= l; i-- ) 00136 { 00137 f = s*e[i]; 00138 b = c*e[i]; 00139 e[i+1] = (r=pythag(f,g)); 00140 if ( r == 0.0 ) 00141 { 00142 d[i+1] -= p; 00143 e[m] = 0.0; 00144 break; 00145 }; 00146 s = f/r; 00147 c = g/r; 00148 g = d[i+1]-p; 00149 r=(d[i]-g)*s+2.0*c*b; 00150 d[i+1]= g +(p=s*r); 00151 g = c*r-b; 00152 for ( k = 0; k < n; k++ )// for ( k = 1; k <= n; k++ ) 00153 { 00154 f = z[k*n+i+1];//f = z[k][i+1]; 00155 z[k*n+i+1] = s*z[k*n+i]+c*f;//z[k][i+1] = s*z[k][i]+c*f; 00156 z[k*n+i] = c*z[k*n+i]-s*f;//z[k][i] = c*z[k][i]-s*f; 00157 }; 00158 } 00159 if ( r == 0.0 && i >= l ) continue; 00160 d[l] -= p; 00161 e[l] = g; 00162 e[m] = 0.0; 00163 }; 00164 } 00165 while( m != l ); 00166 }; 00167 }; 00168 00169 template< class real_t > 00170 void symeigen( real_t * M, real_t * d, int n ) 00171 { 00172 real_t e[n]; 00173 householder( M, d, e, n ); 00174 tqli( d, e, n, M ); 00175 }; 00176 00177 template< typename real_t > 00178 void cpolynom( real_t * p, real_t * M, int n ) 00179 { 00180 real_t cp1[n+1]; 00181 real_t * p1 = cp1; 00182 real_t * p0 = p; 00183 int i,k,c; 00184 00185 for ( k = 0; k <= n; p0[k] = p1[k] = 0, k ++ ) ; 00186 00187 p0[0] = 1; 00188 p1[0] = M[n*(n-1)+(n-1)]; 00189 p1[1] = -1; 00190 00191 c = 0; 00192 00193 for ( i = n-2; i >= 0; i -- ) 00194 { 00195 real_t a = M[i*n+i]; 00196 real_t b = M[i*n+i+1]; 00197 00198 for ( k = 0; k <= n-2-i; k ++ ) p0[k] = -p0[k]*(b*b); 00199 for ( k = 0; k <= n-1-i; k ++ ) p0[k] += a*p1[k]; 00200 for ( k = 1; k <= n-i; k ++ ) 00201 p0[k] -= p1[k-1]; 00202 // for ( k = 0; k <= n-1-i; k ++ ) 00203 // p0[k] += a*p1[k]; 00204 // for ( k = 1; k <= n-i; k ++ ) 00205 // p0[k] -= p1[k-1]; 00206 std::swap(p0,p1); 00207 }; 00208 00209 if ( p1 != p ) 00210 std::copy( p1, p1 + n+1, p ); 00211 // if ( c ) std::copy( p1, p1 + n+1, p0 ); 00212 }; 00213 00214 template< typename real_t > 00215 void caract_polynom( real_t * p, real_t * M, int n ) 00216 { 00217 real_t d[n]; 00218 real_t e[n]; 00219 householder( M, d, e, n ); 00220 cpolynom( p, M, n ); 00221 }; 00222 00223 /* 00224 struct slv1dprm : std::vector< real_t > 00225 { 00226 void output( const real_t& a, const real_t& b ) 00227 { push_back( (a+b)/2.0 ); }; 00228 }; 00229 00230 template< typename real_t > 00231 void symeigen( real_t * vp, real_t * base, real_t * m, 00232 int n, const real_t& eps = 1e-12 ) 00233 { 00234 real_t tmp[n*n]; 00235 std::copy(m,m+n*n,tmp); 00236 real_t p[n+1]; 00237 caract_polynom(p, m, n ); 00238 typedef solvers::descartes_solver< real_t, solvers::bsearch_newton<real_t> > solver_t; 00239 solver_t tmp(n+1,eps); 00240 slv1dprm svp; 00241 tmp.solve( svp, p ); 00242 for ( unsigned i = 0; i < svp.size(); i ++ ) 00243 00244 }; 00245 */ 00246 }; 00247 //-------------------------------------------------------------------- 00248 } //namespace mmx 00249 /********************************************************************/ 00250 #endif //
