#include <base/smartpointer.h>
#include <base/table.h>
+#include <lac/vector_memory.h>
+
#include <vector>
#include <map>
#include <string>
# include <strstream>
#endif
-//TODO:[GK] Obtain aux vector from VectorMemory
template <typename> class BlockVector;
template <typename> class Vector;
* <h3>Requirements on MATRIX</h3>
*
* The template argument <tt>MATRIX</tt> is a class providing the
- * matrix-vector multiplication functions <tt>vmult</tt>,
- * <tt>Tvmult</tt>, <tt>vmult_add</tt> and <tt>Tvmult_add</tt> used in
- * this class, but with arguments of type <tt>VECTOR</tt> instead of
- * <tt>BlockVector<VECTOR></tt>. SparseMatrix is a possible entry
- * type.
+ * 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>. SparseMatrix is a possible entry type.
*
* <h3>Example program</h3>
* We document the relevant parts of <tt>examples/doxygen/block_matrix_array.cc</tt>.
* @skip main
* @until C.fill
*
+ * The BlockMatrixArray needs a VectorMemory<Vector<number>
+ * > object to allocate auxiliary vectors. <tt>number</tt> must
+ * equal the second template argument of BlockMatrixArray and also the
+ * number type of the BlockVector used later. We use the
+ * GrowingVectorMemory type, since it remembers the vector and avoids
+ * reallocating.
+ *
+ * @ line Growing
+ *
* 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
*
* @author Guido Kanschat, 2000 - 2005
*/
-template <class MATRIX>
+template <class MATRIX, typename number = double>
class BlockMatrixArray : public Subscriptor
{
public:
* dimensions.
*/
BlockMatrixArray (const unsigned int n_block_rows,
- const unsigned int n_block_cols);
+ const unsigned int n_block_cols,
+ VectorMemory<Vector<number> >& mem);
/**
* Add a block matrix entry. The
/**
* Matrix-vector multiplication.
*/
- template <class number>
void vmult (BlockVector<number>& dst,
const BlockVector<number>& src) const;
* Matrix-vector multiplication
* adding to <tt>dst</tt>.
*/
- template <class number>
void vmult_add (BlockVector<number>& dst,
const BlockVector<number>& src) const;
* Transposed matrix-vector
* multiplication.
*/
- template <class number>
void Tvmult (BlockVector<number>& dst,
const BlockVector<number>& src) const;
* multiplication adding to
* <tt>dst</tt>.
*/
- template <class number>
void Tvmult_add (BlockVector<number>& dst,
const BlockVector<number>& src) const;
+ /**
+ * Matrix scalar product between
+ * two vectors (at least for a
+ * symmetric matrix).
+ */
+ number matrix_scalar_product (const BlockVector<number>& u,
+ const BlockVector<number>& v) const;
+
+ /**
+ * Square of the matrix norm of a
+ * vector (at least for a
+ * symmetric matrix).
+ */
+ number matrix_norm_square (const BlockVector<number>& u) const;
+
/**
* Print the block structure as a
* LaTeX-array. This output will
* number of blocks per row.
*/
unsigned int block_cols;
+ protected:
+ /**
+ * The memory used for auxiliary
+ * vectors.
+ */
+ mutable SmartPointer<VectorMemory<Vector<number> > > mem;
};
/*@}*/
*
* @author Guido Kanschat, 2001, 2005
*/
-template <class MATRIX>
-class BlockTrianglePrecondition : private BlockMatrixArray<MATRIX>
+template <class MATRIX, typename number = double>
+class BlockTrianglePrecondition
+ : private BlockMatrixArray<MATRIX, number>
{
public:
/**
* insertion instead of forward.
*/
BlockTrianglePrecondition(unsigned int block_rows,
+ VectorMemory<Vector<number> >& mem,
bool backward = false);
/**
* Preconditioning.
*/
- template <class number>
void vmult (BlockVector<number>& dst,
const BlockVector<number>& src) const;
* Preconditioning
* adding to <tt>dst</tt>.
*/
- template <class number>
void vmult_add (BlockVector<number>& dst,
const BlockVector<number>& src) const;
/**
* Transposed preconditioning
*/
- template <class number>
void Tvmult (BlockVector<number>& dst,
const BlockVector<number>& src) const;
* Transposed preconditioning
* adding to <tt>dst</tt>.
*/
- template <class number>
void Tvmult_add (BlockVector<number>& dst,
const BlockVector<number>& src) const;
/**
* Make function of base class available.
*/
- BlockMatrixArray<MATRIX>::enter;
+ BlockMatrixArray<MATRIX, number>::enter;
/**
* Make function of base class available.
*/
- BlockMatrixArray<MATRIX>::print_latex;
+ BlockMatrixArray<MATRIX, number>::print_latex;
/**
* Make function of base class available.
*/
- BlockMatrixArray<MATRIX>::n_block_rows;
+ BlockMatrixArray<MATRIX, number>::n_block_rows;
/**
* Make function of base class available.
*/
- BlockMatrixArray<MATRIX>::n_block_cols;
- BlockMatrixArray<MATRIX>::clear;
- BlockMatrixArray<MATRIX>::Subscriptor::subscribe;
- BlockMatrixArray<MATRIX>::Subscriptor::unsubscribe;
+ BlockMatrixArray<MATRIX, number>::n_block_cols;
+ BlockMatrixArray<MATRIX, number>::clear;
+ BlockMatrixArray<MATRIX, number>::Subscriptor::subscribe;
+ BlockMatrixArray<MATRIX, number>::Subscriptor::unsubscribe;
/** @addtogroup Exceptions
* @{ */
* entry of the diagonal element
* for one row.
*/
- template <class number>
void do_row (BlockVector<number>& dst,
unsigned int row_num) const;
///@if NoDoc
//----------------------------------------------------------------------//
-template <class MATRIX>
+template <class MATRIX, typename number>
inline
-BlockMatrixArray<MATRIX>::Entry::Entry (const MATRIX& matrix,
+BlockMatrixArray<MATRIX,number>::Entry::Entry (const MATRIX& matrix,
unsigned row, unsigned int col,
double prefix, bool transpose)
:
-template <class MATRIX>
+template <class MATRIX, typename number>
inline
-BlockMatrixArray<MATRIX>::BlockMatrixArray (const unsigned int n_block_rows,
- const unsigned int n_block_cols)
+BlockMatrixArray<MATRIX,number>::BlockMatrixArray (
+ const unsigned int n_block_rows,
+ const unsigned int n_block_cols,
+ VectorMemory<Vector<number> >& mem)
: block_rows (n_block_rows),
- block_cols (n_block_cols)
+ block_cols (n_block_cols),
+ mem(&mem)
{}
-template <class MATRIX>
+template <class MATRIX, typename number>
inline
void
-BlockMatrixArray<MATRIX>::enter (const MATRIX& matrix,
+BlockMatrixArray<MATRIX,number>::enter (const MATRIX& matrix,
unsigned row, unsigned int col,
double prefix, bool transpose)
{
-template <class MATRIX>
+template <class MATRIX, typename number>
inline
void
-BlockMatrixArray<MATRIX>::clear ()
+BlockMatrixArray<MATRIX,number>::clear ()
{
entries.clear();
}
-template <class MATRIX>
-template <class number>
+template <class MATRIX, typename number>
inline
void
-BlockMatrixArray<MATRIX>::vmult_add (BlockVector<number>& dst,
+BlockMatrixArray<MATRIX,number>::vmult_add (BlockVector<number>& dst,
const BlockVector<number>& src) const
{
Assert (dst.n_blocks() == block_rows,
Assert (src.n_blocks() == block_cols,
ExcDimensionMismatch(src.n_blocks(), block_cols));
- static Vector<number> aux;
+ Vector<number>* p_aux = mem->alloc();
+ Vector<number>& aux = *p_aux;
typename std::vector<Entry>::const_iterator m = entries.begin();
typename std::vector<Entry>::const_iterator end = entries.end();
m->matrix->vmult(aux, src.block(m->col));
dst.block(m->row).add (m->prefix, aux);
}
+ mem->free(p_aux);
}
-template <class MATRIX>
-template <class number>
+template <class MATRIX, typename number>
inline
void
-BlockMatrixArray<MATRIX>::vmult (BlockVector<number>& dst,
+BlockMatrixArray<MATRIX,number>::vmult (BlockVector<number>& dst,
const BlockVector<number>& src) const
{
dst = 0.;
-template <class MATRIX>
-template <class number>
+template <class MATRIX, typename number>
inline
void
-BlockMatrixArray<MATRIX>::Tvmult_add (BlockVector<number>& dst,
+BlockMatrixArray<MATRIX,number>::Tvmult_add (BlockVector<number>& dst,
const BlockVector<number>& src) const
{
Assert (dst.n_blocks() == block_cols,
typename std::vector<Entry>::const_iterator m = entries.begin();
typename std::vector<Entry>::const_iterator end = entries.end();
- static Vector<number> aux;
+ Vector<number>* p_aux = mem->alloc();
+ Vector<number>& aux = *p_aux;
for (; m != end ; ++m)
{
m->matrix->Tvmult(aux, src.block(m->row));
dst.block(m->col).add (m->prefix, aux);
}
+ mem->free(p_aux);
}
-template <class MATRIX>
-template <class number>
+template <class MATRIX, typename number>
inline
void
-BlockMatrixArray<MATRIX>::Tvmult (BlockVector<number>& dst,
+BlockMatrixArray<MATRIX,number>::Tvmult (BlockVector<number>& dst,
const BlockVector<number>& src) const
{
dst = 0.;
-template <class MATRIX>
+template <class MATRIX, typename number>
+inline
+number
+BlockMatrixArray<MATRIX,number>::matrix_scalar_product (
+ const BlockVector<number>& u,
+ const BlockVector<number>& v) const
+{
+ Assert (u.n_blocks() == block_rows,
+ ExcDimensionMismatch(u.n_blocks(), block_rows));
+ Assert (v.n_blocks() == block_cols,
+ ExcDimensionMismatch(v.n_blocks(), block_cols));
+
+ Vector<number>* p_aux = mem->alloc();
+ Vector<number>& 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;
+ }
+ mem->free(p_aux);
+
+ return result;
+}
+
+
+
+template <class MATRIX, typename number>
+inline
+number
+BlockMatrixArray<MATRIX,number>::matrix_norm_square (
+ const BlockVector<number>& u) const
+{
+ return matrix_scalar_product(u,u);
+}
+
+
+
+template <class MATRIX, typename number>
inline
unsigned int
-BlockMatrixArray<MATRIX>::n_block_rows () const
+BlockMatrixArray<MATRIX,number>::n_block_rows () const
{
return block_rows;
}
-template <class MATRIX>
+template <class MATRIX, typename number>
inline
unsigned int
-BlockMatrixArray<MATRIX>::n_block_cols () const
+BlockMatrixArray<MATRIX,number>::n_block_cols () const
{
return block_cols;
}
-template <class MATRIX>
+template <class MATRIX, typename number>
template <class STREAM>
inline
void
-BlockMatrixArray<MATRIX>::print_latex (STREAM& out) const
+BlockMatrixArray<MATRIX,number>::print_latex (STREAM& out) const
{
out << "\\begin{array}{"
<< std::string(n_block_cols(), 'c')
//----------------------------------------------------------------------//
-template <class MATRIX>
-BlockTrianglePrecondition<MATRIX>::BlockTrianglePrecondition(
+template <class MATRIX, typename number>
+BlockTrianglePrecondition<MATRIX,number>::BlockTrianglePrecondition(
unsigned int block_rows,
+ VectorMemory<Vector<number> >& mem,
bool backward)
:
- BlockMatrixArray<MATRIX> (block_rows, block_rows),
- backward(backward)
+ BlockMatrixArray<MATRIX,number> (block_rows, block_rows, mem),
+ backward(backward)
{}
-template <class MATRIX>
-template <class number>
+template <class MATRIX, typename number>
inline
void
-BlockTrianglePrecondition<MATRIX>::do_row (
+BlockTrianglePrecondition<MATRIX,number>::do_row (
BlockVector<number>& dst,
unsigned int row_num) const
{
- typename std::vector<typename BlockMatrixArray<MATRIX>::Entry>::const_iterator
+ typename std::vector<typename BlockMatrixArray<MATRIX,number>::Entry>::const_iterator
m = this->entries.begin();
- typename std::vector<typename BlockMatrixArray<MATRIX>::Entry>::const_iterator
+ typename std::vector<typename BlockMatrixArray<MATRIX,number>::Entry>::const_iterator
end = this->entries.end();
- typename std::vector<typename BlockMatrixArray<MATRIX>::Entry>::const_iterator
+ typename std::vector<typename BlockMatrixArray<MATRIX,number>::Entry>::const_iterator
diagonal = end;
- static Vector<number> aux;
+ Vector<number>* p_aux = mem->alloc();
+ Vector<number>& aux = *p_aux;
+
aux.reinit(dst.block(row_num), true);
for (; m != end ; ++m)
diagonal->matrix->vmult(aux, dst.block(row_num));
dst.block(row_num).equ (diagonal->prefix, aux);
+ mem->free(p_aux);
}
-template <class MATRIX>
-template <class number>
+template <class MATRIX, typename number>
inline
void
-BlockTrianglePrecondition<MATRIX>::vmult_add (
+BlockTrianglePrecondition<MATRIX,number>::vmult_add (
BlockVector<number>& dst,
const BlockVector<number>& src) const
{
-template <class MATRIX>
-template <class number>
+template <class MATRIX, typename number>
inline
void
-BlockTrianglePrecondition<MATRIX>::vmult (
+BlockTrianglePrecondition<MATRIX,number>::vmult (
BlockVector<number>& dst,
const BlockVector<number>& src) const
{
}
-template <class MATRIX>
-template <class number>
+template <class MATRIX, typename number>
inline
void
-BlockTrianglePrecondition<MATRIX>::Tvmult (
+BlockTrianglePrecondition<MATRIX,number>::Tvmult (
BlockVector<number>&,
const BlockVector<number>&) const
{
}
-template <class MATRIX>
-template <class number>
+template <class MATRIX, typename number>
inline
void
-BlockTrianglePrecondition<MATRIX>::Tvmult_add (
+BlockTrianglePrecondition<MATRIX,number>::Tvmult_add (
BlockVector<number>&,
const BlockVector<number>&) const
{