}
if (this->store_diagonals())
this->diagonal(0) = M_cell;
- this->inverse(0).invert(M_cell);
+ switch (this->inversion)
+ {
+ case PreconditionBlockBase<inverse_type>::gauss_jordan:
+ this->inverse(0).invert(M_cell);
+ break;
+ case PreconditionBlockBase<inverse_type>::householder:
+ this->inverse_householder(0).initialize(M_cell);
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+
+ }
}
else
{
if (this->store_diagonals())
this->diagonal(cell) = M_cell;
- this->inverse(cell).invert(M_cell);
+ switch (this->inversion)
+ {
+ case PreconditionBlockBase<inverse_type>::gauss_jordan:
+ this->inverse(cell).invert(M_cell);
+ break;
+ case PreconditionBlockBase<inverse_type>::householder:
+ this->inverse_householder(cell).initialize(M_cell);
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+
+ }
}
}
this->inverses_computed(true);
if (this->inverses_ready())
{
if (transpose_diagonal)
- this->inverse(cell).Tvmult(x_cell, b_cell);
+ this->inverse_Tvmult(cell, x_cell, b_cell);
else
- this->inverse(cell).vmult(x_cell, b_cell);
+ this->inverse_vmult(cell, x_cell, b_cell);
}
else
{
if (this->inverses_ready())
{
if (transpose_diagonal)
- this->inverse(cell).Tvmult(x_cell, b_cell);
+ this->inverse_Tvmult(cell, x_cell, b_cell);
else
- this->inverse(cell).vmult(x_cell, b_cell);
+ this->inverse_vmult(cell, x_cell, b_cell);
}
else
{
{
b_cell(row_cell)=src(row);
}
- this->inverse(cell).vmult(x_cell, b_cell);
+ this->inverse_vmult(cell, x_cell, b_cell);
// distribute x_cell to dst
for (row=cell*this->blocksize, row_cell=0;
row_cell<this->blocksize;
if (this->inverses_ready())
{
if (transpose_diagonal)
- this->inverse(cell).Tvmult(x_cell, b_cell);
+ this->inverse_Tvmult(cell, x_cell, b_cell);
else
- this->inverse(cell).vmult(x_cell, b_cell);
+ this->inverse_vmult(cell, x_cell, b_cell);
}
else
{
if (this->inverses_ready())
{
if (transpose_diagonal)
- this->inverse(cell).Tvmult(x_cell, b_cell);
+ this->inverse_Tvmult(cell, x_cell, b_cell);
else
- this->inverse(cell).vmult(x_cell, b_cell);
+ this->inverse_vmult(cell, x_cell, b_cell);
}
else
{
#include <base/subscriptor.h>
#include <base/smartpointer.h>
#include <base/memory_consumption.h>
+#include <lac/householder.h>
#include <vector>
class PreconditionBlockBase
{
public:
+ /**
+ * Choose a method for inverting
+ * the blocks, and thus the data
+ * type for the inverses.
+ */
+ enum Inversion
+ {
+ /**
+ * Use the standard
+ * Gauss-Jacobi method
+ * implemented in FullMatrix::inverse().
+ */
+ gauss_jordan,
+ /**
+ * Use QR decomposition of
+ * the Householder class.
+ */
+ householder
+ };
+
/**
* Constructor initializing
* default values.
*/
- PreconditionBlockBase(bool store_diagonals = false);
+ PreconditionBlockBase(bool store_diagonals = false,
+ Inversion method = gauss_jordan);
/**
* The virtual destructor
* then only one block will be
* stored.
*/
- void reinit(unsigned int nblocks, unsigned int blocksize, bool compress);
+ void reinit(unsigned int nblocks, unsigned int blocksize, bool compress,
+ Inversion method = gauss_jordan);
/**
* Tell the class that inverses
* are stored.
*/
number el(unsigned int i, unsigned int j) const;
+
+ /**
+ * Multiply with the inverse
+ * block at position <tt>i</tt>.
+ */
+ template <typename number2>
+ void inverse_vmult(unsigned int i, Vector<number2>& dst, const Vector<number2>& src) const;
+
+ /**
+ * Multiply with the transposed inverse
+ * block at position <tt>i</tt>.
+ */
+ template <typename number2>
+ void inverse_Tvmult(unsigned int i, Vector<number2>& dst, const Vector<number2>& src) const;
/**
* Access to the inverse diagonal
- * blocks.
+ * blocks if Inversion is #gauss_jordan.
*/
FullMatrix<number>& inverse (unsigned int i);
+ /**
+ * Access to the inverse diagonal
+ * blocks if Inversion is #householder.
+ */
+ Householder<number>& inverse_householder (unsigned int i);
+
/**
* Access to the inverse diagonal
* blocks.
*/
DeclException0 (ExcDiagonalsNotStored);
+ protected:
+ /**
+ * The method used for inverting blocks.
+ */
+ Inversion inversion;
+
private:
/**
* The number of (inverse)
*/
unsigned int n_diagonal_blocks;
- /**
+ /**
* Storage of the inverse
* matrices of the diagonal
* blocks matrices as
* <tt>FullMatrix<number></tt>
- * matrices. Using
+ * matrices, if Inversion
+ * #gauss_jordan is used. Using
* <tt>number=float</tt> saves
* memory in comparison with
* <tt>number=double</tt>.
*/
- std::vector<FullMatrix<number> > var_inverse;
+ std::vector<FullMatrix<number> > var_inverse_full;
+
+ /**
+ * Storage of the inverse
+ * matrices of the diagonal
+ * blocks matrices as
+ * <tt>Householder<number></tt>
+ * matrices if Inversion
+ * #householder is used. Using
+ * <tt>number=float</tt> saves
+ * memory in comparison with
+ * <tt>number=double</tt>.
+ */
+ std::vector<Householder<number> > var_inverse_householder;
/**
* Storage of the original diagonal blocks.
template <typename number>
inline
-PreconditionBlockBase<number>::PreconditionBlockBase(bool store)
+PreconditionBlockBase<number>::PreconditionBlockBase(
+ bool store, Inversion method)
:
+ inversion(method),
n_diagonal_blocks(0),
var_store_diagonals(store),
var_same_diagonal(false),
void
PreconditionBlockBase<number>::clear()
{
- if (var_inverse.size()!=0)
- var_inverse.erase(var_inverse.begin(), var_inverse.end());
+ if (var_inverse_full.size()!=0)
+ var_inverse_full.erase(var_inverse_full.begin(), var_inverse_full.end());
+ if (var_inverse_full.size()!=0)
+ var_inverse_householder.erase(var_inverse_householder.begin(), var_inverse_householder.end());
if (var_diagonal.size()!=0)
var_diagonal.erase(var_diagonal.begin(), var_diagonal.end());
var_same_diagonal = false;
template <typename number>
inline
void
-PreconditionBlockBase<number>::reinit(unsigned int n, unsigned int b, bool compress)
+PreconditionBlockBase<number>::reinit(unsigned int n, unsigned int b, bool compress,
+Inversion method)
{
+ inversion = method;
var_same_diagonal = compress;
var_inverses_ready = false;
n_diagonal_blocks = n;
if (compress)
{
- var_inverse.resize(1);
- var_inverse[0].reinit(b,b);
+ switch(inversion)
+ {
+ case gauss_jordan:
+ var_inverse_full.resize(1);
+ var_inverse_full[0].reinit(b,b);
+ break;
+ case householder:
+ var_inverse_householder.resize(1);
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+
if (store_diagonals())
{
var_diagonal.resize(1);
tmp(n, FullMatrix<number>(b));
var_diagonal.swap (tmp);
}
-
- if (true)
+
+ switch(inversion)
{
- std::vector<FullMatrix<number> >
- tmp(n, FullMatrix<number>(b));
- var_inverse.swap (tmp);
- // make sure the tmp object
- // goes out of scope as
- // soon as possible
+ case gauss_jordan:
+ {
+ std::vector<FullMatrix<number> >
+ tmp(n, FullMatrix<number>(b));
+ var_inverse_full.swap (tmp);
+ break;
+ }
+ case householder:
+ var_inverse_householder.resize(n);
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
}
}
}
return n_diagonal_blocks;
}
+template <typename number>
+template <typename number2>
+inline
+void
+PreconditionBlockBase<number>::inverse_vmult(
+ unsigned int i, Vector<number2>& dst, const Vector<number2>& src) const
+{
+ const unsigned int ii = same_diagonal() ? 0U : i;
+
+ switch (inversion)
+ {
+ case gauss_jordan:
+ AssertIndexRange (ii, var_inverse_full.size());
+ var_inverse_full[ii].vmult(dst, src);
+ break;
+ case householder:
+ AssertIndexRange (ii, var_inverse_householder.size());
+ var_inverse_householder[ii].vmult(dst, src);
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+}
+
+
+template <typename number>
+template <typename number2>
+inline
+void
+PreconditionBlockBase<number>::inverse_Tvmult(
+ unsigned int i, Vector<number2>& dst, const Vector<number2>& src) const
+{
+ const unsigned int ii = same_diagonal() ? 0U : i;
+
+ switch (inversion)
+ {
+ case gauss_jordan:
+ AssertIndexRange (ii, var_inverse_full.size());
+ var_inverse_full[ii].Tvmult(dst, src);
+ break;
+ case householder:
+ AssertIndexRange (ii, var_inverse_householder.size());
+ var_inverse_householder[ii].Tvmult(dst, src);
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+}
+
template <typename number>
inline
PreconditionBlockBase<number>::inverse(unsigned int i) const
{
if (same_diagonal())
- return var_inverse[0];
+ return var_inverse_full[0];
- Assert (i < var_inverse.size(), ExcIndexRange(i,0,var_inverse.size()));
- return var_inverse[i];
+ Assert (i < var_inverse_full.size(), ExcIndexRange(i,0,var_inverse_full.size()));
+ return var_inverse_full[i];
}
PreconditionBlockBase<number>::inverse(unsigned int i)
{
if (same_diagonal())
- return var_inverse[0];
+ return var_inverse_full[0];
+
+ Assert (i < var_inverse_full.size(), ExcIndexRange(i,0,var_inverse_full.size()));
+ return var_inverse_full[i];
+}
+
+
+template <typename number>
+inline
+Householder<number>&
+PreconditionBlockBase<number>::inverse_householder(unsigned int i)
+{
+ if (same_diagonal())
+ return var_inverse_householder[0];
- Assert (i < var_inverse.size(), ExcIndexRange(i,0,var_inverse.size()));
- return var_inverse[i];
+ AssertIndexRange (i, var_inverse_householder.size());
+ return var_inverse_householder[i];
}
PreconditionBlockBase<number>::memory_consumption () const
{
unsigned int mem = sizeof(*this);
- for (unsigned int i=0; i<var_inverse.size(); ++i)
- mem += MemoryConsumption::memory_consumption(var_inverse[i]);
+ for (unsigned int i=0; i<var_inverse_full.size(); ++i)
+ mem += MemoryConsumption::memory_consumption(var_inverse_full[i]);
for (unsigned int i=0; i<var_diagonal.size(); ++i)
mem += MemoryConsumption::memory_consumption(var_diagonal[i]);
return mem;
:
relaxation(1.),
invert_diagonal(true),
- same_diagonal(false)
+ same_diagonal(false),
+ inversion(PreconditionBlockBase<inverse_type>::gauss_jordan)
{}
block_list(&bl),
relaxation(relaxation),
invert_diagonal(invert_diagonal),
- same_diagonal(same_diagonal)
+ same_diagonal(same_diagonal),
+ inversion(PreconditionBlockBase<inverse_type>::gauss_jordan)
{}
A = &M;
additional_data = parameters;
- this->reinit(additional_data.block_list->size(), 0, additional_data.same_diagonal);
+ this->reinit(additional_data.block_list->size(), 0, additional_data.same_diagonal,
+ additional_data.inversion);
if (parameters.invert_diagonal)
invert_diagblocks();
this->diagonal(block).reinit(bs, bs);
this->diagonal(block) = M_cell;
}
- this->inverse(block).reinit(bs, bs);
- this->inverse(block).invert(M_cell);
+ switch(this->inversion)
+ {
+ case PreconditionBlockBase<inverse_type>::gauss_jordan:
+ this->inverse(block).reinit(bs, bs);
+ this->inverse(block).invert(M_cell);
+ break;
+ case PreconditionBlockBase<inverse_type>::householder:
+ this->inverse_householder(block).initialize(M_cell);
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
}
}
this->inverses_computed(true);
b_cell(row_cell) -= entry->value() * prev(entry->column());
}
// Apply inverse diagonal
- this->inverse(block).vmult(x_cell, b_cell);
+ this->inverse_vmult(block, x_cell, b_cell);
// Store in result vector
row=additional_data.block_list->begin(block);
for (unsigned int row_cell=0; row_cell<bs; ++row_cell, ++row)