]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add SolverRelaxation and fix (!) SparseMatrix::precondition_SSOR and SparseMatrixEZ...
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Tue, 20 Jul 2010 19:17:02 +0000 (19:17 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Tue, 20 Jul 2010 19:17:02 +0000 (19:17 +0000)
git-svn-id: https://svn.dealii.org/trunk@21546 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/solver_relaxation.h [new file with mode: 0644]
deal.II/lac/include/lac/solver_richardson.h
deal.II/lac/include/lac/sparse_matrix.templates.h
deal.II/lac/include/lac/sparse_matrix_ez.templates.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 (file)
index 0000000..5b2153d
--- /dev/null
@@ -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 <base/config.h>
+#include <base/logstream.h>
+#include <lac/solver.h>
+#include <lac/solver_control.h>
+#include <base/subscriptor.h>
+
+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 VECTOR = Vector<double> >
+class SolverRelaxation : public Solver<VECTOR>
+{
+  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 <i>A</i> itself is only
+                                     * used to compute the residual.
+                                     */
+    template<class MATRIX, class RELAXATION>
+    void
+    solve (const MATRIX& A,
+          VECTOR& x,
+          const VECTOR& b,
+          const RELAXATION& R);
+};
+
+//----------------------------------------------------------------------//
+
+template <class VECTOR>
+SolverRelaxation<VECTOR>::SolverRelaxation(SolverControl &cn,
+                                          const AdditionalData &)
+               :
+               Solver<VECTOR> (cn)
+{}
+
+
+
+template <class VECTOR>
+SolverRelaxation<VECTOR>::~SolverRelaxation()
+{}
+
+
+template <class VECTOR>
+template <class MATRIX, class RELAXATION>
+void
+SolverRelaxation<VECTOR>::solve (
+  const MATRIX& A,
+  VECTOR& x,
+  const VECTOR& b,
+  const RELAXATION &R)
+{
+  GrowingVectorMemory<VECTOR> mem;
+  SolverControl::State conv=SolverControl::iterate;
+
+                                  // Memory allocation
+  typename VectorMemory<VECTOR>::Pointer Vr(mem); VECTOR& r  = *Vr; r.reinit(x);
+  typename VectorMemory<VECTOR>::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
index 13f5494329bb9d10680135f28ec7f3d705f373db..f67afba47a1a14499a46085d3d1b977b48bbb6fb 100644 (file)
@@ -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
index e2c70d60b90da964b58308277d838142f7b3a885..faaf4ea460956008024911150c03aecfefcca26b 100644 (file)
@@ -1357,7 +1357,7 @@ SparseMatrix<number>::precondition_SSOR (Vector<somenumber>              &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<number>::SSOR (Vector<somenumber>& 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(),
index 64f553d6722f0431a65421270cda3e56c0014dbd..da15f1528b0245397fab4a155cc0983e79159fc1 100644 (file)
@@ -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<number>::precondition_SSOR (Vector<somenumber>       &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<RowInfo>::const_reverse_iterator rri;

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.