From 178e30eba5a8adb4b7e7d2d8809b8b71fe2cab2e Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 18 May 2009 07:30:16 +0000 Subject: [PATCH] Provide a correct implemention of TrilinosWrappers::VectorBase::size() even when there is some overlap between the vectors. Move some frequently called Trilinos function with easy code (reference to some Trilinos function only) to the .h file. git-svn-id: https://svn.dealii.org/trunk@18864 0785d39b-7218-0410-832d-ea1e28bc413d --- .../lac/include/lac/trilinos_sparse_matrix.h | 277 ++++++++ .../lac/include/lac/trilinos_vector_base.h | 547 ++++++++++++++++ deal.II/lac/source/trilinos_sparse_matrix.cc | 281 +-------- deal.II/lac/source/trilinos_vector_base.cc | 592 +----------------- 4 files changed, 861 insertions(+), 836 deletions(-) diff --git a/deal.II/lac/include/lac/trilinos_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_sparse_matrix.h index 9447da225d..6f96a1378e 100755 --- a/deal.II/lac/include/lac/trilinos_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_sparse_matrix.h @@ -2203,6 +2203,41 @@ namespace TrilinosWrappers + 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 @@ -2631,6 +2666,248 @@ namespace TrilinosWrappers + // 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 + 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 diff --git a/deal.II/lac/include/lac/trilinos_vector_base.h b/deal.II/lac/include/lac/trilinos_vector_base.h index d86ca5113f..923e32657b 100644 --- a/deal.II/lac/include/lac/trilinos_vector_base.h +++ b/deal.II/lac/include/lac/trilinos_vector_base.h @@ -1105,6 +1105,53 @@ namespace TrilinosWrappers + 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(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 &indices, @@ -1220,6 +1267,24 @@ namespace TrilinosWrappers + 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 VectorBase::local_range () const @@ -1232,6 +1297,488 @@ namespace TrilinosWrappers + 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(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; iUpdate(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 diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index 4ac73a823d..e3ce75b68b 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -705,40 +705,6 @@ namespace TrilinosWrappers - 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) @@ -958,57 +924,6 @@ namespace TrilinosWrappers - 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 - 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 { @@ -1033,191 +948,6 @@ namespace TrilinosWrappers - 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, @@ -1625,15 +1355,14 @@ namespace TrilinosWrappers - // 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()) diff --git a/deal.II/lac/source/trilinos_vector_base.cc b/deal.II/lac/source/trilinos_vector_base.cc index 9039cbb9b2..8fea796f72 100644 --- a/deal.II/lac/source/trilinos_vector_base.cc +++ b/deal.II/lac/source/trilinos_vector_base.cc @@ -106,52 +106,6 @@ namespace TrilinosWrappers - 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(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) { @@ -202,62 +156,11 @@ namespace TrilinosWrappers // particularly frequent operation std::pair local_range = this->local_range (); - std::vector indices (local_range.second - - local_range.first); - std::vector values (local_range.second - - local_range.first); + for (unsigned int i=local_range.first; iMap().NumGlobalElements(); - } - - - - unsigned int - VectorBase::local_size () const - { - return (unsigned int) vector->Map().NumMyElements(); - } @@ -315,110 +218,53 @@ namespace TrilinosWrappers - 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(1./p)); + unsigned int i; + for (i=0; iNormInf (&d); - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + Assert (size() == v.size(), + ExcDimensionMismatch(size(), v.size())); - return d; + return (!(*this==v)); } @@ -487,380 +333,6 @@ namespace TrilinosWrappers - 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; iSumIntoMyValue(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 -- 2.39.5