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
vectorepetra.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):
6  Christophe Prud'homme <christophe.prudhomme@feelpp.org>
7  Klaus Sapelza <klaus.sapelza@epfl.ch>
8  Date: 2006-09-14
9 
10  Copyright (C) 2006,2007 EPFL
11  Copyright (C) 2006,2007,2008 Université Joseph Fourier (Grenoble I)
12 
13  This library is free software; you can redistribute it and/or
14  modify it under the terms of the GNU Lesser General Public
15  License as published by the Free Software Foundation; either
16  version 3.0 of the License, or (at your option) any later version.
17 
18  This library is distributed in the hope that it will be useful,
19  but WITHOUT ANY WARRANTY; without even the implied warranty of
20  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
21  Lesser General Public License for more details.
22 
23  You should have received a copy of the GNU Lesser General Public
24  License along with this library; if not, write to the Free Software
25  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
26 */
33 #ifndef __VectorEpetra_H
34 #define __VectorEpetra_H 1
35 
36 #include <feel/feelconfig.h>
37 
38 #include <feel/feelalg/vector.hpp>
41 #if defined(FEELPP_HAS_TRILINOS_EPETRA)
42 
43 #if defined(FEELPP_HAS_MPI)
44 #include <Epetra_MpiComm.h>
45 #else
46 #include <Epetra_SerialComm.h>
47 #endif /* FEELPP_HAS_MPI */
48 
49 #include <EpetraExt_MultiVectorOut.h>
50 #include <Epetra_FEVector.h>
51 #include <Epetra_Vector.h>
52 
53 
54 
55 namespace Feel
56 {
68 template<typename T>
69 class VectorEpetra : public Vector<T>
70 {
71  typedef Vector<T> super;
72 
73 
74 public:
75 
76 
80 
81  typedef typename super::value_type value_type;
82  typedef typename super::real_type real_type;
83  typedef typename super::clone_ptrtype clone_ptrtype;
84 
85  typedef VectorEpetra<value_type> epetra_vector_type;
86  typedef boost::shared_ptr<epetra_vector_type> epetra_vector_ptrtype;
88 
92 
93  VectorEpetra ();
94 
98  VectorEpetra ( Epetra_BlockMap const& emap );
99 
104  VectorEpetra ( Epetra_FEVector const * v );
105 
110  VectorEpetra ( Epetra_Vector const * v );
111 
116  VectorEpetra ( VectorEpetra const& v );
117 
122  ~VectorEpetra ();
123 
127  Epetra_BlockMap Map() const
128  {
129  return M_vec.Map();
130  }
131 
132 
137  clone_ptrtype clone () const
138  {
139  clone_ptrtype cloned_vector ( new VectorEpetra<T> );
140 
141  *cloned_vector = *this;
142 
143  return cloned_vector;
144  }
145 
146 
159  void init ( Epetra_BlockMap const& emap, const bool fast = false );
160 
161  void init ( const size_type N, const size_type n_local, const bool fast = false );
162 
166  //operator Epetra_Vector() { return M_vec; }
167 
168 
169 
170  value_type operator() ( const size_type i ) const
171  {
172  checkInvariants();
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" );
176 
177  value_type value=0.;
178  int ierr=0, dummy;
179  double* values;
180 
181  ierr = M_vec.ExtractView( &values, &dummy );
182 
183  value = values[i - this->firstLocalIndex()];
184  return static_cast<value_type>( value );
185 
186  }
187 
188  value_type& operator() ( const size_type i )
189  {
190  checkInvariants();
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" );
194 
195  return M_vec[0][ ( int )i-this->firstLocalIndex() ];
196 
197  }
198 
199 
200  VectorEpetra& operator=( VectorEpetra const& v )
201  {
202  if ( &v != this )
203  {
204  super::operator=( v );
205  M_emap = v.Map();
206  M_vec = v.vec();
207 
208  this->M_is_initialized = true;
209  }
210 
211  checkInvariants();
212  return *this;
213  }
214 
219  Vector<T> & operator += ( const Vector<T> &V )
220  {
221 
222  this->add( 1., V );
223 
224  return *this;
225  }
226 
231  super & operator -= ( const super &V )
232  {
233 
234  this->add( -1., V );
235 
236  return *this;
237  }
238 
239 
241 
245 
253  size_type size () const
254  {
255  FEELPP_ASSERT ( this->isInitialized() ).error( "VectorEpetra not initialized" );
256 
257 
258  if ( !this->isInitialized() )
259  return 0;
260 
261  int epetra_size=0;
262  epetra_size = M_vec.GlobalLength();
263  return static_cast<size_type>( epetra_size );
264 
265  return 0;
266  }
271  size_type localSize() const
272  {
273  FEELPP_ASSERT ( this->isInitialized() ).error( "VectorEpetra not initialized" );
274 
275  int epetra_size=0;
276  epetra_size = M_vec.MyLength();
277  DVLOG(2) << "[VectorEpetra::localSize] localSize= " << epetra_size << "\n";
278  return static_cast<size_type>( epetra_size );
279  }
280 
286  Epetra_FEVector const& vec () const
287  {
288  //FEELPP_ASSERT (M_vec != 0).error( "invalid epetra vector" );
289  return M_vec;
290  }
291 
297  Epetra_FEVector& vec ()
298  {
299  //FEELPP_ASSERT (M_vec != 0).error( "invalid epetra vector" );
300  return M_vec;
301  }
302 
310  boost::shared_ptr<Epetra_Vector> epetraVector ()
311  {
312  double** V;
313  M_vec.ExtractView( &V );
314  boost::shared_ptr<Epetra_Vector> EV( new Epetra_Vector( View,M_vec.Map(),V[0] ) );
315  return EV;
316  }
317 
319 
323 
324 
326 
330 
334  void close ()
335  {
336  FEELPP_ASSERT ( this->isInitialized() ).error( "VectorEpetra<> not initialized" );
337 
338  int ierr=0;
339  ierr = M_vec.GlobalAssemble( Add );
340 
341  this->M_is_closed = true;
342  }
343 
348  void zero ()
349  {
350  FEELPP_ASSERT ( this->isInitialized() ).error( "VectorEpetra<> not initialized" );
351 
352  M_vec.PutScalar( 0.0 );
353 
354  }
355 
359  void zero ( size_type /*start*/, size_type /*stop*/ )
360  {
361  this->zero();
362  }
363 
367  void setConstant( value_type v )
368  {
369  this->set( v );
370  }
371 
375  void clear ()
376  {
377  if ( this->isInitialized() ) //&& (this->M_destroy_vec_on_exit))
378  {
379  M_emap = Epetra_BlockMap( -1, 0, 0, M_vec.Comm() );
380  M_vec.ReplaceMap( M_emap );
381  M_vec.PutScalar( 0.0 );
382  }
383 
384  this->M_is_closed = this->M_is_initialized = false;
385  }
386 
390  void set ( const value_type& value );
391 
395  void set ( size_type i, const value_type& value );
396 
400  void add ( size_type i, const value_type& value );
401 
405  void addVector ( int* i, int n, value_type* v );
406 
407 
408  void addVector ( VectorEpetra& v )
409  {
410  M_vec.Update( 1.0, v.vec(), 1.0 );
411  }
412 
413 
418  void addVector ( const Vector<T>& _v,
419  const MatrixSparse<T>& _M );
420 
425  void addVector ( const std::vector<value_type>& v,
426  const std::vector<size_type>& dof_indices )
427  {
428  FEELPP_ASSERT ( v.size() == dof_indices.size() ).error( "invalid dof indices" );
429 
430  for ( size_type i=0; i<v.size(); i++ )
431  this->add ( dof_indices[i], v[i] );
432  }
433 
434 
441  void addVector ( const Vector<value_type>& V,
442  const std::vector<size_type>& dof_indices )
443  {
444  FEELPP_ASSERT ( V.size() == dof_indices.size() ).error( "invalid dof indices" );
445 
446  for ( size_type i=0; i<V.size(); i++ )
447  this->add ( dof_indices[i], V( i ) );
448  }
449  //
450  //
451  // /**
452  // * \f$ U+=A*V\f$, add the product of a \p MatrixSparse \p A
453  // * and a \p Vector \p V to \p this, where \p this=U.
454  // */
455  // void addVector (const Vector<value_type>& V_in,
456  // const MatrixSparse<value_type>& A_in)
457  // {
458  //
459  //
460  // }
461  //
462 
469  void addVector ( const ublas::vector<value_type>& V,
470  const std::vector<size_type>& dof_indices )
471  {
472  FEELPP_ASSERT ( V.size() == dof_indices.size() ).error( "invalid dof indices" );
473 
474  for ( size_type i=0; i<V.size(); i++ )
475  this->add ( dof_indices[i], V( i ) );
476  }
477 
483  void add ( const value_type& v_in )
484  {
485  value_type v = static_cast<value_type>( v_in );
486  const int n = static_cast<int>( this->size() );
487 
488  for ( int i=0 ; i<n ; i++ )
489  {
490  M_vec.SumIntoGlobalValues( 1,&i,&v ); //indices are in global index space
491  //M_vec.SumIntoMyValues(1,&v,&i); //indices are in local index space
492  }
493 
494  }
495 
501  void add ( const Vector<value_type>& v )
502  {
503  VectorEpetra* v_ptr = const_cast<VectorEpetra*>( dynamic_cast< VectorEpetra const*>( &v ) );
504  this->add ( 1., *v_ptr );
505  }
506 
512  void add ( const value_type& a_in, const Vector<value_type>& v_in )
513  {
514 
515  const value_type a = static_cast<value_type>( a_in );
516  const VectorEpetra<T>* v = dynamic_cast<const VectorEpetra<T>*>( &v_in );
517 
518 
519  assert ( v != NULL );
520  assert( this->size() == v->size() );
521 
522  M_vec.Update( a, v->M_vec,1. );
523 
524  }
525 
530  void insert ( const std::vector<T>& v,
531  const std::vector<size_type>& dof_indices )
532  {
533  //to be implemented
534  }
535 
542  void insert ( const Vector<T>& V,
543  const std::vector<size_type>& dof_indices )
544  {
545  //to be implemented
546  }
547 
554  void insert ( const ublas::vector<T>& /*V*/,
555  const std::vector<size_type>& /*dof_indices*/ )
556  {
557  //to be implemented
558  }
559 
564  void scale ( const T /*factor*/ )
565  {
566  //to be implemented
567  }
568 
569 
574  void localize ( std::vector<T>& /*v_local*/ ) const
575  {
576  //to be implemented
577  }
578 
583  void localize ( Vector<T>& /*v_local*/ ) const
584  {
585  //to be implemented
586  }
587 
593  void localize ( Vector<T>& /*v_local*/,
594  const std::vector<size_type>& /*send_list*/ ) const
595  {
596  }
597 
602  void localize ( const size_type /*first_local_idx*/,
603  const size_type /*last_local_idx*/,
604  const std::vector<size_type>& /*send_list*/ )
605  {
606  //to be implemented
607  }
614  void localizeToOneProcessor ( std::vector<T>& /*v_local*/,
615  const size_type /*proc_id*/=0 ) const
616  {
617  //to be implemented
618  }
619 
625  real_type min () const
626  {
627  assert ( this->isInitialized() );
628 
629  double min=0.;
630 
631  M_vec.MinValue( &min );
632 
633  // this return value is correct: VecMin returns a PetscReal
634  return static_cast<double>( min );
635  }
636 
642  real_type max() const
643  {
644  assert ( this->isInitialized() );
645 
646  double max=0.;
647 
648  M_vec.MaxValue( &max );
649 
650  // this return value is correct: VecMin returns a PetscReal
651  return static_cast<double>( max );
652  }
653 
658  real_type l1Norm () const
659  {
660  //assert(this->closed());
661 
662  double value=0.;
663 
664  M_vec.Norm1( &value );
665 
666  return static_cast<Real>( value );
667  }
668 
674  real_type l2Norm () const
675  {
676  //assert(this->closed());
677 
678  double value=0.;
679 
680  M_vec.Norm2( &value );
681 
682  return static_cast<Real>( value );
683 
684  }
685 
691  real_type linftyNorm () const
692  {
693  assert( this->closed() );
694 
695  double value=0.;
696 
697  M_vec.NormInf( &value );
698 
699  return static_cast<Real>( value );
700 
701  }
702 
706  value_type sum() const
707  {
708  //assert(this->closed());
709 
710  double value=0.;
711  double global_sum=0;
712 
713  double const * pointers( M_vec[0] );
714 
715  for ( int i( 0 ); i < M_vec.MyLength(); ++i , ++pointers )
716  value += *pointers;
717 
718  M_vec.Comm().SumAll( &value, &global_sum, 1 );
719 
720  return static_cast<Real>( global_sum );
721 
722 
723  }
724 
725 
732  size_type firstLocalIndex () const
733  {
734 
735  int epetra_first = 0;
736  assert ( this->isInitialized() );
737  //epetra_first = M_vec.Map().MinMyGID();
738  epetra_first = M_vec.Map().MinLID();
739  DVLOG(2) << "[VectorEpetra::firstLocalIndex] firstLocalIndex= " << epetra_first << "\n";
740  return static_cast<size_type>( epetra_first );
741  }
742 
743 
744 
751  size_type lastLocalIndex () const
752  {
753  int epetra_last = 0;
754  assert ( this->isInitialized() );
755  //epetra_last = M_vec.Map().MaxMyGID();
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;
759  }
760 
766  void printMatlab ( const std::string name = "", bool renumber = false ) const;
767 
768  // @}
769 
770 
771 
772 protected:
773 
774 private:
775  void checkInvariants() const;
776 private:
777  Epetra_BlockMap M_emap;
778  Epetra_FEVector M_vec;
779 
784  //const bool M_destroy_vec_on_exit;
785 };
786 
787 template<typename T>
788 DebugStream&
789 operator<<( DebugStream& __os, VectorEpetra<T> const& __n );
790 
791 
792 template<typename T>
793 NdebugStream&
794 operator<<( NdebugStream& __os, VectorEpetra<T> const& __n );
795 
796 DebugStream&
797 operator<<( DebugStream& __os, Epetra_BlockMap const& __n );
798 
799 
800 NdebugStream&
801 operator<<( NdebugStream& __os, Epetra_BlockMap const& __n );
802 
803 
804 } // Feel
805 #endif /* FEELPP_HAS_EPETRA */
806 #endif /* __VectorEpetra_H */
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
Elements & operator=(Elements const &e)
Definition: elements.hpp:335

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