libMesh
trilinos_epetra_vector.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 #ifndef LIBMESH_TRILINOS_EPETRA_VECTOR_H
19 #define LIBMESH_TRILINOS_EPETRA_VECTOR_H
20 
21 
22 #include "libmesh/libmesh_common.h"
23 
24 #ifdef LIBMESH_TRILINOS_HAVE_EPETRA
25 
26 // Local includes
27 #include "libmesh/numeric_vector.h"
28 #include "libmesh/parallel.h"
29 
30 // Trilinos includes
31 #include "libmesh/ignore_warnings.h"
32 #include <Epetra_CombineMode.h>
33 #include <Epetra_Map.h>
34 #include <Epetra_MultiVector.h>
35 #include <Epetra_Vector.h>
36 #include <Epetra_MpiComm.h>
37 #include "libmesh/restore_warnings.h"
38 
39 // C++ includes
40 #include <cstddef>
41 #include <vector>
42 
43 // Forward declarations
44 class Epetra_IntSerialDenseVector;
45 class Epetra_SerialDenseVector;
46 
47 namespace libMesh
48 {
49 
50 // forward declarations
51 template <typename T> class SparseMatrix;
52 
61 template <typename T>
62 class EpetraVector libmesh_final : public NumericVector<T>
63 {
64 public:
65 
69  explicit
70  EpetraVector (const Parallel::Communicator & comm,
71  const ParallelType type = AUTOMATIC);
72 
76  explicit
77  EpetraVector (const Parallel::Communicator & comm,
78  const numeric_index_type n,
79  const ParallelType type = AUTOMATIC);
80 
85  EpetraVector (const Parallel::Communicator & comm,
86  const numeric_index_type n,
87  const numeric_index_type n_local,
88  const ParallelType type = AUTOMATIC);
89 
95  EpetraVector (const Parallel::Communicator & comm,
96  const numeric_index_type N,
97  const numeric_index_type n_local,
98  const std::vector<numeric_index_type> & ghost,
99  const ParallelType type = AUTOMATIC);
100 
108  EpetraVector(Epetra_Vector & v,
109  const Parallel::Communicator & comm
110  LIBMESH_CAN_DEFAULT_TO_COMMWORLD);
111 
116  ~EpetraVector ();
117 
118  virtual void close () libmesh_override;
119 
120  virtual void clear () libmesh_override;
121 
122  virtual void zero () libmesh_override;
123 
124  virtual UniquePtr<NumericVector<T>> zero_clone () const libmesh_override;
125 
126  virtual UniquePtr<NumericVector<T>> clone () const libmesh_override;
127 
128  virtual void init (const numeric_index_type N,
129  const numeric_index_type n_local,
130  const bool fast=false,
131  const ParallelType type=AUTOMATIC) libmesh_override;
132 
133  virtual void init (const numeric_index_type N,
134  const bool fast=false,
135  const ParallelType type=AUTOMATIC) libmesh_override;
136 
137  virtual void init (const numeric_index_type N,
138  const numeric_index_type n_local,
139  const std::vector<numeric_index_type> & ghost,
140  const bool fast = false,
141  const ParallelType = AUTOMATIC) libmesh_override;
142 
143  virtual void init (const NumericVector<T> & other,
144  const bool fast = false) libmesh_override;
145 
146  virtual NumericVector<T> & operator= (const T s) libmesh_override;
147 
148  virtual NumericVector<T> & operator= (const NumericVector<T> & v) libmesh_override;
149 
155  EpetraVector<T> & operator= (const EpetraVector<T> & v);
156 
157  virtual NumericVector<T> & operator= (const std::vector<T> & v) libmesh_override;
158 
159  virtual Real min () const libmesh_override;
160 
161  virtual Real max () const libmesh_override;
162 
163  virtual T sum () const libmesh_override;
164 
165  virtual Real l1_norm () const libmesh_override;
166 
167  virtual Real l2_norm () const libmesh_override;
168 
169  virtual Real linfty_norm () const libmesh_override;
170 
171  virtual numeric_index_type size () const libmesh_override;
172 
173  virtual numeric_index_type local_size() const libmesh_override;
174 
175  virtual numeric_index_type first_local_index() const libmesh_override;
176 
177  virtual numeric_index_type last_local_index() const libmesh_override;
178 
179  virtual T operator() (const numeric_index_type i) const libmesh_override;
180 
181  virtual NumericVector<T> & operator += (const NumericVector<T> & v) libmesh_override;
182 
183  virtual NumericVector<T> & operator -= (const NumericVector<T> & v) libmesh_override;
184 
185  virtual NumericVector<T> & operator /= (NumericVector<T> & v) libmesh_override;
186 
187  virtual void reciprocal() libmesh_override;
188 
189  virtual void conjugate() libmesh_override;
190 
191  virtual void set (const numeric_index_type i, const T value) libmesh_override;
192 
193  virtual void add (const numeric_index_type i, const T value) libmesh_override;
194 
195  virtual void add (const T s) libmesh_override;
196 
197  virtual void add (const NumericVector<T> & v) libmesh_override;
198 
199  virtual void add (const T a, const NumericVector<T> & v) libmesh_override;
200 
205  using NumericVector<T>::add_vector;
206 
207  virtual void add_vector (const T * v,
208  const std::vector<numeric_index_type> & dof_indices) libmesh_override;
209 
210  virtual void add_vector (const NumericVector<T> & v,
211  const SparseMatrix<T> & A) libmesh_override;
212 
213  virtual void add_vector_transpose (const NumericVector<T> & v,
214  const SparseMatrix<T> & A) libmesh_override;
215 
220  using NumericVector<T>::insert;
221 
222  virtual void insert (const T * v,
223  const std::vector<numeric_index_type> & dof_indices) libmesh_override;
224 
225  virtual void scale (const T factor) libmesh_override;
226 
227  virtual void abs() libmesh_override;
228 
229  virtual T dot(const NumericVector<T> & v) const libmesh_override;
230 
231  virtual void localize (std::vector<T> & v_local) const libmesh_override;
232 
233  virtual void localize (NumericVector<T> & v_local) const libmesh_override;
234 
235  virtual void localize (NumericVector<T> & v_local,
236  const std::vector<numeric_index_type> & send_list) const libmesh_override;
237 
238  virtual void localize (std::vector<T> & v_local,
239  const std::vector<numeric_index_type> & indices) const libmesh_override;
240 
241  virtual void localize (const numeric_index_type first_local_idx,
242  const numeric_index_type last_local_idx,
243  const std::vector<numeric_index_type> & send_list) libmesh_override;
244 
245  virtual void localize_to_one (std::vector<T> & v_local,
246  const processor_id_type proc_id=0) const libmesh_override;
247 
248  virtual void pointwise_mult (const NumericVector<T> & vec1,
249  const NumericVector<T> & vec2) libmesh_override;
250 
251  virtual void create_subvector (NumericVector<T> & subvector,
252  const std::vector<numeric_index_type> & rows) const libmesh_override;
253 
254  virtual void swap (NumericVector<T> & v) libmesh_override;
255 
264  Epetra_Vector * vec () { libmesh_assert(_vec); return _vec; }
265 
266 private:
267 
271  Epetra_Vector * _vec;
272 
277 
282  bool _destroy_vec_on_exit;
283 
284  // The following were copied (and slightly modified) from
285  // Epetra_FEVector.h in order to allow us to use a standard
286  // Epetra_Vector... which is more compatible with other Trilinos
287  // packages such as NOX. All of this code is originally under LGPL
288 
293  int SumIntoGlobalValues(int numIDs,
294  const int * GIDs,
295  const double * values);
296 
307  int SumIntoGlobalValues(const Epetra_IntSerialDenseVector & GIDs,
308  const Epetra_SerialDenseVector & values);
309 
314  int ReplaceGlobalValues(int numIDs,
315  const int * GIDs,
316  const double * values);
317 
328  int ReplaceGlobalValues(const Epetra_IntSerialDenseVector & GIDs,
329  const Epetra_SerialDenseVector & values);
330 
331  int SumIntoGlobalValues(int numIDs,
332  const int * GIDs,
333  const int * numValuesPerID,
334  const double * values);
335 
336  int ReplaceGlobalValues(int numIDs,
337  const int * GIDs,
338  const int * numValuesPerID,
339  const double * values);
340 
349  int GlobalAssemble(Epetra_CombineMode mode = Add);
350 
354  void setIgnoreNonLocalEntries(bool flag) {
355  ignoreNonLocalEntries_ = flag;
356  }
357 
358  void FEoperatorequals(const EpetraVector & source);
359 
360  int inputValues(int numIDs,
361  const int * GIDs,
362  const double * values,
363  bool accumulate);
364 
365  int inputValues(int numIDs,
366  const int * GIDs,
367  const int * numValuesPerID,
368  const double * values,
369  bool accumulate);
370 
371  int inputNonlocalValue(int GID,
372  double value,
373  bool accumulate);
374 
375  int inputNonlocalValues(int GID,
376  int numValues,
377  const double * values,
378  bool accumulate);
379 
380  void destroyNonlocalData();
381 
384  double * myCoefs_;
385 
390  double ** nonlocalCoefs_;
391 
397  unsigned char last_edit;
398 
400 };
401 
402 
403 /*----------------------- Inline functions ----------------------------------*/
404 
405 
406 
407 template <typename T>
408 inline
409 EpetraVector<T>::EpetraVector (const Parallel::Communicator & comm,
410  const ParallelType type) :
411  NumericVector<T>(comm, type),
412  _destroy_vec_on_exit(true),
413  myFirstID_(0),
414  myNumIDs_(0),
415  myCoefs_(libmesh_nullptr),
416  nonlocalIDs_(libmesh_nullptr),
417  nonlocalElementSize_(libmesh_nullptr),
418  numNonlocalIDs_(0),
419  allocatedNonlocalLength_(0),
420  nonlocalCoefs_(libmesh_nullptr),
421  last_edit(0),
422  ignoreNonLocalEntries_(false)
423 {
424  this->_type = type;
425 }
426 
427 
428 
429 template <typename T>
430 inline
431 EpetraVector<T>::EpetraVector (const Parallel::Communicator & comm,
432  const numeric_index_type n,
433  const ParallelType type) :
435  _destroy_vec_on_exit(true),
436  myFirstID_(0),
437  myNumIDs_(0),
441  numNonlocalIDs_(0),
444  last_edit(0),
446 
447 {
448  this->init(n, n, false, type);
449 }
450 
451 
452 
453 template <typename T>
454 inline
455 EpetraVector<T>::EpetraVector (const Parallel::Communicator & comm,
456  const numeric_index_type n,
457  const numeric_index_type n_local,
458  const ParallelType type) :
460  _destroy_vec_on_exit(true),
461  myFirstID_(0),
462  myNumIDs_(0),
466  numNonlocalIDs_(0),
469  last_edit(0),
471 {
472  this->init(n, n_local, false, type);
473 }
474 
475 
476 
477 
478 template <typename T>
479 inline
480 EpetraVector<T>::EpetraVector(Epetra_Vector & v,
481  const Parallel::Communicator & comm) :
483  _destroy_vec_on_exit(false),
484  myFirstID_(0),
485  myNumIDs_(0),
489  numNonlocalIDs_(0),
492  last_edit(0),
494 {
495  _vec = &v;
496 
497  this->_type = PARALLEL; // FIXME - need to determine this from v!
498 
499  myFirstID_ = _vec->Map().MinMyGID();
500  myNumIDs_ = _vec->Map().NumMyElements();
501 
502  _map.reset(new Epetra_Map(_vec->GlobalLength(),
503  _vec->MyLength(),
504  0, // IndexBase = 0 for C/C++, 1 for Fortran.
505  Epetra_MpiComm (this->comm().get())));
506 
507  //Currently we impose the restriction that NumVectors==1, so we won't
508  //need the LDA argument when calling ExtractView. Hence the "dummy" arg.
509  int dummy;
510  _vec->ExtractView(&myCoefs_, &dummy);
511 
512  this->_is_closed = true;
513  this->_is_initialized = true;
514 }
515 
516 
517 
518 template <typename T>
519 inline
520 EpetraVector<T>::EpetraVector (const Parallel::Communicator & comm,
521  const numeric_index_type n,
522  const numeric_index_type n_local,
523  const std::vector<numeric_index_type> & ghost,
524  const ParallelType type) :
526  _destroy_vec_on_exit(true),
527  myFirstID_(0),
528  myNumIDs_(0),
532  numNonlocalIDs_(0),
535  last_edit(0),
537 {
538  this->init(n, n_local, ghost, false, type);
539 }
540 
541 
542 
543 // Default implementation for solver packages for which ghosted
544 // vectors are not yet implemented.
545 template <class T>
546 void EpetraVector<T>::init (const NumericVector<T> & other,
547  const bool fast)
548 {
549  this->init(other.size(),other.local_size(),fast,other.type());
550 }
551 
552 
553 
554 template <typename T>
555 inline
556 EpetraVector<T>::~EpetraVector ()
557 {
558  this->clear ();
559 }
560 
561 
562 
563 template <typename T>
564 inline
566  const numeric_index_type n_local,
567  const bool fast,
568  const ParallelType type)
569 {
570  // We default to allocating n_local local storage
571  numeric_index_type my_n_local = n_local;
572 
573  if (type == AUTOMATIC)
574  {
575  if (n == n_local)
576  this->_type = SERIAL;
577  else
578  this->_type = PARALLEL;
579  }
580  else if (type == GHOSTED)
581  {
582  // We don't yet support GHOSTED Epetra vectors, so to get the
583  // same functionality we need a SERIAL vector with local
584  // storage allocated for every entry.
585  this->_type = SERIAL;
586  my_n_local = n;
587  }
588  else
589  this->_type = type;
590 
591  libmesh_assert ((this->_type==SERIAL && n==my_n_local) ||
592  this->_type==PARALLEL);
593 
594  _map.reset(new Epetra_Map(static_cast<int>(n),
595  my_n_local,
596  0,
597  Epetra_MpiComm (this->comm().get())));
598 
599  _vec = new Epetra_Vector(*_map);
600 
601  myFirstID_ = _vec->Map().MinMyGID();
602  myNumIDs_ = _vec->Map().NumMyElements();
603 
604  // Currently we impose the restriction that NumVectors==1, so we won't
605  // need the LDA argument when calling ExtractView. Hence the "dummy" arg.
606  int dummy;
607  _vec->ExtractView(&myCoefs_, &dummy);
608 
609  this->_is_initialized = true;
610  this->_is_closed = true;
611  this->last_edit = 0;
612 
613  if (fast == false)
614  this->zero ();
615 }
616 
617 
618 template <typename T>
619 inline
621  const numeric_index_type n_local,
622  const std::vector<numeric_index_type> & /*ghost*/,
623  const bool fast,
624  const ParallelType type)
625 {
626  // TODO: we shouldn't ignore the ghost sparsity pattern
627  this->init(n, n_local, fast, type);
628 }
629 
630 
631 
632 template <typename T>
633 inline
635  const bool fast,
636  const ParallelType type)
637 {
638  this->init(n,n,fast,type);
639 }
640 
641 
642 
643 template <typename T>
644 inline
645 void EpetraVector<T>::close ()
646 {
647  libmesh_assert (this->initialized());
648 
649  // Are we adding or inserting?
650  unsigned char global_last_edit = last_edit;
651  this->comm().max(global_last_edit);
652  libmesh_assert(!last_edit || last_edit == global_last_edit);
653 
654  if (global_last_edit == 1)
655  this->GlobalAssemble(Insert);
656  else if (global_last_edit == 2)
657  this->GlobalAssemble(Add);
658  else
659  libmesh_assert(!global_last_edit);
660 
661  this->_is_closed = true;
662  this->last_edit = 0;
663 }
664 
665 
666 
667 template <typename T>
668 inline
669 void EpetraVector<T>::clear ()
670 {
671  if (this->initialized())
672  {
673  // We might just be an interface to a user-provided _vec
674  if (this->_destroy_vec_on_exit)
675  {
676  delete _vec;
678  }
679 
680  // But we currently always own our own _map
681  _map.reset();
682  }
683 
684  this->_is_closed = this->_is_initialized = false;
685 }
686 
687 
688 
689 template <typename T>
690 inline
691 void EpetraVector<T>::zero ()
692 {
693  libmesh_assert (this->initialized());
694  libmesh_assert (this->closed());
695 
696  _vec->PutScalar(0.0);
697 }
698 
699 
700 
701 template <typename T>
702 inline
703 UniquePtr<NumericVector<T>> EpetraVector<T>::zero_clone () const
704 {
705  NumericVector<T> * cloned_vector = new EpetraVector<T>(this->comm(), AUTOMATIC);
706  cloned_vector->init(*this);
707  return UniquePtr<NumericVector<T>>(cloned_vector);
708 }
709 
710 
711 
712 template <typename T>
713 inline
714 UniquePtr<NumericVector<T>> EpetraVector<T>::clone () const
715 {
716  NumericVector<T> * cloned_vector = new EpetraVector<T>(this->comm(), AUTOMATIC);
717  cloned_vector->init(*this, true);
718  *cloned_vector = *this;
719  return UniquePtr<NumericVector<T>>(cloned_vector);
720 }
721 
722 
723 
724 template <typename T>
725 inline
726 numeric_index_type EpetraVector<T>::size () const
727 {
728  libmesh_assert (this->initialized());
729 
730  return _vec->GlobalLength();
731 }
732 
733 
734 
735 template <typename T>
736 inline
737 numeric_index_type EpetraVector<T>::local_size () const
738 {
739  libmesh_assert (this->initialized());
740 
741  return _vec->MyLength();
742 }
743 
744 template <typename T>
745 inline
746 numeric_index_type EpetraVector<T>::first_local_index () const
747 {
748  libmesh_assert (this->initialized());
749 
750  return _vec->Map().MinMyGID();
751 }
752 
753 
754 
755 template <typename T>
756 inline
757 numeric_index_type EpetraVector<T>::last_local_index () const
758 {
759  libmesh_assert (this->initialized());
760 
761  return _vec->Map().MaxMyGID()+1;
762 }
763 
764 
765 template <typename T>
766 inline
767 T EpetraVector<T>::operator() (const numeric_index_type i) const
768 {
769  libmesh_assert (this->initialized());
770  libmesh_assert ( ((i >= this->first_local_index()) &&
771  (i < this->last_local_index())) );
772 
773  return (*_vec)[i-this->first_local_index()];
774 }
775 
776 
777 
778 template <typename T>
779 inline
780 Real EpetraVector<T>::min () const
781 {
782  libmesh_assert (this->initialized());
783 
784  T value;
785 
786  _vec->MinValue(&value);
787 
788  return value;
789 }
790 
791 
792 
793 template <typename T>
794 inline
796 {
797  libmesh_assert (this->initialized());
798 
799  T value;
800 
801  _vec->MaxValue(&value);
802 
803  return value;
804 }
805 
806 
807 
808 template <typename T>
809 inline
811 {
812  NumericVector<T>::swap(other);
813 
814  EpetraVector<T> & v = cast_ref<EpetraVector<T> &>(other);
815 
816  std::swap(_vec, v._vec);
817  _map.swap(v._map);
818  std::swap(_destroy_vec_on_exit, v._destroy_vec_on_exit);
819  std::swap(myFirstID_, v.myFirstID_);
820  std::swap(myNumIDs_, v.myNumIDs_);
821  std::swap(myCoefs_, v.myCoefs_);
822  std::swap(nonlocalIDs_, v.nonlocalIDs_);
823  std::swap(nonlocalElementSize_, v.nonlocalElementSize_);
824  std::swap(numNonlocalIDs_, v.numNonlocalIDs_);
825  std::swap(allocatedNonlocalLength_, v.allocatedNonlocalLength_);
826  std::swap(nonlocalCoefs_, v.nonlocalCoefs_);
827  std::swap(last_edit, v.last_edit);
828  std::swap(ignoreNonLocalEntries_, v.ignoreNonLocalEntries_);
829 }
830 
831 } // namespace libMesh
832 
833 
834 #endif // #ifdef LIBMESH_TRILINOS_HAVE_EPETRA
835 #endif // LIBMESH_TRILINOS_EPETRA_VECTOR_H
virtual numeric_index_type n() const libmesh_override
Epetra_Vector * _vec
Actual Epetra vector datatype to hold vector entries.
double abs(double a)
Encapsulates the MPI_Comm object.
Definition: parallel.h:657
virtual numeric_index_type size() const =0
virtual void init(const numeric_index_type, const numeric_index_type, const bool=false, const ParallelType=AUTOMATIC)=0
Change the dimension of the vector to N.
void max(T &r) const
Take a local variable and replace it with the maximum of it&#39;s values on all processors.
void scale(MeshBase &mesh, const Real xs, const Real ys=0., const Real zs=0.)
Scales the mesh.
uint8_t processor_id_type
Definition: id_types.h:99
const class libmesh_nullptr_t libmesh_nullptr
Numeric vector.
Definition: dof_map.h:66
bool _is_initialized
true once init() has been called.
virtual bool initialized() const
unsigned char last_edit
Keep track of whether the last write operation on this vector was nothing (0) or a sum (1) or an add ...
The libMesh namespace provides an interface to certain functionality in the library.
void sum(T &r, const Communicator &comm=Communicator_World)
const Number zero
.
Definition: libmesh.h:178
long double max(long double a, double b)
libmesh_assert(j)
virtual bool closed() const libmesh_override
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
UniquePtr< Epetra_Map > _map
Holds the distributed Map.
DataType _vec
Actual Eigen::SparseVector<> we are wrapping.
ParallelType
Defines an enum for parallel data structure types.
dof_id_type numeric_index_type
Definition: id_types.h:92
bool _destroy_vec_on_exit
This boolean value should only be set to false for the constructor which takes a PETSc Vec object...
Definition: petsc_vector.h:430
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
Epetra_Map * _map
Holds the distributed Map.
void setIgnoreNonLocalEntries(bool flag)
Set whether or not non-local data values should be ignored.
virtual ElemType type() const libmesh_override
Definition: cell_hex20.h:81
virtual void init() libmesh_override
Initialize this matrix using the sparsity structure computed by dof_map.
ParallelType _type
Type of vector.
virtual numeric_index_type local_size() const =0
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void clear() libmesh_override
Restores the NumericVector<T> to a pristine state.
virtual void swap(NumericVector< T > &v)
Swaps the contents of this with v.
X_input swap(X_system)
static PetscErrorCode Mat * A
int GlobalAssemble(Epetra_CombineMode mode=Add)
Gather any overlapping/shared data into the non-overlapping partitioning defined by the Map that was ...
const Parallel::Communicator & comm() const
void swap(Iterator &lhs, Iterator &rhs)
swap, used to implement op=
virtual void zero() libmesh_override
Set all entries to zero.
virtual numeric_index_type last_local_index() const libmesh_override
static const bool value
Definition: xdr_io.C:108
ParallelType type() const
virtual numeric_index_type first_local_index() const libmesh_override
long double min(long double a, double b)
bool _is_closed
Flag which tracks whether the vector&#39;s values are consistent on all processors after insertion or add...