1 /*
2 boost header: numeric/odeint/gram_schmitt.hpp
3
4 Copyright 2011-2013 Karsten Ahnert
5 Copyright 2011 Mario Mulansky
6
7 Distributed under the Boost Software License, Version 1.0.
8 (See accompanying file LICENSE_1_0.txt or
9 copy at http://www.boost.org/LICENSE_1_0.txt)
10 */
11
12 #ifndef BOOST_NUMERIC_ODEINT_GRAM_SCHMITT_HPP_INCLUDED
13 #define BOOST_NUMERIC_ODEINT_GRAM_SCHMITT_HPP_INCLUDED
14
15 #include <boost/throw_exception.hpp>
16 #include <iterator>
17 #include <algorithm>
18 #include <numeric>
19
20 namespace boost {
21 namespace numeric {
22 namespace odeint {
23
24 template< class Iterator , class T >
normalize(Iterator first,Iterator last,T norm)25 void normalize( Iterator first , Iterator last , T norm )
26 {
27 while( first != last ) *first++ /= norm;
28 }
29
30 template< class Iterator , class T >
substract_vector(Iterator first1,Iterator last1,Iterator first2,T val)31 void substract_vector( Iterator first1 , Iterator last1 ,
32 Iterator first2 , T val )
33 {
34 while( first1 != last1 ) *first1++ -= val * ( *first2++ );
35 }
36
37 template< size_t num_of_lyap , class StateType , class LyapType >
gram_schmidt(StateType & x,LyapType & lyap,size_t n)38 void gram_schmidt( StateType &x , LyapType &lyap , size_t n )
39 {
40 if( !num_of_lyap ) return;
41 if( ptrdiff_t( ( num_of_lyap + 1 ) * n ) != std::distance( x.begin() , x.end() ) )
42 BOOST_THROW_EXCEPTION( std::domain_error( "renormalization() : size of state does not match the number of lyapunov exponents." ) );
43
44 typedef typename StateType::value_type value_type;
45 typedef typename StateType::iterator iterator;
46
47 value_type norm[num_of_lyap];
48 value_type tmp[num_of_lyap];
49 iterator first = x.begin() + n;
50 iterator beg1 = first , end1 = first + n ;
51
52 std::fill( norm , norm+num_of_lyap , 0.0 );
53
54 // normalize first vector
55 norm[0] = sqrt( std::inner_product( beg1 , end1 , beg1 , 0.0 ) );
56 normalize( beg1 , end1 , norm[0] );
57
58 beg1 += n;
59 end1 += n;
60
61 for( size_t j=1 ; j<num_of_lyap ; ++j , beg1+=n , end1+=n )
62 {
63 for( size_t k=0 ; k<j ; ++k )
64 {
65 tmp[k] = std::inner_product( beg1 , end1 , first + k*n , 0.0 );
66 // clog << j << " " << k << " " << tmp[k] << "\n";
67 }
68
69
70
71 for( size_t k=0 ; k<j ; ++k )
72 substract_vector( beg1 , end1 , first + k*n , tmp[k] );
73
74 // normalize j-th vector
75 norm[j] = sqrt( std::inner_product( beg1 , end1 , beg1 , 0.0 ) );
76 // clog << j << " " << norm[j] << "\n";
77 normalize( beg1 , end1 , norm[j] );
78 }
79
80 for( size_t j=0 ; j<num_of_lyap ; j++ )
81 lyap[j] += log( norm[j] );
82 }
83
84
85 } // namespace odeint
86 } // namespace numeric
87 } // namespace boost
88
89 #endif //BOOST_NUMERIC_ODEINT_GRAM_SCHMITT_HPP_INCLUDED
90