|
realroot_doc 0.1.1
|
|
realroot_doc 0.1.1
|
| void mmx::linear::caract_polynom | ( | real_t * | p, |
| real_t * | M, | ||
| int | n | ||
| ) |
Definition at line 215 of file linear_householder.hpp.
References cpolynom(), and householder().
{
real_t d[n];
real_t e[n];
householder( M, d, e, n );
cpolynom( p, M, n );
};
| void mmx::linear::cpolynom | ( | real_t * | p, |
| real_t * | M, | ||
| int | n | ||
| ) |
Definition at line 178 of file linear_householder.hpp.
References mmx::sparse::copy(), and mmx::sparse::swap().
Referenced by caract_polynom().
{
real_t cp1[n+1];
real_t * p1 = cp1;
real_t * p0 = p;
int i,k,c;
for ( k = 0; k <= n; p0[k] = p1[k] = 0, k ++ ) ;
p0[0] = 1;
p1[0] = M[n*(n-1)+(n-1)];
p1[1] = -1;
c = 0;
for ( i = n-2; i >= 0; i -- )
{
real_t a = M[i*n+i];
real_t b = M[i*n+i+1];
for ( k = 0; k <= n-2-i; k ++ ) p0[k] = -p0[k]*(b*b);
for ( k = 0; k <= n-1-i; k ++ ) p0[k] += a*p1[k];
for ( k = 1; k <= n-i; k ++ )
p0[k] -= p1[k-1];
// for ( k = 0; k <= n-1-i; k ++ )
// p0[k] += a*p1[k];
// for ( k = 1; k <= n-i; k ++ )
// p0[k] -= p1[k-1];
std::swap(p0,p1);
};
if ( p1 != p )
std::copy( p1, p1 + n+1, p );
// if ( c ) std::copy( p1, p1 + n+1, p0 );
};
| bool mmx::linear::doolittle | ( | real_t const *const | A, |
| real_t * | L, | ||
| real_t * | U, | ||
| unsigned | n | ||
| ) |
Definition at line 41 of file linear_doolittle.hpp.
References mmx::abs(), and mmx::vctops::sum().
Referenced by LUinverse(), and LUsolve().
{
unsigned k,i,p;
for ( k = 0; k < n; k++ )
{
// valueof(k);
for ( i = k; i < n; i ++ )
{
U[ k*n + i ] = A[ k*n + i ];
for ( p = 0; p < k; p ++ ) U[ k*n + i ] -= L[ k*n + p ]*U[ p*n + i ];
};
for ( i = k + 1; i < n; i ++ )
{
if ( std::abs( U[ k*n + k ] ) < numerics::epsilon<real_t>::result*n*n ) return false;
real_t sum = 0.0;
for ( p = 0; p < k; p ++ ) sum += L[i*n+p]*U[p*n+k];
L[i*n+k] = (A[i*n + k] - sum)/U[ k*n + k ];
};
};
return true;
};
| bool mmx::linear::gramschmidt | ( | real_t * | M, |
| unsigned | m, | ||
| unsigned | n | ||
| ) |
Definition at line 17 of file linear_gramschmidt.hpp.
References assert.
{
unsigned i,j,k;
assert(m<=n);
for ( i = 0; i < m; i ++ )
{
real_t s = 0.0;
for ( k = 0; k < n; k ++ )
s += M[i*n+k]*M[i*n+k];
if ( s < numerics::epsilon<real_t>::result ) return false;
for ( j = i+1; j < m; j ++ )
{
real_t d = 0.0;
for ( k = 0; k < n; k ++ )
d += M[i*n+k]*M[j*n+k];
d /= s;
for ( k = 0; k < n; k ++ )
M[j*n+k] -= d*M[i*n+k];
};
};
return true;
};
| void mmx::linear::householder | ( | real_t * | M, |
| real_t * | d, | ||
| real_t * | e, | ||
| int | n | ||
| ) |
Definition at line 22 of file linear_householder.hpp.
References mmx::scale(), and mmx::sqrt().
Referenced by caract_polynom(), and symeigen().
{
int i,j,k,l;
real_t scale, hh, h, g, f;
real_t * mi = M + (n-1)*n;
for ( i = n-1; i >= 1; i--, mi -= n )
{
l = i-1;
h = scale = 0.0;
if ( l > 0 )
{
for ( k = 0; k <= l; scale += fabs(mi[k]), k ++ ) ;
if ( scale == 0.0 ) e[i] = mi[l];
else
{
for ( k = 0; k <= l; mi[k] /= scale, h += mi[k]*mi[k], k ++ ) ;
f = mi[l];
g = (f >= 0.0 ? -std::sqrt(h) : std::sqrt(h));
e[i] = scale*g;
h -= f*g;
mi[l] = f - g;
f = 0.0;
real_t * mj = M;
for ( j = 0; j <= l; j ++, mj += n ) {
mj[i] = mi[j]/h;
g = 0.0;
for ( k = 0; k <= j; g += mj[k]*mi[k], k ++ ) ;
for ( k = j+1; k<= l; g += M[k*n+j]*mi[k], k ++ ) ;
e[j] = g/h;
f += e[j]*mi[j];
};
hh = f /(h+h);
for ( j = 0; j <= l; j ++ )
{
f = mi[j];
e[j] = g = e[j] - hh*f;
for ( k = 0; k <= j; M[j*n+k] -= (f*e[k]+g*mi[k]), k ++ ) ;
};
}
}
else
e[i] = mi[l];
d[i] = h;
};
d[0] = 0;
e[0] = 0;
mi = M;
for ( i = 0; i < n; i ++, mi += n )
{
l = i - 1;
if ( d[i] )
for ( j = 0; j <= l; j ++ )
{
g = 0.0;
for ( k = 0; k <= l; g += mi[k]*M[k*n+j], k ++ ) ;
for ( k = 0; k <= l; M[k*n+j] -= g*M[k*n+i], k ++ ) ;
};
d[i] = mi[i];
mi[i] = 1.0;
for ( j = 0; j <= l; j ++ ) M[j*n+i] = mi[j] = 0.0;
};
};
| void mmx::linear::library | ( | real_t * | M, |
| unsigned | n, | ||
| real_t * | v | ||
| ) |
Definition at line 45 of file linear_gramschmidt.hpp.
References mmx::sparse::copy().
{
int m = n;
int goods[n];
int c = 0;
int i,j,k,p;
int already =0;
real_t sc[ n ];
c = 0;
for ( i = 0; i < n; i ++ )
{
real_t s = 0.0;
for ( k = 0; k < n; s += M[i*n+k]*M[i*n+k], k ++ ) ;
if ( s < numerics::epsilon<real_t>::result ) continue;
goods[c++] = i;
for ( j = i+1; j < n; j ++ )
{
real_t d = 0.0;
for ( k = 0; k < n; k ++ )
d += M[i*n+k]*M[j*n+k];
d /= s;
for ( k = 0; k < n; k ++ )
M[j*n+k] -= d*M[i*n+k];
};
sc[i] = s;
};
if ( c == n ) return;
k = 0;
while( k < c && goods[k] == k ) k ++;
while ( k < c )
{
std::copy( M + goods[k]*n, M + (goods[k]+1)*n,
M + k );
sc[k] = sc[goods[k]];
k++;
};
// vctops::print(goods,c);
// std::cout << std::endl;
// vctops::print(M,n);
// std::cout << std::endl;
// vctops::print(M+n,n);
// std::cout << std::endl;
// std::cout << "c = " << c << std::endl;
// std::cout << "n = " << n << std::endl;
for ( i = c; i < n; i ++ )
for ( k = 0; k < n; k ++ )
M[i*n+k] = ((real_t)rand())/((real_t)RAND_MAX);
// vctops::print(M+c*n,n);
// std::cout << std::endl;
for ( i = c; i < n; i ++ )
{
for ( p = i-1; p >= 0; p-- )
{
// std::cout << sc[p] << std::endl;
real_t d = 0.0;
for ( k = 0; k < n; k ++ )
d += M[p*n+k]*M[i*n+k];
d /= sc[p];
for ( k = 0; k < n; k ++ )
{
// std::cout << M[i*n+k] << std::endl;
M[i*n+k] -= d*M[p*n+k];
// std::cout << M[i*n+k] << std::endl;
};
};
real_t s = 0.0;
for ( k = 0; k < n; k ++ )
s += M[i*n+k]*M[i*n+k];
sc[i] = s;
};
};
| void mmx::linear::Lsolve | ( | real_t * | x, |
| real_t const *const | L, | ||
| real_t const *const | b, | ||
| unsigned | n | ||
| ) |
Definition at line 17 of file linear_doolittle.hpp.
References mmx::vctops::sum().
Referenced by LUinverse(), and LUsolve().
| bool mmx::linear::LUinverse | ( | real_t * | iM, |
| real_t const *const | M, | ||
| unsigned | n | ||
| ) |
Definition at line 74 of file linear_doolittle.hpp.
References doolittle(), mmx::vct::fill(), Lsolve(), and Usolve().
Referenced by mmx::realroot::include3().
{
real_t lu[n*n];
real_t xtmp[n] ;
real_t tmp[n];
real_t _b[n+1];
real_t * b = _b+1;
std::fill( b, b + n , (real_t)0.0 );
if ( ! doolittle( M, lu, lu, n ) ) return false;
for ( unsigned c = 0; c < n; c ++ )
{
if (c > 0) b[c-1] = 0;
b[c] = 1;
Lsolve( xtmp, lu, b, n );
if ( !Usolve( tmp, lu, xtmp, n ) ) return false;
for ( unsigned i = 0; i < n; i ++ ) iM[ i*n + c ] = tmp[i];
};
return true;
};
| bool mmx::linear::LUsolve | ( | real_t * | x, |
| real_t const *const | M, | ||
| real_t const *const | b, | ||
| unsigned | n | ||
| ) |
Definition at line 64 of file linear_doolittle.hpp.
References doolittle(), Lsolve(), and Usolve().
| real_t mmx::linear::pythag | ( | const real_t & | a, |
| const real_t & | b | ||
| ) |
Definition at line 87 of file linear_householder.hpp.
References mmx::abs(), and mmx::sqrt().
Referenced by tqli().
| void mmx::linear::symeigen | ( | real_t * | M, |
| real_t * | d, | ||
| int | n | ||
| ) |
Definition at line 170 of file linear_householder.hpp.
References householder(), and tqli().
{
real_t e[n];
householder( M, d, e, n );
tqli( d, e, n, M );
};
| void mmx::linear::tqli | ( | real_t * | d, |
| real_t * | e, | ||
| int | n, | ||
| real_t * | z | ||
| ) |
Definition at line 106 of file linear_householder.hpp.
References mmx::abs(), assert, and pythag().
Referenced by symeigen().
{
int m,l,iter,i,k;
real_t s,r,p,g,f,dd,c,b;
for ( i = 1; i < n; i ++ ) e[i-1] = e[i];
e[n-1] = 0.0;
for ( l = 0; l < n; l++ )//for ( l = 1; l <= n; l ++ )
{
iter = 0;
do
{
for ( m = l; m < n-1; m++ )// for ( m = 1; m <=n-1; m++ )
{
dd = std::abs(d[m])+std::abs(d[m+1]);
if ( (real_t)(std::abs(e[m])+dd) == dd )
{ break; };
};
if ( m != l )
{
iter++;
assert(iter<30);
g = (d[l+1]-d[l])/(2.0*e[l]);
r = pythag(g,(real_t)1.0);
if ( g >= 0.0 ) g = d[m]-d[l]+e[l]/(g+std::abs(r));
else g = d[m]-d[l]+e[l]/(g-std::abs(r));
s = c = 1.0;
p = 0.0;
for ( i = m-1; i >= l; i-- )
{
f = s*e[i];
b = c*e[i];
e[i+1] = (r=pythag(f,g));
if ( r == 0.0 )
{
d[i+1] -= p;
e[m] = 0.0;
break;
};
s = f/r;
c = g/r;
g = d[i+1]-p;
r=(d[i]-g)*s+2.0*c*b;
d[i+1]= g +(p=s*r);
g = c*r-b;
for ( k = 0; k < n; k++ )// for ( k = 1; k <= n; k++ )
{
f = z[k*n+i+1];//f = z[k][i+1];
z[k*n+i+1] = s*z[k*n+i]+c*f;//z[k][i+1] = s*z[k][i]+c*f;
z[k*n+i] = c*z[k*n+i]-s*f;//z[k][i] = c*z[k][i]-s*f;
};
}
if ( r == 0.0 && i >= l ) continue;
d[l] -= p;
e[l] = g;
e[m] = 0.0;
};
}
while( m != l );
};
};
| bool mmx::linear::Usolve | ( | real_t * | x, |
| real_t const *const | U, | ||
| real_t const *const | b, | ||
| unsigned | n | ||
| ) |
Definition at line 28 of file linear_doolittle.hpp.
References mmx::abs(), and mmx::vctops::sum().
Referenced by LUinverse(), and LUsolve().
