mmx::linear Namespace Reference

Functions

Function Documentation

void mmx::linear::caract_polynom ( real_t *  p,
real_t *  M,
int  n 
) [inline]

Definition at line 215 of file linear_householder.hpp.

References cpolynom(), and householder().

00216   {
00217     real_t d[n];
00218     real_t e[n];
00219     householder( M, d, e, n );
00220     cpolynom( p, M, n );
00221   };

void mmx::linear::cpolynom ( real_t *  p,
real_t *  M,
int  n 
) [inline]

Definition at line 178 of file linear_householder.hpp.

References mmx::sparse::copy(), and mmx::sparse::swap().

Referenced by caract_polynom().

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   };

bool mmx::linear::doolittle ( real_t const *const   A,
real_t *  L,
real_t *  U,
unsigned  n 
) [inline]

Definition at line 41 of file linear_doolittle.hpp.

References mmx::abs(), and mmx::vctops::sum().

Referenced by LUinverse(), and LUsolve().

00042   {
00043     unsigned k,i,p;
00044     for ( k = 0; k < n; k++ )
00045       {
00046         //      valueof(k);
00047         for ( i = k; i < n; i ++ )
00048           {
00049             U[ k*n + i ] = A[ k*n + i ];
00050             for ( p = 0; p < k; p ++ ) U[ k*n + i ] -= L[ k*n + p ]*U[ p*n + i ];
00051           };
00052         for ( i = k + 1; i < n; i ++  )
00053           {
00054             if ( std::abs( U[ k*n + k ] ) < numerics::epsilon<real_t>::result*n*n ) return false;
00055             real_t sum = 0.0;
00056             for ( p = 0; p < k; p ++ ) sum += L[i*n+p]*U[p*n+k];
00057             L[i*n+k] = (A[i*n + k] - sum)/U[ k*n + k ];
00058           };
00059       };
00060     return true;
00061   };

bool mmx::linear::gramschmidt ( real_t *  M,
unsigned  m,
unsigned  n 
) [inline]

Definition at line 17 of file linear_gramschmidt.hpp.

References assert.

00018   {
00019     unsigned i,j,k;
00020     assert(m<=n);
00021     for ( i = 0; i < m; i ++ )
00022       {
00023         real_t s = 0.0;
00024         for ( k = 0; k < n; k ++ )
00025           s += M[i*n+k]*M[i*n+k];
00026         if ( s < numerics::epsilon<real_t>::result ) return false;
00027         for ( j = i+1; j < m; j ++ )
00028           {
00029             real_t d = 0.0;
00030             for ( k = 0; k < n; k ++ )
00031               d += M[i*n+k]*M[j*n+k];
00032             d /= s;
00033             for ( k = 0; k < n; k ++ )
00034               M[j*n+k] -= d*M[i*n+k];
00035           };
00036       };
00037     return true;
00038   };

void mmx::linear::householder ( real_t *  M,
real_t *  d,
real_t *  e,
int  n 
) [inline]

Definition at line 22 of file linear_householder.hpp.

References mmx::scale(), and mmx::sqrt().

Referenced by caract_polynom(), and symeigen().

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   };

void mmx::linear::library ( real_t *  M,
unsigned  n,
real_t *  v 
) [inline]

Definition at line 45 of file linear_gramschmidt.hpp.

References mmx::sparse::copy().

00046   {
00047     int goods[n];
00048     int c = 0;
00049     int i,j,k,p;
00050     real_t sc[ n ];
00051     
00052     c = 0;
00053     for ( i = 0; i < n; i ++ )
00054       {
00055         real_t s = 0.0;
00056         for ( k = 0; k < n; s += M[i*n+k]*M[i*n+k], k ++ ) ;
00057         if ( s < numerics::epsilon<real_t>::result )   continue; 
00058         goods[c++] = i;
00059         for ( j = i+1; j < n; j ++ )
00060           {
00061             real_t d = 0.0;
00062             for ( k = 0; k < n; k ++ )
00063               d += M[i*n+k]*M[j*n+k];
00064             d /= s;
00065             for ( k = 0; k < n; k ++ )
00066               M[j*n+k] -= d*M[i*n+k];
00067           };
00068         sc[i] = s;
00069       };
00070     if ( c == n ) return;
00071     k = 0;
00072     while( k < c && goods[k] == k ) k ++;
00073     while ( k < c )
00074       {
00075         std::copy( M + goods[k]*n, M + (goods[k]+1)*n, 
00076                    M + k );
00077         sc[k] = sc[goods[k]];
00078         k++;
00079       };
00080     //    vctops::print(goods,c);
00081     //    std::cout << std::endl;
00082     //    vctops::print(M,n);
00083     //    std::cout << std::endl;
00084     //    vctops::print(M+n,n);
00085     //    std::cout << std::endl;
00086     //    std::cout << "c = " << c << std::endl;
00087     //    std::cout << "n = " << n << std::endl;
00088     for ( i = c; i < n; i ++ )
00089       for ( k = 0; k < n; k ++ )
00090         M[i*n+k] = ((real_t)rand())/((real_t)RAND_MAX); 
00091     //    vctops::print(M+c*n,n);
00092     //    std::cout << std::endl;
00093     for ( i = c; i  < n; i ++ )
00094       {
00095         for ( p = i-1; p >= 0; p-- )
00096           {
00097             //      std::cout << sc[p] << std::endl;
00098             real_t d = 0.0;
00099             for ( k = 0; k < n; k ++ )
00100               d += M[p*n+k]*M[i*n+k];
00101             d /= sc[p];
00102             for ( k = 0; k < n; k ++ )
00103               {
00104                 //              std::cout << M[i*n+k] << std::endl;
00105                 M[i*n+k] -= d*M[p*n+k];
00106                 //              std::cout << M[i*n+k] << std::endl;
00107               };
00108           };
00109         real_t s = 0.0;
00110         for ( k = 0; k < n; k ++ )
00111           s += M[i*n+k]*M[i*n+k];
00112         sc[i] = s;
00113       };
00114   };

