1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889 |
- /*
- 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
|