30 #ifndef __VectorPetsc_H
31 #define __VectorPetsc_H 1
33 #include <feel/feelconfig.h>
35 #include <feel/feelalg/vector.hpp>
38 #if defined(FEELPP_HAS_PETSC_H)
51 template<
typename T>
class MatrixPetsc;
66 class VectorPetsc :
public Vector<T>
68 typedef Vector<T> super;
72 friend class boost::serialization::access;
78 typedef typename super::value_type value_type;
79 typedef typename super::real_type real_type;
80 typedef typename super::clone_ptrtype clone_ptrtype;
81 typedef typename super::datamap_type datamap_type;
82 typedef typename super::datamap_ptrtype datamap_ptrtype;
96 M_destroy_vec_on_exit( true )
103 VectorPetsc (
const size_type n, WorldComm
const& _worldComm = Environment::worldComm() )
105 super( n, _worldComm ),
106 M_destroy_vec_on_exit( true )
108 this->init( n, n,
false );
117 WorldComm
const& _worldComm = Environment::worldComm() )
119 super( n, n_local, _worldComm ),
120 M_destroy_vec_on_exit( true )
122 this->init( n, n_local,
false );
125 VectorPetsc ( datamap_ptrtype
const& dm,
bool doInit=
true )
128 M_destroy_vec_on_exit( true )
131 this->init( dm->nDof(), dm->nLocalDofWithoutGhost(), false );
145 M_destroy_vec_on_exit( false )
148 this->M_is_initialized =
true;
151 VectorPetsc( Vec v, datamap_ptrtype
const& dm )
154 M_destroy_vec_on_exit( false )
157 this->M_is_initialized =
true;
164 VectorPetsc( VectorPetsc<value_type> &v, IS &is )
167 M_destroy_vec_on_exit( false )
169 #if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2)
173 datamap_ptrtype dm(
new datamap_type(n, n, v.comm()) );
176 VecGetSubVector(v.vec(), is, &this->M_vec);
177 this->M_is_initialized =
true;
183 VectorPetsc( Vector<value_type>
const& v, std::vector<int>
const& index )
187 M_destroy_vec_on_exit( false )
189 #if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2)
191 VectorPetsc<T>
const* V =
dynamic_cast<VectorPetsc<T> const*
> ( &v );
195 int n = index.size();
196 PetscMalloc(n*
sizeof(PetscInt),&map);
197 for (
int i=0; i<n; i++) map[i] = index[i];
198 ierr = ISCreateGeneral(Environment::worldComm(),n,map,PETSC_COPY_VALUES,&is);
199 CHKERRABORT( this->comm(),ierr );
202 datamap_ptrtype dm(
new datamap_type(n, n, V->comm()) );
205 ierr = VecGetSubVector(V->vec(), is, &this->M_vec);
206 CHKERRABORT( this->comm(),ierr );
207 this->M_is_initialized =
true;
226 clone_ptrtype clone ()
const;
243 const bool fast=
false );
249 const bool fast=
false )
251 this->init( N,N,fast );
260 value_type operator() (
const size_type i )
const;
267 Vector<T> & operator += (
const Vector<value_type> &V )
269 FEELPP_ASSERT( this->closed() ).error(
"vector is not closed" );
280 Vector<T> & operator -= (
const Vector<value_type> &V )
282 FEELPP_ASSERT( this->closed() ).error(
"vector is not closed" );
295 bool destroy_vec_on_exit()
const
297 return M_destroy_vec_on_exit;
309 DCHECK( this->isInitialized() ) <<
"VectorPetsc not initialized";
312 if ( !this->isInitialized() )
316 int ierr = VecGetSize( M_vec, &petsc_size );
317 CHKERRABORT( this->comm(),ierr );
318 return static_cast<size_type>( petsc_size );
326 DCHECK( this->isInitialized() ) <<
"VectorPetsc not initialized";
329 int ierr = VecGetLocalSize( M_vec, &petsc_size );
330 CHKERRABORT( this->comm(),ierr );
332 return static_cast<size_type>( petsc_size );
342 FEELPP_ASSERT ( M_vec != 0 ).error(
"invalid petsc vector" );
347 FEELPP_ASSERT ( M_vec != 0 ).error(
"invalid petsc vector" );
369 FEELPP_ASSERT ( this->isInitialized() ).error(
"VectorPetsc<> not initialized" );
373 ierr = VecAssemblyBegin( M_vec );
374 CHKERRABORT( this->comm(),ierr );
375 ierr = VecAssemblyEnd( M_vec );
376 CHKERRABORT( this->comm(),ierr );
378 this->M_is_closed =
true;
398 void setConstant( value_type v )
406 FEELPP_DONT_INLINE
void clear ();
408 void localize(
const Vector<T>& V);
413 void set (
const value_type& value );
418 void set (
size_type i,
const value_type& value );
423 void add (
size_type i,
const value_type& value );
428 void addVector (
int* i,
int n, value_type* v );
435 void addVector (
const std::vector<value_type>& v,
436 const std::vector<size_type>& dof_indices )
438 FEELPP_ASSERT ( v.size() == dof_indices.size() ).error(
"invalid dof indices" );
441 this->add ( dof_indices[i], v[i] );
450 void addVector (
const Vector<value_type>& V,
451 const std::vector<size_type>& dof_indices )
453 FEELPP_ASSERT ( V.size() == dof_indices.size() ).error(
"invalid dof indices" );
456 this->add ( dof_indices[i], V( i ) );
464 void addVector (
const Vector<value_type>& V_in,
465 const MatrixSparse<value_type>& A_in );
474 void addVector (
const ublas::vector<value_type>& V,
475 const std::vector<size_type>& dof_indices )
477 FEELPP_ASSERT ( V.size() == dof_indices.size() ).error(
"invalid dof indices" );
480 this->add ( dof_indices[i], V( i ) );
487 void insert (
const std::vector<T>& ,
488 const std::vector<size_type>& )
490 FEELPP_ASSERT( 0 ).error(
"invalid call, not implemented yet" );
499 void insert (
const Vector<T>& V,
500 const std::vector<size_type>& dof_indices );
509 void insert (
const ublas::vector<T>& V,
510 const std::vector<size_type>& dof_indices );
516 void scale (
const T factor );
524 void add (
const value_type& v_in );
531 void add (
const Vector<value_type>& v );
538 void add (
const value_type& a_in,
const Vector<value_type>& v_in );
545 real_type min ()
const;
552 real_type max()
const;
558 real_type l1Norm ()
const;
565 real_type l2Norm ()
const;
571 value_type sum ()
const;
578 real_type linftyNorm ()
const;
586 assert ( this->isInitialized() );
588 int ierr=0, petsc_first=0, petsc_last=0;
590 ierr = VecGetOwnershipRange ( M_vec, &petsc_first, &petsc_last );
591 CHKERRABORT( this->comm(),ierr );
593 return static_cast<size_type>( petsc_first );
604 assert ( this->isInitialized() );
606 int ierr=0, petsc_first=0, petsc_last=0;
608 ierr = VecGetOwnershipRange ( M_vec, &petsc_first, &petsc_last );
609 CHKERRABORT( this->comm(),ierr );
611 return static_cast<size_type>( petsc_last );
620 void printMatlab(
const std::string name=
"NULL",
bool renumber =
false )
const;
622 value_type dot( Vector<T>
const& __v );
627 template<
class Archive>
628 void serialize(Archive & ar,
const unsigned int version)
const
630 value_type *array = *((value_type**)this->vec()->data);
631 int n = this->vec()->map->n;
632 int N = this->vec()->map->N;
635 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
637 FEELPP_ASSERT( n==N ).error(
"wrong vector type for serialization (!=VECSEQ)" );
639 for (
int i=0; i<n; i++)
651 VectorPetsc( VectorPetsc
const & v )
654 M_destroy_vec_on_exit( true )
656 FEELPP_ASSERT( v.closed() ).error(
"copied vector is not closed" );
658 VecDuplicate( v.M_vec, &M_vec );
659 VecCopy( v.M_vec, M_vec );
660 this->M_is_initialized =
true;
676 bool M_destroy_vec_on_exit;
687 class VectorPetscMPI :
public VectorPetsc<T>
689 typedef VectorPetsc<T> super;
690 typedef typename super::value_type value_type;
691 typedef typename super::datamap_type datamap_type;
692 typedef typename super::datamap_ptrtype datamap_ptrtype;
693 typedef typename super::clone_ptrtype clone_ptrtype;
701 VectorPetscMPI( Vec v, datamap_ptrtype
const& dm );
703 VectorPetscMPI( datamap_ptrtype
const& dm );
709 clone_ptrtype clone ()
const;
712 const bool fast=
false );
714 value_type operator() (
const size_type i )
const;
716 void set(
size_type i,
const value_type& value );
718 void add(
const size_type i,
const value_type& value );
720 void addVector(
int* i,
int n, value_type* v );
731 void duplicateFromOtherPartition( Vector<T>
const& vecInput );
733 value_type dot( Vector<T>
const& __v );
739 void duplicateFromOtherPartition_run( Vector<T>
const& vecInput );
743 VecScatter M_vecScatter;
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319