17 #ifndef __deal2__trilinos_vector_base_h 18 #define __deal2__trilinos_vector_base_h 21 #include <deal.II/base/config.h> 23 #ifdef DEAL_II_WITH_TRILINOS 25 #include <deal.II/base/utilities.h> 26 # include <deal.II/base/std_cxx1x/shared_ptr.h> 27 # include <deal.II/base/subscriptor.h> 28 # include <deal.II/lac/exceptions.h> 29 # include <deal.II/lac/vector.h> 35 # define TrilinosScalar double 36 # include "Epetra_ConfigDefs.h" 37 # ifdef DEAL_II_WITH_MPI // only if MPI is installed 39 # include "Epetra_MpiComm.h" 41 # include "Epetra_SerialComm.h" 43 # include "Epetra_FEVector.h" 48 template <
typename number>
class Vector;
75 typedef ::types::global_dof_index
size_type;
102 VectorReference (VectorBase &vector,
103 const size_type index);
131 const VectorReference &
132 operator = (
const VectorReference &r)
const;
138 const VectorReference &
139 operator = (
const VectorReference &r);
145 const VectorReference &
146 operator = (
const TrilinosScalar &s)
const;
153 const VectorReference &
154 operator += (
const TrilinosScalar &s)
const;
161 const VectorReference &
162 operator -= (
const TrilinosScalar &s)
const;
169 const VectorReference &
170 operator *= (
const TrilinosScalar &s)
const;
177 const VectorReference &
178 operator /= (
const TrilinosScalar &s)
const;
186 operator TrilinosScalar ()
const;
193 <<
"An error with error number " << arg1
194 <<
" occurred while calling a Trilinos function");
200 size_type, size_type, size_type,
201 <<
"You tried to access element " << arg1
202 <<
" of a distributed vector, but it is not stored on " 203 <<
"the current processor. Note: the elements stored " 204 <<
"on the current processor are within the range " 205 << arg2 <<
" through " << arg3
206 <<
" but Trilinos vectors need not store contiguous " 207 <<
"ranges on each processor, and not every element in " 208 <<
"this range may in fact be stored locally.");
221 const size_type index;
229 friend class ::TrilinosWrappers::VectorBase;
281 typedef TrilinosScalar real_type;
282 typedef ::types::global_dof_index size_type;
283 typedef value_type *iterator;
284 typedef const value_type *const_iterator;
285 typedef internal::VectorReference reference;
286 typedef const internal::VectorReference const_reference;
334 const bool fast =
false);
364 void compress (::VectorOperation::values operation);
384 bool is_compressed () const;
410 operator = (const TrilinosScalar s);
440 template <typename Number>
442 operator = (const ::
Vector<Number> &v);
470 size_type size () const;
490 size_type local_size () const;
522 std::pair<size_type, size_type> local_range () const;
532 bool in_local_range (const size_type index) const;
548 IndexSet locally_owned_elements () const;
556 bool has_ghost_elements() const;
564 TrilinosScalar operator * (const
VectorBase &vec) const;
570 real_type norm_sqr () const;
576 TrilinosScalar mean_value () const;
582 TrilinosScalar minimal_value () const;
588 real_type l1_norm () const;
595 real_type l2_norm () const;
604 real_type lp_norm (const TrilinosScalar p) const;
610 real_type linfty_norm () const;
622 bool all_zero () const;
634 bool is_non_negative () const;
648 operator () (const size_type index);
656 operator () (const size_type index) const;
665 operator [] (const size_type index);
675 operator [] (const size_type index) const;
687 void extract_subvector_to (const
std::vector<size_type> &indices,
688 std::vector<TrilinosScalar> &values) const;
694 template <typename ForwardIterator, typename OutputIterator>
695 void extract_subvector_to (ForwardIterator indices_begin,
696 const ForwardIterator indices_end,
697 OutputIterator values_begin) const;
710 TrilinosScalar el (const size_type index) const;
725 const_iterator begin () const;
737 const_iterator end () const;
750 void set (const
std::vector<size_type> &indices,
751 const
std::vector<TrilinosScalar> &values);
760 void set (const
std::vector<size_type> &indices,
761 const ::
Vector<TrilinosScalar> &values);
780 void set (const size_type n_elements,
781 const size_type *indices,
782 const TrilinosScalar *values);
792 void add (const
std::vector<size_type> &indices,
793 const
std::vector<TrilinosScalar> &values);
802 void add (const
std::vector<size_type> &indices,
803 const ::
Vector<TrilinosScalar> &values);
814 void add (const size_type n_elements,
815 const size_type *indices,
816 const TrilinosScalar *values);
822 VectorBase &operator *= (const TrilinosScalar factor);
828 VectorBase &operator /= (const TrilinosScalar factor);
847 void add (const TrilinosScalar s);
860 const
bool allow_different_maps = false);
867 void add (const TrilinosScalar a,
875 void add (const TrilinosScalar a,
877 const TrilinosScalar b,
885 void sadd (const TrilinosScalar s,
893 void sadd (const TrilinosScalar s,
894 const TrilinosScalar a,
900 void sadd (const TrilinosScalar s,
901 const TrilinosScalar a,
903 const TrilinosScalar b,
911 void sadd (const TrilinosScalar s,
912 const TrilinosScalar a,
914 const TrilinosScalar b,
916 const TrilinosScalar c,
928 void scale (const
VectorBase &scaling_factors);
934 void equ (const TrilinosScalar a,
941 void equ (const TrilinosScalar a,
943 const TrilinosScalar b,
978 const Epetra_MultiVector &trilinos_vector () const;
985 Epetra_FEVector &trilinos_vector ();
993 const Epetra_Map &vector_partitioner () const;
1001 void print (const
char *format = 0) const;
1016 void print (
std::ostream &out,
1017 const
unsigned int precision = 3,
1018 const
bool scientific = true,
1019 const
bool across = true) const;
1053 std::
size_t memory_consumption () const;
1060 const MPI_Comm &get_mpi_communicator () const;
1078 << "An error with error number " << arg1
1079 << " occurred while calling a Trilinos function");
1085 size_type, size_type, size_type,
1086 << "You tried to access element " << arg1
1087 << " of a distributed vector, but only entries "
1088 << arg2 << " through " << arg3
1089 << " are stored locally and can be accessed.");
1118 Epetra_CombineMode last_action;
1148 friend class MPI::Vector;
1176 VectorReference::VectorReference (
VectorBase &vector,
1177 const size_type index)
1185 const VectorReference &
1186 VectorReference::operator = (
const VectorReference &r)
const 1192 *
this =
static_cast<TrilinosScalar
> (r);
1200 const VectorReference &
1201 VectorReference::operator = (
const VectorReference &r)
1204 *
this =
static_cast<TrilinosScalar
> (r);
1211 const VectorReference &
1212 VectorReference::operator = (
const TrilinosScalar &value)
const 1214 vector.
set (1, &index, &value);
1221 const VectorReference &
1222 VectorReference::operator += (
const TrilinosScalar &value)
const 1224 vector.
add (1, &index, &value);
1231 const VectorReference &
1232 VectorReference::operator -= (
const TrilinosScalar &value)
const 1234 TrilinosScalar new_value = -value;
1235 vector.
add (1, &index, &new_value);
1242 const VectorReference &
1243 VectorReference::operator *= (
const TrilinosScalar &value)
const 1245 TrilinosScalar new_value =
static_cast<TrilinosScalar
>(*this) * value;
1246 vector.
set (1, &index, &new_value);
1253 const VectorReference &
1254 VectorReference::operator /= (
const TrilinosScalar &value)
const 1256 TrilinosScalar new_value =
static_cast<TrilinosScalar
>(*this) / value;
1257 vector.
set (1, &index, &new_value);
1277 std::pair<size_type, size_type> range = local_range();
1279 return ((index >= range.first) && (index < range.second));
1291 if (vector->Map().LinearMap())
1293 const std::pair<size_type, size_type> x = local_range();
1296 else if (vector->Map().NumMyElements() > 0)
1298 const size_type n_indices = vector->Map().NumMyElements();
1299 #ifndef DEAL_II_USE_LARGE_INDEX_TYPE 1300 unsigned int *vector_indices = (
unsigned int *)vector->Map().MyGlobalElements();
1302 size_type *vector_indices = (size_type *)vector->Map().MyGlobalElements64();
1304 is.
add_indices(vector_indices, vector_indices+n_indices);
1323 internal::VectorReference
1326 return internal::VectorReference (*
this, index);
1332 internal::VectorReference
1335 return operator() (index);
1343 return operator() (index);
1350 std::vector<TrilinosScalar> &values)
const 1352 for (size_type i = 0; i < indices.size(); ++i)
1353 values[i] =
operator()(indices[i]);
1358 template <
typename ForwardIterator,
typename OutputIterator>
1361 const ForwardIterator indices_end,
1362 OutputIterator values_begin)
const 1364 while (indices_begin != indices_end)
1366 *values_begin = operator()(*indices_begin);
1375 VectorBase::iterator
1378 return (*vector)[0];
1384 VectorBase::iterator
1387 return (*vector)[0]+local_size();
1393 VectorBase::const_iterator
1396 return (*vector)[0];
1402 VectorBase::const_iterator
1405 return (*vector)[0]+local_size();
1415 Assert (vector.get() != 0,
1416 ExcMessage(
"Vector has not been constructed properly."));
1418 if (fast ==
false ||
1420 vector.reset (
new Epetra_FEVector(*v.
vector));
1429 ::VectorOperation::values last_action_ =
1430 ::VectorOperation::unknown;
1431 if (last_action == Add)
1432 last_action_ = ::VectorOperation::add;
1433 else if (last_action == Insert)
1434 last_action_ = ::VectorOperation::insert;
1438 compress(last_action_);
1455 Epetra_CombineMode mode = last_action;
1456 if (last_action == Zero)
1458 if (given_last_action==::VectorOperation::add)
1460 else if (given_last_action==::VectorOperation::insert)
1465 # ifdef DEAL_II_WITH_MPI 1471 double double_mode = mode;
1474 dynamic_cast<const Epetra_MpiComm *>
1475 (&vector_partitioner().Comm())->GetMpiComm());
1476 Assert(result.max-result.min<1e-5,
1477 ExcMessage (
"Not all processors agree whether the last operation on " 1478 "this vector was an addition or a set operation. This will " 1479 "prevent the compress() operation from succeeding."));
1486 const int ierr = vector->GlobalAssemble(mode);
1499 compress(VectorOperation::unknown);
1514 const int ierr = vector->PutScalar(s);
1526 const std::vector<TrilinosScalar> &values)
1532 Assert (indices.size() == values.size(),
1535 set (indices.size(), &indices[0], &values[0]);
1543 const ::Vector<TrilinosScalar> &values)
1549 Assert (indices.size() == values.size(),
1552 set (indices.size(), &indices[0], values.begin());
1560 const size_type *indices,
1561 const TrilinosScalar *values)
1567 if (last_action == Add)
1568 vector->GlobalAssemble(Add);
1570 if (last_action != Insert)
1571 last_action = Insert;
1573 for (size_type i=0; i<n_elements; ++i)
1575 const size_type row = indices[i];
1577 if (local_row == -1)
1579 const int ierr = vector->ReplaceGlobalValues (1,
1586 (*vector)[0][local_row] = values[i];
1595 const std::vector<TrilinosScalar> &values)
1600 Assert (indices.size() == values.size(),
1603 add (indices.size(), &indices[0], &values[0]);
1611 const ::Vector<TrilinosScalar> &values)
1616 Assert (indices.size() == values.size(),
1619 add (indices.size(), &indices[0], values.begin());
1627 const size_type *indices,
1628 const TrilinosScalar *values)
1634 if (last_action != Add)
1636 if (last_action == Insert)
1637 vector->GlobalAssemble(Insert);
1641 for (size_type i=0; i<n_elements; ++i)
1643 const size_type row = indices[i];
1645 if (local_row == -1)
1647 const int ierr = vector->SumIntoGlobalValues (1,
1654 (*vector)[0][local_row] += values[i];
1661 VectorBase::size_type
1664 #ifndef DEAL_II_USE_LARGE_INDEX_TYPE 1665 return (size_type) (vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID());
1667 return (size_type) (vector->Map().MaxAllGID64() + 1 - vector->Map().MinAllGID64());
1674 VectorBase::size_type
1677 return (size_type) vector->Map().NumMyElements();
1683 std::pair<VectorBase::size_type, VectorBase::size_type>
1686 #ifndef DEAL_II_USE_LARGE_INDEX_TYPE 1694 Assert (end-begin == vector->Map().NumMyElements(),
1695 ExcMessage (
"This function only makes sense if the elements that this " 1696 "vector stores on the current processor form a contiguous range. " 1697 "This does not appear to be the case for the current vector."));
1699 return std::make_pair (begin, end);
1709 ExcDifferentParallelPartitioning());
1712 TrilinosScalar result;
1714 const int ierr = vector->Dot(*(vec.
vector), &result);
1723 VectorBase::real_type
1726 const TrilinosScalar d = l2_norm();
1738 TrilinosScalar mean;
1739 const int ierr = vector->MeanValue (&mean);
1751 TrilinosScalar min_value;
1752 const int ierr = vector->MinValue (&min_value);
1761 VectorBase::real_type
1767 const int ierr = vector->Norm1 (&d);
1776 VectorBase::real_type
1782 const int ierr = vector->Norm2 (&d);
1791 VectorBase::real_type
1796 TrilinosScalar norm = 0;
1797 TrilinosScalar sum=0;
1798 const size_type n_local = local_size();
1802 for (size_type i=0; i<n_local; i++)
1803 sum += std::pow(std::fabs((*vector)[0][i]), p);
1805 norm = std::pow(sum, static_cast<TrilinosScalar>(1./p));
1813 VectorBase::real_type
1822 const int ierr = vector->NormInf (&d);
1841 const int ierr = vector->Scale(a);
1855 const TrilinosScalar factor = 1./a;
1859 const int ierr = vector->Scale(factor);
1874 ExcDifferentParallelPartitioning());
1876 const int ierr = vector->Update (1.0, *(v.
vector), 1.0);
1891 ExcDifferentParallelPartitioning());
1893 const int ierr = vector->Update (-1.0, *(v.
vector), 1.0);
1910 size_type n_local = local_size();
1911 for (size_type i=0; i<n_local; i++)
1912 (*vector)[0][i] += s;
1930 const int ierr = vector->Update(a, *(v.
vector), 1.);
1940 const TrilinosScalar b,
1954 const int ierr = vector->Update(a, *(v.
vector), b, *(w.
vector), 1.);
1974 const int ierr = vector->Update(1., *(v.
vector), s);
1984 const TrilinosScalar a,
1996 const int ierr = vector->Update(a, *(v.
vector), s);
2006 const TrilinosScalar a,
2008 const TrilinosScalar b,
2023 const int ierr = vector->Update(a, *(v.
vector), b, *(w.
vector), s);
2033 const TrilinosScalar a,
2035 const TrilinosScalar b,
2037 const TrilinosScalar c,
2058 const int ierr = vector->Update(a, *(v.
vector), b, *(w.
vector), s);
2061 const int jerr = vector->Update(c, *(x.
vector), 1.);
2062 Assert (jerr == 0, ExcTrilinosError(jerr));
2078 const int ierr = vector->Multiply (1.0, *(factors.
vector), *vector, 0.0);
2095 if (vector->Map().SameAs(v.
vector->Map())==
false)
2103 int ierr = vector->Update(a, *v.
vector, 0.0);
2117 const TrilinosScalar b,
2130 if (vector->Map().SameAs(v.
vector->Map())==
false)
2143 ExcDifferentParallelPartitioning());
2144 int ierr = vector->Update(a, *v.
vector, b, *w.
vector, 0.0);
2163 const int ierr = vector->ReciprocalMultiply(1.0, *(w.
vector),
2172 const Epetra_MultiVector &
2175 return static_cast<const Epetra_MultiVector &
>(*vector);
2193 return static_cast<const Epetra_Map &
>(vector->Map());
2202 static MPI_Comm comm;
2204 #ifdef DEAL_II_WITH_MPI 2206 const Epetra_MpiComm *mpi_comm
2207 =
dynamic_cast<const Epetra_MpiComm *
>(&vector->Map().Comm());
2208 comm = mpi_comm->Comm();
2212 comm = MPI_COMM_SELF;
2227 DEAL_II_NAMESPACE_CLOSE
2229 #endif // DEAL_II_WITH_TRILINOS
size_type local_size() const
void reinit(const VectorBase &v, const bool fast=false)
const MPI_Comm & get_mpi_communicator() const
reference operator()(const size_type index)
real_type l2_norm() const
IndexSet locally_owned_elements() const
TrilinosScalar operator*(const VectorBase &vec) const
void sadd(const TrilinosScalar s, const VectorBase &V)
::ExceptionBase & ExcMessage(std::string arg1)
bool is_compressed() const
void swap(Vector &u, Vector &v)
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
void scale(const VectorBase &scaling_factors)
reference operator[](const size_type index)
real_type l1_norm() const
#define AssertThrow(cond, exc)
TrilinosScalar minimal_value() const
bool is_finite(const double x)
void set(const std::vector< size_type > &indices, const std::vector< TrilinosScalar > &values)
VectorBase & operator+=(const VectorBase &V)
void add(const std::vector< size_type > &indices, const std::vector< TrilinosScalar > &values)
bool has_ghost_elements() const
::ExceptionBase & ExcGhostsPresent()
#define DeclException1(Exception1, type1, outsequence)
VectorBase & operator/=(const TrilinosScalar factor)
void ratio(const VectorBase &a, const VectorBase &b)
#define Assert(cond, exc)
#define DeclException0(Exception0)
TrilinosScalar value_type
void add_range(const types::global_dof_index begin, const types::global_dof_index end)
BlockCompressedSparsityPattern CompressedBlockSparsityPattern DEAL_II_DEPRECATED
const Epetra_Map & vector_partitioner() const
VectorBase & operator-=(const VectorBase &V)
void compress() DEAL_II_DEPRECATED
std_cxx1x::shared_ptr< Epetra_FEVector > vector
TrilinosScalar mean_value() const
::ExceptionBase & ExcNumberNotFinite()
void extract_subvector_to(const std::vector< size_type > &indices, std::vector< TrilinosScalar > &values) const
VectorBase & operator=(const TrilinosScalar s)
::ExceptionBase & ExcNotImplemented()
real_type linfty_norm() const
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DeclException3(Exception3, type1, type2, type3, outsequence)
const Epetra_MultiVector & trilinos_vector() const
bool in_local_range(const size_type index) const
MinMaxAvg min_max_avg(const double my_value, const MPI_Comm &mpi_communicator)
real_type lp_norm(const TrilinosScalar p) const
VectorBase & operator*=(const TrilinosScalar factor)
std::pair< size_type, size_type > local_range() const
real_type norm_sqr() const
void equ(const TrilinosScalar a, const VectorBase &V)