* This is the operator used by
* LAC iterative solvers.
*/
- void operator() (VECTOR &dst,
- const VECTOR &src) const;
+ void vmult (VECTOR &dst,
+ const VECTOR &src) const;
private:
/**
template<class MG, class VECTOR>
void
-PreconditionMG<MG,VECTOR>::operator() (VECTOR &dst,
- const VECTOR &src) const
+PreconditionMG<MG,VECTOR>::vmult (VECTOR &dst,
+ const VECTOR &src) const
{
multigrid->copy_to_mg(src);
multigrid->vcycle(*pre, *post, *coarse);
PrimitiveVectorMemory<> memory;
SolverCG<> cg(control,memory);
- PreconditionRelaxation<>
- prec(system_matrix,
- &SparseMatrix<double>::template precondition_SSOR<double>,
- 1.2);
+ PreconditionSSOR<> prec;
+ prec.initialize (system_matrix, 1.2);
// solve
cg.solve (system_matrix, solution, right_hand_side, prec);
PrimitiveVectorMemory<> memory;
SolverCG<> cg(control,memory);
- PreconditionRelaxation<>
- prec(mass_matrix,
- &SparseMatrix<double>::template precondition_SSOR<double>,
- 1.2);
+ PreconditionSSOR<> prec;
+ prec.initialize(mass_matrix, 1.2);
// solve
cg.solve (mass_matrix, vec, tmp, prec);
PrimitiveVectorMemory<> memory;
SolverCG<> cg(control,memory);
- PreconditionRelaxation<>
- prec(mass_matrix,
- &SparseMatrix<double>::template precondition_SSOR<double>,
- 1.2);
+ PreconditionSSOR<> prec;
+ prec.initialize(mass_matrix, 1.2);
// solve
cg.solve (mass_matrix, boundary_projection, rhs, prec);
}
memory.free(Vy);
+ memory.free(Vr);
deallog.pop();
// Output
}
memory.free(Vy);
-
+ memory.free(Vr);
+
deallog.pop();
// Output
/**
* Exception
*/
- DeclException2 (ExcInvalidIndex,
- int, int,
- << "The given index " << arg1
- << " should be less than " << arg2 << ".");
- /**
- * Exception
- */
DeclException2 (ExcDimensionMismatch,
int, int,
<< "The two dimensions " << arg1 << " and " << arg2
inline number
FullMatrix<number>::operator() (const unsigned int i, const unsigned int j) const
{
- Assert (i<dim_image, ExcInvalidIndex (i, dim_image));
- Assert (j<dim_range, ExcInvalidIndex (j, dim_range));
+ Assert (i<dim_image, ExcIndexRange (i, 0, dim_image));
+ Assert (j<dim_range, ExcIndexRange (j, 0, dim_range));
return el(i,j);
};
inline number &
FullMatrix<number>::operator() (const unsigned int i, const unsigned int j)
{
- Assert (i<dim_image, ExcInvalidIndex (i, dim_image));
- Assert (j<dim_range, ExcInvalidIndex (j, dim_range));
+ Assert (i<dim_image, ExcIndexRange (i, 0, dim_image));
+ Assert (j<dim_range, ExcIndexRange (j, 0, dim_range));
return el(i,j);
};
#include <iomanip>
#include <algorithm>
-
template <typename number>
FullMatrix<number>::FullMatrix (const unsigned int n) :
val (0),
*/
template<class VECTOR>
void Tvmult (VECTOR&, const VECTOR&) const;
- /**
- * Apply preconditioner.
- */
- template<class VECTOR>
- void operator () (VECTOR&, const VECTOR&) const;
-};
-
-
-/**
- * Jacobi preconditioner using matrix built-in function. The MATRIX
- * class used is required to have a function
- * @p{precondition_Jacobi(VECTOR&, const VECTOR&, double}
- *
- * @author Guido Kanschat, 2000
- */
-template <class MATRIX = SparseMatrix<double> >
-class PreconditionJacobi
-{
- public:
- /**
- * Initialize matrix and
- * relaxation parameter. The
- * matrix is just stored in the
- * preconditioner object. The
- * relaxation parameter should be
- * larger than zero and smaller
- * than 2 for numerical
- * reasons. It defaults to 1.
- */
- void initialize (const MATRIX& A, const double omega = 1.);
- /**
- * Apply preconditioner.
- */
- template<class VECTOR>
- void vmult (VECTOR&, const VECTOR&) const;
- /**
- * Apply transpose
- * preconditioner. Since this is
- * a symmetric preconditioner,
- * this function is the same as
- * @ref{vmult}.
- */
- template<class VECTOR>
- void Tvmult (VECTOR&, const VECTOR&) const;
-
- private:
- /**
- * Pointer to the matrix object.
- */
- SmartPointer<const MATRIX> A;
- /**
- * Relaxation parameter.
- */
- double omega;
};
* of this object with the two
* arguments given here.
*/
- void operator() (VECTOR &dst,
- const VECTOR &src) const;
+ void vmult (VECTOR &dst,
+ const VECTOR &src) const;
private:
/**
/**
- * Preconditioner for builtin relaxation methods.
- * Application of this preconditioner includes
- * use of the #precondition_...# methods of #SparseMatrix#.
+ * Base class for other preconditioners.
+ * Here, only some common features Jacobi, SOR and SSOR preconditioners
+ * are implemented. For preconditioning, refer to derived classes.
*
- * \subsection{Use}
- * You will usually not want to create a named object of this type,
- * although possible. The most common use is like this:
- * \begin{verbatim}
- * SolverGMRES<SparseMatrix<double>,
- * Vector<double> > gmres(control,memory,500);
- *
- * gmres.solve (matrix, solution, right_hand_side,
- * PreconditionRelaxation<SparseMatrix<double>,Vector<double> >
- * (matrix,&SparseMatrix<double>::template precondition_Jacobi,
- * 0.5));
- * \end{verbatim}
- * This creates an unnamed object to be passed as the fourth parameter to
- * the solver function of the #SolverGMRES# class. It assumes that the
- * #SparseMatrix# class has a function #precondition_Jacobi# taking two
- * vectors (source and destination) and a relaxation value as parameters. (Unlike
- * for the #PreconditionUseMatrix# class, this time it should work, with
- * relaxation parameter $0.5$.)
- *
- * Note that due to the default template parameters, the above example
- * could be written shorter as follows:
- * \begin{verbatim}
- * ...
- * gmres.solve (matrix, solution, right_hand_side,
- * PreconditionRelaxation<>
- * (matrix,
- * &SparseMatrix<double>::template precondition_Jacobi,
- * 0.5));
- * \end{verbatim}
- *
- * @author Guido Kanschat, Wolfgang Bangerth, 1999
+ * @author Guido Kanschat, 2000
*/
-template<class MATRIX = SparseMatrix<double>, class VECTOR = Vector<double> >
+template<class MATRIX = SparseMatrix<double> >
class PreconditionRelaxation
{
public:
/**
- * Type of the preconditioning
- * function of the matrix.
+ * Initialize matrix and
+ * relaxation parameter. The
+ * matrix is just stored in the
+ * preconditioner object. The
+ * relaxation parameter should be
+ * larger than zero and smaller
+ * than 2 for numerical
+ * reasons. It defaults to 1.
*/
- typedef void ( MATRIX::* function_ptr)(VECTOR&, const VECTOR&,
- typename MATRIX::value_type) const;
-
+ void initialize (const MATRIX& A, const double omega = 1.);
+ protected:
/**
- * Constructor.
- * This constructor stores a
- * reference to the matrix object
- * for later use and selects a
- * preconditioning method, which
- * must be a member function of
- * that matrix.
+ * Pointer to the matrix object.
*/
- PreconditionRelaxation(const MATRIX &M,
- const function_ptr method,
- const double omega = 1.);
-
+ SmartPointer<const MATRIX> A;
/**
- * Execute preconditioning. Calls the
- * function passed to the constructor
- * of this object with the two
- * arguments given here, and the
- * relaxation parameter passed to the
- * constructor.
+ * Relaxation parameter.
*/
- void operator() (VECTOR&, const VECTOR&) const;
+ double omega;
+};
- private:
+/**
+ * Jacobi preconditioner using matrix built-in function. The MATRIX
+ * class used is required to have a function
+ * @p{precondition_Jacobi(VECTOR&, const VECTOR&, double}
+ *
+ * @author Guido Kanschat, 2000
+ */
+template <class MATRIX = SparseMatrix<double> >
+class PreconditionJacobi : public PreconditionRelaxation<MATRIX>
+{
+ public:
/**
- * Pointer to the matrix in use.
+ * Apply preconditioner.
*/
- const MATRIX& matrix;
-
+ template<class VECTOR>
+ void vmult (VECTOR&, const VECTOR&) const;
/**
- * Pointer to the preconditioning
- * function.
+ * Apply transpose
+ * preconditioner. Since this is
+ * a symmetric preconditioner,
+ * this function is the same as
+ * @ref{vmult}.
*/
- const function_ptr precondition;
-
+ template<class VECTOR>
+ void Tvmult (VECTOR&, const VECTOR&) const;
+};
+
+
+/**
+ * SOR preconditioner using matrix built-in function. The MATRIX
+ * class used is required to have functions
+ * @p{precondition_SOR(VECTOR&, const VECTOR&, double)} and
+ * @p{precondition_TSOR(VECTOR&, const VECTOR&, double)}.
+ *
+ * @author Guido Kanschat, 2000
+ */
+template <class MATRIX = SparseMatrix<double> >
+class PreconditionSOR : public PreconditionRelaxation<MATRIX>
+{
+ public:
/**
- * Relaxation parameter.
+ * Apply preconditioner.
*/
- double omega;
+ template<class VECTOR>
+ void vmult (VECTOR&, const VECTOR&) const;
+ /**
+ * Apply transpose
+ * preconditioner.
+ */
+ template<class VECTOR>
+ void Tvmult (VECTOR&, const VECTOR&) const;
+};
+
+
+/**
+ * SSOR preconditioner using matrix built-in function. The MATRIX
+ * class used is required to have a function
+ * @p{precondition_SSOR(VECTOR&, const VECTOR&, double}
+ *
+ * @author Guido Kanschat, 2000
+ */
+template <class MATRIX = SparseMatrix<double> >
+class PreconditionSSOR : public PreconditionRelaxation<MATRIX>
+{
+ public:
+ /**
+ * Apply preconditioner.
+ */
+ template<class VECTOR>
+ void vmult (VECTOR&, const VECTOR&) const;
+ /**
+ * Apply transpose
+ * preconditioner. Since this is
+ * a symmetric preconditioner,
+ * this function is the same as
+ * @ref{vmult}.
+ */
+ template<class VECTOR>
+ void Tvmult (VECTOR&, const VECTOR&) const;
};
+
/**
* Preconditioner using an iterative solver. This preconditioner uses
* a fully initialized LAC iterative solver for the approximate
* Execute preconditioning.
*/
template<class VECTOR>
- void operator() (VECTOR&, const VECTOR&) const;
+ void vmult (VECTOR&, const VECTOR&) const;
private:
/**
VectorMemory<VECTOR>& mem);
/**
- * Preconditioned matrix-vector-product.
+ * Preconditioned
+ * matrix-vector-product.
*/
void vmult (VECTOR& dst, const VECTOR& src) const;
+ /**
+ * Transposed preconditioned
+ * matrix-vector-product.
+ */
+ void Tvmult (VECTOR& dst, const VECTOR& src) const;
+
/**
* Residual $b-PAx$.
*/
/* ---------------------------------- Inline functions ------------------- */
-template<class VECTOR>
-inline void
-PreconditionIdentity::operator () (VECTOR& dst, const VECTOR& src) const
-{
- dst = src;
-}
-
template<class VECTOR>
inline void
PreconditionIdentity::vmult (VECTOR& dst, const VECTOR& src) const
template <class MATRIX>
inline void
-PreconditionJacobi<MATRIX>::initialize (const MATRIX& rA, double o)
+PreconditionRelaxation<MATRIX>::initialize (const MATRIX& rA, double o)
{
A = &rA;
omega = o;
}
+//----------------------------------------------------------------------//
+
template <class MATRIX>
template<class VECTOR>
inline void
PreconditionJacobi<MATRIX>::vmult (VECTOR& dst, const VECTOR& src) const
{
-//TODO: Assert object not initialized
+ Assert (A!=0, ExcNotInitialized());
+ A->precondition_Jacobi (dst, src, omega);
+}
+
+
+template <class MATRIX>
+template<class VECTOR>
+inline void
+PreconditionJacobi<MATRIX>::Tvmult (VECTOR& dst, const VECTOR& src) const
+{
+ Assert (A!=0, ExcNotInitialized());
A->precondition_Jacobi (dst, src, omega);
}
//----------------------------------------------------------------------//
+template <class MATRIX>
+template<class VECTOR>
+inline void
+PreconditionSOR<MATRIX>::vmult (VECTOR& dst, const VECTOR& src) const
+{
+ Assert (A!=0, ExcNotInitialized());
+ A->precondition_SOR (dst, src, omega);
+}
-template<class MATRIX, class VECTOR>
-PreconditionUseMatrix<MATRIX,VECTOR>::PreconditionUseMatrix(const MATRIX& M,
- function_ptr method)
- :
- matrix(M), precondition(method)
-{}
+template <class MATRIX>
+template<class VECTOR>
+inline void
+PreconditionSOR<MATRIX>::Tvmult (VECTOR& dst, const VECTOR& src) const
+{
+ Assert (A!=0, ExcNotInitialized());
+ A->precondition_TSOR (dst, src, omega);
+}
-template<class MATRIX, class VECTOR>
-void
-PreconditionUseMatrix<MATRIX,VECTOR>::operator() (VECTOR& dst,
- const VECTOR& src) const
+
+//----------------------------------------------------------------------//
+
+template <class MATRIX>
+template<class VECTOR>
+inline void
+PreconditionSSOR<MATRIX>::vmult (VECTOR& dst, const VECTOR& src) const
{
- (matrix.*precondition)(dst, src);
+ Assert (A!=0, ExcNotInitialized());
+ A->precondition_SSOR (dst, src, omega);
+}
+
+
+template <class MATRIX>
+template<class VECTOR>
+inline void
+PreconditionSSOR<MATRIX>::Tvmult (VECTOR& dst, const VECTOR& src) const
+{
+ Assert (A!=0, ExcNotInitialized());
+ A->precondition_SSOR (dst, src, omega);
}
+
+//----------------------------------------------------------------------//
+
+
template<class MATRIX, class VECTOR>
-PreconditionRelaxation<MATRIX,VECTOR>::PreconditionRelaxation(const MATRIX& M,
- function_ptr method,
- double omega)
+PreconditionUseMatrix<MATRIX,VECTOR>::PreconditionUseMatrix(const MATRIX& M,
+ function_ptr method)
:
- matrix(M), precondition(method), omega(omega)
+ matrix(M), precondition(method)
{}
template<class MATRIX, class VECTOR>
void
-PreconditionRelaxation<MATRIX,VECTOR>::operator() (VECTOR& dst,
- const VECTOR& src) const
+PreconditionUseMatrix<MATRIX,VECTOR>::vmult (VECTOR& dst,
+ const VECTOR& src) const
{
- (matrix.*precondition)(dst, src, omega);
+ (matrix.*precondition)(dst, src);
}
//////////////////////////////////////////////////////////////////////
template<class SOLVER, class MATRIX, class PRECONDITION>
template<class VECTOR>
void
-PreconditionLACSolver<SOLVER,MATRIX,PRECONDITION>::operator() (VECTOR& dst,
- const VECTOR& src) const
+PreconditionLACSolver<SOLVER,MATRIX,PRECONDITION>::vmult (VECTOR& dst,
+ const VECTOR& src) const
{
solver.solve(matrix, dst, src, precondition);
}
VECTOR* h = mem.alloc();
h->reinit(src);
A.vmult(*h, src);
- P(dst, *h);
+ P.vmult(dst, *h);
mem.free(h);
}
+
+template<class MATRIX, class PRECOND, class VECTOR>
+inline void
+PreconditionedMatrix<MATRIX, PRECOND, VECTOR>
+::Tvmult (VECTOR& dst,
+ const VECTOR& src) const
+{
+ VECTOR* h = mem.alloc();
+ h->reinit(src);
+ A.Tvmult(*h, src);
+ P.Tvmult(dst, *h);
+ mem.free(h);
+}
+
+
template<class MATRIX, class PRECOND, class VECTOR>
inline double
PreconditionedMatrix<MATRIX, PRECOND, VECTOR>
VECTOR* h = mem.alloc();
h->reinit(src);
A.vmult(*h, src);
- P(dst, *h);
+ P.vmult(dst, *h);
mem.free(h);
dst.sadd(-1.,1.,rhs);
return dst.l2_norm ();
* For all matrices that are empty above and below the diagonal
* blocks (i.e. for all block diagonal matrices) the #BlockJacobi# preconditioner
* is a direct solver. For all matrices that are empty only above the diagonal blocks
- * (e.g. the matrices one gets by the DG method with downstream numbering) the
+ * (e.g. the matrices one gets by the DG method with downstream numbering)
* #BlockSOR# is a direct solver.
*
* This first implementation of the #PreconditionBlock# assumes the
*/
virtual ~PreconditionBlock();
+ /**
+ * Initialize matrix and block
+ * size. We store the matrix and
+ * the block size in the
+ * preconditioner object. In a
+ * second step, the inverses of
+ * the diagonal blocks may be
+ * computed.
+ */
+ void initialize (const SparseMatrix<number>& A,
+ const unsigned int block_size);
+
/**
* Deletes the inverse diagonal block
* matrices if existent, sets the
virtual void clear();
/**
- * Takes the matrix that should be used
- * for the preconditioning. A reference
- * to it is stored within this class,
- * but ownership of the matrix remains
- * with the caller of this function.
+ * Use only the inverse of the
+ * first diagonal block to save
+ * memory and computation time.
+ *
+ * Possible applications:
+ * computing on a cartesian grid,
+ * all diagonal blocks are the
+ * same or all diagonal blocks
+ * are at least similar and
+ * inversion of one of them still
+ * yields a preconditioner.
*/
- void use_matrix(const SparseMatrix<number> &M);
+ void set_same_diagonal ();
- /**
- * Set the right block size before calling
- * #invert_diagblocks# or calling the
- * #operator ()# function
- * of #PreconditionBlockSOR#
- * or #PreconditionBlockJacobi#.
- * If #block_size==1# BlockSOR or
- * BlockJacobi are equal to the
- * standard SOR or Jacobi
- * preconditioner.
- */
- void set_block_size (const unsigned int bsize);
-
/**
* Stores the inverse of
* the diagonal blocks
*/
unsigned int block_size () const;
+ /**
+ * Determine, whether inverses
+ * have been computed.
+ */
+ bool inverses_ready () const;
+
/**
* Exception
*/
*/
DeclException0 (ExcInverseMatricesAlreadyExist);
- /**
- * Exception
- */
- DeclException0 (ExcBlockSizeNotSet);
-
/**
* Exception
*/
DeclException0 (ExcMatrixNotSquare);
-
- /**
- * Exception
- */
- DeclException0 (ExcNoMatrixGivenToUse);
protected:
+ /**
+ * Access to the inverse diagonal blocks.
+ *
+ */
+ const FullMatrix<inverse_type>& inverse (unsigned int i) const;
+
/**
* Size of the blocks. Each diagonal
* block is assumed to be of the
* same size.
*/
unsigned int blocksize;
-
+
+ /**
+ * Pointer to the matrix. Make sure that
+ * the matrix exists as long as this class
+ * needs it, i.e. until calling #invert_diagblocks#,
+ * or (if the inverse matrices should not be
+ * stored) until the last call of the
+ * preconditoining #vmult# function of the
+ * derived classes.
+ */
+ SmartPointer<const SparseMatrix<number> > A;
+
+ private:
/**
* Storage of the inverse matrices of
* the diagonal blocks matrices as
* using #inverse_type=float# saves memory
* in comparison with #inverse_type=double#.
*/
- vector<FullMatrix<inverse_type> > inverse;
+ vector<FullMatrix<inverse_type> > _inverse;
/**
- * Pointer to the matrix. Make sure that
- * the matrix exists as long as this class
- * needs it, i.e. until calling #invert_diagblocks#,
- * or (if the inverse matrices should not be
- * stored) until the last call of the
- * preconditoining #operator()# function of the
- * derived classes.
+ * Flag for diagonal compression.
+ * @see set_same_diagonal()
*/
- SmartPointer<const SparseMatrix<number> > A;
+ bool same_diagonal;
};
public:
/**
* Execute block Jacobi preconditioning.
- * Make sure that the right
- * block size
- * of the matrix is set by
- * #set_block_size#
- * before calling this function.
*
* This function will automatically use the
* inverse matrices if they exist, if not
* matrices in each preconditioning step.
*/
template <typename number2>
- void operator() (Vector<number2>&, const Vector<number2>&) const;
+ void vmult (Vector<number2>&, const Vector<number2>&) const;
+
+ /**
+ * Same as #vmult#, since Jacobi is symmetric.
+ */
+ template <typename number2>
+ void Tvmult (Vector<number2>&, const Vector<number2>&) const;
};
/**
* Execute block SOR preconditioning.
- * Make sure that the right
- * block size
- * of the matrix is set by
- * #set_block_size#
- * before calling this function.
*
* This function will automatically use the
* inverse matrices if they exist, if not
* BlockSOR is a direct solver.
*/
template <typename number2>
- void operator() (Vector<number2>&, const Vector<number2>&) const;
+ void vmult (Vector<number2>&, const Vector<number2>&) const;
private:
/**
#define __deal2__precondition_block_templates_h
+#include <base/exceptions.h>
#include <lac/precondition_block.h>
#include <lac/vector.h>
#include <lac/full_matrix.h>
template <typename number, typename inverse_type>
PreconditionBlock<number,inverse_type>::PreconditionBlock ():
blocksize(0),
- A(0) {};
+ A(0),
+ 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 (_inverse.size()!=0)
+ _inverse.erase(_inverse.begin(), _inverse.end());
}
template <typename number, typename inverse_type>
void PreconditionBlock<number,inverse_type>::clear ()
{
- if (inverse.size()!=0)
- inverse.erase(inverse.begin(), inverse.end());
- blocksize=0;
+ if (_inverse.size()!=0)
+ _inverse.erase(_inverse.begin(), _inverse.end());
+ blocksize = 0;
+ same_diagonal = false;
}
template <typename number, typename inverse_type>
-void PreconditionBlock<number,inverse_type>::use_matrix(
- const SparseMatrix<number> &M)
+void PreconditionBlock<number,inverse_type>
+::initialize (const SparseMatrix<number> &M, unsigned int bsize)
{
+ Assert (M.m() == M.n(), ExcMatrixNotSquare());
A = &M;
+ Assert (bsize>0, ExcIndexRange(bsize, 1, M.m()));
+ Assert (A->m()%bsize==0, ExcWrongBlockSize(bsize, A->m()));
+ blocksize=bsize;
}
template <typename number, typename inverse_type>
-void PreconditionBlock<number,inverse_type>::set_block_size(unsigned int bsize) {
- blocksize=bsize;
+const FullMatrix<inverse_type>&
+PreconditionBlock<number,inverse_type>::inverse(unsigned int i) const
+{
+ if (same_diagonal)
+ return _inverse[0];
+
+ Assert (i < _inverse.size(), ExcIndexRange(i,0,_inverse.size()));
+ return _inverse[i];
}
template <typename number, typename inverse_type>
-unsigned int PreconditionBlock<number,inverse_type>::block_size() const {
+unsigned int PreconditionBlock<number,inverse_type>::block_size() const
+{
return blocksize;
}
+template <typename number, typename inverse_type>
+void
+PreconditionBlock<number,inverse_type>::set_same_diagonal()
+{
+ Assert(_inverse.size()==0, ExcInverseMatricesAlreadyExist());
+ same_diagonal = true;
+}
+
+
+template <typename number, typename inverse_type>
+bool
+PreconditionBlock<number,inverse_type>::inverses_ready() const
+{
+ return (_inverse.size() != 0);
+}
+
+
template <typename number, typename inverse_type>
void PreconditionBlock<number,inverse_type>::invert_diagblocks()
{
- Assert (A!=0, ExcNoMatrixGivenToUse());
- const SparseMatrix<number> &M=*A;
- Assert (M.m() == M.n(), ExcMatrixNotSquare());
- Assert (inverse.size()==0, ExcInverseMatricesAlreadyExist());
+ Assert (A!=0, ExcNotInitialized());
+ Assert (blocksize!=0, ExcNotInitialized());
- Assert (blocksize!=0, ExcBlockSizeNotSet());
- Assert (M.m()%blocksize==0, ExcWrongBlockSize(blocksize, M.m()));
+ const SparseMatrix<number> &M=*A;
+ Assert (_inverse.size()==0, ExcInverseMatricesAlreadyExist());
const unsigned int n_cells = M.m()/blocksize;
- // cell_row, cell_column are the
- // numbering of the blocks (cells).
- // row_cell, column_cell are the local
- // numbering of the unknowns in the
- // blocks.
- // row, column are the global numbering
- // of the unkowns.
-
- // set the #inverse# array to the right
- // size. we could do it like this:
- // inverse = vector<>(n_cells,FullMatrix<>())
- // but this would involve copying many
- // FullMatrix objects.
- //
- // the following is a neat trick which
- // avoids copying
- if (true)
- {
- vector<FullMatrix<inverse_type> > tmp(n_cells,
- FullMatrix<inverse_type>(blocksize));
- inverse.swap (tmp);
- };
-
FullMatrix<inverse_type> M_cell(blocksize);
-
- for (unsigned int cell=0, row=0; cell<n_cells; ++cell)
+
+ if (same_diagonal)
{
- for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell, ++row)
- for (unsigned int column_cell=0, column=cell*blocksize;
- column_cell<blocksize; ++column_cell, ++column)
- M_cell(row_cell,column_cell)=M.el(row,column);
-// try
-// {
+ // Invert only the first block
+ // This is a copy of the code in the
+ // 'else' part, stripped of the outer loop
+ _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);
+
if (blocksize <=4)
{
- inverse[cell].invert(M_cell);
+ _inverse[0].invert(M_cell);
}
else
{
M_cell.gauss_jordan();
- inverse[cell]=M_cell;
+ _inverse[0]=M_cell;
}
+ }
+ else
+ {
+ // cell_row, cell_column are the
+ // numbering of the blocks (cells).
+ // row_cell, column_cell are the local
+ // numbering of the unknowns in the
+ // blocks.
+ // row, column are the global numbering
+ // of the unkowns.
+
+ // set the #_inverse# array to the right
+ // size. we could do it like this:
+ // _inverse = vector<>(n_cells,FullMatrix<>())
+ // but this would involve copying many
+ // FullMatrix objects.
+ //
+ // the following is a neat trick which
+ // avoids copying
+ vector<FullMatrix<inverse_type> > tmp(n_cells,
+ FullMatrix<inverse_type>(blocksize));
+ _inverse.swap (tmp);
+
+ M_cell.clear ();
+
+ for (unsigned int cell=0, row=0; cell<n_cells; ++cell)
+ {
+ for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell, ++row)
+ for (unsigned int column_cell=0, column=cell*blocksize;
+ column_cell<blocksize; ++column_cell, ++column)
+ M_cell(row_cell,column_cell)=M.el(row,column);
+// try
+// {
+ if (blocksize <=4)
+ {
+ _inverse[cell].invert(M_cell);
+ }
+ else
+ {
+ M_cell.gauss_jordan();
+ _inverse[cell]=M_cell;
+ }
// }
// catch (ExcNotImplemented &)
+ }
}
}
template <typename number, typename inverse_type>
template <typename number2>
void PreconditionBlockJacobi<number,inverse_type>
-::operator() (Vector<number2> &dst,
- const Vector<number2> &src) const
+::vmult (Vector<number2> &dst,
+ const Vector<number2> &src) const
{
// introduce the following typedef
// since in the use of exceptions,
// exceptions that do not take
// args...
typedef PreconditionBlock<number,inverse_type> BaseClass;
- Assert(A!=0, typename BaseClass::ExcNoMatrixGivenToUse());
+ Assert(A!=0, ExcNotInitialized());
const SparseMatrix<number> &M=*A;
const unsigned int n_cells=M.m()/blocksize;
- Assert (M.m() == M.n(),
- typename BaseClass::ExcMatrixNotSquare());
- Assert (blocksize!=0,
- typename BaseClass::ExcBlockSizeNotSet());
- Assert (M.m()%blocksize==0,
- typename BaseClass::ExcWrongBlockSize(blocksize, M.m()));
- Assert (inverse.size()==0 || inverse.size()==n_cells,
- typename BaseClass::ExcWrongNumberOfInverses(inverse.size(), n_cells));
-
Vector<number2> b_cell(blocksize), x_cell(blocksize);
// cell_row, cell_column are the
// of the unkowns.
unsigned int row, row_cell, begin_diag_block=0;
- if (inverse.size()==0)
+ if (!inverses_ready())
{
FullMatrix<number> M_cell(blocksize);
for (unsigned int cell=0; cell<n_cells; ++cell)
{
b_cell(row_cell)=src(row);
}
- inverse[cell].vmult(x_cell, b_cell);
+ 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)=x_cell(row_cell);
begin_diag_block+=blocksize;
}
}
+
+
+template <typename number, typename inverse_type>
+template <typename number2>
+void PreconditionBlockJacobi<number,inverse_type>
+::Tvmult (Vector<number2> &dst,
+ const Vector<number2> &src) const
+{
+ vmult(dst, src);
+}
+
+
/*--------------------- PreconditionBlockSOR -----------------------*/
template<typename number, typename inverse_type>
PreconditionBlockSOR<number,inverse_type>::PreconditionBlockSOR(const number omega):
- omega(omega) {}
+ omega(omega)
+{}
template<typename number, typename inverse_type>
-PreconditionBlockSOR<number,inverse_type>::~PreconditionBlockSOR(){}
+PreconditionBlockSOR<number,inverse_type>::~PreconditionBlockSOR()
+{}
}
-
+//TODO: implement Tvmult
template <typename number, typename inverse_type>
template <typename number2>
-void PreconditionBlockSOR<number,inverse_type>::operator() (Vector<number2> &dst,
- const Vector<number2> &src) const
+void PreconditionBlockSOR<number,inverse_type>::vmult (Vector<number2> &dst,
+ const Vector<number2> &src) const
{
// introduce the following typedef
// since in the use of exceptions,
// args...
typedef PreconditionBlock<number,inverse_type> BaseClass;
- Assert (A!=0, typename BaseClass::ExcNoMatrixGivenToUse());
+ Assert (A!=0, ExcNotInitialized());
const SparseMatrix<number> &M=*A;
const unsigned int n_cells=M.m()/blocksize;
- Assert (M.m() == M.n(),
- typename BaseClass::ExcMatrixNotSquare());
- Assert (blocksize!=0,
- typename BaseClass::ExcBlockSizeNotSet());
- Assert (M.m()%blocksize==0,
- typename BaseClass::ExcWrongBlockSize(blocksize, M.m()));
- Assert (inverse.size()==0 || inverse.size()==n_cells,
- typename BaseClass::ExcWrongNumberOfInverses(inverse.size(), n_cells));
-
const SparsityPattern &spars = M.get_sparsity_pattern();
const unsigned int *rowstart = spars.get_rowstart_indices();
const unsigned int *columns = spars.get_column_numbers();
unsigned int row, column, row_cell, begin_diag_block=0;
number2 b_cell_row;
- if (inverse.size()==0)
+ if (!inverses_ready())
{
FullMatrix<number> M_cell(blocksize);
for (unsigned int cell=0; cell<n_cells; ++cell)
}
b_cell(row_cell)=b_cell_row;
}
- inverse[cell].vmult(x_cell, b_cell);
+ 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);
/**
* Solve primal problem only.
*/
- template<class MATRIX, class Preconditioner>
+ template<class MATRIX, class PRECONDITIONER>
typename Solver<VECTOR>::ReturnState
solve (const MATRIX &A,
VECTOR &x,
const VECTOR &b,
- const Preconditioner& precondition);
+ const PRECONDITIONER& precondition);
protected:
/**
beta = rhobar * alpha / (rho * omega);
rho = rhobar;
p.sadd(beta, 1., r, -beta*omega, v);
- precondition(y,p);
+ precondition.vmult(y,p);
A.vmult(v,y);
rhobar = rbar * v;
return typename Solver<VECTOR>::ReturnState(breakdown);
s.equ(1., r, -alpha, v);
- precondition(z,s);
+ precondition.vmult(z,s);
A.vmult(t,z);
rhobar = t*s;
omega = rhobar/(t*t);
/**
* Solver method.
*/
- template<class MATRIX, class Preconditioner>
+ template<class MATRIX, class PRECONDITIONER>
typename Solver<VECTOR>::ReturnState
solve (const MATRIX &A,
VECTOR &x,
const VECTOR &b,
- const Preconditioner& precondition);
+ const PRECONDITIONER& precondition);
protected:
/**
template<class VECTOR>
-template<class MATRIX, class Preconditioner>
+template<class MATRIX, class PRECONDITIONER>
typename Solver<VECTOR>::ReturnState
SolverCG<VECTOR>::solve (const MATRIX &A,
VECTOR &x,
const VECTOR &b,
- const Preconditioner& precondition)
+ const PRECONDITIONER& precondition)
{
SolverControl::State conv=SolverControl::iterate;
};
g.scale(-1.);
- precondition(h,g);
+ precondition.vmult(h,g);
d.equ(-1.,h);
if (conv)
break;
- precondition(h,g);
+ precondition.vmult(h,g);
beta = gh;
gh = g*h;
if (left_precondition)
{
A.residual(p,x,b);
- precondition(v,p);
+ precondition.vmult(v,p);
} else {
A.residual(v,x,b);
};
if (left_precondition)
{
A.vmult(p, *tmp_vectors[inner_iteration]);
- precondition(vv,p);
+ precondition.vmult(vv,p);
} else {
- precondition(p,*tmp_vectors[inner_iteration]);
+ precondition.vmult(p,*tmp_vectors[inner_iteration]);
A.vmult(vv,p);
};
p = 0.;
for (unsigned int i=0; i<dim; ++i)
p.add(h(i), *tmp_vectors[i]);
- precondition(v,p);
+ precondition.vmult(v,p);
x.add(1.,v);
};
/**
* Solver method.
*/
- template<class MATRIX, class Preconditioner>
+ template<class MATRIX, class PRECONDITIONER>
typename Solver<VECTOR>::ReturnState
solve (const MATRIX &A,
VECTOR &x,
const VECTOR &b,
- const Preconditioner& precondition);
+ const PRECONDITIONER& precondition);
/**
* Exception
template<class VECTOR>
-template<class MATRIX, class Preconditioner>
+template<class MATRIX, class PRECONDITIONER>
typename Solver<VECTOR>::ReturnState
SolverMinRes<VECTOR>::solve (const MATRIX &A,
VECTOR &x,
const VECTOR &b,
- const Preconditioner& precondition)
+ const PRECONDITIONER& precondition)
{
SolverControl::State conv=SolverControl::iterate;
// positiv definite and symmetric
// M v = u[1]
- precondition (v,u[1]);
+ precondition.vmult (v,u[1]);
delta[1] = v * u[1];
// Preconditioner positive
// precondition: solve M v = u[2]
// Preconditioner has to be positiv
// definite and symmetric.
- precondition(v,u[2]);
+ precondition.vmult(v,u[2]);
delta[2] = v * u[2];
/**
* Solver method.
*/
- template<class MATRIX, class Preconditioner>
+ template<class MATRIX, class PRECONDITIONER>
typename Solver<VECTOR>::ReturnState
solve (const MATRIX &A,
VECTOR &x,
const VECTOR &b,
- const Preconditioner& precondition);
+ const PRECONDITIONER& precondition);
/**
* Interface for derived class.
/**
* The iteration loop itself.
*/
- template<class MATRIX, class Preconditioner>
+ template<class MATRIX, class PRECONDITIONER>
typename Solver<VECTOR>::ReturnState
- iterate(const MATRIX& A, const Preconditioner& precondition);
+ iterate(const MATRIX& A, const PRECONDITIONER& precondition);
/**
* The current iteration step.
*/
template<class VECTOR>
-template<class MATRIX, class Preconditioner>
+template<class MATRIX, class PRECONDITIONER>
typename Solver<VECTOR>::ReturnState
SolverQMRS<VECTOR>::solve (const MATRIX &A,
VECTOR &x,
const VECTOR &b,
- const Preconditioner& precondition)
+ const PRECONDITIONER& precondition)
{
deallog.push("QMRS");
template<class VECTOR>
-template<class MATRIX, class Preconditioner>
+template<class MATRIX, class PRECONDITIONER>
typename Solver<VECTOR>::ReturnState
SolverQMRS<VECTOR>::iterate(const MATRIX& A,
- const Preconditioner& precondition)
+ const PRECONDITIONER& precondition)
{
/* Remark: the matrix A in the article is the preconditioned matrix.
* Therefore, we have to precondition x before we compute the first residual.
d.reinit(x);
// Apply right preconditioning to x
- precondition(q,x);
+ precondition.vmult(q,x);
// Preconditioned residual
res = A.residual(v, q, b);
p = v;
- precondition(q,p);
+ precondition.vmult(q,p);
tau = v.norm_sqr();
//deallog << "tau:" << tau << endl;
return breakdown;
// Step 7
rho_old = rho;
- precondition(q,v);
+ precondition.vmult(q,v);
rho = q*v;
beta = rho/rho_old;
p.sadd(beta,v);
- precondition(q,p);
+ precondition.vmult(q,p);
}
return exceeded;
}
/**
* Solve $Ax=b$ for $x$.
*/
- template<class MATRIX, class Preconditioner>
+ template<class MATRIX, class PRECONDITIONER>
typename Solver<VECTOR>::ReturnState solve (const MATRIX &A,
VECTOR &x,
const VECTOR &b,
- const Preconditioner& precondition);
+ const PRECONDITIONER& precondition);
/**
* Set the damping-coefficient.
template<class VECTOR>
-template<class MATRIX, class Preconditioner>
+template<class MATRIX, class PRECONDITIONER>
typename Solver<VECTOR>::ReturnState
SolverRichardson<VECTOR>::solve (const MATRIX &A,
VECTOR &x,
const VECTOR &b,
- const Preconditioner& precondition)
+ const PRECONDITIONER& precondition)
{
SolverControl::State conv=SolverControl::iterate;
if (conv != SolverControl::iterate)
break;
- precondition(d,r);
+ precondition.vmult(d,r);
x.add(additional_data.omega,d);
print_vectors(iter,x,r,d);
}
* Same as #apply_decomposition#, format for LAC.
*/
template <typename somenumber>
- void operator() (Vector<somenumber> &dst,
- const Vector<somenumber> &src) const;
+ void vmult (Vector<somenumber> &dst,
+ const Vector<somenumber> &src) const;
/**
* Exception
template <typename number>
template <typename somenumber>
void
-SparseILU<number>::operator() (Vector<somenumber> &dst,
- const Vector<somenumber> &src) const
+SparseILU<number>::vmult (Vector<somenumber> &dst,
+ const Vector<somenumber> &src) const
{
apply_decomposition(dst, src);
}
* update vector in #dst#.
*/
template<typename number2>
- void operator() (Vector<number2> &dst,
- const Vector<number2> &src) const;
+ void vmult (Vector<number2> &dst,
+ const Vector<number2> &src) const;
/**
* Exception
* function with a
* #vector<bool>(n_dofs,true)#.
*
- * The #operator()# of this class
+ * The #vmult# of this class
* of course calls this function
* with a null pointer
*/
* Apply the preconditioner.
*/
template<typename number2>
- void operator() (Vector<number2> &dst,
+ void vmult (Vector<number2> &dst,
const Vector<number2> &src) const;
private:
inverses[row]->gauss_jordan();
};
-
+//TODO: implement Tvmult
template<typename number>
template<typename number2>
void
-SparseVanka<number>::operator ()(Vector<number2> &dst,
- const Vector<number2> &src) const
+SparseVanka<number>::vmult (Vector<number2> &dst,
+ const Vector<number2> &src) const
{
// first set output vector to zero
dst.clear ();
};
};
+//TODO: implement Tvmult
template <typename number>
template <typename number2>
-void SparseBlockVanka<number>::operator() (Vector<number2> &dst,
- const Vector<number2> &src) const
+void SparseBlockVanka<number>::vmult (Vector<number2> &dst,
+ const Vector<number2> &src) const
{
dst.clear ();
// explicit instantiations for "float" PreconditionBlock
template class PreconditionBlockJacobi<float, float>;
-template void PreconditionBlockJacobi<float, float>::operator() (
- Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<float, float>::operator() (
- Vector<double> &, const Vector<double> &) const;
-
+template void PreconditionBlockJacobi<float, float>::vmult
+(Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockJacobi<float, float>::vmult
+(Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockJacobi<float, float>::Tvmult
+(Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockJacobi<float, float>::Tvmult
+(Vector<double> &, const Vector<double> &) const;
template class PreconditionBlockJacobi<double, float>;
-template void PreconditionBlockJacobi<double, float>::operator() (
- Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<double, float>::operator() (
- Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockJacobi<double, float>::vmult
+(Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockJacobi<double, float>::vmult
+(Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockJacobi<double, float>::Tvmult
+(Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockJacobi<double, float>::Tvmult
+(Vector<double> &, const Vector<double> &) const;
template class PreconditionBlockJacobi<double, double>;
-template void PreconditionBlockJacobi<double, double>::operator() (
- Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<double, double>::operator() (
- Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockJacobi<double, double>::vmult
+(Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockJacobi<double, double>::vmult
+(Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockJacobi<double, double>::Tvmult
+(Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockJacobi<double, double>::Tvmult
+(Vector<double> &, const Vector<double> &) const;
/*--------------------- PreconditionBlockGaussSeidel -----------------------*/
// explicit instantiations for "float" PreconditionBlock
template class PreconditionBlockSOR<float, float>;
-template void PreconditionBlockSOR<float, float>::operator() (
+template void PreconditionBlockSOR<float, float>::vmult (
Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSOR<float, float>::operator() (
+template void PreconditionBlockSOR<float, float>::vmult (
Vector<double> &, const Vector<double> &) const;
template class PreconditionBlockSOR<double, float>;
-template void PreconditionBlockSOR<double, float>::operator() (
+template void PreconditionBlockSOR<double, float>::vmult (
Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSOR<double, float>::operator() (
+template void PreconditionBlockSOR<double, float>::vmult (
Vector<double> &, const Vector<double> &) const;
template class PreconditionBlockSOR<double, double>;
-template void PreconditionBlockSOR<double, double>::operator() (
+template void PreconditionBlockSOR<double, double>::vmult (
Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSOR<double, double>::operator() (
+template void PreconditionBlockSOR<double, double>::vmult (
Vector<double> &, const Vector<double> &) const;
template class SparseVanka<float>;
template class SparseVanka<double>;
-template void SparseVanka<double>::operator () (Vector<float> &dst,
- const Vector<float> &src) const;
-template void SparseVanka<double>::operator () (Vector<double> &dst,
- const Vector<double> &src) const;
+template void SparseVanka<double>::vmult (Vector<float> &dst,
+ const Vector<float> &src) const;
+template void SparseVanka<double>::vmult (Vector<double> &dst,
+ const Vector<double> &src) const;
template class SparseBlockVanka<float>;
template class SparseBlockVanka<double>;
-template void SparseBlockVanka<double>::operator () (Vector<float> &dst,
- const Vector<float> &src) const;
-template void SparseBlockVanka<double>::operator () (Vector<double> &dst,
- const Vector<double> &src) const;
+template void SparseBlockVanka<double>::vmult (Vector<float> &dst,
+ const Vector<float> &src) const;
+template void SparseBlockVanka<double>::vmult (Vector<double> &dst,
+ const Vector<double> &src) const;
SolverControl control(2,1.e-300,false,false);
PrimitiveVectorMemory<> mem;
SolverRichardson<> rich(control, mem);
- PreconditionRelaxation<>
- prec((*matrices)[level], &SparseMatrix<double> ::template precondition_SSOR<double>, 1.);
+ PreconditionSSOR<> prec;
+ prec.initialize((*matrices)[level], 1.);
rich.solve((*matrices)[level], u, rhs, prec);
}
PrimitiveVectorMemory<> vector_memory;
SolverCG<> cg (solver_control, vector_memory);
- PreconditionRelaxation<>
- preconditioner(global_system_matrix,
- &SparseMatrix<double>::template precondition_SSOR<double>,
- 1.2);
+ PreconditionSSOR<> preconditioner;
+ preconditioner.initialize(global_system_matrix, 1.2);
solution.clear ();
cg.solve (global_system_matrix, solution, system_rhs,
PrimitiveVectorMemory<Vector> vector_memory;
SolverCG<Vector> cg (solver_control, vector_memory);
- PreconditionRelaxation<Matrix,Vector> preconditioner
- (system_matrix, &Matrix::precondition_Jacobi, 0.8);
+ PreconditionJacobi<Matrix> preconditioner;
+ preconditioner.initialize (system_matrix, 0.8);
cg.solve (system_matrix, solution, system_rhs,
preconditioner);
SolverControl control(2,1.e-300,false,false);
PrimitiveVectorMemory<> mem;
SolverRichardson<> rich(control, mem);
- PreconditionRelaxation<>
- prec((*matrices)[level], &SparseMatrix<double> ::template precondition_SSOR<double>, 1.);
+ PreconditionSSOR<> prec;
+ prec.initialize((*matrices)[level], 1.);
rich.solve((*matrices)[level], u, rhs, prec);
}