From: nauber Date: Tue, 20 Apr 1999 08:39:39 +0000 (+0000) Subject: Implementation of the Gauss-Seidel iterative solver X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8b2f533ace1ac6055ab1213a0e57bf9163414e72;p=dealii-svn.git Implementation of the Gauss-Seidel iterative solver git-svn-id: https://svn.dealii.org/trunk@1184 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/solver_gauss_seidel.h b/deal.II/lac/include/lac/solver_gauss_seidel.h new file mode 100644 index 0000000000..19c8a41c16 --- /dev/null +++ b/deal.II/lac/include/lac/solver_gauss_seidel.h @@ -0,0 +1,136 @@ +/*---------------------------- solver_gauss_seidel.h ---------------------------*/ +/* $Id$ */ +#ifndef __solver_gauss_seidel_H +#define __solver_gauss_seidel_H +/*---------------------------- solver_gauss_seidel.h ---------------------------*/ + +#include +#include +#include + +/** + * Implementation of the Gauss-Seidel method. The stopping criterion + * is the norm of the defect, i.e. if x is the iteration vector, b the + * rhs and A the matrix, the $res = \| A x - b \|$. + */ +template +class SolverGaussSeidel : public Solver { + public: + /** + * Constructor. Takes a SolverControl and some + VectorMemory. VectorMemory is not used. + */ + SolverGaussSeidel (SolverControl &cn, VectorMemory &mem) : + Solver (cn,mem) + {}; + + /** + * Solver method. Just specify the vectors A, x and b and let it + * run. x should contain a reasonable starting value. + */ + virtual ReturnState solve (const Matrix &A, + Vector &x, + const Vector &b); + + protected: + /** + * Implementation of the computation of + * the norm of the residual. + */ + virtual double criterion(); + + /** + * 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. + */ + double res2; +}; + +/*----------------- Implementation of the Gauss-Seidel Method ------------------------*/ + +template +SolverGaussSeidel::ReturnState +SolverGaussSeidel::solve (const Matrix &A, + Vector &x, + const Vector &b) { + + deallog.push(__FILE__); + deallog.push("SolverGaussSeidel::solve"); + + TRACEMSG("Initialize Variables"); + + unsigned int i; + unsigned int n = b.size(); + double res; + + SolverControl::State conv=SolverControl::iterate; + + Vector defect(n); + + TRACEMSG("Entering main loop"); + + deallog.push("Main loop"); + + // Main loop + for(int iter=0; conv==SolverControl::iterate; iter++) + { + // Calculate defect, i.e. defect = A x - b + A.vmult(defect,x); + defect.add(-1,b); + + // Calculate residual + res2 = defect * defect; + + // Apply Gauss-Seidel preconditioner + for (i=0;i +inline double +SolverGaussSeidel::criterion() +{ + return sqrt(res2); +} + +/*---------------------------- solver_gauss_seidel.h ---------------------------*/ +/* end of #ifndef __solver_gauss_seidel_H */ +#endif +/*---------------------------- solver_gauss_seidel.h ---------------------------*/