+ inline
+ void
+ SparseMatrix::compress ()
+ {
+ // flush buffers
+ int ierr;
+ ierr = matrix->GlobalAssemble (col_map, row_map, true);
+
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ ierr = matrix->OptimizeStorage ();
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ last_action = Zero;
+
+ compressed = true;
+ }
+
+
+
+ inline
+ SparseMatrix &
+ SparseMatrix::operator = (const double d)
+ {
+ Assert (d==0, ExcScalarAssignmentOnlyForZeroValue());
+ compress ();
+
+ const int ierr = matrix->PutScalar(d);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return *this;
+ }
+
+
+
// Inline the set() and add()
// functions, since they will be
// called frequently, and the
+ // inline "simple" functions that are
+ // called frequently and do only involve
+ // a call to some Trilinos function.
+ inline
+ unsigned int
+ SparseMatrix::m () const
+ {
+ return matrix -> NumGlobalRows();
+ }
+
+
+
+ inline
+ unsigned int
+ SparseMatrix::n () const
+ {
+ return matrix -> NumGlobalCols();
+ }
+
+
+
+ inline
+ unsigned int
+ SparseMatrix::local_size () const
+ {
+ return matrix -> NumMyRows();
+ }
+
+
+
+ inline
+ std::pair<unsigned int, unsigned int>
+ SparseMatrix::local_range () const
+ {
+ unsigned int begin, end;
+ begin = matrix -> RowMap().MinMyGID();
+ end = matrix -> RowMap().MaxMyGID()+1;
+
+ return std::make_pair (begin, end);
+ }
+
+
+
+ inline
+ unsigned int
+ SparseMatrix::n_nonzero_elements () const
+ {
+ return matrix->NumGlobalNonzeros();
+ }
+
+
+
+ inline
+ TrilinosScalar
+ SparseMatrix::l1_norm () const
+ {
+ if (matrix->Filled() == false)
+ matrix->GlobalAssemble(col_map, row_map, true);
+
+ return matrix->NormOne();
+ }
+
+
+
+ inline
+ TrilinosScalar
+ SparseMatrix::linfty_norm () const
+ {
+ if (matrix->Filled() == false)
+ matrix->GlobalAssemble(col_map, row_map, true);
+
+ return matrix->NormInf();
+ }
+
+
+
+ inline
+ TrilinosScalar
+ SparseMatrix::frobenius_norm () const
+ {
+ if (matrix->Filled() == false)
+ matrix->GlobalAssemble(col_map, row_map, true);
+
+ return matrix->NormFrobenius();
+ }
+
+
+
+ inline
+ SparseMatrix &
+ SparseMatrix::operator *= (const TrilinosScalar a)
+ {
+ const int ierr = matrix->Scale (a);
+ Assert (ierr == 0, ExcTrilinosError(ierr));
+
+ return *this;
+ }
+
+
+
+ inline
+ SparseMatrix &
+ SparseMatrix::operator /= (const TrilinosScalar a)
+ {
+ Assert (a !=0, ExcDivideByZero());
+
+ const TrilinosScalar factor = 1./a;
+
+ const int ierr = matrix->Scale (factor);
+ Assert (ierr == 0, ExcTrilinosError(ierr));
+
+ return *this;
+ }
+
+
+
+ inline
+ void
+ SparseMatrix::vmult (VectorBase &dst,
+ const VectorBase &src) const
+ {
+ Assert (&src != &dst, ExcSourceEqualsDestination());
+
+ if (matrix->Filled() == false)
+ matrix->GlobalAssemble(col_map, row_map, true);
+
+ Assert (src.vector_partitioner().SameAs(matrix->DomainMap()) == true,
+ ExcMessage ("Column map of matrix does not fit with vector map!"));
+ Assert (dst.vector_partitioner().SameAs(matrix->RangeMap()) == true,
+ ExcMessage ("Row map of matrix does not fit with vector map!"));
+
+ const int ierr = matrix->Multiply (false, src.trilinos_vector(),
+ dst.trilinos_vector());
+ Assert (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+
+
+ inline
+ void
+ SparseMatrix::Tvmult (VectorBase &dst,
+ const VectorBase &src) const
+ {
+ Assert (&src != &dst, ExcSourceEqualsDestination());
+
+ if (matrix->Filled() == false)
+ matrix->GlobalAssemble(col_map, row_map, true);
+
+ Assert (src.vector_partitioner().SameAs(matrix->RangeMap()) == true,
+ ExcMessage ("Column map of matrix does not fit with vector map!"));
+ Assert (dst.vector_partitioner().SameAs(matrix->DomainMap()) == true,
+ ExcMessage ("Row map of matrix does not fit with vector map!"));
+
+ const int ierr = matrix->Multiply (true, src.trilinos_vector(),
+ dst.trilinos_vector());
+ Assert (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+
+
+ inline
+ void
+ SparseMatrix::vmult_add (VectorBase &dst,
+ const VectorBase &src) const
+ {
+ Assert (&src != &dst, ExcSourceEqualsDestination());
+
+ // Choose to reinit the vector with fast
+ // argument set, which does not overwrite
+ // the content -- this is what we need
+ // since we're going to overwrite that
+ // anyway in the vmult operation.
+ temp_vector.reinit(dst, true);
+
+ vmult (temp_vector, src);
+ dst += temp_vector;
+ }
+
+
+
+ inline
+ void
+ SparseMatrix::Tvmult_add (VectorBase &dst,
+ const VectorBase &src) const
+ {
+ Assert (&src != &dst, ExcSourceEqualsDestination());
+
+ temp_vector.reinit(dst, true);
+
+ vmult (temp_vector, src);
+ dst += temp_vector;
+ }
+
+
+
+ inline
+ TrilinosScalar
+ SparseMatrix::matrix_norm_square (const VectorBase &v) const
+ {
+ Assert (row_map.SameAs(col_map),
+ ExcDimensionMismatch(row_map.NumGlobalElements(),
+ col_map.NumGlobalElements()));
+
+ temp_vector.reinit(v);
+
+ vmult (temp_vector, v);
+ return temp_vector*v;
+ }
+
+
+
+ inline
+ TrilinosScalar
+ SparseMatrix::matrix_scalar_product (const VectorBase &u,
+ const VectorBase &v) const
+ {
+ Assert (row_map.SameAs(col_map),
+ ExcDimensionMismatch(row_map.NumGlobalElements(),
+ col_map.NumGlobalElements()));
+
+ temp_vector.reinit(v);
+
+ vmult (temp_vector, v);
+ return u*temp_vector;
+ }
+
+
+
+ inline
+ TrilinosScalar
+ SparseMatrix::residual (VectorBase &dst,
+ const VectorBase &x,
+ const VectorBase &b) const
+ {
+ vmult (dst, x);
+ dst -= b;
+ dst *= -1.;
+
+ return dst.l2_norm();
+ }
+
+
inline
const Epetra_CrsMatrix &
SparseMatrix::trilinos_matrix () const
+ inline
+ void
+ VectorBase::reinit (const VectorBase &v,
+ const bool fast)
+ {
+ Assert (&*vector != 0,
+ ExcMessage("Vector has not been constructed properly."));
+
+ if (fast == false || local_range() != v.local_range())
+ vector = std::auto_ptr<Epetra_FEVector>(new Epetra_FEVector(*v.vector));
+ }
+
+
+
+ inline
+ void
+ VectorBase::compress ()
+ {
+ // Now pass over the information about
+ // what we did last to the vector.
+ const int ierr = vector->GlobalAssemble(last_action);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ last_action = Zero;
+
+ compressed = true;
+ }
+
+
+
+ inline
+ VectorBase &
+ VectorBase::operator = (const TrilinosScalar s)
+ {
+
+ Assert (numbers::is_finite(s),
+ ExcMessage("The given value is not finite but either "
+ "infinite or Not A Number (NaN)"));
+
+ const int ierr = vector->PutScalar(s);
+
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return *this;
+ }
+
+
+
inline
void
VectorBase::set (const std::vector<unsigned int> &indices,
+ inline
+ unsigned int
+ VectorBase::size () const
+ {
+ return (unsigned int) vector->Map().MaxAllGID() + 1 - vector->Map().MinAllGID();
+ }
+
+
+
+ inline
+ unsigned int
+ VectorBase::local_size () const
+ {
+ return (unsigned int) vector->Map().NumMyElements();
+ }
+
+
+
inline
std::pair<unsigned int, unsigned int>
VectorBase::local_range () const
+ inline
+ TrilinosScalar
+ VectorBase::operator * (const VectorBase &vec) const
+ {
+ Assert (local_range() == vec.local_range(),
+ ExcDimensionMismatch(size(), vec.size()));
+
+ TrilinosScalar result;
+
+ const int ierr = vector->Dot(*(vec.vector), &result);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return result;
+ }
+
+
+
+ inline
+ VectorBase::real_type
+ VectorBase::norm_sqr () const
+ {
+ const TrilinosScalar d = l2_norm();
+ return d*d;
+ }
+
+
+
+ inline
+ TrilinosScalar
+ VectorBase::mean_value () const
+ {
+ TrilinosScalar mean;
+
+ const int ierr = vector->MeanValue (&mean);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return mean;
+ }
+
+
+
+ inline
+ VectorBase::real_type
+ VectorBase::l1_norm () const
+ {
+ TrilinosScalar d;
+
+ const int ierr = vector->Norm1 (&d);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return d;
+ }
+
+
+
+ inline
+ VectorBase::real_type
+ VectorBase::l2_norm () const
+ {
+ TrilinosScalar d;
+
+ const int ierr = vector->Norm2 (&d);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return d;
+ }
+
+
+
+ inline
+ VectorBase::real_type
+ VectorBase::lp_norm (const TrilinosScalar p) const
+ {
+ // get a representation of the
+ // vector and loop over all
+ // the elements
+ TrilinosScalar *start_ptr;
+ int leading_dimension;
+ int ierr = vector->ExtractView (&start_ptr, &leading_dimension);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ TrilinosScalar norm = 0;
+ TrilinosScalar sum=0;
+
+ const TrilinosScalar * ptr = start_ptr;
+
+ // add up elements
+ // TODO: This
+ // won't work in parallel like
+ // this. Find out a better way to
+ // this in that case.
+ while (ptr != start_ptr+size())
+ sum += std::pow(std::fabs(*ptr++), p);
+
+ norm = std::pow(sum, static_cast<TrilinosScalar>(1./p));
+
+ return norm;
+ }
+
+
+
+ inline
+ VectorBase::real_type
+ VectorBase::linfty_norm () const
+ {
+ TrilinosScalar d;
+
+ const int ierr = vector->NormInf (&d);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return d;
+ }
+
+
+
+ // inline also scalar products, vector
+ // additions etc. since they are all
+ // representable by a single Trilinos
+ // call. This reduces the overhead of the
+ // wrapper class.
+ inline
+ VectorBase &
+ VectorBase::operator *= (const TrilinosScalar a)
+ {
+
+ Assert (numbers::is_finite(a),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+
+ const int ierr = vector->Scale(a);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return *this;
+ }
+
+
+
+ inline
+ VectorBase &
+ VectorBase::operator /= (const TrilinosScalar a)
+ {
+
+ Assert (numbers::is_finite(a),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+
+ const TrilinosScalar factor = 1./a;
+
+ Assert (numbers::is_finite(factor),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+
+ const int ierr = vector->Scale(factor);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return *this;
+ }
+
+
+
+ inline
+ VectorBase &
+ VectorBase::operator += (const VectorBase &v)
+ {
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
+
+ const int ierr = vector->Update (1.0, *(v.vector), 1.0);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return *this;
+ }
+
+
+
+ inline
+ VectorBase &
+ VectorBase::operator -= (const VectorBase &v)
+ {
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
+
+ const int ierr = vector->Update (-1.0, *(v.vector), 1.0);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return *this;
+ }
+
+
+
+ inline
+ void
+ VectorBase::add (const TrilinosScalar s)
+ {
+
+ Assert (numbers::is_finite(s),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+
+ unsigned int n_local = local_size();
+ for (unsigned int i=0; i<n_local; i++)
+ (*vector)[0][i] += s;
+ }
+
+
+
+ inline
+ void
+ VectorBase::add (const TrilinosScalar a,
+ const VectorBase &v)
+ {
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
+
+ Assert (numbers::is_finite(a),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+
+ const int ierr = vector->Update(a, *(v.vector), 1.);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+
+
+ inline
+ void
+ VectorBase::add (const TrilinosScalar a,
+ const VectorBase &v,
+ const TrilinosScalar b,
+ const VectorBase &w)
+ {
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
+ Assert (size() == w.size(),
+ ExcDimensionMismatch(size(), w.size()));
+
+ Assert (numbers::is_finite(a),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+ Assert (numbers::is_finite(b),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+
+ const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
+
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+
+
+ inline
+ void
+ VectorBase::sadd (const TrilinosScalar s,
+ const VectorBase &v)
+ {
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
+
+ Assert (numbers::is_finite(s),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+
+ const int ierr = vector->Update(1., *(v.vector), s);
+
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+
+
+ inline
+ void
+ VectorBase::sadd (const TrilinosScalar s,
+ const TrilinosScalar a,
+ const VectorBase &v)
+ {
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
+
+ Assert (numbers::is_finite(s),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+ Assert (numbers::is_finite(a),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+
+ const int ierr = vector->Update(a, *(v.vector), s);
+
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+
+
+ inline
+ void
+ VectorBase::sadd (const TrilinosScalar s,
+ const TrilinosScalar a,
+ const VectorBase &v,
+ const TrilinosScalar b,
+ const VectorBase &w)
+ {
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
+ Assert (size() == w.size(),
+ ExcDimensionMismatch(size(), w.size()));
+
+
+ Assert (numbers::is_finite(s),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+ Assert (numbers::is_finite(a),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+ Assert (numbers::is_finite(b),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+
+ const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s);
+
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+
+
+ inline
+ void
+ VectorBase::sadd (const TrilinosScalar s,
+ const TrilinosScalar a,
+ const VectorBase &v,
+ const TrilinosScalar b,
+ const VectorBase &w,
+ const TrilinosScalar c,
+ const VectorBase &x)
+ {
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
+ Assert (size() == w.size(),
+ ExcDimensionMismatch(size(), w.size()));
+ Assert (size() == x.size(),
+ ExcDimensionMismatch(size(), x.size()));
+
+ Assert (numbers::is_finite(s),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+ Assert (numbers::is_finite(a),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+ Assert (numbers::is_finite(b),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+ Assert (numbers::is_finite(c),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+
+ // Update member can only
+ // input two other vectors so
+ // do it in two steps
+ const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ const int jerr = vector->Update(c, *(x.vector), 1.);
+ Assert (jerr == 0, ExcTrilinosError(jerr));
+ }
+
+
+
+ inline
+ void
+ VectorBase::scale (const VectorBase &factors)
+ {
+ Assert (size() == factors.size(),
+ ExcDimensionMismatch(size(), factors.size()));
+
+ const int ierr = vector->Multiply (1.0, *(factors.vector), *vector, 0.0);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+
+
+ inline
+ void
+ VectorBase::equ (const TrilinosScalar a,
+ const VectorBase &v)
+ {
+
+ Assert (numbers::is_finite(a),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+
+ // If we don't have the same map, copy.
+ if (local_range() != v.local_range())
+ {
+ vector.reset();
+ last_action = Zero;
+ *vector = *v.vector;
+ *this *= a;
+ }
+ else
+ {
+ // Otherwise, just update
+ Assert (vector->Map().SameAs(v.vector->Map()) == true,
+ ExcMessage ("The Epetra maps in the assignment operator ="
+ " do not match, even though the local_range "
+ " seems to be the same. Check vector setup!"));
+ int ierr;
+ ierr = vector->GlobalAssemble(last_action);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ ierr = vector->Update(a, *v.vector, 0.0);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ last_action = Zero;
+ }
+
+ }
+
+
+
+ inline
+ void
+ VectorBase::equ (const TrilinosScalar a,
+ const VectorBase &v,
+ const TrilinosScalar b,
+ const VectorBase &w)
+ {
+
+ Assert (v.size() == w.size(),
+ ExcDimensionMismatch (v.size(), w.size()));
+
+ Assert (numbers::is_finite(a),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+ Assert (numbers::is_finite(b),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+
+ // If we don't have the same map, copy.
+ if (local_range() != v.local_range())
+ {
+ vector.reset();
+ last_action = Zero;
+ *vector = *v.vector;
+ sadd (a, b, w);
+ }
+ else
+ {
+ // Otherwise, just update
+ Assert (vector->Map().SameAs(v.vector->Map()) == true,
+ ExcMessage ("The Epetra maps in the assignment operator ="
+ " do not match, even though the local_range "
+ " seems to be the same. Check vector setup!"));
+ int ierr;
+ ierr = vector->GlobalAssemble(last_action);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ ierr = vector->Update(a, *v.vector, b, *w.vector, 0.0);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ last_action = Zero;
+ }
+ }
+
+
+
+ inline
+ void
+ VectorBase::ratio (const VectorBase &v,
+ const VectorBase &w)
+ {
+ Assert (v.size() == w.size(),
+ ExcDimensionMismatch (v.size(), w.size()));
+
+ Assert (size() == w.size(),
+ ExcDimensionMismatch (size(), w.size()));
+
+ const int ierr = vector->ReciprocalMultiply(1.0, *(w.vector),
+ *(v.vector), 0.0);
+
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+
+
inline
const Epetra_MultiVector &
VectorBase::trilinos_vector () const
- void
- SparseMatrix::compress ()
- {
- // flush buffers
- int ierr;
- ierr = matrix->GlobalAssemble (col_map, row_map, true);
-
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- ierr = matrix->OptimizeStorage ();
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- last_action = Zero;
-
- compressed = true;
- }
-
-
-
- SparseMatrix &
- SparseMatrix::operator = (const double d)
- {
- Assert (d==0, ExcScalarAssignmentOnlyForZeroValue());
-
- compress ();
-
- const int ierr = matrix->PutScalar(d);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return *this;
- }
-
-
-
void
SparseMatrix::clear_row (const unsigned int row,
const TrilinosScalar new_diag_value)
- unsigned int
- SparseMatrix::m () const
- {
- int n_rows = matrix -> NumGlobalRows();
-
- return n_rows;
- }
-
-
-
- unsigned int
- SparseMatrix::n () const
- {
- unsigned int n_cols = matrix -> NumGlobalCols();
- return n_cols;
- }
-
-
-
- unsigned int
- SparseMatrix::local_size () const
- {
- unsigned int n_rows = matrix -> NumMyRows();
-
- return n_rows;
- }
-
-
-
- std::pair<unsigned int, unsigned int>
- SparseMatrix::local_range () const
- {
- unsigned int begin, end;
- begin = matrix -> RowMap().MinMyGID();
- end = matrix -> RowMap().MaxMyGID()+1;
-
- return std::make_pair (begin, end);
- }
-
-
-
- unsigned int
- SparseMatrix::n_nonzero_elements () const
- {
- unsigned int nnz = matrix->NumGlobalNonzeros();
-
- return nnz;
- }
-
-
-
unsigned int
SparseMatrix::row_length (const unsigned int row) const
{
- TrilinosScalar
- SparseMatrix::l1_norm () const
- {
- if (matrix->Filled() == false)
- matrix->GlobalAssemble(col_map, row_map, true);
-
- TrilinosScalar result = matrix->NormOne();
-
- return result;
- }
-
-
-
- TrilinosScalar
- SparseMatrix::linfty_norm () const
- {
- if (matrix->Filled() == false)
- matrix->GlobalAssemble(col_map, row_map, true);
-
- TrilinosScalar result = matrix->NormInf();
-
- return result;
- }
-
-
-
- TrilinosScalar
- SparseMatrix::frobenius_norm () const
- {
- if (matrix->Filled() == false)
- matrix->GlobalAssemble(col_map, row_map, true);
-
- TrilinosScalar result = matrix->NormFrobenius();
-
- return result;
- }
-
-
-
- SparseMatrix &
- SparseMatrix::operator *= (const TrilinosScalar a)
- {
- const int ierr = matrix->Scale (a);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return *this;
- }
-
-
-
- SparseMatrix &
- SparseMatrix::operator /= (const TrilinosScalar a)
- {
- Assert (a !=0, ExcDivideByZero());
-
- const TrilinosScalar factor = 1./a;
-
- const int ierr = matrix->Scale (factor);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return *this;
- }
-
-
-
- void
- SparseMatrix::vmult (VectorBase &dst,
- const VectorBase &src) const
- {
- Assert (&src != &dst, ExcSourceEqualsDestination());
-
- if (matrix->Filled() == false)
- matrix->GlobalAssemble(col_map, row_map, true);
-
- Assert (src.vector_partitioner().SameAs(matrix->DomainMap()) == true,
- ExcMessage ("Column map of matrix does not fit with vector map!"));
- Assert (dst.vector_partitioner().SameAs(matrix->RangeMap()) == true,
- ExcMessage ("Row map of matrix does not fit with vector map!"));
-
- const int ierr = matrix->Multiply (false, src.trilinos_vector(),
- dst.trilinos_vector());
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
-
-
- void
- SparseMatrix::Tvmult (VectorBase &dst,
- const VectorBase &src) const
- {
- Assert (&src != &dst, ExcSourceEqualsDestination());
-
- if (matrix->Filled() == false)
- matrix->GlobalAssemble(col_map, row_map, true);
-
- Assert (src.vector_partitioner().SameAs(matrix->RangeMap()) == true,
- ExcMessage ("Column map of matrix does not fit with vector map!"));
- Assert (dst.vector_partitioner().SameAs(matrix->DomainMap()) == true,
- ExcMessage ("Row map of matrix does not fit with vector map!"));
-
- const int ierr = matrix->Multiply (true, src.trilinos_vector(),
- dst.trilinos_vector());
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
-
-
- void
- SparseMatrix::vmult_add (VectorBase &dst,
- const VectorBase &src) const
- {
- Assert (&src != &dst, ExcSourceEqualsDestination());
-
- // Choose to reinit the vector with fast
- // argument set, which does not overwrite
- // the content -- this is what we need
- // since we're going to overwrite that
- // anyway in the vmult operation.
- temp_vector.reinit(dst, true);
-
- vmult (temp_vector, src);
- dst += temp_vector;
- }
-
-
-
- void
- SparseMatrix::Tvmult_add (VectorBase &dst,
- const VectorBase &src) const
- {
- Assert (&src != &dst, ExcSourceEqualsDestination());
-
- temp_vector.reinit(dst, true);
-
- vmult (temp_vector, src);
- dst += temp_vector;
- }
-
-
-
- TrilinosScalar
- SparseMatrix::matrix_norm_square (const VectorBase &v) const
- {
- Assert (row_map.SameAs(col_map),
- ExcDimensionMismatch(row_map.NumGlobalElements(),
- col_map.NumGlobalElements()));
-
- temp_vector.reinit(v);
-
- vmult (temp_vector, v);
- return temp_vector*v;
- }
-
-
-
- TrilinosScalar
- SparseMatrix::matrix_scalar_product (const VectorBase &u,
- const VectorBase &v) const
- {
- Assert (row_map.SameAs(col_map),
- ExcDimensionMismatch(row_map.NumGlobalElements(),
- col_map.NumGlobalElements()));
-
- temp_vector.reinit(v);
-
- vmult (temp_vector, v);
- return u*temp_vector;
- }
-
-
-
- TrilinosScalar
- SparseMatrix::residual (VectorBase &dst,
- const VectorBase &x,
- const VectorBase &b) const
- {
- vmult (dst, x);
- dst -= b;
- dst *= -1.;
-
- return dst.l2_norm();
- }
-
-
-
void
SparseMatrix::mmult (SparseMatrix &C,
const SparseMatrix &B,
- // TODO: Currently this only flips a
- // flag that tells Trilinos that any
- // application should be done with
- // the transpose. However, the
- // matrix structure is not
- // reset. Can we leave it like this?
void
SparseMatrix::transpose ()
{
+ // This only flips a flag that tells
+ // Trilinos that any vmult operation
+ // should be done with the
+ // transpose. However, the matrix
+ // structure is not reset.
int ierr;
if (!matrix->UseTranspose())
- void
- VectorBase::reinit (const VectorBase &v,
- const bool fast)
- {
- Assert (&*vector != 0,
- ExcMessage("Vector has not been constructed properly."));
-
- if (fast == false || local_range() != v.local_range())
- vector = std::auto_ptr<Epetra_FEVector>(new Epetra_FEVector(*v.vector));
- }
-
-
-
- void
- VectorBase::compress ()
- {
- // Now pass over the information about
- // what we did last to the vector.
- const int ierr = vector->GlobalAssemble(last_action);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- last_action = Zero;
-
- compressed = true;
- }
-
-
-
- VectorBase &
- VectorBase::operator = (const TrilinosScalar s)
- {
-
- Assert (numbers::is_finite(s),
- ExcMessage("The given value is not finite but either "
- "infinite or Not A Number (NaN)"));
-
- compress();
-
- const int ierr = vector->PutScalar(s);
-
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return *this;
- }
-
-
-
VectorBase &
VectorBase::operator = (const VectorBase &v)
{
// particularly frequent operation
std::pair<unsigned int, unsigned int>
local_range = this->local_range ();
- std::vector<unsigned int> indices (local_range.second -
- local_range.first);
- std::vector<TrilinosScalar> values (local_range.second -
- local_range.first);
+ for (unsigned int i=local_range.first; i<local_range.second; ++i)
+ (*vector)[0][i-local_range.first] = v(i);
- for (unsigned int i=0; i<local_range.second-local_range.first; ++i)
- {
- indices[i] = i + local_range.first;
- values[i] = v(i + local_range.first);
- }
-
- set (indices.size(), &indices[0], &values[0]);
-
return *this;
}
-
-
-
- bool
- VectorBase::operator == (const VectorBase &v) const
- {
- Assert (size() == v.size(),
- ExcDimensionMismatch(size(), v.size()));
- unsigned int i;
- for (i=0; i<size(); i++)
- if ((*(v.vector))[0][i]!=(*vector)[0][i]) return false;
-
- return true;
- }
-
-
-
- bool
- VectorBase::operator != (const VectorBase &v) const
- {
- Assert (size() == v.size(),
- ExcDimensionMismatch(size(), v.size()));
-
- return (!(*this==v));
- }
-
-
-
- unsigned int
- VectorBase::size () const
- {
- return (unsigned int) vector->Map().NumGlobalElements();
- }
-
-
-
- unsigned int
- VectorBase::local_size () const
- {
- return (unsigned int) vector->Map().NumMyElements();
- }
- TrilinosScalar
- VectorBase::operator * (const VectorBase &vec) const
- {
- Assert (local_range() == vec.local_range(),
- ExcDimensionMismatch(size(), vec.size()));
-
- TrilinosScalar result;
-
- const int ierr = vector->Dot(*(vec.vector), &result);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return result;
- }
-
-
-
- VectorBase::real_type
- VectorBase::norm_sqr () const
- {
- const TrilinosScalar d = l2_norm();
- return d*d;
- }
-
-
-
- TrilinosScalar
- VectorBase::mean_value () const
- {
- TrilinosScalar mean;
-
- const int ierr = vector->MeanValue (&mean);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return mean;
- }
-
-
-
- VectorBase::real_type
- VectorBase::l1_norm () const
+ inline
+ void
+ VectorBase::add (const VectorBase &v,
+ const bool allow_different_maps)
{
- TrilinosScalar d;
-
- const int ierr = vector->Norm1 (&d);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return d;
- }
-
-
+ if (allow_different_maps == false)
+ *this += v;
+ else
+ {
+ AssertThrow (size() == v.size(),
+ ExcDimensionMismatch (size(), v.size()));
- VectorBase::real_type
- VectorBase::l2_norm () const
- {
- TrilinosScalar d;
+ Epetra_Import data_exchange (vector->Map(), v.vector->Map());
- const int ierr = vector->Norm2 (&d);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ int ierr = vector->Import(*v.vector, data_exchange, Add);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- return d;
+ last_action = Insert;
+ }
}
- VectorBase::real_type
- VectorBase::lp_norm (const TrilinosScalar p) const
+ bool
+ VectorBase::operator == (const VectorBase &v) const
{
- // get a representation of the
- // vector and loop over all
- // the elements
- TrilinosScalar *start_ptr;
- int leading_dimension;
- int ierr = vector->ExtractView (&start_ptr, &leading_dimension);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- TrilinosScalar norm = 0;
- TrilinosScalar sum=0;
-
- const TrilinosScalar * ptr = start_ptr;
-
- // add up elements
- // TODO: This
- // won't work in parallel like
- // this. Find out a better way to
- // this in that case.
- while (ptr != start_ptr+size())
- sum += std::pow(std::fabs(*ptr++), p);
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
+ if (local_size() != v.local_size())
+ return false;
- norm = std::pow(sum, static_cast<TrilinosScalar>(1./p));
+ unsigned int i;
+ for (i=0; i<local_size(); i++)
+ if ((*(v.vector))[0][i]!=(*vector)[0][i]) return false;
- return norm;
+ return true;
}
- VectorBase::real_type
- VectorBase::linfty_norm () const
+ bool
+ VectorBase::operator != (const VectorBase &v) const
{
- TrilinosScalar d;
-
- const int ierr = vector->NormInf (&d);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
- return d;
+ return (!(*this==v));
}
- VectorBase &
- VectorBase::operator *= (const TrilinosScalar a)
- {
-
- Assert (numbers::is_finite(a),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
-
- const int ierr = vector->Scale(a);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return *this;
- }
-
-
-
- VectorBase &
- VectorBase::operator /= (const TrilinosScalar a)
- {
-
- Assert (numbers::is_finite(a),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
-
- const TrilinosScalar factor = 1./a;
-
- Assert (numbers::is_finite(factor),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
-
- const int ierr = vector->Scale(factor);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return *this;
- }
-
-
-
- VectorBase &
- VectorBase::operator += (const VectorBase &v)
- {
- Assert (size() == v.size(),
- ExcDimensionMismatch(size(), v.size()));
-
- const int ierr = vector->Update (1.0, *(v.vector), 1.0);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return *this;
- }
-
-
-
- VectorBase &
- VectorBase::operator -= (const VectorBase &v)
- {
- Assert (size() == v.size(),
- ExcDimensionMismatch(size(), v.size()));
-
- const int ierr = vector->Update (-1.0, *(v.vector), 1.0);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return *this;
- }
-
-
-
- void
- VectorBase::add (const TrilinosScalar s)
- {
-
- Assert (numbers::is_finite(s),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
-
- unsigned int n_local = local_size();
- int ierr;
-
- for (unsigned int i=0; i<n_local; i++)
- {
- ierr = vector->SumIntoMyValue(i,0,s);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
- }
-
-
-
- void
- VectorBase::add (const VectorBase &v,
- const bool allow_different_maps)
- {
- if (allow_different_maps == false)
- *this += v;
- else
- {
- AssertThrow (size() == v.size(),
- ExcDimensionMismatch (size(), v.size()));
-
- Epetra_Import data_exchange (vector->Map(), v.vector->Map());
-
- int ierr = vector->Import(*v.vector, data_exchange, Add);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- last_action = Insert;
- }
- }
-
-
-
- void
- VectorBase::add (const TrilinosScalar a,
- const VectorBase &v)
- {
- Assert (size() == v.size(),
- ExcDimensionMismatch(size(), v.size()));
-
- Assert (numbers::is_finite(a),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
-
- const int ierr = vector->Update(a, *(v.vector), 1.);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
-
-
- void
- VectorBase::add (const TrilinosScalar a,
- const VectorBase &v,
- const TrilinosScalar b,
- const VectorBase &w)
- {
- Assert (size() == v.size(),
- ExcDimensionMismatch(size(), v.size()));
- Assert (size() == w.size(),
- ExcDimensionMismatch(size(), w.size()));
-
- Assert (numbers::is_finite(a),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
- Assert (numbers::is_finite(b),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
-
- const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
-
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
-
-
- void
- VectorBase::sadd (const TrilinosScalar s,
- const VectorBase &v)
- {
- Assert (size() == v.size(),
- ExcDimensionMismatch(size(), v.size()));
-
- Assert (numbers::is_finite(s),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
-
- const int ierr = vector->Update(1., *(v.vector), s);
-
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
-
-
- void
- VectorBase::sadd (const TrilinosScalar s,
- const TrilinosScalar a,
- const VectorBase &v)
- {
- Assert (size() == v.size(),
- ExcDimensionMismatch(size(), v.size()));
-
- Assert (numbers::is_finite(s),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
- Assert (numbers::is_finite(a),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
-
- const int ierr = vector->Update(a, *(v.vector), s);
-
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
-
-
- void
- VectorBase::sadd (const TrilinosScalar s,
- const TrilinosScalar a,
- const VectorBase &v,
- const TrilinosScalar b,
- const VectorBase &w)
- {
- Assert (size() == v.size(),
- ExcDimensionMismatch(size(), v.size()));
- Assert (size() == w.size(),
- ExcDimensionMismatch(size(), w.size()));
-
-
- Assert (numbers::is_finite(s),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
- Assert (numbers::is_finite(a),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
- Assert (numbers::is_finite(b),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
-
- const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s);
-
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
-
-
- void
- VectorBase::sadd (const TrilinosScalar s,
- const TrilinosScalar a,
- const VectorBase &v,
- const TrilinosScalar b,
- const VectorBase &w,
- const TrilinosScalar c,
- const VectorBase &x)
- {
- Assert (size() == v.size(),
- ExcDimensionMismatch(size(), v.size()));
- Assert (size() == w.size(),
- ExcDimensionMismatch(size(), w.size()));
- Assert (size() == x.size(),
- ExcDimensionMismatch(size(), x.size()));
-
- Assert (numbers::is_finite(s),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
- Assert (numbers::is_finite(a),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
- Assert (numbers::is_finite(b),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
- Assert (numbers::is_finite(c),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
-
- // Update member can only
- // input two other vectors so
- // do it in two steps
- const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- const int jerr = vector->Update(c, *(x.vector), 1.);
- Assert (jerr == 0, ExcTrilinosError(jerr));
- }
-
-
-
- void
- VectorBase::scale (const VectorBase &factors)
- {
- Assert (size() == factors.size(),
- ExcDimensionMismatch(size(), factors.size()));
-
- const int ierr = vector->Multiply (1.0, *(factors.vector), *vector, 0.0);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
-
-
- void
- VectorBase::equ (const TrilinosScalar a,
- const VectorBase &v)
- {
-
- Assert (numbers::is_finite(a),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
-
- // If we don't have the same map, copy.
- if (local_range() != v.local_range())
- {
- vector.reset();
- last_action = Zero;
- *vector = *v.vector;
- *this *= a;
- }
- else
- {
- // Otherwise, just update
- Assert (vector->Map().SameAs(v.vector->Map()) == true,
- ExcMessage ("The Epetra maps in the assignment operator ="
- " do not match, even though the local_range "
- " seems to be the same. Check vector setup!"));
- int ierr;
- ierr = vector->GlobalAssemble(last_action);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- ierr = vector->Update(a, *v.vector, 0.0);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- last_action = Zero;
- }
-
- }
-
-
-
- void
- VectorBase::equ (const TrilinosScalar a,
- const VectorBase &v,
- const TrilinosScalar b,
- const VectorBase &w)
- {
-
- Assert (v.size() == w.size(),
- ExcDimensionMismatch (v.size(), w.size()));
-
- Assert (numbers::is_finite(a),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
- Assert (numbers::is_finite(b),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
-
- // If we don't have the same map, copy.
- if (local_range() != v.local_range())
- {
- vector.reset();
- last_action = Zero;
- *vector = *v.vector;
- sadd (a, b, w);
- }
- else
- {
- // Otherwise, just update
- Assert (vector->Map().SameAs(v.vector->Map()) == true,
- ExcMessage ("The Epetra maps in the assignment operator ="
- " do not match, even though the local_range "
- " seems to be the same. Check vector setup!"));
- int ierr;
- ierr = vector->GlobalAssemble(last_action);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- ierr = vector->Update(a, *v.vector, b, *w.vector, 0.0);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- last_action = Zero;
- }
- }
-
-
-
- void
- VectorBase::ratio (const VectorBase &v,
- const VectorBase &w)
- {
- Assert (v.size() == w.size(),
- ExcDimensionMismatch (v.size(), w.size()));
-
- Assert (size() == w.size(),
- ExcDimensionMismatch (size(), w.size()));
-
- const int ierr = vector->ReciprocalMultiply(1.0, *(w.vector),
- *(v.vector), 0.0);
-
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
-
-
// TODO: up to now only local
// data printed out! Find a
// way to neatly output