diff options
Diffstat (limited to 'src/boost/libs/numeric/odeint/examples/gram_schmidt.hpp')
-rw-r--r-- | src/boost/libs/numeric/odeint/examples/gram_schmidt.hpp | 89 |
1 files changed, 89 insertions, 0 deletions
diff --git a/src/boost/libs/numeric/odeint/examples/gram_schmidt.hpp b/src/boost/libs/numeric/odeint/examples/gram_schmidt.hpp new file mode 100644 index 000000000..f5f56808b --- /dev/null +++ b/src/boost/libs/numeric/odeint/examples/gram_schmidt.hpp @@ -0,0 +1,89 @@ +/* + boost header: numeric/odeint/gram_schmitt.hpp + + Copyright 2011-2013 Karsten Ahnert + Copyright 2011 Mario Mulansky + + Distributed under the Boost Software License, Version 1.0. + (See accompanying file LICENSE_1_0.txt or + copy at http://www.boost.org/LICENSE_1_0.txt) + */ + +#ifndef BOOST_NUMERIC_ODEINT_GRAM_SCHMITT_HPP_INCLUDED +#define BOOST_NUMERIC_ODEINT_GRAM_SCHMITT_HPP_INCLUDED + +#include <boost/throw_exception.hpp> +#include <iterator> +#include <algorithm> +#include <numeric> + +namespace boost { +namespace numeric { +namespace odeint { + +template< class Iterator , class T > +void normalize( Iterator first , Iterator last , T norm ) +{ + while( first != last ) *first++ /= norm; +} + +template< class Iterator , class T > +void substract_vector( Iterator first1 , Iterator last1 , + Iterator first2 , T val ) +{ + while( first1 != last1 ) *first1++ -= val * ( *first2++ ); +} + +template< size_t num_of_lyap , class StateType , class LyapType > +void gram_schmidt( StateType &x , LyapType &lyap , size_t n ) +{ + if( !num_of_lyap ) return; + if( ptrdiff_t( ( num_of_lyap + 1 ) * n ) != std::distance( x.begin() , x.end() ) ) + BOOST_THROW_EXCEPTION( std::domain_error( "renormalization() : size of state does not match the number of lyapunov exponents." ) ); + + typedef typename StateType::value_type value_type; + typedef typename StateType::iterator iterator; + + value_type norm[num_of_lyap]; + value_type tmp[num_of_lyap]; + iterator first = x.begin() + n; + iterator beg1 = first , end1 = first + n ; + + std::fill( norm , norm+num_of_lyap , 0.0 ); + + // normalize first vector + norm[0] = sqrt( std::inner_product( beg1 , end1 , beg1 , 0.0 ) ); + normalize( beg1 , end1 , norm[0] ); + + beg1 += n; + end1 += n; + + for( size_t j=1 ; j<num_of_lyap ; ++j , beg1+=n , end1+=n ) + { + for( size_t k=0 ; k<j ; ++k ) + { + tmp[k] = std::inner_product( beg1 , end1 , first + k*n , 0.0 ); + // clog << j << " " << k << " " << tmp[k] << "\n"; + } + + + + for( size_t k=0 ; k<j ; ++k ) + substract_vector( beg1 , end1 , first + k*n , tmp[k] ); + + // normalize j-th vector + norm[j] = sqrt( std::inner_product( beg1 , end1 , beg1 , 0.0 ) ); + // clog << j << " " << norm[j] << "\n"; + normalize( beg1 , end1 , norm[j] ); + } + + for( size_t j=0 ; j<num_of_lyap ; j++ ) + lyap[j] += log( norm[j] ); +} + + +} // namespace odeint +} // namespace numeric +} // namespace boost + +#endif //BOOST_NUMERIC_ODEINT_GRAM_SCHMITT_HPP_INCLUDED |