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
vectorublas.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: 2006-11-13
7 
8  Copyright (C) 2005,2006 EPFL
9  Copyright (C) 2007-2010 Universit� Joseph Fourier (Grenoble I)
10 
11  This library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU Lesser General Public
13  License as published by the Free Software Foundation; either
14  version 3.0 of the License, or (at your option) any later version.
15 
16  This library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  Lesser General Public License for more details.
20 
21  You should have received a copy of the GNU Lesser General Public
22  License along with this library; if not, write to the Free Software
23  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24 */
30 #ifndef __VectorUblas_H
31 #define __VectorUblas_H 1
32 
33 #include <set>
34 #include <boost/operators.hpp>
35 #include <boost/numeric/ublas/vector.hpp>
36 #include <boost/numeric/ublas/vector_proxy.hpp>
38 #include <feel/feelalg/vector.hpp>
39 
40 
41 
42 namespace Feel
43 {
44 namespace ublas = boost::numeric::ublas;
56 template<typename T, typename Storage = ublas::vector<T> >
57 class VectorUblas
58  : public Vector<T>
59  , boost::addable<VectorUblas<T,Storage> >
60  , boost::subtractable<VectorUblas<T,Storage> >
61  , boost::multipliable<VectorUblas<T,Storage>, T >
62 {
63  typedef Vector<T> super1;
64 public:
65 
66 
70 
71  typedef T value_type;
72  typedef typename type_traits<value_type>::real_type real_type;
73 
74  typedef Storage vector_type;
75  typedef typename vector_type::difference_type difference_type;
76  typedef ublas::basic_range<size_type, difference_type> range_type;
77  typedef ublas::basic_slice<size_type, difference_type> slice_type;
78  typedef Vector<value_type> clone_type;
79  typedef boost::shared_ptr<clone_type> clone_ptrtype;
80  typedef VectorUblas<value_type, Storage> this_type;
81 
82  typedef typename vector_type::iterator iterator;
83  typedef typename vector_type::const_iterator const_iterator;
84 
85  struct range
86  {
87  typedef ublas::vector_range<ublas::vector<value_type> > subtype;
88  typedef VectorUblas<value_type,subtype> type;
89  };
90 
91  struct slice
92  {
93  typedef ublas::vector_slice<ublas::vector<value_type> > subtype;
94  typedef VectorUblas<value_type,subtype> type;
95  };
96 
97 
98  typedef typename super1::datamap_type datamap_type;
99  typedef typename super1::datamap_ptrtype datamap_ptrtype;
100 
102 
106 
107  VectorUblas();
108 
109  VectorUblas( size_type __s );
110 
111  VectorUblas( datamap_ptrtype const& dm );
112 
113  VectorUblas( size_type __s, size_type __n_local );
114 
115  VectorUblas( VectorUblas const & m );
116 
117  VectorUblas( VectorUblas<value_type>& m, range_type const& range, datamap_ptrtype const& dm );
118 
119  VectorUblas( ublas::vector<value_type>& m, range_type const& range );
120 
121  VectorUblas( VectorUblas<value_type>& m, slice_type const& slice );
122 
123  VectorUblas( ublas::vector<value_type>& m, slice_type const& slice );
124 
125 
126  ~VectorUblas();
127 
140  void init ( const size_type N,
141  const size_type n_local,
142  const bool fast=false );
143 
147  void init ( const size_type n,
148  const bool fast=false );
149 
153  void init( datamap_ptrtype const& dm );
154 
155 
160  clone_ptrtype clone() const
161  {
162  return clone_ptrtype( new this_type( *this ) );
163  }
164 
166 
170 
175 
176  template<typename AE>
177  VectorUblas<value_type, Storage>& operator=( ublas::vector_expression<AE> const& e )
178  {
179  this->outdateGlobalValues();
180  M_vec.operator=( e );
181  return *this;
182  }
183 
187  T operator()( size_type i ) const
188  {
189  FEELPP_ASSERT ( this->isInitialized() ).error( "vector not initialized" );
190  FEELPP_ASSERT ( ( i >= this->firstLocalIndex() ) &&
191  ( i < this->lastLocalIndex() ) )
192  ( i )
193  ( this->firstLocalIndex() )
194  ( this->lastLocalIndex() ).error( "vector invalid index" );
195 
196  return M_vec.operator()( i-this->firstLocalIndex() );
197  }
198 
203  {
204  FEELPP_ASSERT ( this->isInitialized() ).error( "vector not initialized" );
205  FEELPP_ASSERT ( ( i >= this->firstLocalIndex() ) &&
206  ( i < this->lastLocalIndex() ) )
207  ( i )
208  ( this->firstLocalIndex() )
209  ( this->lastLocalIndex() ).error( "vector invalid index" );
210  this->outdateGlobalValues();
211  return M_vec.operator()( i-this->firstLocalIndex() );
212  }
213 
217  T operator[]( size_type i ) const
218  {
219  FEELPP_ASSERT ( this->isInitialized() ).error( "vector not initialized" );
220  FEELPP_ASSERT ( ( i >= this->firstLocalIndex() ) &&
221  ( i < this->lastLocalIndex() ) )
222  ( i )
223  ( this->firstLocalIndex() )
224  ( this->lastLocalIndex() ).error( "vector invalid index" );
225 
226  return M_vec.operator()( i-this->firstLocalIndex() );
227  }
228 
233  {
234  FEELPP_ASSERT ( this->isInitialized() ).error( "vector not initialized" );
235  FEELPP_ASSERT ( ( i >= this->firstLocalIndex() ) &&
236  ( i <= this->lastLocalIndex() ) )
237  ( i )
238  ( this->firstLocalIndex() )
239  ( this->lastLocalIndex() ).error( "vector invalid index" );
240  this->outdateGlobalValues();
241  return M_vec.operator()( i-this->firstLocalIndex() );
242  }
243 
249  {
250  checkInvariant();
251  this->outdateGlobalValues();
252  add( 1., v );
253  return *this;
254  }
255 
261  {
262  checkInvariant();
263  this->outdateGlobalValues();
264  add( -1., v );
265 
266  return *this;
267  }
268 
272  Vector<T>& operator*=( T const& v )
273  {
274  checkInvariant();
275  this->outdateGlobalValues();
276  this->scale( v );
277 
278  return *this;
279  }
281 
285 
286  iterator begin()
287  {
288  return M_vec.begin();
289  }
290  const_iterator begin() const
291  {
292  return M_vec.begin();
293  }
294 
295  iterator end()
296  {
297  return M_vec.end();
298  }
299  const_iterator end() const
300  {
301  return M_vec.end();
302  }
303 
308  size_type start() const;
309 
314  unsigned int rowStart () const
315  {
316  checkInvariant();
317  return 0;
318  }
319 
325  {
326  checkInvariant();
327  return 0;
328  }
329 
333  bool isInitialized() const
334  {
335  return true;
336  }
337 
342  void close () const;
343 
344 
349  bool closed() const
350  {
351  return true;
352  }
353 
354 
358  vector_type const& vec () const
359  {
360  return M_vec;
361  }
362 
366  vector_type & vec ()
367  {
368  return M_vec;
369  }
370 
375  {
376  return M_global_values_updated;
377  }
378 
382  void updateGlobalValues() const
383  {
384  //this->localize( M_global_values );
385  M_global_values_updated = true;
386  }
387 
391  value_type globalValue( size_type i ) const
392  {
393  return this->operator()( i );
394  }
395  //{ return M_global_values( i ); }
396 
397  //@
398 
399 
403 
408  {
409  M_global_values_updated = false;
410  }
411 
415  void setGlobalValue( size_type i, value_type v ) const
416  {
417  //M_global_values( i ) = v;
418  }
419 
423  void setConstant( value_type v )
424  {
425  M_vec = ublas::scalar_vector<double>( M_vec.size(), v );
426  }
427 
429 
433 
437  void resize( size_type n, bool preserve = true );
438 
445  void clear ();
446 
451  void zero ()
452  {
453  //M_vec.clear();
454  this->outdateGlobalValues();
455  std::fill( this->begin(), this->end(), value_type( 0 ) );
456  }
457 
458  void zero ( size_type /*start1*/, size_type /*stop1*/ )
459  {
460  //ublas::project( (*this), ublas::range( start1, stop1 ) ) = ublas::zero_vector<value_type>( stop1 );
461  }
462 
466  void add ( const size_type i, const value_type& value )
467  {
468  checkInvariant();
469  this->outdateGlobalValues();
470  ( *this )( i ) += value;
471  }
472 
476  void addVector ( int* i, int n, value_type* v )
477  {
478  for ( int j = 0; j < n; ++j )
479  ( *this )( i[j] ) += v[j];
480  }
481 
485  void set ( size_type i, const value_type& value )
486  {
487  checkInvariant();
488  this->outdateGlobalValues();
489  ( *this )( i ) = value;
490  }
491 
497  void addVector ( const std::vector<value_type>& v,
498  const std::vector<size_type>& dof_indices )
499  {
500  FEELPP_ASSERT ( v.size() == dof_indices.size() ).error( "invalid dof indices" );
501  this->outdateGlobalValues();
502 
503  for ( size_type i=0; i<v.size(); i++ )
504  this->add ( dof_indices[i], v[i] );
505  }
506 
513  void addVector ( const Vector<value_type>& V,
514  const std::vector<size_type>& dof_indices )
515  {
516  FEELPP_ASSERT ( V.size() == dof_indices.size() ).error( "invalid dof indices" );
517  this->outdateGlobalValues();
518 
519  for ( size_type i=0; i<V.size(); i++ )
520  this->add ( dof_indices[i], V( i ) );
521  }
522 
523 
528  void addVector ( const Vector<value_type>& /*V_in*/,
529  const MatrixSparse<value_type>& /*A_in*/ )
530  {
531  FEELPP_ASSERT( 0 ).error( "invalid call, not implemented yet" );
532  }
533 
540  void addVector ( const ublas::vector<value_type>& V,
541  const std::vector<size_type>& dof_indices )
542  {
543  FEELPP_ASSERT ( V.size() == dof_indices.size() ).error( "invalid dof indices" );
544  this->outdateGlobalValues();
545 
546  for ( size_type i=0; i<V.size(); i++ )
547  this->add ( dof_indices[i], V( i ) );
548  }
549 
554  void insert ( const std::vector<T>& /*v*/,
555  const std::vector<size_type>& /*dof_indices*/ )
556  {
557  FEELPP_ASSERT( 0 ).error( "invalid call, not implemented yet" );
558  }
559 
566  void insert ( const Vector<T>& /*V*/,
567  const std::vector<size_type>& /*dof_indices*/ )
568  {
569  FEELPP_ASSERT( 0 ).error( "invalid call, not implemented yet" );
570  }
571 
572 
579  void insert ( const ublas::vector<T>& /*V*/,
580  const std::vector<size_type>& /*dof_indices*/ )
581  {
582  FEELPP_ASSERT( 0 ).error( "invalid call, not implemented yet" );
583  }
584 
589  void scale ( const T factor )
590  {
591  this->outdateGlobalValues();
592  M_vec.operator *=( factor );
593  }
594 
601  void printMatlab( const std::string name="NULL", bool renumber = false ) const;
602 
603  void close() {}
604 
609  real_type min() const
610  {
611  checkInvariant();
612 
613  real_type local_min = *std::min_element( M_vec.begin(), M_vec.end() );
614 
615  real_type global_min = local_min;
616 
617 
618 
619 #ifdef FEELPP_HAS_MPI
620 
621  if ( this->comm().size() > 1 )
622  {
623  MPI_Allreduce ( &local_min, &global_min, 1,
624  MPI_DOUBLE, MPI_MIN, this->comm() );
625  }
626 
627 #endif
628 
629  return global_min;
630  }
635  real_type max() const
636  {
637  checkInvariant();
638 
639  real_type local_max = *std::max_element( M_vec.begin(), M_vec.end() );
640 
641  real_type global_max = local_max;
642 
643 
644 #ifdef FEELPP_HAS_MPI
645 
646  if ( this->comm().size() > 1 )
647  {
648  MPI_Allreduce ( &local_max, &global_max, 1,
649  MPI_DOUBLE, MPI_MAX, this->comm() );
650  }
651 
652 #endif
653 
654  return global_max;
655  }
656 
661  real_type l1Norm() const
662  {
663  checkInvariant();
664  double local_l1 = ublas::norm_1( M_vec );
665 
666  double global_l1 = local_l1;
667 
668 
669 #ifdef FEELPP_HAS_MPI
670 
671  if ( this->comm().size() > 1 )
672  {
673  mpi::all_reduce( this->comm(), local_l1, global_l1, std::plus<double>() );
674  }
675 
676 #endif
677 
678  return global_l1;
679 
680  }
681 
686  real_type l2Norm() const
687  {
688  checkInvariant();
689  real_type local_norm2 = 0, global_norm2=0;
690 
691  if ( this->comm().size() == 1 )
692  {
693  local_norm2 = ublas::inner_prod( M_vec, M_vec );
694  global_norm2 = local_norm2;
695  }
696  else
697  {
698  size_type s = this->localSize();
699  size_type start = this->firstLocalIndex();
700  for ( size_type i = 0; i < s; ++i )
701  {
702  if ( !this->localIndexIsGhost( start + i ) )
703  local_norm2 += std::pow(M_vec.operator()( start + i ),2);
704  }
705 #ifdef FEELPP_HAS_MPI
706  mpi::all_reduce( this->comm(), local_norm2, global_norm2, std::plus<real_type>() );
707 #endif
708  }
709 
710  return math::sqrt( global_norm2 );
711  }
712 
717  real_type linftyNorm() const
718  {
719  checkInvariant();
720  real_type local_norminf = ublas::norm_inf( M_vec );
721  real_type global_norminf = local_norminf;
722 
723 
724 #ifdef FEELPP_HAS_MPI
725 
726  if ( this->comm().size() > 1 )
727  {
728  mpi::all_reduce( this->comm(), local_norminf, global_norminf, mpi::maximum<real_type>() );
729  }
730 
731 #endif
732  return global_norminf;
733  }
734 
735 
739  value_type sum() const
740  {
741  checkInvariant();
742  double local_sum = ublas::sum( M_vec );
743 
744  double global_sum = local_sum;
745 
746 
747 #ifdef FEELPP_HAS_MPI
748 
749  if ( this->comm().size() > 1 )
750  {
751  mpi::all_reduce( this->comm(), local_sum, global_sum, std::plus<value_type>() );
752  }
753 
754 #endif
755 
756  return global_sum;
757 
758  }
759 
763  FEELPP_DONT_INLINE this_type sqrt() const;
764 
765 
769  this_type pow( int n ) const;
770 
771 
777  void add( const T& a )
778  {
779  checkInvariant();
780  this->outdateGlobalValues();
781 
782  for ( size_type i = 0; i < this->localSize(); ++i )
783  M_vec.operator[]( i ) += a;
784 
785  return;
786  }
787 
792  void add( const Vector<T>& v )
793  {
794  checkInvariant();
795  this->outdateGlobalValues();
796  add( 1., v );
797  return;
798  }
799 
805  void add( const T& a, const Vector<T>& v )
806  {
807  checkInvariant();
808  this->outdateGlobalValues();
809 
810  for ( size_type i = 0; i < this->localSize(); ++i )
811  M_vec.operator()( i ) += a*v( v.firstLocalIndex() + i );
812 
813  return;
814  }
815 
820  void localize ( std::vector<value_type>& /*v_local*/ ) const
821  {
822  FEELPP_ASSERT( 0 ).error( "invalid call, not implemented yet" );
823  }
824 
829  void localize ( ublas::vector<value_type>& v_local ) const;
830 
835  void localize ( ublas::vector_range<ublas::vector<value_type> >& v_local ) const;
836 
841  void localize ( ublas::vector_slice<ublas::vector<value_type> >& v_local ) const;
842 
847  void localize ( Vector<T>& v_local ) const;
848 
854  void localize ( Vector<T>& v_local,
855  const std::vector<size_type>& send_list ) const;
856 
861  void localize ( const size_type first_local_idx,
862  const size_type last_local_idx,
863  const std::vector<size_type>& send_list );
864 
871  void localizeToOneProcessor ( ublas::vector<T>& v_local,
872  const size_type proc_id = 0 ) const;
873 
880  void localizeToOneProcessor ( std::vector<T>& v_local,
881  const size_type proc_id = 0 ) const;
882 
883 
884  value_type dot( Vector<T> const& __v );
886 
887 
888 
889 protected:
890 
891 private:
892 
896  void checkInvariant() const;
897 
898 private:
899  vector_type M_vec;
900  mutable bool M_global_values_updated;
901  //mutable ublas::vector<value_type> M_global_values;
902 };
903 
911 template <typename T>
912 VectorUblas<T>
914 {
915  FEELPP_ASSERT( v1.localSize() == v2.localSize() &&
916  v1.size() == v2.size() )
917  ( v1.localSize() )( v2.localSize() )
918  ( v1.size() )( v2.size() ).error( "incompatible vector sizes" );
919 
920  typedef typename type_traits<T>::real_type real_type;
921 
922  VectorUblas<real_type> _t( v1.mapPtr() );
923  size_type s = v1.localSize();
924  size_type start = v1.firstLocalIndex();
925 
926  for ( size_type i = 0; i < s; ++i )
927  _t.operator()( start+i ) = v1.operator()( start + i )* v2.operator()( start + i );
928 
929  return _t;
930 }
931 
939 template <typename T>
940 VectorUblas<T>
941 element_product( boost::shared_ptr<VectorUblas<T> > const& v1,
942  boost::shared_ptr<VectorUblas<T> > const& v2 )
943 {
944  return element_product( *v1, *v2 );
945 }
946 
952 #if !defined( FEELPP_INSTANTIATE_VECTORUBLAS )
953 extern template class VectorUblas<double,ublas::vector<double> >;
954 extern template class VectorUblas<double,ublas::vector_range<ublas::vector<double> > >;
955 extern template class VectorUblas<double,ublas::vector_slice<ublas::vector<double> > >;
956 #endif
957 
958 } // Feel
959 
960 
961 #endif /* __VectorUblas_H */

Generated on Sun Oct 20 2013 08:25:06 for Feel++ by doxygen 1.8.4