Merge branch 'cherry-pick-650cbeac' into 'master'

all of FETI_Jump in a single commit See merge request jschoeberl/ngsolve!322
NGSolve · Mar 19, 2018 · eb4a85a · eb4a85a
2 parents 992c8cc + 79ff365
commit eb4a85a
Show file tree

Hide file tree

Showing 2 changed files with 97 additions and 22 deletions.
diff --git a/parallel/parallel_matrices.cpp b/parallel/parallel_matrices.cpp
@@ -484,15 +484,82 @@ namespace ngla
 
 
 
+  FETI_Jump_Matrix :: FETI_Jump_Matrix (shared_ptr<ParallelDofs> apardofs)
+      : BaseMatrix(apardofs)
+  {
+
+    size_t njs = 0;
+    for(auto p:paralleldofs->GetDistantProcs())
+      njs += paralleldofs->GetExchangeDofs(p).Size();
+
+    Array<size_t> ones(njs);
+    ones = 1;
+    Table<int>* dps = new Table<int>(ones);
+    njs = 0;
+    for(auto p:paralleldofs->GetDistantProcs()) {
+      for(auto d:paralleldofs->GetExchangeDofs(p)) {
+	(*dps)[njs++][0] = p;
+      }
+    }    
 
+    this->jump_paralleldofs = make_shared<ParallelDofs>(paralleldofs->GetCommunicator(), dps);
+
+    return;
+  }
 
 
   void FETI_Jump_Matrix :: MultAdd (double s, const BaseVector & x, BaseVector & y) const
-  { ; }
+  {
+    y.Distribute();
+    size_t count = 0;
+    for(auto p:paralleldofs->GetDistantProcs()) {
+      auto exdofs = paralleldofs->GetExchangeDofs(p);
+      if(p<MyMPI_GetId(paralleldofs->GetCommunicator())) {
+	for(auto k:Range(exdofs.Size())) {
+	  y.FVDouble()[count++] -= s*x.FVDouble()[exdofs[k]];
+	}
+      }
+      else {
+	for(auto k:Range(exdofs.Size())) {
+	  y.FVDouble()[count++] += s*x.FVDouble()[exdofs[k]];
+	}
+      }
+    }
+    return;
+  }
 
   void FETI_Jump_Matrix :: MultTransAdd (double s, const BaseVector & x, BaseVector & y) const
-  { ; }
+  {
+    x.Cumulate();
+    size_t count = 0;
+    for(auto p:paralleldofs->GetDistantProcs()) {
+      auto exdofs = paralleldofs->GetExchangeDofs(p);
+      if(p<MyMPI_GetId(paralleldofs->GetCommunicator())) {
+	for(auto k:Range(exdofs.Size())) {
+	  y.FVDouble()[exdofs[k]] -= s*x.FVDouble()[count++];
+	}
+      }
+      else {
+	for(auto k:Range(exdofs.Size())) {
+	  y.FVDouble()[exdofs[k]] += s*x.FVDouble()[count++];
+	}
+      }
+    }
+    return;
+  }
 
+  AutoVector FETI_Jump_Matrix :: CreateRowVector () const
+  {
+    return make_shared<VVector<double>> (paralleldofs->GetNDofLocal());
+  }
+
+  AutoVector FETI_Jump_Matrix :: CreateColVector () const
+  {
+    return make_shared<ParallelVVector<double>> (jump_paralleldofs->GetNDofLocal(),
+						 jump_paralleldofs);
+  }
+
+
 }
 
 #endif
diff --git a/parallel/parallel_matrices.hpp b/parallel/parallel_matrices.hpp
@@ -43,45 +43,53 @@ namespace ngla
     // : mat(*amat), pardofs(*apardofs) 
     // {const_cast<BaseMatrix&>(mat).SetParallelDofs (apardofs);}
 
-    virtual ~ParallelMatrix () override;
-    virtual bool IsComplex() const override { return mat->IsComplex(); } 
-    virtual void MultAdd (double s, const BaseVector & x, BaseVector & y) const override;
-    virtual void MultTransAdd (double s, const BaseVector & x, BaseVector & y) const override;
+    virtual ~ParallelMatrix ();
+    virtual bool IsComplex() const { return mat->IsComplex(); } 
+    virtual void MultAdd (double s, const BaseVector & x, BaseVector & y) const;
+    virtual void MultTransAdd (double s, const BaseVector & x, BaseVector & y) const;
 
-    virtual BaseVector & AsVector() override { return mat->AsVector(); }
-    virtual const BaseVector & AsVector() const override { return mat->AsVector(); }
+    virtual BaseVector & AsVector() { return mat->AsVector(); }
+    virtual const BaseVector & AsVector() const { return mat->AsVector(); }
 
     shared_ptr<BaseMatrix> GetMatrix() const { return mat; }
-    virtual shared_ptr<BaseMatrix> CreateMatrix () const override;
-    virtual AutoVector CreateVector () const override;
+    virtual shared_ptr<BaseMatrix> CreateMatrix () const;
+    virtual AutoVector CreateVector () const;
 
-    virtual ostream & Print (ostream & ost) const override;
+    virtual ostream & Print (ostream & ost) const;
 
-    virtual int VHeight() const override;
-    virtual int VWidth() const override;
+    virtual int VHeight() const;
+    virtual int VWidth() const;
 
     // virtual const ParallelDofs * GetParallelDofs () const {return &pardofs;}
 
 
-    virtual shared_ptr<BaseMatrix> InverseMatrix (shared_ptr<BitArray> subset = 0) const override;
+    virtual shared_ptr<BaseMatrix> InverseMatrix (shared_ptr<BitArray> subset = 0) const;
     template <typename TM>
     shared_ptr<BaseMatrix> InverseMatrixTM (shared_ptr<BitArray> subset = 0) const;
 
     virtual shared_ptr<BaseMatrix> InverseMatrix (shared_ptr<const Array<int>> clusters) const override;
-    virtual INVERSETYPE SetInverseType ( INVERSETYPE ainversetype ) const override;
-    virtual INVERSETYPE SetInverseType ( string ainversetype ) const override;
-    virtual INVERSETYPE  GetInverseType () const override;
+    virtual INVERSETYPE SetInverseType ( INVERSETYPE ainversetype ) const;
+    virtual INVERSETYPE SetInverseType ( string ainversetype ) const;
+    virtual INVERSETYPE  GetInverseType () const;
   };
 
+
   class FETI_Jump_Matrix : public BaseMatrix
   {
   public:
-    FETI_Jump_Matrix (shared_ptr<ParallelDofs> pardofs)
-      : BaseMatrix(pardofs) { ; }
+    FETI_Jump_Matrix (shared_ptr<ParallelDofs> pardofs);
 
-    virtual bool IsComplex() const { return false; }
-    virtual void MultAdd (double s, const BaseVector & x, BaseVector & y) const;
-    virtual void MultTransAdd (double s, const BaseVector & x, BaseVector & y) const;
+    virtual bool IsComplex() const override { return false; }
+    virtual void MultAdd (double s, const BaseVector & x, BaseVector & y) const override;
+    virtual void MultTransAdd (double s, const BaseVector & x, BaseVector & y) const override;
+
+    virtual AutoVector CreateRowVector () const override;
+    virtual AutoVector CreateColVector () const override;
+
+  protected:
+
+    shared_ptr<ParallelDofs> jump_paralleldofs;
+
   };
 
 #endif