From: Guido Kanschat Date: Fri, 2 Jul 2010 18:17:17 +0000 (+0000) Subject: IterativeInverse replaces PreconditionLACSolver X-Git-Tag: v8.0.0~5856 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e1aa1666434971c7dea5ddc60e7d919f4f5fd0d1;p=dealii.git IterativeInverse replaces PreconditionLACSolver git-svn-id: https://svn.dealii.org/trunk@21447 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/iterative_inverse.h b/deal.II/lac/include/lac/iterative_inverse.h new file mode 100644 index 0000000000..0bf78ff0c9 --- /dev/null +++ b/deal.II/lac/include/lac/iterative_inverse.h @@ -0,0 +1,183 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2009, 2010 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- +#ifndef __deal2__iterative_inverse_h +#define __deal2__iterative_inverse_h + +#include +#include +//#include +//#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + + +/** + * Implementation of the inverse of a matrix, using an iterative + * method. + * + * The function vmult() of this class starts an iterative solver in + * order to approximate the action of the inverse matrix. + * + * Krylov space methods like SolverCG or SolverBicgstab + * become inefficient if soution down to machine accuracy is + * needed. This is due to the fact, that round-off errors spoil the + * orthogonality of the vector sequences. Therefore, a nested + * iteration of two methods is proposed: The outer method is + * SolverRichardson, since it is robust with respect to round-of + * errors. The inner loop is an appropriate Krylov space method, since + * it is fast. + * + * @code + * // Declare related objects + * + * SparseMatrix A; + * Vector x; + * Vector b; + * GrowingVectorMemory > mem; + + * ReductionControl inner_control (10, 1.e-30, 1.e-2) + * PreconditionSSOR > inner_precondition; + * inner_precondition.initialize (A, 1.2); + * + * IterativeInverse > precondition; + * precondition.initialize (A, inner_precondition); + * precondition.solver.select("cg"); + * precondition.solver.control = inner_control; + * + * SolverControl outer_control(100, 1.e-16); + * SolverRichardson > outer_iteration; + * + * outer_iteration.solve (A, x, b, precondition); + * @endcode + * + * Each time we call the inner loop, reduction of the residual by a + * factor 1.e-2 is attempted. Since the right hand side vector of + * the inner iteration is the residual of the outer loop, the relative + * errors are far from machine accuracy, even if the errors of the + * outer loop are in the range of machine accuracy. + * + * @ingroup Matrix2 + * @author Guido Kanschat + * @date 2010 + */ +template +class IterativeInverse : public Subscriptor +{ + public: + /** + * Initialization + * function. Provide a matrix and + * preconditioner for th solve in + * vmult(). + */ + template + void initialize (const MATRIX&, const PRECONDITION&); + + /** + * Delete the pointers to matrix + * and preconditioner. + */ + void clear(); + + /** + * Solve for right hand side src. + */ + void vmult (VECTOR& dst, const VECTOR& src) const; + + /** + * Solve the transposed system + * for right hand side + * src. + */ +// void Tvmult (VECTOR& dst, const VECTOR& src) const; + + /** + * The solver, which allows + * selection of the actual solver + * as well as adjuxtment of + * parameters. + */ + SolverSelector solver; + + private: + /** + * The matrix in use. + */ + std_cxx1x::shared_ptr > matrix; + + /** + * The preconditioner to use. + */ + std_cxx1x::shared_ptr > preconditioner; + /** + * The transpose of the matrix in use. + */ + std_cxx1x::shared_ptr > transpose_matrix; + + /** + * The transpose of the preconditioner to use. + */ + std_cxx1x::shared_ptr > transpose_preconditioner; +}; + + +template +template +inline +void +IterativeInverse::initialize(const MATRIX& m, const PRECONDITION& p) +{ + // dummy variable + VECTOR* v; + matrix = std_cxx1x::shared_ptr > (new_pointer_matrix_base(m, *v)); + preconditioner = std_cxx1x::shared_ptr > (new_pointer_matrix_base(p, *v)); +} + + +template +inline +void +IterativeInverse::clear() +{ + matrix = 0; + preconditioner = 0; +} + + +template +inline void +IterativeInverse::vmult (VECTOR& dst, const VECTOR& src) const +{ + Assert(matrix != 0, ExcNotInitialized()); + Assert(preconditioner != 0, ExcNotInitialized()); + solver.solve(*matrix, dst, src, *preconditioner); +} + + +// template +// inline void +// IterativeInverse::Tvmult (VECTOR& dst, const VECTOR& src) const +// { +// TransposeMatrix< +// solver.solve(*matrix, dst, src, *preconditioner); +// } + + +DEAL_II_NAMESPACE_CLOSE + +#endif + + diff --git a/deal.II/lac/include/lac/precondition.h b/deal.II/lac/include/lac/precondition.h index 62e532ad5e..8aa289fd24 100644 --- a/deal.II/lac/include/lac/precondition.h +++ b/deal.II/lac/include/lac/precondition.h @@ -710,6 +710,9 @@ class PreconditionPSOR : public PreconditionRelaxation /** + * @deprecated This class has been superseded by IterativeInverse, + * which is more flexible and easier to use. + * * Preconditioner using an iterative solver. This preconditioner uses * a fully initialized LAC iterative solver for the approximate * inverse of the matrix. Naturally, this solver needs another