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"
45 DEAL_II_NAMESPACE_OPEN
48 template <
typename number>
class Vector;
55 namespace TrilinosWrappers
75 typedef ::types::global_dof_index
size_type;
102 VectorReference (VectorBase &vector,
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");
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.");
229 friend class ::TrilinosWrappers::VectorBase;
281 typedef TrilinosScalar real_type;
282 typedef ::types::global_dof_index size_type;
334 const bool fast =
false);
364 void compress (::VectorOperation::values operation);
410 operator = (const TrilinosScalar s);
440 template <typename Number>
442 operator = (const ::
Vector<Number> &v);
470 size_type
size () const;
522 std::pair<size_type, size_type>
local_range () const;
564 TrilinosScalar operator * (const
VectorBase &vec) const;
604 real_type
lp_norm (const TrilinosScalar p) 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;
688 std::
vector<TrilinosScalar> &values) const;
694 template <typename ForwardIterator, typename OutputIterator>
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,
934 void equ (const TrilinosScalar a,
941 void equ (const TrilinosScalar a,
943 const TrilinosScalar b,
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;
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.");
1139 std_cxx1x::shared_ptr<Epetra_FEVector> vector;
1146 friend class internal::VectorReference;
1148 friend class MPI::Vector;
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);
1266 VectorBase::is_compressed ()
const
1275 VectorBase::in_local_range (
const size_type index)
const
1277 std::pair<size_type, size_type> range = local_range();
1279 return ((index >= range.first) && (index < range.second));
1286 VectorBase::locally_owned_elements()
const
1291 if (vector->Map().LinearMap())
1293 const std::pair<size_type, size_type> x = local_range();
1294 is.add_range (x.first, x.second);
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();
1304 is.add_indices(vector_indices, vector_indices+n_indices);
1315 VectorBase::has_ghost_elements()
const
1323 internal::VectorReference
1324 VectorBase::operator () (
const size_type index)
1326 return internal::VectorReference (*
this, index);
1332 internal::VectorReference
1333 VectorBase::operator [] (
const size_type index)
1335 return operator() (index);
1341 VectorBase::operator [] (
const size_type index)
const
1343 return operator() (index);
1349 void VectorBase::extract_subvector_to (
const std::vector<size_type> &indices,
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>
1360 void VectorBase::extract_subvector_to (ForwardIterator indices_begin,
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
1394 VectorBase::begin()
const
1396 return (*vector)[0];
1402 VectorBase::const_iterator
1403 VectorBase::end()
const
1405 return (*vector)[0]+local_size();
1412 VectorBase::reinit (
const VectorBase &v,
1415 Assert (vector.get() != 0,
1416 ExcMessage(
"Vector has not been constructed properly."));
1418 if (fast ==
false ||
1419 vector_partitioner().SameAs(v.vector_partitioner())==
false)
1420 vector.reset (
new Epetra_FEVector(*v.vector));
1427 VectorBase::compress (
const Epetra_CombineMode last_action)
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_);
1445 VectorBase::compress (::VectorOperation::values given_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);
1497 VectorBase::compress ()
1499 compress(VectorOperation::unknown);
1506 VectorBase::operator = (
const TrilinosScalar s)
1514 const int ierr = vector->PutScalar(s);
1525 VectorBase::set (
const std::vector<size_type> &indices,
1526 const std::vector<TrilinosScalar> &values)
1532 Assert (indices.size() == values.size(),
1535 set (indices.size(), &indices[0], &values[0]);
1542 VectorBase::set (
const std::vector<size_type> &indices,
1543 const ::Vector<TrilinosScalar> &values)
1549 Assert (indices.size() == values.size(),
1552 set (indices.size(), &indices[0], values.begin());
1559 VectorBase::set (
const size_type n_elements,
1561 const TrilinosScalar *values)
1567 if (last_action == Add)
1568 vector->GlobalAssemble(Add);
1570 if (last_action != Insert)
1571 last_action = Insert;
1577 if (local_row == -1)
1579 const int ierr = vector->ReplaceGlobalValues (1,
1586 (*vector)[0][local_row] = values[i];
1594 VectorBase::add (
const std::vector<size_type> &indices,
1595 const std::vector<TrilinosScalar> &values)
1600 Assert (indices.size() == values.size(),
1603 add (indices.size(), &indices[0], &values[0]);
1610 VectorBase::add (
const std::vector<size_type> &indices,
1611 const ::Vector<TrilinosScalar> &values)
1616 Assert (indices.size() == values.size(),
1619 add (indices.size(), &indices[0], values.begin());
1626 VectorBase::add (
const size_type n_elements,
1628 const TrilinosScalar *values)
1634 if (last_action != Add)
1636 if (last_action == Insert)
1637 vector->GlobalAssemble(Insert);
1645 if (local_row == -1)
1647 const int ierr = vector->SumIntoGlobalValues (1,
1654 (*vector)[0][local_row] += values[i];
1661 VectorBase::size_type
1662 VectorBase::size ()
const
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
1675 VectorBase::local_size ()
const
1677 return (
size_type) vector->Map().NumMyElements();
1683 std::pair<VectorBase::size_type, VectorBase::size_type>
1684 VectorBase::local_range ()
const
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);
1706 VectorBase::operator * (
const VectorBase &vec)
const
1708 Assert (vector->Map().SameAs(vec.vector->Map()),
1709 ExcDifferentParallelPartitioning());
1712 TrilinosScalar result;
1714 const int ierr = vector->Dot(*(vec.vector), &result);
1723 VectorBase::real_type
1724 VectorBase::norm_sqr ()
const
1726 const TrilinosScalar d = l2_norm();
1734 VectorBase::mean_value ()
const
1738 TrilinosScalar
mean;
1739 const int ierr = vector->MeanValue (&mean);
1749 VectorBase::minimal_value ()
const
1751 TrilinosScalar min_value;
1752 const int ierr = vector->MinValue (&min_value);
1761 VectorBase::real_type
1762 VectorBase::l1_norm ()
const
1767 const int ierr = vector->Norm1 (&d);
1776 VectorBase::real_type
1777 VectorBase::l2_norm ()
const
1782 const int ierr = vector->Norm2 (&d);
1791 VectorBase::real_type
1792 VectorBase::lp_norm (
const TrilinosScalar p)
const
1796 TrilinosScalar
norm = 0;
1797 TrilinosScalar
sum=0;
1803 sum += std::pow(std::fabs((*vector)[0][i]), p);
1805 norm = std::pow(sum, static_cast<TrilinosScalar>(1./p));
1813 VectorBase::real_type
1814 VectorBase::linfty_norm ()
const
1822 const int ierr = vector->NormInf (&d);
1837 VectorBase::operator *= (
const TrilinosScalar a)
1841 const int ierr = vector->Scale(a);
1851 VectorBase::operator /= (
const TrilinosScalar a)
1855 const TrilinosScalar factor = 1./a;
1859 const int ierr = vector->Scale(factor);
1869 VectorBase::operator += (
const VectorBase &v)
1871 Assert (size() == v.size(),
1873 Assert (vector->Map().SameAs(v.vector->Map()),
1874 ExcDifferentParallelPartitioning());
1876 const int ierr = vector->Update (1.0, *(v.vector), 1.0);
1886 VectorBase::operator -= (
const VectorBase &v)
1888 Assert (size() == v.size(),
1890 Assert (vector->Map().SameAs(v.vector->Map()),
1891 ExcDifferentParallelPartitioning());
1893 const int ierr = vector->Update (-1.0, *(v.vector), 1.0);
1903 VectorBase::add (
const TrilinosScalar s)
1912 (*vector)[0][i] += s;
1919 VectorBase::add (
const TrilinosScalar a,
1920 const VectorBase &v)
1925 Assert (local_size() == v.local_size(),
1930 const int ierr = vector->Update(a, *(v.vector), 1.);
1938 VectorBase::add (
const TrilinosScalar a,
1939 const VectorBase &v,
1940 const TrilinosScalar b,
1941 const VectorBase &w)
1946 Assert (local_size() == v.local_size(),
1948 Assert (local_size() == w.local_size(),
1954 const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
1963 VectorBase::sadd (
const TrilinosScalar s,
1964 const VectorBase &v)
1969 Assert (local_size() == v.local_size(),
1974 const int ierr = vector->Update(1., *(v.vector), s);
1983 VectorBase::sadd (
const TrilinosScalar s,
1984 const TrilinosScalar a,
1985 const VectorBase &v)
1990 Assert (local_size() == v.local_size(),
1996 const int ierr = vector->Update(a, *(v.vector), s);
2005 VectorBase::sadd (
const TrilinosScalar s,
2006 const TrilinosScalar a,
2007 const VectorBase &v,
2008 const TrilinosScalar b,
2009 const VectorBase &w)
2014 Assert (local_size() == v.local_size(),
2016 Assert (local_size() == w.local_size(),
2023 const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s);
2032 VectorBase::sadd (
const TrilinosScalar s,
2033 const TrilinosScalar a,
2034 const VectorBase &v,
2035 const TrilinosScalar b,
2036 const VectorBase &w,
2037 const TrilinosScalar c,
2038 const VectorBase &x)
2043 Assert (local_size() == v.local_size(),
2045 Assert (local_size() == w.local_size(),
2047 Assert (local_size() == x.local_size(),
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));
2070 VectorBase::scale (
const VectorBase &factors)
2075 Assert (local_size() == factors.local_size(),
2078 const int ierr = vector->Multiply (1.0, *(factors.vector), *vector, 0.0);
2086 VectorBase::equ (
const TrilinosScalar a,
2087 const VectorBase &v)
2095 if (vector->Map().SameAs(v.vector->Map())==
false)
2097 *vector = *v.vector;
2103 int ierr = vector->Update(a, *v.vector, 0.0);
2115 VectorBase::equ (
const TrilinosScalar a,
2116 const VectorBase &v,
2117 const TrilinosScalar b,
2118 const VectorBase &w)
2123 Assert (v.local_size() == w.local_size(),
2130 if (vector->Map().SameAs(v.vector->Map())==
false)
2132 *vector = *v.vector;
2142 Assert (vector->Map().SameAs(w.vector->Map()),
2143 ExcDifferentParallelPartitioning());
2144 int ierr = vector->Update(a, *v.vector, b, *w.vector, 0.0);
2155 VectorBase::ratio (
const VectorBase &v,
2156 const VectorBase &w)
2158 Assert (v.local_size() == w.local_size(),
2160 Assert (local_size() == w.local_size(),
2163 const int ierr = vector->ReciprocalMultiply(1.0, *(w.vector),
2172 const Epetra_MultiVector &
2173 VectorBase::trilinos_vector ()
const
2175 return static_cast<const Epetra_MultiVector &
>(*vector);
2182 VectorBase::trilinos_vector ()
2191 VectorBase::vector_partitioner ()
const
2193 return static_cast<const Epetra_Map &
>(vector->Map());
2200 VectorBase::get_mpi_communicator ()
const
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
DeclException0(ExcGhostsPresent)
size_type local_size() const
void reinit(const VectorBase &v, const bool fast=false)
const MPI_Comm & get_mpi_communicator() const
real_type l2_norm() const
IndexSet locally_owned_elements() const
void sadd(const TrilinosScalar s, const VectorBase &V)
::ExceptionBase & ExcMessage(std::string arg1)
bool is_compressed() const
void scale(const VectorBase &scaling_factors)
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)
std::size_t memory_consumption() const
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)
void print(const char *format=0) const
T sum(const T &t, const MPI_Comm &mpi_communicator)
void ratio(const VectorBase &a, const VectorBase &b)
#define Assert(cond, exc)
TrilinosScalar value_type
BlockCompressedSparsityPattern CompressedBlockSparsityPattern DEAL_II_DEPRECATED
const Epetra_Map & vector_partitioner() const
void compress() DEAL_II_DEPRECATED
bool is_non_negative() const
double norm(const FEValuesBase< dim > &fe, const VectorSlice< const std::vector< std::vector< Tensor< 1, dim > > > > &Du)
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
Epetra_CombineMode last_action
::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
DeclException3(ExcAccessToNonlocalElement, size_type, size_type, size_type,<< "You tried to access element "<< arg1<< " of a distributed vector, but only entries "<< arg2<< " through "<< arg3<< " are stored locally and can be accessed.")
friend class internal::VectorReference
MinMaxAvg min_max_avg(const double my_value, const MPI_Comm &mpi_communicator)
real_type lp_norm(const TrilinosScalar p) const
DeclException1(ExcTrilinosError, int,<< "An error with error number "<< arg1<< " occurred while calling a Trilinos function")
std::pair< size_type, size_type > local_range() const
TrilinosScalar el(const size_type index) const
real_type norm_sqr() const
void equ(const TrilinosScalar a, const VectorBase &V)