From 3cf7ed19d707eb4d2a00539fefdcbd727584f1b4 Mon Sep 17 00:00:00 2001 From: richter Date: Wed, 8 Mar 2000 14:03:49 +0000 Subject: [PATCH] Add minres git-svn-id: https://svn.dealii.org/trunk@2567 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/solver_minres.h | 310 ++++++++++++++++++++++++ 1 file changed, 310 insertions(+) create mode 100644 deal.II/lac/include/lac/solver_minres.h diff --git a/deal.II/lac/include/lac/solver_minres.h b/deal.II/lac/include/lac/solver_minres.h new file mode 100644 index 0000000000..2e50c305a2 --- /dev/null +++ b/deal.II/lac/include/lac/solver_minres.h @@ -0,0 +1,310 @@ +//---------------------------- solver_minres.h --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000 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. +// +//---------------------------- solver_minres.h ----------------------- +#ifndef __deal2__solver_minres_h +#define __deal2__solver_minres_h + + +#include +#include +#include +#include + + + +/** + * Preconditioned MinRes method. + * + * Like all other solver classes, this class has a local structure called + * #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 #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. + * + * However, since the MinRes method does not need additional data, the respective + * structure is empty and does not offer any functionality. The constructor + * has a default argument, so you may call it without the additional + * parameter. + * + * The preconditioner has to be positive definite and symmetric + * + * @author Thomas Richter, 2000 + */ +template , class Vector = Vector > +class SolverMinRes : public Solver +{ + public: + /** + * Standardized data struct to + * pipe additional data to the + * solver. This solver does not + * need additional data yet. + */ + struct AdditionalData {}; + + /** + * Constructor. + */ + SolverMinRes (SolverControl &cn, + VectorMemory &mem, + const AdditionalData &data=AdditionalData()); + + /** + * Solver method. + */ + template + typename Solver::ReturnState + solve (const Matrix &A, + Vector &x, + const Vector &b, + const Preconditioner& precondition); + + /** + * Exception + */ + DeclException0 (ExcPreconditionerNotDefinite); + + + protected: + /** + * Implementation of the computation of + * the norm of the residual. + */ + virtual long double criterion(); + + /** + * Temporary vectors, allocated through + * the #VectorMemory# object at the start + * of the actual solution process and + * deallocated at the end. + */ + Vector *Vu0, *Vu1, *Vu2; + Vector *Vm0, *Vm1, *Vm2; + Vector *Vv; + + /** + * Within the iteration loop, the + * square of the residual vector is + * stored in this variable. The + * function #criterion# uses this + * variable to compute the convergence + * value, which in this class is the + * norm of the residual vector and thus + * the square root of the #res2# value. + */ + long double res2; +}; + + +/*------------------------- Implementation ----------------------------*/ + + +template +SolverMinRes::SolverMinRes(SolverControl &cn, + VectorMemory &mem, + const AdditionalData &) : + Solver(cn,mem) {}; + + +template +long double +SolverMinRes::criterion() +{ + return res2; +}; + + +template +template +typename Solver::ReturnState +SolverMinRes::solve (const Matrix &A, + Vector &x, + const Vector &b, + const Preconditioner& precondition) +{ + SolverControl::State conv=SolverControl::iterate; + + deallog.push("minres"); + + + unsigned int VS = b.size(); + + + // Memory allocation + Vu0 = memory.alloc(); + Vu1 = memory.alloc(); + Vu2 = memory.alloc(); + Vv = memory.alloc(); + Vm0 = memory.alloc(); + Vm1 = memory.alloc(); + Vm2 = memory.alloc(); + // define some aliases for simpler access + typedef Vector vecref; + vecref u[3] = {*Vu0, *Vu1, *Vu2}; + vecref m[3] = {*Vm0, *Vm1, *Vm2}; + vecref v = *Vv; + // resize the vectors, but do not set + // the values since they'd be overwritten + // soon anyway. + u[0].reinit(VS,true); + u[1].reinit(VS,true); + u[2].reinit(VS,true); + m[0].reinit(VS,true); + m[1].reinit(VS,true); + m[2].reinit(VS,true); + v.reinit(VS,true); + + // some values needed + vector delta(3); + vector f(2); + vector e(2); + + double r_l2 = 0; + double r0 = 0; + double tau = 0; + double c = 0; + double gamma = 0; + double s = 0; + double d_ = 0; + double d = 0; + + // The iteration step. + int j = 1; + + + // Start of the solving process + + // The algorithm is taken from + // Astrid Battermann, Master thesis + // with some changes + + A.vmult(m[0],x); + u[1] = b; + u[1].add(-1.,m[0]); + // Precondition is applied. + // The preconditioner has to be + // positiv definite and symmetric + + // M v = u[1] + precondition (v,u[1]); + + delta[1] = v * u[1]; + Assert (delta[1]>=0, ExcPreconditionerNotDefinite()); + + r0 = sqrt(delta[1]); + r_l2 = r0; + + + u[0].reinit(VS,0); + delta[0] = 1.; + m[0].reinit(VS,0); + m[1].reinit(VS,0); + m[2].reinit(VS,0); + + conv = control().check(0,r_l2); + + while (conv==SolverControl::iterate) + { + + if (delta[1]!=0) + v.scale(1./sqrt(delta[1])); + else + v.reinit(VS,0); + + A.vmult(u[2],v); + u[2].add(-sqrt(delta[1]/delta[0]),u[0]); + + gamma = u[2] * v; + u[2].add (-gamma / sqrt(delta[1]), u[1]); + m[0] = v; + + // precondition: solve M v = u[2] + // Preconditioner has to be positiv + // definite and symmetric. + precondition(v,u[2]); + + delta[2] = v * u[2]; + + Assert (delta[2]>=0, ExcPreconditionerNotDefinite()); + + if (j==1) + { + d_ = gamma; + e[1] = sqrt(delta[2]); + } + if (j>1) + { + d_ = s * e[0] - c * gamma; + e[0] = c * e[0] + s * gamma; + f[1] = s * sqrt(delta[2]); + e[1] = -c * sqrt(delta[2]); + } + + d = sqrt (d_*d_ + delta[2]); + + if (j>1) tau *= s / c; + c = d_ / d; + tau *= c; + + s = sqrt(delta[2]) / d; + + if (j==1) + tau = r0 * c; + + m[0].add (-e[0], m[1]); + if (j>1) + m[0].add (-f[0],m[2]); + m[0].scale(1./d); + x.add (tau, m[0]); + r_l2 *= fabs(s); + + conv = control().check(j,r_l2); + + // next iteration step + ++j; + // All vectors have to be shifted + // one iteration step. + // This should be changed one time. + m[2] = m[1]; + m[1] = m[0]; + + u[0] = u[1]; + u[1] = u[2]; + f[0] = f[1]; + e[0] = e[1]; + delta[0] = delta[1]; + delta[1] = delta[2]; + } + + // Deallocation of Memory + memory.free(Vu0); + memory.free(Vu1); + memory.free(Vu2); + memory.free(Vv); + memory.free(Vm0); + memory.free(Vm1); + memory.free(Vm2); + // Output + deallog.pop (); + + if (conv == SolverControl::failure) + return exceeded; + + return success; +}; + + +#endif + -- 2.39.5