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
vectoreigen.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: 2012-06-20
7 
8  Copyright (C) 2012 Universit� Joseph Fourier (Grenoble I)
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 __VectorEigen_H
30 #define __VectorEigen_H 1
31 
32 #include <set>
33 #include <boost/operators.hpp>
34 #include <Eigen/Core>
36 #include <feel/feelalg/vector.hpp>
37 
38 
39 
40 namespace Feel
41 {
42 template<typename T>
43 struct get_real_type
44 {};
45 
46 template<>
47 struct get_real_type<double>
48 {
49  typedef double value_type;
50 };
51 
52 template<>
53 struct get_real_type<std::complex<double> >
54 {
55  typedef double value_type;
56 };
57 
69 template<typename T >
71  : public Vector<T>
72  , boost::addable<VectorEigen<T> >
73  , boost::subtractable<VectorEigen<T> >
74 
75 {
76  typedef Vector<T> super1;
77 public:
78 
79 
83 
84  typedef T value_type;
85  typedef typename get_real_type<T>::value_type real_type;
86 
87  typedef Eigen::Matrix<value_type,Eigen::Dynamic,1> vector_type;
89  typedef typename super1::clone_ptrtype clone_ptrtype;
90 
91  typedef typename super1::datamap_type datamap_type;
92  typedef typename super1::datamap_ptrtype datamap_ptrtype;
93 
95 
99 
100  VectorEigen();
101 
102  VectorEigen( size_type __s );
103 
104  VectorEigen( datamap_ptrtype const& dm );
105 
106  VectorEigen( size_type __s, size_type __n_local );
107 
108  VectorEigen( VectorEigen const & m );
109 
110  ~VectorEigen();
111 
116  clone_ptrtype clone () const ;
117 
130  void init ( const size_type N,
131  const size_type n_local,
132  const bool fast=false );
133 
137  void init ( const size_type n,
138  const bool fast=false );
139 
143  void init( datamap_ptrtype const& dm );
144 
145 
147 
150 
155 
159  T operator()( size_type i ) const
160  {
161  FEELPP_ASSERT ( this->isInitialized() ).error( "vector not initialized" );
162  FEELPP_ASSERT ( ( i >= this->firstLocalIndex() ) &&
163  ( i < this->lastLocalIndex() ) )
164  ( i )
165  ( this->firstLocalIndex() )
166  ( this->lastLocalIndex() ).error( "vector invalid index" );
167 
168  return M_vec.operator()( i-this->firstLocalIndex() );
169  }
170 
175  {
176  FEELPP_ASSERT ( this->isInitialized() ).error( "vector not initialized" );
177  FEELPP_ASSERT ( ( i >= this->firstLocalIndex() ) &&
178  ( i < this->lastLocalIndex() ) )
179  ( i )
180  ( this->firstLocalIndex() )
181  ( this->lastLocalIndex() ).error( "vector invalid index" );
182  return M_vec.operator()( i-this->firstLocalIndex() );
183  }
184 
188  T operator[]( size_type i ) const
189  {
190  FEELPP_ASSERT ( this->isInitialized() ).error( "vector not initialized" );
191  FEELPP_ASSERT ( ( i >= this->firstLocalIndex() ) &&
192  ( i < this->lastLocalIndex() ) )
193  ( i )
194  ( this->firstLocalIndex() )
195  ( this->lastLocalIndex() ).error( "vector invalid index" );
196 
197  return M_vec.operator()( i-this->firstLocalIndex() );
198  }
199 
204  {
205  FEELPP_ASSERT ( this->isInitialized() ).error( "vector not initialized" );
206  FEELPP_ASSERT ( ( i >= this->firstLocalIndex() ) &&
207  ( i <= this->lastLocalIndex() ) )
208  ( i )
209  ( this->firstLocalIndex() )
210  ( this->lastLocalIndex() ).error( "vector invalid index" );
211  return M_vec.operator()( i-this->firstLocalIndex() );
212  }
213 
219  {
220  checkInvariant();
221  add( 1., v );
222  return *this;
223  }
224 
230  {
231  checkInvariant();
232  add( -1., v );
233 
234  return *this;
235  }
236 
238 
242 
243 
248  size_type start() const;
249 
254  unsigned int rowStart () const
255  {
256  checkInvariant();
257  return 0;
258  }
259 
265  {
266  checkInvariant();
267  return 0;
268  }
269 
273  bool isInitialized() const
274  {
275  return true;
276  }
277 
282  void close () const;
283 
284 
289  bool closed() const
290  {
291  return true;
292  }
293 
294 
298  vector_type const& vec () const
299  {
300  return M_vec;
301  }
302 
306  vector_type & vec ()
307  {
308  return M_vec;
309  }
310 
311 
312  //@
313 
314 
318 
319 
323  void setConstant( value_type v )
324  {
325  M_vec.setConstant( v );
326  }
327 
329 
333 
337  void resize( size_type n, bool preserve = true );
338 
345  void clear ();
346 
351  void zero ()
352  {
353  M_vec.setZero( M_vec.size() );
354  }
355 
356  void zero ( size_type /*start1*/, size_type /*stop1*/ )
357  {
358  //eigen::project( (*this), eigen::range( start1, stop1 ) ) = eigen::zero_vector<value_type>( stop1 );
359  }
360 
364  void add ( const size_type i, const value_type& value )
365  {
366  checkInvariant();
367  M_vec( i ) += value;
368  }
369 
373  void addVector ( int* i, int n, value_type* v )
374  {
375  for ( int j = 0; j < n; ++j )
376  M_vec( i[j] ) += v[j];
377  }
378 
382  void set ( size_type i, const value_type& value )
383  {
384  checkInvariant();
385  M_vec( i ) = value;
386  }
387 
393  void addVector ( const std::vector<value_type>& v,
394  const std::vector<size_type>& dof_indices )
395  {
396  FEELPP_ASSERT ( v.size() == dof_indices.size() ).error( "invalid dof indices" );
397 
398  for ( size_type i=0; i<v.size(); i++ )
399  this->add ( dof_indices[i], v[i] );
400  }
401 
408  void addVector ( const Vector<value_type>& V,
409  const std::vector<size_type>& dof_indices )
410  {
411  FEELPP_ASSERT ( V.size() == dof_indices.size() ).error( "invalid dof indices" );
412 
413  for ( size_type i=0; i<V.size(); i++ )
414  this->add ( dof_indices[i], V( i ) );
415  }
416 
417 
422  void addVector ( const Vector<value_type>& /*V_in*/,
423  const MatrixSparse<value_type>& /*A_in*/ )
424  {
425  FEELPP_ASSERT( 0 ).error( "invalid call, not implemented yet" );
426  }
427 
434  void addVector ( const vector_type& V,
435  const std::vector<size_type>& dof_indices )
436  {
437  FEELPP_ASSERT ( V.size() == dof_indices.size() ).error( "invalid dof indices" );
438 
439  for ( size_type i=0; i<V.size(); i++ )
440  this->add ( dof_indices[i], V( i ) );
441  }
442 
447  void insert ( const std::vector<T>& /*v*/,
448  const std::vector<size_type>& /*dof_indices*/ )
449  {
450  FEELPP_ASSERT( 0 ).error( "invalid call, not implemented yet" );
451  }
452 
459  void insert ( const Vector<T>& /*V*/,
460  const std::vector<size_type>& /*dof_indices*/ )
461  {
462  FEELPP_ASSERT( 0 ).error( "invalid call, not implemented yet" );
463  }
464 
465 
472  void insert ( const vector_type& /*V*/,
473  const std::vector<size_type>& /*dof_indices*/ )
474  {
475  FEELPP_ASSERT( 0 ).error( "invalid call, not implemented yet" );
476  }
477 
484  void insert ( const ublas::vector<T>& V,
485  const std::vector<size_type>& dof_indices );
486 
487 
492  void scale ( const T factor )
493  {
494  M_vec *= factor;
495  }
496 
503  void printMatlab( const std::string name="NULL", bool renumber = false ) const;
504 
505  void close() {}
506 
511  real_type min() const
512  {
513  checkInvariant();
514 
515  real_type local_min = 0;//M_vec.minCoeff();
516 
517  real_type global_min = local_min;
518 
519 
520 
521 #ifdef FEELPP_HAS_MPI
522 
523  if ( this->comm().size() > 1 )
524  {
525  MPI_Allreduce ( &local_min, &global_min, 1,
526  MPI_DOUBLE, MPI_MIN, this->comm() );
527  }
528 
529 #endif
530 
531  return global_min;
532  }
537  real_type max() const
538  {
539  checkInvariant();
540 
541  real_type local_max = 0;//M_vec.maxCoeff();
542 
543  real_type global_max = local_max;
544 
545 
546 #ifdef FEELPP_HAS_MPI
547 
548  if ( this->comm().size() > 1 )
549  {
550  MPI_Allreduce ( &local_max, &global_max, 1,
551  MPI_DOUBLE, MPI_MAX, this->comm() );
552  }
553 
554 #endif
555 
556  return global_max;
557  }
558 
563  real_type l1Norm() const
564  {
565  checkInvariant();
566  double local_l1 = M_vec.array().abs().sum();
567 
568  double global_l1 = local_l1;
569 
570 
571 #ifdef FEELPP_HAS_MPI
572 
573  if ( this->comm().size() > 1 )
574  {
575  mpi::all_reduce( this->comm(), local_l1, global_l1, std::plus<double>() );
576  }
577 
578 #endif
579 
580  return global_l1;
581 
582  }
583 
588  real_type l2Norm() const
589  {
590  checkInvariant();
591  real_type local_norm2 = M_vec.squaredNorm();
592  real_type global_norm2 = local_norm2;
593 
594 
595 #ifdef FEELPP_HAS_MPI
596 
597  if ( this->comm().size() > 1 )
598  {
599  mpi::all_reduce( this->comm(), local_norm2, global_norm2, std::plus<real_type>() );
600  }
601 
602 #endif
603  return math::sqrt( global_norm2 );
604  }
605 
610  real_type linftyNorm() const
611  {
612  checkInvariant();
613  real_type local_norminf = M_vec.array().abs().maxCoeff();
614  real_type global_norminf = local_norminf;
615 
616 
617 #ifdef FEELPP_HAS_MPI
618 
619  if ( this->comm().size() > 1 )
620  {
621  mpi::all_reduce( this->comm(), local_norminf, global_norminf, mpi::maximum<real_type>() );
622  }
623 
624 #endif
625  return global_norminf;
626  }
627 
628 
632  value_type sum() const
633  {
634  checkInvariant();
635  value_type local_sum = M_vec.array().sum();
636 
637  value_type global_sum = local_sum;
638 
639 
640 #ifdef FEELPP_HAS_MPI
641 
642  if ( this->comm().size() > 1 )
643  {
644  mpi::all_reduce( this->comm(), local_sum, global_sum, std::plus<value_type>() );
645  }
646 
647 #endif
648 
649  return global_sum;
650 
651  }
652 
656  FEELPP_DONT_INLINE this_type sqrt() const;
657 
658 
662  this_type pow( int n ) const;
663 
664 
670  void add( const T& a )
671  {
672  checkInvariant();
673 
674  for ( size_type i = 0; i < this->localSize(); ++i )
675  M_vec.operator[]( i ) += a;
676 
677  return;
678  }
679 
684  void add( const Vector<T>& v )
685  {
686  checkInvariant();
687  add( 1., v );
688  return;
689  }
690 
696  void add( const T& a, const Vector<T>& v )
697  {
698  checkInvariant();
699 
700  for ( size_type i = 0; i < this->localSize(); ++i )
701  M_vec.operator()( i ) += a*v( v.firstLocalIndex() + i );
702 
703  return;
704  }
705 
710  void localize ( std::vector<value_type>& /*v_local*/ ) const
711  {
712  FEELPP_ASSERT( 0 ).error( "invalid call, not implemented yet" );
713  }
714 
719  void localize ( const size_type first_local_idx,
720  const size_type last_local_idx,
721  const std::vector<size_type>& send_list );
726  void localize ( Vector<T>& v_local ) const;
727  void localize ( vector_type& v_local ) const;
728 
734  void localize ( Vector<T>& v_local,
735  const std::vector<size_type>& send_list ) const;
736 
737 
738 
745  void localizeToOneProcessor ( std::vector<T>& v_local,
746  const size_type proc_id = 0 ) const;
747  void localizeToOneProcessor ( vector_type& v_local,
748  const size_type proc_id = 0 ) const;
749 
750 
751  value_type dot( Vector<T> const& __v )
752  {
753  throw std::logic_error( "[vetor eigen] ERROR dot function not yet implemented" );
754  return 0;
755  }
757 
758 
759 
760 protected:
761 
762 private:
763 
767  void checkInvariant() const;
768 
769 private:
770  vector_type M_vec;
771 };
772 
780 template <typename T>
781 VectorEigen<T>
783 {
784  return v1.vec().array()*v2.vec().array();
785 }
786 
794 template <typename T>
795 VectorEigen<T>
796 element_product( boost::shared_ptr<VectorEigen<T> > const& v1,
797  boost::shared_ptr<VectorEigen<T> > const& v2 )
798 {
799  return v1->vec().array()*v2->vec().array();
800 }
801 
807 #if !defined( FEELPP_INSTANTIATE_VECTOREIGEN )
808 extern template class VectorEigen<double>;
809 #endif
810 
811 } // Feel
812 
813 
814 #endif /* __VectorEigen_H */

Generated on Fri Oct 25 2013 14:24:26 for Feel++ by doxygen 1.8.4