petsc_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 
19 
20 
21 #ifndef LIBMESH_PETSC_VECTOR_H
22 #define LIBMESH_PETSC_VECTOR_H
23 
24 
25 #include "libmesh/libmesh_config.h"
26 
27 #ifdef LIBMESH_HAVE_PETSC
28 
29 // Local includes
30 #include "libmesh/numeric_vector.h"
31 #include "libmesh/petsc_macro.h"
32 
36 EXTERN_C_FOR_PETSC_BEGIN
37 # include <petscvec.h>
38 EXTERN_C_FOR_PETSC_END
39 
40 // C++ includes
41 #include <cstddef>
42 #include <cstring>
43 #include <vector>
44 
45 namespace libMesh
46 {
47 
48 
49 
50 // forward declarations
51 template <typename T> class SparseMatrix;
52 
60 template <typename T>
61 class PetscVector : public NumericVector<T>
62 {
63 public:
64 
68  explicit
70  const ParallelType type = AUTOMATIC);
71 
75  explicit
77  const numeric_index_type n,
78  const ParallelType type = AUTOMATIC);
79 
85  const numeric_index_type n,
87  const ParallelType type = AUTOMATIC);
88 
95  const numeric_index_type N,
96  const numeric_index_type n_local,
97  const std::vector<numeric_index_type>& ghost,
98  const ParallelType type = AUTOMATIC);
99 
107  PetscVector(Vec v,
108  const Parallel::Communicator &comm
109  LIBMESH_CAN_DEFAULT_TO_COMMWORLD);
110 
115  ~PetscVector ();
116 
120  void close ();
121 
125  void clear ();
126 
131  void zero ();
132 
138  virtual AutoPtr<NumericVector<T> > zero_clone () const;
139 
143  AutoPtr<NumericVector<T> > clone () const;
144 
158  void init (const numeric_index_type N,
159  const numeric_index_type n_local,
160  const bool fast=false,
161  const ParallelType type=AUTOMATIC);
162 
166  void init (const numeric_index_type N,
167  const bool fast=false,
168  const ParallelType type=AUTOMATIC);
169 
174  virtual void init (const numeric_index_type /*N*/,
175  const numeric_index_type /*n_local*/,
176  const std::vector<numeric_index_type>& /*ghost*/,
177  const bool /*fast*/ = false,
178  const ParallelType = AUTOMATIC);
179 
184  virtual void init (const NumericVector<T>& other,
185  const bool fast = false);
186 
187  // /**
188  // * Change the dimension to that of the
189  // * vector \p V. The same applies as for
190  // * the other \p init function.
191  // *
192  // * The elements of \p V are not copied, i.e.
193  // * this function is the same as calling
194  // * \p init(V.size(),fast).
195  // */
196  // void init (const NumericVector<T>& V,
197  // const bool fast=false);
198 
202  NumericVector<T> & operator= (const T s);
203 
208 
213 
217  NumericVector<T> & operator= (const std::vector<T> &v);
218 
224  Real min () const;
225 
231  Real max () const;
232 
236  T sum () const;
237 
242  Real l1_norm () const;
243 
249  Real l2_norm () const;
250 
256  Real linfty_norm () const;
257 
265  numeric_index_type size () const;
266 
272 
278 
284 
292 
296  T operator() (const numeric_index_type i) const;
297 
303  virtual void get(const std::vector<numeric_index_type>& index, std::vector<T>& values) const;
304 
310 
316 
320  virtual void reciprocal();
321 
326  virtual void conjugate();
327 
331  void set (const numeric_index_type i, const T value);
332 
336  void add (const numeric_index_type i, const T value);
337 
343  void add (const T s);
344 
350  void add (const NumericVector<T>& V);
351 
357  void add (const T a, const NumericVector<T>& v);
358 
364  void add_vector (const std::vector<T>& v,
365  const std::vector<numeric_index_type>& dof_indices);
366 
373  void add_vector (const NumericVector<T>& V,
374  const std::vector<numeric_index_type>& dof_indices);
375 
376 
381  void add_vector (const NumericVector<T> &V,
382  const SparseMatrix<T> &A);
383 
390  void add_vector (const DenseVector<T>& V,
391  const std::vector<numeric_index_type>& dof_indices);
392 
398  void add_vector_transpose (const NumericVector<T> &V,
399  const SparseMatrix<T> &A_trans);
400 
407  const SparseMatrix<T> &A_trans);
408 
413  virtual void insert (const std::vector<T>& v,
414  const std::vector<numeric_index_type>& dof_indices);
415 
422  virtual void insert (const NumericVector<T>& V,
423  const std::vector<numeric_index_type>& dof_indices);
424 
430  virtual void insert (const DenseVector<T>& V,
431  const std::vector<numeric_index_type>& dof_indices);
432 
438  virtual void insert (const DenseSubVector<T>& V,
439  const std::vector<numeric_index_type>& dof_indices);
440 
441 
446  void scale (const T factor);
447 
452 
457  virtual void abs();
458 
463  virtual T dot(const NumericVector<T>&) const;
464 
469  virtual T indefinite_dot(const NumericVector<T>&) const;
470 
475  void localize (std::vector<T>& v_local) const;
476 
481  void localize (NumericVector<T>& v_local) const;
482 
488  void localize (NumericVector<T>& v_local,
489  const std::vector<numeric_index_type>& send_list) const;
490 
495  void localize (const numeric_index_type first_local_idx,
496  const numeric_index_type last_local_idx,
497  const std::vector<numeric_index_type>& send_list);
498 
505  void localize_to_one (std::vector<T>& v_local,
506  const processor_id_type proc_id=0) const;
507 
512  virtual void pointwise_mult (const NumericVector<T>& vec1,
513  const NumericVector<T>& vec2);
514 
521  void print_matlab(const std::string name="NULL") const;
522 
527  virtual void create_subvector(NumericVector<T>& subvector,
528  const std::vector<numeric_index_type>& rows) const;
529 
533  virtual void swap (NumericVector<T> &v);
534 
540  Vec vec () { libmesh_assert (_vec); return _vec; }
541 
542 
543 
544 private:
545 
550  Vec _vec;
551 
557  mutable bool _array_is_present;
558 
559 #ifndef NDEBUG
560 
566 #endif
567 
573  mutable Vec _local_form;
574 
579  mutable PetscScalar* _values;
580 
585  void _get_array(void) const;
586 
591  void _restore_array(void) const;
592 
596  typedef std::map<numeric_index_type,numeric_index_type> GlobalToLocalMap;
597 
603 
609 };
610 
611 
612 /*----------------------- Inline functions ----------------------------------*/
613 
614 
615 
616 template <typename T>
617 inline
619  : NumericVector<T>(comm, ptype),
620  _array_is_present(false),
621  _local_form(NULL),
622  _values(NULL),
623  _global_to_local_map(),
624  _destroy_vec_on_exit(true)
625 {
626  this->_type = ptype;
627 }
628 
629 
630 
631 template <typename T>
632 inline
634  const numeric_index_type n,
635  const ParallelType ptype)
636  : NumericVector<T>(comm, ptype),
637  _array_is_present(false),
638  _local_form(NULL),
639  _values(NULL),
640  _global_to_local_map(),
641  _destroy_vec_on_exit(true)
642 {
643  this->init(n, n, false, ptype);
644 }
645 
646 
647 
648 template <typename T>
649 inline
651  const numeric_index_type n,
653  const ParallelType ptype)
654  : NumericVector<T>(comm, ptype),
655  _array_is_present(false),
656  _local_form(NULL),
657  _values(NULL),
658  _global_to_local_map(),
659  _destroy_vec_on_exit(true)
660 {
661  this->init(n, n_local, false, ptype);
662 }
663 
664 
665 
666 template <typename T>
667 inline
669  const numeric_index_type n,
671  const std::vector<numeric_index_type>& ghost,
672  const ParallelType ptype)
673  : NumericVector<T>(comm, ptype),
674  _array_is_present(false),
675  _local_form(NULL),
676  _values(NULL),
677  _global_to_local_map(),
678  _destroy_vec_on_exit(true)
679 {
680  this->init(n, n_local, ghost, false, ptype);
681 }
682 
683 
684 
685 
686 
687 template <typename T>
688 inline
691  : NumericVector<T>(comm, AUTOMATIC),
692  _array_is_present(false),
693  _local_form(NULL),
694  _values(NULL),
695  _global_to_local_map(),
696  _destroy_vec_on_exit(false)
697 {
698  this->_vec = v;
699  this->_is_closed = true;
700  this->_is_initialized = true;
701 
702  /* We need to ask PETSc about the (local to global) ghost value
703  mapping and create the inverse mapping out of it. */
705  PetscInt petsc_local_size=0;
706  ierr = VecGetLocalSize(_vec, &petsc_local_size);
707  LIBMESH_CHKERRABORT(ierr);
708 
709  // Get the vector type from PETSc.
710  // As of Petsc 3.0.0, the VecType #define lost its const-ness, so we
711  // need to have it in the code
712 #if PETSC_VERSION_LESS_THAN(3,0,0) || !PETSC_VERSION_LESS_THAN(3,4,0)
713  // Pre-3.0 and petsc-dev (as of October 2012) use non-const versions
714  VecType ptype;
715 #else
716  const VecType ptype;
717 #endif
718  ierr = VecGetType(_vec, &ptype);
719  LIBMESH_CHKERRABORT(ierr);
720 
721  if((std::strcmp(ptype,VECSHARED) == 0) || (std::strcmp(ptype,VECMPI) == 0))
722  {
723 #if PETSC_RELEASE_LESS_THAN(3,1,1)
724  ISLocalToGlobalMapping mapping = _vec->mapping;
725 #else
726  ISLocalToGlobalMapping mapping;
727  ierr = VecGetLocalToGlobalMapping(_vec, &mapping);
728  LIBMESH_CHKERRABORT(ierr);
729 #endif
730 
731  // If is a sparsely stored vector, set up our new mapping
732  if (mapping)
733  {
734  const numeric_index_type my_local_size = static_cast<numeric_index_type>(petsc_local_size);
735  const numeric_index_type ghost_begin = static_cast<numeric_index_type>(petsc_local_size);
736 #if PETSC_RELEASE_LESS_THAN(3,4,0)
737  const numeric_index_type ghost_end = static_cast<numeric_index_type>(mapping->n);
738 #else
739  PetscInt n;
740  ierr = ISLocalToGlobalMappingGetSize(mapping, &n);
741  LIBMESH_CHKERRABORT(ierr);
742  const numeric_index_type ghost_end = static_cast<numeric_index_type>(n);
743 #endif
744 #if PETSC_RELEASE_LESS_THAN(3,1,1)
745  const PetscInt *indices = mapping->indices;
746 #else
747  const PetscInt *indices;
748  ierr = ISLocalToGlobalMappingGetIndices(mapping,&indices);
749  LIBMESH_CHKERRABORT(ierr);
750 #endif
751  for(numeric_index_type i=ghost_begin; i<ghost_end; i++)
752  _global_to_local_map[indices[i]] = i-my_local_size;
753  this->_type = GHOSTED;
754 #if !PETSC_RELEASE_LESS_THAN(3,1,1)
755  ierr = ISLocalToGlobalMappingRestoreIndices(mapping, &indices);
756  LIBMESH_CHKERRABORT(ierr);
757 #endif
758  }
759  else
760  this->_type = PARALLEL;
761  }
762  else
763  this->_type = SERIAL;
764 
765  this->close();
766 }
767 
768 
769 
770 
771 template <typename T>
772 inline
774 {
775  this->clear ();
776 }
777 
778 
779 
780 template <typename T>
781 inline
784  const bool fast,
785  const ParallelType ptype)
786 {
787  PetscErrorCode ierr=0;
788  PetscInt petsc_n=static_cast<PetscInt>(n);
789  PetscInt petsc_n_local=static_cast<PetscInt>(n_local);
790 
791 
792  // Clear initialized vectors
793  if (this->initialized())
794  this->clear();
795 
796  if (ptype == AUTOMATIC)
797  {
798  if (n == n_local)
799  this->_type = SERIAL;
800  else
801  this->_type = PARALLEL;
802  }
803  else
804  this->_type = ptype;
805 
806  libmesh_assert ((this->_type==SERIAL && n==n_local) ||
807  this->_type==PARALLEL);
808 
809  // create a sequential vector if on only 1 processor
810  if (this->_type == SERIAL)
811  {
812  ierr = VecCreateSeq (PETSC_COMM_SELF, petsc_n, &_vec);
813  CHKERRABORT(PETSC_COMM_SELF,ierr);
814 
815  ierr = VecSetFromOptions (_vec);
816  CHKERRABORT(PETSC_COMM_SELF,ierr);
817  }
818  // otherwise create an MPI-enabled vector
819  else if (this->_type == PARALLEL)
820  {
821 #ifdef LIBMESH_HAVE_MPI
822  libmesh_assert_less_equal (n_local, n);
823  ierr = VecCreateMPI (this->comm().get(), petsc_n_local, petsc_n,
824  &_vec);
825  LIBMESH_CHKERRABORT(ierr);
826 #else
827  libmesh_assert_equal_to (n_local, n);
828  ierr = VecCreateSeq (PETSC_COMM_SELF, petsc_n, &_vec);
829  CHKERRABORT(PETSC_COMM_SELF,ierr);
830 #endif
831 
832  ierr = VecSetFromOptions (_vec);
833  LIBMESH_CHKERRABORT(ierr);
834  }
835  else
836  libmesh_error();
837 
838  this->_is_initialized = true;
839  this->_is_closed = true;
840 
841 
842  if (fast == false)
843  this->zero ();
844 }
845 
846 
847 
848 template <typename T>
849 inline
851  const bool fast,
852  const ParallelType ptype)
853 {
854  this->init(n,n,fast,ptype);
855 }
856 
857 
858 
859 template <typename T>
860 inline
863  const std::vector<numeric_index_type>& ghost,
864  const bool fast,
865  const ParallelType libmesh_dbg_var(ptype))
866 {
867  PetscErrorCode ierr=0;
868  PetscInt petsc_n=static_cast<PetscInt>(n);
869  PetscInt petsc_n_local=static_cast<PetscInt>(n_local);
870  PetscInt petsc_n_ghost=static_cast<PetscInt>(ghost.size());
871 
872  // If the mesh is not disjoint, every processor will either have
873  // all the dofs, none of the dofs, or some non-zero dofs at the
874  // boundary between processors.
875  //
876  // However we can't assert this, because someone might want to
877  // construct a GHOSTED vector which doesn't include neighbor element
878  // dofs. Boyce tried to do so in user code, and we're going to want
879  // to do so in System::project_vector().
880  //
881  // libmesh_assert(n_local == 0 || n_local == n || !ghost.empty());
882 
883  libmesh_assert(sizeof(PetscInt) == sizeof(numeric_index_type));
884 
885  PetscInt* petsc_ghost = ghost.empty() ? PETSC_NULL :
886  const_cast<PetscInt*>(reinterpret_cast<const PetscInt*>(&ghost[0]));
887 
888  // Clear initialized vectors
889  if (this->initialized())
890  this->clear();
891 
892  libmesh_assert(ptype == AUTOMATIC || ptype == GHOSTED);
893  this->_type = GHOSTED;
894 
895  /* Make the global-to-local ghost cell map. */
896  for(numeric_index_type i=0; i<ghost.size(); i++)
897  {
898  _global_to_local_map[ghost[i]] = i;
899  }
900 
901  /* Create vector. */
902  ierr = VecCreateGhost (this->comm().get(), petsc_n_local, petsc_n,
903  petsc_n_ghost, petsc_ghost,
904  &_vec);
905  LIBMESH_CHKERRABORT(ierr);
906 
907  ierr = VecSetFromOptions (_vec);
908  LIBMESH_CHKERRABORT(ierr);
909 
910  this->_is_initialized = true;
911  this->_is_closed = true;
912 
913  if (fast == false)
914  this->zero ();
915 }
916 
917 
918 
919 template <typename T>
920 inline
922  const bool fast)
923 {
924  // Clear initialized vectors
925  if (this->initialized())
926  this->clear();
927 
928  const PetscVector<T>& v = libmesh_cast_ref<const PetscVector<T>&>(other);
929 
930  // Other vector should restore array.
931  if(v.initialized())
932  {
933  v._restore_array();
934  }
935 
936  this->_global_to_local_map = v._global_to_local_map;
937 
938  // Even if we're initializeing sizes based on an uninitialized or
939  // unclosed vector, *this* vector is being initialized now and is
940  // initially closed.
941  this->_is_closed = true; // v._is_closed;
942  this->_is_initialized = true; // v._is_initialized;
943 
944  this->_type = v._type;
945 
946  if (v.size() != 0)
947  {
948  PetscErrorCode ierr = 0;
949 
950  ierr = VecDuplicate (v._vec, &this->_vec);
951  LIBMESH_CHKERRABORT(ierr);
952  }
953 
954  if (fast == false)
955  this->zero ();
956 }
957 
958 
959 
960 template <typename T>
961 inline
963 {
964  this->_restore_array();
965 
966  PetscErrorCode ierr=0;
967 
968  ierr = VecAssemblyBegin(_vec);
969  LIBMESH_CHKERRABORT(ierr);
970  ierr = VecAssemblyEnd(_vec);
971  LIBMESH_CHKERRABORT(ierr);
972 
973  if(this->type() == GHOSTED)
974  {
975  ierr = VecGhostUpdateBegin(_vec,INSERT_VALUES,SCATTER_FORWARD);
976  LIBMESH_CHKERRABORT(ierr);
977  ierr = VecGhostUpdateEnd(_vec,INSERT_VALUES,SCATTER_FORWARD);
978  LIBMESH_CHKERRABORT(ierr);
979  }
980 
981  this->_is_closed = true;
982 }
983 
984 
985 
986 template <typename T>
987 inline
989 {
990  if (this->initialized())
991  this->_restore_array();
992 
993  if ((this->initialized()) && (this->_destroy_vec_on_exit))
994  {
995  PetscErrorCode ierr=0;
996 
997  ierr = LibMeshVecDestroy(&_vec);
998  LIBMESH_CHKERRABORT(ierr);
999  }
1000 
1001  this->_is_closed = this->_is_initialized = false;
1002 
1003  _global_to_local_map.clear();
1004 }
1005 
1006 
1007 
1008 template <typename T>
1009 inline
1011 {
1012  libmesh_assert(this->closed());
1013 
1014  this->_restore_array();
1015 
1016  PetscErrorCode ierr=0;
1017 
1018  PetscScalar z=0.;
1019 
1020  if(this->type() != GHOSTED)
1021  {
1022 #if PETSC_VERSION_LESS_THAN(2,3,0)
1023  // 2.2.x & earlier style
1024  ierr = VecSet (&z, _vec);
1025  LIBMESH_CHKERRABORT(ierr);
1026 #else
1027  // 2.3.x & newer
1028  ierr = VecSet (_vec, z);
1029  LIBMESH_CHKERRABORT(ierr);
1030 #endif
1031  }
1032  else
1033  {
1034  /* Vectors that include ghost values require a special
1035  handling. */
1036  Vec loc_vec;
1037  ierr = VecGhostGetLocalForm (_vec,&loc_vec);
1038  LIBMESH_CHKERRABORT(ierr);
1039 #if PETSC_VERSION_LESS_THAN(2,3,0)
1040  // 2.2.x & earlier style
1041  ierr = VecSet (&z, loc_vec);
1042  LIBMESH_CHKERRABORT(ierr);
1043 #else
1044  // 2.3.x & newer
1045  ierr = VecSet (loc_vec, z);
1046  LIBMESH_CHKERRABORT(ierr);
1047 #endif
1048  ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
1049  LIBMESH_CHKERRABORT(ierr);
1050  }
1051 }
1052 
1053 
1054 
1055 template <typename T>
1056 inline
1058 {
1059  AutoPtr<NumericVector<T> > cloned_vector
1060  (new PetscVector<T>(this->comm(), this->type()));
1061 
1062  cloned_vector->init(*this);
1063 
1064  return cloned_vector;
1065 }
1066 
1067 
1068 
1069 template <typename T>
1070 inline
1072 {
1073  AutoPtr<NumericVector<T> > cloned_vector
1074  (new PetscVector<T>(this->comm(), this->type()));
1075 
1076  cloned_vector->init(*this, true);
1077 
1078  *cloned_vector = *this;
1079 
1080  return cloned_vector;
1081 }
1082 
1083 
1084 
1085 template <typename T>
1086 inline
1088 {
1089  libmesh_assert (this->initialized());
1090 
1091  PetscErrorCode ierr=0;
1092  PetscInt petsc_size=0;
1093 
1094  if (!this->initialized())
1095  return 0;
1096 
1097  ierr = VecGetSize(_vec, &petsc_size);
1098  LIBMESH_CHKERRABORT(ierr);
1099 
1100  return static_cast<numeric_index_type>(petsc_size);
1101 }
1102 
1103 
1104 
1105 template <typename T>
1106 inline
1108 {
1109  libmesh_assert (this->initialized());
1110 
1111  PetscErrorCode ierr=0;
1112  PetscInt petsc_size=0;
1113 
1114  ierr = VecGetLocalSize(_vec, &petsc_size);
1115  LIBMESH_CHKERRABORT(ierr);
1116 
1117  return static_cast<numeric_index_type>(petsc_size);
1118 }
1119 
1120 
1121 
1122 template <typename T>
1123 inline
1125 {
1126  libmesh_assert (this->initialized());
1127 
1128  PetscErrorCode ierr=0;
1129  PetscInt petsc_first=0, petsc_last=0;
1130 
1131  ierr = VecGetOwnershipRange (_vec, &petsc_first, &petsc_last);
1132  LIBMESH_CHKERRABORT(ierr);
1133 
1134  return static_cast<numeric_index_type>(petsc_first);
1135 }
1136 
1137 
1138 
1139 template <typename T>
1140 inline
1142 {
1143  libmesh_assert (this->initialized());
1144 
1145  PetscErrorCode ierr=0;
1146  PetscInt petsc_first=0, petsc_last=0;
1147 
1148  ierr = VecGetOwnershipRange (_vec, &petsc_first, &petsc_last);
1149  LIBMESH_CHKERRABORT(ierr);
1150 
1151  return static_cast<numeric_index_type>(petsc_last);
1152 }
1153 
1154 
1155 
1156 template <typename T>
1157 inline
1159 {
1160  libmesh_assert (this->initialized());
1161 
1162  PetscErrorCode ierr=0;
1163  PetscInt petsc_first=0, petsc_last=0;
1164  ierr = VecGetOwnershipRange (_vec, &petsc_first, &petsc_last);
1165  LIBMESH_CHKERRABORT(ierr);
1166  const numeric_index_type first = static_cast<numeric_index_type>(petsc_first);
1167  const numeric_index_type last = static_cast<numeric_index_type>(petsc_last);
1168 
1169  if((i>=first) && (i<last))
1170  {
1171  return i-first;
1172  }
1173 
1174  GlobalToLocalMap::const_iterator it = _global_to_local_map.find(i);
1175 #ifndef NDEBUG
1176  const GlobalToLocalMap::const_iterator end = _global_to_local_map.end();
1177  if (it == end)
1178  {
1179  std::ostringstream error_message;
1180  error_message << "No index " << i << " in ghosted vector.\n"
1181  << "Vector contains [" << first << ',' << last << ")\n";
1182  GlobalToLocalMap::const_iterator b = _global_to_local_map.begin();
1183  if (b == end)
1184  {
1185  error_message << "And empty ghost array.\n";
1186  }
1187  else
1188  {
1189  error_message << "And ghost array {" << b->first;
1190  for (b++; b != end; b++)
1191  error_message << ',' << b->first;
1192  error_message << "}\n";
1193  }
1194 
1195  libMesh::err << error_message.str();
1196 
1197  libmesh_error();
1198  }
1199  libmesh_assert (it != _global_to_local_map.end());
1200 #endif
1201  return it->second+last-first;
1202 }
1203 
1204 
1205 
1206 template <typename T>
1207 inline
1209 {
1210  this->_get_array();
1211 
1212  const numeric_index_type local_index = this->map_global_to_local_index(i);
1213 
1214 #ifndef NDEBUG
1215  if(this->type() == GHOSTED)
1216  {
1217  libmesh_assert_less (local_index, _local_size);
1218  }
1219 #endif
1220 
1221  return static_cast<T>(_values[local_index]);
1222 }
1223 
1224 
1225 
1226 template <typename T>
1227 inline
1228 void PetscVector<T>::get(const std::vector<numeric_index_type>& index, std::vector<T>& values) const
1229 {
1230  this->_get_array();
1231 
1232  const std::size_t num = index.size();
1233  values.resize(num);
1234 
1235  for(std::size_t i=0; i<num; i++)
1236  {
1237  const numeric_index_type local_index = this->map_global_to_local_index(index[i]);
1238 #ifndef NDEBUG
1239  if(this->type() == GHOSTED)
1240  {
1241  libmesh_assert_less (local_index, _local_size);
1242  }
1243 #endif
1244  values[i] = static_cast<T>(_values[local_index]);
1245  }
1246 }
1247 
1248 
1249 
1250 template <typename T>
1251 inline
1253 {
1254  this->_restore_array();
1255 
1256  PetscErrorCode ierr=0;
1257  PetscInt index=0;
1258  PetscReal returnval=0.;
1259 
1260  ierr = VecMin (_vec, &index, &returnval);
1261  LIBMESH_CHKERRABORT(ierr);
1262 
1263  // this return value is correct: VecMin returns a PetscReal
1264  return static_cast<Real>(returnval);
1265 }
1266 
1267 
1268 
1269 template <typename T>
1270 inline
1272 {
1273  this->_restore_array();
1274 
1275  PetscErrorCode ierr=0;
1276  PetscInt index=0;
1277  PetscReal returnval=0.;
1278 
1279  ierr = VecMax (_vec, &index, &returnval);
1280  LIBMESH_CHKERRABORT(ierr);
1281 
1282  // this return value is correct: VecMax returns a PetscReal
1283  return static_cast<Real>(returnval);
1284 }
1285 
1286 
1287 
1288 template <typename T>
1289 inline
1291 {
1292  NumericVector<T>::swap(other);
1293 
1294  PetscVector<T>& v = libmesh_cast_ref<PetscVector<T>&>(other);
1295 
1296  std::swap(_vec, v._vec);
1297  std::swap(_destroy_vec_on_exit, v._destroy_vec_on_exit);
1298  std::swap(_global_to_local_map, v._global_to_local_map);
1299  std::swap(_array_is_present, v._array_is_present);
1300  std::swap(_local_form, v._local_form);
1301  std::swap(_values, v._values);
1302 }
1303 
1304 
1305 
1306 template <typename T>
1307 inline
1309 {
1310  libmesh_assert (this->initialized());
1311  if(!_array_is_present)
1312  {
1313  PetscErrorCode ierr=0;
1314  if(this->type() != GHOSTED)
1315  {
1316  ierr = VecGetArray(_vec, &_values);
1317  LIBMESH_CHKERRABORT(ierr);
1318  }
1319  else
1320  {
1321  ierr = VecGhostGetLocalForm (_vec,&_local_form);
1322  LIBMESH_CHKERRABORT(ierr);
1323  ierr = VecGetArray(_local_form, &_values);
1324  LIBMESH_CHKERRABORT(ierr);
1325 #ifndef NDEBUG
1326  PetscInt my_local_size = 0;
1327  ierr = VecGetLocalSize(_local_form, &my_local_size);
1328  LIBMESH_CHKERRABORT(ierr);
1329  _local_size = static_cast<numeric_index_type>(my_local_size);
1330 #endif
1331  }
1332  _array_is_present = true;
1333  }
1334 }
1335 
1336 
1337 
1338 template <typename T>
1339 inline
1341 {
1342  libmesh_assert (this->initialized());
1343  if(_array_is_present)
1344  {
1345  PetscErrorCode ierr=0;
1346  if(this->type() != GHOSTED)
1347  {
1348  ierr = VecRestoreArray (_vec, &_values);
1349  LIBMESH_CHKERRABORT(ierr);
1350  _values = NULL;
1351  }
1352  else
1353  {
1354  ierr = VecRestoreArray (_local_form, &_values);
1355  LIBMESH_CHKERRABORT(ierr);
1356  _values = NULL;
1357  ierr = VecGhostRestoreLocalForm (_vec,&_local_form);
1358  LIBMESH_CHKERRABORT(ierr);
1359  _local_form = NULL;
1360 #ifndef NDEBUG
1361  _local_size = 0;
1362 #endif
1363  }
1364  _array_is_present = false;
1365  }
1366 }
1367 
1368 
1369 } // namespace libMesh
1370 
1371 
1372 #endif // #ifdef LIBMESH_HAVE_PETSC
1373 #endif // LIBMESH_PETSC_VECTOR_H

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

Hosted By:
SourceForge.net Logo