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
vectorpetsc.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: 2005-10-18
7 
8  Copyright (C) 2005,2006 EPFL
9  Copyright (C) 2007 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 __VectorPetsc_H
31 #define __VectorPetsc_H 1
32 
33 #include <feel/feelconfig.h>
34 
35 #include <feel/feelalg/vector.hpp>
37 
38 #if defined(FEELPP_HAS_PETSC_H)
40 
41 
42 extern "C"
43 {
44 #include <petscmat.h>
45 }
46 
47 
48 
49 namespace Feel
50 {
51 template<typename T> class MatrixPetsc;
52 
65 template<typename T>
66 class VectorPetsc : public Vector<T>
67 {
68  typedef Vector<T> super;
69 
70 public:
71 
72  friend class boost::serialization::access;
73 
77 
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;
83 
85 
89 
93  VectorPetsc ()
94  :
95  super(),
96  M_destroy_vec_on_exit( true )
97  {
98  }
99 
103  VectorPetsc ( const size_type n, WorldComm const& _worldComm = Environment::worldComm() )
104  :
105  super( n, _worldComm ),
106  M_destroy_vec_on_exit( true )
107  {
108  this->init( n, n, false );
109  }
110 
115  VectorPetsc ( const size_type n,
116  const size_type n_local,
117  WorldComm const& _worldComm = Environment::worldComm() )
118  :
119  super( n, n_local, _worldComm ),
120  M_destroy_vec_on_exit( true )
121  {
122  this->init( n, n_local, false );
123  }
124 
125  VectorPetsc ( datamap_ptrtype const& dm, bool doInit=true )
126  :
127  super( dm ),
128  M_destroy_vec_on_exit( true )
129  {
130  if ( doInit )
131  this->init( dm->nDof(), dm->nLocalDofWithoutGhost(), false );
132  }
133 
134 
142  VectorPetsc( Vec v )
143  :
144  super(),
145  M_destroy_vec_on_exit( false )
146  {
147  this->M_vec = v;
148  this->M_is_initialized = true;
149  }
150 
151  VectorPetsc( Vec v, datamap_ptrtype const& dm )
152  :
153  super( dm ),
154  M_destroy_vec_on_exit( false )
155  {
156  this->M_vec = v;
157  this->M_is_initialized = true;
158  }
159 
164  VectorPetsc( VectorPetsc<value_type> &v, IS &is )
165  :
166  super(),
167  M_destroy_vec_on_exit( false )
168  {
169 #if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2)
170  /* map */
171  PetscInt n;
172  ISGetSize(is,&n);
173  datamap_ptrtype dm( new datamap_type(n, n, v.comm()) );
174  this->setMap(dm);
175  /* init */
176  VecGetSubVector(v.vec(), is, &this->M_vec);
177  this->M_is_initialized = true;
178  /* close */
179  this->close(); /* no // assembly required */
180 #endif
181  }
182 
183  VectorPetsc( Vector<value_type> const& v, std::vector<int> const& index )
184  :
185  super(),
186  //super(v,index),
187  M_destroy_vec_on_exit( false )
188  {
189 #if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2)
190 
191  VectorPetsc<T> const* V = dynamic_cast<VectorPetsc<T> const*> ( &v );
192  int ierr=0;
193  IS is;
194  PetscInt *map;
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 );
200  PetscFree(map);
201 
202  datamap_ptrtype dm( new datamap_type(n, n, V->comm()) );
203  this->setMap(dm);
204  /* init */
205  ierr = VecGetSubVector(V->vec(), is, &this->M_vec);
206  CHKERRABORT( this->comm(),ierr );
207  this->M_is_initialized = true;
208  /* close */
209  this->close(); /* no // assembly required */
210 #endif
211  }
212 
217  ~VectorPetsc ()
218  {
219  this->clear();
220  }
221 
226  clone_ptrtype clone () const;
227 
228 
241  void init ( const size_type N,
242  const size_type n_local,
243  const bool fast=false );
244 
248  void init ( const size_type N,
249  const bool fast=false )
250  {
251  this->init( N,N,fast );
252  }
253 
255 
259 
260  value_type operator() ( const size_type i ) const;
261 
262 
267  Vector<T> & operator += ( const Vector<value_type> &V )
268  {
269  FEELPP_ASSERT( this->closed() ).error( "vector is not closed" );
270 
271  this->add( 1., V );
272 
273  return *this;
274  }
275 
280  Vector<T> & operator -= ( const Vector<value_type> &V )
281  {
282  FEELPP_ASSERT( this->closed() ).error( "vector is not closed" );
283 
284  this->add( -1., V );
285 
286  return *this;
287  }
288 
290 
294 
295  bool destroy_vec_on_exit() const
296  {
297  return M_destroy_vec_on_exit;
298  }
299 
307  size_type size () const
308  {
309  DCHECK( this->isInitialized() ) << "VectorPetsc not initialized";
310 
311 
312  if ( !this->isInitialized() )
313  return 0;
314 
315  int petsc_size=0;
316  int ierr = VecGetSize( M_vec, &petsc_size );
317  CHKERRABORT( this->comm(),ierr );
318  return static_cast<size_type>( petsc_size );
319  }
324  size_type localSize() const
325  {
326  DCHECK( this->isInitialized() ) << "VectorPetsc not initialized";
327 
328  int petsc_size=0;
329  int ierr = VecGetLocalSize( M_vec, &petsc_size );
330  CHKERRABORT( this->comm(),ierr );
331 
332  return static_cast<size_type>( petsc_size );
333  }
334 
340  Vec vec () const
341  {
342  FEELPP_ASSERT ( M_vec != 0 ).error( "invalid petsc vector" );
343  return M_vec;
344  }
345  Vec& vec ()
346  {
347  FEELPP_ASSERT ( M_vec != 0 ).error( "invalid petsc vector" );
348  return M_vec;
349  }
350 
352 
356 
357 
359 
363 
367  void close ()
368  {
369  FEELPP_ASSERT ( this->isInitialized() ).error( "VectorPetsc<> not initialized" );
370 
371  int ierr=0;
372 
373  ierr = VecAssemblyBegin( M_vec );
374  CHKERRABORT( this->comm(),ierr );
375  ierr = VecAssemblyEnd( M_vec );
376  CHKERRABORT( this->comm(),ierr );
377 
378  this->M_is_closed = true;
379  }
380 
385  void zero ();
386 
390  void zero ( size_type /*start*/, size_type /*stop*/ )
391  {
392  this->zero();
393  }
394 
398  void setConstant( value_type v )
399  {
400  this->set( v );
401  }
402 
406  FEELPP_DONT_INLINE void clear ();
407 
408  void localize( const Vector<T>& V);
409 
413  void set ( const value_type& value );
414 
418  void set ( size_type i, const value_type& value );
419 
423  void add ( size_type i, const value_type& value );
424 
428  void addVector ( int* i, int n, value_type* v );
429 
435  void addVector ( const std::vector<value_type>& v,
436  const std::vector<size_type>& dof_indices )
437  {
438  FEELPP_ASSERT ( v.size() == dof_indices.size() ).error( "invalid dof indices" );
439 
440  for ( size_type i=0; i<v.size(); i++ )
441  this->add ( dof_indices[i], v[i] );
442  }
443 
450  void addVector ( const Vector<value_type>& V,
451  const std::vector<size_type>& dof_indices )
452  {
453  FEELPP_ASSERT ( V.size() == dof_indices.size() ).error( "invalid dof indices" );
454 
455  for ( size_type i=0; i<V.size(); i++ )
456  this->add ( dof_indices[i], V( i ) );
457  }
458 
459 
464  void addVector ( const Vector<value_type>& V_in,
465  const MatrixSparse<value_type>& A_in );
466 
467 
474  void addVector ( const ublas::vector<value_type>& V,
475  const std::vector<size_type>& dof_indices )
476  {
477  FEELPP_ASSERT ( V.size() == dof_indices.size() ).error( "invalid dof indices" );
478 
479  for ( size_type i=0; i<V.size(); i++ )
480  this->add ( dof_indices[i], V( i ) );
481  }
482 
487  void insert ( const std::vector<T>& /*v*/,
488  const std::vector<size_type>& /*dof_indices*/ )
489  {
490  FEELPP_ASSERT( 0 ).error( "invalid call, not implemented yet" );
491  }
492 
499  void insert ( const Vector<T>& V,
500  const std::vector<size_type>& dof_indices );
501 
502 
509  void insert ( const ublas::vector<T>& V,
510  const std::vector<size_type>& dof_indices );
511 
516  void scale ( const T factor );
517 
518 
524  void add ( const value_type& v_in );
525 
531  void add ( const Vector<value_type>& v );
532 
538  void add ( const value_type& a_in, const Vector<value_type>& v_in );
539 
545  real_type min () const;
546 
552  real_type max() const;
553 
558  real_type l1Norm () const;
559 
565  real_type l2Norm () const;
566 
571  value_type sum () const;
572 
578  real_type linftyNorm () const;
579 
584  size_type firstLocalIndex () const
585  {
586  assert ( this->isInitialized() );
587 
588  int ierr=0, petsc_first=0, petsc_last=0;
589 
590  ierr = VecGetOwnershipRange ( M_vec, &petsc_first, &petsc_last );
591  CHKERRABORT( this->comm(),ierr );
592 
593  return static_cast<size_type>( petsc_first );
594  }
595 
596 
597 
602  size_type lastLocalIndex () const
603  {
604  assert ( this->isInitialized() );
605 
606  int ierr=0, petsc_first=0, petsc_last=0;
607 
608  ierr = VecGetOwnershipRange ( M_vec, &petsc_first, &petsc_last );
609  CHKERRABORT( this->comm(),ierr );
610 
611  return static_cast<size_type>( petsc_last );
612  }
613 
614 
620  void printMatlab( const std::string name="NULL", bool renumber = false ) const;
621 
622  value_type dot( Vector<T> const& __v );
623 
627  template<class Archive>
628  void serialize(Archive & ar, const unsigned int version) const
629  {
630  value_type *array = *((value_type**)this->vec()->data);
631  int n = this->vec()->map->n;
632  int N = this->vec()->map->N;
633 
634  int rank;
635  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
636 
637  FEELPP_ASSERT( n==N ).error( "wrong vector type for serialization (!=VECSEQ)" );
638 
639  for (int i=0; i<n; i++)
640  ar & array[i];
641  }
643 
644 
645 
646 protected:
647 
648 public:
649 
650  // disable
651  VectorPetsc( VectorPetsc const & v )
652  :
653  super( v ),
654  M_destroy_vec_on_exit( true )
655  {
656  FEELPP_ASSERT( v.closed() ).error( "copied vector is not closed" );
657 
658  VecDuplicate( v.M_vec, &M_vec );
659  VecCopy( v.M_vec, M_vec );
660  this->M_is_initialized = true;
661  this->close();
662  }
663 
664 
665 protected:
666 
670  Vec M_vec;
671 
676  bool M_destroy_vec_on_exit;
677 };
678 
679 
680 //----------------------------------------------------------------------------------------------------//
681 //----------------------------------------------------------------------------------------------------//
682 //----------------------------------------------------------------------------------------------------//
683 //----------------------------------------------------------------------------------------------------//
684 //----------------------------------------------------------------------------------------------------//
685 
686 template<typename T>
687 class VectorPetscMPI : public VectorPetsc<T>
688 {
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;
694 public:
695 
696  VectorPetscMPI()
697  :
698  super()
699  {}
700 
701  VectorPetscMPI( Vec v, datamap_ptrtype const& dm );
702 
703  VectorPetscMPI( datamap_ptrtype const& dm );
704 
705  ~VectorPetscMPI()
706  {
707  this->clear();
708  }
709  clone_ptrtype clone () const;
710  void init( const size_type N,
711  const size_type n_local,
712  const bool fast=false );
713 
714  value_type operator() ( const size_type i ) const;
715 
716  void set( size_type i, const value_type& value );
717 
718  void add( const size_type i, const value_type& value );
719 
720  void addVector( int* i, int n, value_type* v );
721 
722  void clear();
723 
724  void localize();
725 
726  void close();
727 
728  size_type firstLocalIndex() const;
729  size_type lastLocalIndex() const;
730 
731  void duplicateFromOtherPartition( Vector<T> const& vecInput );
732 
733  value_type dot( Vector<T> const& __v );
734 
735  size_type localSize() const;
736 
737 private :
738 
739  void duplicateFromOtherPartition_run( Vector<T> const& vecInput );
740 
741  Vec M_vecLocal;
742 
743  VecScatter M_vecScatter;
744 
745 };
746 
747 } // Feel
748 #endif /* FEELPP_HAS_PETSC */
749 #endif /* __VectorPetsc_H */
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319

Generated on Sun Dec 22 2013 13:11:15 for Feel++ by doxygen 1.8.5