33 #ifndef __VectorEpetra_H
34 #define __VectorEpetra_H 1
36 #include <feel/feelconfig.h>
38 #include <feel/feelalg/vector.hpp>
41 #if defined(FEELPP_HAS_TRILINOS_EPETRA)
43 #if defined(FEELPP_HAS_MPI)
44 #include <Epetra_MpiComm.h>
46 #include <Epetra_SerialComm.h>
49 #include <EpetraExt_MultiVectorOut.h>
50 #include <Epetra_FEVector.h>
51 #include <Epetra_Vector.h>
69 class VectorEpetra :
public Vector<T>
71 typedef Vector<T> super;
81 typedef typename super::value_type value_type;
82 typedef typename super::real_type real_type;
83 typedef typename super::clone_ptrtype clone_ptrtype;
85 typedef VectorEpetra<value_type> epetra_vector_type;
86 typedef boost::shared_ptr<epetra_vector_type> epetra_vector_ptrtype;
98 VectorEpetra ( Epetra_BlockMap
const& emap );
104 VectorEpetra ( Epetra_FEVector
const * v );
110 VectorEpetra ( Epetra_Vector
const * v );
116 VectorEpetra ( VectorEpetra
const& v );
127 Epetra_BlockMap Map()
const
137 clone_ptrtype clone ()
const
139 clone_ptrtype cloned_vector (
new VectorEpetra<T> );
141 *cloned_vector = *
this;
143 return cloned_vector;
159 void init ( Epetra_BlockMap
const& emap,
const bool fast =
false );
170 value_type operator() (
const size_type i )
const
173 FEELPP_ASSERT ( this->isInitialized() ).error(
"vector not initialized" );
174 FEELPP_ASSERT ( ( ( i >= this->firstLocalIndex() ) &&
175 ( i < this->lastLocalIndex() ) ) )( i )( this->firstLocalIndex() )( this->lastLocalIndex() ).warn(
"invalid vector index" );
181 ierr = M_vec.ExtractView( &values, &dummy );
183 value = values[i - this->firstLocalIndex()];
184 return static_cast<value_type
>( value );
188 value_type& operator() (
const size_type i )
191 FEELPP_ASSERT ( this->isInitialized() ).error(
"vector not initialized" );
192 FEELPP_ASSERT ( ( ( i >= this->firstLocalIndex() ) &&
193 ( i < this->lastLocalIndex() ) ) )( i )( this->firstLocalIndex() )( this->lastLocalIndex() ).warn(
"invalid vector index" );
195 return M_vec[0][ ( int )i-this->firstLocalIndex() ];
200 VectorEpetra&
operator=( VectorEpetra
const& v )
208 this->M_is_initialized =
true;
219 Vector<T> & operator += (
const Vector<T> &V )
231 super & operator -= (
const super &V )
255 FEELPP_ASSERT ( this->isInitialized() ).error(
"VectorEpetra not initialized" );
258 if ( !this->isInitialized() )
262 epetra_size = M_vec.GlobalLength();
263 return static_cast<size_type>( epetra_size );
273 FEELPP_ASSERT ( this->isInitialized() ).error(
"VectorEpetra not initialized" );
276 epetra_size = M_vec.MyLength();
277 DVLOG(2) <<
"[VectorEpetra::localSize] localSize= " << epetra_size <<
"\n";
278 return static_cast<size_type>( epetra_size );
286 Epetra_FEVector
const& vec ()
const
297 Epetra_FEVector& vec ()
310 boost::shared_ptr<Epetra_Vector> epetraVector ()
313 M_vec.ExtractView( &V );
314 boost::shared_ptr<Epetra_Vector> EV(
new Epetra_Vector( View,M_vec.Map(),V[0] ) );
336 FEELPP_ASSERT ( this->isInitialized() ).error(
"VectorEpetra<> not initialized" );
339 ierr = M_vec.GlobalAssemble( Add );
341 this->M_is_closed =
true;
350 FEELPP_ASSERT ( this->isInitialized() ).error(
"VectorEpetra<> not initialized" );
352 M_vec.PutScalar( 0.0 );
367 void setConstant( value_type v )
377 if ( this->isInitialized() )
379 M_emap = Epetra_BlockMap( -1, 0, 0, M_vec.Comm() );
380 M_vec.ReplaceMap( M_emap );
381 M_vec.PutScalar( 0.0 );
384 this->M_is_closed = this->M_is_initialized =
false;
390 void set (
const value_type& value );
395 void set (
size_type i,
const value_type& value );
400 void add (
size_type i,
const value_type& value );
405 void addVector (
int* i,
int n, value_type* v );
408 void addVector ( VectorEpetra& v )
410 M_vec.Update( 1.0, v.vec(), 1.0 );
418 void addVector (
const Vector<T>& _v,
419 const MatrixSparse<T>& _M );
425 void addVector (
const std::vector<value_type>& v,
426 const std::vector<size_type>& dof_indices )
428 FEELPP_ASSERT ( v.size() == dof_indices.size() ).error(
"invalid dof indices" );
431 this->add ( dof_indices[i], v[i] );
441 void addVector (
const Vector<value_type>& V,
442 const std::vector<size_type>& dof_indices )
444 FEELPP_ASSERT ( V.size() == dof_indices.size() ).error(
"invalid dof indices" );
447 this->add ( dof_indices[i], V( i ) );
469 void addVector (
const ublas::vector<value_type>& V,
470 const std::vector<size_type>& dof_indices )
472 FEELPP_ASSERT ( V.size() == dof_indices.size() ).error(
"invalid dof indices" );
475 this->add ( dof_indices[i], V( i ) );
483 void add (
const value_type& v_in )
485 value_type v =
static_cast<value_type
>( v_in );
486 const int n =
static_cast<int>( this->size() );
488 for (
int i=0 ; i<n ; i++ )
490 M_vec.SumIntoGlobalValues( 1,&i,&v );
501 void add (
const Vector<value_type>& v )
503 VectorEpetra* v_ptr =
const_cast<VectorEpetra*
>(
dynamic_cast< VectorEpetra const*
>( &v ) );
504 this->add ( 1., *v_ptr );
512 void add (
const value_type& a_in,
const Vector<value_type>& v_in )
515 const value_type a =
static_cast<value_type
>( a_in );
516 const VectorEpetra<T>* v =
dynamic_cast<const VectorEpetra<T>*
>( &v_in );
519 assert ( v != NULL );
520 assert( this->size() == v->size() );
522 M_vec.Update( a, v->M_vec,1. );
530 void insert (
const std::vector<T>& v,
531 const std::vector<size_type>& dof_indices )
542 void insert (
const Vector<T>& V,
543 const std::vector<size_type>& dof_indices )
554 void insert (
const ublas::vector<T>& ,
555 const std::vector<size_type>& )
564 void scale (
const T )
574 void localize ( std::vector<T>& )
const
583 void localize ( Vector<T>& )
const
593 void localize ( Vector<T>& ,
594 const std::vector<size_type>& )
const
604 const std::vector<size_type>& )
614 void localizeToOneProcessor ( std::vector<T>& ,
625 real_type min ()
const
627 assert ( this->isInitialized() );
631 M_vec.MinValue( &min );
634 return static_cast<double>( min );
642 real_type max()
const
644 assert ( this->isInitialized() );
648 M_vec.MaxValue( &max );
651 return static_cast<double>( max );
658 real_type l1Norm ()
const
664 M_vec.Norm1( &value );
666 return static_cast<Real
>( value );
674 real_type l2Norm ()
const
680 M_vec.Norm2( &value );
682 return static_cast<Real
>( value );
691 real_type linftyNorm ()
const
693 assert( this->closed() );
697 M_vec.NormInf( &value );
699 return static_cast<Real
>( value );
706 value_type sum()
const
713 double const * pointers( M_vec[0] );
715 for (
int i( 0 ); i < M_vec.MyLength(); ++i , ++pointers )
718 M_vec.Comm().SumAll( &value, &global_sum, 1 );
720 return static_cast<Real
>( global_sum );
735 int epetra_first = 0;
736 assert ( this->isInitialized() );
738 epetra_first = M_vec.Map().MinLID();
739 DVLOG(2) <<
"[VectorEpetra::firstLocalIndex] firstLocalIndex= " << epetra_first <<
"\n";
740 return static_cast<size_type>( epetra_first );
754 assert ( this->isInitialized() );
756 epetra_last = M_vec.Map().MaxLID();
757 DVLOG(2) <<
"[VectorEpetra::lastLocalIndex] lastLocalIndex= " << epetra_last+1 <<
"\n";
758 return static_cast<size_type>( epetra_last )+1;
766 void printMatlab (
const std::string name =
"",
bool renumber =
false )
const;
775 void checkInvariants()
const;
777 Epetra_BlockMap M_emap;
778 Epetra_FEVector M_vec;
789 operator<<( DebugStream& __os, VectorEpetra<T>
const& __n );
794 operator<<( NdebugStream& __os, VectorEpetra<T>
const& __n );
797 operator<<( DebugStream& __os, Epetra_BlockMap
const& __n );
801 operator<<( NdebugStream& __os, Epetra_BlockMap
const& __n );