trilinos_epetra_vector.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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 
25 #ifdef LIBMESH_HAVE_TRILINOS
26 
27 // Local includes
28 #include "libmesh/numeric_vector.h"
29 #include "libmesh/parallel.h"
30 
31 // 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>
37 
38 // C++ includes
39 #include <cstddef>
40 #include <vector>
41 
42 // Forward declarations
43 class Epetra_IntSerialDenseVector;
44 class Epetra_SerialDenseVector;
45 
46 namespace libMesh
47 {
48 
49 // forward declarations
50 template <typename T> class SparseMatrix;
51 
59 template <typename T>
60 class EpetraVector : public NumericVector<T>
61 {
62 public:
63 
67  explicit
69  const ParallelType type = AUTOMATIC);
70 
74  explicit
76  const numeric_index_type n,
77  const ParallelType type = AUTOMATIC);
78 
84  const numeric_index_type n,
86  const ParallelType type = AUTOMATIC);
87 
94  const numeric_index_type N,
95  const numeric_index_type n_local,
96  const std::vector<numeric_index_type>& ghost,
97  const ParallelType type = AUTOMATIC);
98 
106  EpetraVector(Epetra_Vector & v,
107  const Parallel::Communicator &comm
108  LIBMESH_CAN_DEFAULT_TO_COMMWORLD);
109 
114  ~EpetraVector ();
115 
119  void close ();
120 
124  void clear ();
125 
130  void zero ();
131 
137  virtual AutoPtr<NumericVector<T> > zero_clone () const;
138 
142  AutoPtr<NumericVector<T> > clone () const;
143 
157  void init (const numeric_index_type N,
158  const numeric_index_type n_local,
159  const bool fast=false,
160  const ParallelType type=AUTOMATIC);
161 
165  void init (const numeric_index_type N,
166  const bool fast=false,
167  const ParallelType type=AUTOMATIC);
168 
173  void init (const numeric_index_type /*N*/,
174  const numeric_index_type /*n_local*/,
175  const std::vector<numeric_index_type>& /*ghost*/,
176  const bool /*fast*/ = false,
177  const ParallelType = AUTOMATIC);
178 
183  virtual void init (const NumericVector<T>& other,
184  const bool fast = false);
185 
186  // /**
187  // * Change the dimension to that of the
188  // * vector \p V. The same applies as for
189  // * the other \p init function.
190  // *
191  // * The elements of \p V are not copied, i.e.
192  // * this function is the same as calling
193  // * \p init(V.size(),fast).
194  // */
195  // void init (const NumericVector<T>& V,
196  // const bool fast=false);
197 
201  NumericVector<T> & operator= (const T s);
202 
207 
212 
216  NumericVector<T> & operator= (const std::vector<T> &v);
217 
223  Real min () const;
224 
230  Real max () const;
231 
235  T sum () const;
236 
241  Real l1_norm () const;
242 
248  Real l2_norm () const;
249 
255  Real linfty_norm () const;
256 
264  numeric_index_type size () const;
265 
271 
277 
283 
287  T operator() (const numeric_index_type i) const;
288 
294 
300 
305 
309  virtual void reciprocal();
310 
316  virtual void conjugate();
317 
321  void set (const numeric_index_type i, const T value);
322 
326  void add (const numeric_index_type i, const T value);
327 
333  void add (const T s);
334 
340  void add (const NumericVector<T>& V);
341 
347  void add (const T a, const NumericVector<T>& v);
348 
354  void add_vector (const std::vector<T>& v,
355  const std::vector<numeric_index_type>& dof_indices);
356 
363  void add_vector (const NumericVector<T>& V,
364  const std::vector<numeric_index_type>& dof_indices);
365 
366 
371  void add_vector (const NumericVector<T> &V,
372  const SparseMatrix<T> &A);
373 
380  void add_vector (const DenseVector<T>& V,
381  const std::vector<numeric_index_type>& dof_indices);
382 
387  void add_vector_transpose (const NumericVector<T> &V,
388  const SparseMatrix<T> &A_trans);
389 
394  virtual void insert (const std::vector<T>& v,
395  const std::vector<numeric_index_type>& dof_indices);
396 
403  virtual void insert (const NumericVector<T>& V,
404  const std::vector<numeric_index_type>& dof_indices);
405 
411  virtual void insert (const DenseVector<T>& V,
412  const std::vector<numeric_index_type>& dof_indices);
413 
419  virtual void insert (const DenseSubVector<T>& V,
420  const std::vector<numeric_index_type>& dof_indices);
421 
426  void scale (const T factor);
427 
432  virtual void abs();
433 
434 
438  virtual T dot(const NumericVector<T>& V) const;
439 
444  void localize (std::vector<T>& v_local) const;
445 
450  void localize (NumericVector<T>& v_local) const;
451 
457  void localize (NumericVector<T>& v_local,
458  const std::vector<numeric_index_type>& send_list) const;
459 
464  void localize (const numeric_index_type first_local_idx,
465  const numeric_index_type last_local_idx,
466  const std::vector<numeric_index_type>& send_list);
467 
474  void localize_to_one (std::vector<T>& v_local,
475  const processor_id_type proc_id=0) const;
476 
481  virtual void pointwise_mult (const NumericVector<T>& vec1,
482  const NumericVector<T>& vec2);
483 
490  void print_matlab (const std::string name="NULL") const;
491 
496  virtual void create_subvector (NumericVector<T>& subvector,
497  const std::vector<numeric_index_type>& rows) const;
498 
502  virtual void swap (NumericVector<T> &v);
503 
509  Epetra_Vector * vec () { libmesh_assert(_vec); return _vec; }
510 
511 private:
512 
517  Epetra_Vector * _vec;
518 
522  Epetra_Map * _map;
523 
529 
530 
531 
532  /*********************************************************************
533  * The following were copied (and slightly modified) from
534  * Epetra_FEVector.h in order to allow us to use a standard
535  * Epetra_Vector... which is more compatible with other Trilinos
536  * packages such as NOX. All of this code is originally under LGPL
537  *********************************************************************/
538 
542  int SumIntoGlobalValues(int numIDs, const int* GIDs, const double* values);
543 
553  int SumIntoGlobalValues(const Epetra_IntSerialDenseVector& GIDs,
554  const Epetra_SerialDenseVector& values);
555 
559  int ReplaceGlobalValues(int numIDs, const int* GIDs, const double* values);
560 
570  int ReplaceGlobalValues(const Epetra_IntSerialDenseVector& GIDs,
571  const Epetra_SerialDenseVector& values);
572 
573  int SumIntoGlobalValues(int numIDs, const int* GIDs,
574  const int* numValuesPerID,
575  const double* values);
576 
577  int ReplaceGlobalValues(int numIDs, const int* GIDs,
578  const int* numValuesPerID,
579  const double* values);
580 
588  int GlobalAssemble(Epetra_CombineMode mode = Add);
589 
592  void setIgnoreNonLocalEntries(bool flag) {
593  ignoreNonLocalEntries_ = flag;
594  }
595 
596  void FEoperatorequals(const EpetraVector& source);
597 
598  int inputValues(int numIDs,
599  const int* GIDs, const double* values,
600  bool accumulate);
601 
602  int inputValues(int numIDs,
603  const int* GIDs, const int* numValuesPerID,
604  const double* values,
605  bool accumulate);
606 
607  int inputNonlocalValue(int GID, double value, bool accumulate);
608 
609  int inputNonlocalValues(int GID, int numValues, const double* values,
610  bool accumulate);
611 
612  void destroyNonlocalData();
613 
616  double* myCoefs_;
617 
622  double** nonlocalCoefs_;
623 
629  unsigned char last_edit;
630 
632 };
633 
634 
635 /*----------------------- Inline functions ----------------------------------*/
636 
637 
638 
639 template <typename T>
640 inline
642  const ParallelType type)
643 : NumericVector<T>(comm, type),
644  _destroy_vec_on_exit(true),
645  myFirstID_(0),
646  myNumIDs_(0),
647  myCoefs_(NULL),
648  nonlocalIDs_(NULL),
649  nonlocalElementSize_(NULL),
650  numNonlocalIDs_(0),
651  allocatedNonlocalLength_(0),
652  nonlocalCoefs_(NULL),
653  last_edit(0),
654  ignoreNonLocalEntries_(false)
655 {
656  this->_type = type;
657 }
658 
659 
660 
661 template <typename T>
662 inline
664  const numeric_index_type n,
665  const ParallelType type)
666 : NumericVector<T>(comm, type),
667  _destroy_vec_on_exit(true),
668  myFirstID_(0),
669  myNumIDs_(0),
670  myCoefs_(NULL),
671  nonlocalIDs_(NULL),
672  nonlocalElementSize_(NULL),
673  numNonlocalIDs_(0),
674  allocatedNonlocalLength_(0),
675  nonlocalCoefs_(NULL),
676  last_edit(0),
677  ignoreNonLocalEntries_(false)
678 
679 {
680  this->init(n, n, false, type);
681 }
682 
683 
684 
685 template <typename T>
686 inline
688  const numeric_index_type n,
690  const ParallelType type)
691 : NumericVector<T>(comm, type),
692  _destroy_vec_on_exit(true),
693  myFirstID_(0),
694  myNumIDs_(0),
695  myCoefs_(NULL),
696  nonlocalIDs_(NULL),
697  nonlocalElementSize_(NULL),
698  numNonlocalIDs_(0),
699  allocatedNonlocalLength_(0),
700  nonlocalCoefs_(NULL),
701  last_edit(0),
702  ignoreNonLocalEntries_(false)
703 {
704  this->init(n, n_local, false, type);
705 }
706 
707 
708 
709 
710 template <typename T>
711 inline
712 EpetraVector<T>::EpetraVector(Epetra_Vector & v,
714  : NumericVector<T>(comm, AUTOMATIC),
715  _destroy_vec_on_exit(false),
716  myFirstID_(0),
717  myNumIDs_(0),
718  myCoefs_(NULL),
719  nonlocalIDs_(NULL),
720  nonlocalElementSize_(NULL),
721  numNonlocalIDs_(0),
722  allocatedNonlocalLength_(0),
723  nonlocalCoefs_(NULL),
724  last_edit(0),
725  ignoreNonLocalEntries_(false)
726 {
727  _vec = &v;
728 
729  this->_type = PARALLEL; // FIXME - need to determine this from v!
730 
731  myFirstID_ = _vec->Map().MinMyGID();
732  myNumIDs_ = _vec->Map().NumMyElements();
733 
734  _map = new Epetra_Map(_vec->GlobalLength(),
735  _vec->MyLength(),
736  0,
737  Epetra_MpiComm (this->comm().get()));
738 
739  //Currently we impose the restriction that NumVectors==1, so we won't
740  //need the LDA argument when calling ExtractView. Hence the "dummy" arg.
741  int dummy;
742  _vec->ExtractView(&myCoefs_, &dummy);
743 
744  this->_is_closed = true;
745  this->_is_initialized = true;
746 }
747 
748 
749 
750 template <typename T>
751 inline
753  const numeric_index_type n,
755  const std::vector<numeric_index_type>& ghost,
756  const ParallelType type)
757 : NumericVector<T>(comm, AUTOMATIC),
758  _destroy_vec_on_exit(true),
759  myFirstID_(0),
760  myNumIDs_(0),
761  myCoefs_(NULL),
762  nonlocalIDs_(NULL),
763  nonlocalElementSize_(NULL),
764  numNonlocalIDs_(0),
765  allocatedNonlocalLength_(0),
766  nonlocalCoefs_(NULL),
767  last_edit(0),
768  ignoreNonLocalEntries_(false)
769 {
770  this->init(n, n_local, ghost, false, type);
771 }
772 
773 
774 
775 /* Default implementation for solver packages for which ghosted
776  vectors are not yet implemented. */
777 template <class T>
779  const bool fast)
780 {
781  this->init(other.size(),other.local_size(),fast,other.type());
782 }
783 
784 
785 
786 template <typename T>
787 inline
789 {
790  this->clear ();
791 }
792 
793 
794 
795 template <typename T>
796 inline
799  const bool fast,
800  const ParallelType type)
801 {
802  // We default to allocating n_local local storage
803  numeric_index_type my_n_local = n_local;
804 
805  if (type == AUTOMATIC)
806  {
807  if (n == n_local)
808  this->_type = SERIAL;
809  else
810  this->_type = PARALLEL;
811  }
812  else if (type == GHOSTED)
813  {
814  // We don't yet support GHOSTED Epetra vectors, so to get the
815  // same functionality we need a SERIAL vector with local
816  // storage allocated for every entry.
817  this->_type = SERIAL;
818  my_n_local = n;
819  }
820  else
821  this->_type = type;
822 
823  libmesh_assert ((this->_type==SERIAL && n==my_n_local) ||
824  this->_type==PARALLEL);
825 
826  _map = new Epetra_Map(static_cast<int>(n),
827  my_n_local,
828  0,
829  Epetra_MpiComm (this->comm().get()));
830 
831  _vec = new Epetra_Vector(*_map);
832 
833  myFirstID_ = _vec->Map().MinMyGID();
834  myNumIDs_ = _vec->Map().NumMyElements();
835 
836  //Currently we impose the restriction that NumVectors==1, so we won't
837  //need the LDA argument when calling ExtractView. Hence the "dummy" arg.
838  int dummy;
839  _vec->ExtractView(&myCoefs_, &dummy);
840 
841  this->_is_initialized = true;
842  this->_is_closed = true;
843  this->last_edit = 0;
844 
845  if (fast == false)
846  this->zero ();
847 }
848 
849 
850 template <typename T>
851 inline
854  const std::vector<numeric_index_type>& /*ghost*/,
855  const bool fast,
856  const ParallelType type)
857 {
858  // TODO: we shouldn't ignore the ghost sparsity pattern
859  this->init(n, n_local, fast, type);
860 }
861 
862 
863 
864 template <typename T>
865 inline
867  const bool fast,
868  const ParallelType type)
869 {
870  this->init(n,n,fast,type);
871 }
872 
873 
874 
875 template <typename T>
876 inline
878 {
879  libmesh_assert (this->initialized());
880 
881  // Are we adding or inserting?
882  unsigned char global_last_edit = last_edit;
883  this->comm().max(global_last_edit);
884  libmesh_assert(!last_edit || last_edit == global_last_edit);
885 
886  if (global_last_edit == 1)
887  this->GlobalAssemble(Insert);
888  else if (global_last_edit == 2)
889  this->GlobalAssemble(Add);
890  else
891  libmesh_assert(!global_last_edit);
892 
893  this->_is_closed = true;
894  this->last_edit = 0;
895 }
896 
897 
898 
899 template <typename T>
900 inline
902 {
903  if (this->initialized())
904  {
905  // We might just be an interface to a user-provided _vec
906  if (this->_destroy_vec_on_exit)
907  {
908  delete _vec;
909  _vec = NULL;
910  }
911 
912  // But we currently always own our own _map
913  delete _map;
914  _map = NULL;
915  }
916 
917  this->_is_closed = this->_is_initialized = false;
918 }
919 
920 
921 
922 template <typename T>
923 inline
925 {
926  libmesh_assert (this->initialized());
927  libmesh_assert (this->closed());
928 
929  _vec->PutScalar(0.0);
930 }
931 
932 
933 
934 template <typename T>
935 inline
937 {
938  AutoPtr<NumericVector<T> > cloned_vector
939  (new EpetraVector<T>(this->comm(), AUTOMATIC));
940 
941  cloned_vector->init(*this);
942 
943  return cloned_vector;
944 }
945 
946 
947 
948 template <typename T>
949 inline
951 {
952  AutoPtr<NumericVector<T> > cloned_vector
953  (new EpetraVector<T>(this->comm(), AUTOMATIC));
954 
955  cloned_vector->init(*this, true);
956 
957  *cloned_vector = *this;
958 
959  return cloned_vector;
960 }
961 
962 
963 
964 template <typename T>
965 inline
967 {
968  libmesh_assert (this->initialized());
969 
970  return _vec->GlobalLength();
971 }
972 
973 
974 
975 template <typename T>
976 inline
978 {
979  libmesh_assert (this->initialized());
980 
981  return _vec->MyLength();
982 }
983 
984 template <typename T>
985 inline
987 {
988  libmesh_assert (this->initialized());
989 
990  return _vec->Map().MinMyGID();
991 }
992 
993 
994 
995 template <typename T>
996 inline
998 {
999  libmesh_assert (this->initialized());
1000 
1001  return _vec->Map().MaxMyGID()+1;
1002 }
1003 
1004 
1005 template <typename T>
1006 inline
1008 {
1009  libmesh_assert (this->initialized());
1010  libmesh_assert ( ((i >= this->first_local_index()) &&
1011  (i < this->last_local_index())) );
1012 
1013  return (*_vec)[i-this->first_local_index()];
1014 }
1015 
1016 
1017 
1018 template <typename T>
1019 inline
1021 {
1022  libmesh_assert (this->initialized());
1023 
1024  T value;
1025 
1026  _vec->MinValue(&value);
1027 
1028  return value;
1029 }
1030 
1031 
1032 
1033 template <typename T>
1034 inline
1036 {
1037  libmesh_assert (this->initialized());
1038 
1039  T value;
1040 
1041  _vec->MaxValue(&value);
1042 
1043  return value;
1044 }
1045 
1046 
1047 
1048 template <typename T>
1049 inline
1051 {
1052  NumericVector<T>::swap(other);
1053 
1054  EpetraVector<T>& v = libmesh_cast_ref<EpetraVector<T>&>(other);
1055 
1056  std::swap(_vec, v._vec);
1057  std::swap(_map, v._map);
1058  std::swap(_destroy_vec_on_exit, v._destroy_vec_on_exit);
1059  std::swap(myFirstID_, v.myFirstID_);
1060  std::swap(myNumIDs_, v.myNumIDs_);
1061  std::swap(myCoefs_, v.myCoefs_);
1062  std::swap(nonlocalIDs_, v.nonlocalIDs_);
1063  std::swap(nonlocalElementSize_, v.nonlocalElementSize_);
1064  std::swap(numNonlocalIDs_, v.numNonlocalIDs_);
1065  std::swap(allocatedNonlocalLength_, v.allocatedNonlocalLength_);
1066  std::swap(nonlocalCoefs_, v.nonlocalCoefs_);
1067  std::swap(last_edit, v.last_edit);
1068  std::swap(ignoreNonLocalEntries_, v.ignoreNonLocalEntries_);
1069 }
1070 
1071 } // namespace libMesh
1072 
1073 
1074 #endif // #ifdef HAVE_EPETRA
1075 #endif // LIBMESH_TRILINOS_EPETRA_VECTOR_H

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:07 UTC

Hosted By:
SourceForge.net Logo