We have a parallel solver :-) Hurra !

[[Imported from SVN: r3582]]

istl/ownercopy.hh

@@ -10,6 +10,8 @@
 #include <map>
 #include <set>
 
+#include "math.h"
+
 #include "mpi.h"
 
 #include "dune/common/tripel.hh"
@@ -81,7 +83,7 @@ namespace Dune {
       }
     };
 
-    void buildOwnerOverlapToAllInterface ()
+    void buildOwnerOverlapToAllInterface () const
     {
       if (OwnerOverlapToAllInterfaceBuilt)
         OwnerOverlapToAllInterface.free();
@@ -93,7 +95,7 @@ namespace Dune {
       OwnerOverlapToAllInterfaceBuilt = true;
     }
 
-    void buildOwnerToAllInterface ()
+    void buildOwnerToAllInterface () const
     {
       if (OwnerToAllInterfaceBuilt)
         OwnerToAllInterface.free();
@@ -109,7 +111,7 @@ namespace Dune {
   public:
 
     template<class T>
-    void copyOwnerToAll (const T& source, T& dest)
+    void copyOwnerToAll (const T& source, T& dest) const
     {
       if (!OwnerToAllInterfaceBuilt)
         buildOwnerToAllInterface ();
@@ -120,7 +122,7 @@ namespace Dune {
     }
 
     template<class T>
-    void addOwnerOverlapToAll (const T& source, T& dest)
+    void addOwnerOverlapToAll (const T& source, T& dest) const
     {
       if (!OwnerOverlapToAllInterfaceBuilt)
         buildOwnerOverlapToAllInterface ();
@@ -131,7 +133,7 @@ namespace Dune {
     }
 
     template<class T1, class T2>
-    void dot (const T1& x, const T1& y, T2& result)
+    void dot (const T1& x, const T1& y, T2& result) const
     {
       // set up mask vector
       if (mask.size()!=x.size())
@@ -139,41 +141,57 @@ namespace Dune {
         mask.resize(x.size());
         for (int i=0; i<mask.size(); i++) mask[i] = 1;
         for (typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i)
-          if (i->local().attribute()==copy)
+          if (i->local().attribute()!=owner)
             mask[i->local().local()] = 0;
       }
       result = 0;
-      for (int i=0; i<mask.size(); i++)
+      for (int i=0; i<x.size(); i++)
         result += x[i].operator*(y[i])*mask[i];
       int procs;
       MPI_Comm_size(comm,&procs);
       if (procs==1) return;
       double res;     // assumes that result is double \todo: template magick to treat complex<...>
       MPI_Allreduce(&result,&res,1,MPI_DOUBLE,MPI_SUM,comm);
+      result = res;
       return;
     }
 
     template<class T1>
-    double norm (const T1& x)
+    double norm (const T1& x) const
     {
+      int rank;
+      MPI_Comm_rank(comm,&rank);
+
       // set up mask vector
       if (mask.size()!=x.size())
       {
         mask.resize(x.size());
         for (int i=0; i<mask.size(); i++) mask[i] = 1;
         for (typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i)
-          if (i->local().attribute()==copy)
+          if (i->local().attribute()!=owner)
             mask[i->local().local()] = 0;
       }
       double result = 0;
-      for (int i=0; i<mask.size(); i++)
-        result += x[i].two_norm()*mask[i];
+      for (int i=0; i<x.size(); i++)
+      {
+        result += x[i].two_norm2()*mask[i];
+        //                if (mask[i]==1)
+        //                      std::cout << rank << ": " << "index=" << i << " value=" << x[i].two_norm2() << std::endl;
+      }
       int procs;
       MPI_Comm_size(comm,&procs);
-      if (procs==1) return result;
+      if (procs==1) return sqrt(result);
       double res;
       MPI_Allreduce(&result,&res,1,MPI_DOUBLE,MPI_SUM,comm);
-      return res;
+      return sqrt(res);
+    }
+
+    template<class T1>
+    void project (T1& x) const
+    {
+      for (typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i)
+        if (i->local().attribute()==copy)
+          x[i->local().local()] = 0;
     }
 
 
@@ -200,12 +218,12 @@ namespace Dune {
           pis.add(i->first,LI(i->second,overlap,true));
         if (i->third==copy)
           pis.add(i->first,LI(i->second,copy,true));
-        std::cout << rank << ": adding index " << i->first << " " << i->second << " " << i->third << std::endl;
+        //                std::cout << rank << ": adding index " << i->first << " " << i->second << " " << i->third << std::endl;
       }
       pis.endResize();
 
       // build remote indices WITHOUT communication
-      std::cout << rank << ": build remote indices" << std::endl;
+      //          std::cout << rank << ": build remote indices" << std::endl;
       ri.setIndexSets(pis,pis,comm);
       if (othersindices.size()>0)
       {
@@ -230,7 +248,7 @@ namespace Dune {
             DUNE_THROW(ISTLError,"OwnerOverlapCopyCommunication: global index not in index set");
 
           // insert entry
-          std::cout << rank << ": adding remote index " << i->first << " " << i->second << " " << i->third << std::endl;
+          //                      std::cout << rank << ": adding remote index " << i->first << " " << i->second << " " << i->third << std::endl;
           if (i->third==owner)
             modifier.insert(RX(owner,&(*pi)));
           if (i->third==overlap)
@@ -253,11 +271,11 @@ namespace Dune {
     MPI_Comm comm;
     PIS pis;
     RI ri;
-    IF OwnerToAllInterface;
-    bool OwnerToAllInterfaceBuilt;
-    IF OwnerOverlapToAllInterface;
-    bool OwnerOverlapToAllInterfaceBuilt;
-    std::vector<double> mask;
+    mutable IF OwnerToAllInterface;
+    mutable bool OwnerToAllInterfaceBuilt;
+    mutable IF OwnerOverlapToAllInterface;
+    mutable bool OwnerOverlapToAllInterfaceBuilt;
+    mutable std::vector<double> mask;
   };
 
 

istl/schwarz.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_SCHWARZ_HH
+#define DUNE_SCHWARZ_HH
+
+#include <iostream>               // for input/output to shell
+#include <fstream>                // for input/output to files
+#include <vector>                 // STL vector class
+#include <sstream>
+
+#include <math.h>                 // Yes, we do some math here
+#include <stdio.h>                // There is nothing better than sprintf
+#include <sys/times.h>           // for timing measurements
+
+#include "dune/common/timer.hh"
+
+#include "io.hh"
+#include "bvector.hh"
+#include "vbvector.hh"
+#include "bcrsmatrix.hh"
+#include "io.hh"
+#include "gsetc.hh"
+#include "ilu.hh"
+#include "operators.hh"
+#include "solvers.hh"
+#include "preconditioners.hh"
+#include "scalarproducts.hh"
+
+namespace Dune {
+
+  /**
+     @addtogroup ISTL
+     @{
+   */
+
+  /**
+   * @brief Categories for the solvers.
+   */
+
+  template<class M, class X, class Y, class C>
+  class OverlappingSchwarzOperator : public AssembledLinearOperator<M,X,Y>
+  {
+  public:
+    //! export types
+    typedef M matrix_type;
+    typedef X domain_type;
+    typedef Y range_type;
+    typedef typename X::field_type field_type;
+
+    //! define the category
+    enum {category=SolverCategory::overlapping};
+
+    //! constructor: just store a reference to a matrix
+    OverlappingSchwarzOperator (const M& A, const C& com)
+      : _A_(A), communication(com)
+    {}
+
+    //! apply operator to x:  \f$ y = A(x) \f$
+    virtual void apply (const X& x, Y& y) const
+    {
+      y = 0;
+      _A_.umv(x,y);     // result is consistent on interior+border
+      communication.project(y);
+    }
+
+    //! apply operator to x, scale and add:  \f$ y = y + \alpha A(x) \f$
+    virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
+    {
+      _A_.usmv(alpha,x,y);     // result is consistent on interior+border
+      communication.project(y);
+    }
+
+    //! get matrix via *
+    virtual const M& getmat () const
+    {
+      return _A_;
+    }
+
+  private:
+    const M& _A_;
+    const C& communication;
+  };
+
+
+
+  //! Scalar product assuming consistent vectors in interior+border
+  template<class X, class C>
+  class OverlappingSchwarzScalarProduct : public ScalarProduct<X>
+  {
+  public:
+    //! export types
+    typedef X domain_type;
+    typedef typename X::field_type field_type;
+
+    //! define the category
+    enum {category=SolverCategory::overlapping};
+
+    /*! \brief Constructor needs to know the grid
+     */
+    OverlappingSchwarzScalarProduct (const C& com)
+      : communication(com)
+    {}
+
+    /*! \brief Dot product of two vectors.
+       It is assumed that the vectors are consistent on the interior+border
+       partition.
+     */
+    virtual field_type dot (const X& x, const X& y)
+    {
+      field_type result;
+      communication.dot(x,y,result);
+      return result;
+    }
+
+    /*! \brief Norm of a right-hand side vector.
+       The vector must be consistent on the interior+border partition
+     */
+    virtual double norm (const X& x)
+    {
+      return communication.norm(x);
+    }
+
+  private:
+    const C& communication;
+  };
+
+
+
+
+
+  template<class M, class X, class Y, class C>
+  class ParSSOR : public Preconditioner<X,Y> {
+  public:
+    //! \brief The matrix type the preconditioner is for.
+    typedef M matrix_type;
+    //! \brief The domain type of the preconditioner.
+    typedef X domain_type;
+    //! \brief The range type of the preconditioner.
+    typedef Y range_type;
+    //! \brief The field type of the preconditioner.
+    typedef typename X::field_type field_type;
+
+    // define the category
+    enum {
+      //! \brief The category the precondtioner is part of.
+      category=SolverCategory::overlapping
+    };
+
+    /*! \brief Constructor.
+
+       constructor gets all parameters to operate the prec.
+       \param A The matrix to operate on.
+       \param n The number of iterations to perform.
+       \param w The relaxation factor.
+     */
+    ParSSOR (const M& A, int n, field_type w, const C& c)
+      : _A_(A), _n(n), _w(w), communication(c)
+    {   }
+
+    /*!
+       \brief Prepare the preconditioner.
+
+       \copydoc Preconditioner::pre(X&,Y&)
+     */
+    virtual void pre (X& x, Y& b) {}
+
+    /*!
+       \brief Apply the precondtioner
+
+       \copydoc Preconditioner::apply(X&,Y&)
+     */
+    virtual void apply (X& v, const Y& d)
+    {
+      for (int i=0; i<_n; i++) {
+        bsorf(_A_,v,d,_w);
+        bsorb(_A_,v,d,_w);
+      }
+      communication.copyOwnerToAll(v,v);
+    }
+
+    /*!
+       \brief Clean up.
+
+       \copydoc Preconditioner::post(X&)
+     */
+    virtual void post (X& x) {}
+
+  private:
+    //! \brief The matrix we operate on.
+    const M& _A_;
+    //! \brief The number of steps to do in apply
+    int _n;
+    //! \brief The relaxation factor to use
+    field_type _w;
+    //! \brief the communication object
+    const C& communication;
+  };
+
+  /** @} end documentation */
+
+} // end namespace
+
+#endif