From b9ca22fab9c26a600adf0d8f4d7731e20556efaf Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 2 Jun 2004 18:11:40 +0000 Subject: [PATCH] Merge the results from the post_5_0 development branch. git-svn-id: https://svn.dealii.org/trunk@9365 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/aclocal.m4 | 50 + deal.II/base/include/base/config.h.in | 4 + deal.II/base/include/base/logstream.h | 13 + deal.II/configure.in | 1 + .../deal.II/include/dofs/dof_renumbering.h | 2 +- deal.II/deal.II/include/dofs/dof_tools.h | 2 +- deal.II/deal.II/source/dofs/dof_accessor.cc | 75 +- .../deal.II/source/dofs/dof_constraints.cc | 14 +- deal.II/deal.II/source/fe/fe_values.cc | 37 +- deal.II/deal.II/source/numerics/data_out.cc | 43 + .../numerics/derivative_approximation.cc | 2 + .../source/numerics/error_estimator.cc | 2 + deal.II/examples/step-17/step-17.cc | 22 +- deal.II/lac/include/lac/block_matrix_base.h | 1668 +++++++++++++ deal.II/lac/include/lac/block_sparse_matrix.h | 1572 +----------- .../lac/block_sparse_matrix.templates.h | 101 +- deal.II/lac/include/lac/block_vector.h | 1678 ++----------- .../lac/include/lac/block_vector.templates.h | 574 +---- deal.II/lac/include/lac/block_vector_base.h | 2150 +++++++++++++++++ .../include/lac/petsc_block_sparse_matrix.h | 201 ++ deal.II/lac/include/lac/petsc_block_vector.h | 489 ++++ deal.II/lac/include/lac/petsc_matrix_base.h | 60 +- .../lac/petsc_parallel_block_sparse_matrix.h | 198 ++ .../include/lac/petsc_parallel_block_vector.h | 485 ++++ .../lac/petsc_parallel_sparse_matrix.h | 258 +- .../lac/include/lac/petsc_parallel_vector.h | 41 +- deal.II/lac/include/lac/petsc_solver.h | 159 +- deal.II/lac/include/lac/petsc_sparse_matrix.h | 109 + deal.II/lac/include/lac/petsc_vector.h | 33 +- deal.II/lac/include/lac/petsc_vector_base.h | 14 +- deal.II/lac/include/lac/vector.h | 43 +- deal.II/lac/include/lac/vector.templates.h | 34 +- deal.II/lac/source/block_sparse_matrix.cc | 6 +- deal.II/lac/source/block_vector.cc | 29 +- .../lac/source/petsc_block_sparse_matrix.cc | 97 + deal.II/lac/source/petsc_matrix_base.cc | 25 + .../petsc_parallel_block_sparse_matrix.cc | 109 + .../lac/source/petsc_parallel_block_vector.cc | 45 + .../source/petsc_parallel_sparse_matrix.cc | 126 +- deal.II/lac/source/petsc_solver.cc | 141 +- deal.II/lac/source/petsc_sparse_matrix.cc | 70 + tests/bits/Makefile | 1 + tests/bits/block_sparse_matrix_iterator_04.cc | 2 +- tests/bits/block_vector_iterator_01.cc | 112 + tests/bits/block_vector_iterator_02.cc | 106 + tests/bits/metis_02.cc | 122 + tests/bits/petsc_block_vector_iterator_01.cc | 116 + tests/bits/petsc_block_vector_iterator_02.cc | 110 + tests/bits/petsc_block_vector_iterator_03.cc | 365 +++ 49 files changed, 7823 insertions(+), 3893 deletions(-) create mode 100644 deal.II/lac/include/lac/block_matrix_base.h create mode 100644 deal.II/lac/include/lac/block_vector_base.h create mode 100644 deal.II/lac/include/lac/petsc_block_sparse_matrix.h create mode 100644 deal.II/lac/include/lac/petsc_block_vector.h create mode 100644 deal.II/lac/include/lac/petsc_parallel_block_sparse_matrix.h create mode 100644 deal.II/lac/include/lac/petsc_parallel_block_vector.h create mode 100644 deal.II/lac/source/petsc_block_sparse_matrix.cc create mode 100644 deal.II/lac/source/petsc_parallel_block_sparse_matrix.cc create mode 100644 deal.II/lac/source/petsc_parallel_block_vector.cc create mode 100644 tests/bits/block_vector_iterator_01.cc create mode 100644 tests/bits/block_vector_iterator_02.cc create mode 100644 tests/bits/metis_02.cc create mode 100644 tests/bits/petsc_block_vector_iterator_01.cc create mode 100644 tests/bits/petsc_block_vector_iterator_02.cc create mode 100644 tests/bits/petsc_block_vector_iterator_03.cc diff --git a/deal.II/aclocal.m4 b/deal.II/aclocal.m4 index 748d09fbfc..9b9522d0c9 100644 --- a/deal.II/aclocal.m4 +++ b/deal.II/aclocal.m4 @@ -2986,6 +2986,56 @@ AC_DEFUN(DEAL_II_CHECK_TEMPL_OP_DISAMBIGUATION_BUG, dnl +dnl ------------------------------------------------------------- +dnl Some older versions of gcc compile this, despite the 'explicit' +dnl keyword: +dnl --------------------------------- +dnl struct X { +dnl template +dnl explicit X(T); +dnl }; +dnl +dnl void f(X); +dnl +dnl int main () { f(1); } +dnl --------------------------------- +dnl +dnl Check for this misfeature. +dnl +dnl Usage: DEAL_II_CHECK_EXPLICIT_CONSTRUCTOR_BUG +dnl +dnl -------------------------------------------------------------- +AC_DEFUN(DEAL_II_CHECK_EXPLICIT_CONSTRUCTOR_BUG, dnl +[ + AC_MSG_CHECKING(for explicit template constructor bug) + AC_LANG(C++) + CXXFLAGS="$CXXFLAGSG" + AC_TRY_COMPILE( + [ + struct X { + template + explicit X(T); + }; + + void f(X); + ], + [ + f(1); + ], + [ + AC_MSG_RESULT(yes. disabling some functions) + AC_DEFINE(DEAL_II_EXPLICIT_CONSTRUCTOR_BUG, 1, + [Defined if the compiler does not honor the explicit + keyword on template constructors.]) + ], + [ + AC_MSG_RESULT(no) + ]) +]) + + + + dnl ------------------------------------------------------------- dnl The boost::shared_ptr class has a templated assignment operator dnl but no assignment operator matching the default operator diff --git a/deal.II/base/include/base/config.h.in b/deal.II/base/include/base/config.h.in index c7105dad31..b5e05bfd0b 100644 --- a/deal.II/base/include/base/config.h.in +++ b/deal.II/base/include/base/config.h.in @@ -58,6 +58,10 @@ */ #undef DEAL_II_COMPAT_MAPPING +/* Defined if the compiler does not honor the explicit keyword on template + constructors. */ +#undef DEAL_II_EXPLICIT_CONSTRUCTOR_BUG + /* Defined if the compiler needs a workaround for certain problems with taking the address of template template functions. For the details, look at aclocal.m4 in the top-level directory. */ diff --git a/deal.II/base/include/base/logstream.h b/deal.II/base/include/base/logstream.h index 56b1659776..ac60e3ced0 100644 --- a/deal.II/base/include/base/logstream.h +++ b/deal.II/base/include/base/logstream.h @@ -203,6 +203,19 @@ class LogStream * a line head the next time * something is written by this * stream. + * + * An overload of this function is needed + * anyway, since the compiler can't bind + * manipulators like @p std::endl + * directly to template arguments @p T + * like in the previous general + * template. This is due to the fact that + * @p std::endl is actually an overloaded + * set of functions for @p std::ostream, + * @p std::wostream, and potentially more + * of this kind. This function is + * therefore necessary to pick one + * element from this overload set. */ LogStream & operator<< (std::ostream& (*p) (std::ostream&)); diff --git a/deal.II/configure.in b/deal.II/configure.in index 09d9497c46..38aed7e99d 100644 --- a/deal.II/configure.in +++ b/deal.II/configure.in @@ -184,6 +184,7 @@ DEAL_II_CHECK_ANON_NAMESPACE_BUG3 DEAL_II_CHECK_SFINAE_BUG DEAL_II_CHECK_TEMPL_OP_DISAMBIGUATION_BUG DEAL_II_CHECK_ARRAY_CONDITIONAL_DECAY_BUG +DEAL_II_CHECK_EXPLICIT_CONSTRUCTOR_BUG DEAL_II_CHECK_BOOST_SHARED_PTR_ASSIGNMENT DEAL_II_HAVE_PRETTY_FUNCTION DEAL_II_HAVE_STD_ITERATOR diff --git a/deal.II/deal.II/include/dofs/dof_renumbering.h b/deal.II/deal.II/include/dofs/dof_renumbering.h index 4dfd3fdc84..0b6dddfd56 100644 --- a/deal.II/deal.II/include/dofs/dof_renumbering.h +++ b/deal.II/deal.II/include/dofs/dof_renumbering.h @@ -284,7 +284,7 @@ class DoFRenumbering * components are ordered in a * different way than suggested * by the @p FESystem object you - * use. To this end, Set up the + * use. To this end, set up the * vector @p target_component * such that the entry at index * @p i denotes the number of diff --git a/deal.II/deal.II/include/dofs/dof_tools.h b/deal.II/deal.II/include/dofs/dof_tools.h index dffee2e726..fc1f71d337 100644 --- a/deal.II/deal.II/include/dofs/dof_tools.h +++ b/deal.II/deal.II/include/dofs/dof_tools.h @@ -613,7 +613,7 @@ class DoFTools * @p dof. The size of * @p selected_dofs shall equal * dof_handler.n_dofs(). Previous - * contents of this array or + * contents of this array are * overwritten. * * If the finite element under diff --git a/deal.II/deal.II/source/dofs/dof_accessor.cc b/deal.II/deal.II/source/dofs/dof_accessor.cc index a3a2f11be7..9239f8e4eb 100644 --- a/deal.II/deal.II/source/dofs/dof_accessor.cc +++ b/deal.II/deal.II/source/dofs/dof_accessor.cc @@ -13,8 +13,9 @@ #include -#include #include +#include +#include #include #include @@ -694,6 +695,23 @@ template void DoFObjectAccessor<1,deal_II_dimension>::set_dof_values (const Vector&, PETScWrappers::Vector&) const; + +template +void +DoFObjectAccessor<1,deal_II_dimension>::get_dof_values +(const PETScWrappers::BlockVector &, Vector&) const; +template +void +DoFObjectAccessor<1,deal_II_dimension>::get_dof_values +(const PETScWrappers::BlockVector &, Vector&) const; +template +void +DoFObjectAccessor<1,deal_II_dimension>::set_dof_values +(const Vector &, PETScWrappers::BlockVector&) const; +template +void +DoFObjectAccessor<1,deal_II_dimension>::set_dof_values +(const Vector&, PETScWrappers::BlockVector&) const; #endif #if deal_II_dimension >= 2 @@ -783,6 +801,23 @@ template void DoFObjectAccessor<2,deal_II_dimension>::set_dof_values (const Vector&, PETScWrappers::Vector&) const; + +template +void +DoFObjectAccessor<2,deal_II_dimension>::get_dof_values +(const PETScWrappers::BlockVector &, Vector&) const; +template +void +DoFObjectAccessor<2,deal_II_dimension>::get_dof_values +(const PETScWrappers::BlockVector &, Vector&) const; +template +void +DoFObjectAccessor<2,deal_II_dimension>::set_dof_values +(const Vector &, PETScWrappers::BlockVector&) const; +template +void +DoFObjectAccessor<2,deal_II_dimension>::set_dof_values +(const Vector&, PETScWrappers::BlockVector&) const; #endif #endif @@ -876,6 +911,23 @@ template void DoFObjectAccessor<3,deal_II_dimension>::set_dof_values (const Vector&, PETScWrappers::Vector&) const; + +template +void +DoFObjectAccessor<3,deal_II_dimension>::get_dof_values +(const PETScWrappers::BlockVector &, Vector&) const; +template +void +DoFObjectAccessor<3,deal_II_dimension>::get_dof_values +(const PETScWrappers::BlockVector &, Vector&) const; +template +void +DoFObjectAccessor<3,deal_II_dimension>::set_dof_values +(const Vector&, PETScWrappers::BlockVector &) const; +template +void +DoFObjectAccessor<3,deal_II_dimension>::set_dof_values +(const Vector&, PETScWrappers::BlockVector&) const; #endif #endif @@ -996,6 +1048,27 @@ DoFCellAccessor:: set_dof_values_by_interpolation (const Vector&, PETScWrappers::Vector&) const; +template +void +DoFCellAccessor:: +get_interpolated_dof_values +(const PETScWrappers::BlockVector&, Vector&) const; +template +void +DoFCellAccessor:: +set_dof_values_by_interpolation +(const Vector&, PETScWrappers::BlockVector&) const; + +template +void +DoFCellAccessor:: +get_interpolated_dof_values +(const PETScWrappers::BlockVector&, Vector&) const; +template +void +DoFCellAccessor:: +set_dof_values_by_interpolation +(const Vector&, PETScWrappers::BlockVector&) const; #endif diff --git a/deal.II/deal.II/source/dofs/dof_constraints.cc b/deal.II/deal.II/source/dofs/dof_constraints.cc index ec70e24e9e..47a861a5d9 100644 --- a/deal.II/deal.II/source/dofs/dof_constraints.cc +++ b/deal.II/deal.II/source/dofs/dof_constraints.cc @@ -19,13 +19,17 @@ #include #include #include +#include +#include +#include #include +#include #include +#include #include +#include #include -#include -#include -#include +#include #include #include @@ -1272,7 +1276,9 @@ VECTOR_FUNCTIONS(BlockVector); #ifdef DEAL_II_USE_PETSC VECTOR_FUNCTIONS(PETScWrappers::Vector); +VECTOR_FUNCTIONS(PETScWrappers::BlockVector); VECTOR_FUNCTIONS(PETScWrappers::MPI::Vector); +VECTOR_FUNCTIONS(PETScWrappers::MPI::BlockVector); #endif @@ -1319,5 +1325,7 @@ MATRIX_FUNCTIONS(BlockSparseMatrix); #ifdef DEAL_II_USE_PETSC MATRIX_FUNCTIONS(PETScWrappers::SparseMatrix); +MATRIX_FUNCTIONS(PETScWrappers::BlockSparseMatrix); MATRIX_FUNCTIONS(PETScWrappers::MPI::SparseMatrix); +MATRIX_FUNCTIONS(PETScWrappers::MPI::BlockSparseMatrix); #endif diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc index 913b37d8a0..fdead32927 100644 --- a/deal.II/deal.II/source/fe/fe_values.cc +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -17,13 +17,14 @@ #include #include #include +#include +#include +#include +#include #include #include #include #include -#include -#include -#include #include @@ -911,6 +912,14 @@ template void FEValuesBase::get_function_values (const PETScWrappers::Vector &, std::vector&) const; +template +void FEValuesBase::get_function_values +(const PETScWrappers::BlockVector &, + std::vector&) const; +template +void FEValuesBase::get_function_values +(const PETScWrappers::BlockVector &, + std::vector&) const; #endif //----------------------------------------------------------------------------- @@ -941,6 +950,10 @@ template void FEValuesBase::get_function_values (const PETScWrappers::Vector &, std::vector > &) const; +template +void FEValuesBase::get_function_values +(const PETScWrappers::BlockVector &, + std::vector > &) const; #endif //----------------------------------------------------------------------------- @@ -967,6 +980,10 @@ template void FEValuesBase::get_function_grads (const PETScWrappers::Vector &, std::vector > &) const; +template +void FEValuesBase::get_function_grads +(const PETScWrappers::BlockVector &, + std::vector > &) const; #endif //----------------------------------------------------------------------------- @@ -993,6 +1010,11 @@ template void FEValuesBase::get_function_grads (const PETScWrappers::Vector &, std::vector > > &) const; + +template +void FEValuesBase::get_function_grads +(const PETScWrappers::BlockVector &, + std::vector > > &) const; #endif //----------------------------------------------------------------------------- @@ -1015,6 +1037,10 @@ template void FEValuesBase::get_function_2nd_derivatives (const PETScWrappers::Vector &, std::vector > &) const; +template +void FEValuesBase::get_function_2nd_derivatives +(const PETScWrappers::BlockVector &, + std::vector > &) const; #endif //----------------------------------------------------------------------------- @@ -1037,4 +1063,9 @@ template void FEValuesBase::get_function_2nd_derivatives (const PETScWrappers::Vector &, std::vector > > &) const; + +template +void FEValuesBase::get_function_2nd_derivatives +(const PETScWrappers::BlockVector &, + std::vector > > &) const; #endif diff --git a/deal.II/deal.II/source/numerics/data_out.cc b/deal.II/deal.II/source/numerics/data_out.cc index 1c57fb71ac..6f39903eb3 100644 --- a/deal.II/deal.II/source/numerics/data_out.cc +++ b/deal.II/deal.II/source/numerics/data_out.cc @@ -18,6 +18,7 @@ #include #include #include +#include #include #include #include @@ -721,6 +722,13 @@ add_data_vector ( const PETScWrappers::Vector &vec, const std::vector &name, const DataVectorType type); + +template void +DataOut_DoFData:: +add_data_vector ( + const PETScWrappers::BlockVector &vec, + const std::vector &name, + const DataVectorType type); #endif template void @@ -744,6 +752,13 @@ add_data_vector ( const PETScWrappers::Vector &vec, const std::string &name, const DataVectorType type); + +template void +DataOut_DoFData:: +add_data_vector ( + const PETScWrappers::BlockVector &vec, + const std::string &name, + const DataVectorType type); #endif @@ -771,6 +786,13 @@ add_data_vector ( const PETScWrappers::Vector &vec, const std::vector &name, const DataVectorType type); + +template void +DataOut_DoFData:: +add_data_vector ( + const PETScWrappers::BlockVector &vec, + const std::vector &name, + const DataVectorType type); #endif template void @@ -794,6 +816,13 @@ add_data_vector ( const PETScWrappers::Vector &vec, const std::string &name, const DataVectorType type); + +template void +DataOut_DoFData:: +add_data_vector ( + const PETScWrappers::BlockVector &vec, + const std::string &name, + const DataVectorType type); #endif @@ -826,6 +855,13 @@ add_data_vector ( const PETScWrappers::Vector &vec, const std::vector &name, const DataVectorType type); + +template void +DataOut_DoFData:: +add_data_vector ( + const PETScWrappers::BlockVector &vec, + const std::vector &name, + const DataVectorType type); #endif template void @@ -849,6 +885,13 @@ add_data_vector ( const PETScWrappers::Vector &vec, const std::string &name, const DataVectorType type); + +template void +DataOut_DoFData:: +add_data_vector ( + const PETScWrappers::BlockVector &vec, + const std::string &name, + const DataVectorType type); #endif #endif diff --git a/deal.II/deal.II/source/numerics/derivative_approximation.cc b/deal.II/deal.II/source/numerics/derivative_approximation.cc index d569fcc80d..71005d19da 100644 --- a/deal.II/deal.II/source/numerics/derivative_approximation.cc +++ b/deal.II/deal.II/source/numerics/derivative_approximation.cc @@ -18,6 +18,7 @@ #include #include #include +#include #include #include #include @@ -756,6 +757,7 @@ INSTANTIATE(BlockVector); #ifdef DEAL_II_USE_PETSC INSTANTIATE(PETScWrappers::Vector); +INSTANTIATE(PETScWrappers::BlockVector); #endif diff --git a/deal.II/deal.II/source/numerics/error_estimator.cc b/deal.II/deal.II/source/numerics/error_estimator.cc index e84bdde0ca..724ab95752 100644 --- a/deal.II/deal.II/source/numerics/error_estimator.cc +++ b/deal.II/deal.II/source/numerics/error_estimator.cc @@ -18,6 +18,7 @@ #include #include #include +#include #include #include #include @@ -1338,4 +1339,5 @@ INSTANTIATE(BlockVector); #ifdef DEAL_II_USE_PETSC INSTANTIATE(PETScWrappers::Vector); +INSTANTIATE(PETScWrappers::BlockVector); #endif diff --git a/deal.II/examples/step-17/step-17.cc b/deal.II/examples/step-17/step-17.cc index ce3a2ace8d..69e8f18ae8 100644 --- a/deal.II/examples/step-17/step-17.cc +++ b/deal.II/examples/step-17/step-17.cc @@ -446,13 +446,29 @@ void ElasticProblem::setup_system () // vectors. Since they all need to work in // parallel, we have to pass them an MPI // communication object, as well as their - // global sizes, and also how many rows out - // of this global size are to be stored - // locally: + // global sizes (both dimensions are equal + // to the number of degrees of freedom), + // and also how many rows out of this + // global size are to be stored locally + // (``n_local_dofs''). In addition, PETSc + // needs to know how to partition the + // columns in the chunk of the matrix that + // is stored locally; for square matrices, + // the columns should be partitioned in the + // same way as the rows (indicated by the + // second ``n_local_dofs'' in the call) but + // in the case of rectangular matrices one + // has to partition the columns in the same + // way as vectors are partitioned with + // which the matrix is multiplied, while + // rows have to partitioned in the same way + // as destination vectors of matrix-vector + // multiplications: system_matrix.reinit (mpi_communicator, dof_handler.n_dofs(), dof_handler.n_dofs(), n_local_dofs, + n_local_dofs, dof_handler.max_couplings_between_dofs()); solution.reinit (mpi_communicator, dof_handler.n_dofs(), n_local_dofs); diff --git a/deal.II/lac/include/lac/block_matrix_base.h b/deal.II/lac/include/lac/block_matrix_base.h new file mode 100644 index 0000000000..2fe4402d2b --- /dev/null +++ b/deal.II/lac/include/lac/block_matrix_base.h @@ -0,0 +1,1668 @@ +//---------------------------- block_matrix_base.h --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2000, 2001, 2002, 2003, 2004 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- block_matrix_base.h --------------------------- +#ifndef __deal2__block_matrix_base_h +#define __deal2__block_matrix_base_h + + +#include +#include +#include +#include + +#include + + + +/*! @addtogroup Matrix1 + *@{ + */ + + +namespace internal +{ + +/** + * Namespace in which iterators in block matrices are implemented. + * + * @author Wolfgang Bangerth, 2004 + */ + namespace BlockMatrixIterators + { + template class ConstIterator; + + /** + * Accessor class for iterators + */ + template + class Accessor + { + public: + /** + * Typedef the value type of the + * matrix we point into. + */ + typedef typename BlockMatrix::value_type value_type; + + /** + * Constructor. Since we use + * accessors only for read + * access, a const matrix + * pointer is sufficient. + */ + Accessor (const BlockMatrix *m, + const unsigned int row, + const unsigned int index); + + /** + * Row number of the element + * represented by this + * object. + */ + unsigned int row() const; + + /** + * Index in row of the element + * represented by this + * object. + */ + unsigned int index() const; + + /** + * Column number of the + * element represented by + * this object. + */ + unsigned int column() const; + + /** + * Value of this matrix entry. + */ + value_type value() const; + + /** + * Block row of the + * element represented by + * this object. + */ + unsigned int block_row() const; + + /** + * Block column of the + * element represented by + * this object. + */ + unsigned int block_column() const; + + protected: + /** + * The matrix accessed. + */ + const BlockMatrix * matrix; + + /** + * Iterator of the underlying matrix + * class. + */ + typename BlockMatrix::BlockType::const_iterator base_iterator; + + /** + * Number of block where row lies in. + */ + unsigned int row_block; + + /** + * First row of block. + */ + unsigned int row_start; + + /** + * Number of block column where + * column lies in. + */ + unsigned int col_block; + + /** + * First column of block. + */ + unsigned int col_start; + + /** + * Index in whole row. + */ + unsigned int a_index; + + /** + * Let the iterator class be a + * friend. + */ + template + friend class ConstIterator; + }; + + /** + * STL conforming iterator. + */ + template + class ConstIterator + { + public: + /** + * Constructor. + */ + ConstIterator(const BlockMatrix *matrix, + const unsigned int row, + const unsigned int index); + + /** + * Prefix increment. + */ + ConstIterator& operator++ (); + + /** + * Postfix increment. + */ + ConstIterator& operator++ (int); + + /** + * Dereferencing operator. + */ + const Accessor & operator* () const; + + /** + * Dereferencing operator. + */ + const Accessor * operator-> () const; + + /** + * Comparison. True, if + * both iterators point to + * the same matrix + * position. + */ + bool operator == (const ConstIterator&) const; + /** + * Inverse of operator==(). + */ + bool operator != (const ConstIterator&) const; + + /** + * Comparison + * operator. Result is true + * if either the first row + * number is smaller or if + * the row numbers are + * equal and the first + * index is smaller. + */ + bool operator < (const ConstIterator&) const; + + private: + /** + * The accessor with which we work. + */ + Accessor accessor; + }; + } +} + + +/** + * Blocked matrix class. The behaviour of objects of this type is almost as + * for the usual matrix objects, with most of the functions being implemented + * in both classes. The main difference is that the matrix represented by this + * object is composed of an array of matrices (e.g. of type + * SparseMatrix) and all accesses to the elements of this object are + * relayed to accesses of the base matrices. The actual type of the individual + * blocks of this matrix is the type of the template argument, and can, for + * example be the usual SparseMatrix or PETScWrappers::SparseMatrix. + * + * In addition to the usual matrix access and linear algebra + * functions, there are functions block() which allow access to the + * different blocks of the matrix. This may, for example, be of help + * when you want to implement Schur complement methods, or block + * preconditioners, where each block belongs to a specific component + * of the equation you are presently discretizing. + * + * Note that the numbers of blocks and rows are implicitly determined + * by the sparsity pattern objects used. + * + * Objects of this type are frequently used when a system of differential + * equations has solutions with variables that fall into different + * classes. For example, solutions of the Stokes or Navier-Stokes equations + * have @p dim velocity components and one pressure component. In this case, + * it may make sense to consider the linear system of equations as a system of + * 2x2 blocks, and one can construct preconditioners or solvers based on this + * 2x2 block structure. This class can help you in these cases, as it allows + * to view the matrix alternatively as one big matrix, or as a number of + * individual blocks. + * + * + * @section Inheritance Inheriting from this class + * + * Since this class simply forwards its calls to the subobjects (if necessary + * after adjusting indices denoting which subobject is meant), this class is + * completely independent of the actual type of the subobject. The functions + * that set up block matrices and destroy them, however, have to be + * implemented in derived classes. These functions also have to fill the data + * members provided by this base class, as they are only used passively in + * this class. + * + * + * @ref Instantiations: some (@ @). Most of the + * functions take a vector or block vector argument. These functions can, in + * general, only successfully be compiled if the individual blocks of this + * matrix implement the respective functions operating on the vector type in + * question. For example, if you have a block sparse matrix over deal.II + * SparseMatrix objects, then you will likely not be able to form the + * matrix-vector multiplication with a block vector over + * PETScWrappers::SparseMatrix objects. If you attempt anyway, you will likely + * get a number of compiler errors. + * + * @author Wolfgang Bangerth, 2000, 2004 + */ +template +class BlockMatrixBase : public Subscriptor +{ + public: + /** + * Typedef the type of the underlying + * matrix. + */ + typedef MatrixType BlockType; + + /** + * Type of matrix entries. In analogy to + * the STL container classes. + */ + typedef typename BlockType::value_type value_type; + typedef value_type *pointer; + typedef const value_type *const_pointer; + typedef value_type &reference; + typedef const value_type &const_reference; + typedef size_t size_type; + + typedef + internal::BlockMatrixIterators::ConstIterator + iterator; + + typedef + internal::BlockMatrixIterators::ConstIterator + const_iterator; + + + /** + * Default constructor. + */ + BlockMatrixBase (); + + /** + * Copy the given matrix to this + * one. The operation throws an + * error if the sparsity patterns + * of the two involved matrices + * do not point to the same + * object, since in this case the + * copy operation is + * cheaper. Since this operation + * is notheless not for free, we + * do not make it available + * through operator=(), since + * this may lead to unwanted + * usage, e.g. in copy arguments + * to functions, which should + * really be arguments by + * reference. + * + * The source matrix may be a + * matrix of arbitrary type, as + * long as its data type is + * convertible to the data type + * of this matrix. + * + * The function returns a + * reference to this. + */ + template + BlockMatrixBase & + copy_from (const BlockMatrixType &source); + + /** + * Release all memory and return + * to a state just like after + * having called the default + * constructor. It also forgets + * the sparsity pattern it was + * previously tied to. + * + * This calls SparseMatrix::clear on all + * sub-matrices and then resets this + * object to have no blocks at all. + */ + void clear (); + + /** + * Access the block with the + * given coordinates. + */ + BlockType & + block (const unsigned int row, + const unsigned int column); + + + /** + * Access the block with the + * given coordinates. Version for + * constant objects. + */ + const BlockType & + block (const unsigned int row, + const unsigned int column) const; + + /** + * Return the dimension of the + * image space. To remember: the + * matrix is of dimension + * $m \times n$. + */ + unsigned int m () const; + + /** + * Return the dimension of the + * range space. To remember: the + * matrix is of dimension + * $m \times n$. + */ + unsigned int n () const; + + + /** + * Return the number of blocks in + * a column. Returns zero if no + * sparsity pattern is presently + * associated to this matrix. + */ + unsigned int n_block_rows () const; + + /** + * Return the number of blocks in + * a row. Returns zero if no + * sparsity pattern is presently + * associated to this matrix. + */ + unsigned int n_block_cols () const; + + /** + * Set the element (i,j) + * to value. Throws an + * error if the entry does not + * exist. Still, it is allowed to + * store zero values in + * non-existent fields. + */ + void set (const unsigned int i, + const unsigned int j, + const value_type value); + + /** + * Add value to the element + * (i,j). Throws an error if + * the entry does not + * exist. Still, it is allowed to + * store zero values in + * non-existent fields. + */ + void add (const unsigned int i, + const unsigned int j, + const value_type value); + + /** + * Return the value of the entry + * (i,j). This may be an + * expensive operation and you + * should always take care where + * to call this function. In + * order to avoid abuse, this + * function throws an exception + * if the wanted element does not + * exist in the matrix. + */ + value_type operator () (const unsigned int i, + const unsigned int j) const; + + /** + * This function is mostly like + * operator()() in that it + * returns the value of the + * matrix entry (i,j). The only + * difference is that if this + * entry does not exist in the + * sparsity pattern, then instead + * of raising an exception, zero + * is returned. While this may be + * convenient in some cases, note + * that it is simple to write + * algorithms that are slow + * compared to an optimal + * solution, since the sparsity + * of the matrix is not used. + */ + value_type el (const unsigned int i, + const unsigned int j) const; + + /** + * Call the compress() function on all + * the subblocks of the matrix. + */ + void compress (); + + /** + * Multiply the entire matrix by a + * fixed factor. + */ + BlockMatrixBase & operator *= (const value_type factor); + + /** + * Divide the entire matrix by a + * fixed factor. + */ + BlockMatrixBase & operator /= (const value_type factor); + + /** + * Add matrix scaled by + * factor to this matrix, + * i.e. the matrix factor*matrix + * is added to this. This + * function throws an error if the + * sparsity patterns of the two involved + * matrices do not point to the same + * object, since in this case the + * operation is cheaper. + * + * The source matrix may be a sparse + * matrix over an arbitrary underlying + * scalar type, as long as its data type + * is convertible to the data type of + * this matrix. + */ + template + void add_scaled (const value_type factor, + const BlockMatrixType &matrix); + + /** + * Matrix-vector multiplication: + * let $dst = M*src$ with $M$ + * being this matrix. + */ + template + void vmult (BlockVectorType &dst, + const BlockVectorType &src) const; + + /** + * Matrix-vector + * multiplication. Just like the + * previous function, but only + * applicable if the matrix has + * only one block column. + */ + template + void vmult (BlockVectorType &dst, + const VectorType &src) const; + + /** + * Matrix-vector + * multiplication. Just like the + * previous function, but only + * applicable if the matrix has + * only one block row. + */ + template + void vmult (VectorType &dst, + const BlockVectorType &src) const; + + /** + * Matrix-vector + * multiplication. Just like the + * previous function, but only + * applicable if the matrix has + * only one block. + */ + template + void vmult (VectorType &dst, + const VectorType &src) const; + + /** + * Matrix-vector multiplication: + * let $dst = M^T*src$ with $M$ + * being this matrix. This + * function does the same as + * vmult() but takes the + * transposed matrix. + */ + template + void Tvmult (BlockVectorType &dst, + const BlockVectorType &src) const; + + /** + * Matrix-vector + * multiplication. Just like the + * previous function, but only + * applicable if the matrix has + * only one block row. + */ + template + void Tvmult (BlockVectorType &dst, + const VectorType &src) const; + + /** + * Matrix-vector + * multiplication. Just like the + * previous function, but only + * applicable if the matrix has + * only one block column. + */ + template + void Tvmult (VectorType &dst, + const BlockVectorType &src) const; + + /** + * Matrix-vector + * multiplication. Just like the + * previous function, but only + * applicable if the matrix has + * only one block. + */ + template + void Tvmult (VectorType &dst, + const VectorType &src) const; + + /** + * Adding Matrix-vector + * multiplication. Add $M*src$ on + * $dst$ with $M$ being this + * matrix. + */ + template + void vmult_add (BlockVectorType &dst, + const BlockVectorType &src) const; + + /** + * Adding Matrix-vector + * multiplication. Add + * MTsrc to + * dst with M being + * this matrix. This function + * does the same as vmult_add() + * but takes the transposed + * matrix. + */ + template + void Tvmult_add (BlockVectorType &dst, + const BlockVectorType &src) const; + + /** + * Return the norm of the vector + * v with respect to the + * norm induced by this matrix, + * i.e. vTMv). This + * is useful, e.g. in the finite + * element context, where the + * LT-norm of a + * function equals the matrix + * norm with respect to the mass + * matrix of the vector + * representing the nodal values + * of the finite element + * function. Note that even + * though the function's name + * might suggest something + * different, for historic + * reasons not the norm but its + * square is returned, as defined + * above by the scalar product. + * + * Obviously, the matrix needs to + * be square for this operation. + */ + template + value_type + matrix_norm_square (const BlockVectorType &v) const; + + /** + * Compute the matrix scalar + * product $\left(u,Mv\right)$. + */ + template + value_type + matrix_scalar_product (const BlockVectorType &u, + const BlockVectorType &v) const; + + /** + * Compute the residual + * r=b-Ax. Write the + * residual into dst. + */ + template + value_type residual (BlockVectorType &dst, + const BlockVectorType &x, + const BlockVectorType &b) const; + + /** + * STL-like iterator with the + * first entry. + */ + const_iterator begin () const; + + /** + * Final iterator. + */ + const_iterator end () const; + + /** + * STL-like iterator with the + * first entry of row r. + */ + const_iterator begin (const unsigned int r) const; + + /** + * Final iterator of row r. + */ + const_iterator end (const unsigned int r) const; + + + /** + * Exception + */ + DeclException0 (ExcMatrixNotBlockSquare); + + /** + * Exception + */ + DeclException0 (ExcMatrixNotInitialized); + /** + * Exception + */ + DeclException4 (ExcIncompatibleRowNumbers, + int, int, int, int, + << "The blocks [" << arg1 << ',' << arg2 << "] and [" + << arg3 << ',' << arg4 << "] have differing row numbers."); + /** + * Exception + */ + DeclException4 (ExcIncompatibleColNumbers, + int, int, int, int, + << "The blocks [" << arg1 << ',' << arg2 << "] and [" + << arg3 << ',' << arg4 << "] have differing column numbers."); + + protected: + /** + * Index arrays for rows and columns. + */ + BlockIndices row_block_indices; + BlockIndices column_block_indices; + + /** + * Array of sub-matrices. + */ + Table<2,SmartPointer > sub_objects; + + /** + * This function collects the + * sizes of the sub-objects and + * stores them in internal + * arrays, in order to be able to + * relay global indices into the + * matrix to indices into the + * subobjects. You *must* call + * this function each time after + * you have changed the size of + * the sub-objects. + * + * Derived classes should call this + * function whenever the size of the + * sub-objects has changed and the @p + * X_block_indices arrays need to be + * updated. + * + * Note that this function is not public + * since not all derived classes need to + * export its interface. For example, for + * the usual deal.II SparseMatrix class, + * the sizes are implicitly determined + * whenever reinit() is called, and + * individual blocks cannot be + * resized. For that class, this function + * therefore does not have to be + * public. On the other hand, for the + * PETSc classes, there is no associated + * sparsity pattern object that + * determines the block sizes, and for + * these the function needs to be + * publicly available. These classes + * therefore export this function. + */ + void collect_sizes (); + + /** + * Make the iterator class a + * friend. We have to work around + * a compiler bug here again. + */ +#ifndef DEAL_II_NAMESP_TEMPL_FRIEND_BUG + template + friend class internal::BlockMatrixIterators::Accessor; + + template + friend class internal::BlockMatrixIterators::ConstIterator; +#else + typedef internal::BlockMatrixIterators::Accessor Accessor; + friend class Accessor; + + friend class const_iterator; +#endif +}; + + +/*@}*/ +/* ------------------------- Template functions ---------------------- */ + + +namespace internal +{ + namespace BlockMatrixIterators + { + template + inline + Accessor:: + Accessor (const BlockMatrix *matrix, + const unsigned int r, + const unsigned int i) + : + matrix(matrix), + base_iterator(matrix->block(0,0).begin()), + row_block(0), + row_start(0), + col_block(0), + col_start(0), + a_index(0) + { + Assert (i==0, ExcNotImplemented()); + + if (r < matrix->m()) + { + std::pair indices + = matrix->row_block_indices.global_to_local(r); + row_block = indices.first; + base_iterator = matrix->block(indices.first, 0).begin(indices.second); + row_start = matrix->row_block_indices.local_to_global(row_block, 0); + } + else + { + row_block = matrix->n_block_rows(); + base_iterator = matrix->block(0, 0).begin(); + } + } + + + template + inline + unsigned int + Accessor::row() const + { + return row_start + base_iterator->row(); + } + + + template + inline + unsigned int + Accessor::index() const + { + return a_index; + } + + + template + inline + unsigned int + Accessor::column() const + { + return col_start + base_iterator->column(); + } + + + template + inline + unsigned int + Accessor::block_row() const + { + return row_block; + } + + + template + inline + unsigned int + Accessor::block_column() const + { + return col_block; + } + + + template + inline + typename Accessor::value_type + Accessor::value () const + { + return base_iterator->value(); + } + + +//----------------------------------------------------------------------// + + + template + inline + ConstIterator:: + ConstIterator (const BlockMatrix *m, + const unsigned int r, + const unsigned int i) + : + accessor (m, r, i) + {} + + + + template + inline + ConstIterator & + ConstIterator::operator++ () + { + Assert (accessor.row_block < accessor.matrix->n_block_rows(), + ExcIteratorPastEnd()); + + // Remember current row inside block + unsigned int local_row = accessor.base_iterator->row(); + + // Advance inside block + ++accessor.base_iterator; + ++accessor.a_index; + + // If end of row inside block, + // advance to next block + if (accessor.base_iterator == + accessor.matrix->block(accessor.row_block, + accessor.col_block).end(local_row)) + { + if (accessor.col_block < accessor.matrix->n_block_cols()-1) + { + // Advance to next block in + // row +//TODO: what if this row in that block is empty? + ++accessor.col_block; + accessor.col_start = + accessor.matrix->column_block_indices + .local_to_global(accessor.col_block, 0); + } + else + { + // Advance to first block + // in next row + accessor.col_block = 0; + accessor.col_start = 0; + accessor.a_index = 0; + ++local_row; + if (local_row >= + accessor.matrix->block(accessor.row_block,0).m()) + { + // If final row in + // block, go to next + // block row + local_row = 0; + ++accessor.row_block; + if (accessor.row_block < accessor.matrix->n_block_rows()) + accessor.row_start = + accessor.matrix->row_block_indices + .local_to_global(accessor.row_block, 0); + } + } + // Finally, set base_iterator + // to start of row determined + // above + if (accessor.row_block < accessor.matrix->n_block_rows()) + accessor.base_iterator = + accessor.matrix->block(accessor.row_block, accessor.col_block) + .begin(local_row); + else + // Set base_iterator to a + // defined state for + // comparison. This is the + // end() state. +//TODO: this is a particularly bad choice, since it is actually a valid iterator. it should rather be something like the end iterator of the last block! + accessor.base_iterator = accessor.matrix->block(0, 0).begin(); + } + return *this; + } + + + +// template +// inline +// const_iterator& +// ConstIterator::operator++ (int) +// { +// Assert (false, ExcNotImplemented()); +// } + + + + + template + inline + const Accessor & + ConstIterator::operator* () const + { + return accessor; + } + + + template + inline + const Accessor * + ConstIterator::operator-> () const + { + return &accessor; + } + + + + template + inline + bool + ConstIterator:: + operator == (const ConstIterator& i) const + { + if (accessor.matrix != i.accessor.matrix) + return false; + + if (accessor.row_block == i.accessor.row_block + && accessor.col_block == i.accessor.col_block + && accessor.base_iterator == i.accessor.base_iterator) + return true; + return false; + } + + + + template + inline + bool + ConstIterator:: + operator != (const ConstIterator& i) const + { + return !(*this == i); + } + + + + template + inline + bool + ConstIterator:: + operator < (const ConstIterator& i) const + { + if (accessor.row_block < i.accessor.row_block) + return true; + if (accessor.row_block == i.accessor.row_block) + { + if (accessor.base_iterator->row() < i.accessor.base_iterator->row()) + return true; + if (accessor.base_iterator->row() == i.accessor.base_iterator->row()) + if (accessor.a_index < i.accessor.a_index) + return true; + } + return false; + } + + } +} + +//----------------------------------------------------------------------// + + +template +inline +BlockMatrixBase::BlockMatrixBase () +{} + + + +template +template +inline +BlockMatrixBase & +BlockMatrixBase:: +copy_from (const BlockMatrixType &source) +{ + for (unsigned int r=0; r +inline +void +BlockMatrixBase::clear () +{ + for (unsigned int r=0; r +inline +typename BlockMatrixBase::BlockType & +BlockMatrixBase::block (const unsigned int row, + const unsigned int column) +{ + Assert (row +inline +const typename BlockMatrixBase::BlockType & +BlockMatrixBase::block (const unsigned int row, + const unsigned int column) const +{ + Assert (row +inline +unsigned int +BlockMatrixBase::m () const +{ + return row_block_indices.total_size(); +} + + + +template +inline +unsigned int +BlockMatrixBase::n () const +{ + return column_block_indices.total_size(); +} + + + +template +inline +unsigned int +BlockMatrixBase::n_block_cols () const +{ + return column_block_indices.size(); +} + + + +template +inline +unsigned int +BlockMatrixBase::n_block_rows () const +{ + return row_block_indices.size(); +} + + + +template +inline +void +BlockMatrixBase::set (const unsigned int i, + const unsigned int j, + const value_type value) +{ + const std::pair + row_index = row_block_indices.global_to_local (i), + col_index = column_block_indices.global_to_local (j); + block(row_index.first,col_index.first).set (row_index.second, + col_index.second, + value); +} + + + +template +inline +void +BlockMatrixBase::add (const unsigned int i, + const unsigned int j, + const value_type value) +{ + const std::pair + row_index = row_block_indices.global_to_local (i), + col_index = column_block_indices.global_to_local (j); + block(row_index.first,col_index.first).add (row_index.second, + col_index.second, + value); +} + + + +template +inline +typename BlockMatrixBase::value_type +BlockMatrixBase::operator () (const unsigned int i, + const unsigned int j) const +{ + const std::pair + row_index = row_block_indices.global_to_local (i), + col_index = column_block_indices.global_to_local (j); + return block(row_index.first,col_index.first) (row_index.second, + col_index.second); +} + + + +template +inline +typename BlockMatrixBase::value_type +BlockMatrixBase::el (const unsigned int i, + const unsigned int j) const +{ + const std::pair + row_index = row_block_indices.global_to_local (i), + col_index = column_block_indices.global_to_local (j); + return block(row_index.first,col_index.first).el (row_index.second, + col_index.second); +} + + + +template +inline +void +BlockMatrixBase::compress () +{ + for (unsigned int r=0; r +inline +BlockMatrixBase & +BlockMatrixBase::operator *= (const value_type factor) +{ + Assert (n_block_cols() != 0, ExcMatrixNotInitialized()); + Assert (n_block_rows() != 0, ExcMatrixNotInitialized()); + + for (unsigned int r=0; r +inline +BlockMatrixBase & +BlockMatrixBase::operator /= (const value_type factor) +{ + Assert (n_block_cols() != 0, ExcMatrixNotInitialized()); + Assert (n_block_rows() != 0, ExcMatrixNotInitialized()); + Assert (factor !=0, ExcDivideByZero()); + + const value_type factor_inv = 1. / factor; + + for (unsigned int r=0; r +template +void +BlockMatrixBase:: +add_scaled (const value_type factor, + const BlockMatrixType &matrix) +{ + for (unsigned int r=0; r +template +void +BlockMatrixBase::vmult (BlockVectorType &dst, + const BlockVectorType &src) const +{ + Assert (dst.n_blocks() == n_block_rows(), + ExcDimensionMismatch(dst.n_blocks(), n_block_rows())); + Assert (src.n_blocks() == n_block_cols(), + ExcDimensionMismatch(src.n_blocks(), n_block_cols())); + + for (unsigned int row=0; row +template +void +BlockMatrixBase::vmult (VectorType &dst, + const BlockVectorType &src) const +{ + Assert (n_block_rows() == 1, + ExcDimensionMismatch(1, n_block_rows())); + Assert (src.n_blocks() == n_block_cols(), + ExcDimensionMismatch(src.n_blocks(), n_block_cols())); + + block(0,0).vmult (dst, src.block(0)); + for (unsigned int col=1; col +template +void +BlockMatrixBase::vmult (BlockVectorType &dst, + const VectorType &src) const +{ + Assert (dst.n_blocks() == n_block_rows(), + ExcDimensionMismatch(dst.n_blocks(), n_block_rows())); + Assert (1 == n_block_cols(), + ExcDimensionMismatch(1, n_block_cols())); + + for (unsigned int row=0; row +template +void +BlockMatrixBase::vmult (VectorType &dst, + const VectorType &src) const +{ + Assert (1 == n_block_rows(), + ExcDimensionMismatch(1, n_block_rows())); + Assert (1 == n_block_cols(), + ExcDimensionMismatch(1, n_block_cols())); + + block(0,0).vmult (dst, src); +} + + + +template +template +void +BlockMatrixBase::vmult_add (BlockVectorType &dst, + const BlockVectorType &src) const +{ + Assert (dst.n_blocks() == n_block_rows(), + ExcDimensionMismatch(dst.n_blocks(), n_block_rows())); + Assert (src.n_blocks() == n_block_cols(), + ExcDimensionMismatch(src.n_blocks(), n_block_cols())); + + for (unsigned int row=0; row +template +void +BlockMatrixBase::Tvmult (BlockVectorType &dst, + const BlockVectorType &src) const +{ + Assert (dst.n_blocks() == n_block_cols(), + ExcDimensionMismatch(dst.n_blocks(), n_block_cols())); + Assert (src.n_blocks() == n_block_rows(), + ExcDimensionMismatch(src.n_blocks(), n_block_rows())); + + dst = 0.; + + for (unsigned int row=0; row +template +void +BlockMatrixBase::Tvmult (BlockVectorType &dst, + const VectorType &src) const +{ + Assert (dst.n_blocks() == n_block_cols(), + ExcDimensionMismatch(dst.n_blocks(), n_block_cols())); + Assert (1 == n_block_rows(), + ExcDimensionMismatch(1, n_block_rows())); + + dst = 0.; + + for (unsigned int col=0; col +template +void +BlockMatrixBase::Tvmult (VectorType &dst, + const BlockVectorType &src) const +{ + Assert (1 == n_block_cols(), + ExcDimensionMismatch(1, n_block_cols())); + Assert (src.n_blocks() == n_block_rows(), + ExcDimensionMismatch(src.n_blocks(), n_block_rows())); + + block(0,0).Tvmult (dst, src.block(0)); + + for (unsigned int row=1; row +template +void +BlockMatrixBase::Tvmult (VectorType &dst, + const VectorType &src) const +{ + Assert (1 == n_block_cols(), + ExcDimensionMismatch(1, n_block_cols())); + Assert (1 == n_block_rows(), + ExcDimensionMismatch(1, n_block_rows())); + + block(0,0).Tvmult (dst, src); +} + + + +template +template +void +BlockMatrixBase::Tvmult_add (BlockVectorType& dst, + const BlockVectorType& src) const +{ + Assert (dst.n_blocks() == n_block_cols(), + ExcDimensionMismatch(dst.n_blocks(), n_block_cols())); + Assert (src.n_blocks() == n_block_rows(), + ExcDimensionMismatch(src.n_blocks(), n_block_rows())); + + for (unsigned int row=0; row +template +typename BlockMatrixBase::value_type +BlockMatrixBase::matrix_norm_square (const BlockVectorType &v) const +{ + Assert (n_block_rows() == n_block_cols(), + ExcMatrixNotBlockSquare()); + Assert (v.n_blocks() == n_block_rows(), + ExcDimensionMismatch(v.n_blocks(), n_block_rows())); + + value_type norm_sqr = 0; + for (unsigned int row=0; row +template +typename BlockMatrixBase::value_type +BlockMatrixBase:: +matrix_scalar_product (const BlockVectorType &u, + const BlockVectorType &v) const +{ + Assert (u.n_blocks() == n_block_rows(), + ExcDimensionMismatch(u.n_blocks(), n_block_rows())); + Assert (v.n_blocks() == n_block_cols(), + ExcDimensionMismatch(v.n_blocks(), n_block_cols())); + + value_type result = 0; + for (unsigned int row=0; row +template +typename BlockMatrixBase::value_type +BlockMatrixBase:: +residual (BlockVectorType &dst, + const BlockVectorType &x, + const BlockVectorType &b) const +{ + Assert (dst.n_blocks() == n_block_rows(), + ExcDimensionMismatch(dst.n_blocks(), n_block_rows())); + Assert (b.n_blocks() == n_block_rows(), + ExcDimensionMismatch(b.n_blocks(), n_block_rows())); + Assert (x.n_blocks() == n_block_cols(), + ExcDimensionMismatch(x.n_blocks(), n_block_cols())); + // in block notation, the residual is + // r_i = b_i - \sum_j A_ij x_j. + // this can be written as + // r_i = b_i - A_i0 x_0 - \sum_{j>0} A_ij x_j. + // + // for the first two terms, we can + // call the residual function of + // A_i0. for the other terms, we + // use vmult_add. however, we want + // to subtract, so in order to + // avoid a temporary vector, we + // perform a sign change of the + // first two term before, and after + // adding up + for (unsigned int row=0; row +inline +typename BlockMatrixBase::const_iterator +BlockMatrixBase::begin () const +{ + return const_iterator(this, 0, 0); +} + + + +template +inline +typename BlockMatrixBase::const_iterator +BlockMatrixBase::end () const +{ + return const_iterator(this, m(), 0); +} + + + +template +inline +typename BlockMatrixBase::const_iterator +BlockMatrixBase::begin (const unsigned int r) const +{ + Assert (r +inline +typename BlockMatrixBase::const_iterator +BlockMatrixBase::end (const unsigned int r) const +{ + Assert (r +void +BlockMatrixBase::collect_sizes () +{ + std::vector row_sizes (this->n_block_rows()); + std::vector col_sizes (this->n_block_cols()); + + // first find out the row sizes + // from the first block column + for (unsigned int r=0; rn_block_rows(); ++r) + row_sizes[r] = sub_objects[r][0]->m(); + // then check that the following + // block columns have the same + // sizes + for (unsigned int c=1; cn_block_cols(); ++c) + for (unsigned int r=0; rn_block_rows(); ++r) + Assert (row_sizes[r] == sub_objects[r][c]->m(), + ExcIncompatibleRowNumbers (r,0,r,c)); + + // finally initialize the row + // indices with this array + this->row_block_indices.reinit (row_sizes); + + + // then do the same with the columns + for (unsigned int c=0; cn_block_cols(); ++c) + col_sizes[c] = sub_objects[0][c]->n(); + for (unsigned int r=1; rn_block_rows(); ++r) + for (unsigned int c=0; cn_block_cols(); ++c) + Assert (col_sizes[c] == sub_objects[r][c]->n(), + ExcIncompatibleRowNumbers (0,c,r,c)); + + // finally initialize the row + // indices with this array + this->column_block_indices.reinit (col_sizes); +} + + + +#endif // __deal2__block_matrix_base_h diff --git a/deal.II/lac/include/lac/block_sparse_matrix.h b/deal.II/lac/include/lac/block_sparse_matrix.h index 7c2b4346a5..9b834ce760 100644 --- a/deal.II/lac/include/lac/block_sparse_matrix.h +++ b/deal.II/lac/include/lac/block_sparse_matrix.h @@ -16,267 +16,57 @@ #include #include +#include #include #include #include -template class Vector; -template class BlockVector; -template class BlockSparseMatrix; - /*! @addtogroup Matrix1 *@{ */ -namespace internal -{ - -/** - * Namespace in which iterators in block matrices are implemented. - * - * @author Wolfgang Bangerth, 2004 - */ - namespace BlockMatrixIterators - { - template class ConstIterator; - - /** - * Accessor class for iterators - */ - template - class Accessor - { - public: - /** - * Typedef the value type of the - * matrix we point into. - */ - typedef typename BlockMatrix::value_type value_type; - - /** - * Constructor. Since we use - * accessors only for read - * access, a const matrix - * pointer is sufficient. - */ - Accessor (const BlockMatrix *m, - const unsigned int row, - const unsigned int index); - - /** - * Row number of the element - * represented by this - * object. - */ - unsigned int row() const; - - /** - * Index in row of the element - * represented by this - * object. - */ - unsigned int index() const; - - /** - * Column number of the - * element represented by - * this object. - */ - unsigned int column() const; - - /** - * Value of this matrix entry. - */ - value_type value() const; - - /** - * Block row of the - * element represented by - * this object. - */ - unsigned int block_row() const; - - /** - * Block column of the - * element represented by - * this object. - */ - unsigned int block_column() const; - - protected: - /** - * The matrix accessed. - */ - const BlockMatrix * matrix; - - /** - * Iterator of the underlying matrix - * class. - */ - typename BlockMatrix::BlockType::const_iterator base_iterator; - - /** - * Number of block where row lies in. - */ - unsigned int row_block; - - /** - * First row of block. - */ - unsigned int row_start; - - /** - * Number of block column where - * column lies in. - */ - unsigned int col_block; - - /** - * First column of block. - */ - unsigned int col_start; - - /** - * Index in whole row. - */ - unsigned int a_index; - - /** - * Let the iterator class be a - * friend. - */ - template - friend class ConstIterator; - }; - - /** - * STL conforming iterator. - */ - template - class ConstIterator : private Accessor - { - public: - /** - * Constructor. - */ - ConstIterator(const BlockMatrix *matrix, - const unsigned int row, - const unsigned int index); - - /** - * Prefix increment. - */ - ConstIterator& operator++ (); - - /** - * Postfix increment. - */ - ConstIterator& operator++ (int); - - /** - * Dereferencing operator. - */ - const Accessor & operator* () const; - - /** - * Dereferencing operator. - */ - const Accessor * operator-> () const; - - /** - * Comparison. True, if - * both iterators point to - * the same matrix - * position. - */ - bool operator == (const ConstIterator&) const; - /** - * Inverse of operator==(). - */ - bool operator != (const ConstIterator&) const; - - /** - * Comparison - * operator. Result is true - * if either the first row - * number is smaller or if - * the row numbers are - * equal and the first - * index is smaller. - */ - bool operator < (const ConstIterator&) const; - }; - } -} - - /** - * Blocked sparse matrix. The behaviour of objects of this type is - * almost as for the SparseMatrix objects, with most of the - * functions being implemented in both classes. The main difference is - * that the matrix represented by this object is composed of an array - * of sparse matrices (i.e. of type SparseMatrix) and all - * accesses to the elements of this object are relayed to accesses of - * the base matrices. - * - * In addition to the usual matrix access and linear algebra - * functions, there are functions block() which allow access to the - * different blocks of the matrix. This may, for example, be of help - * when you want to implement Schur complement methods, or block - * preconditioners, where each block belongs to a specific component - * of the equation you are presently discretizing. - * - * Note that the numbers of blocks and rows are implicitly determined - * by the sparsity pattern objects used. - * - * Objects of this type are frequently used when a system of differential - * equations has solutions with variables that fall into different - * classes. For example, solutions of the Stokes or Navier-Stokes equations - * have @p dim velocity components and one pressure component. In this case, - * it may make sense to consider the linear system of equations as a system of - * 2x2 blocks, and one can construct preconditioners or solvers based on this - * 2x2 block structure. This class can help you in these cases, as it allows - * to view the matrix alternatively as one big matrix, or as a number of - * individual blocks. - * - * @ref Instantiations: some (@ @) + * Blocked sparse matrix based on the SparseMatrix class. This class + * implements the functions that are specific to the SparseMatrix base objects + * for a blocked sparse matrix, and leaves the actual work relaying most of + * the calls to the individual blocks to the functions implemented in the base + * class. See there also for a description of when this class is useful. * - * @author Wolfgang Bangerth, 2000 + * @author Wolfgang Bangerth, 2000, 2004 */ template -class BlockSparseMatrix : public Subscriptor +class BlockSparseMatrix : public BlockMatrixBase > { public: + /** + * Typedef the base class for simpler + * access to its own typedefs. + */ + typedef BlockMatrixBase > BaseClass; + /** * Typedef the type of the underlying * matrix. */ - typedef SparseMatrix BlockType; + typedef typename BaseClass::BlockType BlockType; /** - * Type of matrix entries. In analogy to - * the STL container classes. + * Import the typedefs from the base + * class. */ - typedef typename BlockType::value_type value_type; - typedef value_type *pointer; - typedef const value_type *const_pointer; - typedef value_type &reference; - typedef const value_type &const_reference; - typedef size_t size_type; - - typedef - internal::BlockMatrixIterators::ConstIterator - iterator; - - typedef - internal::BlockMatrixIterators::ConstIterator - const_iterator; - + typedef typename BaseClass::value_type value_type; + typedef typename BaseClass::pointer pointer; + typedef typename BaseClass::const_pointer const_pointer; + typedef typename BaseClass::reference reference; + typedef typename BaseClass::const_reference const_reference; + typedef typename BaseClass::size_type size_type; + typedef typename BaseClass::iterator iterator; + typedef typename BaseClass::const_iterator const_iterator; - /** + /** * Constructor; initializes the * matrix to be empty, without * any structure, i.e. the @@ -330,13 +120,12 @@ class BlockSparseMatrix : public Subscriptor /** - * Pseudo operator only copying - * empty objects. The sizes of - * the block matrices need to be - * the same. + * Pseudo copy operator only copying + * empty objects. The sizes of the block + * matrices need to be the same. */ - BlockSparseMatrix & - operator = (const BlockSparseMatrix &); + BlockSparseMatrix & + operator = (const BlockSparseMatrix &); /** * This operator assigns a scalar to a @@ -375,54 +164,6 @@ class BlockSparseMatrix : public Subscriptor */ virtual void reinit (const BlockSparsityPattern &sparsity); - - /** - * Access the block with the - * given coordinates. - */ - BlockType & - block (const unsigned int row, - const unsigned int column); - - - /** - * Access the block with the - * given coordinates. Version for - * constant objects. - */ - const BlockType & - block (const unsigned int row, - const unsigned int column) const; - - /** - * Release all memory and return - * to a state just like after - * having called the default - * constructor. It also forgets - * the sparsity pattern it was - * previously tied to. - * - * This calls SparseMatrix::clear on all - * sub-matrices. - */ - virtual void clear (); - - /** - * Return the number of blocks in - * a column. Returns zero if no - * sparsity pattern is presently - * associated to this matrix. - */ - unsigned int n_block_rows () const; - - /** - * Return the number of blocks in - * a row. Returns zero if no - * sparsity pattern is presently - * associated to this matrix. - */ - unsigned int n_block_cols () const; - /** * Return whether the object is * empty. It is empty if either @@ -432,22 +173,6 @@ class BlockSparseMatrix : public Subscriptor */ bool empty () const; - /** - * Return the dimension of the - * image space. To remember: the - * matrix is of dimension - * $m \times n$. - */ - unsigned int m () const; - - /** - * Return the dimension of the - * range space. To remember: the - * matrix is of dimension - * $m \times n$. - */ - unsigned int n () const; - /** * Return the number of nonzero * elements of this @@ -465,287 +190,7 @@ class BlockSparseMatrix : public Subscriptor * the number of actually nonzero * elements of all the blocks. */ - unsigned int n_actually_nonzero_elements () const; - - /** - * Set the element (i,j) - * to value. Throws an - * error if the entry does not - * exist. Still, it is allowed to - * store zero values in - * non-existent fields. - */ - void set (const unsigned int i, - const unsigned int j, - const number value); - - /** - * Add value to the element - * (i,j). Throws an error if - * the entry does not - * exist. Still, it is allowed to - * store zero values in - * non-existent fields. - */ - void add (const unsigned int i, - const unsigned int j, - const number value); - - /** - * Multiply the entire matrix by a - * fixed factor. - */ - BlockSparseMatrix & operator *= (const number factor); - - /** - * Divide the entire matrix by a - * fixed factor. - */ - BlockSparseMatrix & operator /= (const number factor); - - /** - * Copy the given matrix to this - * one. The operation throws an - * error if the sparsity patterns - * of the two involved matrices - * do not point to the same - * object, since in this case the - * copy operation is - * cheaper. Since this operation - * is notheless not for free, we - * do not make it available - * through operator=(), since - * this may lead to unwanted - * usage, e.g. in copy arguments - * to functions, which should - * really be arguments by - * reference. - * - * The source matrix may be a - * matrix of arbitrary type, as - * long as its data type is - * convertible to the data type - * of this matrix. - * - * The function returns a - * reference to this. - */ - template - BlockSparseMatrix & - copy_from (const BlockSparseMatrix &source); - - /** - * Add matrix scaled by - * factor to this matrix, - * i.e. the matrix factor*matrix - * is added to this. This - * function throws an error if the - * sparsity patterns of the two involved - * matrices do not point to the same - * object, since in this case the - * operation is cheaper. - * - * The source matrix may be a sparse - * matrix over an arbitrary underlying - * scalar type, as long as its data type - * is convertible to the data type of - * this matrix. - */ - template - void add_scaled (const number factor, - const BlockSparseMatrix &matrix); - - /** - * Return the value of the entry - * (i,j). This may be an - * expensive operation and you - * should always take care where - * to call this function. In - * order to avoid abuse, this - * function throws an exception - * if the wanted element does not - * exist in the matrix. - */ - number operator () (const unsigned int i, - const unsigned int j) const; - - /** - * This function is mostly like - * operator()() in that it - * returns the value of the - * matrix entry (i,j). The only - * difference is that if this - * entry does not exist in the - * sparsity pattern, then instead - * of raising an exception, zero - * is returned. While this may be - * convenient in some cases, note - * that it is simple to write - * algorithms that are slow - * compared to an optimal - * solution, since the sparsity - * of the matrix is not used. - */ - number el (const unsigned int i, - const unsigned int j) const; - - /** - * Matrix-vector multiplication: - * let $dst = M*src$ with $M$ - * being this matrix. - */ - template - void vmult (BlockVector &dst, - const BlockVector &src) const; - - /** - * Matrix-vector - * multiplication. Just like the - * previous function, but only - * applicable if the matrix has - * only one block column. - */ - template - void vmult (BlockVector &dst, - const Vector &src) const; - - /** - * Matrix-vector - * multiplication. Just like the - * previous function, but only - * applicable if the matrix has - * only one block row. - */ - template - void vmult (Vector &dst, - const BlockVector &src) const; - - /** - * Matrix-vector - * multiplication. Just like the - * previous function, but only - * applicable if the matrix has - * only one block. - */ - template - void vmult (Vector &dst, - const Vector &src) const; - - /** - * Matrix-vector multiplication: - * let $dst = M^T*src$ with $M$ - * being this matrix. This - * function does the same as - * vmult() but takes the - * transposed matrix. - */ - template - void Tvmult (BlockVector &dst, - const BlockVector &src) const; - - /** - * Matrix-vector - * multiplication. Just like the - * previous function, but only - * applicable if the matrix has - * only one block row. - */ - template - void Tvmult (BlockVector &dst, - const Vector &src) const; - - /** - * Matrix-vector - * multiplication. Just like the - * previous function, but only - * applicable if the matrix has - * only one block column. - */ - template - void Tvmult (Vector &dst, - const BlockVector &src) const; - - /** - * Matrix-vector - * multiplication. Just like the - * previous function, but only - * applicable if the matrix has - * only one block. - */ - template - void Tvmult (Vector &dst, - const Vector &src) const; - - /** - * Adding Matrix-vector - * multiplication. Add $M*src$ on - * $dst$ with $M$ being this - * matrix. - */ - template - void vmult_add (BlockVector &dst, - const BlockVector &src) const; - - /** - * Adding Matrix-vector - * multiplication. Add - * MTsrc to - * dst with M being - * this matrix. This function - * does the same as vmult_add() - * but takes the transposed - * matrix. - */ - template - void Tvmult_add (BlockVector &dst, - const BlockVector &src) const; - - /** - * Return the norm of the vector - * v with respect to the - * norm induced by this matrix, - * i.e. vTMv). This - * is useful, e.g. in the finite - * element context, where the - * LT-norm of a - * function equals the matrix - * norm with respect to the mass - * matrix of the vector - * representing the nodal values - * of the finite element - * function. Note that even - * though the function's name - * might suggest something - * different, for historic - * reasons not the norm but its - * square is returned, as defined - * above by the scalar product. - * - * Obviously, the matrix needs to - * be square for this operation. - */ - template - somenumber - matrix_norm_square (const BlockVector &v) const; - - /** - * Compute the matrix scalar - * product $\left(u,Mv\right)$. - */ - template - somenumber - matrix_scalar_product (const BlockVector &u, - const BlockVector &v) const; - - /** - * Compute the residual - * r=b-Ax. Write the - * residual into dst. - */ - template - somenumber residual (BlockVector &dst, - const BlockVector &x, - const BlockVector &b) const; + unsigned int n_actually_nonzero_elements () const; /** * Apply the Jacobi @@ -762,39 +207,11 @@ class BlockSparseMatrix : public Subscriptor * square matrices for this * operation. */ - template - void precondition_Jacobi (BlockVector &dst, - const BlockVector &src, - const number omega = 1.) const; + template + void precondition_Jacobi (BlockVectorType &dst, + const BlockVectorType &src, + const number omega = 1.) const; - /** - * Apply the Jacobi - * preconditioner, which - * multiplies every element of - * the src vector by the - * inverse of the respective - * diagonal element and - * multiplies the result with the - * relaxation parameter - * omega. - * - * All diagonal blocks must be - * square matrices for this - * operation. - * - * In contrast to the previous - * function, input and output - * vectors are not block - * vectors. This operation is - * thus only allowed if the block - * matrix consists of only one - * block. - */ - template - void precondition_Jacobi (Vector &dst, - const Vector &src, - const number omega = 1.) const; - /** * Print the matrix in the usual * format, i.e. as a matrix and @@ -857,28 +274,6 @@ class BlockSparseMatrix : public Subscriptor */ const BlockSparsityPattern & get_sparsity_pattern () const; - - /** - * STL-like iterator with the - * first entry. - */ - const_iterator begin () const; - - /** - * Final iterator. - */ - const_iterator end () const; - - /** - * STL-like iterator with the - * first entry of row r. - */ - const_iterator begin (const unsigned int r) const; - - /** - * Final iterator of row r. - */ - const_iterator end (const unsigned int r) const; /** * Determine an estimate for the @@ -887,37 +282,12 @@ class BlockSparseMatrix : public Subscriptor */ unsigned int memory_consumption () const; - /** - * Exception - */ - DeclException0 (ExcMatrixNotBlockSquare); - - /** - * Exception - */ - DeclException0 (ExcMatrixNotInitialized); + /** + * Exception + */ + DeclException0 (ExcBlockDimensionMismatch); private: - /** - * Number of block rows in this - * matrix. Is set by default to - * zero, and is only changed if a - * sparsity pattern is given to - * the constructor or the - * reinit() function. - */ - unsigned int rows; - - /** - * Number of block columns in this - * matrix. Is set by default to - * zero, and is only changed if a - * sparsity pattern is given to - * the constructor or the - * reinit() function. - */ - unsigned int columns; - /** * Pointer to the block sparsity * pattern used for this @@ -928,29 +298,6 @@ class BlockSparseMatrix : public Subscriptor * class. */ SmartPointer sparsity_pattern; - - /** - * Array of sub-matrices. - */ - Table<2,SmartPointer > sub_objects; - - /** - * Make the iterator class a - * friend. We have to work around - * a compiler bug here again. - */ -#ifndef DEAL_II_NAMESP_TEMPL_FRIEND_BUG - template - friend class internal::BlockMatrixIterators::Accessor; - - template - friend class internal::BlockMatrixIterators::ConstIterator; -#else - typedef internal::BlockMatrixIterators::Accessor Accessor; - friend class Accessor; - - friend class const_iterator; -#endif }; @@ -959,837 +306,30 @@ class BlockSparseMatrix : public Subscriptor /* ------------------------- Template functions ---------------------- */ -namespace internal -{ - namespace BlockMatrixIterators - { - template - inline - Accessor:: - Accessor (const BlockMatrix *matrix, - const unsigned int r, - const unsigned int i) - : - matrix(matrix), - base_iterator(matrix->block(0,0).begin()), - row_block(0), - row_start(0), - col_block(0), - col_start(0), - a_index(0) - { - Assert (i==0, ExcNotImplemented()); - - if (r < matrix->m()) - { - std::pair indices - = matrix->sparsity_pattern->get_row_indices().global_to_local(r); - row_block = indices.first; - base_iterator = matrix->block(indices.first, 0).begin(indices.second); - row_start = matrix->sparsity_pattern - ->get_row_indices().local_to_global(row_block, 0); - } - else - { - row_block = matrix->n_block_rows(); - base_iterator = matrix->block(0, 0).begin(); - } - } - - - template - inline - unsigned int - Accessor::row() const - { - return row_start + base_iterator->row(); - } - - - template - inline - unsigned int - Accessor::index() const - { - return a_index; - } - - - template - inline - unsigned int - Accessor::column() const - { - return col_start + base_iterator->column(); - } - - - template - inline - unsigned int - Accessor::block_row() const - { - return row_block; - } - - - template - inline - unsigned int - Accessor::block_column() const - { - return col_block; - } - - - template - inline - typename Accessor::value_type - Accessor::value () const - { - return base_iterator->value(); - } - - -//----------------------------------------------------------------------// - - - template - inline - ConstIterator:: - ConstIterator (const BlockMatrix *m, - const unsigned int r, - const unsigned int i) - : - Accessor(m, r, i) - {} - - - - template - inline - ConstIterator & - ConstIterator::operator++ () - { - Assert (this->row_blockmatrix->n_block_rows(), ExcIteratorPastEnd()); - - // Remember current row inside block - unsigned int local_row = this->base_iterator->row(); - - // Advance inside block - ++this->base_iterator; - ++this->a_index; - - // If end of row inside block, - // advance to next block - if (this->base_iterator == - this->matrix->block(this->row_block, this->col_block).end(local_row)) - { - if (this->col_block < this->matrix->n_block_cols()-1) - { - // Advance to next block in - // row -//TODO: what if this row in that block is empty? - ++this->col_block; - this->col_start = this->matrix->sparsity_pattern - ->get_column_indices().local_to_global(this->col_block, 0); - } - else - { - // Advance to first block - // in next row - this->col_block = 0; - this->col_start = 0; - this->a_index = 0; - ++local_row; - if (local_row >= this->matrix->block(this->row_block,0).m()) - { - // If final row in - // block, go to next - // block row - local_row = 0; - ++this->row_block; - if (this->row_block < this->matrix->n_block_rows()) - this->row_start = this->matrix->sparsity_pattern - ->get_row_indices().local_to_global(this->row_block, 0); - } - } - // Finally, set base_iterator - // to start of row determined - // above - if (this->row_block < this->matrix->n_block_rows()) - this->base_iterator = this->matrix->block(this->row_block, this->col_block).begin(local_row); - else - // Set base_iterator to a - // defined state for - // comparison. This is the - // end() state. -//TODO: this is a particularly bad choice, since it is actually a valid iterator. it should rather be something like the end iterator of the last block! - this->base_iterator = this->matrix->block(0, 0).begin(); - } - return *this; - } - - - -// template -// inline -// const_iterator& -// ConstIterator::operator++ (int) -// { -// Assert (false, ExcNotImplemented()); -// } - - - - - template - inline - const Accessor & - ConstIterator::operator* () const - { - return *this; - } - - - template - inline - const Accessor * - ConstIterator::operator-> () const - { - return this; - } - - - - template - inline - bool - ConstIterator:: - operator == (const ConstIterator& i) const - { - if (this->matrix != i->matrix) - return false; - - if (this->row_block == i->row_block - && this->col_block == i->col_block - && this->base_iterator == i->base_iterator) - return true; - return false; - } - - - - template - inline - bool - ConstIterator:: - operator != (const ConstIterator& i) const - { - return !(*this == i); - } - - - - template - inline - bool - ConstIterator:: - operator < (const ConstIterator& i) const - { - if (this->row_block < i->row_block) - return true; - if (this->row_block == i->row_block) - { - if (this->base_iterator->row() < i->base_iterator->row()) - return true; - if (this->base_iterator->row() == i->base_iterator->row()) - { - if (this->a_index < i->a_index) - return true; - } - } - return false; - } - - } -} - -//----------------------------------------------------------------------// template -inline -unsigned int -BlockSparseMatrix::n_block_cols () const -{ - return columns; -} - - - -template -inline -unsigned int -BlockSparseMatrix::n_block_rows () const -{ - return rows; -} - - -template -inline -typename BlockSparseMatrix::BlockType & -BlockSparseMatrix::block (const unsigned int row, - const unsigned int column) -{ - Assert (row -inline -const typename BlockSparseMatrix::BlockType & -BlockSparseMatrix::block (const unsigned int row, - const unsigned int column) const -{ - Assert (row -inline -unsigned int -BlockSparseMatrix::m () const -{ - return sparsity_pattern->n_rows(); -} - - - -template -inline -unsigned int -BlockSparseMatrix::n () const -{ - return sparsity_pattern->n_cols(); -} - - - -template -inline -void -BlockSparseMatrix::set (const unsigned int i, - const unsigned int j, - const number value) -{ - const std::pair - row_index = sparsity_pattern->row_indices.global_to_local (i), - col_index = sparsity_pattern->column_indices.global_to_local (j); - block(row_index.first,col_index.first).set (row_index.second, - col_index.second, - value); -} - - - -template -inline -BlockSparseMatrix & -BlockSparseMatrix::operator *= (const number factor) -{ - Assert (columns != 0, ExcMatrixNotInitialized()); - Assert (rows != 0, ExcMatrixNotInitialized()); - - for (unsigned int r=0; r -inline -BlockSparseMatrix & -BlockSparseMatrix::operator /= (const number factor) -{ - Assert (columns != 0, ExcMatrixNotInitialized()); - Assert (rows != 0, ExcMatrixNotInitialized()); - Assert (factor !=0, ExcDivideByZero()); - - const number factor_inv = 1. / factor; - - for (unsigned int r=0; r -inline -void -BlockSparseMatrix::add (const unsigned int i, - const unsigned int j, - const number value) -{ - const std::pair - row_index = sparsity_pattern->row_indices.global_to_local (i), - col_index = sparsity_pattern->column_indices.global_to_local (j); - block(row_index.first,col_index.first).add (row_index.second, - col_index.second, - value); -} - - - -template -inline -number -BlockSparseMatrix::operator () (const unsigned int i, - const unsigned int j) const -{ - const std::pair - row_index = sparsity_pattern->row_indices.global_to_local (i), - col_index = sparsity_pattern->column_indices.global_to_local (j); - return block(row_index.first,col_index.first) (row_index.second, - col_index.second); -} - - - -template -inline -number -BlockSparseMatrix::el (const unsigned int i, - const unsigned int j) const -{ - const std::pair - row_index = sparsity_pattern->row_indices.global_to_local (i), - col_index = sparsity_pattern->column_indices.global_to_local (j); - return block(row_index.first,col_index.first).el (row_index.second, - col_index.second); -} - - - - -template -template -BlockSparseMatrix & -BlockSparseMatrix:: -copy_from (const BlockSparseMatrix &source) -{ - for (unsigned int r=0; r -template +template void BlockSparseMatrix:: -add_scaled (const number factor, - const BlockSparseMatrix &matrix) -{ - for (unsigned int r=0; r -template -void -BlockSparseMatrix::vmult (BlockVector &dst, - const BlockVector &src) const -{ - Assert (dst.n_blocks() == rows, - ExcDimensionMismatch(dst.n_blocks(), rows)); - Assert (src.n_blocks() == columns, - ExcDimensionMismatch(src.n_blocks(), columns)); - - for (unsigned int row=0; row -template -void -BlockSparseMatrix::vmult (Vector &dst, - const BlockVector &src) const -{ - Assert (rows == 1, - ExcDimensionMismatch(1, rows)); - Assert (src.n_blocks() == columns, - ExcDimensionMismatch(src.n_blocks(), columns)); - - block(0,0).vmult (dst, src.block(0)); - for (unsigned int col=1; col -template -void -BlockSparseMatrix::vmult (BlockVector &dst, - const Vector &src) const -{ - Assert (dst.n_blocks() == rows, - ExcDimensionMismatch(dst.n_blocks(), rows)); - Assert (1 == columns, - ExcDimensionMismatch(1, columns)); - - for (unsigned int row=0; row -template -void -BlockSparseMatrix::vmult (Vector &dst, - const Vector &src) const +precondition_Jacobi (BlockVectorType &dst, + const BlockVectorType &src, + const number omega) const { - Assert (1 == rows, - ExcDimensionMismatch(1, rows)); - Assert (1 == columns, - ExcDimensionMismatch(1, columns)); + Assert (this->n_block_rows() == this->n_block_cols(), + typename BaseClass::ExcMatrixNotBlockSquare()); + Assert (dst.n_blocks() == this->n_block_rows(), + ExcDimensionMismatch(dst.n_blocks(), this->n_block_rows())); + Assert (src.n_blocks() == this->n_block_cols(), + ExcDimensionMismatch(src.n_blocks(), this->n_block_cols())); - block(0,0).vmult (dst, src); + // do a diagonal preconditioning. uses only + // the diagonal blocks of the matrix + for (unsigned int i=0; in_block_rows(); ++i) + this->block(i,i).precondition_Jacobi (dst.block(i), + src.block(i), + omega); } - -template -template -void -BlockSparseMatrix::vmult_add (BlockVector &dst, - const BlockVector &src) const -{ - Assert (dst.n_blocks() == rows, - ExcDimensionMismatch(dst.n_blocks(), rows)); - Assert (src.n_blocks() == columns, - ExcDimensionMismatch(src.n_blocks(), columns)); - - for (unsigned int row=0; row -template -void -BlockSparseMatrix::Tvmult (BlockVector &dst, - const BlockVector &src) const -{ - Assert (dst.n_blocks() == n_block_cols(), - ExcDimensionMismatch(dst.n_blocks(), n_block_cols())); - Assert (src.n_blocks() == n_block_rows(), - ExcDimensionMismatch(src.n_blocks(), n_block_rows())); - - dst = 0.; - - for (unsigned int row=0; row -template -void -BlockSparseMatrix::Tvmult (BlockVector &dst, - const Vector &src) const -{ - Assert (dst.n_blocks() == n_block_cols(), - ExcDimensionMismatch(dst.n_blocks(), n_block_cols())); - Assert (1 == n_block_rows(), - ExcDimensionMismatch(1, n_block_rows())); - - dst = 0.; - - for (unsigned int col=0; col -template -void -BlockSparseMatrix::Tvmult (Vector &dst, - const BlockVector &src) const -{ - Assert (1 == n_block_cols(), - ExcDimensionMismatch(1, n_block_cols())); - Assert (src.n_blocks() == n_block_rows(), - ExcDimensionMismatch(src.n_blocks(), n_block_rows())); - - block(0,0).Tvmult (dst, src.block(0)); - - for (unsigned int row=1; row -template -void -BlockSparseMatrix::Tvmult (Vector &dst, - const Vector &src) const -{ - Assert (1 == n_block_cols(), - ExcDimensionMismatch(1, n_block_cols())); - Assert (1 == n_block_rows(), - ExcDimensionMismatch(1, n_block_rows())); - - block(0,0).Tvmult (dst, src); -} - - - -template -template -void -BlockSparseMatrix::Tvmult_add (BlockVector& dst, - const BlockVector& src) const -{ - Assert (dst.n_blocks() == n_block_cols(), - ExcDimensionMismatch(dst.n_blocks(), n_block_cols())); - Assert (src.n_blocks() == n_block_rows(), - ExcDimensionMismatch(src.n_blocks(), n_block_rows())); - - for (unsigned int row=0; row -template -somenumber -BlockSparseMatrix::matrix_norm_square (const BlockVector &v) const -{ - Assert (rows == columns, ExcMatrixNotBlockSquare()); - Assert (v.n_blocks() == rows, - ExcDimensionMismatch(v.n_blocks(), rows)); - - somenumber norm_sqr = 0; - for (unsigned int row=0; row -template -somenumber -BlockSparseMatrix:: -matrix_scalar_product (const BlockVector &u, - const BlockVector &v) const -{ - Assert (u.n_blocks() == rows, - ExcDimensionMismatch(u.n_blocks(), rows)); - Assert (v.n_blocks() == columns, - ExcDimensionMismatch(v.n_blocks(), columns)); - - somenumber result = 0; - for (unsigned int row=0; row -template -somenumber -BlockSparseMatrix:: -residual (BlockVector &dst, - const BlockVector &x, - const BlockVector &b) const -{ - Assert (dst.n_blocks() == rows, - ExcDimensionMismatch(dst.n_blocks(), rows)); - Assert (b.n_blocks() == rows, - ExcDimensionMismatch(b.n_blocks(), rows)); - Assert (x.n_blocks() == columns, - ExcDimensionMismatch(x.n_blocks(), columns)); - // in block notation, the residual is - // r_i = b_i - \sum_j A_ij x_j. - // this can be written as - // r_i = b_i - A_i0 x_0 - \sum_{j>0} A_ij x_j. - // - // for the first two terms, we can - // call the residual function of - // A_i0. for the other terms, we - // use vmult_add. however, we want - // to subtract, so in order to - // avoid a temporary vector, we - // perform a sign change of the - // first two term before, and after - // adding up - for (unsigned int row=0; row -template -void -BlockSparseMatrix:: -precondition_Jacobi (BlockVector &dst, - const BlockVector &src, - const number omega) const -{ - Assert (rows == columns, ExcMatrixNotBlockSquare()); - Assert (dst.n_blocks() == rows, - ExcDimensionMismatch(dst.n_blocks(), rows)); - Assert (src.n_blocks() == columns, - ExcDimensionMismatch(src.n_blocks(), columns)); - - for (unsigned int i=0; i -template -void -BlockSparseMatrix:: -precondition_Jacobi (Vector &dst, - const Vector &src, - const number omega) const -{ - Assert (rows == columns, ExcMatrixNotBlockSquare()); - Assert (1 == rows, - ExcDimensionMismatch(1, rows)); - Assert (1 == columns, - ExcDimensionMismatch(1, columns)); - - block(0,0).precondition_Jacobi (dst, src, omega); -} - - - -template -inline -typename BlockSparseMatrix::const_iterator -BlockSparseMatrix::begin () const -{ - return const_iterator(this, 0, 0); -} - - - -template -inline -typename BlockSparseMatrix::const_iterator -BlockSparseMatrix::end () const -{ - return const_iterator(this, m(), 0); -} - - - -template -inline -typename BlockSparseMatrix::const_iterator -BlockSparseMatrix::begin (const unsigned int r) const -{ - Assert (r -inline -typename BlockSparseMatrix::const_iterator -BlockSparseMatrix::end (const unsigned int r) const -{ - Assert (r BlockSparseMatrix::BlockSparseMatrix () : - rows (0), - columns (0), - sparsity_pattern (0) + sparsity_pattern (0) {} @@ -34,8 +32,7 @@ template BlockSparseMatrix:: BlockSparseMatrix (const BlockSparsityPattern &sparsity) : - rows (0), - columns (0) + sparsity_pattern (0) { reinit (sparsity); } @@ -47,13 +44,15 @@ BlockSparseMatrix::~BlockSparseMatrix () { // delete previous content of // the subobjects array - for (unsigned int r=0; rn_block_rows(); ++r) + for (unsigned int c=0; cn_block_cols(); ++c) { - BlockType *p = sub_objects[r][c]; - sub_objects[r][c] = 0; + BlockType *p = this->sub_objects[r][c]; + this->sub_objects[r][c] = 0; delete p; } + + sparsity_pattern = 0; } @@ -63,15 +62,18 @@ BlockSparseMatrix & BlockSparseMatrix:: operator = (const BlockSparseMatrix &m) { - Assert (rows == m.rows, ExcDimensionMismatch(rows, m.rows)); - Assert (columns == m.columns, ExcDimensionMismatch(columns, m.columns)); - // this operator does not do + Assert (this->row_block_indices == m.row_block_indices, + ExcBlockDimensionMismatch()); + Assert (this->column_block_indices == m.column_block_indices, + ExcBlockDimensionMismatch()); + + // this operator does not do // anything except than checking // whether the base objects want to // do something - for (unsigned int r=0; rn_block_rows(); ++r) + for (unsigned int c=0; cn_block_cols(); ++c) + this->block(r,c) = m.block(r,c); return *this; } @@ -84,9 +86,9 @@ BlockSparseMatrix::operator = (const double d) { Assert (d==0, ExcScalarAssignmentOnlyForZeroValue()); - for (unsigned int r=0; rn_block_rows(); ++r) + for (unsigned int c=0; cn_block_cols(); ++c) + block(r,c) = d; return *this; } @@ -101,11 +103,11 @@ reinit (const BlockSparsityPattern &sparsity) // first delete previous content of // the subobjects array and delete // the table completely - for (unsigned int r=0; rn_block_rows(); ++r) + for (unsigned int c=0; cn_block_cols(); ++c) { - BlockType *p = sub_objects[r][c]; - sub_objects[r][c] = 0; + BlockType *p = this->sub_objects[r][c]; + this->sub_objects[r][c] = 0; delete p; } sub_objects.reinit (0,0); @@ -113,41 +115,32 @@ reinit (const BlockSparsityPattern &sparsity) // then associate new sparsity // pattern and resize sparsity_pattern = &sparsity; - rows = sparsity.n_block_rows(); - columns = sparsity.n_block_cols(); - sub_objects.reinit (rows, columns); + + this->row_block_indices = sparsity.row_indices; + this->column_block_indices = sparsity.column_indices; + + this->sub_objects.reinit (sparsity.n_block_rows(), + sparsity.n_block_cols()); // and reinitialize the blocks - for (unsigned int r=0; rn_block_rows(); ++r) + for (unsigned int c=0; cn_block_cols(); ++c) { BlockType *p = new SparseMatrix(); p->reinit (sparsity.block(r,c)); - sub_objects[r][c] = p; + this->sub_objects[r][c] = p; } } -template -void -BlockSparseMatrix::clear () -{ - sparsity_pattern = 0; - for (unsigned int r=0; r bool BlockSparseMatrix::empty () const { - for (unsigned int r=0; rn_block_rows(); ++r) + for (unsigned int c=0; cn_block_cols(); ++c) + if (this->block(r,c).empty () == false) return false; return true; @@ -170,9 +163,9 @@ unsigned int BlockSparseMatrix::n_actually_nonzero_elements () const { unsigned int count = 0; - for (unsigned int i=0; in_actually_nonzero_elements (); + for (unsigned int i=0; in_block_rows(); ++i) + for (unsigned int j=0; jn_block_cols(); ++j) + count += this->sub_objects[i][j]->n_actually_nonzero_elements (); return count; } @@ -198,12 +191,12 @@ print_formatted (std::ostream &out, const char *zero_string, const double denominator) const { - for (unsigned int r=0;rn_block_rows();++r) + for (unsigned int c=0;cn_block_cols();++c) { out << "Component (" << r << "," << c << ")" << std::endl; - block(r,c).print_formatted (out, precision, scientific, - width, zero_string, denominator); + this->block(r,c).print_formatted (out, precision, scientific, + width, zero_string, denominator); } } @@ -214,10 +207,10 @@ unsigned int BlockSparseMatrix::memory_consumption () const { unsigned int mem = sizeof(*this); - mem += MemoryConsumption::memory_consumption (sub_objects); - for (unsigned int r=0; rsub_objects); + for (unsigned int r=0; rn_block_rows(); ++r) + for (unsigned int c=0; cn_block_cols(); ++c) + mem += MemoryConsumption::memory_consumption(*this->sub_objects[r][c]); return mem; } diff --git a/deal.II/lac/include/lac/block_vector.h b/deal.II/lac/include/lac/block_vector.h index ca4f4ce557..6c38c645e6 100644 --- a/deal.II/lac/include/lac/block_vector.h +++ b/deal.II/lac/include/lac/block_vector.h @@ -15,647 +15,66 @@ #include +#include #include #include -#include +#include + #include #include -#include -// forward declaration -template -class BlockVector; /*! @addtogroup Vectors *@{ */ -namespace internal -{ /** - * Namespace in which iterators in block vectors are implemented. - * - * @author Wolfgang Bangerth, 2001 - */ - namespace BlockVectorIterators - { - /** - * Declaration of the general - * template of a structure which is - * used to determine some types - * based on the template arguments - * of other classes. - */ - template - struct Types - { - }; - - - - /** - * Declaration of a specialized - * template of a structure which is - * used to determine some types - * based on the template arguments - * of other classes. - * - * This is for the use of non-const - * iterators. - */ - template - struct Types - { - /** - * Type of the vector - * underlying the block vector - * used in non-const - * iterators. There, the - * vector must not be constant. - */ - typedef ::Vector Vector; - - /** - * Type of the block vector - * used in non-const - * iterators. There, the block - * vector must not be constant. - */ - typedef ::BlockVector BlockVector; - - /** - * Type of the numbers we point - * to. Here, they are not - * constant. - */ - typedef number NumberType; - }; - - - - /** - * Declaration of a specialized - * template of a structure which is - * used to determine some types - * based on the template arguments - * of other classes. - * - * This is for the use of - * const_iterator. - */ - template - struct Types - { - /** - * Type of the vector - * underlying the block vector - * used in - * const_iterator. There, - * the vector must be - * constant. - */ - typedef const ::Vector Vector; - - /** - * Type of the block vector - * used in - * const_iterator. There, - * the block vector must be - * constant. - */ - typedef const ::BlockVector BlockVector; - - /** - * Type of the numbers we point - * to. Here, they are constant - * since the block vector we - * use is constant. - */ - typedef const number NumberType; - }; - - - /** - * General random-access iterator - * class for block vectors. Since - * we do not want to have two - * classes for non-const - * iterator and - * const_iterator, we take a - * second template argument which - * denotes whether the vector we - * point into is a constant object - * or not. The first template - * argument is always the number - * type of the block vector in use. - * - * This class satisfies all - * requirements of random access - * iterators defined in the C++ - * standard. Operations on these - * iterators are constant in the - * number of elements in the block - * vector. However, they are - * sometimes linear in the number - * of blocks in the vector, but - * since that does rarely change - * dynamically within an - * application, this is a constant - * and we again have that the - * iterator satisfies the - * requirements of a random access - * iterator. - * - * The implementation of this class - * has to work around some problems - * in compilers and standard - * libraries. One of these requires - * us to write all comparison - * operators twice, once comparison - * with iterators of the same type - * and once with iterators pointing - * to numbers of opposite constness - * specification. The reason is - * that if we would have written - * the comparison operators as a - * template on the constness of the - * right hand side, then gcc2.95 - * signals an error that these - * operators ambiguate operators - * declared somewhere within the - * standard library. Likewise, we - * have to work around some - * problems with granting other - * iterators friendship. This makes - * the implementation somewhat - * non-optimal at places, but at - * least everything works. - * - * @author Wolfgang Bangerth, 2001 - */ - template - class Iterator : -#ifdef HAVE_STD_ITERATOR_CLASS - public std::iterator::NumberType> -#else - random_access_iterator::NumberType,int> -#endif - { - private: - /** - * Typedef an iterator with - * opposite constness - * requirements on the elements - * it points to. - */ - typedef Iterator InverseConstnessIterator; - - public: - /** - * Declare some typedefs which - * are standard for iterators - * and are used by algorithms - * to enquire about the - * specifics of the iterators - * they work on. - */ - typedef std::random_access_iterator_tag iterator_type; - typedef typename Types::NumberType value_type; - typedef ptrdiff_t difference_type; - typedef value_type &reference; - typedef value_type *pointer; - - /** - * Typedef the type of the - * block vector (which differs - * in constness, depending on - * the second template - * parameter). - */ - typedef typename Types::BlockVector BlockVectorType; - - /** - * Type of the number this - * iterator points - * to. Depending on the value - * of the second template - * parameter, this is either a - * constant or non-const - * number. - */ - typedef typename Types::NumberType NumberType; - - /** - * Construct an iterator from - * a vector to which we point - * and the global index of - * the element pointed to. - * - * Depending on the value of - * the constness template - * argument of this class, - * the first argument of this - * constructor is either is a - * const or non-const - * reference. - */ - Iterator (BlockVectorType &parent, - const unsigned int global_index); - - /** - * Copy constructor. - */ - Iterator (const Iterator &c); - - /** - * Copy constructor for - * conversion between iterators - * with different constness - * requirements. This - * constructor throws an error - * if an attempt is made at - * converting a constant to a - * non-constant iterator. - */ - Iterator (const InverseConstnessIterator &c); - - private: - /** - * Constructor used internally - * in this class. The arguments - * match exactly the values of - * the respective member - * variables. - */ - Iterator (BlockVectorType &parent, - const unsigned int global_index, - const unsigned int current_block, - const unsigned int index_within_block, - const unsigned int next_break_forward, - const unsigned int next_break_backward); - - public: - - /** - * Copy operator. - */ - Iterator & operator = (const Iterator &c); - - /** - * Dereferencing operator. If the - * template argument - * constness is - * true, then no writing to - * the result is possible, making - * this a const_iterator. - */ - reference operator * () const; - - /** - * Dereferencing operator. If the - * template argument - * constness is - * true, then no writing to - * the result is possible, making - * this a const_iterator. - */ - pointer operator -> () const; - - /** - * Random access operator, - * grant access to arbitrary - * elements relative to the one - * presently pointed to. - */ - reference operator [] (const difference_type d) const; - - /** - * Prefix increment operator. This - * operator advances the iterator to - * the next element and returns a - * reference to *this. - */ - Iterator & operator ++ (); - - /** - * Postfix increment - * operator. This operator - * advances the iterator to - * the next element and - * returns a copy of the old - * value of this iterator. - */ - Iterator operator ++ (int); - - /** - * Prefix decrement operator. This - * operator retracts the iterator to - * the previous element and returns a - * reference to *this. - */ - Iterator & operator -- (); - - /** - * Postfix decrement - * operator. This operator - * retracts the iterator to - * the previous element and - * returns a copy of the old - * value of this iterator. - */ - Iterator operator -- (int); - - /** - * Compare for equality of - * iterators. This operator - * checks whether the vectors - * pointed to are the same, - * and if not it throws an - * exception. - */ - bool operator == (const Iterator &i) const; - - /** - * Same, but compare with an - * iterator of different - * constness. - */ - bool operator == (const InverseConstnessIterator &i) const; - - /** - * Compare for inequality of - * iterators. This operator - * checks whether the vectors - * pointed to are the same, - * and if not it throws an - * exception. - */ - bool operator != (const Iterator &i) const; - - /** - * Same, but compare with an - * iterator of different - * constness. - */ - bool operator != (const InverseConstnessIterator &i) const; - - /** - * Check whether this - * iterators points to an - * element previous to the - * one pointed to by the - * given argument. This - * operator checks whether - * the vectors pointed to are - * the same, and if not it - * throws an exception. - */ - bool operator < (const Iterator &i) const; - - /** - * Same, but compare with an - * iterator of different - * constness. - */ - bool operator < (const InverseConstnessIterator &i) const; - - /** - * Comparison operator alike - * to the one above. - */ - bool operator <= (const Iterator &i) const; - - /** - * Same, but compare with an - * iterator of different - * constness. - */ - bool operator <= (const InverseConstnessIterator &i) const; - - /** - * Comparison operator alike - * to the one above. - */ - bool operator > (const Iterator &i) const; - - /** - * Same, but compare with an - * iterator of different - * constness. - */ - bool operator > (const InverseConstnessIterator &i) const; - - /** - * Comparison operator alike - * to the one above. - */ - bool operator >= (const Iterator &i) const; - - /** - * Same, but compare with an - * iterator of different - * constness. - */ - bool operator >= (const InverseConstnessIterator &i) const; - - /** - * Return the distance between - * the two iterators, in - * elements. - */ - difference_type operator - (const Iterator &i) const; - - /** - * Same, but for iterators of - * opposite constness. - */ - difference_type operator - (const InverseConstnessIterator &i) const; - - /** - * Return an iterator which is - * the given number of elements - * in front of the present one. - */ - Iterator operator + (const difference_type &d) const; - - /** - * Return an iterator which is - * the given number of elements - * behind the present one. - */ - Iterator operator - (const difference_type &d) const; - - /** - * Move the iterator d - * elements forward at once, - * and return the result. - */ - Iterator & operator += (const difference_type &d); - - /** - * Move the iterator d - * elements backward at once, - * and return the result. - */ - Iterator & operator -= (const difference_type &d); - - /** - * Exception. - */ - DeclException0 (ExcPointerToDifferentVectors); - /** - * Exception. - */ - DeclException0 (ExcCastingAwayConstness); - - private: - /** - * Pointer to the block - * vector object to which - * this iterator - * points. Depending on the - * value of the constness - * template argument of this - * class, this is a const - * or non-const pointer. - */ - BlockVectorType *parent; - - /** - * Global index of the - * element to which we - * presently point. - */ - unsigned int global_index; - - /** - * Current block and index - * within this block of the - * element presently pointed - * to. - */ - unsigned int current_block; - unsigned int index_within_block; - - /** - * Indices of the global - * element address at which - * we have to move on to - * another block when moving - * forward and - * backward. These indices - * are kept as a cache since - * this is much more - * efficient than always - * asking the parent object. - */ - unsigned int next_break_forward; - unsigned int next_break_backward; - - /** - * Move forward one element. - */ - void move_forward (); - - /** - * Move backward one element. - */ - void move_backward (); - - -#ifndef DEAL_II_NAMESP_TEMPL_FRIEND_BUG - /** - * Mark all other instances of - * this template as friends. In - * fact, we only need the - * inverse constness iterator - * as friend, but this is - * something that ISO C++ does - * not allow to specify. If we - * have detected a compiler bug - * during configuration of the - * library, use a workaround - * that works for this - * particular compiler, but is - * not ISO C++ conforming. - */ - template - friend class Iterator; -#else - friend class InverseConstnessIterator; -#endif - }; - } // namespace BlockVectorIterators -} // namespace internal - - - -/** - * A vector composed of several blocks each representing a vector of - * its own. - * - * The BlockVector is a collection of normal LAC-Vectors. Each of - * the vectors inside can have a different size. The special case of a - * block vector with constant block size is supported by constructor - * and reinit() functions. - * - * The functionality of BlockVector includes everything a Vector can do, plus - * the access to a single Vector inside the BlockVector by - * block(i). It also has a complete random access iterator, just as - * the LAC Vector class or the standard C++ library template - * std::vector. Therefore, all algorithms working on iterators also - * work with objects of this class. - * - * - * @section BlockVectorAccess Accessing individual blocks, and resizing vectors - * - * Apart from using this object as a whole, you can use each block - * separately as a Vector, using the block() function. There - * is a single caveat: if you have changed the size of one of several - * block, you must call the function collect_sizes() of the block - * vector to update its internal structures. - * - * @attention Warning: If you change the sizes of single blocks - * without calling collect_sizes(), results may be unpredictable. The - * debug version does not check consistency here for performance - * reasons! + * An implementation of block vectors based on deal.II vectors. While the base + * class provides for most of the interface, this class handles the actual + * allocation of vectors and provides functions that are specific to the + * underlying vector type. * * @ref Instantiations: some (@ @) * - * @author Wolfgang Bangerth, Guido Kanschat, 1999, 2000, 2001, 2002 + * @author Wolfgang Bangerth, Guido Kanschat, 1999, 2000, 2001, 2002, 2004 */ template -class BlockVector +class BlockVector : public BlockVectorBase > { public: - /* - * Declare standard types used in - * all containers. These types - * parallel those in the - * C++ standard - * libraries - * std::vector<...> - * class. This includes iterator - * types. + /** + * Typedef the base class for simpler + * access to its own typedefs. + */ + typedef BlockVectorBase > BaseClass; + + /** + * Typedef the type of the underlying + * vector. + */ + typedef typename BaseClass::BlockType BlockType; + + /** + * Import the typedefs from the base + * class. */ - typedef Number value_type; - typedef value_type *pointer; - typedef const value_type *const_pointer; - typedef internal::BlockVectorIterators::Iterator iterator; - typedef internal::BlockVectorIterators::Iterator const_iterator; - typedef value_type &reference; - typedef const value_type &const_reference; - typedef std::size_t size_type; + typedef typename BaseClass::value_type value_type; + typedef typename BaseClass::pointer pointer; + typedef typename BaseClass::const_pointer const_pointer; + typedef typename BaseClass::reference reference; + typedef typename BaseClass::const_reference const_reference; + typedef typename BaseClass::size_type size_type; + typedef typename BaseClass::iterator iterator; + typedef typename BaseClass::const_iterator const_iterator; /** * Constructor. There are three * ways to use this * constructor. First, without * any arguments, it generates - * an objetct with no + * an object with no * blocks. Given one argument, * it initializes num_blocks * blocks, but these blocks have @@ -669,8 +88,8 @@ class BlockVector * use blocks of different * sizes. */ - explicit BlockVector (unsigned int num_blocks = 0, - unsigned int block_size = 0); + explicit BlockVector (const unsigned int num_blocks = 0, + const unsigned int block_size = 0); /** * Copy-Constructor. Dimension set to @@ -680,30 +99,41 @@ class BlockVector BlockVector (const BlockVector& V); -// note: I disabled this function for the time being, since egcs1.1.2 -// does not respect the "explicit" keyword for template constructors. -// this leads to unwanted conversions and in some places to automatically -// generated temporaries, where this is not a good idea. [WB] -// /** -// * Copy constructor taking a BlockVector of -// * another data type. This will fail if -// * there is no conversion path from -// * OtherNumber to Number. Note that -// * you may lose accuracy when copying -// * to a BlockVector with data elements with -// * less accuracy. -// */ -// template -// explicit -// BlockVector (const BlockVector &v); - +#ifndef DEAL_II_EXPLICIT_CONSTRUCTOR_BUG /** - * Constructor. Set the number of - * blocks to n.size() and - * initialize each block with - * n[i] zero elements. - */ - BlockVector (const std::vector &n); + * Copy constructor taking a BlockVector of + * another data type. This will fail if + * there is no conversion path from + * OtherNumber to Number. Note that + * you may lose accuracy when copying + * to a BlockVector with data elements with + * less accuracy. + * + * Older versions of gcc did not honor + * the @p explicit keyword on template + * constructors. In such cases, it is + * easy to accidentally write code that + * can be very inefficient, since the + * compiler starts performing hidden + * conversions. To avoid this, this + * function is disabled if we have + * detected a broken compiler during + * configuration. + */ + template + explicit + BlockVector (const BlockVector &v); +#endif + + /** + * Constructor. Set the number of + * blocks to + * block_sizes.size() and + * initialize each block with + * block_sizes[i] zero + * elements. + */ + BlockVector (const std::vector &block_sizes); /** * Constructor. Set the number of @@ -734,6 +164,28 @@ class BlockVector ~BlockVector (); /** + * Copy operator: fill all components of + * the vector with the given scalar + * value. + */ + BlockVector & operator = (const value_type s); + + /** + * Copy operator for arguments of the + * same type. + */ + BlockVector & + operator= (const BlockVector &V); + + /** + * Copy operator for template arguments + * of different types. + */ + template + BlockVector & + operator= (const BlockVector &V); + + /** * Reinitialize the BlockVector to * contain num_blocks blocks of * size block_size each. @@ -746,22 +198,17 @@ class BlockVector const bool fast = false); /** - * Reinitialize the BlockVector - * such that it contains - * N.size() blocks. Each - * Block is reinitialized to - * dimension N[i]. + * Reinitialize the BlockVector such that + * it contains + * block_sizes.size() + * blocks. Each block is reinitialized to + * dimension block_sizes[i]. * * If the number of blocks is the * same as before this function * was called, all vectors remain * the same and reinit() is - * called for each vector. While - * reinitailizing a usual vector - * can consume a lot of time, - * this function here definitely - * has a potential to slow down a - * program considerably. + * called for each vector. * * If fast==false, the vector * is filled with zeros. @@ -774,9 +221,9 @@ class BlockVector * block, to allow the block * vector to update its caches of * vector sizes. If you call - * reinit() of one of the + * reinit() on one of the * blocks, then subsequent - * actions of this object may + * actions on this object may * yield unpredictable results * since they may be routed to * the wrong block. @@ -813,19 +260,29 @@ class BlockVector template void reinit (const BlockVector &V, const bool fast=false); + + /** + * Scale each element of the + * vector by the given factor. + * + * This function is deprecated + * and will be removed in a + * future version. Use + * operator *= and + * operator /= instead. + * + * @deprecated Use operator*= + * instead. + */ + void scale (const value_type factor); /** - * Update internal structures - * after resizing - * vectors. Whenever you reinited - * a block of a block vector, the - * internal data structures are - * corrupted. Therefore, you - * should call this function - * after al blocks got their new - * size. + * Multiply each element of this + * vector by the corresponding + * element of v. */ - void collect_sizes (); + template + void scale (const BlockVector2 &v); /** * Swap the contents of this @@ -857,307 +314,8 @@ class BlockVector * u.swap(v), again in analogy * to standard functions. */ - void swap (BlockVector &v); - - /** - * Access to a single block. - */ - Vector & block (const unsigned int i); - - /** - * Read-only access to a single block. - */ - const Vector & - block (const unsigned int i) const; - - /** - * Return a reference on the - * object that describes the - * mapping between block and - * global indices. The use of - * this function is highly - * deprecated and it should - * vanish in one of the next - * versions - */ - const BlockIndices & - get_block_indices () const; - - /** - * Set all components of the vector to - * the given number @p s. Simply pass - * this down to the individual block - * objects, but we still need to declare - * this function to make the example - * given in the discussion about making - * the constructor explicit work. - * - * - * Since the semantics of assigning a - * scalar to a vector are not - * immediately clear, this operator - * should really only be used if you - * want to set the entire vector to - * zero. This allows the intuitive - * notation v=0. Assigning - * other values is deprecated and may - * be disallowed in the future. - */ - BlockVector & operator = (const Number s); - - /** - * $U = V$: copy all components. - */ - BlockVector & - operator= (const BlockVector& V); - - /** - * $U = V$ for different types. - */ - template - BlockVector & - operator= (const BlockVector& V); - - /** - * $U = U * V$: scalar product. - */ - Number operator* (const BlockVector& V) const; - - /** - * Return square of the $l_2$-norm. - */ - Number norm_sqr () const; - - /** - * Return the mean value of the elements - * of this vector. - */ - Number mean_value () const; - - /** - * Return the $l_1$-norm of the vector, - * i.e. the sum of the absolute values. - */ - Number l1_norm () const; - - /** - * Return the $l_2$-norm of the vector, - * i.e. the square root of the sum of - * the squares of the elements. - */ - Number l2_norm () const; - - /** - * Return the maximum absolute value of - * the elements of this vector, which is - * the $l_\infty$-norm of a vector. - */ - Number linfty_norm () const; - - /** - * Number of blocks. - */ - unsigned int n_blocks () const; - - /** - * Return dimension of the vector. This - * is the sum of the dimensions of all - * components. - */ - unsigned int size () const; - - /** - * Return whether the vector contains only - * elements with value zero. This function - * is mainly for internal consistency - * check and should seldomly be used when - * not in debug mode since it uses quite - * some time. - */ - bool all_zero () const; - - /** - * Return @p true if the vector has no - * negative entries, i.e. all entries are - * zero or positive. This function is - * used, for example, to check whether - * refinement indicators are really all - * positive (or zero). - */ - bool is_non_negative () const; - - /** - * @name 2: Data-Access - */ - //@{ - /** - * Access components, returns U(i). - */ - Number operator() (const unsigned int i) const; - - /** - * Access components, returns U(i) - * as a writeable reference. - */ - Number& operator() (const unsigned int i); - - /** - * Return an iterator pointing to - * the first element. - */ - iterator begin (); - - /** - * Return an iterator pointing to - * the first element of a - * constant block vector. - */ - const_iterator begin () const; - - /** - * Return an iterator pointing to - * the element past the end. - */ - iterator end (); - - /** - * Return an iterator pointing to - * the element past the end of a - * constant block vector. - */ - const_iterator end () const; - //@} - - - /** - * @name 3: Modification of vectors - */ - //@{ - /** - * Addition operator. Fast equivalent to - * U.add(1, V). - */ - BlockVector & - operator += (const BlockVector &V); - - /** - * Subtraction operator. Fast equivalent - * to U.add(-1, V). - */ - BlockVector & - operator -= (const BlockVector &V); + void swap (BlockVector &v); - /** - * $U(0-DIM)+=s$. Addition of s - * to all components. Note that - * s is a scalar and not a - * vector. - */ - void add (const Number s); - - /** - * U+=V. - * Simple vector addition, equal to the - * operator +=. - */ - void add (const BlockVector& V); - - /** - * U+=a*V. - * Simple addition of a scaled vector. - */ - void add (const Number a, const BlockVector& V); - - /** - * U+=a*V+b*W. - * Multiple addition of scaled vectors. - */ - void add (const Number a, const BlockVector& V, - const Number b, const BlockVector& W); - - /** - * U=s*U+V. - * Scaling and simple vector addition. - */ - void sadd (const Number s, const BlockVector& V); - - /** - * U=s*U+a*V. - * Scaling and simple addition. - */ - void sadd (const Number s, const Number a, const BlockVector& V); - - /** - * U=s*U+a*V+b*W. - * Scaling and multiple addition. - */ - void sadd (const Number s, const Number a, - const BlockVector& V, - const Number b, const BlockVector& W); - - /** - * U=s*U+a*V+b*W+c*X. - * Scaling and multiple addition. - */ - void sadd (const Number s, const Number a, - const BlockVector& V, - const Number b, const BlockVector& W, - const Number c, const BlockVector& X); - - /** - * Scale each element of the - * vector by the given factor. - * - * This function is deprecated - * and will be removed in a - * future version. Use - * operator *= and - * operator /= instead. - */ - void scale (const Number factor); - - /** - * Scale each element of the - * vector by a constant - * value. - */ - BlockVector & operator *= (const Number factor); - - /** - * Scale each element of the - * vector by the inverse of the - * given value. - */ - BlockVector & operator /= (const Number factor); - - /** - * Multiply each element of this - * vector by the corresponding - * element of v. - */ - template - void scale (const BlockVector& v); - - /** - * U=a*V. Assignment. - */ - template - void equ (const Number a, const BlockVector& V); - - /** - * U=a*V+b*W. - * Replacing by sum. - */ - void equ (const Number a, const BlockVector& V, - const Number b, const BlockVector& W); - - //@} - - - /** - * @name 5: Mixed stuff - */ - //@{ /** * Output of vector in user-defined * format. @@ -1206,50 +364,10 @@ class BlockVector */ unsigned int memory_consumption () const; - //@} - /** * Exception */ DeclException0 (ExcIteratorRangeDoesNotMatchVectorSize); - - protected: - /** - * Pointer to the array of components. - */ - std::vector > components; - - /** - * Object managing the - * transformation between global - * indices and indices within the - * different blocks. - */ - BlockIndices block_indices; - - private: - /** - * The number of blocks. This - * number is redundant to - * components.size() and stored - * here for convenience. - */ - unsigned int num_blocks; - - /** - * Make the iterator class a - * friend. We have to work around - * a compiler bug here again. - */ -#ifndef DEAL_II_NAMESP_TEMPL_FRIEND_BUG - template - friend class internal::BlockVectorIterators::Iterator; -#else - friend class iterator; - friend class const_iterator; -#endif - - template friend class BlockVector; }; /*@}*/ @@ -1280,42 +398,13 @@ BlockVector::BlockVector (const std::vector &n, - -template -inline -unsigned int BlockVector::size () const -{ - return block_indices.total_size(); -} - - -template -inline -unsigned int BlockVector::n_blocks () const -{ - return num_blocks; -} - - - -template -inline -Number BlockVector::operator() (const unsigned int i) const -{ - const std::pair local_index - = block_indices.global_to_local (i); - return components[local_index.first](local_index.second); -} - - - template inline -Number& BlockVector::operator() (const unsigned int i) +BlockVector & +BlockVector::operator = (const value_type s) { - const std::pair local_index - = block_indices.global_to_local (i); - return components[local_index.first](local_index.second); + BaseClass::operator = (s); + return *this; } @@ -1323,540 +412,43 @@ Number& BlockVector::operator() (const unsigned int i) template inline BlockVector & -BlockVector::operator *= (const Number factor) +BlockVector::operator = (const BlockVector &v) { - scale (factor); + BaseClass::operator = (v); return *this; } template +template inline BlockVector & -BlockVector::operator /= (const Number factor) +BlockVector::operator = (const BlockVector &v) { - scale (1./factor); + BaseClass::operator = (v); return *this; } - template -inline -Vector & -BlockVector::block (const unsigned int i) +void BlockVector::scale (const value_type factor) { - Assert(i -inline -const Vector & -BlockVector::block (const unsigned int i) const +template +void BlockVector::scale (const BlockVector2 &v) { - Assert(i -inline -const BlockIndices& -BlockVector::get_block_indices () const -{ - return block_indices; -} - - -namespace internal -{ - namespace BlockVectorIterators - { - - template - inline - Iterator:: - Iterator (const Iterator &c) - : - parent (c.parent), - global_index (c.global_index), - current_block (c.current_block), - index_within_block (c.index_within_block), - next_break_forward (c.next_break_forward), - next_break_backward (c.next_break_backward) - {} - - - - template - inline - Iterator:: - Iterator (const InverseConstnessIterator &c) - : - parent (const_cast(c.parent)), - global_index (c.global_index), - current_block (c.current_block), - index_within_block (c.index_within_block), - next_break_forward (c.next_break_forward), - next_break_backward (c.next_break_backward) - { - // if constness==false, then the - // constness of the iterator we - // got is true and we are trying - // to cast away the - // constness. disallow this - Assert (constness==true, ExcCastingAwayConstness()); - } - - - - template - inline - Iterator:: - Iterator (BlockVectorType &parent, - const unsigned int global_index, - const unsigned int current_block, - const unsigned int index_within_block, - const unsigned int next_break_forward, - const unsigned int next_break_backward) - : - parent (&parent), - global_index (global_index), - current_block (current_block), - index_within_block (index_within_block), - next_break_forward (next_break_forward), - next_break_backward (next_break_backward) - { - } - - - - template - inline - Iterator & - Iterator:: - operator = (const Iterator &c) - { - parent = c.parent; - global_index = c.global_index; - index_within_block = c.index_within_block; - current_block = c.current_block; - next_break_forward = c.next_break_forward; - next_break_backward = c.next_break_backward; - - return *this; - } - - - - template - inline - typename Iterator::reference - Iterator::operator * () const - { - // we might want to return the - // following directly, but if we do - // that we run into a gcc bug where - // it complains that we return a - // reference to a temporary - reference p = parent->block(current_block)(index_within_block); - return p; - } - - - - template - inline - typename Iterator::pointer - Iterator::operator -> () const - { - // we might want to return the - // following directly, but if we do - // that we run into a gcc bug where - // it complains that we return a - // reference to a temporary - reference p = parent->block(current_block)(index_within_block); - return &p; - } - - - - template - inline - typename Iterator::reference - Iterator::operator [] (const difference_type d) const - { - // if the index pointed to is - // still within the block we - // currently point into, then we - // can save the computation of - // the block - if ((global_index+d >= next_break_backward) && - (global_index+d <= next_break_forward)) - { - reference p = parent->block(current_block)(index_within_block + d); - return p; - }; - - - // if the index is not within the - // block of the block vector into - // which we presently point, then - // there is no way: we have to - // search for the block. this can - // be done through the parent - // class as well. - reference p = (*parent)(global_index+d); - return p; - } - - - - template - inline - Iterator & - Iterator::operator ++ () - { - move_forward (); - return *this; - } - - - - template - inline - Iterator - Iterator::operator ++ (int) - { - const Iterator old_value = *this; - move_forward (); - return old_value; - } - - - - template - inline - Iterator & - Iterator::operator -- () - { - move_backward (); - return *this; - } - - - - template - inline - Iterator - Iterator::operator -- (int) - { - const Iterator old_value = *this; - move_backward (); - return old_value; - } - - - - template - inline - bool - Iterator:: - operator == (const Iterator &i) const - { - Assert (parent == i.parent, ExcPointerToDifferentVectors()); - - return (global_index == i.global_index); - } - - - - template - inline - bool - Iterator:: - operator == (const InverseConstnessIterator &i) const - { - Assert (parent == i.parent, ExcPointerToDifferentVectors()); - - return (global_index == i.global_index); - } - - - - template - inline - bool - Iterator:: - operator != (const Iterator &i) const - { - Assert (parent == i.parent, ExcPointerToDifferentVectors()); - - return (global_index != i.global_index); - } - - - - template - inline - bool - Iterator:: - operator != (const InverseConstnessIterator &i) const - { - Assert (parent == i.parent, ExcPointerToDifferentVectors()); - - return (global_index != i.global_index); - } - - - - template - inline - bool - Iterator:: - operator < (const Iterator &i) const - { - Assert (parent == i.parent, ExcPointerToDifferentVectors()); - - return (global_index < i.global_index); - } - - - - template - inline - bool - Iterator:: - operator < (const InverseConstnessIterator &i) const - { - Assert (parent == i.parent, ExcPointerToDifferentVectors()); - - return (global_index < i.global_index); - } - - - - template - inline - bool - Iterator:: - operator <= (const Iterator &i) const - { - Assert (parent == i.parent, ExcPointerToDifferentVectors()); - - return (global_index <= i.global_index); - } - - - - template - inline - bool - Iterator:: - operator <= (const InverseConstnessIterator &i) const - { - Assert (parent == i.parent, ExcPointerToDifferentVectors()); - - return (global_index <= i.global_index); - } - - - - template - inline - bool - Iterator:: - operator > (const Iterator &i) const - { - Assert (parent == i.parent, ExcPointerToDifferentVectors()); - - return (global_index > i.global_index); - } - - - - template - inline - bool - Iterator:: - operator > (const InverseConstnessIterator &i) const - { - Assert (parent == i.parent, ExcPointerToDifferentVectors()); - - return (global_index > i.global_index); - } - - - - template - inline - bool - Iterator:: - operator >= (const Iterator &i) const - { - Assert (parent == i.parent, ExcPointerToDifferentVectors()); - - return (global_index >= i.global_index); - } - - - - template - inline - bool - Iterator:: - operator >= (const InverseConstnessIterator &i) const - { - Assert (parent == i.parent, ExcPointerToDifferentVectors()); - - return (global_index >= i.global_index); - } - - - - template - inline - typename Iterator::difference_type - Iterator:: - operator - (const Iterator &i) const - { - Assert (parent == i.parent, ExcPointerToDifferentVectors()); - - return (static_cast(global_index) - - static_cast(i.global_index)); - } - - - - template - inline - typename Iterator::difference_type - Iterator:: - operator - (const InverseConstnessIterator &i) const - { - Assert (parent == i.parent, ExcPointerToDifferentVectors()); - - return (static_cast(global_index) - - static_cast(i.global_index)); - } - - - - template - inline - Iterator - Iterator:: - operator + (const difference_type &d) const - { - // if the index pointed to is - // still within the block we - // currently point into, then we - // can save the computation of - // the block - if ((global_index+d >= next_break_backward) && - (global_index+d <= next_break_forward)) - return Iterator (*parent, global_index+d, current_block, - index_within_block+d, - next_break_forward, next_break_backward); - else - // outside present block, so - // have to seek new block - // anyway - return Iterator (*parent, global_index+d); - } - - - - template - inline - Iterator - Iterator:: - operator - (const difference_type &d) const - { - // if the index pointed to is - // still within the block we - // currently point into, then we - // can save the computation of - // the block - if ((global_index-d >= next_break_backward) && - (global_index-d <= next_break_forward)) - return Iterator (*parent, global_index-d, current_block, - index_within_block-d, - next_break_forward, next_break_backward); - else - // outside present block, so - // have to seek new block - // anyway - return Iterator (*parent, global_index-d); - } - - - - template - inline - Iterator & - Iterator:: - operator += (const difference_type &d) - { - // if the index pointed to is - // still within the block we - // currently point into, then we - // can save the computation of - // the block - if ((global_index+d >= next_break_backward) && - (global_index+d <= next_break_forward)) - { - global_index += d; - index_within_block += d; - } - else - // outside present block, so - // have to seek new block - // anyway - *this = Iterator (*parent, global_index+d); - - return *this; - } - - - - template - inline - Iterator & - Iterator:: - operator -= (const difference_type &d) - { - // if the index pointed to is - // still within the block we - // currently point into, then we - // can save the computation of - // the block - if ((global_index-d >= next_break_backward) && - (global_index-d <= next_break_forward)) - { - global_index -= d; - index_within_block -= d; - } - else - // outside present block, so - // have to seek new block - // anyway - *this = Iterator (*parent, global_index-d); - - return *this; - } - - } // namespace BlockVectorIterators - -} //namespace internal - /*! @addtogroup Vectors *@{ diff --git a/deal.II/lac/include/lac/block_vector.templates.h b/deal.II/lac/include/lac/block_vector.templates.h index ded153a7e1..5013e26b4b 100644 --- a/deal.II/lac/include/lac/block_vector.templates.h +++ b/deal.II/lac/include/lac/block_vector.templates.h @@ -22,8 +22,8 @@ template -BlockVector::BlockVector (unsigned int n_blocks, - unsigned int block_size) +BlockVector::BlockVector (const unsigned int n_blocks, + const unsigned int block_size) { reinit (n_blocks, block_size); } @@ -39,32 +39,29 @@ BlockVector::BlockVector (const std::vector &n) template BlockVector::BlockVector (const BlockVector& v) - : components(v.num_blocks), - block_indices(v.block_indices), - num_blocks(v.num_blocks) + : + BlockVectorBase > () { + components.resize (v.num_blocks); + block_indices = v.block_indices; + num_blocks = v.num_blocks; + for (unsigned int i=0; i -// template -// BlockVector::Vector (const BlockVector< OtherNumber>& v) : -// dim(v.size()), -// maxdim(v.size()), -// val(0) -// { -// if (dim) -// { -// val = new Number[maxdim]; -// Assert (val != 0, ExcOutOfMemory()); -// copy (v.begin(), v.end(), begin()); -// } -// } +#ifndef DEAL_II_EXPLICIT_CONSTRUCTOR_BUG + +template +template +BlockVector::BlockVector (const BlockVector& v) +{ + reinit (v, true); + *this = v; +} +#endif template void BlockVector::reinit (const unsigned int n_bl, @@ -93,28 +90,15 @@ void BlockVector::reinit (const std::vector &n, template template void BlockVector::reinit (const BlockVector& v, - const bool fast) + const bool fast) { - block_indices = v.block_indices; - num_blocks = v.num_blocks; + block_indices = v.get_block_indices(); + num_blocks = v.n_blocks(); if (components.size() != num_blocks) components.resize(num_blocks); for (unsigned int i=0;i -void -BlockVector::collect_sizes () -{ - std::vector sizes (num_blocks); - - for (unsigned int i=0; i::swap (BlockVector &v) -template -bool -BlockVector::all_zero () const -{ - for (unsigned int i=0;i -bool -BlockVector::is_non_negative () const -{ - for (unsigned int i=0;i -Number BlockVector::operator * (const BlockVector& v) const -{ - Assert (num_blocks == v.num_blocks, - ExcDimensionMismatch(num_blocks, v.num_blocks)); - - Number sum = 0.; - for (unsigned int i=0;i -typename BlockVector::iterator -BlockVector::begin() -{ - return iterator(*this, 0U); -} - - - -template -typename BlockVector::const_iterator -BlockVector::begin() const -{ - return const_iterator(*this, 0U); -} - - -template -typename BlockVector::iterator -BlockVector::end() -{ - return iterator(*this, size()); -} - - - -template -typename BlockVector::const_iterator -BlockVector::end() const -{ - return const_iterator(*this, size()); -} - - - -template -Number BlockVector::norm_sqr () const -{ - Number sum = 0.; - for (unsigned int i=0;i -Number BlockVector::mean_value () const -{ - Number sum = 0.; - for (unsigned int i=0;i -Number BlockVector::l1_norm () const -{ - Number sum = 0.; - for (unsigned int i=0;i -Number BlockVector::l2_norm () const -{ - return std::sqrt(norm_sqr()); -} - - - -template -Number BlockVector::linfty_norm () const -{ - Number sum = 0.; - for (unsigned int i=0;i -BlockVector& -BlockVector::operator += (const BlockVector& v) -{ - add (v); - return *this; -} - - - -template -BlockVector& -BlockVector::operator -= (const BlockVector& v) -{ - Assert (num_blocks == v.num_blocks, - ExcDimensionMismatch(num_blocks, v.num_blocks)); - - for (unsigned int i=0;i -void BlockVector::add (const Number a) -{ - for (unsigned int i=0;i -void BlockVector::add (const BlockVector& v) -{ - Assert (num_blocks == v.num_blocks, - ExcDimensionMismatch(num_blocks, v.num_blocks)); - - for (unsigned int i=0;i -void BlockVector::add (const Number a, - const BlockVector& v) -{ - Assert (num_blocks == v.num_blocks, - ExcDimensionMismatch(num_blocks, v.num_blocks)); - - for (unsigned int i=0;i -void BlockVector::add (const Number a, - const BlockVector& v, - const Number b, - const BlockVector& w) -{ - Assert (num_blocks == v.num_blocks, - ExcDimensionMismatch(num_blocks, v.num_blocks)); - Assert (num_blocks == w.num_blocks, - ExcDimensionMismatch(num_blocks, w.num_blocks)); - - - for (unsigned int i=0;i -void BlockVector::sadd (const Number x, - const BlockVector& v) -{ - Assert (num_blocks == v.num_blocks, - ExcDimensionMismatch(num_blocks, v.num_blocks)); - - for (unsigned int i=0;i -void BlockVector::sadd (const Number x, const Number a, - const BlockVector& v) -{ - Assert (num_blocks == v.num_blocks, - ExcDimensionMismatch(num_blocks, v.num_blocks)); - - for (unsigned int i=0;i -void BlockVector::sadd (const Number x, const Number a, - const BlockVector& v, - const Number b, - const BlockVector& w) -{ - Assert (num_blocks == v.num_blocks, - ExcDimensionMismatch(num_blocks, v.num_blocks)); - Assert (num_blocks == w.num_blocks, - ExcDimensionMismatch(num_blocks, w.num_blocks)); - - for (unsigned int i=0;i -void BlockVector::sadd (const Number x, const Number a, - const BlockVector& v, - const Number b, - const BlockVector& w, - const Number c, - const BlockVector& y) -{ - Assert (num_blocks == v.num_blocks, - ExcDimensionMismatch(num_blocks, v.num_blocks)); - Assert (num_blocks == w.num_blocks, - ExcDimensionMismatch(num_blocks, w.num_blocks)); - Assert (num_blocks == y.num_blocks, - ExcDimensionMismatch(num_blocks, y.num_blocks)); - - for (unsigned int i=0;i -void BlockVector::scale (const Number factor) -{ - for (unsigned int i=0;i -template -void BlockVector::scale (const BlockVector& v) -{ - Assert (num_blocks == v.num_blocks, - ExcDimensionMismatch(num_blocks, v.num_blocks)); - for (unsigned int i=0;i -void BlockVector::equ (const Number a, - const BlockVector& v, - const Number b, - const BlockVector& w) -{ - Assert (num_blocks == v.num_blocks, - ExcDimensionMismatch(num_blocks, v.num_blocks)); - Assert (num_blocks == w.num_blocks, - ExcDimensionMismatch(num_blocks, w.num_blocks)); - - for (unsigned int i=0;i -template -void BlockVector::equ (const Number a, - const BlockVector& v) -{ - Assert (num_blocks == v.num_blocks, - ExcDimensionMismatch(num_blocks, v.num_blocks)); - - for (unsigned int i=0;i -BlockVector& BlockVector::operator = (const Number s) -{ - for (unsigned int i=0;i -BlockVector& -BlockVector::operator = (const BlockVector& v) -{ - Assert (num_blocks == v.num_blocks, - ExcDimensionMismatch(num_blocks, v.num_blocks)); - - for (unsigned int i=0;i -template -BlockVector& -BlockVector::operator = (const BlockVector& v) -{ - Assert (num_blocks == v.num_blocks, - ExcDimensionMismatch(num_blocks, v.num_blocks)); - - for (unsigned int i=0;i void BlockVector::print (std::ostream &out, const unsigned int precision, @@ -585,122 +173,4 @@ BlockVector::memory_consumption () const -namespace internal -{ - namespace BlockVectorIterators - { - - template - Iterator:: - Iterator (BlockVectorType &parent, - const unsigned global_index) - : - parent (&parent), - global_index (global_index) - { - // find which block we are - // in. for this, take into - // account that it happens at - // times that people want to - // initialize iterators - // past-the-end - if (global_index < parent.size()) - { - const std::pair - indices = parent.block_indices.global_to_local(global_index); - current_block = indices.first; - index_within_block = indices.second; - - next_break_backward - = parent.block_indices.local_to_global (current_block, 0); - next_break_forward - = (parent.block_indices.local_to_global (current_block, 0) - +parent.block_indices.block_size(current_block)-1); - } - else - // past the end. only have one - // value for this - { - this->global_index = parent.size (); - current_block = parent.n_blocks(); - index_within_block = 0; - next_break_backward = global_index; - next_break_forward = deal_II_numbers::invalid_unsigned_int; - }; - } - - - - template - void - Iterator::move_forward () - { - if (global_index != next_break_forward) - ++index_within_block; - else - { - // ok, we traverse a boundary - // between blocks: - index_within_block = 0; - ++current_block; - - // break backwards is now old - // break forward - next_break_backward = next_break_forward+1; - - // compute new break forward - if (current_block < parent->block_indices.size()) - next_break_forward - += parent->block_indices.block_size(current_block); - else - // if we are beyond the end, - // then move the next - // boundary arbitrarily far - // away - next_break_forward = deal_II_numbers::invalid_unsigned_int; - }; - - ++global_index; - } - - - - template - void - Iterator::move_backward () - { - if (global_index != next_break_backward) - --index_within_block; - else - if (current_block != 0) - { - // ok, we traverse a boundary - // between blocks: - --current_block; - index_within_block = parent->block_indices.block_size(current_block)-1; - - // break forwards is now old - // break backward - next_break_forward = next_break_backward-1; - - // compute new break forward - next_break_backward - -= parent->block_indices.block_size (current_block); - } - else - // current block was 0, we now - // get into unspecified terrain - { - --current_block; - index_within_block = deal_II_numbers::invalid_unsigned_int; - next_break_forward = 0; - next_break_backward = 0; - }; - - --global_index; - } - } -} - - #endif diff --git a/deal.II/lac/include/lac/block_vector_base.h b/deal.II/lac/include/lac/block_vector_base.h new file mode 100644 index 0000000000..6eca3e5d28 --- /dev/null +++ b/deal.II/lac/include/lac/block_vector_base.h @@ -0,0 +1,2150 @@ +//---------------------------- block_vector_base.h --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- block_vector_base.h --------------------------- +#ifndef __deal2__block_vector_base_h +#define __deal2__block_vector_base_h + + +#include +#include +#include + +#include +#include + + + +/*! @addtogroup Vectors + *@{ + */ + +namespace internal +{ + +/** + * Namespace in which iterators in block vectors are implemented. + * + * @author Wolfgang Bangerth, 2001 + */ + namespace BlockVectorIterators + { + /** + * Declaration of the general + * template of a structure which is + * used to determine some types + * based on the template arguments + * of other classes. + */ + template + struct Types + { + }; + + + + /** + * Declaration of a specialized + * template of a structure which is + * used to determine some types + * based on the template arguments + * of other classes. + * + * This is for the use of non-const + * iterators. + */ + template + struct Types + { + /** + * Type of the vector + * underlying the block vector + * used in non-const + * iterators. There, the + * vector must not be constant. + */ + typedef typename BlockVectorType::BlockType Vector; + + /** + * Type of the block vector + * used in non-const + * iterators. There, the block + * vector must not be constant. + */ + typedef BlockVectorType BlockVector; + + /** + * Type of the numbers we point + * to. Here, they are not + * constant. + */ + typedef typename BlockVector::value_type value_type; + + /** + * Typedef the result of a + * dereferencing operation for an + * iterator of the underlying + * iterator. + */ + typedef typename Vector::reference dereference_type; + }; + + + + /** + * Declaration of a specialized + * template of a structure which is + * used to determine some types + * based on the template arguments + * of other classes. + * + * This is for the use of + * const_iterator. + */ + template + struct Types + { + /** + * Type of the vector + * underlying the block vector + * used in + * const_iterator. There, + * the vector must be + * constant. + */ + typedef const typename BlockVectorType::BlockType Vector; + + /** + * Type of the block vector + * used in + * const_iterator. There, + * the block vector must be + * constant. + */ + typedef const BlockVectorType BlockVector; + + /** + * Type of the numbers we point + * to. Here, they are constant + * since the block vector we + * use is constant. + */ + typedef const typename BlockVector::value_type value_type; + + /** + * Typedef the result of a + * dereferencing operation for an + * iterator of the underlying + * iterator. Since this is for + * constant iterators, we can only + * return values, no actual + * references. + */ + typedef value_type dereference_type; + }; + + + /** + * General random-access iterator + * class for block vectors. Since + * we do not want to have two + * classes for non-const + * iterator and + * const_iterator, we take a + * second template argument which + * denotes whether the vector we + * point into is a constant object + * or not. The first template + * argument is always the number + * type of the block vector in use. + * + * This class satisfies all + * requirements of random access + * iterators defined in the C++ + * standard. Operations on these + * iterators are constant in the + * number of elements in the block + * vector. However, they are + * sometimes linear in the number + * of blocks in the vector, but + * since that does rarely change + * dynamically within an + * application, this is a constant + * and we again have that the + * iterator satisfies the + * requirements of a random access + * iterator. + * + * The implementation of this class + * has to work around some problems + * in compilers and standard + * libraries. One of these requires + * us to write all comparison + * operators twice, once comparison + * with iterators of the same type + * and once with iterators pointing + * to numbers of opposite constness + * specification. The reason is + * that if we would have written + * the comparison operators as a + * template on the constness of the + * right hand side, then gcc2.95 + * signals an error that these + * operators ambiguate operators + * declared somewhere within the + * standard library. Likewise, we + * have to work around some + * problems with granting other + * iterators friendship. This makes + * the implementation somewhat + * non-optimal at places, but at + * least everything works. + * + * @author Wolfgang Bangerth, 2001 + */ + template + class Iterator : +#ifdef HAVE_STD_ITERATOR_CLASS + public std::iterator::value_type> +#else + random_access_iterator::value_type,int> +#endif + { + private: + /** + * Typedef an iterator with + * opposite constness + * requirements on the elements + * it points to. + */ + typedef Iterator InverseConstnessIterator; + + public: + /** + * Declare some typedefs which + * are standard for iterators + * and are used by algorithms + * to enquire about the + * specifics of the iterators + * they work on. + */ + typedef std::random_access_iterator_tag iterator_type; + typedef typename Types::value_type value_type; + typedef ptrdiff_t difference_type; + typedef typename BlockVectorType::reference reference; + typedef value_type *pointer; + + typedef + typename Types::dereference_type + dereference_type; + + /** + * Typedef the type of the + * block vector (which differs + * in constness, depending on + * the second template + * parameter). + */ + typedef + typename Types::BlockVector + BlockVector; + + /** + * Type of the number this + * iterator points + * to. Depending on the value + * of the second template + * parameter, this is either a + * constant or non-const + * number. + */ + typedef + typename Types::value_type + value_type; + + /** + * Construct an iterator from + * a vector to which we point + * and the global index of + * the element pointed to. + * + * Depending on the value of + * the constness template + * argument of this class, + * the first argument of this + * constructor is either is a + * const or non-const + * reference. + */ + Iterator (BlockVector &parent, + const unsigned int global_index); + + /** + * Copy constructor. + */ + Iterator (const Iterator &c); + + /** + * Copy constructor for + * conversion between iterators + * with different constness + * requirements. This + * constructor throws an error + * if an attempt is made at + * converting a constant to a + * non-constant iterator. + */ + Iterator (const InverseConstnessIterator &c); + + private: + /** + * Constructor used internally + * in this class. The arguments + * match exactly the values of + * the respective member + * variables. + */ + Iterator (BlockVector &parent, + const unsigned int global_index, + const unsigned int current_block, + const unsigned int index_within_block, + const unsigned int next_break_forward, + const unsigned int next_break_backward); + + public: + + /** + * Copy operator. + */ + Iterator & operator = (const Iterator &c); + + /** + * Dereferencing operator. If the + * template argument + * constness is + * true, then no writing to + * the result is possible, making + * this a const_iterator. + */ + dereference_type operator * () const; + + /** + * Random access operator, + * grant access to arbitrary + * elements relative to the one + * presently pointed to. + */ + dereference_type operator [] (const difference_type d) const; + + /** + * Prefix increment operator. This + * operator advances the iterator to + * the next element and returns a + * reference to *this. + */ + Iterator & operator ++ (); + + /** + * Postfix increment + * operator. This operator + * advances the iterator to + * the next element and + * returns a copy of the old + * value of this iterator. + */ + Iterator operator ++ (int); + + /** + * Prefix decrement operator. This + * operator retracts the iterator to + * the previous element and returns a + * reference to *this. + */ + Iterator & operator -- (); + + /** + * Postfix decrement + * operator. This operator + * retracts the iterator to + * the previous element and + * returns a copy of the old + * value of this iterator. + */ + Iterator operator -- (int); + + /** + * Compare for equality of + * iterators. This operator + * checks whether the vectors + * pointed to are the same, + * and if not it throws an + * exception. + */ + bool operator == (const Iterator &i) const; + + /** + * Same, but compare with an + * iterator of different + * constness. + */ + bool operator == (const InverseConstnessIterator &i) const; + + /** + * Compare for inequality of + * iterators. This operator + * checks whether the vectors + * pointed to are the same, + * and if not it throws an + * exception. + */ + bool operator != (const Iterator &i) const; + + /** + * Same, but compare with an + * iterator of different + * constness. + */ + bool operator != (const InverseConstnessIterator &i) const; + + /** + * Check whether this + * iterators points to an + * element previous to the + * one pointed to by the + * given argument. This + * operator checks whether + * the vectors pointed to are + * the same, and if not it + * throws an exception. + */ + bool operator < (const Iterator &i) const; + + /** + * Same, but compare with an + * iterator of different + * constness. + */ + bool operator < (const InverseConstnessIterator &i) const; + + /** + * Comparison operator alike + * to the one above. + */ + bool operator <= (const Iterator &i) const; + + /** + * Same, but compare with an + * iterator of different + * constness. + */ + bool operator <= (const InverseConstnessIterator &i) const; + + /** + * Comparison operator alike + * to the one above. + */ + bool operator > (const Iterator &i) const; + + /** + * Same, but compare with an + * iterator of different + * constness. + */ + bool operator > (const InverseConstnessIterator &i) const; + + /** + * Comparison operator alike + * to the one above. + */ + bool operator >= (const Iterator &i) const; + + /** + * Same, but compare with an + * iterator of different + * constness. + */ + bool operator >= (const InverseConstnessIterator &i) const; + + /** + * Return the distance between + * the two iterators, in + * elements. + */ + difference_type operator - (const Iterator &i) const; + + /** + * Same, but for iterators of + * opposite constness. + */ + difference_type operator - (const InverseConstnessIterator &i) const; + + /** + * Return an iterator which is + * the given number of elements + * in front of the present one. + */ + Iterator operator + (const difference_type &d) const; + + /** + * Return an iterator which is + * the given number of elements + * behind the present one. + */ + Iterator operator - (const difference_type &d) const; + + /** + * Move the iterator d + * elements forward at once, + * and return the result. + */ + Iterator & operator += (const difference_type &d); + + /** + * Move the iterator d + * elements backward at once, + * and return the result. + */ + Iterator & operator -= (const difference_type &d); + + /** + * Exception. + */ + DeclException0 (ExcPointerToDifferentVectors); + /** + * Exception. + */ + DeclException0 (ExcCastingAwayConstness); + + private: + /** + * Pointer to the block + * vector object to which + * this iterator + * points. Depending on the + * value of the constness + * template argument of this + * class, this is a const + * or non-const pointer. + */ + BlockVector *parent; + + /** + * Global index of the + * element to which we + * presently point. + */ + unsigned int global_index; + + /** + * Current block and index + * within this block of the + * element presently pointed + * to. + */ + unsigned int current_block; + unsigned int index_within_block; + + /** + * Indices of the global + * element address at which + * we have to move on to + * another block when moving + * forward and + * backward. These indices + * are kept as a cache since + * this is much more + * efficient than always + * asking the parent object. + */ + unsigned int next_break_forward; + unsigned int next_break_backward; + + /** + * Move forward one element. + */ + void move_forward (); + + /** + * Move backward one element. + */ + void move_backward (); + + +#ifndef DEAL_II_NAMESP_TEMPL_FRIEND_BUG + /** + * Mark all other instances of + * this template as friends. In + * fact, we only need the + * inverse constness iterator + * as friend, but this is + * something that ISO C++ does + * not allow to specify. If we + * have detected a compiler bug + * during configuration of the + * library, use a workaround + * that works for this + * particular compiler, but is + * not ISO C++ conforming. + */ + template + friend class Iterator; +#else + friend class InverseConstnessIterator; +#endif + }; + } // namespace BlockVectorIterators +} // namespace internal + + + +/** + * A vector composed of several blocks each representing a vector of + * its own. + * + * The BlockVector is a collection of Vectors (e.g. of either deal.II Vector + * objects or PETScWrappers::Vector object). Each of the vectors inside can + * have a different size. + * + * The functionality of BlockVector includes everything a Vector can do, plus + * the access to a single Vector inside the BlockVector by + * block(i). It also has a complete random access iterator, just as + * the other Vector classes or the standard C++ library template + * std::vector. Therefore, all algorithms working on iterators also + * work with objects of this class. + * + * While this base class implements most of the functionality by dispatching + * calls to its member functions to the respective functions on each of the + * individual blocks, this class does not actually allocate some memory or + * change the size of vectors. For this, the constructors, assignment + * operators and reinit() functions of derived classes are responsible. This + * class only handles the common part that is independent of the actual vector + * type the block vector is built on. + * + * + * @section BlockVectorAccess Accessing individual blocks, and resizing vectors + * + * Apart from using this object as a whole, you can use each block + * separately as a vector, using the block() function. There + * is a single caveat: if you have changed the size of one or several + * blocks, you must call the function collect_sizes() of the block + * vector to update its internal structures. + * + * @attention Warning: If you change the sizes of single blocks + * without calling collect_sizes(), results may be unpredictable. The + * debug version does not check consistency here for performance + * reasons! + * + * @author Wolfgang Bangerth, Guido Kanschat, 1999, 2000, 2001, 2002, 2004 + */ +template +class BlockVectorBase +{ + public: + /** + * Typedef the type of the underlying + * vector. + */ + typedef VectorType BlockType; + + /* + * Declare standard types used in + * all containers. These types + * parallel those in the + * C++ standard + * libraries + * std::vector<...> + * class. This includes iterator + * types. + */ + typedef typename BlockType::value_type value_type; + typedef value_type *pointer; + typedef const value_type *const_pointer; + typedef internal::BlockVectorIterators::Iterator iterator; + typedef internal::BlockVectorIterators::Iterator const_iterator; + typedef typename BlockType::reference reference; + typedef const typename BlockType::reference const_reference; + + typedef std::size_t size_type; + + /** + * Default constructor. + */ + BlockVectorBase (); + + /** + * Update internal structures + * after resizing + * vectors. Whenever you reinited + * a block of a block vector, the + * internal data structures are + * corrupted. Therefore, you + * should call this function + * after al blocks got their new + * size. + */ + void collect_sizes (); + + /** + * Call the compress() function on all + * the subblocks of the matrix. + */ + void compress (); + + /** + * Set all entries to zero. Equivalent to + * v=0, but more obvious and + * faster. Note that this function does + * not change the size of the vector, + * unlike the STL's + * vector<>::clear function. + */ + void clear (); + + /** + * Access to a single block. + */ + BlockType & + block (const unsigned int i); + + /** + * Read-only access to a single block. + */ + const BlockType & + block (const unsigned int i) const; + + /** + * Return a reference on the + * object that describes the + * mapping between block and + * global indices. The use of + * this function is highly + * deprecated and it should + * vanish in one of the next + * versions + */ + const BlockIndices & + get_block_indices () const; + + /** + * Number of blocks. + */ + unsigned int n_blocks () const; + + /** + * Return dimension of the vector. This + * is the sum of the dimensions of all + * components. + */ + unsigned int size () const; + + /** + * Return an iterator pointing to + * the first element. + */ + iterator begin (); + + /** + * Return an iterator pointing to + * the first element of a + * constant block vector. + */ + const_iterator begin () const; + + /** + * Return an iterator pointing to + * the element past the end. + */ + iterator end (); + + /** + * Return an iterator pointing to + * the element past the end of a + * constant block vector. + */ + const_iterator end () const; + + /** + * Access components, returns U(i). + */ + value_type operator() (const unsigned int i) const; + + /** + * Access components, returns U(i) + * as a writeable reference. + */ + reference operator() (const unsigned int i); + + /** + * Copy operator: fill all components of + * the vector with the given scalar + * value. + */ + BlockVectorBase & operator = (const value_type s); + + /** + * Copy operator for arguments of the + * same type. + */ + BlockVectorBase & + operator= (const BlockVectorBase& V); + + /** + * Copy operator for template arguments + * of different types. + */ + template + BlockVectorBase & + operator= (const BlockVectorBase &V); + + /** + * Check for equality of two block vector + * types. This operation is only allowed + * if the two vectors already have the + * same block structure. + */ + template + bool + operator == (const BlockVectorBase &v) const; + + /** + * $U = U * V$: scalar product. + */ + value_type operator* (const BlockVectorBase& V) const; + + /** + * Return square of the $l_2$-norm. + */ + value_type norm_sqr () const; + + /** + * Return the mean value of the elements + * of this vector. + */ + value_type mean_value () const; + + /** + * Return the $l_1$-norm of the vector, + * i.e. the sum of the absolute values. + */ + value_type l1_norm () const; + + /** + * Return the $l_2$-norm of the vector, + * i.e. the square root of the sum of + * the squares of the elements. + */ + value_type l2_norm () const; + + /** + * Return the maximum absolute value of + * the elements of this vector, which is + * the $l_\infty$-norm of a vector. + */ + value_type linfty_norm () const; + + /** + * Return whether the vector contains only + * elements with value zero. This function + * is mainly for internal consistency + * check and should seldomly be used when + * not in debug mode since it uses quite + * some time. + */ + bool all_zero () const; + + /** + * Return @p true if the vector has no + * negative entries, i.e. all entries are + * zero or positive. This function is + * used, for example, to check whether + * refinement indicators are really all + * positive (or zero). + */ + bool is_non_negative () const; + + /** + * Addition operator. Fast equivalent to + * U.add(1, V). + */ + BlockVectorBase & + operator += (const BlockVectorBase &V); + + /** + * Subtraction operator. Fast equivalent + * to U.add(-1, V). + */ + BlockVectorBase & + operator -= (const BlockVectorBase &V); + + /** + * $U(0-DIM)+=s$. Addition of s + * to all components. Note that + * s is a scalar and not a + * vector. + */ + void add (const value_type s); + + /** + * U+=V. + * Simple vector addition, equal to the + * operator +=. + */ + void add (const BlockVectorBase& V); + + /** + * U+=a*V. + * Simple addition of a scaled vector. + */ + void add (const value_type a, const BlockVectorBase& V); + + /** + * U+=a*V+b*W. + * Multiple addition of scaled vectors. + */ + void add (const value_type a, const BlockVectorBase& V, + const value_type b, const BlockVectorBase& W); + + /** + * U=s*U+V. + * Scaling and simple vector addition. + */ + void sadd (const value_type s, const BlockVectorBase& V); + + /** + * U=s*U+a*V. + * Scaling and simple addition. + */ + void sadd (const value_type s, const value_type a, const BlockVectorBase& V); + + /** + * U=s*U+a*V+b*W. + * Scaling and multiple addition. + */ + void sadd (const value_type s, const value_type a, + const BlockVectorBase& V, + const value_type b, const BlockVectorBase& W); + + /** + * U=s*U+a*V+b*W+c*X. + * Scaling and multiple addition. + */ + void sadd (const value_type s, const value_type a, + const BlockVectorBase& V, + const value_type b, const BlockVectorBase& W, + const value_type c, const BlockVectorBase& X); + + /** + * Scale each element of the + * vector by a constant + * value. + */ + BlockVectorBase & operator *= (const value_type factor); + + /** + * Scale each element of the + * vector by the inverse of the + * given value. + */ + BlockVectorBase & operator /= (const value_type factor); + + /** + * Multiply each element of this + * vector by the corresponding + * element of v. + */ + template + void scale (const BlockVector2 &v); + + /** + * U=a*V. Assignment. + */ + template + void equ (const value_type a, const BlockVector2 &V); + + /** + * U=a*V+b*W. + * Replacing by sum. + */ + void equ (const value_type a, const BlockVectorBase& V, + const value_type b, const BlockVectorBase& W); + + /** + * Exception + */ + DeclException0 (ExcNotmatchingBlockSizes); + + protected: + /** + * Pointer to the array of components. + */ + std::vector components; + + /** + * Object managing the + * transformation between global + * indices and indices within the + * different blocks. + */ + BlockIndices block_indices; + + /** + * The number of blocks. This + * number is redundant to + * components.size() and stored + * here for convenience. + */ + unsigned int num_blocks; + + /** + * Make the iterator class a + * friend. We have to work around + * a compiler bug here again. + */ +#ifndef DEAL_II_NAMESP_TEMPL_FRIEND_BUG + template + friend class internal::BlockVectorIterators::Iterator; +#else + friend class iterator; + friend class const_iterator; +#endif + + template friend class BlockVectorBase; +}; + + +/*@}*/ + +/*----------------------- Inline functions ----------------------------------*/ + + +namespace internal +{ + namespace BlockVectorIterators + { + + template + inline + Iterator:: + Iterator (const Iterator &c) + : + parent (c.parent), + global_index (c.global_index), + current_block (c.current_block), + index_within_block (c.index_within_block), + next_break_forward (c.next_break_forward), + next_break_backward (c.next_break_backward) + {} + + + + template + inline + Iterator:: + Iterator (const InverseConstnessIterator &c) + : + parent (const_cast(c.parent)), + global_index (c.global_index), + current_block (c.current_block), + index_within_block (c.index_within_block), + next_break_forward (c.next_break_forward), + next_break_backward (c.next_break_backward) + { + // if constness==false, then the + // constness of the iterator we + // got is true and we are trying + // to cast away the + // constness. disallow this + Assert (constness==true, ExcCastingAwayConstness()); + } + + + + template + inline + Iterator:: + Iterator (BlockVector &parent, + const unsigned int global_index, + const unsigned int current_block, + const unsigned int index_within_block, + const unsigned int next_break_forward, + const unsigned int next_break_backward) + : + parent (&parent), + global_index (global_index), + current_block (current_block), + index_within_block (index_within_block), + next_break_forward (next_break_forward), + next_break_backward (next_break_backward) + { + } + + + + template + inline + Iterator & + Iterator:: + operator = (const Iterator &c) + { + parent = c.parent; + global_index = c.global_index; + index_within_block = c.index_within_block; + current_block = c.current_block; + next_break_forward = c.next_break_forward; + next_break_backward = c.next_break_backward; + + return *this; + } + + + + template + inline + typename Iterator::dereference_type + Iterator::operator * () const + { + return parent->block(current_block)(index_within_block); + } + + + + template + inline + typename Iterator::dereference_type + Iterator::operator [] (const difference_type d) const + { + // if the index pointed to is + // still within the block we + // currently point into, then we + // can save the computation of + // the block + if ((global_index+d >= next_break_backward) && + (global_index+d <= next_break_forward)) + return parent->block(current_block)(index_within_block + d); + + // if the index is not within the + // block of the block vector into + // which we presently point, then + // there is no way: we have to + // search for the block. this can + // be done through the parent + // class as well. + return (*parent)(global_index+d); + } + + + + template + inline + Iterator & + Iterator::operator ++ () + { + move_forward (); + return *this; + } + + + + template + inline + Iterator + Iterator::operator ++ (int) + { + const Iterator old_value = *this; + move_forward (); + return old_value; + } + + + + template + inline + Iterator & + Iterator::operator -- () + { + move_backward (); + return *this; + } + + + + template + inline + Iterator + Iterator::operator -- (int) + { + const Iterator old_value = *this; + move_backward (); + return old_value; + } + + + + template + inline + bool + Iterator:: + operator == (const Iterator &i) const + { + Assert (parent == i.parent, ExcPointerToDifferentVectors()); + + return (global_index == i.global_index); + } + + + + template + inline + bool + Iterator:: + operator == (const InverseConstnessIterator &i) const + { + Assert (parent == i.parent, ExcPointerToDifferentVectors()); + + return (global_index == i.global_index); + } + + + + template + inline + bool + Iterator:: + operator != (const Iterator &i) const + { + Assert (parent == i.parent, ExcPointerToDifferentVectors()); + + return (global_index != i.global_index); + } + + + + template + inline + bool + Iterator:: + operator != (const InverseConstnessIterator &i) const + { + Assert (parent == i.parent, ExcPointerToDifferentVectors()); + + return (global_index != i.global_index); + } + + + + template + inline + bool + Iterator:: + operator < (const Iterator &i) const + { + Assert (parent == i.parent, ExcPointerToDifferentVectors()); + + return (global_index < i.global_index); + } + + + + template + inline + bool + Iterator:: + operator < (const InverseConstnessIterator &i) const + { + Assert (parent == i.parent, ExcPointerToDifferentVectors()); + + return (global_index < i.global_index); + } + + + + template + inline + bool + Iterator:: + operator <= (const Iterator &i) const + { + Assert (parent == i.parent, ExcPointerToDifferentVectors()); + + return (global_index <= i.global_index); + } + + + + template + inline + bool + Iterator:: + operator <= (const InverseConstnessIterator &i) const + { + Assert (parent == i.parent, ExcPointerToDifferentVectors()); + + return (global_index <= i.global_index); + } + + + + template + inline + bool + Iterator:: + operator > (const Iterator &i) const + { + Assert (parent == i.parent, ExcPointerToDifferentVectors()); + + return (global_index > i.global_index); + } + + + + template + inline + bool + Iterator:: + operator > (const InverseConstnessIterator &i) const + { + Assert (parent == i.parent, ExcPointerToDifferentVectors()); + + return (global_index > i.global_index); + } + + + + template + inline + bool + Iterator:: + operator >= (const Iterator &i) const + { + Assert (parent == i.parent, ExcPointerToDifferentVectors()); + + return (global_index >= i.global_index); + } + + + + template + inline + bool + Iterator:: + operator >= (const InverseConstnessIterator &i) const + { + Assert (parent == i.parent, ExcPointerToDifferentVectors()); + + return (global_index >= i.global_index); + } + + + + template + inline + typename Iterator::difference_type + Iterator:: + operator - (const Iterator &i) const + { + Assert (parent == i.parent, ExcPointerToDifferentVectors()); + + return (static_cast(global_index) - + static_cast(i.global_index)); + } + + + + template + inline + typename Iterator::difference_type + Iterator:: + operator - (const InverseConstnessIterator &i) const + { + Assert (parent == i.parent, ExcPointerToDifferentVectors()); + + return (static_cast(global_index) - + static_cast(i.global_index)); + } + + + + template + inline + Iterator + Iterator:: + operator + (const difference_type &d) const + { + // if the index pointed to is + // still within the block we + // currently point into, then we + // can save the computation of + // the block + if ((global_index+d >= next_break_backward) && + (global_index+d <= next_break_forward)) + return Iterator (*parent, global_index+d, current_block, + index_within_block+d, + next_break_forward, next_break_backward); + else + // outside present block, so + // have to seek new block + // anyway + return Iterator (*parent, global_index+d); + } + + + + template + inline + Iterator + Iterator:: + operator - (const difference_type &d) const + { + // if the index pointed to is + // still within the block we + // currently point into, then we + // can save the computation of + // the block + if ((global_index-d >= next_break_backward) && + (global_index-d <= next_break_forward)) + return Iterator (*parent, global_index-d, current_block, + index_within_block-d, + next_break_forward, next_break_backward); + else + // outside present block, so + // have to seek new block + // anyway + return Iterator (*parent, global_index-d); + } + + + + template + inline + Iterator & + Iterator:: + operator += (const difference_type &d) + { + // if the index pointed to is + // still within the block we + // currently point into, then we + // can save the computation of + // the block + if ((global_index+d >= next_break_backward) && + (global_index+d <= next_break_forward)) + { + global_index += d; + index_within_block += d; + } + else + // outside present block, so + // have to seek new block + // anyway + *this = Iterator (*parent, global_index+d); + + return *this; + } + + + + template + inline + Iterator & + Iterator:: + operator -= (const difference_type &d) + { + // if the index pointed to is + // still within the block we + // currently point into, then we + // can save the computation of + // the block + if ((global_index-d >= next_break_backward) && + (global_index-d <= next_break_forward)) + { + global_index -= d; + index_within_block -= d; + } + else + // outside present block, so + // have to seek new block + // anyway + *this = Iterator (*parent, global_index-d); + + return *this; + } + + + template + Iterator:: + Iterator (BlockVector &parent, + const unsigned global_index) + : + parent (&parent), + global_index (global_index) + { + // find which block we are + // in. for this, take into + // account that it happens at + // times that people want to + // initialize iterators + // past-the-end + if (global_index < parent.size()) + { + const std::pair + indices = parent.block_indices.global_to_local(global_index); + current_block = indices.first; + index_within_block = indices.second; + + next_break_backward + = parent.block_indices.local_to_global (current_block, 0); + next_break_forward + = (parent.block_indices.local_to_global (current_block, 0) + +parent.block_indices.block_size(current_block)-1); + } + else + // past the end. only have one + // value for this + { + this->global_index = parent.size (); + current_block = parent.n_blocks(); + index_within_block = 0; + next_break_backward = global_index; + next_break_forward = deal_II_numbers::invalid_unsigned_int; + }; + } + + + + template + void + Iterator::move_forward () + { + if (global_index != next_break_forward) + ++index_within_block; + else + { + // ok, we traverse a boundary + // between blocks: + index_within_block = 0; + ++current_block; + + // break backwards is now old + // break forward + next_break_backward = next_break_forward+1; + + // compute new break forward + if (current_block < parent->block_indices.size()) + next_break_forward + += parent->block_indices.block_size(current_block); + else + // if we are beyond the end, + // then move the next + // boundary arbitrarily far + // away + next_break_forward = deal_II_numbers::invalid_unsigned_int; + }; + + ++global_index; + } + + + + template + void + Iterator::move_backward () + { + if (global_index != next_break_backward) + --index_within_block; + else + if (current_block != 0) + { + // ok, we traverse a boundary + // between blocks: + --current_block; + index_within_block = parent->block_indices.block_size(current_block)-1; + + // break forwards is now old + // break backward + next_break_forward = next_break_backward-1; + + // compute new break forward + next_break_backward + -= parent->block_indices.block_size (current_block); + } + else + // current block was 0, we now + // get into unspecified terrain + { + --current_block; + index_within_block = deal_II_numbers::invalid_unsigned_int; + next_break_forward = 0; + next_break_backward = 0; + }; + + --global_index; + } + + + } // namespace BlockVectorIterators + +} //namespace internal + + +template +inline +BlockVectorBase::BlockVectorBase () +{} + + + +template +inline +unsigned int +BlockVectorBase::size () const +{ + return block_indices.total_size(); +} + + + +template +inline +unsigned int +BlockVectorBase::n_blocks () const +{ + return num_blocks; +} + + +template +inline +typename BlockVectorBase::BlockType & +BlockVectorBase::block (const unsigned int i) +{ + Assert(i +inline +const typename BlockVectorBase::BlockType & +BlockVectorBase::block (const unsigned int i) const +{ + Assert(i +inline +const BlockIndices& +BlockVectorBase::get_block_indices () const +{ + return block_indices; +} + + +template +inline +void +BlockVectorBase::collect_sizes () +{ + std::vector sizes (num_blocks); + + for (unsigned int i=0; i +inline +void +BlockVectorBase::compress () +{ + for (unsigned int i=0; i +inline +void +BlockVectorBase::clear () +{ + for (unsigned int i=0;i +inline +typename BlockVectorBase::iterator +BlockVectorBase::begin() +{ + return iterator(*this, 0U); +} + + + +template +inline +typename BlockVectorBase::const_iterator +BlockVectorBase::begin() const +{ + return const_iterator(*this, 0U); +} + + +template +inline +typename BlockVectorBase::iterator +BlockVectorBase::end() +{ + return iterator(*this, size()); +} + + + +template +inline +typename BlockVectorBase::const_iterator +BlockVectorBase::end() const +{ + return const_iterator(*this, size()); +} + + + +template +bool +BlockVectorBase::all_zero () const +{ + for (unsigned int i=0;i +bool +BlockVectorBase::is_non_negative () const +{ + for (unsigned int i=0;i +typename BlockVectorBase::value_type +BlockVectorBase:: +operator * (const BlockVectorBase& v) const +{ + Assert (num_blocks == v.num_blocks, + ExcDimensionMismatch(num_blocks, v.num_blocks)); + + value_type sum = 0.; + for (unsigned int i=0;i +typename BlockVectorBase::value_type +BlockVectorBase::norm_sqr () const +{ + value_type sum = 0.; + for (unsigned int i=0;i +typename BlockVectorBase::value_type +BlockVectorBase::mean_value () const +{ + value_type sum = 0.; + for (unsigned int i=0;i +typename BlockVectorBase::value_type +BlockVectorBase::l1_norm () const +{ + value_type sum = 0.; + for (unsigned int i=0;i +typename BlockVectorBase::value_type +BlockVectorBase::l2_norm () const +{ + return std::sqrt(norm_sqr()); +} + + + +template +typename BlockVectorBase::value_type +BlockVectorBase::linfty_norm () const +{ + value_type sum = 0.; + for (unsigned int i=0;i +BlockVectorBase & +BlockVectorBase::operator += (const BlockVectorBase& v) +{ + add (v); + return *this; +} + + + +template +BlockVectorBase & +BlockVectorBase::operator -= (const BlockVectorBase& v) +{ + Assert (num_blocks == v.num_blocks, + ExcDimensionMismatch(num_blocks, v.num_blocks)); + + for (unsigned int i=0;i +void BlockVectorBase::add (const value_type a) +{ + for (unsigned int i=0;i +void BlockVectorBase::add (const BlockVectorBase& v) +{ + Assert (num_blocks == v.num_blocks, + ExcDimensionMismatch(num_blocks, v.num_blocks)); + + for (unsigned int i=0;i +void BlockVectorBase::add (const value_type a, + const BlockVectorBase& v) +{ + Assert (num_blocks == v.num_blocks, + ExcDimensionMismatch(num_blocks, v.num_blocks)); + + for (unsigned int i=0;i +void BlockVectorBase::add (const value_type a, + const BlockVectorBase& v, + const value_type b, + const BlockVectorBase& w) +{ + Assert (num_blocks == v.num_blocks, + ExcDimensionMismatch(num_blocks, v.num_blocks)); + Assert (num_blocks == w.num_blocks, + ExcDimensionMismatch(num_blocks, w.num_blocks)); + + + for (unsigned int i=0;i +void BlockVectorBase::sadd (const value_type x, + const BlockVectorBase& v) +{ + Assert (num_blocks == v.num_blocks, + ExcDimensionMismatch(num_blocks, v.num_blocks)); + + for (unsigned int i=0;i +void BlockVectorBase::sadd (const value_type x, const value_type a, + const BlockVectorBase& v) +{ + Assert (num_blocks == v.num_blocks, + ExcDimensionMismatch(num_blocks, v.num_blocks)); + + for (unsigned int i=0;i +void BlockVectorBase::sadd (const value_type x, const value_type a, + const BlockVectorBase& v, + const value_type b, + const BlockVectorBase& w) +{ + Assert (num_blocks == v.num_blocks, + ExcDimensionMismatch(num_blocks, v.num_blocks)); + Assert (num_blocks == w.num_blocks, + ExcDimensionMismatch(num_blocks, w.num_blocks)); + + for (unsigned int i=0;i +void BlockVectorBase::sadd (const value_type x, const value_type a, + const BlockVectorBase& v, + const value_type b, + const BlockVectorBase& w, + const value_type c, + const BlockVectorBase& y) +{ + Assert (num_blocks == v.num_blocks, + ExcDimensionMismatch(num_blocks, v.num_blocks)); + Assert (num_blocks == w.num_blocks, + ExcDimensionMismatch(num_blocks, w.num_blocks)); + Assert (num_blocks == y.num_blocks, + ExcDimensionMismatch(num_blocks, y.num_blocks)); + + for (unsigned int i=0;i +template +void BlockVectorBase::scale (const BlockVector2 &v) +{ + Assert (num_blocks == v.num_blocks, + ExcDimensionMismatch(num_blocks, v.num_blocks)); + for (unsigned int i=0;i +void BlockVectorBase::equ (const value_type a, + const BlockVectorBase& v, + const value_type b, + const BlockVectorBase& w) +{ + Assert (num_blocks == v.num_blocks, + ExcDimensionMismatch(num_blocks, v.num_blocks)); + Assert (num_blocks == w.num_blocks, + ExcDimensionMismatch(num_blocks, w.num_blocks)); + + for (unsigned int i=0;i +template +void BlockVectorBase::equ (const value_type a, + const BlockVector2 &v) +{ + Assert (num_blocks == v.num_blocks, + ExcDimensionMismatch(num_blocks, v.num_blocks)); + + for (unsigned int i=0;i +BlockVectorBase& +BlockVectorBase::operator = (const value_type s) +{ + for (unsigned int i=0;i +BlockVectorBase& +BlockVectorBase::operator = (const BlockVectorBase& v) +{ + Assert (num_blocks == v.num_blocks, + ExcDimensionMismatch(num_blocks, v.num_blocks)); + + for (unsigned int i=0;i +template +BlockVectorBase& +BlockVectorBase::operator = (const BlockVectorBase &v) +{ + Assert (num_blocks == v.num_blocks, + ExcDimensionMismatch(num_blocks, v.num_blocks)); + + for (unsigned int i=0;i +template +inline +bool +BlockVectorBase:: +operator == (const BlockVectorBase &v) const +{ + Assert (block_indices == v.block_indices, + ExcNotmatchingBlockSizes()); + + for (unsigned int i=0;i +inline +BlockVectorBase & +BlockVectorBase::operator *= (const value_type factor) +{ + for (unsigned int i=0;i +inline +BlockVectorBase & +BlockVectorBase::operator /= (const value_type factor) +{ + for (unsigned int i=0;i +inline +typename BlockVectorBase::value_type +BlockVectorBase::operator() (const unsigned int i) const +{ + const std::pair local_index + = block_indices.global_to_local (i); + return components[local_index.first](local_index.second); +} + + + +template +inline +typename BlockVectorBase::reference +BlockVectorBase::operator() (const unsigned int i) +{ + const std::pair local_index + = block_indices.global_to_local (i); + return components[local_index.first](local_index.second); +} + + +#endif diff --git a/deal.II/lac/include/lac/petsc_block_sparse_matrix.h b/deal.II/lac/include/lac/petsc_block_sparse_matrix.h new file mode 100644 index 0000000000..42e4b4ee9d --- /dev/null +++ b/deal.II/lac/include/lac/petsc_block_sparse_matrix.h @@ -0,0 +1,201 @@ +//---------------------------- petsc_block_sparse_matrix.h --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- petsc_block_sparse_matrix.h --------------------------- +#ifndef __deal2__petsc_block_sparse_matrix_h +#define __deal2__petsc_block_sparse_matrix_h + + +#include +#include +#include +#include +#include + + + +/*! @addtogroup PETSc + *@{ + */ + +namespace PETScWrappers +{ + + +/** + * Blocked sparse matrix based on the PETScWrappers::SparseMatrix class. This + * class implements the functions that are specific to the PETSc SparseMatrix + * base objects for a blocked sparse matrix, and leaves the actual work + * relaying most of the calls to the individual blocks to the functions + * implemented in the base class. See there also for a description of when + * this class is useful. + * + * In contrast to the deal.II-type SparseMatrix class, the PETSc matrices do + * not have external objects for the sparsity patterns. Thus, one does not + * determine the size of the individual blocks of a block matrix of this type + * by attaching a block sparsity pattern, but by calling reinit() to set the + * number of blocks and then by setting the size of each block separately. In + * order to fix the data structures of the block matrix, it is then necessary + * to let it know that we have changed the sizes of the underlying + * matrices. For this, one has to call the collect_sizes() function, for much + * the same reason as is documented with the BlockSparsityPattern class. + * + * @author Wolfgang Bangerth, 2004 + */ + class BlockSparseMatrix : public BlockMatrixBase + { + public: + /** + * Typedef the base class for simpler + * access to its own typedefs. + */ + typedef BlockMatrixBase BaseClass; + + /** + * Typedef the type of the underlying + * matrix. + */ + typedef BaseClass::BlockType BlockType; + + /** + * Import the typedefs from the base + * class. + */ + typedef BaseClass::value_type value_type; + typedef BaseClass::pointer pointer; + typedef BaseClass::const_pointer const_pointer; + typedef BaseClass::reference reference; + typedef BaseClass::const_reference const_reference; + typedef BaseClass::size_type size_type; + typedef BaseClass::iterator iterator; + typedef BaseClass::const_iterator const_iterator; + + /** + * Constructor; initializes the + * matrix to be empty, without + * any structure, i.e. the + * matrix is not usable at + * all. This constructor is + * therefore only useful for + * matrices which are members of + * a class. All other matrices + * should be created at a point + * in the data flow where all + * necessary information is + * available. + * + * You have to initialize the + * matrix before usage with + * reinit(BlockSparsityPattern). The + * number of blocks per row and + * column are then determined by + * that function. + */ + BlockSparseMatrix (); + + /** + * Destructor. + */ + ~BlockSparseMatrix (); + + /** + * Pseudo copy operator only copying + * empty objects. The sizes of the block + * matrices need to be the same. + */ + BlockSparseMatrix & + operator = (const BlockSparseMatrix &); + + /** + * This operator assigns a scalar to a + * matrix. Since this does usually not + * make much sense (should we set all + * matrix entries to this value? Only + * the nonzero entries of the sparsity + * pattern?), this operation is only + * allowed if the actual value to be + * assigned is zero. This operator only + * exists to allow for the obvious + * notation matrix=0, which + * sets all elements of the matrix to + * zero, but keep the sparsity pattern + * previously used. + */ + BlockSparseMatrix & + operator = (const double d); + + /** + * Resize the matrix, by setting + * the number of block rows and + * columns. This deletes all + * blocks and replaces them by + * unitialized ones, i.e. ones + * for which also the sizes are + * not yet set. You have to do + * that by calling the @p reinit + * functions of the blocks + * themselves. Do not forget to + * call @p collect_sizes after + * that on this object. + * + * The reason that you have to + * set sizes of the blocks + * yourself is that the sizes may + * be varying, the maximum number + * of elements per row may be + * varying, etc. It is simpler + * not to reproduce the interface + * of the @p SparsityPattern + * class here but rather let the + * user call whatever function + * she desires. + */ + void reinit (const unsigned int n_block_rows, + const unsigned int n_block_columns); + + + /** + * This function collects the + * sizes of the sub-objects and + * stores them in internal + * arrays, in order to be able to + * relay global indices into the + * matrix to indices into the + * subobjects. You *must* call + * this function each time after + * you have changed the size of + * the sub-objects. + */ + void collect_sizes (); + + /** + * Exception + */ + DeclException4 (ExcIncompatibleRowNumbers, + int, int, int, int, + << "The blocks [" << arg1 << ',' << arg2 << "] and [" + << arg3 << ',' << arg4 << "] have differing row numbers."); + /** + * Exception + */ + DeclException4 (ExcIncompatibleColNumbers, + int, int, int, int, + << "The blocks [" << arg1 << ',' << arg2 << "] and [" + << arg3 << ',' << arg4 << "] have differing column numbers."); + }; + + + +/*@}*/ +} + + +#endif // __deal2__petsc_block_sparse_matrix_h diff --git a/deal.II/lac/include/lac/petsc_block_vector.h b/deal.II/lac/include/lac/petsc_block_vector.h new file mode 100644 index 0000000000..651a381206 --- /dev/null +++ b/deal.II/lac/include/lac/petsc_block_vector.h @@ -0,0 +1,489 @@ +//---------------------------- petsc_block_vector.h --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- petsc_block_vector.h --------------------------- +#ifndef __deal2__petsc_block_vector_h +#define __deal2__petsc_block_vector_h + + +#include +#include +#include +#include +#include +#include + +#ifdef DEAL_II_USE_PETSC + + +/*! @addtogroup PETScWrappers + *@{ + */ + +namespace PETScWrappers +{ + +/** + * An implementation of block vectors based on the vector class implemented in + * PETScWrappers. While the base class provides for most of the interface, + * this class handles the actual allocation of vectors and provides functions + * that are specific to the underlying vector type. + * + * @author Wolfgang Bangerth, 2004 + */ + class BlockVector : public BlockVectorBase + { + public: + /** + * Typedef the base class for simpler + * access to its own typedefs. + */ + typedef BlockVectorBase BaseClass; + + /** + * Typedef the type of the underlying + * vector. + */ + typedef BaseClass::BlockType BlockType; + + /** + * Import the typedefs from the base + * class. + */ + typedef BaseClass::value_type value_type; + typedef BaseClass::pointer pointer; + typedef BaseClass::const_pointer const_pointer; + typedef BaseClass::reference reference; + typedef BaseClass::const_reference const_reference; + typedef BaseClass::size_type size_type; + typedef BaseClass::iterator iterator; + typedef BaseClass::const_iterator const_iterator; + + /** + * Constructor. There are three + * ways to use this + * constructor. First, without + * any arguments, it generates + * an object with no + * blocks. Given one argument, + * it initializes num_blocks + * blocks, but these blocks have + * size zero. The third variant + * finally initializes all + * blocks to the same size + * block_size. + * + * Confer the other constructor + * further down if you intend to + * use blocks of different + * sizes. + */ + explicit BlockVector (const unsigned int num_blocks = 0, + const unsigned int block_size = 0); + + /** + * Copy-Constructor. Dimension set to + * that of V, all components are copied + * from V + */ + BlockVector (const BlockVector &V); + + /** + * Copy-constructor: copy the values + * from a PETSc wrapper parallel block + * vector class. + * + * + * Note that due to the communication + * model of MPI, @em all processes have + * to actually perform this operation, + * even if they do not use the + * result. It is not sufficient if only + * one processor tries to copy the + * elements from the other processors + * over to its own process space. + */ + explicit BlockVector (const MPI::BlockVector &v); + + /** + * Constructor. Set the number of + * blocks to n.size() and + * initialize each block with + * n[i] zero elements. + */ + BlockVector (const std::vector &n); + + /** + * Constructor. Set the number of + * blocks to + * n.size(). Initialize the + * vector with the elements + * pointed to by the range of + * iterators given as second and + * third argument. Apart from the + * first argument, this + * constructor is in complete + * analogy to the respective + * constructor of the + * std::vector class, but the + * first argument is needed in + * order to know how to subdivide + * the block vector into + * different blocks. + */ + template + BlockVector (const std::vector &n, + const InputIterator first, + const InputIterator end); + + /** + * Destructor. Clears memory + */ + ~BlockVector (); + + /** + * Copy operator: fill all components of + * the vector with the given scalar + * value. + */ + BlockVector & operator = (const value_type s); + + /** + * Copy operator for arguments of the + * same type. + */ + BlockVector & + operator= (const BlockVector &V); + + /** + * Copy all the elements of the + * parallel block vector @p v into this + * local vector. Note that due to the + * communication model of MPI, @em all + * processes have to actually perform + * this operation, even if they do not + * use the result. It is not sufficient + * if only one processor tries to copy + * the elements from the other + * processors over to its own process + * space. + */ + BlockVector & + operator = (const MPI::BlockVector &v); + + /** + * Reinitialize the BlockVector to + * contain num_blocks blocks of + * size block_size each. + * + * If fast==false, the vector + * is filled with zeros. + */ + void reinit (const unsigned int num_blocks, + const unsigned int block_size, + const bool fast = false); + + /** + * Reinitialize the BlockVector such + * that it contains + * block_sizes.size() + * blocks. Each block is reinitialized + * to dimension + * block_sizes[i]. + * + * If the number of blocks is the + * same as before this function + * was called, all vectors remain + * the same and reinit() is + * called for each vector. + * + * If fast==false, the vector + * is filled with zeros. + * + * Note that you must call this + * (or the other reinit() + * functions) function, rather + * than calling the reinit() + * functions of an individual + * block, to allow the block + * vector to update its caches of + * vector sizes. If you call + * reinit() on one of the + * blocks, then subsequent + * actions on this object may + * yield unpredictable results + * since they may be routed to + * the wrong block. + */ + void reinit (const std::vector &N, + const bool fast=false); + + /** + * Change the dimension to that + * of the vector V. The same + * applies as for the other + * reinit() function. + * + * The elements of V are not + * copied, i.e. this function is + * the same as calling reinit + * (V.size(), fast). + * + * Note that you must call this + * (or the other reinit() + * functions) function, rather + * than calling the reinit() + * functions of an individual + * block, to allow the block + * vector to update its caches of + * vector sizes. If you call + * reinit() of one of the + * blocks, then subsequent + * actions of this object may + * yield unpredictable results + * since they may be routed to + * the wrong block. + */ + void reinit (const BlockVector &V, + const bool fast=false); + + /** + * Swap the contents of this + * vector and the other vector + * v. One could do this + * operation with a temporary + * variable and copying over the + * data elements, but this + * function is significantly more + * efficient since it only swaps + * the pointers to the data of + * the two vectors and therefore + * does not need to allocate + * temporary storage and move + * data around. + * + * Limitation: right now this + * function only works if both + * vectors have the same number + * of blocks. If needed, the + * numbers of blocks should be + * exchanged, too. + * + * This function is analog to the + * the swap() function of all C++ + * standard containers. Also, + * there is a global function + * swap(u,v) that simply calls + * u.swap(v), again in analogy + * to standard functions. + */ + void swap (BlockVector &v); + + /** + * Exception + */ + DeclException0 (ExcIteratorRangeDoesNotMatchVectorSize); + }; + +/*@}*/ + +/*----------------------- Inline functions ----------------------------------*/ + + + + inline + BlockVector::BlockVector (const unsigned int n_blocks, + const unsigned int block_size) + { + reinit (n_blocks, block_size); + } + + + + inline + BlockVector::BlockVector (const std::vector &n) + { + reinit (n, false); + } + + + inline + BlockVector::BlockVector (const BlockVector& v) + : + BlockVectorBase () + { + components.resize (v.num_blocks); + block_indices = v.block_indices; + num_blocks = v.num_blocks; + + for (unsigned int i=0; i () + { + components.resize (v.get_block_indices().size()); + block_indices = v.get_block_indices(); + num_blocks = block_indices.size(); + + for (unsigned int i=0; i + BlockVector::BlockVector (const std::vector &n, + const InputIterator first, + const InputIterator end) + { + // first set sizes of blocks, but + // don't initialize them as we will + // copy elements soon + reinit (n, true); + InputIterator start = first; + for (unsigned int b=0; b(n[b])); + std::copy (start, end, this->block(b).begin()); + start = end; + }; + Assert (start == end, ExcIteratorRangeDoesNotMatchVectorSize()); + } + + + + inline + BlockVector & + BlockVector::operator = (const value_type s) + { + BaseClass::operator = (s); + return *this; + } + + + + inline + BlockVector & + BlockVector::operator = (const BlockVector &v) + { + BaseClass::operator = (v); + return *this; + } + + + + inline + BlockVector & + BlockVector::operator = (const MPI::BlockVector &v) + { + BaseClass::operator = (v); + return *this; + } + + + + inline + BlockVector::~BlockVector () + {} + + + inline + void + BlockVector::reinit (const unsigned int n_bl, + const unsigned int bl_sz, + const bool fast) + { + std::vector n(n_bl, bl_sz); + reinit(n, fast); + } + + + + inline + void + BlockVector::reinit (const std::vector &n, + const bool fast) + { + block_indices.reinit (n); + num_blocks = n.size(); + if (components.size() != num_blocks) + components.resize(num_blocks); + + for (unsigned int i=0; i #include +#include #ifdef DEAL_II_USE_PETSC @@ -280,7 +281,7 @@ namespace PETScWrappers * @ingroup PETScWrappers * @author Wolfgang Bangerth, 2004 */ - class MatrixBase + class MatrixBase : public Subscriptor { public: /** @@ -431,6 +432,47 @@ namespace PETScWrappers */ unsigned int n () const; + /** + * Return the local dimension of the + * matrix, i.e. the number of rows + * stored on the present MPI + * process. For sequential matrices, + * this number is the same as m(), + * but for parallel matrices it may be + * smaller. + * + * To figure out which elements + * exactly are stored locally, + * use local_range(). + */ + unsigned int local_size () const; + + /** + * Return a pair of indices + * indicating which rows of + * this matrix are stored + * locally. The first number is + * the index of the first + * row stored, the second + * the index of the one past + * the last one that is stored + * locally. If this is a + * sequential matrix, then the + * result will be the pair + * (0,m()), otherwise it will be + * a pair (i,i+n), where + * n=local_size(). + */ + std::pair + local_range () const; + + /** + * Return whether @p index is + * in the local range or not, + * see also local_range(). + */ + bool in_local_range (const unsigned int index) const; + /** * Return the number of nonzero * elements of this @@ -898,6 +940,22 @@ namespace PETScWrappers Assert (r < m(), ExcIndexRange(r, 0, m())); return const_iterator(this, r+1, 0); } + + + + inline + bool + MatrixBase::in_local_range (const unsigned int index) const + { + int begin, end; + const int ierr = MatGetOwnershipRange (static_cast(matrix), + &begin, &end); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + return ((index >= static_cast(begin)) && + (index < static_cast(end))); + } + /// @endif } diff --git a/deal.II/lac/include/lac/petsc_parallel_block_sparse_matrix.h b/deal.II/lac/include/lac/petsc_parallel_block_sparse_matrix.h new file mode 100644 index 0000000000..733d19fe9f --- /dev/null +++ b/deal.II/lac/include/lac/petsc_parallel_block_sparse_matrix.h @@ -0,0 +1,198 @@ +//---------------------------- petsc_parallel_block_sparse_matrix.h --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- petsc_parallel_block_sparse_matrix.h --------------------------- +#ifndef __deal2__petsc_parallel_block_sparse_matrix_h +#define __deal2__petsc_parallel_block_sparse_matrix_h + + +#include +#include +#include +#include +#include + + + +/*! @addtogroup PETSc + *@{ + */ + +namespace PETScWrappers +{ + namespace MPI + { + + +/** + * Blocked sparse matrix based on the PETScWrappers::SparseMatrix class. This + * class implements the functions that are specific to the PETSc SparseMatrix + * base objects for a blocked sparse matrix, and leaves the actual work + * relaying most of the calls to the individual blocks to the functions + * implemented in the base class. See there also for a description of when + * this class is useful. + * + * In contrast to the deal.II-type SparseMatrix class, the PETSc matrices do + * not have external objects for the sparsity patterns. Thus, one does not + * determine the size of the individual blocks of a block matrix of this type + * by attaching a block sparsity pattern, but by calling reinit() to set the + * number of blocks and then by setting the size of each block separately. In + * order to fix the data structures of the block matrix, it is then necessary + * to let it know that we have changed the sizes of the underlying + * matrices. For this, one has to call the collect_sizes() function, for much + * the same reason as is documented with the BlockSparsityPattern class. + * + * @author Wolfgang Bangerth, 2004 + */ + class BlockSparseMatrix : public BlockMatrixBase + { + public: + /** + * Typedef the base class for simpler + * access to its own typedefs. + */ + typedef BlockMatrixBase BaseClass; + + /** + * Typedef the type of the underlying + * matrix. + */ + typedef BaseClass::BlockType BlockType; + + /** + * Import the typedefs from the base + * class. + */ + typedef BaseClass::value_type value_type; + typedef BaseClass::pointer pointer; + typedef BaseClass::const_pointer const_pointer; + typedef BaseClass::reference reference; + typedef BaseClass::const_reference const_reference; + typedef BaseClass::size_type size_type; + typedef BaseClass::iterator iterator; + typedef BaseClass::const_iterator const_iterator; + + /** + * Constructor; initializes the + * matrix to be empty, without + * any structure, i.e. the + * matrix is not usable at + * all. This constructor is + * therefore only useful for + * matrices which are members of + * a class. All other matrices + * should be created at a point + * in the data flow where all + * necessary information is + * available. + * + * You have to initialize the + * matrix before usage with + * reinit(BlockSparsityPattern). The + * number of blocks per row and + * column are then determined by + * that function. + */ + BlockSparseMatrix (); + + /** + * Destructor. + */ + ~BlockSparseMatrix (); + + /** + * Pseudo copy operator only copying + * empty objects. The sizes of the + * block matrices need to be the + * same. + */ + BlockSparseMatrix & + operator = (const BlockSparseMatrix &); + + /** + * This operator assigns a scalar to + * a matrix. Since this does usually + * not make much sense (should we set + * all matrix entries to this value? + * Only the nonzero entries of the + * sparsity pattern?), this operation + * is only allowed if the actual + * value to be assigned is zero. This + * operator only exists to allow for + * the obvious notation + * matrix=0, which sets all + * elements of the matrix to zero, + * but keep the sparsity pattern + * previously used. + */ + BlockSparseMatrix & + operator = (const PetscScalar d); + + /** + * Resize the matrix, by setting + * the number of block rows and + * columns. This deletes all + * blocks and replaces them by + * unitialized ones, i.e. ones + * for which also the sizes are + * not yet set. You have to do + * that by calling the @p reinit + * functions of the blocks + * themselves. Do not forget to + * call @p collect_sizes after + * that on this object. + * + * The reason that you have to + * set sizes of the blocks + * yourself is that the sizes may + * be varying, the maximum number + * of elements per row may be + * varying, etc. It is simpler + * not to reproduce the interface + * of the @p SparsityPattern + * class here but rather let the + * user call whatever function + * she desires. + */ + void reinit (const unsigned int n_block_rows, + const unsigned int n_block_columns); + + /** + * This function collects the + * sizes of the sub-objects and + * stores them in internal + * arrays, in order to be able to + * relay global indices into the + * matrix to indices into the + * subobjects. You *must* call + * this function each time after + * you have changed the size of + * the sub-objects. + */ + void collect_sizes (); + + /** + * Return a reference to the MPI + * communicator object in use with + * this vector. + */ + const MPI_Comm & get_mpi_communicator () const; + }; + + + +/*@}*/ + } + +} + + +#endif // __deal2__petsc_parallel_block_sparse_matrix_h diff --git a/deal.II/lac/include/lac/petsc_parallel_block_vector.h b/deal.II/lac/include/lac/petsc_parallel_block_vector.h new file mode 100644 index 0000000000..4fb6d07795 --- /dev/null +++ b/deal.II/lac/include/lac/petsc_parallel_block_vector.h @@ -0,0 +1,485 @@ +//---------------------------- petsc_parallel_block_vector.h --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- petsc_parallel_block_vector.h --------------------------- +#ifndef __deal2__petsc_parallel_block_vector_h +#define __deal2__petsc_parallel_block_vector_h + + +#include +#include +#include +#include +#include + +#ifdef DEAL_II_USE_PETSC + + +/*! @addtogroup PETScWrappers + *@{ + */ + +namespace PETScWrappers +{ + // forward declaration + class BlockVector; + + namespace MPI + { + +/** + * An implementation of block vectors based on the parallel vector class + * implemented in PETScWrappers. While the base class provides for most of the + * interface, this class handles the actual allocation of vectors and provides + * functions that are specific to the underlying vector type. + * + * The model of distribution of data is such that each of the blocks is + * distributed across all MPI processes named in the MPI communicator. I.e. we + * don't just distribute the whole vector, but each component. In the + * constructors and reinit() functions, one therefore not only has to specify + * the sizes of the individual blocks, but also the number of elements of each + * of these blocks to be stored on the local process. + * + * @author Wolfgang Bangerth, 2004 + */ + class BlockVector : public BlockVectorBase + { + public: + /** + * Typedef the base class for simpler + * access to its own typedefs. + */ + typedef BlockVectorBase BaseClass; + + /** + * Typedef the type of the underlying + * vector. + */ + typedef BaseClass::BlockType BlockType; + + /** + * Import the typedefs from the base + * class. + */ + typedef BaseClass::value_type value_type; + typedef BaseClass::pointer pointer; + typedef BaseClass::const_pointer const_pointer; + typedef BaseClass::reference reference; + typedef BaseClass::const_reference const_reference; + typedef BaseClass::size_type size_type; + typedef BaseClass::iterator iterator; + typedef BaseClass::const_iterator const_iterator; + + /** + * Default constructor. Generate an + * empty vector without any blocks. + */ + BlockVector (); + + /** + * Constructor. Generate a block + * vector with @p n_blocks blocks, + * each of which is a parallel + * vector across @p communicator + * with @p block_size elements of + * which @p local_size elements are + * stored on the present process. + */ + explicit BlockVector (const unsigned int n_blocks, + const MPI_Comm &communicator, + const unsigned int block_size, + const unsigned int local_size); + + /** + * Copy-Constructor. Set all the + * properties of the parallel vector + * to those of the given argument and + * copy the elements. + */ + BlockVector (const BlockVector &V); + + /** + * Constructor. Set the number of + * blocks to + * block_sizes.size() and + * initialize each block with + * block_sizes[i] zero + * elements. The individual blocks + * are distributed across the given + * communicator, and each store + * local_elements[i] + * elements on the present process. + */ + BlockVector (const std::vector &block_sizes, + const MPI_Comm &communicator, + const std::vector &local_elements); + + /** + * Destructor. Clears memory + */ + ~BlockVector (); + + /** + * Copy operator: fill all components + * of the vector that are locally + * stored with the given scalar value. + */ + BlockVector & operator = (const value_type s); + + /** + * Copy operator for arguments of the + * same type. + */ + BlockVector & + operator= (const BlockVector &V); + + /** + * Copy the given sequential + * (non-distributed) block vector + * into the present parallel block + * vector. It is assumed that they + * have the same size, and this + * operation does not change the + * partitioning of the parallel + * vectors by which its elements are + * distributed across several MPI + * processes. What this operation + * therefore does is to copy that + * chunk of the given vector @p v + * that corresponds to elements of + * the target vector that are stored + * locally, and copies them, for each + * of the individual blocks of this + * object. Elements that are not + * stored locally are not touched. + * + * This being a parallel vector, you + * must make sure that @em all + * processes call this function at + * the same time. It is not possible + * to change the local part of a + * parallel vector on only one + * process, independent of what other + * processes do, with this function. + */ + BlockVector & + operator = (const PETScWrappers::BlockVector &v); + + /** + * Reinitialize the BlockVector to + * contain @p n_blocks of size @p + * block_size, each of which stores + * @p local_size elements + * locally. The @p communicator + * argument denotes which MPI channel + * each of these blocks shall + * communicate. + * + * If fast==false, the vector + * is filled with zeros. + */ + void reinit (const unsigned int n_blocks, + const MPI_Comm &communicator, + const unsigned int block_size, + const unsigned int local_size, + const bool fast = false); + + /** + * Reinitialize the BlockVector such + * that it contains + * block_sizes.size() + * blocks. Each block is + * reinitialized to dimension + * block_sizes[i]. Each of + * them stores + * local_sizes[i] elements + * on the present process. + * + * If the number of blocks is the + * same as before this function + * was called, all vectors remain + * the same and reinit() is + * called for each vector. + * + * If fast==false, the vector + * is filled with zeros. + * + * Note that you must call this + * (or the other reinit() + * functions) function, rather + * than calling the reinit() + * functions of an individual + * block, to allow the block + * vector to update its caches of + * vector sizes. If you call + * reinit() of one of the + * blocks, then subsequent + * actions on this object may + * yield unpredictable results + * since they may be routed to + * the wrong block. + */ + void reinit (const std::vector &block_sizes, + const MPI_Comm &communicator, + const std::vector &local_sizes, + const bool fast=false); + + /** + * Change the dimension to that + * of the vector V. The same + * applies as for the other + * reinit() function. + * + * The elements of V are not + * copied, i.e. this function is + * the same as calling reinit + * (V.size(), fast). + * + * Note that you must call this + * (or the other reinit() + * functions) function, rather + * than calling the reinit() + * functions of an individual + * block, to allow the block + * vector to update its caches of + * vector sizes. If you call + * reinit() on one of the + * blocks, then subsequent + * actions on this object may + * yield unpredictable results + * since they may be routed to + * the wrong block. + */ + void reinit (const BlockVector &V, + const bool fast=false); + + /** + * Return a reference to the MPI + * communicator object in use with + * this vector. + */ + const MPI_Comm & get_mpi_communicator () const; + + /** + * Swap the contents of this + * vector and the other vector + * v. One could do this + * operation with a temporary + * variable and copying over the + * data elements, but this + * function is significantly more + * efficient since it only swaps + * the pointers to the data of + * the two vectors and therefore + * does not need to allocate + * temporary storage and move + * data around. + * + * Limitation: right now this + * function only works if both + * vectors have the same number + * of blocks. If needed, the + * numbers of blocks should be + * exchanged, too. + * + * This function is analog to the + * the swap() function of all C++ + * standard containers. Also, + * there is a global function + * swap(u,v) that simply calls + * u.swap(v), again in analogy + * to standard functions. + */ + void swap (BlockVector &v); + + /** + * Exception + */ + DeclException0 (ExcIteratorRangeDoesNotMatchVectorSize); + /** + * Exception + */ + DeclException0 (ExcNonMatchingBlockVectors); + }; + +/*@}*/ + +/*----------------------- Inline functions ----------------------------------*/ + + + inline + BlockVector::BlockVector () + {} + + + + inline + BlockVector::BlockVector (const unsigned int n_blocks, + const MPI_Comm &communicator, + const unsigned int block_size, + const unsigned int local_size) + { + reinit (n_blocks, communicator, block_size, local_size); + } + + + + inline + BlockVector::BlockVector (const std::vector &block_sizes, + const MPI_Comm &communicator, + const std::vector &local_elements) + { + reinit (block_sizes, communicator, local_elements, false); + } + + + inline + BlockVector::BlockVector (const BlockVector& v) + : + BlockVectorBase () + { + components.resize (v.num_blocks); + block_indices = v.block_indices; + num_blocks = v.num_blocks; + + for (unsigned int i=0; i(n_blocks, block_size), + communicator, + std::vector(n_blocks, local_size), + fast); + } + + + + inline + void + BlockVector::reinit (const std::vector &block_sizes, + const MPI_Comm &communicator, + const std::vector &local_sizes, + const bool fast) + { + block_indices.reinit (block_sizes); + num_blocks = block_sizes.size(); + if (components.size() != num_blocks) + components.resize(num_blocks); + + for (unsigned int i=0; i +#include #include @@ -43,6 +44,30 @@ namespace PETScWrappers * There are a number of comments on the communication model as well as access * to individual elements in the documentation to the parallel vector * class. These comments apply here as well. + * + * + * @section Partitioning Partitioning of matrices + * + * PETSc partitions parallel matrices so that each MPI process "owns" a + * certain number of rows (i.e. only this process stores the respective + * entries in these rows). The number of rows each process owns has to be + * passed to the constructors and reinit() functions via the argument @p + * local_rows. The individual values passed as @p local_rows on all the MPI + * processes of course have to add up to the global number of rows of the + * matrix. + * + * In addition to this, PETSc also partitions the rectangular chunk of the + * matrix it owns (i.e. the @p local_rows times n() elements in the matrix), + * so that matrix vector multiplications can be performed efficiently. This + * column-partitioning therefore has to match the partitioning of the vectors + * with which the matrix is multiplied, just as the row-partitioning has to + * match the partitioning of destination vectors. This partitioning is passed + * to the constructors and reinit() functions through the @p local_columns + * variable, which again has to add up to the global number of columns in the + * matrix. The name @p local_columns may be named inappropriately since it + * does not reflect that only these columns are stored locally, but it + * reflects the fact that these are the columns for which the elements of + * incoming vectors are stored locally. * * @ingroup PETScWrappers * @author Wolfgang Bangerth, 2004 @@ -51,23 +76,28 @@ namespace PETScWrappers { public: /** - * Default constructor. Create an empty - * matrix. + * Default constructor. Create an + * empty matrix. */ SparseMatrix (); /** - * Create a sparse matrix of dimensions - * @p m times @p n, with an - * initial guess of - * @p n_nonzero_per_row nonzero - * elements per row. PETSc is able to - * cope with the situation that more - * than this number of elements is - * later allocated for a row, but this + * Create a sparse matrix of + * dimensions @p m times @p n, with + * an initial guess of @p + * n_nonzero_per_row nonzero elements + * per row. PETSc is able to cope + * with the situation that more than + * this number of elements are later + * allocated for a row, but this * involves copying data, and is thus * expensive. * + * For the meaning of the @p + * local_row and @p local_columns + * parameters, see the class + * documentation. + * * The @p is_symmetric flag determines * whether we should tell PETSc that * the matrix is going to be symmetric @@ -85,43 +115,111 @@ namespace PETScWrappers const unsigned int m, const unsigned int n, const unsigned int local_rows, + const unsigned int local_columns, const unsigned int n_nonzero_per_row, const bool is_symmetric = false); /** - * Initialize a rectangular matrix with - * @p m rows and @p n - * columns. The maximal number of - * nonzero entries for each row - * separately is given by the - * @p row_lengths array. + * Initialize a rectangular matrix + * with @p m rows and @p n columns. + * The maximal number of nonzero + * entries for each row separately is + * given by the @p row_lengths array. * - * Just as for the other constructors: - * PETSc is able to cope with the - * situation that more than this number - * of elements is later allocated for a - * row, but this involves copying data, - * and is thus expensive. + * For the meaning of the @p + * local_row and @p local_columns + * parameters, see the class + * documentation. + * + * Just as for the other + * constructors: PETSc is able to + * cope with the situation that more + * than this number of elements are + * later allocated for a row, but + * this involves copying data, and is + * thus expensive. * - * The @p is_symmetric flag determines - * whether we should tell PETSc that - * the matrix is going to be symmetric - * (as indicated by the call - * MatSetOption(mat, + * The @p is_symmetric flag + * determines whether we should tell + * PETSc that the matrix is going to + * be symmetric (as indicated by the + * call MatSetOption(mat, * MAT_SYMMETRIC). Note that the - * PETSc documentation states that one - * cannot form an ILU decomposition of - * a matrix for which this flag has - * been set to @p true, only an - * ICC. The default value of this flag - * is @p false. + * PETSc documentation states that + * one cannot form an ILU + * decomposition of a matrix for + * which this flag has been set to @p + * true, only an ICC. The default + * value of this flag is @p false. */ SparseMatrix (const MPI_Comm &communicator, const unsigned int m, const unsigned int n, const unsigned int local_rows, + const unsigned int local_columns, const std::vector &row_lengths, const bool is_symmetric = false); + + /** + * Initialize using the given + * sparsity pattern with + * communication happening over the + * provided @p communicator. + * + * For the meaning of the @p + * local_row and @p local_columns + * parameters, see the class + * documentation. + * + * Note that PETSc can be very slow + * if you do not provide it with a + * good estimate of the lengths of + * rows. Using the present function + * is a very efficient way to do + * this, as it uses the exact number + * of nonzero entries for each row of + * the matrix by using the given + * sparsity pattern argument. If the + * @p preset_nonzero_locations flag + * is @p true, this function in + * addition not only sets the correct + * row sizes up front, but also + * pre-allocated the correct nonzero + * entries in the matrix. + * + * PETsc allows to later add + * additional nonzero entries to a + * matrix, by simply writing to these + * elements. However, this will then + * lead to additional memory + * allocations which are very + * inefficient and will greatly slow + * down your program. It is therefore + * significantly more efficient to + * get memory allocation right from + * the start. + * + * Despite the fact that it would + * seem to be an obvious win, setting + * the @p preset_nonzero_locations + * flag to @p true doesn't seem to + * accelerate program. Rather on the + * contrary, it seems to be able to + * slow down entire programs + * somewhat. This is suprising, since + * we can use efficient function + * calls into PETSc that allow to + * create multiple entries at once; + * nevertheless, given the fact that + * it is inefficient, the respective + * flag has a default value equal to + * @p false. + */ + SparseMatrix (const MPI_Comm &communicator, + const CompressedSparsityPattern &sparsity_pattern, + const unsigned int local_rows, + const unsigned int local_columns, + const bool preset_nonzero_locations = false); /** * This operator assigns a scalar to @@ -153,6 +251,7 @@ namespace PETScWrappers const unsigned int m, const unsigned int n, const unsigned int local_rows, + const unsigned int local_columns, const unsigned int n_nonzero_per_row, const bool is_symmetric = false); @@ -168,9 +267,73 @@ namespace PETScWrappers const unsigned int m, const unsigned int n, const unsigned int local_rows, + const unsigned int local_columns, const std::vector &row_lengths, const bool is_symmetric = false); + /** + * Initialize using the given + * sparsity pattern with + * communication happening over the + * provided @p communicator. + * + * Note that PETSc can be very slow + * if you do not provide it with a + * good estimate of the lengths of + * rows. Using the present function + * is a very efficient way to do + * this, as it uses the exact number + * of nonzero entries for each row of + * the matrix by using the given + * sparsity pattern argument. If the + * @p preset_nonzero_locations flag + * is @p true, this function in + * addition not only sets the correct + * row sizes up front, but also + * pre-allocated the correct nonzero + * entries in the matrix. + * + * PETsc allows to later add + * additional nonzero entries to a + * matrix, by simply writing to these + * elements. However, this will then + * lead to additional memory + * allocations which are very + * inefficient and will greatly slow + * down your program. It is therefore + * significantly more efficient to + * get memory allocation right from + * the start. + * + * Despite the fact that it would + * seem to be an obvious win, setting + * the @p preset_nonzero_locations + * flag to @p true doesn't seem to + * accelerate program. Rather on the + * contrary, it seems to be able to + * slow down entire programs + * somewhat. This is suprising, since + * we can use efficient function + * calls into PETSc that allow to + * create multiple entries at once; + * nevertheless, given the fact that + * it is inefficient, the respective + * flag has a default value equal to + * @p false. + */ + void reinit (const MPI_Comm &communicator, + const CompressedSparsityPattern &sparsity_pattern, + const unsigned int local_rows, + const unsigned int local_columns, + const bool preset_nonzero_locations = false); + + /** + * Return a reference to the MPI + * communicator object in use with + * this matrix. + */ + const MPI_Comm & get_mpi_communicator () const; + /** * Exception */ @@ -189,14 +352,16 @@ namespace PETScWrappers /** * Do the actual work for the - * respective reinit() function and the - * matching constructor, i.e. create a - * matrix. Getting rid of the previous - * matrix is left to the caller. + * respective reinit() function and + * the matching constructor, + * i.e. create a matrix. Getting rid + * of the previous matrix is left to + * the caller. */ void do_reinit (const unsigned int m, const unsigned int n, const unsigned int local_rows, + const unsigned int local_columns, const unsigned int n_nonzero_per_row, const bool is_symmetric = false); @@ -206,10 +371,29 @@ namespace PETScWrappers void do_reinit (const unsigned int m, const unsigned int n, const unsigned int local_rows, + const unsigned int local_columns, const std::vector &row_lengths, const bool is_symmetric = false); + + /** + * Same as previous functions. + */ + void do_reinit (const CompressedSparsityPattern &sparsity_pattern, + const unsigned int local_rows, + const unsigned int local_columns, + const bool preset_nonzero_locations); }; + + +// -------- template and inline functions ---------- + + inline + const MPI_Comm & + SparseMatrix::get_mpi_communicator () const + { + return communicator; + } } } diff --git a/deal.II/lac/include/lac/petsc_parallel_vector.h b/deal.II/lac/include/lac/petsc_parallel_vector.h index 566f7d96a7..4905fcf282 100644 --- a/deal.II/lac/include/lac/petsc_parallel_vector.h +++ b/deal.II/lac/include/lac/petsc_parallel_vector.h @@ -24,6 +24,13 @@ #include +// forward declaration +template class Vector; + + +/*! @addtogroup PETSc + *@{ + */ namespace PETScWrappers { /** @@ -335,6 +342,13 @@ namespace PETScWrappers void reinit (const Vector &v, const bool fast = false); + /** + * Return a reference to the MPI + * communicator object in use with + * this vector. + */ + const MPI_Comm & get_mpi_communicator () const; + protected: /** * Create a vector of length @@ -361,6 +375,21 @@ namespace PETScWrappers // ------------------ template and inline functions ------------- +/** + * Global function @p swap which overloads the default implementation + * of the C++ standard library which uses a temporary object. The + * function simply exchanges the data of the two vectors. + * + * @author Wolfgang Bangerth, 2004 + */ + inline + void swap (Vector &u, Vector &v) + { + u.swap (v); + } + + + template Vector::Vector (const MPI_Comm &communicator, const ::Vector &v, @@ -459,7 +488,17 @@ namespace PETScWrappers compress (); return *this; - } + } + + + + inline + const MPI_Comm & + Vector::get_mpi_communicator () const + { + return communicator; + } + /// @endif } } diff --git a/deal.II/lac/include/lac/petsc_solver.h b/deal.II/lac/include/lac/petsc_solver.h index c527127e26..b96f313c17 100644 --- a/deal.II/lac/include/lac/petsc_solver.h +++ b/deal.II/lac/include/lac/petsc_solver.h @@ -19,6 +19,7 @@ #ifdef DEAL_II_USE_PETSC +#include #include @@ -37,6 +38,28 @@ namespace PETScWrappers * classes simply set the right flags to select one solver or another, or to * set certain parameters for individual solvers. * + * One of the gotchas of PETSc is that -- in particular in MPI mode -- it + * often does not produce very helpful error messages. In order to save + * other users some time in searching a hard to track down error, here is + * one situation and the error message one gets there: + * when you don't specify an MPI communicator to your solver's constructor. In + * this case, you will get an error of the following form from each of your + * parallel processes: + * @verbatim + * [1]PETSC ERROR: PCSetVector() line 1173 in src/ksp/pc/interface/precon.c + * [1]PETSC ERROR: Arguments must have same communicators! + * [1]PETSC ERROR: Different communicators in the two objects: Argument # 1 and 2! + * [1]PETSC ERROR: KSPSetUp() line 195 in src/ksp/ksp/interface/itfunc.c + * @endverbatim + * + * This error, on which one can spend a very long time figuring out what + * exactly goes wrong, results from not specifying an MPI communicator. Note + * that the communicator @em must match that of the matrix and all vectors in + * the linear system which we want to solve. Aggravating the situation is the + * fact that the default argument to the solver classes, @p PETSC_COMM_SELF, + * is the appropriate argument for the sequential case (which is why it is the + * default argument), so this error only shows up in parallel mode. + * * @ingroup PETScWrappers * @author Wolfgang Bangerth, 2004 */ @@ -55,10 +78,12 @@ namespace PETScWrappers * right hand side object of the solve * to be done with this * solver. Otherwise, PETSc will - * generate hard to track down errors. + * generate hard to track down errors, + * see the documentation of the + * SolverBase class. */ - SolverBase (SolverControl &cn, - MPI_Comm &mpi_communicator); + SolverBase (SolverControl &cn, + const MPI_Comm &mpi_communicator); /** * Destructor. @@ -78,7 +103,7 @@ namespace PETScWrappers solve (const MatrixBase &A, VectorBase &x, const VectorBase &b, - const PreconditionerBase &preconditioner) const; + const PreconditionerBase &preconditioner); /** @@ -113,7 +138,7 @@ namespace PETScWrappers * Copy of the MPI communicator object * to be used for the solver. */ - MPI_Comm mpi_communicator; + const MPI_Comm mpi_communicator; /** * Function that takes a Krylov @@ -140,6 +165,64 @@ namespace PETScWrappers const double residual_norm, KSPConvergedReason *reason, void *solver_control); + + /** + * A structure that contains the PETSc + * solver and preconditioner + * objects. This object is preserved + * between subsequent calls to the + * solver if the same preconditioner is + * used as in the previous solver + * step. This may save some computation + * time, if setting up a preconditioner + * is expensive, such as in the case of + * an ILU for example. + * + * The actual declaration of this class + * is complicated by the fact that + * PETSc changed its solver interface + * completely and incompatibly between + * versions 2.1.6 and 2.2.0 :-( + * + * Objects of this type are explicitly + * created, but are destroyed when the + * surrounding solver object goes out + * of scope, or when we assign a new + * value to the pointer to this + * object. The respective *Destroy + * functions are therefore written into + * the destructor of this object, even + * though the object does not have a + * constructor. + */ + struct SolverData + { + /** + * Destructor + */ + ~SolverData (); + +#if (PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR < 2) + /** + * A PETSc solver object. + */ + SLES sles; +#endif + /** + * Objects for Krylov subspace + * solvers and preconditioners. + */ + KSP ksp; + PC pc; + }; + + /** + * Pointer to an object that stores the + * solver context. This is recreated in + * the main solver routine if + * necessary. + */ + boost::shared_ptr solver_data; }; @@ -202,10 +285,12 @@ namespace PETScWrappers * right hand side object of the solve * to be done with this * solver. Otherwise, PETSc will - * generate hard to track down errors. + * generate hard to track down errors, + * see the documentation of the + * SolverBase class. */ SolverRichardson (SolverControl &cn, - MPI_Comm &mpi_communicator = PETSC_COMM_SELF, + const MPI_Comm &mpi_communicator = PETSC_COMM_SELF, const AdditionalData &data = AdditionalData(1)); protected: @@ -268,10 +353,12 @@ namespace PETScWrappers * right hand side object of the solve * to be done with this * solver. Otherwise, PETSc will - * generate hard to track down errors. + * generate hard to track down errors, + * see the documentation of the + * SolverBase class. */ SolverChebychev (SolverControl &cn, - MPI_Comm &mpi_communicator = PETSC_COMM_SELF, + const MPI_Comm &mpi_communicator = PETSC_COMM_SELF, const AdditionalData &data = AdditionalData()); protected: @@ -334,10 +421,12 @@ namespace PETScWrappers * right hand side object of the solve * to be done with this * solver. Otherwise, PETSc will - * generate hard to track down errors. + * generate hard to track down errors, + * see the documentation of the + * SolverBase class. */ SolverCG (SolverControl &cn, - MPI_Comm &mpi_communicator = PETSC_COMM_SELF, + const MPI_Comm &mpi_communicator = PETSC_COMM_SELF, const AdditionalData &data = AdditionalData()); protected: @@ -400,10 +489,12 @@ namespace PETScWrappers * right hand side object of the solve * to be done with this * solver. Otherwise, PETSc will - * generate hard to track down errors. + * generate hard to track down errors, + * see the documentation of the + * SolverBase class. */ SolverBiCG (SolverControl &cn, - MPI_Comm &mpi_communicator = PETSC_COMM_SELF, + const MPI_Comm &mpi_communicator = PETSC_COMM_SELF, const AdditionalData &data = AdditionalData()); protected: @@ -483,10 +574,12 @@ namespace PETScWrappers * right hand side object of the solve * to be done with this * solver. Otherwise, PETSc will - * generate hard to track down errors. + * generate hard to track down errors, + * see the documentation of the + * SolverBase class. */ SolverGMRES (SolverControl &cn, - MPI_Comm &mpi_communicator = PETSC_COMM_SELF, + const MPI_Comm &mpi_communicator = PETSC_COMM_SELF, const AdditionalData &data = AdditionalData(30)); protected: @@ -549,10 +642,12 @@ namespace PETScWrappers * right hand side object of the solve * to be done with this * solver. Otherwise, PETSc will - * generate hard to track down errors. + * generate hard to track down errors, + * see the documentation of the + * SolverBase class. */ SolverBicgstab (SolverControl &cn, - MPI_Comm &mpi_communicator = PETSC_COMM_SELF, + const MPI_Comm &mpi_communicator = PETSC_COMM_SELF, const AdditionalData &data = AdditionalData()); protected: @@ -615,10 +710,12 @@ namespace PETScWrappers * right hand side object of the solve * to be done with this * solver. Otherwise, PETSc will - * generate hard to track down errors. + * generate hard to track down errors, + * see the documentation of the + * SolverBase class. */ SolverCGS (SolverControl &cn, - MPI_Comm &mpi_communicator = PETSC_COMM_SELF, + const MPI_Comm &mpi_communicator = PETSC_COMM_SELF, const AdditionalData &data = AdditionalData()); protected: @@ -681,10 +778,12 @@ namespace PETScWrappers * right hand side object of the solve * to be done with this * solver. Otherwise, PETSc will - * generate hard to track down errors. + * generate hard to track down errors, + * see the documentation of the + * SolverBase class. */ SolverTFQMR (SolverControl &cn, - MPI_Comm &mpi_communicator = PETSC_COMM_SELF, + const MPI_Comm &mpi_communicator = PETSC_COMM_SELF, const AdditionalData &data = AdditionalData()); protected: @@ -752,10 +851,12 @@ namespace PETScWrappers * right hand side object of the solve * to be done with this * solver. Otherwise, PETSc will - * generate hard to track down errors. + * generate hard to track down errors, + * see the documentation of the + * SolverBase class. */ SolverTCQMR (SolverControl &cn, - MPI_Comm &mpi_communicator = PETSC_COMM_SELF, + const MPI_Comm &mpi_communicator = PETSC_COMM_SELF, const AdditionalData &data = AdditionalData()); protected: @@ -818,10 +919,12 @@ namespace PETScWrappers * right hand side object of the solve * to be done with this * solver. Otherwise, PETSc will - * generate hard to track down errors. + * generate hard to track down errors, + * see the documentation of the + * SolverBase class. */ SolverCR (SolverControl &cn, - MPI_Comm &mpi_communicator = PETSC_COMM_SELF, + const MPI_Comm &mpi_communicator = PETSC_COMM_SELF, const AdditionalData &data = AdditionalData()); protected: @@ -884,10 +987,12 @@ namespace PETScWrappers * right hand side object of the solve * to be done with this * solver. Otherwise, PETSc will - * generate hard to track down errors. + * generate hard to track down errors, + * see the documentation of the + * SolverBase class. */ SolverLSQR (SolverControl &cn, - MPI_Comm &mpi_communicator = PETSC_COMM_SELF, + const MPI_Comm &mpi_communicator = PETSC_COMM_SELF, const AdditionalData &data = AdditionalData()); protected: diff --git a/deal.II/lac/include/lac/petsc_sparse_matrix.h b/deal.II/lac/include/lac/petsc_sparse_matrix.h index 62668e18bb..76291d259e 100644 --- a/deal.II/lac/include/lac/petsc_sparse_matrix.h +++ b/deal.II/lac/include/lac/petsc_sparse_matrix.h @@ -20,6 +20,7 @@ #ifdef DEAL_II_USE_PETSC #include +#include #include @@ -109,6 +110,57 @@ namespace PETScWrappers const std::vector &row_lengths, const bool is_symmetric = false); + /** + * Initialize a sparse matrix using the + * given sparsity pattern. + * + * Note that PETSc can be very slow + * if you do not provide it with a + * good estimate of the lengths of + * rows. Using the present function + * is a very efficient way to do + * this, as it uses the exact number + * of nonzero entries for each row of + * the matrix by using the given + * sparsity pattern argument. If the + * @p preset_nonzero_locations flag + * is @p true, this function in + * addition not only sets the correct + * row sizes up front, but also + * pre-allocated the correct nonzero + * entries in the matrix. + * + * PETsc allows to later add + * additional nonzero entries to a + * matrix, by simply writing to these + * elements. However, this will then + * lead to additional memory + * allocations which are very + * inefficient and will greatly slow + * down your program. It is therefore + * significantly more efficient to + * get memory allocation right from + * the start. + * + * Despite the fact that it would + * seem to be an obvious win, setting + * the @p preset_nonzero_locations + * flag to @p true doesn't seem to + * accelerate program. Rather on the + * contrary, it seems to be able to + * slow down entire programs + * somewhat. This is suprising, since + * we can use efficient function + * calls into PETSc that allow to + * create multiple entries at once; + * nevertheless, given the fact that + * it is inefficient, the respective + * flag has a default value equal to + * @p false. + */ + SparseMatrix (const CompressedSparsityPattern &sparsity_pattern, + const bool preset_nonzero_locations = false); + /** * This operator assigns a scalar to * a matrix. Since this does usually @@ -153,6 +205,57 @@ namespace PETScWrappers const std::vector &row_lengths, const bool is_symmetric = false); + /** + * Initialize a sparse matrix using the + * given sparsity pattern. + * + * Note that PETSc can be very slow + * if you do not provide it with a + * good estimate of the lengths of + * rows. Using the present function + * is a very efficient way to do + * this, as it uses the exact number + * of nonzero entries for each row of + * the matrix by using the given + * sparsity pattern argument. If the + * @p preset_nonzero_locations flag + * is @p true, this function in + * addition not only sets the correct + * row sizes up front, but also + * pre-allocated the correct nonzero + * entries in the matrix. + * + * PETsc allows to later add + * additional nonzero entries to a + * matrix, by simply writing to these + * elements. However, this will then + * lead to additional memory + * allocations which are very + * inefficient and will greatly slow + * down your program. It is therefore + * significantly more efficient to + * get memory allocation right from + * the start. + * + * Despite the fact that it would + * seem to be an obvious win, setting + * the @p preset_nonzero_locations + * flag to @p true doesn't seem to + * accelerate program. Rather on the + * contrary, it seems to be able to + * slow down entire programs + * somewhat. This is suprising, since + * we can use efficient function + * calls into PETSc that allow to + * create multiple entries at once; + * nevertheless, given the fact that + * it is inefficient, the respective + * flag has a default value equal to + * @p false. + */ + void reinit (const CompressedSparsityPattern &sparsity_pattern, + const bool preset_nonzero_locations = false); + private: /** @@ -174,6 +277,12 @@ namespace PETScWrappers const unsigned int n, const std::vector &row_lengths, const bool is_symmetric = false); + + /** + * Same as previous function. + */ + void do_reinit (const CompressedSparsityPattern &sparsity_pattern, + const bool preset_nonzero_locations); }; } diff --git a/deal.II/lac/include/lac/petsc_vector.h b/deal.II/lac/include/lac/petsc_vector.h index 6b6168b9c3..e4d1ede1ba 100644 --- a/deal.II/lac/include/lac/petsc_vector.h +++ b/deal.II/lac/include/lac/petsc_vector.h @@ -77,11 +77,21 @@ namespace PETScWrappers * Copy-constructor the values from a * PETSc wrapper vector class. */ - explicit Vector (const Vector &v); + Vector (const Vector &v); /** - * Copy-constructor the values from a - * PETSc wrapper parallel vector class. + * Copy-constructor: copy the values + * from a PETSc wrapper parallel vector + * class. + * + * Note that due to the communication + * model of MPI, @em all processes have + * to actually perform this operation, + * even if they do not use the + * result. It is not sufficient if only + * one processor tries to copy the + * elements from the other processors + * over to its own process space. */ explicit Vector (const MPI::Vector &v); @@ -93,7 +103,7 @@ namespace PETScWrappers /** * Copy all the elements of the - * parallel vector @arg v into this + * parallel vector @p v into this * local vector. Note that due to the * communication model of MPI, @em all * processes have to actually perform @@ -185,6 +195,21 @@ namespace PETScWrappers // ------------------ template and inline functions ------------- +/** + * Global function @p swap which overloads the default implementation + * of the C++ standard library which uses a temporary object. The + * function simply exchanges the data of the two vectors. + * + * @author Wolfgang Bangerth, 2004 + */ + inline + void swap (Vector &u, Vector &v) + { + u.swap (v); + } + + + template Vector::Vector (const ::Vector &v) { diff --git a/deal.II/lac/include/lac/petsc_vector_base.h b/deal.II/lac/include/lac/petsc_vector_base.h index af3a6cb8f0..fa917458dc 100644 --- a/deal.II/lac/include/lac/petsc_vector_base.h +++ b/deal.II/lac/include/lac/petsc_vector_base.h @@ -213,8 +213,9 @@ namespace PETScWrappers * standard libraries vector<...> * class. */ - typedef PetscScalar value_type; - typedef size_t size_type; + typedef PetscScalar value_type; + typedef size_t size_type; + typedef internal::VectorReference reference; /** * Default constructor. It doesn't do @@ -325,7 +326,8 @@ namespace PETScWrappers * a pair (i,i+n), where * n=local_size(). */ - std::pair local_range () const; + std::pair + local_range () const; /** * Return whether @p index is @@ -338,7 +340,7 @@ namespace PETScWrappers * Provide access to a given element, * both read and write. */ - internal::VectorReference + reference operator () (const unsigned int index); /** @@ -585,7 +587,6 @@ namespace PETScWrappers */ void swap (VectorBase &v); - /** * Conversion operator to gain access * to the underlying PETSc type. If you @@ -865,7 +866,8 @@ namespace PETScWrappers &begin, &end); AssertThrow (ierr == 0, ExcPETScError(ierr)); - return (index >= static_cast(begin)) && (index < static_cast(end)); + return ((index >= static_cast(begin)) && + (index < static_cast(end))); } diff --git a/deal.II/lac/include/lac/vector.h b/deal.II/lac/include/lac/vector.h index b147e99251..194fcc5ae2 100644 --- a/deal.II/lac/include/lac/vector.h +++ b/deal.II/lac/include/lac/vector.h @@ -104,23 +104,32 @@ class Vector Vector (const Vector &v); -// note: I disabled this function for the time being, since gcc2.95 -// does not respect the "explicit" keyword for template constructors. -// this leads to unwanted conversions and in some places to automatically -// generated temporaries, where this is not a good idea. [WB] -// /** -// * Copy constructor taking a vector of -// * another data type. This will fail if -// * there is no conversion path from -// * @p OtherNumber to @p Number. Note that -// * you may lose accuracy when copying -// * to a vector with data elements with -// * less accuracy. -// */ -// template -// explicit -// Vector (const Vector &v); - +#ifndef DEAL_II_EXPLICIT_CONSTRUCTOR_BUG + /** + * Copy constructor taking a vector of + * another data type. This will fail if + * there is no conversion path from + * @p OtherNumber to @p Number. Note that + * you may lose accuracy when copying + * to a vector with data elements with + * less accuracy. + * + * Older versions of gcc did not honor + * the @p explicit keyword on template + * constructors. In such cases, it is + * easy to accidentally write code that + * can be very inefficient, since the + * compiler starts performing hidden + * conversions. To avoid this, this + * function is disabled if we have + * detected a broken compiler during + * configuration. + */ + template + explicit + Vector (const Vector &v); +#endif + #ifdef DEAL_II_USE_PETSC /** * Another copy constructor: copy the diff --git a/deal.II/lac/include/lac/vector.templates.h b/deal.II/lac/include/lac/vector.templates.h index 07ea7fd246..0710bbcd04 100644 --- a/deal.II/lac/include/lac/vector.templates.h +++ b/deal.II/lac/include/lac/vector.templates.h @@ -60,23 +60,25 @@ Vector::Vector (const Vector& v) } -// see the .h file for why this function was disabled -// -// template -// template -// Vector::Vector (const Vector& v) : -// dim(v.size()), -// maxdim(v.size()), -// val(0) -// { -// if (dim != 0) -// { -// val = new Number[maxdim]; -// Assert (val != 0, ExcOutOfMemory()); -// copy (v.begin(), v.end(), begin()); -// } -// } +#ifndef DEAL_II_EXPLICIT_CONSTRUCTOR_BUG +template +template +Vector::Vector (const Vector& v) + : + dim(v.size()), + maxdim(v.size()), + val(0) +{ + if (dim != 0) + { + val = new Number[maxdim]; + Assert (val != 0, ExcOutOfMemory()); + copy (v.begin(), v.end(), begin()); + } +} + +#endif #ifdef DEAL_II_USE_PETSC diff --git a/deal.II/lac/source/block_sparse_matrix.cc b/deal.II/lac/source/block_sparse_matrix.cc index a2d4133770..31e41e3dd2 100644 --- a/deal.II/lac/source/block_sparse_matrix.cc +++ b/deal.II/lac/source/block_sparse_matrix.cc @@ -1,15 +1,15 @@ -//---------------------------- block_block_sparse_matrix.cc --------------------------- +//---------------------------- block_sparse_matrix.cc --------------------------- // $Id$ // Version: $Name$ // -// Copyright (C) 2000, 2001, 2002 by the deal authors +// Copyright (C) 2000, 2001, 2002, 2004 by the deal authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer // to the file deal.II/doc/license.html for the text and // further information on this license. // -//---------------------------- block_block_sparse_matrix.cc --------------------------- +//---------------------------- block_sparse_matrix.cc --------------------------- #include diff --git a/deal.II/lac/source/block_vector.cc b/deal.II/lac/source/block_vector.cc index 54ae82e39b..1e9433b3fb 100644 --- a/deal.II/lac/source/block_vector.cc +++ b/deal.II/lac/source/block_vector.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2000, 2001, 2002, 2003 by the deal authors +// Copyright (C) 2000, 2001, 2002, 2003, 2004 by the deal authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -15,42 +15,27 @@ // explicit instantiations template class BlockVector; -template BlockVector& BlockVector::DEAL_II_MEMBER_OP_TEMPLATE_INST operator=( - const BlockVector&); +template class BlockVector; + template void BlockVector::reinit(const BlockVector&, const bool); template void BlockVector::reinit(const BlockVector&, const bool); -template void BlockVector::equ(const double, - const BlockVector&); -template void BlockVector::equ(const double, - const BlockVector&); -template void BlockVector::scale(const BlockVector&); -template void BlockVector::scale(const BlockVector&); -template class BlockVector; -template BlockVector& BlockVector::DEAL_II_MEMBER_OP_TEMPLATE_INST operator=( - const BlockVector&); template void BlockVector::reinit(const BlockVector&, const bool); template void BlockVector::reinit(const BlockVector&, const bool); -template void BlockVector::equ(const float, - const BlockVector&); -template void BlockVector::equ(const float, - const BlockVector&); -template void BlockVector::scale(const BlockVector&); -template void BlockVector::scale(const BlockVector&); namespace internal { namespace BlockVectorIterators { - template class Iterator; - template class Iterator; + template class Iterator,false>; + template class Iterator,true>; - template class Iterator; - template class Iterator; + template class Iterator,false>; + template class Iterator,true>; } } diff --git a/deal.II/lac/source/petsc_block_sparse_matrix.cc b/deal.II/lac/source/petsc_block_sparse_matrix.cc new file mode 100644 index 0000000000..d2354b2297 --- /dev/null +++ b/deal.II/lac/source/petsc_block_sparse_matrix.cc @@ -0,0 +1,97 @@ +//---------------------------- petsc_block_sparse_matrix.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004 by the deal authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- petsc_block_sparse_matrix.cc --------------------------- + +#include + + + + +namespace PETScWrappers +{ + BlockSparseMatrix::BlockSparseMatrix () + {} + + + + BlockSparseMatrix::~BlockSparseMatrix () + {} + + + + BlockSparseMatrix & + BlockSparseMatrix::operator = (const BlockSparseMatrix &m) + { + BaseClass::operator = (m); + + return *this; + } + + + + BlockSparseMatrix & + BlockSparseMatrix::operator = (const double d) + { + for (unsigned int r=0; rn_block_rows(); ++r) + for (unsigned int c=0; cn_block_cols(); ++c) + this->block(r,c) = d; + + return *this; + } + + + void + BlockSparseMatrix:: + reinit (const unsigned int n_block_rows, + const unsigned int n_block_columns) + { + // first delete previous content of + // the subobjects array + for (unsigned int r=0; rn_block_rows(); ++r) + for (unsigned int c=0; cn_block_cols(); ++c) + { + BlockType *p = this->sub_objects[r][c]; + this->sub_objects[r][c] = 0; + delete p; + } + + this->sub_objects.reinit (0,0); + + // then resize. set sizes of blocks to + // zero. user will later have to call + // collect_sizes for this + this->sub_objects.reinit (n_block_rows, + n_block_columns); + this->row_block_indices.reinit (n_block_rows, 0); + this->column_block_indices.reinit (n_block_columns, 0); + + // and reinitialize the blocks + for (unsigned int r=0; rn_block_rows(); ++r) + for (unsigned int c=0; cn_block_cols(); ++c) + { + BlockType *p = new BlockType(); + this->sub_objects[r][c] = p; + } + } + + + + void + BlockSparseMatrix::collect_sizes () + { + BaseClass::collect_sizes (); + } + +} + + + diff --git a/deal.II/lac/source/petsc_matrix_base.cc b/deal.II/lac/source/petsc_matrix_base.cc index 6f6ae79eb4..d8fb80b9c3 100644 --- a/deal.II/lac/source/petsc_matrix_base.cc +++ b/deal.II/lac/source/petsc_matrix_base.cc @@ -238,6 +238,31 @@ namespace PETScWrappers + unsigned int + MatrixBase::local_size () const + { + int n_rows, n_cols; + int ierr = MatGetLocalSize (matrix, &n_rows, &n_cols); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + return n_rows; + } + + + + std::pair + MatrixBase::local_range () const + { + int begin, end; + const int ierr = MatGetOwnershipRange (static_cast(matrix), + &begin, &end); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + return std::make_pair (begin, end); + } + + + unsigned int MatrixBase::n_nonzero_elements () const { diff --git a/deal.II/lac/source/petsc_parallel_block_sparse_matrix.cc b/deal.II/lac/source/petsc_parallel_block_sparse_matrix.cc new file mode 100644 index 0000000000..2f7113cf92 --- /dev/null +++ b/deal.II/lac/source/petsc_parallel_block_sparse_matrix.cc @@ -0,0 +1,109 @@ +//---------------------------- petsc_parallel_block_sparse_matrix.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004 by the deal authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- petsc_parallel_block_sparse_matrix.cc --------------------------- + +#include + + + + +namespace PETScWrappers +{ + namespace MPI + { + + BlockSparseMatrix::BlockSparseMatrix () + {} + + + + BlockSparseMatrix::~BlockSparseMatrix () + {} + + + + BlockSparseMatrix & + BlockSparseMatrix::operator = (const BlockSparseMatrix &m) + { + BaseClass::operator = (m); + + return *this; + } + + + + BlockSparseMatrix & + BlockSparseMatrix::operator = (const PetscScalar d) + { + for (unsigned int r=0; rn_block_rows(); ++r) + for (unsigned int c=0; cn_block_cols(); ++c) + this->block(r,c) = d; + + return *this; + } + + + void + BlockSparseMatrix:: + reinit (const unsigned int n_block_rows, + const unsigned int n_block_columns) + { + // first delete previous content of + // the subobjects array + for (unsigned int r=0; rn_block_rows(); ++r) + for (unsigned int c=0; cn_block_cols(); ++c) + { + BlockType *p = this->sub_objects[r][c]; + this->sub_objects[r][c] = 0; + delete p; + } + + this->sub_objects.reinit (0,0); + + // then resize. set sizes of blocks to + // zero. user will later have to call + // collect_sizes for this + this->sub_objects.reinit (n_block_rows, + n_block_columns); + this->row_block_indices.reinit (n_block_rows, 0); + this->column_block_indices.reinit (n_block_columns, 0); + + // and reinitialize the blocks + for (unsigned int r=0; rn_block_rows(); ++r) + for (unsigned int c=0; cn_block_cols(); ++c) + { + BlockType *p = new BlockType(); + this->sub_objects[r][c] = p; + } + } + + + + void + BlockSparseMatrix::collect_sizes () + { + BaseClass::collect_sizes (); + } + + + + const MPI_Comm & + BlockSparseMatrix::get_mpi_communicator () const + { + return block(0,0).get_mpi_communicator(); + } + + } +} + + + diff --git a/deal.II/lac/source/petsc_parallel_block_vector.cc b/deal.II/lac/source/petsc_parallel_block_vector.cc new file mode 100644 index 0000000000..8817d718b9 --- /dev/null +++ b/deal.II/lac/source/petsc_parallel_block_vector.cc @@ -0,0 +1,45 @@ +//---------------------------- petsc_parallel_block_vector.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004 by the deal authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- petsc_parallel_block_vector.cc --------------------------- + + +#include +#include + + +#ifdef DEAL_II_USE_PETSC + + +namespace PETScWrappers +{ + namespace MPI + { + BlockVector & + BlockVector::operator = (const PETScWrappers::BlockVector &v) + { + Assert (v.get_block_indices() == this->get_block_indices(), + ExcNonMatchingBlockVectors()); + + for (unsigned int i=0; i &row_lengths, const bool is_symmetric) : communicator (communicator) { - do_reinit (m, n, local_rows, row_lengths, is_symmetric); + do_reinit (m, n, local_rows, local_columns, + row_lengths, is_symmetric); } + SparseMatrix:: + SparseMatrix (const MPI_Comm &communicator, + const CompressedSparsityPattern &sparsity_pattern, + const unsigned int local_rows, + const unsigned int local_columns, + const bool preset_nonzero_locations) + : + communicator (communicator) + { + do_reinit (sparsity_pattern, local_rows, local_columns, + preset_nonzero_locations); + } + + + SparseMatrix & SparseMatrix::operator = (const double d) { @@ -79,6 +98,7 @@ namespace PETScWrappers const unsigned int m, const unsigned int n, const unsigned int local_rows, + const unsigned int local_columns, const unsigned int n_nonzero_per_row, const bool is_symmetric) { @@ -89,7 +109,8 @@ namespace PETScWrappers const int ierr = MatDestroy (matrix); AssertThrow (ierr == 0, ExcPETScError(ierr)); - do_reinit (m, n, local_rows, n_nonzero_per_row, is_symmetric); + do_reinit (m, n, local_rows, local_columns, + n_nonzero_per_row, is_symmetric); } @@ -99,6 +120,7 @@ namespace PETScWrappers const unsigned int m, const unsigned int n, const unsigned int local_rows, + const unsigned int local_columns, const std::vector &row_lengths, const bool is_symmetric) { @@ -109,15 +131,37 @@ namespace PETScWrappers const int ierr = MatDestroy (matrix); AssertThrow (ierr == 0, ExcPETScError(ierr)); - do_reinit (m, n, local_rows, row_lengths, is_symmetric); + do_reinit (m, n, local_rows, local_columns, row_lengths, is_symmetric); } + void + SparseMatrix:: + reinit (const MPI_Comm &communicator, + const CompressedSparsityPattern &sparsity_pattern, + const unsigned int local_rows, + const unsigned int local_columns, + const bool preset_nonzero_locations) + { + this->communicator = communicator; + + // get rid of old matrix and generate a + // new one + const int ierr = MatDestroy (matrix); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + do_reinit (sparsity_pattern, local_rows, local_columns, + preset_nonzero_locations); + } + + + void SparseMatrix::do_reinit (const unsigned int m, const unsigned int n, const unsigned int local_rows, + const unsigned int local_columns, const unsigned int n_nonzero_per_row, const bool is_symmetric) { @@ -126,16 +170,10 @@ namespace PETScWrappers // use the call sequence indicating only // a maximal number of elements per row // for all rows globally - // - // note that we partition the columns - // in the same way as we partition the - // rows, i.e. we assume that users - // specify the same number for - // local_rows as they do for local_size - // on parallel vectors +//TODO[WB]: We should do better by providing ways to tell PETSc how to partition the columns const int ierr = MatCreateMPIAIJ(communicator, - local_rows, std::min(local_rows, n), + local_rows, local_columns, m, n, n_nonzero_per_row, 0, 0, 0, &matrix); @@ -156,6 +194,7 @@ namespace PETScWrappers SparseMatrix::do_reinit (const unsigned int m, const unsigned int n, const unsigned int local_rows, + const unsigned int local_columns, const std::vector &row_lengths, const bool is_symmetric) { @@ -175,16 +214,11 @@ namespace PETScWrappers const std::vector int_row_lengths (row_lengths.begin(), row_lengths.end()); - // note that we partition the columns - // in the same way as we partition the - // rows, i.e. we assume that users - // specify the same number for - // local_rows as they do for local_size - // on parallel vectors +//TODO[WB]: We should do better by providing ways to tell PETSc how to partition the columns const int ierr = MatCreateMPIAIJ(communicator, - local_rows, local_rows, - m, std::min(local_rows, n), + local_rows, local_columns, + m, n, 0, &int_row_lengths[0], 0, 0, &matrix); AssertThrow (ierr == 0, ExcPETScError(ierr)); @@ -197,6 +231,58 @@ namespace PETScWrappers AssertThrow (ierr == 0, ExcPETScError(ierr)); } } + + + + void + SparseMatrix::do_reinit (const CompressedSparsityPattern &sparsity_pattern, + const unsigned int local_rows, + const unsigned int local_columns, + const bool preset_nonzero_locations) + { + std::vector row_lengths (sparsity_pattern.n_rows()); + for (unsigned int i=0; i row_entries; + std::vector row_values; + for (unsigned int i=0; iksp); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + // set the matrices involved. the + // last argument is irrelevant here, + // since we use the solver only once + // anyway + ierr = KSPSetOperators (solver_data->ksp, A, preconditioner, + SAME_PRECONDITIONER); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + // let derived classes set the solver + // type, and the preconditioning + // object set the type of + // preconditioner + set_solver_type (solver_data->ksp); - PC pc; - ierr = KSPGetPC (ksp, &pc); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - preconditioner.set_preconditioner_type (pc); - - // in the deal.II solvers, we always - // honor the initial guess in the - // solution vector. do so here as well: - KSPSetInitialGuessNonzero (ksp, PETSC_TRUE); - - // then a convergence monitor - // function. that function simply checks - // with the solver_control object we have - // in this object for convergence - KSPSetConvergenceTest (ksp, &convergence_test, - reinterpret_cast(&solver_control)); + ierr = KSPGetPC (solver_data->ksp, &solver_data->pc); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + preconditioner.set_preconditioner_type (solver_data->pc); + + // in the deal.II solvers, we always + // honor the initial guess in the + // solution vector. do so here as + // well: + KSPSetInitialGuessNonzero (solver_data->ksp, PETSC_TRUE); + + // then a convergence monitor + // function. that function simply + // checks with the solver_control + // object we have in this object for + // convergence + KSPSetConvergenceTest (solver_data->ksp, &convergence_test, + reinterpret_cast(&solver_control)); + } // then do the real work: set up solver // internal data and solve the // system. - ierr = KSPSetRhs(ksp,b); + ierr = KSPSetRhs(solver_data->ksp,b); AssertThrow (ierr == 0, ExcPETScError(ierr)); - ierr = KSPSetSolution(ksp,x); + ierr = KSPSetSolution(solver_data->ksp, x); AssertThrow (ierr == 0, ExcPETScError(ierr)); - ierr = KSPSetUp (ksp); + ierr = KSPSetUp (solver_data->ksp); AssertThrow (ierr == 0, ExcPETScError(ierr)); - ierr = KSPSolve (ksp); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - // and destroy the solver object - ierr = KSPDestroy (ksp); + ierr = KSPSolve (solver_data->ksp); AssertThrow (ierr == 0, ExcPETScError(ierr)); #endif + // destroy solver object + solver_data.reset (); + // in case of failure: throw // exception if (solver_control.last_check() != SolverControl::success) throw SolverControl::NoConvergence (solver_control.last_step(), solver_control.last_value()); - // otherwise exit as normal + // otherwise exit as normal } @@ -244,7 +263,7 @@ namespace PETScWrappers SolverRichardson::SolverRichardson (SolverControl &cn, - MPI_Comm &mpi_communicator, + const MPI_Comm &mpi_communicator, const AdditionalData &data) : SolverBase (cn, mpi_communicator), @@ -272,7 +291,7 @@ namespace PETScWrappers /* ---------------------- SolverChebychev ------------------------ */ SolverChebychev::SolverChebychev (SolverControl &cn, - MPI_Comm &mpi_communicator, + const MPI_Comm &mpi_communicator, const AdditionalData &data) : SolverBase (cn, mpi_communicator), @@ -296,7 +315,7 @@ namespace PETScWrappers /* ---------------------- SolverCG ------------------------ */ SolverCG::SolverCG (SolverControl &cn, - MPI_Comm &mpi_communicator, + const MPI_Comm &mpi_communicator, const AdditionalData &data) : SolverBase (cn, mpi_communicator), @@ -320,7 +339,7 @@ namespace PETScWrappers /* ---------------------- SolverBiCG ------------------------ */ SolverBiCG::SolverBiCG (SolverControl &cn, - MPI_Comm &mpi_communicator, + const MPI_Comm &mpi_communicator, const AdditionalData &data) : SolverBase (cn, mpi_communicator), @@ -352,7 +371,7 @@ namespace PETScWrappers SolverGMRES::SolverGMRES (SolverControl &cn, - MPI_Comm &mpi_communicator, + const MPI_Comm &mpi_communicator, const AdditionalData &data) : SolverBase (cn, mpi_communicator), @@ -405,7 +424,7 @@ namespace PETScWrappers /* ---------------------- SolverBicgstab ------------------------ */ SolverBicgstab::SolverBicgstab (SolverControl &cn, - MPI_Comm &mpi_communicator, + const MPI_Comm &mpi_communicator, const AdditionalData &data) : SolverBase (cn, mpi_communicator), @@ -429,7 +448,7 @@ namespace PETScWrappers /* ---------------------- SolverCGS ------------------------ */ SolverCGS::SolverCGS (SolverControl &cn, - MPI_Comm &mpi_communicator, + const MPI_Comm &mpi_communicator, const AdditionalData &data) : SolverBase (cn, mpi_communicator), @@ -453,7 +472,7 @@ namespace PETScWrappers /* ---------------------- SolverTFQMR ------------------------ */ SolverTFQMR::SolverTFQMR (SolverControl &cn, - MPI_Comm &mpi_communicator, + const MPI_Comm &mpi_communicator, const AdditionalData &data) : SolverBase (cn, mpi_communicator), @@ -477,7 +496,7 @@ namespace PETScWrappers /* ---------------------- SolverTCQMR ------------------------ */ SolverTCQMR::SolverTCQMR (SolverControl &cn, - MPI_Comm &mpi_communicator, + const MPI_Comm &mpi_communicator, const AdditionalData &data) : SolverBase (cn, mpi_communicator), @@ -501,7 +520,7 @@ namespace PETScWrappers /* ---------------------- SolverCR ------------------------ */ SolverCR::SolverCR (SolverControl &cn, - MPI_Comm &mpi_communicator, + const MPI_Comm &mpi_communicator, const AdditionalData &data) : SolverBase (cn, mpi_communicator), @@ -525,7 +544,7 @@ namespace PETScWrappers /* ---------------------- SolverLSQR ------------------------ */ SolverLSQR::SolverLSQR (SolverControl &cn, - MPI_Comm &mpi_communicator, + const MPI_Comm &mpi_communicator, const AdditionalData &data) : SolverBase (cn, mpi_communicator), diff --git a/deal.II/lac/source/petsc_sparse_matrix.cc b/deal.II/lac/source/petsc_sparse_matrix.cc index 9f8de5208e..701b7bdb84 100644 --- a/deal.II/lac/source/petsc_sparse_matrix.cc +++ b/deal.II/lac/source/petsc_sparse_matrix.cc @@ -51,6 +51,15 @@ namespace PETScWrappers + SparseMatrix:: + SparseMatrix (const CompressedSparsityPattern &sparsity_pattern, + const bool preset_nonzero_locations) + { + do_reinit (sparsity_pattern, preset_nonzero_locations); + } + + + SparseMatrix & SparseMatrix::operator = (const double d) { @@ -91,7 +100,21 @@ namespace PETScWrappers } + void + SparseMatrix:: + reinit (const CompressedSparsityPattern &sparsity_pattern, + const bool preset_nonzero_locations) + { + // get rid of old matrix and generate a + // new one + const int ierr = MatDestroy (matrix); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + do_reinit (sparsity_pattern, preset_nonzero_locations); + } + + + void SparseMatrix::do_reinit (const unsigned int m, const unsigned int n, @@ -149,6 +172,53 @@ namespace PETScWrappers AssertThrow (ierr == 0, ExcPETScError(ierr)); } } + + + + void + SparseMatrix::do_reinit (const CompressedSparsityPattern &sparsity_pattern, + const bool preset_nonzero_locations) + { + std::vector row_lengths (sparsity_pattern.n_rows()); + for (unsigned int i=0; i row_entries; + std::vector row_values; + for (unsigned int i=0; i +#include +#include + + +void test () +{ + BlockVector v(2,1); + v(0) = 1; + v(1) = 2; + + // first check reading through a const + // iterator + { + BlockVector::const_iterator i=v.begin(); + Assert (*i == 1, ExcInternalError()); + ++i; + Assert (*i == 2, ExcInternalError()); + } + + // same, but create iterator in a different + // way + { + BlockVector::const_iterator + i=const_cast&>(v).begin(); + Assert (*i == 1, ExcInternalError()); + ++i; + Assert (*i == 2, ExcInternalError()); + } + + // read through a read-write iterator + { + BlockVector::iterator i = v.begin(); + Assert (*i == 1, ExcInternalError()); + ++i; + Assert (*i == 2, ExcInternalError()); + } + + // write through a read-write iterator + { + BlockVector::iterator i = v.begin(); + + *i = 2; + ++i; + *i = 3; + } + + // and read again + { + BlockVector::iterator i = v.begin(); + Assert (*i == 2, ExcInternalError()); + ++i; + Assert (*i == 3, ExcInternalError()); + } + + deallog << "OK" << std::endl; +} + + + +int main () +{ + std::ofstream logfile("block_vector_iterator_01.output"); + deallog.attach(logfile); + deallog.depth_console(0); + + try + { + test (); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/bits/block_vector_iterator_02.cc b/tests/bits/block_vector_iterator_02.cc new file mode 100644 index 0000000000..b0ff564a2d --- /dev/null +++ b/tests/bits/block_vector_iterator_02.cc @@ -0,0 +1,106 @@ +//---------------------------- block_vector_iterator_02.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- block_vector_iterator_02.cc --------------------------- + + +// like _01, except that we use operator[] instead of operator* + +#include "../tests.h" +#include +#include +#include + + +void test () +{ + BlockVector v(2,1); + v(0) = 1; + v(1) = 2; + + // first check reading through a const + // iterator + { + BlockVector::const_iterator i=v.begin(); + Assert (i[0] == 1, ExcInternalError()); + Assert (i[1] == 2, ExcInternalError()); + } + + // same, but create iterator in a different + // way + { + BlockVector::const_iterator + i=const_cast&>(v).begin(); + Assert (i[0] == 1, ExcInternalError()); + Assert (i[1] == 2, ExcInternalError()); + } + + // read through a read-write iterator + { + BlockVector::iterator i = v.begin(); + Assert (i[0] == 1, ExcInternalError()); + Assert (i[1] == 2, ExcInternalError()); + } + + // write through a read-write iterator + { + BlockVector::iterator i = v.begin(); + i[0] = 2; + i[1] = 3; + } + + // and read again + { + BlockVector::iterator i = v.begin(); + Assert (i[0] == 2, ExcInternalError()); + Assert (i[1] == 3, ExcInternalError()); + } + + deallog << "OK" << std::endl; +} + + + +int main () +{ + std::ofstream logfile("block_vector_iterator_02.output"); + deallog.attach(logfile); + deallog.depth_console(0); + + try + { + test (); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/bits/metis_02.cc b/tests/bits/metis_02.cc new file mode 100644 index 0000000000..40e8f07f89 --- /dev/null +++ b/tests/bits/metis_02.cc @@ -0,0 +1,122 @@ +//---------------------------- metis_02.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- metis_02.cc --------------------------- + + +// check GridTools::partition_triangulation. generate some output + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +template +void test () +{ + deallog << dim << "D" << std::endl; + + Triangulation triangulation; + GridGenerator::hyper_cube (triangulation); + triangulation.refine_global (4-dim); + for (unsigned int i=0; i<11-2*dim; ++i) + { + typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + for (unsigned int index=0; cell != triangulation.end(); ++cell, ++index) + if (index % (3*dim) == 0) + cell->set_refine_flag(); + triangulation.execute_coarsening_and_refinement (); + } + + + // subdivide into 5 subdomains + GridTools::partition_triangulation (5, triangulation); + + // generate a field where the value equals + // the subdomain number, and output it + FE_Q fe(2); + DoFHandler dof_handler (triangulation); + dof_handler.distribute_dofs (fe); + + Vector partitions (triangulation.n_active_cells()); + { + typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + for (unsigned int index=0; cell != triangulation.end(); ++cell, ++index) + partitions(index) = cell->subdomain_id(); + } + Vector dof_part (dof_handler.n_dofs()); + std::vector assoc (dof_handler.n_dofs()); + DoFTools::get_subdomain_association (dof_handler, assoc); + + deallog << "Cell association:" << std::endl; + for (unsigned int i=0; i (); + test<2> (); + test<3> (); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/bits/petsc_block_vector_iterator_01.cc b/tests/bits/petsc_block_vector_iterator_01.cc new file mode 100644 index 0000000000..71567befcf --- /dev/null +++ b/tests/bits/petsc_block_vector_iterator_01.cc @@ -0,0 +1,116 @@ +//---------------------------- petsc_block_vector_iterator_01.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- petsc_block_vector_iterator_01.cc --------------------------- + + +// make sure that block vector iterator allows reading and writing correctly + +#include "../tests.h" +#include +#include +#include + + +void test () +{ + PETScWrappers::BlockVector v(2,1); + v(0) = 1; + v(1) = 2; + + // first check reading through a const + // iterator + { + PETScWrappers::BlockVector::const_iterator i=v.begin(); + Assert (*i == 1, ExcInternalError()); + ++i; + Assert (*i == 2, ExcInternalError()); + } + + // same, but create iterator in a different + // way + { + PETScWrappers::BlockVector::const_iterator + i=const_cast(v).begin(); + Assert (*i == 1, ExcInternalError()); + ++i; + Assert (*i == 2, ExcInternalError()); + } + + // read through a read-write iterator + { + PETScWrappers::BlockVector::iterator i = v.begin(); + Assert (*i == 1, ExcInternalError()); + ++i; + Assert (*i == 2, ExcInternalError()); + } + + // write through a read-write iterator + { + PETScWrappers::BlockVector::iterator i = v.begin(); + + *i = 2; + ++i; + *i = 3; + } + + // and read again + { + PETScWrappers::BlockVector::iterator i = v.begin(); + Assert (*i == 2, ExcInternalError()); + ++i; + Assert (*i == 3, ExcInternalError()); + } + + deallog << "OK" << std::endl; +} + + + +int main (int argc,char **argv) +{ + std::ofstream logfile("petsc_block_vector_iterator_01.output"); + deallog.attach(logfile); + deallog.depth_console(0); + + try + { + PetscInitialize(&argc,&argv,0,0); + { + test (); + } + PetscFinalize(); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/bits/petsc_block_vector_iterator_02.cc b/tests/bits/petsc_block_vector_iterator_02.cc new file mode 100644 index 0000000000..d395838c3b --- /dev/null +++ b/tests/bits/petsc_block_vector_iterator_02.cc @@ -0,0 +1,110 @@ +//---------------------------- petsc_block_vector_iterator_02.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- petsc_block_vector_iterator_02.cc --------------------------- + + +// like _01, except that we use operator[] instead of operator* + +#include "../tests.h" +#include +#include +#include + + +void test () +{ + PETScWrappers::BlockVector v(2,1); + v(0) = 1; + v(1) = 2; + + // first check reading through a const + // iterator + { + PETScWrappers::BlockVector::const_iterator i=v.begin(); + Assert (i[0] == 1, ExcInternalError()); + Assert (i[1] == 2, ExcInternalError()); + } + + // same, but create iterator in a different + // way + { + PETScWrappers::BlockVector::const_iterator + i=const_cast(v).begin(); + Assert (i[0] == 1, ExcInternalError()); + Assert (i[1] == 2, ExcInternalError()); + } + + // read through a read-write iterator + { + PETScWrappers::BlockVector::iterator i = v.begin(); + Assert (i[0] == 1, ExcInternalError()); + Assert (i[1] == 2, ExcInternalError()); + } + + // write through a read-write iterator + { + PETScWrappers::BlockVector::iterator i = v.begin(); + i[0] = 2; + i[1] = 3; + } + + // and read again + { + PETScWrappers::BlockVector::iterator i = v.begin(); + Assert (i[0] == 2, ExcInternalError()); + Assert (i[1] == 3, ExcInternalError()); + } + + deallog << "OK" << std::endl; +} + + + +int main (int argc,char **argv) +{ + std::ofstream logfile("petsc_block_vector_iterator_02.output"); + deallog.attach(logfile); + deallog.depth_console(0); + + try + { + PetscInitialize(&argc,&argv,0,0); + { + test (); + } + PetscFinalize(); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/bits/petsc_block_vector_iterator_03.cc b/tests/bits/petsc_block_vector_iterator_03.cc new file mode 100644 index 0000000000..5b07026229 --- /dev/null +++ b/tests/bits/petsc_block_vector_iterator_03.cc @@ -0,0 +1,365 @@ +//-------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2000, 2001, 2002, 2003, 2004 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//-------------------------------------------------------------------- + +// this test is an adaptation of lac/block_vector_iterator for PETSc block +// vectors + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include + +template +bool operator == (const PETScWrappers::BlockVector &v1, + const PETScWrappers::BlockVector &v2) +{ + if (v1.size() != v2.size()) + return false; + for (unsigned int i=0; i ivector(4); + ivector[0] = 2; + ivector[1] = 4; + ivector[2] = 3; + ivector[3] = 5; + + // Check 1: initialization via + // iterators + if (true) + { + PETScWrappers::BlockVector v1(ivector); + PETScWrappers::BlockVector v2(ivector); + + // initialize first vector with + // simple loop + for (unsigned int i=0; i(v1.size()) ? + "true" : "false") + << std::endl; + + // check std::copy + PETScWrappers::BlockVector v2(ivector); + std::copy (v1.begin(), v1.end(), v2.begin()); + deallog << "Check 6: " << (v1 == v2 ? "true" : "false") << std::endl; + + // check std::transform + std::transform (v1.begin(), v1.end(), v2.begin(), + std::bind2nd (std::multiplies(), + 2.0)); + v2 *= 1./2.; + deallog << "Check 7: " << (v1 == v2 ? "true" : "false") << std::endl; + + + // check operators +/-, +=/-= + deallog << "Check 8: " + << (std::distance(v1.begin(), v1.begin()+3) == 3 ? + "true" : "false") + << std::endl; + deallog << "Check 9: " + << (std::distance(v1.end()-6, v1.end()) == 6 ? + "true" : "false") + << std::endl; + deallog << "Check 10: " + << (std::distance(v1.begin(), v1.end()) == (signed)v1.size() ? + "true" : "false") + << std::endl; + deallog << "Check 11: " + << (std::distance(v1.begin(), (v1.begin()+=7)) == 7 ? + "true" : "false") + << std::endl; + deallog << "Check 12: " + << (std::distance((v1.end()-=4), v1.end()) == 4 ? + "true" : "false") + << std::endl; + + // check advance + PETScWrappers::BlockVector::iterator p2 = v1.begin(); + std::advance (p2, v1.size()); + deallog << "Check 13: " << (p2 == v1.end() ? "true" : "false") << std::endl; + + PETScWrappers::BlockVector::const_iterator p3 = v1.begin(); + std::advance (p3, v1.size()); + deallog << "Check 14: " << (p3 == v1.end() ? "true" : "false") << std::endl; + }; + + // Check 15: initialization through + // iterators + if (true) + { + PETScWrappers::BlockVector v1(ivector); + // initialize first vector with + // simple loop + for (unsigned int i=0; i v2(v1.begin(), v1.end()); + + // and reverse way + PETScWrappers::BlockVector v3(ivector, v2.begin(), v2.end()); + deallog << "Check 15: " << (v1==v3 ? "true" : "false") << std::endl; + }; + + // Check 16: initialization through + // iterators. variant with one + // constant object + if (true) + { + PETScWrappers::BlockVector v1(ivector); + // initialize first vector with + // simple loop + for (unsigned int i=0; i v2(v1.begin(), v1.end()); + + // and reverse way + PETScWrappers::BlockVector v3(ivector, v2.begin(), v2.end()); + deallog << "Check 16: " << (v1==v3 ? "true" : "false") << std::endl; + }; + + // Check 17: initialization through + // iterators. variant with two + // constant object + if (true) + { + PETScWrappers::BlockVector v1(ivector); + // initialize first vector with + // simple loop + for (unsigned int i=0; i v2(v1.begin(), v1.end()); + + // and reverse way + const PETScWrappers::BlockVector v3(ivector, v2.begin(), v2.end()); + deallog << "Check 17: " << (v1==v3 ? "true" : "false") << std::endl; + }; + + // Check 18: initialization through + // iterators. variant with three + // constant object + if (true) + { + PETScWrappers::BlockVector v0(ivector); + // initialize first vector with + // simple loop + for (unsigned int i=0; i v2(v1.begin(), v1.end()); + + // and reverse way + const PETScWrappers::BlockVector v3(ivector, v2.begin(), v2.end()); + deallog << "Check 18: " << (v1==v3 ? "true" : "false") << std::endl; + }; + + // Check 19: operator[] + if (true) + { + PETScWrappers::BlockVector v1(ivector); + for (unsigned int i=0; i