avoid gcollect in load caculation

git-svn-id: https://dune.mathematik.uni-freiburg.de/svn/alugrid/trunk@1609 0d966ed9-3843-0410-af09-ebfb50bd7c74

src/parallel/gitter_pll_sti.cc

@@ -153,6 +153,7 @@ void GitterPll :: printsize ()
     cout << endl ;
   }
 
+  // could better use MPI_gather here 
   vector < vector < int > > in = mpAccess ().gcollect (n) ;
   assert (static_cast<int> (in.size ()) == np) ;
 
@@ -1099,16 +1100,34 @@ bool GitterPll :: checkPartitioning( LoadBalancer :: DataBase& db )
     // mean  - mittlere ElementLast
     // nload - Lastverh"altnis
 
+    // number of leaf elements 
     const double load = db.accVertexLoad () ;
 
+    // get:  min(load), max(load), sum(load)
+    const double (&minmaxsum)[3] = mpAccess ().minmaxsum( load ); 
+
+    // get mean value of leaf elements 
+    const double mean = minmaxsum[ 2 ] / double( np );
+
+    const double minload = minmaxsum[ 0 ];
+    const double maxload = minmaxsum[ 1 ];
+
+    /*
+    // old version using Allgather 
     vector < double > v (mpAccess ().gcollect (load)) ;
     const vector < double > :: iterator iEnd =  v.end () ;
 
     // sum up values and devide by number of cores 
     const double mean = accumulate (v.begin (), v.end (), 0.0) / double (np) ;
+    std::cout << mean << " mean value " << std::endl;
 
     for (vector < double > :: iterator i = v.begin () ; i != iEnd ; ++i)
       neu |= (*i > mean ? (*i > (_ldbOver * mean) ? true : false) : (*i < (_ldbUnder * mean) ? true : false)) ;
+    */
+    std::cout << mean << " mean value " << minload << "  minload  " << maxload << std::endl;
+
+    if( maxload > (_ldbOver * mean) || minload < (_ldbUnder * mean) ) 
+      neu = true ;
   }
 
 #ifndef NDEBUG

src/parallel/mpAccess.h

