* second step, the inverses of
* the diagonal blocks may be
* computed.
+ *
+ * Additionally, a relaxation
+ * parameter for derived classes
+ * may be provided.
*/
void initialize (const SparseMatrix<number>& A,
- const unsigned int block_size);
+ const unsigned int block_size,
+ const double relaxation = 1.);
/**
* Deletes the inverse diagonal
* Exception
*/
DeclException0 (ExcMatrixNotSquare);
+
+ /**
+ * Exception
+ */
+ DeclException0 (ExcDiagonalsNotStored);
protected:
/**
*/
const FullMatrix<inverse_type>& inverse (unsigned int i) const;
+ /**
+ * Access to the diagonal
+ * blocks.
+ */
+ const FullMatrix<inverse_type>& diagonal (unsigned int i) const;
+
/**
* Size of the blocks. Each
* diagonal block is assumed to
* derived classes.
*/
SmartPointer<const SparseMatrix<number> > A;
+ /**
+ * Relaxation parameter to be
+ * used by derived classes.
+ */
+ double relaxation;
+
+
+ /**
+ * Flag for storing the diagonal
+ * blocks of the matrix.
+ */
+ bool store_diagonals;
private:
/**
* matrices of the diagonal
* blocks matrices as
* @p{FullMatrix<inverse_type>}
- * matrices. For BlockSOR as
- * preconditioning using
+ * matrices. Using
* @p{inverse_type=float} saves
* memory in comparison with
* @p{inverse_type=double}.
*/
- vector<FullMatrix<inverse_type> > _inverse;
+ vector<FullMatrix<inverse_type> > var_inverse;
+ /**
+ * Storage of the original diagonal blocks.
+ * These are only filled if @p{store_diagonals}
+ * is @p{true}.
+ *
+ * Used by the blocked SSOR method.
+ */
+ vector<FullMatrix<inverse_type> > var_diagonal;
+
/**
* Flag for diagonal compression.
* @see set_same_diagonal()
/**
* Block SOR preconditioning.
*
- * The diagonal blocks and the elements
- * (of arbitray structure) below the diagonal blocks are used
- * in the @p{operator ()} function of this class.
+ * The functions @p{vmult} and @p{Tvmult} execute a (transposed)
+ * block-SOR step, based on the blocks in @ref{PreconditionBlock}. The
+ * elements outside the diagonal blocks may be distributed
+ * arbitrarily.
*
* See @ref{PreconditionBlock} for requirements on the matrix.
* @author Ralf Hartmann, Guido Kanschat, 1999, 2000
*/
template<typename number, typename inverse_type = number>
class PreconditionBlockSOR : public Subscriptor,
- private PreconditionBlock<number,inverse_type>
+ protected PreconditionBlock<number,inverse_type>
{
public:
- /**
- * Constructor. Takes the damping
- * Parameter as
- */
- PreconditionBlockSOR (const number omega=1.);
-
- /**
- * Destructor.
- */
- ~PreconditionBlockSOR();
-
/**
* Make initialization function
* publicly available.
*/
PreconditionBlock<number,inverse_type>::invert_diagblocks;
- /**
- * Set the relaxation parameter.
- */
- void set_omega(number omega);
-
/**
* Execute block SOR
* preconditioning.
*/
template <typename number2>
void Tvmult (Vector<number2>&, const Vector<number2>&) const;
+};
- private:
+
+/**
+ * Block SSOR preconditioning.
+ *
+ * The functions @p{vmult} and @p{Tvmult} execute a block-SSOR step,
+ * based on the implementation in @ref{PreconditionBlockSOR}. This
+ * class requires storage of the diagonal blocks and their inverses.
+ *
+ * See @ref{PreconditionBlock} for requirements on the matrix.
+ * @author Ralf Hartmann, Guido Kanschat, 1999, 2000
+ */
+template<typename number, typename inverse_type = number>
+class PreconditionBlockSSOR : private PreconditionBlockSOR<number,inverse_type>
+{
+ public:
/**
- * Relaxation parameter.
+ * Constructor.
+ */
+ PreconditionBlockSSOR ();
+ /**
+ * Make initialization function
+ * publicly available.
+ */
+ PreconditionBlock<number,inverse_type>::initialize;
+
+ /**
+ * Make function of base class public again.
+ */
+ PreconditionBlock<number,inverse_type>::clear;
+
+ /**
+ * Make function of base class public again.
+ */
+ PreconditionBlock<number,inverse_type>::set_same_diagonal;
+
+ /**
+ * Make function of base class public again.
*/
- number omega;
+ PreconditionBlock<number,inverse_type>::invert_diagblocks;
+
+ /**
+ * Execute block SSOR
+ * preconditioning.
+ *
+ * This function will
+ * automatically use the inverse
+ * matrices if they exist, if not
+ * then BlockSOR will waste much
+ * time inverting the diagonal
+ * block matrices in each
+ * preconditioning step.
+ */
+ template <typename number2>
+ void vmult (Vector<number2>&, const Vector<number2>&) const;
+
+ /**
+ * Same as @ref{vmult}
+ */
+ template <typename number2>
+ void Tvmult (Vector<number2>&, const Vector<number2>&) const;
};
PreconditionBlock<number,inverse_type>::PreconditionBlock ():
blocksize(0),
A(0),
+ store_diagonals(false),
same_diagonal(false)
{};
template <typename number, typename inverse_type>
PreconditionBlock<number,inverse_type>::~PreconditionBlock ()
{
- if (_inverse.size()!=0)
- _inverse.erase(_inverse.begin(), _inverse.end());
+ if (var_inverse.size()!=0)
+ var_inverse.erase(var_inverse.begin(), var_inverse.end());
+ if (var_diagonal.size()!=0)
+ var_diagonal.erase(var_diagonal.begin(), var_diagonal.end());
}
template <typename number, typename inverse_type>
void PreconditionBlock<number,inverse_type>::clear ()
{
- if (_inverse.size()!=0)
- _inverse.erase(_inverse.begin(), _inverse.end());
+ if (var_inverse.size()!=0)
+ var_inverse.erase(var_inverse.begin(), var_inverse.end());
+ if (var_diagonal.size()!=0)
+ var_diagonal.erase(var_diagonal.begin(), var_diagonal.end());
blocksize = 0;
same_diagonal = false;
A = 0;
template <typename number, typename inverse_type>
void PreconditionBlock<number,inverse_type>::initialize (const SparseMatrix<number> &M,
- unsigned int bsize)
+ unsigned int bsize,
+ const double relax)
{
clear();
Assert (M.m() == M.n(), ExcMatrixNotSquare());
Assert (bsize>0, ExcIndexRange(bsize, 1, M.m()));
Assert (A->m()%bsize==0, ExcWrongBlockSize(bsize, A->m()));
blocksize=bsize;
+ relaxation = relax;
}
PreconditionBlock<number,inverse_type>::inverse(unsigned int i) const
{
if (same_diagonal)
- return _inverse[0];
+ return var_inverse[0];
- Assert (i < _inverse.size(), ExcIndexRange(i,0,_inverse.size()));
- return _inverse[i];
+ Assert (i < var_inverse.size(), ExcIndexRange(i,0,var_inverse.size()));
+ return var_inverse[i];
+}
+
+
+template <typename number, typename inverse_type>
+const FullMatrix<inverse_type>&
+PreconditionBlock<number,inverse_type>::diagonal(unsigned int i) const
+{
+ Assert(store_diagonals, ExcDiagonalsNotStored());
+
+ if (same_diagonal)
+ return var_diagonal[0];
+
+ Assert (i < var_diagonal.size(), ExcIndexRange(i,0,var_diagonal.size()));
+ return var_diagonal[i];
}
void
PreconditionBlock<number,inverse_type>::set_same_diagonal()
{
- Assert(_inverse.size()==0, ExcInverseMatricesAlreadyExist());
+ Assert(var_inverse.size()==0, ExcInverseMatricesAlreadyExist());
same_diagonal = true;
}
bool
PreconditionBlock<number,inverse_type>::inverses_ready() const
{
- return (_inverse.size() != 0);
+ return (var_inverse.size() != 0);
}
Assert (blocksize!=0, ExcNotInitialized());
const SparseMatrix<number> &M=*A;
- Assert (_inverse.size()==0, ExcInverseMatricesAlreadyExist());
+ Assert (var_inverse.size()==0, ExcInverseMatricesAlreadyExist());
const unsigned int n_cells = M.m()/blocksize;
// Invert only the first block
// This is a copy of the code in the
// 'else' part, stripped of the outer loop
- _inverse.resize(1);
+ if (store_diagonals)
+ var_diagonal.resize(1);
+ var_inverse.resize(1);
for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell)
for (unsigned int column_cell=0; column_cell<blocksize; ++column_cell)
M_cell(row_cell,column_cell)=M.el(row_cell,column_cell);
- _inverse[0].invert(M_cell);
+ if (store_diagonals)
+ var_diagonal[0] = M_cell;
+ var_inverse[0].invert(M_cell);
}
else
{
// row, column are the global numbering
// of the unkowns.
- // set the @p{_inverse} array to the right
+ // set the @p{var_inverse} array to the right
// size. we could do it like this:
- // _inverse = vector<>(n_cells,FullMatrix<>())
+ // var_inverse = vector<>(n_cells,FullMatrix<>())
// but this would involve copying many
// FullMatrix objects.
//
// the following is a neat trick which
// avoids copying
+ if (store_diagonals)
+ {
+ vector<FullMatrix<inverse_type> > tmp(n_cells,
+ FullMatrix<inverse_type>(blocksize));
+ var_diagonal.swap (tmp);
+ }
+
vector<FullMatrix<inverse_type> > tmp(n_cells,
FullMatrix<inverse_type>(blocksize));
- _inverse.swap (tmp);
+ var_inverse.swap (tmp);
M_cell.clear ();
column_cell<blocksize; ++column_cell, ++column)
M_cell(row_cell,column_cell)=M.el(row,column);
- _inverse[cell].invert(M_cell);
+ if (store_diagonals)
+ var_diagonal[cell] = M_cell;
+ var_inverse[cell].invert(M_cell);
}
}
}
begin_diag_block+=blocksize;
}
+ dst.scale(relaxation);
}
/*--------------------- PreconditionBlockSOR -----------------------*/
-template<typename number, typename inverse_type>
-PreconditionBlockSOR<number,inverse_type>::PreconditionBlockSOR(const number omega):
- omega(omega)
-{}
-
-
-template<typename number, typename inverse_type>
-PreconditionBlockSOR<number,inverse_type>::~PreconditionBlockSOR()
-{}
-
-
-
-template<typename number, typename inverse_type>
-void
-PreconditionBlockSOR<number,inverse_type>::set_omega(number om)
-{
- omega = om;
-}
-
-
-
template <typename number, typename inverse_type>
template <typename number2>
void PreconditionBlockSOR<number,inverse_type>::vmult (Vector<number2> &dst,
M_cell.backward(x_cell,b_cell);
// distribute x_cell to dst
for (row=cell*blocksize, row_cell=0; row_cell<blocksize; ++row_cell, ++row)
- dst(row)=omega*x_cell(row_cell);
+ dst(row)=relaxation*x_cell(row_cell);
begin_diag_block+=blocksize;
}
inverse(cell).vmult(x_cell, b_cell);
// distribute x_cell to dst
for (row=cell*blocksize, row_cell=0; row_cell<blocksize; ++row_cell, ++row)
- dst(row)=omega*x_cell(row_cell);
+ dst(row)=relaxation*x_cell(row_cell);
begin_diag_block+=blocksize;
}
M_cell.backward(x_cell,b_cell);
// distribute x_cell to dst
for (row=cell*blocksize, row_cell=0; row_cell<blocksize; ++row_cell, ++row)
- dst(row)=omega*x_cell(row_cell);
+ dst(row)=relaxation*x_cell(row_cell);
end_diag_block-=blocksize;
}
inverse(cell).vmult(x_cell, b_cell);
// distribute x_cell to dst
for (row=cell*blocksize, row_cell=0; row_cell<blocksize; ++row_cell, ++row)
- dst(row)=omega*x_cell(row_cell);
+ dst(row)=relaxation*x_cell(row_cell);
end_diag_block-=blocksize;
}
}
+//----------------------------------------------------------------------//
+
+
+template <typename number, typename inverse_type>
+PreconditionBlockSSOR<number,inverse_type>::PreconditionBlockSSOR ()
+{
+ store_diagonals = 1;
+}
+
+
+template <typename number, typename inverse_type>
+template <typename number2>
+void PreconditionBlockSSOR<number,inverse_type>::vmult (Vector<number2> &dst,
+ const Vector<number2> &src) const
+{
+ Vector<number2> help;
+ help.reinit(dst);
+
+ PreconditionBlockSOR<number,inverse_type>::vmult(help, src);
+
+ Vector<inverse_type> cell_src(blocksize);
+ Vector<inverse_type> cell_dst(blocksize);
+
+ const unsigned int n_cells = A->m()/blocksize;
+
+ // Multiply with diagonal blocks
+ for (unsigned int cell=0; cell<n_cells; ++cell)
+ {
+ unsigned int row = cell*blocksize;
+
+ for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell)
+ cell_src(row_cell)=help(row+row_cell);
+
+ diagonal(cell).vmult(cell_dst, cell_src);
+
+ for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell)
+ help(row+row_cell)=relaxation*cell_dst(row_cell);
+ }
+
+ PreconditionBlockSOR<number,inverse_type>::Tvmult(dst, help);
+}
+
+template <typename number, typename inverse_type>
+template <typename number2>
+void PreconditionBlockSSOR<number,inverse_type>::Tvmult (Vector<number2> &dst,
+ const Vector<number2> &src) const
+{
+ Vector<number2> help;
+ help.reinit(dst);
+
+ PreconditionBlockSOR<number,inverse_type>::Tvmult(help, src);
+
+ Vector<inverse_type> cell_src(blocksize);
+ Vector<inverse_type> cell_dst(blocksize);
+
+ const unsigned int n_cells = A->m()/blocksize;
+
+ // Multiply with diagonal blocks
+ for (unsigned int cell=0; cell<n_cells; ++cell)
+ {
+ unsigned int row = cell*blocksize;
+
+ for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell)
+ cell_src(row_cell)=help(row+row_cell);
+
+ diagonal(cell).vmult(cell_dst, cell_src);
+
+ for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell)
+ help(row+row_cell)=relaxation*cell_dst(row_cell);
+ }
+
+ PreconditionBlockSOR<number,inverse_type>::vmult(dst, help);
+}
#endif
Vector<double> &, const Vector<double> &) const;
+/*--------------------- PreconditionBlockSSOR -----------------------*/
+
+
+// explicit instantiations for "float" PreconditionBlock
+template class PreconditionBlockSSOR<float, float>;
+
+template void PreconditionBlockSSOR<float, float>::vmult (
+ Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockSSOR<float, float>::vmult (
+ Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockSSOR<float, float>::Tvmult (
+ Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockSSOR<float, float>::Tvmult (
+ Vector<double> &, const Vector<double> &) const;
+
+
+// the instantiation for class PreconditionBlockSSOR<float, double> is skipped
+// because it does not make sense to have inverse block matrices with
+// higher precision than the matrix itself
+
+
+// explicit instantiations for "double" PreconditionBlockSSOR
+template class PreconditionBlockSSOR<double, float>;
+
+
+template void PreconditionBlockSSOR<double, float>::vmult (
+ Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockSSOR<double, float>::vmult (
+ Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockSSOR<double, float>::Tvmult (
+ Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockSSOR<double, float>::Tvmult (
+ Vector<double> &, const Vector<double> &) const;
+
+template class PreconditionBlockSSOR<double, double>;
+
+template void PreconditionBlockSSOR<double, double>::vmult (
+ Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockSSOR<double, double>::vmult (
+ Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockSSOR<double, double>::Tvmult (
+ Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockSSOR<double, double>::Tvmult (
+ Vector<double> &, const Vector<double> &) const;
+
+
projects / dealii-svn.git / commitdiff