--- /dev/null
+Removed: The deprecated classes BlockMatrixArray (and derived classes) and MeanValueFilter
+have been removed. Use BlockLinearOperator instead.
+<br>
+(Daniel Arndt, 2020/03/28)
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2005 - 2018 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE.md at
-// the top level directory of deal.II.
-//
-// ---------------------------------------------------------------------
-
-
-// See documentation of BlockMatrixArray for documentation of this example
-
-#include <deal.II/base/logstream.h>
-#include <deal.II/lac/block_matrix_array.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/block_vector.h>
-#include <deal.II/lac/precondition.h>
-#include <deal.II/lac/solver_cg.h>
-#include <deal.II/lac/solver_gmres.h>
-
-#include <iostream>
-#include <fstream>
-
-using namespace dealii;
-
-double Adata[4][4] = {{4., .5, .1, 0.},
- {.5, 4., .5, .1},
- {.1, .5, 4., .5},
- {0., .1, .5, 4.}};
-
-double B1data[4][2] = {{.5, .1}, {.4, .2}, {.3, .3}, {.2, .4}};
-
-double B2data[2][4] = {{.3, 0., -.3, 0.}, {-.3, 0., .3, 0.}};
-
-double Cdata[2][2] = {{8., 1.}, {1., 8.}};
-
-int main()
-{
- FullMatrix<float> A(4, 4);
- FullMatrix<float> B1(4, 2);
- FullMatrix<float> B2(2, 4);
- FullMatrix<float> C(2, 2);
-
- A.fill(&Adata[0][0]);
- B1.fill(&B1data[0][0]);
- B2.fill(&B2data[0][0]);
- C.fill(&Cdata[0][0]);
-
- BlockMatrixArray<double> matrix(2, 2);
-
- matrix.enter(A, 0, 0, 2.);
- matrix.enter(B1, 0, 1, -1.);
- matrix.enter(B2, 0, 1, 1., true);
- matrix.enter(B2, 1, 0, 1.);
- matrix.enter(B1, 1, 0, -1., true);
- matrix.enter(C, 1, 1);
- matrix.print_latex(deallog);
-
- std::vector<unsigned int> block_sizes(2);
- block_sizes[0] = 4;
- block_sizes[1] = 2;
-
- BlockVector<double> result(block_sizes);
- BlockVector<double> x(block_sizes);
- BlockVector<double> y(block_sizes);
- for (unsigned int i = 0; i < result.size(); ++i)
- result(i) = i;
-
- matrix.vmult(y, result);
-
- SolverControl control(100, 1.e-10);
- PreconditionIdentity id;
-
- SolverCG<BlockVector<double>> cg(control);
- cg.solve(matrix, x, y, id);
- x.add(-1., result);
- deallog << "Error " << x.l2_norm() << std::endl;
-
- deallog << "Error A-norm " << std::sqrt(matrix.matrix_norm_square(x))
- << std::endl;
-
- FullMatrix<float> Ainv(4, 4);
- Ainv.invert(A);
- FullMatrix<float> Cinv(2, 2);
- Cinv.invert(C);
-
- BlockTrianglePrecondition<double> precondition(2);
- precondition.enter(Ainv, 0, 0, .5);
- precondition.enter(Cinv, 1, 1);
-
- cg.solve(matrix, x, y, precondition);
- x.add(-1., result);
- deallog << "Error " << x.l2_norm() << std::endl;
-
- precondition.enter(B1, 1, 0, -1., true);
- precondition.enter(B2, 1, 0, 1.);
-
- SolverGMRES<BlockVector<double>> gmres(control);
- gmres.solve(matrix, x, y, precondition);
- x.add(-1., result);
- deallog << "Error " << x.l2_norm() << std::endl;
-
- return 0;
-}
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2001 - 2019 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE.md at
-// the top level directory of deal.II.
-//
-// ---------------------------------------------------------------------
-
-#ifndef dealii_block_matrix_array_h
-#define dealii_block_matrix_array_h
-
-#include <deal.II/base/config.h>
-
-#include <deal.II/base/subscriptor.h>
-#include <deal.II/base/table.h>
-
-#include <deal.II/lac/block_vector.h>
-#include <deal.II/lac/pointer_matrix.h>
-#include <deal.II/lac/vector_memory.h>
-
-#include <map>
-#include <memory>
-#include <sstream>
-#include <string>
-#include <vector>
-
-DEAL_II_NAMESPACE_OPEN
-
-/*! @addtogroup Matrix2
- *@{
- */
-
-
-/**
- * Block matrix composed of different single matrices; these matrices may even
- * be of different types.
- *
- * Given a set of arbitrary matrices <i>A<sub>i</sub></i>, this class
- * implements a block matrix with block entries of the form <i>M<sub>jk</sub>
- * = s<sub>jk</sub>A<sub>i</sub></i>. Each <i>A<sub>i</sub></i> may be used
- * several times with different prefix. The matrices are not copied into the
- * BlockMatrixArray object, but rather a PointerMatrix referencing each of
- * them will be stored along with factors and transposition flags.
- *
- * Non-zero entries are registered by the function enter(), zero entries are
- * not stored at all. Using enter() with the same location <tt>(i,j)</tt>
- * several times will add the corresponding matrices in matrix-vector
- * multiplications. These matrices will not be actually added, but the
- * multiplications with them will be summed up.
- *
- * @note This mechanism makes it impossible to access single entries of
- * BlockMatrixArray. In particular, (block) relaxation preconditioners based
- * on PreconditionRelaxation or PreconditionBlock <b>cannot</b> be used with
- * this class. If you need a preconditioner for a BlockMatrixArray object, use
- * BlockTrianglePrecondition.
- *
- * <h3>Requirements on MatrixType</h3>
- *
- * The template argument <tt>MatrixType</tt> is a class providing the matrix-
- * vector multiplication functions vmult(), Tvmult(), vmult_add() and
- * Tvmult_add() used in this class, but with arguments of type
- * Vector<number> instead of BlockVector<number>. Every matrix
- * which can be used by PointerMatrix is allowed, in particular SparseMatrix
- * is a possible entry type.
- *
- * <h3>Example program</h3> We document the relevant parts of
- * <tt>examples/doxygen/block_matrix_array.cc</tt>.
- *
- * @dontinclude block_matrix_array.cc
- *
- * Obviously, we have to include the header file containing the definition of
- * BlockMatrixArray:
- * @skipline block_matrix_array.h
- *
- * First, we set up some matrices to be entered into the blocks.
- * @skip main
- * @until C.fill
- *
- * Now, we are ready to build a <i>2x2</i> BlockMatrixArray.
- * @line Block
- * First, we enter the matrix <tt>A</tt> multiplied by 2 in the upper left
- * block
- * @line enter
- * Now -1 times <tt>B1</tt> in the upper right block.
- * @line enter
- * We add the transpose of <tt>B2</tt> to the upper right block and continue
- * in a similar fashion. In the end, the block matrix structure is printed
- * into an LaTeX table.
- * @until latex
- *
- * Now, we set up vectors to be multiplied with this matrix and do a
- * multiplication.
- * @until vmult
- *
- * Finally, we solve a linear system with BlockMatrixArray, using no
- * preconditioning and the conjugate gradient method.
- * @until Error
- *
- * The remaining code of this sample program concerns preconditioning and is
- * described in the documentation of BlockTrianglePrecondition.
- *
- * @deprecated This class has been superseded by BlockLinearOperator.
- *
- * @see
- * @ref GlossBlockLA "Block (linear algebra)"
- * @author Guido Kanschat
- * @date 2000-2005, 2010
- */
-template <typename number = double,
- typename BlockVectorType = BlockVector<number>>
-class DEAL_II_DEPRECATED BlockMatrixArray : public Subscriptor
-{
-public:
- /**
- * Declare the type for container size.
- */
- using size_type = types::global_dof_index;
-
- /**
- * Default constructor creating a useless object. initialize() must be
- * called before using it.
- */
- BlockMatrixArray();
-
- /**
- * Constructor fixing the dimensions.
- */
- BlockMatrixArray(const unsigned int n_block_rows,
- const unsigned int n_block_cols);
-
- /**
- * Initialize object completely. This is the function to call for an object
- * created by the default constructor.
- */
- void
- initialize(const unsigned int n_block_rows, const unsigned int n_block_cols);
-
- /**
- * Adjust the matrix to a new size and delete all blocks.
- */
- void
- reinit(const unsigned int n_block_rows, const unsigned int n_block_cols);
-
- /**
- * Add a block matrix entry. The <tt>matrix</tt> is entered into a list of
- * blocks for multiplication, together with its coordinates <tt>row</tt> and
- * <tt>col</tt> as well as optional multiplication factor <tt>prefix</tt>
- * and transpose flag <tt>transpose</tt>.
- *
- * @note No check for consistency of block sizes is made. Therefore,
- * entering a block of wrong dimension here will only lead to a
- * ExcDimensionMismatch in one of the multiplication functions.
- */
- template <typename MatrixType>
- void
- enter(const MatrixType & matrix,
- const unsigned int row,
- const unsigned int col,
- const number prefix = 1.,
- const bool transpose = false);
-
- /**
- * Delete all entries, i.e. reset the matrix to an empty state.
- */
- void
- clear();
-
- /**
- * Number of block-entries per column.
- */
- unsigned int
- n_block_rows() const;
-
- /**
- * Number of block-entries per row.
- */
- unsigned int
- n_block_cols() const;
-
- /**
- * Matrix-vector multiplication.
- */
- void
- vmult(BlockVectorType &dst, const BlockVectorType &src) const;
-
- /**
- * Matrix-vector multiplication adding to <tt>dst</tt>.
- */
- void
- vmult_add(BlockVectorType &dst, const BlockVectorType &src) const;
-
- /**
- * Transposed matrix-vector multiplication.
- */
- void
- Tvmult(BlockVectorType &dst, const BlockVectorType &src) const;
-
- /**
- * Transposed matrix-vector multiplication adding to <tt>dst</tt>.
- */
- void
- Tvmult_add(BlockVectorType &dst, const BlockVectorType &src) const;
-
- /**
- * Matrix scalar product between two vectors (at least for a symmetric
- * matrix).
- */
- number
- matrix_scalar_product(const BlockVectorType &u,
- const BlockVectorType &v) const;
-
- /**
- * Compute $u^T M u$. This is the square of the norm induced by the matrix
- * assuming the matrix is symmetric positive definitive.
- */
- number
- matrix_norm_square(const BlockVectorType &u) const;
-
- /**
- * Print the block structure as a LaTeX-array. This output will not be very
- * intuitive, since the current object lacks knowledge about what the
- * individual blocks represent or how they should be named. Instead, what
- * you will see is an entry for each block showing all the matrices with
- * their multiplication factors and possibly transpose marks. The matrices
- * itself are named successively as they are encountered. If the same matrix
- * is entered several times, it will be listed with different names.
- *
- * As an example, consider the following code:
- * @code
- * FullMatrix<double> A1(4,4);
- * FullMatrix<double> A2(4,4);
- * FullMatrix<double> B(4,3);
- * FullMatrix<double> C(3,3);
- *
- * BlockMatrixArray<double> block(2,2);
- *
- * block.enter(A1,0,0);
- * block.enter(A2,0,0,2,true);
- * block.enter(B,0,1,-3.);
- * block.enter(B,0,1,-3.,true);
- * block.enter(C,1,1,1.,true);
- *
- * block.print_latex(std::cout);
- * @endcode
- * The current function will then produce output of the following kind:
- * @code
- * \begin{array}{cc}
- * M0+2xM1^T & -3xM2-3xM3^T\\
- * & M4^T
- * \end{array}
- * @endcode
- * Note how the individual blocks here are just numbered successively as
- * <code>M0</code> to <code>M4</code> and that the output misses the fact
- * that <code>M2</code> and <code>M3</code> are, in fact, the same matrix.
- * Nevertheless, the output at least gives some kind of idea of the block
- * structure of this matrix.
- */
- template <class StreamType>
- void
- print_latex(StreamType &out) const;
-
-protected:
- /**
- * Internal data structure.
- *
- * For each entry of a BlockMatrixArray, its position, matrix, prefix and
- * optional transposition must be stored. This structure encapsulates all of
- * them.
- *
- * @author Guido Kanschat, 2000, 2001
- */
- class Entry
- {
- public:
- /**
- * Constructor initializing all data fields. A PointerMatrix object is
- * generated for <tt>matrix</tt>.
- */
- template <typename MatrixType>
- Entry(const MatrixType &matrix,
- size_type row,
- size_type col,
- number prefix,
- bool transpose);
-
- /**
- * Copy constructor invalidating the old object. Since it is only used for
- * entering temporary objects into a vector, this is ok.
- *
- * For a deep copy, we would need a reproduction operator in
- * PointerMatixBase.
- */
- Entry(const Entry &);
-
- /**
- * Destructor, where we delete the PointerMatrix created by the
- * constructor.
- */
- ~Entry();
-
- /**
- * Row number in the block matrix.
- */
- size_type row;
-
- /**
- * Column number in the block matrix.
- */
- size_type col;
-
- /**
- * Factor in front of the matrix block.
- */
- number prefix;
-
- /**
- * Indicates that matrix block must be transposed for multiplication.
- */
- bool transpose;
-
- /**
- * The matrix block itself.
- */
- PointerMatrixBase<typename BlockVectorType::BlockType> *matrix;
-
- /**
- * Assignment operator.
- *
- * @note Since the copy constructor is destructive (see its documentation)
- * and only exists for convenience there is no reasonable way to implement
- * this, so it is explicitly deleted.
- */
- Entry &
- operator=(const Entry &) = delete;
- };
-
- /**
- * Array of block entries in the matrix.
- */
- std::vector<Entry> entries;
-
-private:
- /**
- * Number of blocks per column.
- */
- unsigned int block_rows;
- /**
- * number of blocks per row.
- */
- unsigned int block_cols;
-};
-
-/*@}*/
-
-
-/**
- * Inversion of a block-triangular matrix.
- *
- * In this block matrix, the inverses of the diagonal blocks are stored
- * together with the off-diagonal blocks of a block matrix. Then, forward or
- * backward insertion is performed block-wise. The diagonal blocks are NOT
- * inverted for this purpose!
- *
- * Like for all preconditioners, the preconditioning operation is performed by
- * the vmult() member function.
- *
- * @note While block indices may be duplicated (see BlockMatrixArray) to add
- * blocks, this has to be used with caution, since summing up the inverse of
- * two blocks does not yield the inverse of the sum. While the latter would be
- * desirable, we can only perform the first.
- *
- * The implementation may be a little clumsy, but it should be sufficient as
- * long as the block sizes are much larger than the number of blocks.
- *
- * <h3>Example</h3> Here, we document the second part of
- * <tt>examples/doxygen/block_matrix_array.cc</tt>. For the beginning of this
- * file, see BlockMatrixArray.
- *
- * In order to set up the preconditioner, we have to compute the inverses of
- * the diagonal blocks ourselves. Since we used FullMatrix objects, this is
- * fairly easy.
- * @dontinclude block_matrix_array.cc
- * @skip Error
- * @until Cinv.invert
- *
- * After creating a <i>2x2</i> BlockTrianglePrecondition object, we only fill
- * its diagonals. The scaling factor <i>1/2</i> used for <tt>A</tt> is the
- * reciprocal of the scaling factor used for the <tt>matrix</tt> itself.
- * Remember, this preconditioner actually <b>multiplies</b> with the diagonal
- * blocks.
- * @until Cinv
- *
- * Now, we have a block Jacobi preconditioner, which is still symmetric, since
- * the blocks are symmetric. Therefore, we can still use the preconditioned
- * conjugate gradient method.
- * @until Error
- *
- * Now, we enter the subdiagonal block. This is the same as in
- * <tt>matrix</tt>.
- * @until B2
- *
- * Since the preconditioner is not symmetric anymore, we use the GMRES method
- * for solving.
- * @until Error
- *
- * @deprecated This class has been superseded by block_back_substitution and
- * block_forward_substitution, which build on BlockLinearOperator.
- *
- * @ingroup Preconditioners
- * @author Guido Kanschat, 2001, 2005
- */
-template <typename number = double,
- typename BlockVectorType = BlockVector<number>>
-class DEAL_II_DEPRECATED BlockTrianglePrecondition
- : private BlockMatrixArray<number, BlockVectorType>
-{
-public:
- /**
- * Declare type for container size.
- */
- using size_type = types::global_dof_index;
-
- /**
- * Default constructor creating a useless object. initialize() must be
- * called before using it.
- */
- BlockTrianglePrecondition();
-
- /**
- * Constructor. This matrix must be block-quadratic, and <tt>n_blocks</tt>
- * is the number of blocks in each direction.
- */
- BlockTrianglePrecondition(const unsigned int n_blocks);
-
- /**
- * Resize preconditioner to a new size and clear all blocks.
- */
- void
- reinit(const unsigned int n_block_rows);
-
-
- /**
- * Enter a block. This calls BlockMatrixArray::enter(). Remember that the
- * diagonal blocks should actually be inverse matrices or preconditioners.
- */
- template <typename MatrixType>
- void
- enter(const MatrixType &matrix,
- const size_type row,
- const size_type col,
- const number prefix = 1.,
- const bool transpose = false);
-
- /**
- * Preconditioning.
- */
- void
- vmult(BlockVectorType &dst, const BlockVectorType &src) const;
-
- /**
- * Preconditioning adding to <tt>dst</tt>.
- */
- void
- vmult_add(BlockVectorType &dst, const BlockVectorType &src) const;
-
- /**
- * Transposed preconditioning
- */
- void
- Tvmult(BlockVectorType &dst, const BlockVectorType &src) const;
-
- /**
- * Transposed preconditioning adding to <tt>dst</tt>.
- */
- void
- Tvmult_add(BlockVectorType &dst, const BlockVectorType &src) const;
-
- /**
- * Make function of base class available.
- */
- using BlockMatrixArray<number, BlockVectorType>::print_latex;
-
- /**
- * Make function of base class available.
- */
- using BlockMatrixArray<number, BlockVectorType>::n_block_rows;
-
- /**
- * Make function of base class available.
- */
- using BlockMatrixArray<number, BlockVectorType>::n_block_cols;
- using BlockMatrixArray<number, BlockVectorType>::clear;
- using BlockMatrixArray<number, BlockVectorType>::Subscriptor::subscribe;
- using BlockMatrixArray<number, BlockVectorType>::Subscriptor::unsubscribe;
-
- /**
- * @addtogroup Exceptions
- * @{
- */
-
- /**
- * Each diagonal block must contain one and only one matrix. If this
- * exception is thrown, you did not enter a matrix here.
- */
- DeclException1(ExcNoDiagonal,
- size_type,
- << "No diagonal entry was added for block " << arg1);
-
- /**
- * Each diagonal block must contain one and only one matrix. If this
- * exception is thrown, you entered a second matrix here.
- */
- DeclException1(ExcMultipleDiagonal,
- size_type,
- << "Inverse diagonal entries may not be added in block "
- << arg1);
- //@}
-private:
- /**
- * Add all off-diagonal contributions and return the entry of the diagonal
- * element for one row.
- */
- void
- do_row(BlockVectorType &dst, size_type row_num) const;
-
- /**
- * Flag for backward insertion.
- */
- bool backward;
-};
-
-
-#ifndef DOXYGEN
-//---------------------------------------------------------------------------
-
-template <typename number, typename BlockVectorType>
-template <typename MatrixType>
-inline BlockMatrixArray<number, BlockVectorType>::Entry::Entry(
- const MatrixType &m,
- size_type row,
- size_type col,
- number prefix,
- bool transpose)
- : row(row)
- , col(col)
- , prefix(prefix)
- , transpose(transpose)
- , matrix(new_pointer_matrix_base(m,
- typename BlockVectorType::BlockType(),
- typeid(*this).name()))
-{}
-
-
-
-template <typename number, typename BlockVectorType>
-template <typename MatrixType>
-inline void
-BlockMatrixArray<number, BlockVectorType>::enter(const MatrixType &matrix,
- unsigned int row,
- unsigned int col,
- number prefix,
- bool transpose)
-{
- AssertIndexRange(row, n_block_rows());
- AssertIndexRange(col, n_block_cols());
- entries.push_back(Entry(matrix, row, col, prefix, transpose));
-}
-
-
-template <typename number, typename BlockVectorType>
-template <class StreamType>
-inline void
-BlockMatrixArray<number, BlockVectorType>::print_latex(StreamType &out) const
-{
- out << "\\begin{array}{" << std::string(n_block_cols(), 'c') << "}"
- << std::endl;
-
- Table<2, std::string> array(n_block_rows(), n_block_cols());
-
- using NameMap =
- std::map<const PointerMatrixBase<typename BlockVectorType::BlockType> *,
- std::string>;
- NameMap matrix_names;
-
- typename std::vector<Entry>::const_iterator m = entries.begin();
- typename std::vector<Entry>::const_iterator end = entries.end();
-
- size_type matrix_number = 0;
- for (; m != end; ++m)
- {
- if (matrix_names.find(m->matrix) == matrix_names.end())
- {
- std::pair<typename NameMap::iterator, bool> x = matrix_names.insert(
- std::pair<
- const PointerMatrixBase<typename BlockVectorType::BlockType> *,
- std::string>(m->matrix, std::string("M")));
- std::ostringstream stream;
- stream << matrix_number++;
-
- x.first->second += stream.str();
- }
-
- std::ostringstream stream;
-
- if (!array(m->row, m->col).empty() && m->prefix >= 0)
- stream << "+";
- if (m->prefix != 1.)
- stream << m->prefix << 'x';
- stream << matrix_names.find(m->matrix)->second;
- // stream << '(' << m->matrix << ')';
- if (m->transpose)
- stream << "^T";
-
- array(m->row, m->col) += stream.str();
- }
- for (unsigned int i = 0; i < n_block_rows(); ++i)
- for (unsigned int j = 0; j < n_block_cols(); ++j)
- {
- out << '\t' << array(i, j);
- if (j == n_block_cols() - 1)
- {
- if (i != n_block_rows() - 1)
- out << "\\\\" << std::endl;
- else
- out << std::endl;
- }
- else
- out << " &";
- }
- out << "\\end{array}" << std::endl;
-}
-
-template <typename number, typename BlockVectorType>
-template <typename MatrixType>
-inline void
-BlockTrianglePrecondition<number, BlockVectorType>::enter(
- const MatrixType &matrix,
- size_type row,
- size_type col,
- number prefix,
- bool transpose)
-{
- BlockMatrixArray<number, BlockVectorType>::enter(
- matrix, row, col, prefix, transpose);
-}
-
-
-
-#endif // DOXYGEN
-
-DEAL_II_NAMESPACE_CLOSE
-
-#endif
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2002 - 2018 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE.md at
-// the top level directory of deal.II.
-//
-// ---------------------------------------------------------------------
-
-#ifndef dealii_matrix_lib_h
-#define dealii_matrix_lib_h
-
-#include <deal.II/base/config.h>
-
-#include <deal.II/base/subscriptor.h>
-
-#include <deal.II/lac/solver_richardson.h>
-#include <deal.II/lac/vector_memory.h>
-
-DEAL_II_NAMESPACE_OPEN
-
-// Forward declarations
-#ifndef DOXYGEN
-template <typename number>
-class Vector;
-template <typename number>
-class BlockVector;
-template <typename number>
-class SparseMatrix;
-#endif
-
-/*! @addtogroup Matrix2
- *@{
- */
-
-
-/**
- * Mean value filter. The vmult() functions of this matrix filter out mean
- * values of the vector. If the vector is of type BlockVector, then an
- * additional parameter selects a single component for this operation.
- *
- * In mathematical terms, this class acts as if it was the matrix $I-\frac
- * 1n{\mathbf 1}_n{\mathbf 1}_n^T$ where ${\mathbf 1}_n$ is a vector of size
- * $n$ that has only ones as its entries. Thus, taking the dot product between
- * a vector $\mathbf v$ and $\frac 1n {\mathbf 1}_n$ yields the <i>mean
- * value</i> of the entries of ${\mathbf v}$. Consequently, $ \left[I-\frac
- * 1n{\mathbf 1}_n{\mathbf 1}_n^T\right] \mathbf v = \mathbf v - \left[\frac
- * 1n {\mathbf v} \cdot {\mathbf 1}_n\right]{\mathbf 1}_n$ subtracts from every
- * vector element the mean value of all elements.
- *
- * @deprecated Use a LinearOperator, or a BlockLinearOperator instead. you
- * can construct such a filter by using mean_value_filter, or
- * block_diagonal_operator (with a mean_value_filter block), respectively.
- *
- * @author Guido Kanschat, 2002, 2003
- */
-class DEAL_II_DEPRECATED MeanValueFilter : public Subscriptor
-{
-public:
- /**
- * Declare type for container size.
- */
- using size_type = types::global_dof_index;
-
- /**
- * Constructor, optionally selecting a component.
- */
- MeanValueFilter(const size_type component = numbers::invalid_size_type);
-
- /**
- * Subtract mean value from @p v.
- */
- template <typename number>
- void
- filter(Vector<number> &v) const;
-
- /**
- * Subtract mean value from @p v.
- */
- template <typename number>
- void
- filter(BlockVector<number> &v) const;
-
- /**
- * Return the source vector with subtracted mean value.
- */
- template <typename number>
- void
- vmult(Vector<number> &dst, const Vector<number> &src) const;
-
- /**
- * Add source vector with subtracted mean value to dest.
- */
- template <typename number>
- void
- vmult_add(Vector<number> &dst, const Vector<number> &src) const;
-
- /**
- * Return the source vector with subtracted mean value in selected
- * component.
- */
- template <typename number>
- void
- vmult(BlockVector<number> &dst, const BlockVector<number> &src) const;
-
- /**
- * Add a source to dest, where the mean value in the selected component is
- * subtracted.
- */
- template <typename number>
- void
- vmult_add(BlockVector<number> &dst, const BlockVector<number> &src) const;
-
-
- /**
- * Not implemented.
- */
- template <typename VectorType>
- void
- Tvmult(VectorType &, const VectorType &) const;
-
- /**
- * Not implemented.
- */
- template <typename VectorType>
- void
- Tvmult_add(VectorType &, const VectorType &) const;
-
-private:
- /**
- * Component for filtering block vectors.
- */
- const size_type component;
-};
-
-
-
-/*@}*/
-//---------------------------------------------------------------------------
-
-
-template <typename VectorType>
-inline void
-MeanValueFilter::Tvmult(VectorType &, const VectorType &) const
-{
- Assert(false, ExcNotImplemented());
-}
-
-
-template <typename VectorType>
-inline void
-MeanValueFilter::Tvmult_add(VectorType &, const VectorType &) const
-{
- Assert(false, ExcNotImplemented());
-}
-
-
-DEAL_II_NAMESPACE_CLOSE
-
-#endif
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2002 - 2018 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE.md at
-// the top level directory of deal.II.
-//
-// ---------------------------------------------------------------------
-
-#ifndef dealii_matrix_lib_templates_h
-#define dealii_matrix_lib_templates_h
-
-#include <deal.II/base/config.h>
-
-#include <deal.II/lac/block_vector.h>
-#include <deal.II/lac/matrix_lib.h>
-#include <deal.II/lac/vector.h>
-
-DEAL_II_NAMESPACE_OPEN
-
-template <typename number>
-void
-MeanValueFilter::filter(Vector<number> &v) const
-{
- number mean = v.mean_value();
-
- for (size_type i = 0; i < v.size(); ++i)
- v(i) -= mean;
-}
-
-
-
-template <typename number>
-void
-MeanValueFilter::vmult(Vector<number> &dst, const Vector<number> &src) const
-{
- Assert(dst.size() == src.size(),
- ExcDimensionMismatch(dst.size(), src.size()));
-
- number mean = src.mean_value();
-
- for (size_type i = 0; i < dst.size(); ++i)
- dst(i) = src(i) - mean;
-}
-
-
-
-template <typename number>
-void
-MeanValueFilter::vmult_add(Vector<number> &dst, const Vector<number> &src) const
-{
- Assert(dst.size() == src.size(),
- ExcDimensionMismatch(dst.size(), src.size()));
-
- number mean = src.mean_value();
-
- for (size_type i = 0; i < dst.size(); ++i)
- dst(i) += src(i) - mean;
-}
-
-
-
-template <typename number>
-void
-MeanValueFilter::filter(BlockVector<number> &v) const
-{
- Assert(component != numbers::invalid_unsigned_int, ExcNotInitialized());
-
- for (unsigned int i = 0; i < v.n_blocks(); ++i)
- if (i == component)
- vmult(v.block(i), v.block(i));
-}
-
-
-
-template <typename number>
-void
-MeanValueFilter::vmult(BlockVector<number> & dst,
- const BlockVector<number> &src) const
-{
- Assert(component != numbers::invalid_unsigned_int, ExcNotInitialized());
-
- Assert(dst.n_blocks() == src.n_blocks(),
- ExcDimensionMismatch(dst.n_blocks(), src.n_blocks()));
-
- for (unsigned int i = 0; i < dst.n_blocks(); ++i)
- if (i == component)
- vmult(dst.block(i), src.block(i));
- else
- dst.block(i) = src.block(i);
-}
-
-
-
-template <typename number>
-void
-MeanValueFilter::vmult_add(BlockVector<number> & dst,
- const BlockVector<number> &src) const
-{
- Assert(component != numbers::invalid_unsigned_int, ExcNotInitialized());
-
- Assert(dst.n_blocks() == src.n_blocks(),
- ExcDimensionMismatch(dst.n_blocks(), src.n_blocks()));
-
- for (unsigned int i = 0; i < dst.n_blocks(); ++i)
- if (i == component)
- vmult_add(dst.block(i), src.block(i));
- else
- dst.block(i) += src.block(i);
-}
-
-
-DEAL_II_NAMESPACE_CLOSE
-
-#endif
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/householder.h>
#include <deal.II/lac/linear_operator.h>
-#include <deal.II/lac/matrix_lib.h>
#include <deal.II/multigrid/mg_base.h>
#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/lac/block_matrix_array.h>
#include <deal.II/lac/block_sparsity_pattern.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/sparsity_pattern.h>
SET(_unity_include_src
affine_constraints.cc
- block_matrix_array.cc
block_sparse_matrix.cc
block_sparse_matrix_ez.cc
block_sparsity_pattern.cc
la_vector.cc
la_parallel_vector.cc
la_parallel_block_vector.cc
- matrix_lib.cc
matrix_out.cc
precondition_block.cc
precondition_block_ez.cc
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2005 - 2018 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE.md at
-// the top level directory of deal.II.
-//
-// ---------------------------------------------------------------------
-
-
-#include <deal.II/lac/block_matrix_array.h>
-#include <deal.II/lac/block_vector.h>
-#include <deal.II/lac/petsc_block_vector.h>
-#include <deal.II/lac/trilinos_parallel_block_vector.h>
-#include <deal.II/lac/vector.h>
-
-DEAL_II_NAMESPACE_OPEN
-
-
-template <typename number, typename BlockVectorType>
-BlockMatrixArray<number, BlockVectorType>::Entry::Entry(const Entry &e)
- : row(e.row)
- , col(e.col)
- , prefix(e.prefix)
- , transpose(e.transpose)
- , matrix(e.matrix)
-{
- Entry &e2 = const_cast<Entry &>(e);
- e2.matrix = nullptr;
-}
-
-
-
-template <typename number, typename BlockVectorType>
-BlockMatrixArray<number, BlockVectorType>::Entry::~Entry()
-{
- if (matrix)
- delete matrix;
-}
-
-
-template <typename number, typename BlockVectorType>
-BlockMatrixArray<number, BlockVectorType>::BlockMatrixArray()
- : block_rows(0)
- , block_cols(0)
-{}
-
-
-
-template <typename number, typename BlockVectorType>
-BlockMatrixArray<number, BlockVectorType>::BlockMatrixArray(
- const unsigned int n_block_rows,
- const unsigned int n_block_cols)
- : block_rows(n_block_rows)
- , block_cols(n_block_cols)
-{}
-
-
-template <typename number, typename BlockVectorType>
-void
-BlockMatrixArray<number, BlockVectorType>::initialize(
- const unsigned int n_block_rows,
- const unsigned int n_block_cols)
-{
- block_rows = n_block_rows;
- block_cols = n_block_cols;
-}
-
-
-
-template <typename number, typename BlockVectorType>
-void
-BlockMatrixArray<number, BlockVectorType>::reinit(
- const unsigned int n_block_rows,
- const unsigned int n_block_cols)
-{
- clear();
- block_rows = n_block_rows;
- block_cols = n_block_cols;
-}
-
-
-
-template <typename number, typename BlockVectorType>
-void
-BlockMatrixArray<number, BlockVectorType>::clear()
-{
- entries.clear();
-}
-
-
-template <typename number, typename BlockVectorType>
-void
-BlockMatrixArray<number, BlockVectorType>::vmult_add(
- BlockVectorType & dst,
- const BlockVectorType &src) const
-{
- GrowingVectorMemory<typename BlockVectorType::BlockType> mem;
- Assert(dst.n_blocks() == block_rows,
- ExcDimensionMismatch(dst.n_blocks(), block_rows));
- Assert(src.n_blocks() == block_cols,
- ExcDimensionMismatch(src.n_blocks(), block_cols));
-
- typename VectorMemory<typename BlockVectorType::BlockType>::Pointer p_aux(
- mem);
- typename BlockVectorType::BlockType &aux = *p_aux;
-
- typename std::vector<Entry>::const_iterator m = entries.begin();
- typename std::vector<Entry>::const_iterator end = entries.end();
-
- for (; m != end; ++m)
- {
- aux.reinit(dst.block(m->row));
- if (m->transpose)
- m->matrix->Tvmult(aux, src.block(m->col));
- else
- m->matrix->vmult(aux, src.block(m->col));
- dst.block(m->row).add(m->prefix, aux);
- }
-}
-
-
-
-template <typename number, typename BlockVectorType>
-void
-BlockMatrixArray<number, BlockVectorType>::vmult(
- BlockVectorType & dst,
- const BlockVectorType &src) const
-{
- dst = 0.;
- vmult_add(dst, src);
-}
-
-
-
-template <typename number, typename BlockVectorType>
-void
-BlockMatrixArray<number, BlockVectorType>::Tvmult_add(
- BlockVectorType & dst,
- const BlockVectorType &src) const
-{
- GrowingVectorMemory<typename BlockVectorType::BlockType> mem;
- Assert(dst.n_blocks() == block_cols,
- ExcDimensionMismatch(dst.n_blocks(), block_cols));
- Assert(src.n_blocks() == block_rows,
- ExcDimensionMismatch(src.n_blocks(), block_rows));
-
- typename std::vector<Entry>::const_iterator m = entries.begin();
- typename std::vector<Entry>::const_iterator end = entries.end();
-
- typename VectorMemory<typename BlockVectorType::BlockType>::Pointer p_aux(
- mem);
- typename BlockVectorType::BlockType &aux = *p_aux;
-
- for (; m != end; ++m)
- {
- aux.reinit(dst.block(m->col));
- if (m->transpose)
- m->matrix->vmult(aux, src.block(m->row));
- else
- m->matrix->Tvmult(aux, src.block(m->row));
- dst.block(m->col).add(m->prefix, aux);
- }
-}
-
-
-
-template <typename number, typename BlockVectorType>
-void
-BlockMatrixArray<number, BlockVectorType>::Tvmult(
- BlockVectorType & dst,
- const BlockVectorType &src) const
-{
- dst = 0.;
- Tvmult_add(dst, src);
-}
-
-
-
-template <typename number, typename BlockVectorType>
-number
-BlockMatrixArray<number, BlockVectorType>::matrix_scalar_product(
- const BlockVectorType &u,
- const BlockVectorType &v) const
-{
- GrowingVectorMemory<typename BlockVectorType::BlockType> mem;
- Assert(u.n_blocks() == block_rows,
- ExcDimensionMismatch(u.n_blocks(), block_rows));
- Assert(v.n_blocks() == block_cols,
- ExcDimensionMismatch(v.n_blocks(), block_cols));
-
- typename VectorMemory<typename BlockVectorType::BlockType>::Pointer p_aux(
- mem);
- typename BlockVectorType::BlockType &aux = *p_aux;
-
- typename std::vector<Entry>::const_iterator m;
- typename std::vector<Entry>::const_iterator end = entries.end();
-
- number result = 0.;
-
- for (unsigned int i = 0; i < block_rows; ++i)
- {
- aux.reinit(u.block(i));
- for (m = entries.begin(); m != end; ++m)
- {
- if (m->row != i)
- continue;
- if (m->transpose)
- m->matrix->Tvmult_add(aux, v.block(m->col));
- else
- m->matrix->vmult(aux, v.block(m->col));
- }
- result += u.block(i) * aux;
- }
-
- return result;
-}
-
-
-
-template <typename number, typename BlockVectorType>
-number
-BlockMatrixArray<number, BlockVectorType>::matrix_norm_square(
- const BlockVectorType &u) const
-{
- return matrix_scalar_product(u, u);
-}
-
-
-
-template <typename number, typename BlockVectorType>
-unsigned int
-BlockMatrixArray<number, BlockVectorType>::n_block_rows() const
-{
- return block_rows;
-}
-
-
-
-template <typename number, typename BlockVectorType>
-unsigned int
-BlockMatrixArray<number, BlockVectorType>::n_block_cols() const
-{
- return block_cols;
-}
-
-
-
-//---------------------------------------------------------------------------
-
-template <typename number, typename BlockVectorType>
-BlockTrianglePrecondition<number, BlockVectorType>::BlockTrianglePrecondition()
- : BlockMatrixArray<number, BlockVectorType>()
- , backward(false)
-{}
-
-
-template <typename number, typename BlockVectorType>
-BlockTrianglePrecondition<number, BlockVectorType>::BlockTrianglePrecondition(
- const unsigned int block_rows)
- : BlockMatrixArray<number, BlockVectorType>(block_rows, block_rows)
- , backward(false)
-{}
-
-
-template <typename number, typename BlockVectorType>
-void
-BlockTrianglePrecondition<number, BlockVectorType>::reinit(const unsigned int n)
-{
- BlockMatrixArray<number, BlockVectorType>::reinit(n, n);
-}
-
-
-template <typename number, typename BlockVectorType>
-void
-BlockTrianglePrecondition<number, BlockVectorType>::do_row(
- BlockVectorType &dst,
- size_type row_num) const
-{
- GrowingVectorMemory<typename BlockVectorType::BlockType> mem;
- typename std::vector<
- typename BlockMatrixArray<number, BlockVectorType>::Entry>::const_iterator
- m = this->entries.begin();
- typename std::vector<
- typename BlockMatrixArray<number, BlockVectorType>::Entry>::const_iterator
- end = this->entries.end();
- std::vector<typename std::vector<
- typename BlockMatrixArray<number, BlockVectorType>::Entry>::const_iterator>
- diagonals;
-
- typename VectorMemory<typename BlockVectorType::BlockType>::Pointer p_aux(
- mem);
- typename BlockVectorType::BlockType &aux = *p_aux;
-
- aux.reinit(dst.block(row_num), true);
-
- // Loop over all entries, since
- // they are not ordered by rows.
- for (; m != end; ++m)
- {
- const size_type i = m->row;
- // Ignore everything not in
- // this row
- if (i != row_num)
- continue;
- const size_type j = m->col;
- // Only use the lower (upper)
- // triangle for forward
- // (backward) substitution
- if (((j > i) && !backward) || ((j < i) && backward))
- continue;
- if (j == i)
- {
- diagonals.push_back(m);
- }
- else
- {
- if (m->transpose)
- m->matrix->Tvmult(aux, dst.block(j));
- else
- m->matrix->vmult(aux, dst.block(j));
- dst.block(i).add(-1.0 * m->prefix, aux);
- }
- }
- Assert(diagonals.size() != 0, ExcNoDiagonal(row_num));
-
- // Inverting the diagonal block is
- // simple, if there is only one
- // matrix
- if (diagonals.size() == 1)
- {
- if (diagonals[0]->transpose)
- diagonals[0]->matrix->Tvmult(aux, dst.block(row_num));
- else
- diagonals[0]->matrix->vmult(aux, dst.block(row_num));
- dst.block(row_num).equ(diagonals[0]->prefix, aux);
- }
- else
- {
- aux = 0.;
- for (size_type i = 0; i < diagonals.size(); ++i)
- {
- m = diagonals[i];
- // First, divide by the current
- // factor, such that we can
- // multiply by it later.
- aux /= m->prefix;
- if (m->transpose)
- m->matrix->Tvmult_add(aux, dst.block(row_num));
- else
- m->matrix->vmult_add(aux, dst.block(row_num));
- aux *= m->prefix;
- }
- dst.block(row_num) = aux;
- }
-}
-
-
-
-template <typename number, typename BlockVectorType>
-void
-BlockTrianglePrecondition<number, BlockVectorType>::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()));
-
- BlockVectorType aux;
- aux.reinit(dst);
- vmult(aux, src);
- dst += aux;
-}
-
-
-
-template <typename number, typename BlockVectorType>
-void
-BlockTrianglePrecondition<number, BlockVectorType>::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()));
-
- dst.equ(1., src);
-
- if (backward)
- {
- for (unsigned int i = n_block_rows(); i > 0;)
- do_row(dst, --i);
- }
- else
- {
- for (unsigned int i = 0; i < n_block_rows(); ++i)
- do_row(dst, i);
- }
-}
-
-template <typename number, typename BlockVectorType>
-void
-BlockTrianglePrecondition<number, BlockVectorType>::Tvmult(
- BlockVectorType &,
- const BlockVectorType &) const
-{
- Assert(false, ExcNotImplemented());
-}
-
-
-template <typename number, typename BlockVectorType>
-void
-BlockTrianglePrecondition<number, BlockVectorType>::Tvmult_add(
- BlockVectorType &,
- const BlockVectorType &) const
-{
- Assert(false, ExcNotImplemented());
-}
-
-template class BlockMatrixArray<float>;
-template class BlockMatrixArray<double>;
-template class BlockTrianglePrecondition<float>;
-template class BlockTrianglePrecondition<double>;
-
-#ifdef DEAL_II_WITH_TRILINOS
-template class BlockMatrixArray<float, TrilinosWrappers::MPI::BlockVector>;
-template class BlockMatrixArray<double, TrilinosWrappers::MPI::BlockVector>;
-template class BlockTrianglePrecondition<float,
- TrilinosWrappers::MPI::BlockVector>;
-template class BlockTrianglePrecondition<double,
- TrilinosWrappers::MPI::BlockVector>;
-#endif
-
-#ifdef DEAL_II_WITH_PETSC
-# ifdef PETSC_USE_COMPLEX
-template class BlockMatrixArray<std::complex<float>,
- PETScWrappers::MPI::BlockVector>;
-template class BlockMatrixArray<std::complex<double>,
- PETScWrappers::MPI::BlockVector>;
-template class BlockTrianglePrecondition<std::complex<float>,
- PETScWrappers::MPI::BlockVector>;
-template class BlockTrianglePrecondition<std::complex<double>,
- PETScWrappers::MPI::BlockVector>;
-# else
-template class BlockMatrixArray<float, PETScWrappers::MPI::BlockVector>;
-template class BlockMatrixArray<double, PETScWrappers::MPI::BlockVector>;
-template class BlockTrianglePrecondition<float,
- PETScWrappers::MPI::BlockVector>;
-template class BlockTrianglePrecondition<double,
- PETScWrappers::MPI::BlockVector>;
-# endif
-#endif
-
-DEAL_II_NAMESPACE_CLOSE
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2002 - 2018 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE.md at
-// the top level directory of deal.II.
-//
-// ---------------------------------------------------------------------
-
-#include <deal.II/lac/matrix_lib.templates.h>
-#include <deal.II/lac/sparse_matrix.h>
-
-DEAL_II_NAMESPACE_OPEN
-
-MeanValueFilter::MeanValueFilter(size_type component)
- : component(component)
-{}
-
-template void
-MeanValueFilter::filter(Vector<float> &) const;
-template void
-MeanValueFilter::filter(Vector<double> &) const;
-template void
-MeanValueFilter::filter(BlockVector<float> &) const;
-template void
-MeanValueFilter::filter(BlockVector<double> &) const;
-template void
-MeanValueFilter::vmult(Vector<float> &, const Vector<float> &) const;
-template void
-MeanValueFilter::vmult(Vector<double> &, const Vector<double> &) const;
-template void
-MeanValueFilter::vmult(BlockVector<float> &, const BlockVector<float> &) const;
-template void
-MeanValueFilter::vmult(BlockVector<double> &,
- const BlockVector<double> &) const;
-
-template void
-MeanValueFilter::vmult_add(Vector<float> &, const Vector<float> &) const;
-template void
-MeanValueFilter::vmult_add(Vector<double> &, const Vector<double> &) const;
-template void
-MeanValueFilter::vmult_add(BlockVector<float> &,
- const BlockVector<float> &) const;
-template void
-MeanValueFilter::vmult_add(BlockVector<double> &,
- const BlockVector<double> &) const;
-
-DEAL_II_NAMESPACE_CLOSE
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2005 - 2018 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE.md at
-// the top level directory of deal.II.
-//
-// ---------------------------------------------------------------------
-
-
-
-// This tests the construction of a BlockMatrixArray and outputs the
-// entered blocks using print_latex.
-
-#include <deal.II/lac/block_matrix_array.h>
-#include <deal.II/lac/full_matrix.h>
-
-#include "../tests.h"
-
-
-int
-main()
-{
- initlog();
-
- FullMatrix<double> A1(4, 4);
- FullMatrix<double> A2(4, 4);
- FullMatrix<double> B(4, 3);
- FullMatrix<double> C(3, 3);
-
- BlockMatrixArray<double> block(2, 2);
-
- block.enter(A1, 0, 0);
- block.enter(A2, 0, 0, 2, true);
- block.enter(B, 0, 1, -3.);
- block.enter(B, 0, 1, -3., true);
- block.enter(C, 1, 1, 1., true);
-
- block.print_latex(deallog);
-
- return 0;
-}
+++ /dev/null
-
-DEAL::\begin{array}{cc}
-DEAL:: M0+2xM1^T & -3xM2-3xM3^T\\
-DEAL:: & M4^T
-DEAL::\end{array}
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria.h>
-#include <deal.II/lac/block_matrix_array.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/block_sparsity_pattern.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/filtered_matrix.h>
#include <deal.II/lac/linear_operator.h>
-#include <deal.II/lac/matrix_lib.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/vector.h>
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2000 - 2018 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE.md at
-// the top level directory of deal.II.
-//
-// ---------------------------------------------------------------------
-
-
-
-#include <deal.II/base/mg_level_object.h>
-
-#include <deal.II/lac/block_matrix_array.h>
-#include <deal.II/lac/block_vector.h>
-#include <deal.II/lac/vector.h>
-
-#include <deal.II/multigrid/mg_block_smoother.h>
-
-#include <algorithm>
-
-#include "../tests.h"
-
-using namespace std;
-
-template <typename number>
-class ScalingMatrix : public Subscriptor
-{
-public:
- /**
- * Constructor setting the
- * scaling factor. Default is
- * constructing the identity
- * matrix.
- */
- ScalingMatrix(number scaling_factor = 1.);
- /**
- * Apply preconditioner.
- */
- template <typename VectorType>
- void
- vmult(VectorType &, const VectorType &) const;
-
- /**
- * Apply transpose
- * preconditioner. Since this is
- * the identity, this function is
- * the same as
- * vmult().
- */
- template <typename VectorType>
- void
- Tvmult(VectorType &, const VectorType &) const;
- /**
- * Apply preconditioner, adding to the previous value.
- */
- template <typename VectorType>
- void
- vmult_add(VectorType &, const VectorType &) const;
-
- /**
- * Apply transpose
- * preconditioner, adding. Since this is
- * the identity, this function is
- * the same as
- * vmult_add().
- */
- template <typename VectorType>
- void
- Tvmult_add(VectorType &, const VectorType &) const;
-
-private:
- number factor;
-};
-
-
-//----------------------------------------------------------------------//
-
-template <typename number>
-ScalingMatrix<number>::ScalingMatrix(number factor)
- : factor(factor)
-{}
-
-
-template <typename number>
-template <typename VectorType>
-inline void
-ScalingMatrix<number>::vmult(VectorType &dst, const VectorType &src) const
-{
- dst.equ(factor, src);
-}
-
-template <typename number>
-template <typename VectorType>
-inline void
-ScalingMatrix<number>::Tvmult(VectorType &dst, const VectorType &src) const
-{
- dst.equ(factor, src);
-}
-
-template <typename number>
-template <typename VectorType>
-inline void
-ScalingMatrix<number>::vmult_add(VectorType &dst, const VectorType &src) const
-{
- dst.add(factor, src);
-}
-
-
-
-template <typename number>
-template <typename VectorType>
-inline void
-ScalingMatrix<number>::Tvmult_add(VectorType &dst, const VectorType &src) const
-{
- dst.add(factor, src);
-}
-
-//----------------------------------------------------------------------//
-
-template <typename MatrixType, class RELAX>
-void
-check_smoother(const MGLevelObject<MatrixType> &m,
- const MGLevelObject<RELAX> & r)
-{
- GrowingVectorMemory<BlockVector<double>> mem;
- MGSmootherBlock<MatrixType, RELAX, double> smoother;
-
- smoother.initialize(m, r);
-
- for (unsigned int l = m.min_level(); l <= m.max_level(); ++l)
- {
- deallog << "Level " << l << std::endl;
-
- BlockVector<double> &u = *mem.alloc();
- BlockVector<double> &f = *mem.alloc();
- u.reinit(m[l].n_block_rows(), 3);
- f.reinit(u);
- for (unsigned int b = 0; b < f.n_blocks(); ++b)
- for (unsigned int i = 0; i < f.block(b).size(); ++i)
- f.block(b)(i) = (b + 1) * (i + l);
-
- deallog << "First step" << std::endl;
- smoother.set_steps(1);
- smoother.smooth(l, u, f);
-
- for (unsigned int b = 0; b < u.n_blocks(); ++b)
- {
- for (unsigned int i = 0; i < u.block(b).size(); ++i)
- deallog << '\t' << (int)(u.block(b)(i) + .5);
- deallog << std::endl;
- }
-
- deallog << "Second step" << std::endl;
- smoother.smooth(l, u, f);
-
- for (unsigned int b = 0; b < u.n_blocks(); ++b)
- {
- for (unsigned int i = 0; i < u.block(b).size(); ++i)
- deallog << '\t' << (int)(u.block(b)(i) + .5);
- deallog << std::endl;
- }
-
- deallog << "Two steps" << std::endl;
- u = 0.;
- smoother.set_steps(2);
- smoother.smooth(l, u, f);
-
- for (unsigned int b = 0; b < u.n_blocks(); ++b)
- {
- for (unsigned int i = 0; i < u.block(b).size(); ++i)
- deallog << '\t' << (int)(u.block(b)(i) + .5);
- deallog << std::endl;
- }
-
- mem.free(&u);
- mem.free(&f);
- }
-}
-
-void
-check()
-{
- ScalingMatrix<double> s1(-1.);
- ScalingMatrix<double> s2(2.);
- ScalingMatrix<double> s8(8.);
-
- GrowingVectorMemory<Vector<double>> mem;
- MGLevelObject<BlockMatrixArray<double>> A(2, 4);
- MGLevelObject<BlockTrianglePrecondition<double>> P(2, 4);
-
- for (unsigned int l = A.min_level(); l <= A.max_level(); ++l)
- {
- A[l].initialize(3, 3);
- P[l].reinit(3);
- for (unsigned int b = 0; b < A[l].n_block_rows(); ++b)
- {
- P[l].enter(s2, b, b, A[l].n_block_rows() - b);
- A[l].enter(s8, b, b, 1);
- for (unsigned int b2 = 0; b2 < A[l].n_block_rows(); ++b2)
- A[l].enter(s1, b, b2, 1.);
- }
- }
-
- check_smoother(A, P);
-}
-
-
-int
-main()
-{
- initlog();
- deallog << std::setprecision(3);
-
- check();
-}
+++ /dev/null
-
-DEAL::Level 2
-DEAL::First step
-DEAL:: 12 18 24
-DEAL:: 16 24 32
-DEAL:: 12 18 24
-DEAL::Second step
-DEAL:: -311 -467 -623
-DEAL:: -319 -479 -639
-DEAL:: -87 -131 -175
-DEAL::Two steps
-DEAL:: -311 -467 -623
-DEAL:: -319 -479 -639
-DEAL:: -87 -131 -175
-DEAL::Level 3
-DEAL::First step
-DEAL:: 18 24 30
-DEAL:: 24 32 40
-DEAL:: 18 24 30
-DEAL::Second step
-DEAL:: -467 -623 -779
-DEAL:: -479 -639 -799
-DEAL:: -131 -175 -219
-DEAL::Two steps
-DEAL:: -467 -623 -779
-DEAL:: -479 -639 -799
-DEAL:: -131 -175 -219
-DEAL::Level 4
-DEAL::First step
-DEAL:: 24 30 36
-DEAL:: 32 40 48
-DEAL:: 24 30 36
-DEAL::Second step
-DEAL:: -623 -779 -935
-DEAL:: -639 -799 -959
-DEAL:: -175 -219 -263
-DEAL::Two steps
-DEAL:: -623 -779 -935
-DEAL:: -639 -799 -959
-DEAL:: -175 -219 -263