void mmx::linear::Lsolve ( real_t *  x,
real_t const *const   L,
real_t const *const   b,
unsigned  n 
) [inline]

Definition at line 17 of file linear_doolittle.hpp.

References mmx::vctops::sum().

Referenced by LUinverse(), and LUsolve().

00018   {
00019     for ( unsigned i = 0; i  < n ; i ++ )
00020       {
00021         real_t sum = 0.0;
00022         for ( unsigned j = 0; j < i; j ++ ) sum += L[i*n+j]*x[j];
00023         x[i] = (b[i]-sum);
00024       };
00025   };

bool mmx::linear::LUinverse ( real_t *  iM,
real_t const *const   M,
unsigned  n 
) [inline]

Definition at line 74 of file linear_doolittle.hpp.

References doolittle(), mmx::vct::fill(), Lsolve(), and Usolve().

Referenced by mmx::realroot::include3().

00075   {
00076     real_t lu[n*n];
00077     real_t xtmp[n] ;
00078     real_t tmp[n];
00079     real_t _b[n+1];
00080     real_t * b = _b+1;
00081     std::fill( b, b + n , (real_t)0.0 );
00082     if ( ! doolittle( M, lu, lu, n ) ) return false;
00083     for ( unsigned c = 0; c < n; c ++ )
00084       {
00085         if (c > 0) b[c-1] = 0;
00086         b[c] = 1;
00087         Lsolve( xtmp, lu, b, n );
00088         if ( !Usolve( tmp, lu, xtmp, n ) ) return false;
00089         for ( unsigned i = 0; i < n; i ++ ) iM[ i*n + c ] = tmp[i]; 
00090       };
00091     return true;
00092   };

bool mmx::linear::LUsolve ( real_t *  x,
real_t const *const   M,
real_t const *const   b,
unsigned  n 
) [inline]

Definition at line 64 of file linear_doolittle.hpp.

References doolittle(), Lsolve(), and Usolve().

00065   {
00066     real_t tmp[ n*n ];
00067     real_t xtmp[n] ;
00068     if ( !doolittle( M, tmp, tmp, n ) ) return false;
00069     Lsolve( xtmp, tmp, b, n );
00070     return Usolve( x, tmp, xtmp, n );
00071   };

real_t mmx::linear::pythag ( const real_t &  a,
const real_t &  b 
) [inline]

Definition at line 87 of file linear_householder.hpp.

References mmx::abs(), and mmx::sqrt().

Referenced by tqli().

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    };

void mmx::linear::symeigen ( real_t *  M,
real_t *  d,
int  n 
) [inline]

Definition at line 170 of file linear_householder.hpp.

References householder(), and tqli().

00171   {
00172     real_t e[n];
00173     householder( M, d, e, n );
00174     tqli( d, e, n, M );
00175   };

void mmx::linear::tqli ( real_t *  d,
real_t *  e,
int  n,
real_t *  z 
) [inline]

Definition at line 106 of file linear_householder.hpp.

References mmx::abs(), assert, and pythag().

Referenced by symeigen().

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   };

bool mmx::linear::Usolve ( real_t *  x,
real_t const *const   U,
real_t const *const   b,
unsigned  n 
) [inline]

Definition at line 28 of file linear_doolittle.hpp.

References mmx::abs(), and mmx::vctops::sum().

Referenced by LUinverse(), and LUsolve().

00029   {
00030     for ( int  i = n-1; i  >=0  ; i -- )
00031       {
00032         real_t sum = 0.0;
00033         for ( int j = n-1; j > i; j -- ) sum += U[i*n+j]*x[j];
00034         if ( std::abs(U[i*n+i]) < numerics::epsilon<real_t>::result*n*n) return false;
00035         x[i] = (b[i]-sum)/U[i*n+i];
00036       };
00037     return true;
00038   };

Generated on 6 Dec 2012 for realroot by  doxygen 1.6.1