@@ -34,6 +34,7 @@ class MpAccessGlobal {
     virtual void gmax (int*,int,int*) const = 0 ;
     virtual void gmin (int*,int,int*) const = 0 ;
     virtual void gsum (int*,int,int*) const = 0 ;
+    virtual const double (& minmaxsum( double ) const)[3] = 0;
     virtual pair<double,double> gmax (pair<double,double>) const = 0 ;
     virtual pair<double,double> gmin (pair<double,double>) const = 0 ;
     virtual pair<double,double> gsum (pair<double,double>) const = 0 ;

src/parallel/mpAccess_MPI.cc

@@ -51,12 +51,15 @@ int MpAccessMPI :: getRank()
 
 MpAccessMPI :: MpAccessMPI (const MpAccessMPI & a)
 : _mpiCommPtr( a._mpiCommPtr->clone() ),
+  _minmaxsum( 0 ),
   _psize( getSize() ) , _myrank( getRank() )
 {
+  initMinMaxSum();
 }
 
 MpAccessMPI :: ~MpAccessMPI ()
 {
+  delete _minmaxsum;
   delete _mpiCommPtr;
   _mpiCommPtr = 0;
 }
@@ -299,6 +302,72 @@ void MpAccessMPI :: gsum (int* a,int size,int *x) const {
   assert (test == MPI_SUCCESS) ;
 }
 
+struct MinMaxSumOp : public MpAccessMPI :: MinMaxSumIF 
+{
+  MinMaxSumOp ( const MpAccessMPI& mpAccess ) 
+    : _mpAccess( mpAccess )
+  {
+    MPI_Type_contiguous (3, MPI_DOUBLE, &_minmaxsum_t);
+    MPI_Type_commit (&_minmaxsum_t);
+    MPI_Op_create((MPI_User_function *)&min_max_sum, TRUE, &_op);
+  }
+
+  ~MinMaxSumOp() 
+  {
+    MPI_Op_free (&_op);
+  }
+
+  const MpAccessMPI& _mpAccess ;
+
+  MPI_Op _op; 
+  MPI_Datatype _minmaxsum_t;
+
+  mutable double recvbuf[ 3 ];
+
+  static void 
+  min_max_sum( double* in, double* inout, 
+               int* len, MPI_Datatype* datatype ) 
+  {
+    const int size = *len ;
+    for( int i = 0; i < size; i += 3 )
+    {
+      inout[ i   ] = std::min( in[ i   ], inout[ i   ] ); 
+      inout[ i+1 ] = std::max( in[ i+1 ], inout[ i+1 ] ); 
+      inout[ i+2 ] = inout[ i+2 ] + in[ i+2 ] ;
+    }
+  } 
+
+  const double (& minmaxsum( double value ) const)[3] 
+  {
+    const MpAccessMPI :: CommIF* _mpiCommPtr = _mpAccess.mpiCommPtr();
+    // get mpi communicator (use define, see above)
+    MPI_Comm comm = _mpiComm ;
+
+    double sendbuf[ 3 ] = { value, value, value };
+    MPI_Allreduce( &sendbuf[ 0 ], &recvbuf[ 0 ], 3, _minmaxsum_t, _op, comm );
+    return recvbuf ;
+  }
+};
+
+void MpAccessMPI :: initMinMaxSum() 
+{
+  if( ! _minmaxsum ) 
+    _minmaxsum = new MinMaxSumOp( *this );
+}
+
+const double (& MpAccessMPI :: minmaxsum( double value ) const )[3]
+{
+  static double recvbuf[ 3 ];
+  recvbuf[ 0 ] = gmin( value );
+  recvbuf[ 1 ] = gmax( value );
+  recvbuf[ 2 ] = gsum( value );
+
+  // does not work yet due to some double free corruption 
+  //assert( _minmaxsum );
+  //return _minmaxsum->minmaxsum( value );
+  return recvbuf;
+}
+
 pair<double,double> MpAccessMPI :: gmax (pair<double,double> p) const {
   double x[2] ;
   double a[2]={p.first,p.second};

src/parallel/mpAccess_MPI.h

@@ -6,6 +6,15 @@
 class MpAccessMPI : public MpAccessLocal 
 {
 public:
+  class MinMaxSumIF 
+  {
+  protected:  
+      MinMaxSumIF () {}
+  public:  
+    virtual ~MinMaxSumIF() {}
+    virtual const double (& minmaxsum( double ) const )[3] = 0;
+  };
+
   typedef MpAccessGlobal :: CommIF CommIF;
 
   template <class MPICommunicator>
@@ -37,9 +46,10 @@ public:
   const CommIF* mpiCommPtr() const { return _mpiCommPtr; }
 
 protected:  
-
   // class holding the MPI communicator 
   const CommIF* _mpiCommPtr;
+  // pointer to minmaxsum communication 
+  const MinMaxSumIF* _minmaxsum;
   // number of processors
   const int _psize; 
   // my processor number  
@@ -48,14 +58,18 @@ protected:
   int mpi_allgather (int *, int , int *, int) const ;
   int mpi_allgather (char *, int, char *, int) const ;
   int mpi_allgather (double *, int, double *, int ) const ;
+
+  void initMinMaxSum() ;
 public :
   // constructor taking MPI_Comm 
   // to avoid MPI types here this is a template constructor 
   template <class MPICommunicator>  
   inline MpAccessMPI (MPICommunicator mpicomm ) 
     : _mpiCommPtr( new Comm<MPICommunicator> ( mpicomm ) ), 
+      _minmaxsum( 0 ),
       _psize( getSize() ), _myrank( getRank() )
   {
+    initMinMaxSum();
   }
 
   // copy constructor 
@@ -85,6 +99,7 @@ public:
   void gmax (int*,int,int*) const ;
   void gmin (int*,int,int*) const ;
   void gsum (int*,int,int*) const ;
+  const double (& minmaxsum( double ) const)[3] ;
   pair<double,double> gmax (pair<double,double>) const ;
   pair<double,double> gmin (pair<double,double>) const ;
   pair<double,double> gsum (pair<double,double>) const ;