Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
glas.hpp
Go to the documentation of this file.
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2005-08-17
7 
8  Copyright (C) 2005,2006 EPFL
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 3.0 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
29 #ifndef __FEELPP_GLAS_HPP
30 #define __FEELPP_GLAS_HPP 1
31 
32 #include <string>
33 #include <sstream>
34 #include <fstream>
35 
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>
40 
41 #include <boost/lambda/lambda.hpp>
42 #include <boost/lambda/algorithm.hpp>
43 #include <boost/lambda/bind.hpp>
44 #include <boost/lambda/if.hpp>
45 
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>
53 
54 #include <feel/feelcore/feel.hpp>
55 #include <feel/feelcore/debug.hpp>
56 #include <feel/feelcore/traits.hpp>
57 
58 namespace Feel
59 {
60 const double Pi = 3.14159265358979323846264338328;
61 const double TGV = 1e20;
62 
63 template <class T>
64 inline T Min ( const T &a, const T &b )
65 {
66  return a < b ? a : b;
67 }
68 template <class T>
69 inline T Max ( const T &a, const T & b )
70 {
71  return a > b ? a : b;
72 }
73 template <class T>
74 inline T Abs ( const T &a )
75 {
76  return a < 0 ? -a : a;
77 }
78 
79 template <class T>
80 inline void Exchange ( T& a, T& b )
81 {
82  T c = a;
83  a = b;
84  b = c;
85 }
86 template <class T>
87 inline T Max ( const T &a, const T & b, const T & c )
88 {
89  return Max( Max( a, b ), c );
90 }
91 template <class T>
92 inline T Min ( const T &a, const T & b, const T & c )
93 {
94  return Min( Min( a, b ), c );
95 }
96 
97 namespace ublas = boost::numeric::ublas;
98 
99 
100 
101 /* reduce operations */
102 using ublas::sum;
103 using ublas::norm_1;
104 using ublas::norm_2;
105 using ublas::norm_inf;
106 using ublas::index_norm_inf;
107 
108 /* binary operations */
109 using ublas::outer_prod;
110 using ublas::inner_prod;
111 
112 struct norm_inf_adaptor
113 {
114  template<typename E>
115  Real operator()( E const& __v ) const
116  {
117  return norm_inf( __v );
118  }
119 };
120 
124 template <typename T = double, uint16_type S = 3>
125 struct node
126 {
127  //typedef ublas::vector<T, ublas::bounded_array<T, S> > type;
128  typedef ublas::vector<T> type;
129 };
130 #if 0
131 
135 typedef node<double, 3>::type node_type;
136 
140 typedef node<double, 3>::type gradient_node_type;
141 #endif
142 
145 typedef ublas::symmetric_matrix<double, ublas::lower, ublas::row_major, ublas::bounded_array<double, 9> > hessian_node_type;
146 
147 
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;
150 
151 typedef lapack_matrix_type transformation_matrix_type;
152 
156 template <typename T = double, uint16_type S = 256>
158 {
159  //typedef ublas::matrix<T, ublas::column_major, ublas::bounded_array<T, S> > type;
160  typedef ublas::matrix<T, ublas::column_major> type;
161 };
162 
163 #if 0
164 typedef matrix_node<double>::type matrix_node_type;
165 #endif
166 
167 inline
168 DebugStream&
169 operator<<( DebugStream& __os, node<real64_type>::type const& __n )
170 {
171  if ( __os.doPrint() )
172  {
173  std::ostringstream __str;
174 
175  __str << __n;
176 
177  __os << __str.str() << "\n";
178  }
179 
180  return __os;
181 }
182 inline
183 NdebugStream&
184 operator<<( NdebugStream& os, node<real64_type>::type const& /*n*/ )
185 {
186  return os;
187 }
188 
189 
190 #if defined( FEELPP_HAS_QD_QD_H )
191 inline
192 DebugStream&
193 operator<<( DebugStream& __os, node<dd_real>::type const& __n )
194 {
195  if ( __os.doPrint() )
196  {
197  std::ostringstream __str;
198 
199  __str << __n;
200 
201  __os << __str.str() << "\n";
202  }
203 
204  return __os;
205 }
206 inline
207 NdebugStream&
208 operator<<( NdebugStream& __os, node<dd_real>::type const& __n )
209 {
210  return __os;
211 }
212 
213 inline
214 DebugStream&
215 operator<<( DebugStream& __os, node<qd_real>::type const& __n )
216 {
217  if ( __os.doPrint() )
218  {
219  std::ostringstream __str;
220 
221  __str << __n;
222 
223  __os << __str.str() << "\n";
224  }
225 
226  return __os;
227 }
228 inline
229 NdebugStream&
230 operator<<( NdebugStream& __os, node<qd_real>::type const& __n )
231 {
232  return __os;
233 }
234 
235 
236 #endif /* FEELPP_HAS_QD_QD_H */
237 
238 template<typename T>
239 inline
240 DebugStream&
241 operator<<( DebugStream& __os, ublas::vector<T> const& __n )
242 {
243  if ( __os.doPrint() )
244  {
245  std::ostringstream __str;
246 
247  __str << __n;
248 
249  __os << __str.str() << "\n";
250  }
251 
252  return __os;
253 }
254 template<typename T>
255 inline
256 NdebugStream&
257 operator<<( NdebugStream& __os, ublas::vector<T> const& /*__n*/ )
258 {
259  return __os;
260 }
261 template<typename T, typename Orient>
262 inline
263 DebugStream&
264 operator<<( DebugStream& __os, ublas::matrix<T, Orient> const& __n )
265 {
266  if ( __os.doPrint() )
267  {
268  std::ostringstream __str;
269 
270  __str << __n;
271 
272  __os << __str.str() << "\n";
273  }
274 
275  return __os;
276 }
277 template<typename T,typename Orient>
278 inline
279 NdebugStream&
280 operator<<( NdebugStream& __os, ublas::matrix<T,Orient> const& /*__n*/ )
281 {
282  return __os;
283 }
284 
285 
286 //
287 // sparse matrices
288 //
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;
294 #else
295 typedef ublas::compressed_matrix<double,
296  ublas::row_major, 0,
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;
303 #endif
304 
309 template<typename MatrixType>
310 void spy( MatrixType const& __m, std::string const &filename )
311 {
312  std::string name = filename;
313  std::string separator = " , ";
314 
315  // check on the file name
316  int i = filename.find( "." );
317 
318  if ( i <= 0 )
319  name = filename + ".m";
320 
321  else
322  {
323  if ( ( unsigned int ) i != filename.size() - 2 ||
324  filename[ i + 1 ] != 'm' )
325  {
326  std::cerr << "Wrong file name ";
327  name = filename + ".m";
328  }
329  }
330 
331  std::ofstream file_out( name.c_str() );
332 
333  FEELPP_ASSERT( file_out )( filename ).error( "[Feel::spy] ERROR: File cannot be opened for writing." );
334 
335  file_out << "S = [ ";
336 
337  for ( typename MatrixType::const_iterator1 i1=__m.begin1();
338  i1!=__m.end1(); ++i1 )
339  {
340  for ( typename MatrixType::const_iterator2 i2=i1.begin();
341  i2!=i1.end(); ++i2 )
342  file_out << i2.index1() + 1 << separator
343  << i2.index2() + 1 << separator
344  << *i2 << std::endl;
345  }
346 
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;
350 }
351 
352 namespace glas
353 {
354 namespace ublas = boost::numeric::ublas;
355 
356 template<typename T, typename Orien>
357 inline
358 ublas::matrix<T,Orien>
359 average( ublas::matrix<T, Orien> const& m )
360 {
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() );
364 
365  for ( size_type i = 0; i < v.size1(); ++i )
366  v( i, 0 ) = ublas::inner_prod( ublas::row( m, i ), avg )/n_val;
367 
368  return v;
369 }
370 template<typename T>
371 inline
372 void
373 clean( T& t,
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 ) )
376 {
377  std::for_each( t.data().begin(),
378  t.data().end(),
379  lambda::if_then( lambda::_1 < lambda::constant( treshold ) && lambda::_1 > -lambda::constant( treshold ),
380  lambda::_1 = lambda::constant( new_value ) ) );
381 }
382 
383 template<typename T>
384 inline
385 ublas::vector<T>
386 linspace( T const& __a, T const& __b, size_type __N, int interior = 0 )
387 {
388  size_type N = __N-2*interior;
389  ublas::vector<T> v( N );
390  T h = ( __b-__a )/T( __N-1 );
391  T a = __a+T( interior )*h;
392  //T b = __b-T(interior)*h;
393 #if 0
394  size_type i = 0;
395  std::for_each( v.begin(),
396  v.end(),
397  ( lambda::_1 = a+lambda::var( i )*h, ++lambda::var( i ) ) );
398 #else
399 
400  for ( size_type i = 0; i < N; ++i )
401  v[i]=a+T( i )*h;
402 
403 #endif
404  return v;
405 }
406 
407 template<typename T>
408 inline
409 void
410 randomize( T& t )
411 {
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 );
418 
419  std::for_each( t.data().begin(),
420  t.data().end(),
421  lambda::_1 = lambda::bind<value_type>( uni ) );
422 
423 }
424 
425 //
426 //
427 //
428 } // namespace glas
429 } // namespace Feel
430 
431 #endif /* __FEELPP_GLAS_HPP */

Generated on Sun Oct 20 2013 08:24:59 for Feel++ by doxygen 1.8.4