parallel mtl distribution

AMDiS/src/MTL4Solver.h

@@ -23,6 +23,12 @@
 #include "OEMSolver.h"
 #include "MatrixStreams.h"
 #include <iostream>
+#include <boost/mpl/bool.hpp>
+#include <boost/numeric/mtl/utility/is_distributed.hpp>
+
+#if defined(HAVE_PARALLEL_DOMAIN_AMDIS) && defined(HAVE_PARALLEL_MTL4)
+# include <boost/numeric/mtl/par/distribution.hpp>
+#endif
 
 namespace AMDiS {
 
@@ -32,23 +38,52 @@ namespace AMDiS {
   protected:    
     Worker worker;
 
+    template< typename Mapper, bool IsDist>
+    void init(Mapper& mapper, boost::mpl::bool_<IsDist>)
+    {
+        matrix.change_dim(mapper.getNumRows(), mapper.getNumCols());
+    }
+
+#if defined(HAVE_PARALLEL_DOMAIN_AMDIS) && defined(HAVE_PARALLEL_MTL4)
+    void init(ParallelMapper& mapper, boost::mpl::true_)
+    {
+	mtl::par::block_distribution dist(mapper.getNumRows() /* , communicator */ );
+	dist.setup_from_local_size(mapper.meshDistributor().getNumberRankDofs()*mapper.getNumComponents());
+        matrix.init_distribution(dist, dist, mapper.getNumRows(), mapper.getNumRows());
+    }
+#endif
+
+    template< bool IsDist>
+    void init(MTLVector& v, boost::mpl::bool_<IsDist>)
+    {
+      v.change_dim(num_rows(matrix));
+    }
+
+    void init(MTLVector& v, boost::mpl::true_)
+    {
+      v.init_distribution(row_distribution(matrix), num_rows(matrix));
+    }
+
     template< typename Matrix, typename Vector, typename Mapper >
     int solve(const Matrix& A, Vector& x, Vector& b, Mapper& mapper) 
     {
       if(num_rows(matrix) == 0 || !getMultipleRhs() ) {
-        matrix.change_dim(mapper.nRow(), mapper.nCol());
+	init(mapper, mtl::traits::is_distributed<MTLMatrix>());
+
         set_to_zero(matrix);
         MatMap< const Matrix, Mapper > matMap(A,mapper);
         matrix << matMap;
 	worker.init(matrix);
       }
 
-      MTLVector mtl_x(mapper.nRow());
+      MTLVector mtl_x;
+      init(mtl_x,mtl::traits::is_distributed<MTLVector>());
       set_to_zero(mtl_x);
       VecMap< Vector, Mapper > xVecMap(x, mapper);
       mtl_x << xVecMap;
 
-      MTLVector mtl_b(mapper.nRow());
+      MTLVector mtl_b;
+      init(mtl_b,mtl::traits::is_distributed<MTLVector>());
       set_to_zero(mtl_b);
       VecMap< Vector, Mapper> bVecMap(b, mapper);
       mtl_b << bVecMap;

AMDiS/src/Mapper.h

@@ -51,9 +51,11 @@ namespace AMDiS {
 
     inline int col(int c) const { return c + colOffset; }
 
-    inline int nCol() const { return ncol; }
+    inline int getNumCols() const { return ncol; }
 
-    inline int nRow() const { return nrow; }
+    inline int getNumRows() const { return nrow; }
+
+    inline int getNumComponents() { return nComp; }
 
   protected:
 

AMDiS/src/SolverMatrix.cc

@@ -6,7 +6,7 @@ namespace AMDiS {
   void SolverMatrix< Matrix<DOFMatrix*> >::buildMatrix() const
   {
     VectorialMapper mapper(*this);
-    matrix.change_dim(mapper.nRow(), mapper.nCol());
+    matrix.change_dim(mapper.getNumRows(), mapper.getNumCols());
     set_to_zero(matrix);
     MatMap< const SolverMatrix< Matrix<DOFMatrix*> >, VectorialMapper > matMap(*this,mapper);
     matrix << matMap;

AMDiS/src/parallel/ParallelMapper.h

@@ -22,21 +22,24 @@ namespace AMDiS {
 		}
 		size_type row(size_type r_) 
 		{
-			//size_type ret = md.mapLocalToGlobal(r_)*nComponents+r;
-			size_type ret = md.mapLocalToGlobal(r_)+r*md.getNumberOverallDofs() ;
+			size_type ret = md.mapLocalToGlobal(r_)*nComponents+r;
+			//size_type ret = md.mapLocalToGlobal(r_)+r*md.getNumberOverallDofs() ;
 			return ret;
 		}
 		size_type col(size_type c_) 
 		{
-			//size_type ret = md.mapLocalToGlobal(c_)*nComponents+c;
-			size_type ret = md.mapLocalToGlobal(c_)+c*md.getNumberOverallDofs() ;
+			size_type ret = md.mapLocalToGlobal(c_)*nComponents+c;
+			//size_type ret = md.mapLocalToGlobal(c_)+c*md.getNumberOverallDofs() ;
 			return ret ;
 		}
 
+		inline MeshDistributor& meshDistributor() { return md; }
+		inline int getNumComponents() { return nComponents; }
+
 		inline void setRow(size_type r_) { r = r_ ;}
 		inline void setCol(size_type c_) { c = c_; }
-		inline size_type nRow() { return md.getNumberOverallDofs() * nComponents; }
-		inline size_type nCol() { return md.getNumberOverallDofs() * nComponents; }
+		inline size_type getNumRows() { return md.getNumberOverallDofs() * nComponents; }
+		inline size_type getNumCols() { return md.getNumberOverallDofs() * nComponents; }
 	};
 
 	typedef ParallelMapper_base< unsigned int > ParallelMapper;