29 #ifndef __FEELPP_GLAS_HPP
30 #define __FEELPP_GLAS_HPP 1
36 #include <boost/random/linear_congruential.hpp>
37 #include <boost/random/uniform_int.hpp>
38 #include <boost/random/uniform_real.hpp>
39 #include <boost/random/variate_generator.hpp>
41 #include <boost/lambda/lambda.hpp>
42 #include <boost/lambda/algorithm.hpp>
43 #include <boost/lambda/bind.hpp>
44 #include <boost/lambda/if.hpp>
46 #include <boost/numeric/ublas/vector.hpp>
47 #include <boost/numeric/ublas/matrix.hpp>
48 #include <boost/numeric/ublas/matrix_sparse.hpp>
49 #include <boost/numeric/ublas/symmetric.hpp>
50 #include <boost/numeric/ublas/operation.hpp>
51 #include <boost/numeric/ublas/matrix_proxy.hpp>
52 #include <boost/numeric/ublas/io.hpp>
54 #include <feel/feelcore/feel.hpp>
56 #include <feel/feelcore/traits.hpp>
60 const double Pi = 3.14159265358979323846264338328;
61 const double TGV = 1e20;
64 inline T Min (
const T &a,
const T &b )
69 inline T Max (
const T &a,
const T & b )
74 inline T Abs (
const T &a )
76 return a < 0 ? -a : a;
80 inline void Exchange ( T& a, T& b )
87 inline T Max (
const T &a,
const T & b,
const T & c )
89 return Max( Max( a, b ), c );
92 inline T Min (
const T &a,
const T & b,
const T & c )
94 return Min( Min( a, b ), c );
97 namespace ublas = boost::numeric::ublas;
105 using ublas::norm_inf;
106 using ublas::index_norm_inf;
109 using ublas::outer_prod;
110 using ublas::inner_prod;
112 struct norm_inf_adaptor
115 Real operator()( E
const& __v )
const
117 return norm_inf( __v );
124 template <
typename T =
double, u
int16_type S = 3>
128 typedef ublas::vector<T> type;
135 typedef node<double, 3>::type node_type;
140 typedef node<double, 3>::type gradient_node_type;
145 typedef ublas::symmetric_matrix<double, ublas::lower, ublas::row_major, ublas::bounded_array<double, 9> >
hessian_node_type;
148 typedef ublas::matrix<double, ublas::column_major, ublas::bounded_array<double, 9> > lapack_matrix_type;
149 typedef ublas::symmetric_adaptor<lapack_matrix_type, ublas::lower> symmetric_matrix_type;
151 typedef lapack_matrix_type transformation_matrix_type;
156 template <
typename T =
double, u
int16_type S = 256>
160 typedef ublas::matrix<T, ublas::column_major> type;
164 typedef matrix_node<double>::type matrix_node_type;
169 operator<<( DebugStream& __os, node<real64_type>::type
const& __n )
171 if ( __os.doPrint() )
173 std::ostringstream __str;
177 __os << __str.str() <<
"\n";
184 operator<<( NdebugStream& os, node<real64_type>::type
const& )
190 #if defined( FEELPP_HAS_QD_QD_H )
193 operator<<( DebugStream& __os, node<dd_real>::type
const& __n )
195 if ( __os.doPrint() )
197 std::ostringstream __str;
201 __os << __str.str() <<
"\n";
208 operator<<( NdebugStream& __os, node<dd_real>::type
const& __n )
215 operator<<( DebugStream& __os, node<qd_real>::type
const& __n )
217 if ( __os.doPrint() )
219 std::ostringstream __str;
223 __os << __str.str() <<
"\n";
230 operator<<( NdebugStream& __os, node<qd_real>::type
const& __n )
241 operator<<( DebugStream& __os, ublas::vector<T>
const& __n )
243 if ( __os.doPrint() )
245 std::ostringstream __str;
249 __os << __str.str() <<
"\n";
257 operator<<( NdebugStream& __os, ublas::vector<T>
const& )
261 template<
typename T,
typename Orient>
264 operator<<( DebugStream& __os, ublas::matrix<T, Orient>
const& __n )
266 if ( __os.doPrint() )
268 std::ostringstream __str;
272 __os << __str.str() <<
"\n";
277 template<
typename T,
typename Orient>
280 operator<<( NdebugStream& __os, ublas::matrix<T,Orient>
const& )
289 #if defined( FEELPP_SIZET_SAME_AS_UINT )
290 typedef ublas::compressed_matrix<double,
291 ublas::row_major > csr_matrix_type;
292 typedef ublas::compressed_matrix<double,
293 ublas::column_major > csc_matrix_type;
295 typedef ublas::compressed_matrix<double,
297 ublas::unbounded_array<int>,
298 ublas::unbounded_array<double> > csr_matrix_type;
299 typedef ublas::compressed_matrix<double,
300 ublas::column_major, 0,
301 ublas::unbounded_array<int>,
302 ublas::unbounded_array<double> > csc_matrix_type;
309 template<
typename MatrixType>
310 void spy( MatrixType
const& __m, std::string
const &filename )
312 std::string name = filename;
313 std::string separator =
" , ";
316 int i = filename.find(
"." );
319 name = filename +
".m";
323 if ( (
unsigned int ) i != filename.size() - 2 ||
324 filename[ i + 1 ] !=
'm' )
326 std::cerr <<
"Wrong file name ";
327 name = filename +
".m";
331 std::ofstream file_out( name.c_str() );
333 FEELPP_ASSERT( file_out )( filename ).error(
"[Feel::spy] ERROR: File cannot be opened for writing." );
335 file_out <<
"S = [ ";
337 for (
typename MatrixType::const_iterator1 i1=__m.begin1();
338 i1!=__m.end1(); ++i1 )
340 for (
typename MatrixType::const_iterator2 i2=i1.begin();
342 file_out << i2.index1() + 1 << separator
343 << i2.index2() + 1 << separator
347 file_out <<
"];" << std::endl;
348 file_out <<
"I=S(:,1); J=S(:,2); S=S(:,3);" << std::endl;
349 file_out <<
"A=sparse(I,J,S); spy(A);" << std::endl;
354 namespace ublas = boost::numeric::ublas;
356 template<
typename T,
typename Orien>
358 ublas::matrix<T,Orien>
359 average( ublas::matrix<T, Orien>
const& m )
361 ublas::matrix<T, Orien> v( m.size1(), 1 );
362 ublas::scalar_vector<T> avg( m.size2(), T( 1 ) );
363 T n_val = int( m.size2() );
365 for (
size_type i = 0; i < v.size1(); ++i )
366 v( i, 0 ) = ublas::inner_prod( ublas::row( m, i ), avg )/n_val;
374 typename T::value_type
const& treshold = type_traits<typename T::value_type>::epsilon(),
375 typename T::value_type
const& new_value =
typename T::value_type( 0.0 ) )
377 std::for_each( t.data().begin(),
379 lambda::if_then( lambda::_1 < lambda::constant( treshold ) && lambda::_1 > -lambda::constant( treshold ),
380 lambda::_1 = lambda::constant( new_value ) ) );
386 linspace( T
const& __a, T
const& __b,
size_type __N,
int interior = 0 )
389 ublas::vector<T> v( N );
390 T h = ( __b-__a )/T( __N-1 );
391 T a = __a+T( interior )*h;
395 std::for_each( v.begin(),
397 ( lambda::_1 = a+lambda::var( i )*h, ++lambda::var( i ) ) );
412 typedef typename T::value_type value_type;
413 typedef boost::minstd_rand base_generator_type;
414 base_generator_type generator( 42u );
415 boost::uniform_real<value_type> uni_dist( 0.0,1.0 );
416 typedef boost::variate_generator<base_generator_type&, boost::uniform_real<value_type> > rand_type;
417 rand_type uni( generator, uni_dist );
419 std::for_each( t.data().begin(),
421 lambda::_1 = lambda::bind<value_type>( uni ) );