From 2ea558d1011dd6d9b63fe397a6ab4a72b69eb827 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Tue, 20 Jul 2010 19:17:02 +0000 Subject: [PATCH] add SolverRelaxation and fix (!) SparseMatrix::precondition_SSOR and SparseMatrixEZ::precondition_SSOR git-svn-id: https://svn.dealii.org/trunk@21546 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/solver_relaxation.h | 154 ++++++++++++++++++ deal.II/lac/include/lac/solver_richardson.h | 5 +- .../lac/include/lac/sparse_matrix.templates.h | 5 +- .../include/lac/sparse_matrix_ez.templates.h | 4 +- 4 files changed, 162 insertions(+), 6 deletions(-) create mode 100644 deal.II/lac/include/lac/solver_relaxation.h diff --git a/deal.II/lac/include/lac/solver_relaxation.h b/deal.II/lac/include/lac/solver_relaxation.h new file mode 100644 index 0000000000..5b2153d634 --- /dev/null +++ b/deal.II/lac/include/lac/solver_relaxation.h @@ -0,0 +1,154 @@ +//--------------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 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__solver_relaxation_h +#define __deal2__solver_relaxation_h + + +#include +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +/** + * Implementation of an iterative solver based on relaxation + * methods. The stopping criterion is the norm of the residual. + * + * For the requirements on matrices and vectors in order to work with + * this class, see the documentation of the Solver base class. + * + * Like all other solver classes, this class has a local structure called + * @p AdditionalData which is used to pass additional parameters to the + * solver, like damping parameters or the number of temporary vectors. We + * use this additional structure instead of passing these values directly + * to the constructor because this makes the use of the @p SolverSelector and + * other classes much easier and guarantees that these will continue to + * work even if number or type of the additional parameters for a certain + * solver changes. AdditionalData of this class currently does not + * contain any data. + * + * @author Guido Kanschat + * @date 2010 + */ +template > +class SolverRelaxation : public Solver +{ + public: + /** + * Standardized data struct to + * pipe additional data to the + * solver. There is no data in + * here for relaxation methods. + */ + struct AdditionalData {}; + + /** + * Constructor. + */ + SolverRelaxation (SolverControl &cn, + const AdditionalData &data=AdditionalData()); + + /** + * Virtual destructor. + */ + virtual ~SolverRelaxation (); + + /** + * Solve the system $Ax = b$ + * using the relaxation method + * $x_{k+1} = R(x_k,b)$. The + * amtrix A itself is only + * used to compute the residual. + */ + template + void + solve (const MATRIX& A, + VECTOR& x, + const VECTOR& b, + const RELAXATION& R); +}; + +//----------------------------------------------------------------------// + +template +SolverRelaxation::SolverRelaxation(SolverControl &cn, + const AdditionalData &) + : + Solver (cn) +{} + + + +template +SolverRelaxation::~SolverRelaxation() +{} + + +template +template +void +SolverRelaxation::solve ( + const MATRIX& A, + VECTOR& x, + const VECTOR& b, + const RELAXATION &R) +{ + GrowingVectorMemory mem; + SolverControl::State conv=SolverControl::iterate; + + // Memory allocation + typename VectorMemory::Pointer Vr(mem); VECTOR& r = *Vr; r.reinit(x); + typename VectorMemory::Pointer Vd(mem); VECTOR& d = *Vd; d.reinit(x); + + deallog.push("Relaxation"); + + try + { + // Main loop + for(int iter=0; conv==SolverControl::iterate; iter++) + { + // Compute residual + A.vmult(r,x); + r.sadd(-1.,1.,b); + + // The required norm of the + // (preconditioned) + // residual is computed in + // criterion() and stored + // in res. + conv = this->control().check (iter, r.l2_norm()); + if (conv != SolverControl::iterate) + break; + R.step(x,b); + } + } + catch (...) + { + deallog.pop(); + throw; + } + deallog.pop(); + + // in case of failure: throw + // exception + if (this->control().last_check() != SolverControl::success) + throw SolverControl::NoConvergence (this->control().last_step(), + this->control().last_value()); + // otherwise exit as normal +} + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/deal.II/lac/include/lac/solver_richardson.h b/deal.II/lac/include/lac/solver_richardson.h index 13f5494329..f67afba47a 100644 --- a/deal.II/lac/include/lac/solver_richardson.h +++ b/deal.II/lac/include/lac/solver_richardson.h @@ -1,8 +1,7 @@ //--------------------------------------------------------------------------- // $Id$ -// Version: $Name$ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -26,7 +25,7 @@ DEAL_II_NAMESPACE_OPEN /*@{*/ /** - * Implementation of the richardson iteration method. The stopping criterion + * Implementation of the preconditioned Richardson iteration method. The stopping criterion * is the norm of the residual. * * For the requirements on matrices and vectors in order to work with diff --git a/deal.II/lac/include/lac/sparse_matrix.templates.h b/deal.II/lac/include/lac/sparse_matrix.templates.h index e2c70d60b9..faaf4ea460 100644 --- a/deal.II/lac/include/lac/sparse_matrix.templates.h +++ b/deal.II/lac/include/lac/sparse_matrix.templates.h @@ -1357,7 +1357,7 @@ SparseMatrix::precondition_SSOR (Vector &dst, rowstart_ptr = &cols->rowstart[0]; dst_ptr = &dst(0); for ( ; rowstart_ptr!=&cols->rowstart[n]; ++rowstart_ptr, ++dst_ptr) - *dst_ptr *= (2.-om)*val[*rowstart_ptr]; + *dst_ptr *= om*(2.-om)*val[*rowstart_ptr]; // backward sweep rowstart_ptr = &cols->rowstart[n-1]; @@ -1709,6 +1709,9 @@ void SparseMatrix::SSOR (Vector& dst, const number om) const { +//TODO: Is this called anywhere? If so, multiplication with om(2-om)D is missing + Assert(false, ExcNotImplemented()); + Assert (cols != 0, ExcNotInitialized()); Assert (val != 0, ExcNotInitialized()); Assert (cols->optimize_diagonal(), diff --git a/deal.II/lac/include/lac/sparse_matrix_ez.templates.h b/deal.II/lac/include/lac/sparse_matrix_ez.templates.h index 64f553d672..da15f1528b 100644 --- a/deal.II/lac/include/lac/sparse_matrix_ez.templates.h +++ b/deal.II/lac/include/lac/sparse_matrix_ez.templates.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors +// Copyright (C) 2002, 2003, 2004, 2005, 2006, 2007, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -353,7 +353,7 @@ SparseMatrixEZ::precondition_SSOR (Vector &dst, // Diagonal dst_ptr = dst.begin(); for (ri = row_info.begin(); ri != end; ++dst_ptr, ++ri) - *dst_ptr *= (2.-om) * data[ri->start + ri->diagonal].value; + *dst_ptr *= om*(2.-om) * data[ri->start + ri->diagonal].value; // Backward typename std::vector::const_reverse_iterator rri; -- 2.39.5