30 #ifndef __VectorUblas_H
31 #define __VectorUblas_H 1
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>
44 namespace ublas = boost::numeric::ublas;
56 template<
typename T,
typename Storage = ublas::vector<T> >
59 , boost::addable<VectorUblas<T,Storage> >
60 , boost::subtractable<VectorUblas<T,Storage> >
61 , boost::multipliable<VectorUblas<T,Storage>, T >
63 typedef Vector<T> super1;
72 typedef typename type_traits<value_type>::real_type real_type;
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;
82 typedef typename vector_type::iterator iterator;
83 typedef typename vector_type::const_iterator const_iterator;
87 typedef ublas::vector_range<ublas::vector<value_type> > subtype;
88 typedef VectorUblas<value_type,subtype> type;
93 typedef ublas::vector_slice<ublas::vector<value_type> > subtype;
94 typedef VectorUblas<value_type,subtype> type;
98 typedef typename super1::datamap_type datamap_type;
99 typedef typename super1::datamap_ptrtype datamap_ptrtype;
111 VectorUblas( datamap_ptrtype
const& dm );
115 VectorUblas( VectorUblas
const & m );
117 VectorUblas( VectorUblas<value_type>& m, range_type
const& range, datamap_ptrtype
const& dm );
119 VectorUblas( ublas::vector<value_type>& m, range_type
const& range );
121 VectorUblas( VectorUblas<value_type>& m, slice_type
const& slice );
123 VectorUblas( ublas::vector<value_type>& m, slice_type
const& slice );
142 const bool fast=
false );
148 const bool fast=
false );
153 void init( datamap_ptrtype
const& dm );
162 return clone_ptrtype(
new this_type( *
this ) );
176 template<
typename AE>
180 M_vec.operator=( e );
189 FEELPP_ASSERT ( this->
isInitialized() ).error(
"vector not initialized" );
204 FEELPP_ASSERT ( this->
isInitialized() ).error(
"vector not initialized" );
219 FEELPP_ASSERT ( this->
isInitialized() ).error(
"vector not initialized" );
234 FEELPP_ASSERT ( this->
isInitialized() ).error(
"vector not initialized" );
288 return M_vec.begin();
290 const_iterator begin()
const
292 return M_vec.begin();
299 const_iterator end()
const
358 vector_type
const&
vec ()
const
376 return M_global_values_updated;
385 M_global_values_updated =
true;
409 M_global_values_updated =
false;
425 M_vec = ublas::scalar_vector<double>( M_vec.size(), v );
455 std::fill( this->begin(), this->end(), value_type( 0 ) );
470 ( *this )( i ) += value;
478 for (
int j = 0; j < n; ++j )
479 ( *
this )( i[j] ) += v[j];
489 ( *this )( i ) = value;
498 const std::vector<size_type>& dof_indices )
500 FEELPP_ASSERT ( v.size() == dof_indices.size() ).error(
"invalid dof indices" );
504 this->
add ( dof_indices[i], v[i] );
514 const std::vector<size_type>& dof_indices )
516 FEELPP_ASSERT ( V.
size() == dof_indices.size() ).error(
"invalid dof indices" );
520 this->
add ( dof_indices[i], V( i ) );
531 FEELPP_ASSERT( 0 ).error(
"invalid call, not implemented yet" );
541 const std::vector<size_type>& dof_indices )
543 FEELPP_ASSERT ( V.size() == dof_indices.size() ).error(
"invalid dof indices" );
547 this->
add ( dof_indices[i], V( i ) );
555 const std::vector<size_type>& )
557 FEELPP_ASSERT( 0 ).error(
"invalid call, not implemented yet" );
567 const std::vector<size_type>& )
569 FEELPP_ASSERT( 0 ).error(
"invalid call, not implemented yet" );
580 const std::vector<size_type>& )
582 FEELPP_ASSERT( 0 ).error(
"invalid call, not implemented yet" );
592 M_vec.operator *=( factor );
601 void printMatlab(
const std::string name=
"NULL",
bool renumber =
false )
const;
613 real_type local_min = *std::min_element( M_vec.begin(), M_vec.end() );
615 real_type global_min = local_min;
619 #ifdef FEELPP_HAS_MPI
621 if ( this->
comm().size() > 1 )
623 MPI_Allreduce ( &local_min, &global_min, 1,
624 MPI_DOUBLE, MPI_MIN, this->
comm() );
639 real_type local_max = *std::max_element( M_vec.begin(), M_vec.end() );
641 real_type global_max = local_max;
644 #ifdef FEELPP_HAS_MPI
646 if ( this->
comm().size() > 1 )
648 MPI_Allreduce ( &local_max, &global_max, 1,
649 MPI_DOUBLE, MPI_MAX, this->
comm() );
664 double local_l1 = ublas::norm_1( M_vec );
666 double global_l1 = local_l1;
669 #ifdef FEELPP_HAS_MPI
673 mpi::all_reduce( this->
comm(), local_l1, global_l1, std::plus<double>() );
689 real_type local_norm2 = 0, global_norm2=0;
693 local_norm2 = ublas::inner_prod( M_vec, M_vec );
694 global_norm2 = local_norm2;
702 if ( !this->localIndexIsGhost( start + i ) )
703 local_norm2 += std::pow(M_vec.operator()( start + i ),2);
705 #ifdef FEELPP_HAS_MPI
706 mpi::all_reduce( this->
comm(), local_norm2, global_norm2, std::plus<real_type>() );
710 return math::sqrt( global_norm2 );
720 real_type local_norminf = ublas::norm_inf( M_vec );
721 real_type global_norminf = local_norminf;
724 #ifdef FEELPP_HAS_MPI
728 mpi::all_reduce( this->
comm(), local_norminf, global_norminf, mpi::maximum<real_type>() );
732 return global_norminf;
742 double local_sum = ublas::sum( M_vec );
744 double global_sum = local_sum;
747 #ifdef FEELPP_HAS_MPI
751 mpi::all_reduce( this->
comm(), local_sum, global_sum, std::plus<value_type>() );
763 FEELPP_DONT_INLINE this_type
sqrt()
const;
769 this_type
pow(
int n )
const;
783 M_vec.operator[]( i ) += a;
822 FEELPP_ASSERT( 0 ).error(
"invalid call, not implemented yet" );
829 void localize ( ublas::vector<value_type>& v_local )
const;
835 void localize ( ublas::vector_range<ublas::vector<value_type> >& v_local )
const;
841 void localize ( ublas::vector_slice<ublas::vector<value_type> >& v_local )
const;
855 const std::vector<size_type>& send_list )
const;
863 const std::vector<size_type>& send_list );
896 void checkInvariant()
const;
900 mutable bool M_global_values_updated;
911 template <
typename T>
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" );
920 typedef typename type_traits<T>::real_type real_type;
927 _t.operator()( start+i ) = v1.operator()( start + i )* v2.operator()( start + i );
939 template <
typename T>
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> > >;