trilinos_epetra_vector.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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
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>
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 final : public NumericVector<T>
63 {
64 public:
65 
69  explicit
71  const ParallelType type = AUTOMATIC);
72 
76  explicit
78  const numeric_index_type n,
79  const ParallelType type = AUTOMATIC);
80 
86  const numeric_index_type n,
87  const numeric_index_type n_local,
88  const ParallelType type = AUTOMATIC);
89 
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 
116  EpetraVector (EpetraVector &&) = delete;
117  EpetraVector (const EpetraVector &) = delete;
118  EpetraVector & operator= (const EpetraVector &) = delete;
119  EpetraVector & operator= (EpetraVector &&) = delete;
120  virtual ~EpetraVector ();
121 
122  virtual void close () override;
123 
124  virtual void clear () override;
125 
126  virtual void zero () override;
127 
128  virtual std::unique_ptr<NumericVector<T>> zero_clone () const override;
129 
130  virtual std::unique_ptr<NumericVector<T>> clone () const override;
131 
132  virtual void init (const numeric_index_type N,
133  const numeric_index_type n_local,
134  const bool fast=false,
135  const ParallelType type=AUTOMATIC) override;
136 
137  virtual void init (const numeric_index_type N,
138  const bool fast=false,
139  const ParallelType type=AUTOMATIC) override;
140 
141  virtual void init (const numeric_index_type N,
142  const numeric_index_type n_local,
143  const std::vector<numeric_index_type> & ghost,
144  const bool fast = false,
145  const ParallelType = AUTOMATIC) override;
146 
147  virtual void init (const NumericVector<T> & other,
148  const bool fast = false) override;
149 
150  virtual NumericVector<T> & operator= (const T s) override;
151 
152  virtual NumericVector<T> & operator= (const NumericVector<T> & v) override;
153 
154  virtual NumericVector<T> & operator= (const std::vector<T> & v) override;
155 
156  virtual Real min () const override;
157 
158  virtual Real max () const override;
159 
160  virtual T sum () const override;
161 
162  virtual Real l1_norm () const override;
163 
164  virtual Real l2_norm () const override;
165 
166  virtual Real linfty_norm () const override;
167 
168  virtual numeric_index_type size () const override;
169 
170  virtual numeric_index_type local_size() const override;
171 
172  virtual numeric_index_type first_local_index() const override;
173 
174  virtual numeric_index_type last_local_index() const override;
175 
176  virtual T operator() (const numeric_index_type i) const override;
177 
178  virtual NumericVector<T> & operator += (const NumericVector<T> & v) override;
179 
180  virtual NumericVector<T> & operator -= (const NumericVector<T> & v) override;
181 
182  virtual NumericVector<T> & operator /= (const NumericVector<T> & v) override;
183 
184  virtual void reciprocal() override;
185 
186  virtual void conjugate() override;
187 
188  virtual void set (const numeric_index_type i, const T value) override;
189 
190  virtual void add (const numeric_index_type i, const T value) override;
191 
192  virtual void add (const T s) override;
193 
194  virtual void add (const NumericVector<T> & v) override;
195 
196  virtual void add (const T a, const NumericVector<T> & v) override;
197 
203 
204  virtual void add_vector (const T * v,
205  const std::vector<numeric_index_type> & dof_indices) override;
206 
207  virtual void add_vector (const NumericVector<T> & v,
208  const SparseMatrix<T> & A) override;
209 
210  virtual void add_vector_transpose (const NumericVector<T> & v,
211  const SparseMatrix<T> & A) override;
212 
218 
219  virtual void insert (const T * v,
220  const std::vector<numeric_index_type> & dof_indices) override;
221 
222  virtual void scale (const T factor) override;
223 
224  virtual void abs() override;
225 
226  virtual T dot(const NumericVector<T> & v) const override;
227 
228  virtual void localize (std::vector<T> & v_local) const override;
229 
230  virtual void localize (NumericVector<T> & v_local) const override;
231 
232  virtual void localize (NumericVector<T> & v_local,
233  const std::vector<numeric_index_type> & send_list) const override;
234 
235  virtual void localize (std::vector<T> & v_local,
236  const std::vector<numeric_index_type> & indices) const override;
237 
238  virtual void localize (const numeric_index_type first_local_idx,
239  const numeric_index_type last_local_idx,
240  const std::vector<numeric_index_type> & send_list) override;
241 
242  virtual void localize_to_one (std::vector<T> & v_local,
243  const processor_id_type proc_id=0) const override;
244 
245  virtual void pointwise_mult (const NumericVector<T> & vec1,
246  const NumericVector<T> & vec2) override;
247 
248  virtual void create_subvector (NumericVector<T> & subvector,
249  const std::vector<numeric_index_type> & rows) const override;
250 
251  virtual void swap (NumericVector<T> & v) override;
252 
261  Epetra_Vector * vec () { libmesh_assert(_vec); return _vec; }
262 
263 private:
264 
268  Epetra_Vector * _vec;
269 
273  std::unique_ptr<Epetra_Map> _map;
274 
280 
281  // The following were copied (and slightly modified) from
282  // Epetra_FEVector.h in order to allow us to use a standard
283  // Epetra_Vector... which is more compatible with other Trilinos
284  // packages such as NOX. All of this code is originally under LGPL
285 
290  int SumIntoGlobalValues(int numIDs,
291  const int * GIDs,
292  const double * values);
293 
304  int SumIntoGlobalValues(const Epetra_IntSerialDenseVector & GIDs,
305  const Epetra_SerialDenseVector & values);
306 
311  int ReplaceGlobalValues(int numIDs,
312  const int * GIDs,
313  const double * values);
314 
325  int ReplaceGlobalValues(const Epetra_IntSerialDenseVector & GIDs,
326  const Epetra_SerialDenseVector & values);
327 
328  int SumIntoGlobalValues(int numIDs,
329  const int * GIDs,
330  const int * numValuesPerID,
331  const double * values);
332 
333  int ReplaceGlobalValues(int numIDs,
334  const int * GIDs,
335  const int * numValuesPerID,
336  const double * values);
337 
346  int GlobalAssemble(Epetra_CombineMode mode = Add);
347 
351  void setIgnoreNonLocalEntries(bool flag) {
352  ignoreNonLocalEntries_ = flag;
353  }
354 
355  void FEoperatorequals(const EpetraVector & source);
356 
357  int inputValues(int numIDs,
358  const int * GIDs,
359  const double * values,
360  bool accumulate);
361 
362  int inputValues(int numIDs,
363  const int * GIDs,
364  const int * numValuesPerID,
365  const double * values,
366  bool accumulate);
367 
368  int inputNonlocalValue(int GID,
369  double value,
370  bool accumulate);
371 
372  int inputNonlocalValues(int GID,
373  int numValues,
374  const double * values,
375  bool accumulate);
376 
377  void destroyNonlocalData();
378 
381  double * myCoefs_;
382 
387  double ** nonlocalCoefs_;
388 
394  unsigned char last_edit;
395 
397 };
398 
399 
400 /*----------------------- Inline functions ----------------------------------*/
401 
402 
403 
404 template <typename T>
405 inline
407  const ParallelType type) :
408  NumericVector<T>(comm, type),
409  _destroy_vec_on_exit(true),
410  myFirstID_(0),
411  myNumIDs_(0),
412  myCoefs_(nullptr),
413  nonlocalIDs_(nullptr),
414  nonlocalElementSize_(nullptr),
415  numNonlocalIDs_(0),
416  allocatedNonlocalLength_(0),
417  nonlocalCoefs_(nullptr),
418  last_edit(0),
419  ignoreNonLocalEntries_(false)
420 {
421  this->_type = type;
422 }
423 
424 
425 
426 template <typename T>
427 inline
429  const numeric_index_type n,
430  const ParallelType type) :
431  NumericVector<T>(comm, type),
432  _destroy_vec_on_exit(true),
433  myFirstID_(0),
434  myNumIDs_(0),
435  myCoefs_(nullptr),
436  nonlocalIDs_(nullptr),
437  nonlocalElementSize_(nullptr),
438  numNonlocalIDs_(0),
439  allocatedNonlocalLength_(0),
440  nonlocalCoefs_(nullptr),
441  last_edit(0),
442  ignoreNonLocalEntries_(false)
443 
444 {
445  this->init(n, n, false, type);
446 }
447 
448 
449 
450 template <typename T>
451 inline
453  const numeric_index_type n,
454  const numeric_index_type n_local,
455  const ParallelType type) :
456  NumericVector<T>(comm, type),
457  _destroy_vec_on_exit(true),
458  myFirstID_(0),
459  myNumIDs_(0),
460  myCoefs_(nullptr),
461  nonlocalIDs_(nullptr),
462  nonlocalElementSize_(nullptr),
463  numNonlocalIDs_(0),
464  allocatedNonlocalLength_(0),
465  nonlocalCoefs_(nullptr),
466  last_edit(0),
467  ignoreNonLocalEntries_(false)
468 {
469  this->init(n, n_local, false, type);
470 }
471 
472 
473 
474 
475 template <typename T>
476 inline
478  const Parallel::Communicator & comm) :
479  NumericVector<T>(comm, AUTOMATIC),
480  _destroy_vec_on_exit(false),
481  myFirstID_(0),
482  myNumIDs_(0),
483  myCoefs_(nullptr),
484  nonlocalIDs_(nullptr),
485  nonlocalElementSize_(nullptr),
486  numNonlocalIDs_(0),
487  allocatedNonlocalLength_(0),
488  nonlocalCoefs_(nullptr),
489  last_edit(0),
490  ignoreNonLocalEntries_(false)
491 {
492  _vec = &v;
493 
494  this->_type = PARALLEL; // FIXME - need to determine this from v!
495 
496  myFirstID_ = _vec->Map().MinMyGID();
497  myNumIDs_ = _vec->Map().NumMyElements();
498 
499  _map.reset(new Epetra_Map(_vec->GlobalLength(),
500  _vec->MyLength(),
501  0, // IndexBase = 0 for C/C++, 1 for Fortran.
502  Epetra_MpiComm (this->comm().get())));
503 
504  //Currently we impose the restriction that NumVectors==1, so we won't
505  //need the LDA argument when calling ExtractView. Hence the "dummy" arg.
506  int dummy;
507  _vec->ExtractView(&myCoefs_, &dummy);
508 
509  this->_is_closed = true;
510  this->_is_initialized = true;
511 }
512 
513 
514 
515 template <typename T>
516 inline
518  const numeric_index_type n,
519  const numeric_index_type n_local,
520  const std::vector<numeric_index_type> & ghost,
521  const ParallelType type) :
522  NumericVector<T>(comm, AUTOMATIC),
523  _destroy_vec_on_exit(true),
524  myFirstID_(0),
525  myNumIDs_(0),
526  myCoefs_(nullptr),
527  nonlocalIDs_(nullptr),
528  nonlocalElementSize_(nullptr),
529  numNonlocalIDs_(0),
530  allocatedNonlocalLength_(0),
531  nonlocalCoefs_(nullptr),
532  last_edit(0),
533  ignoreNonLocalEntries_(false)
534 {
535  this->init(n, n_local, ghost, false, type);
536 }
537 
538 
539 
540 // Default implementation for solver packages for which ghosted
541 // vectors are not yet implemented.
542 template <class T>
544  const bool fast)
545 {
546  this->init(other.size(),other.local_size(),fast,other.type());
547 }
548 
549 
550 
551 template <typename T>
552 inline
554 {
555  this->clear ();
556 }
557 
558 
559 
560 template <typename T>
561 inline
563  const numeric_index_type n_local,
564  const bool fast,
565  const ParallelType type)
566 {
567  // We default to allocating n_local local storage
568  numeric_index_type my_n_local = n_local;
569 
570  if (type == AUTOMATIC)
571  {
572  if (n == n_local)
573  this->_type = SERIAL;
574  else
575  this->_type = PARALLEL;
576  }
577  else if (type == GHOSTED)
578  {
579  // We don't yet support GHOSTED Epetra vectors, so to get the
580  // same functionality we need a SERIAL vector with local
581  // storage allocated for every entry.
582  this->_type = SERIAL;
583  my_n_local = n;
584  }
585  else
586  this->_type = type;
587 
588  libmesh_assert ((this->_type==SERIAL && n==my_n_local) ||
589  this->_type==PARALLEL);
590 
591  _map.reset(new Epetra_Map(static_cast<int>(n),
592  my_n_local,
593  0,
594  Epetra_MpiComm (this->comm().get())));
595 
596  _vec = new Epetra_Vector(*_map);
597 
598  myFirstID_ = _vec->Map().MinMyGID();
599  myNumIDs_ = _vec->Map().NumMyElements();
600 
601  // Currently we impose the restriction that NumVectors==1, so we won't
602  // need the LDA argument when calling ExtractView. Hence the "dummy" arg.
603  int dummy;
604  _vec->ExtractView(&myCoefs_, &dummy);
605 
606  this->_is_initialized = true;
607  this->_is_closed = true;
608  this->last_edit = 0;
609 
610  if (fast == false)
611  this->zero ();
612 }
613 
614 
615 template <typename T>
616 inline
618  const numeric_index_type n_local,
619  const std::vector<numeric_index_type> & /*ghost*/,
620  const bool fast,
621  const ParallelType type)
622 {
623  // TODO: we shouldn't ignore the ghost sparsity pattern
624  this->init(n, n_local, fast, type);
625 }
626 
627 
628 
629 template <typename T>
630 inline
632  const bool fast,
633  const ParallelType type)
634 {
635  this->init(n,n,fast,type);
636 }
637 
638 
639 
640 template <typename T>
641 inline
643 {
644  libmesh_assert (this->initialized());
645 
646  // Are we adding or inserting?
647  unsigned char global_last_edit = last_edit;
648  this->comm().max(global_last_edit);
649  libmesh_assert(!last_edit || last_edit == global_last_edit);
650 
651  if (global_last_edit == 1)
652  this->GlobalAssemble(Insert);
653  else if (global_last_edit == 2)
654  this->GlobalAssemble(Add);
655  else
656  libmesh_assert(!global_last_edit);
657 
658  this->_is_closed = true;
659  this->last_edit = 0;
660 }
661 
662 
663 
664 template <typename T>
665 inline
667 {
668  if (this->initialized())
669  {
670  // We might just be an interface to a user-provided _vec
671  if (this->_destroy_vec_on_exit)
672  {
673  delete _vec;
674  _vec = nullptr;
675  }
676 
677  // But we currently always own our own _map
678  _map.reset();
679  }
680 
681  this->_is_closed = this->_is_initialized = false;
682 }
683 
684 
685 
686 template <typename T>
687 inline
689 {
690  libmesh_assert (this->initialized());
691  libmesh_assert (this->closed());
692 
693  _vec->PutScalar(0.0);
694 }
695 
696 
697 
698 template <typename T>
699 inline
700 std::unique_ptr<NumericVector<T>> EpetraVector<T>::zero_clone () const
701 {
702  NumericVector<T> * cloned_vector = new EpetraVector<T>(this->comm(), AUTOMATIC);
703  cloned_vector->init(*this);
704  return std::unique_ptr<NumericVector<T>>(cloned_vector);
705 }
706 
707 
708 
709 template <typename T>
710 inline
711 std::unique_ptr<NumericVector<T>> EpetraVector<T>::clone () const
712 {
713  NumericVector<T> * cloned_vector = new EpetraVector<T>(this->comm(), AUTOMATIC);
714  cloned_vector->init(*this, true);
715  *cloned_vector = *this;
716  return std::unique_ptr<NumericVector<T>>(cloned_vector);
717 }
718 
719 
720 
721 template <typename T>
722 inline
724 {
725  libmesh_assert (this->initialized());
726 
727  return _vec->GlobalLength();
728 }
729 
730 
731 
732 template <typename T>
733 inline
735 {
736  libmesh_assert (this->initialized());
737 
738  return _vec->MyLength();
739 }
740 
741 template <typename T>
742 inline
744 {
745  libmesh_assert (this->initialized());
746 
747  return _vec->Map().MinMyGID();
748 }
749 
750 
751 
752 template <typename T>
753 inline
755 {
756  libmesh_assert (this->initialized());
757 
758  return _vec->Map().MaxMyGID()+1;
759 }
760 
761 
762 template <typename T>
763 inline
765 {
766  libmesh_assert (this->initialized());
767  libmesh_assert ( ((i >= this->first_local_index()) &&
768  (i < this->last_local_index())) );
769 
770  return (*_vec)[i-this->first_local_index()];
771 }
772 
773 
774 
775 template <typename T>
776 inline
778 {
779  libmesh_assert (this->initialized());
780 
781  T value;
782 
783  _vec->MinValue(&value);
784 
785  return value;
786 }
787 
788 
789 
790 template <typename T>
791 inline
793 {
794  libmesh_assert (this->initialized());
795 
796  T value;
797 
798  _vec->MaxValue(&value);
799 
800  return value;
801 }
802 
803 
804 
805 template <typename T>
806 inline
808 {
809  NumericVector<T>::swap(other);
810 
811  EpetraVector<T> & v = cast_ref<EpetraVector<T> &>(other);
812 
813  std::swap(_vec, v._vec);
814  _map.swap(v._map);
815  std::swap(_destroy_vec_on_exit, v._destroy_vec_on_exit);
816  std::swap(myFirstID_, v.myFirstID_);
817  std::swap(myNumIDs_, v.myNumIDs_);
818  std::swap(myCoefs_, v.myCoefs_);
819  std::swap(nonlocalIDs_, v.nonlocalIDs_);
820  std::swap(nonlocalElementSize_, v.nonlocalElementSize_);
821  std::swap(numNonlocalIDs_, v.numNonlocalIDs_);
822  std::swap(allocatedNonlocalLength_, v.allocatedNonlocalLength_);
823  std::swap(nonlocalCoefs_, v.nonlocalCoefs_);
824  std::swap(last_edit, v.last_edit);
825  std::swap(ignoreNonLocalEntries_, v.ignoreNonLocalEntries_);
826 }
827 
828 
829 // Trilinos only got serious about const in version 10.4
830 inline
832 {
833  return reinterpret_cast<int *>(const_cast<numeric_index_type *>(p));
834 }
835 
836 } // namespace libMesh
837 
838 
839 #endif // #ifdef LIBMESH_TRILINOS_HAVE_EPETRA
840 #endif // LIBMESH_TRILINOS_EPETRA_VECTOR_H
virtual NumericVector< T > & operator/=(const NumericVector< T > &v) override
bool closed()
Definition: libmesh.C:265
virtual T operator()(const numeric_index_type i) const override
std::unique_ptr< Epetra_Map > _map
virtual void pointwise_mult(const NumericVector< T > &vec1, const NumericVector< T > &vec2) override
virtual void reciprocal() override
virtual numeric_index_type last_local_index() const override
virtual void add(const numeric_index_type i, const T value) override
virtual void close() override
virtual std::unique_ptr< NumericVector< T > > zero_clone() const override
virtual numeric_index_type size() const =0
virtual void localize_to_one(std::vector< T > &v_local, const processor_id_type proc_id=0) const override
uint8_t processor_id_type
Definition: id_types.h:99
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
Definition: diff_context.h:40
const Parallel::Communicator & comm() const
virtual void insert(const T *v, const std::vector< numeric_index_type > &dof_indices) override
int inputNonlocalValues(int GID, int numValues, const double *values, bool accumulate)
virtual void conjugate() override
virtual Real max() const override
const Number zero
Definition: libmesh.h:239
virtual void zero() override
virtual T sum() const override
virtual void init(const numeric_index_type n, const numeric_index_type n_local, const bool fast=false, const ParallelType ptype=AUTOMATIC)=0
int GlobalAssemble(Epetra_CombineMode mode=Add)
virtual numeric_index_type size() const override
virtual void add_vector_transpose(const NumericVector< T > &v, const SparseMatrix< T > &A) override
int ReplaceGlobalValues(int numIDs, const int *GIDs, const double *values)
int inputValues(int numIDs, const int *GIDs, const double *values, bool accumulate)
void setIgnoreNonLocalEntries(bool flag)
dof_id_type numeric_index_type
Definition: id_types.h:92
EpetraVector(const Parallel::Communicator &comm, const ParallelType type=AUTOMATIC)
virtual NumericVector< T > & operator+=(const NumericVector< T > &v) override
virtual Real l1_norm() const override
void init(triangulateio &t)
virtual void abs() override
virtual NumericVector< T > & operator-=(const NumericVector< T > &v) override
int inputNonlocalValue(int GID, double value, bool accumulate)
virtual Real min() const override
virtual Real linfty_norm() const override
virtual void init(const numeric_index_type N, const numeric_index_type n_local, const bool fast=false, const ParallelType type=AUTOMATIC) override
ParallelType type() const
virtual void localize(std::vector< T > &v_local) const override
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void swap(NumericVector< T > &v)
virtual T dot(const NumericVector< T > &v) const override
int SumIntoGlobalValues(int numIDs, const int *GIDs, const double *values)
static PetscErrorCode Mat * A
virtual Real l2_norm() const override
void swap(Iterator &lhs, Iterator &rhs)
virtual numeric_index_type local_size() const =0
virtual void clear() override
static const bool value
Definition: xdr_io.C:109
virtual void create_subvector(NumericVector< T > &subvector, const std::vector< numeric_index_type > &rows) const override
virtual numeric_index_type first_local_index() const override
void FEoperatorequals(const EpetraVector &source)
virtual void swap(NumericVector< T > &v) override
virtual std::unique_ptr< NumericVector< T > > clone() const override
bool initialized()
Definition: libmesh.C:258
int * numeric_trilinos_cast(const numeric_index_type *p)
virtual void scale(const T factor) override
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices) override
virtual numeric_index_type local_size() const override
EpetraVector & operator=(const EpetraVector &)=delete