linear_algebra_distribution.cc

Go to the documentation of this file.

//LIC// ====================================================================
//LIC// This file forms part of oomph-lib, the object-oriented, 
//LIC// multi-physics finite-element library, available 
//LIC// at http://www.oomph-lib.org.
//LIC// 
//LIC//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
//LIC// 
//LIC// This library is free software; you can redistribute it and/or
//LIC// modify it under the terms of the GNU Lesser General Public
//LIC// License as published by the Free Software Foundation; either
//LIC// version 2.1 of the License, or (at your option) any later version.
//LIC// 
//LIC// This library is distributed in the hope that it will be useful,
//LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
//LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//LIC// Lesser General Public License for more details.
//LIC// 
//LIC// You should have received a copy of the GNU Lesser General Public
//LIC// License along with this library; if not, write to the Free Software
//LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
//LIC// 02110-1301  USA.
//LIC// 
//LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
//LIC// 
//LIC//====================================================================
#include "linear_algebra_distribution.h"

namespace oomph
{


 //============================================================================
 /// Sets the distribution. Takes first_row, local_nrow and 
 /// global_nrow as arguments. If global_nrow is not provided or equal to
 /// 0 then it is computed automatically
 //============================================================================
 void LinearAlgebraDistribution::build(const OomphCommunicator* 
                                       const comm_pt,
                                       const unsigned& first_row, 
                                       const unsigned& local_nrow,
                                       const unsigned& global_nrow)
 {
  // copy the communicator
  delete Comm_pt;
  Comm_pt = new OomphCommunicator(*comm_pt);

  // get the rank and the number of processors
  int my_rank = Comm_pt->my_rank();
  int nproc = Comm_pt->nproc();

  // resize the storage
  First_row.clear();
  First_row.resize(nproc);
  Nrow_local.clear();
  Nrow_local.resize(nproc);
  
  // set first row and local nrow for this processor
  First_row[my_rank] = first_row;
  Nrow_local[my_rank] = local_nrow;
  
#ifdef OOMPH_HAS_MPI
  // gather the First_row vector
  unsigned my_nrow_local = Nrow_local[my_rank];
  MPI_Allgather(&my_nrow_local,1,MPI_UNSIGNED,
                &Nrow_local[0],1,MPI_UNSIGNED,
                Comm_pt->mpi_comm());

  // gather the Nrow_local vector
  unsigned my_first_row = First_row[my_rank];
  MPI_Allgather(&my_first_row,1,MPI_UNSIGNED,
                &First_row[0],1,MPI_UNSIGNED,
                Comm_pt->mpi_comm());    
#endif

  // if global nrow is not provided then compute by summing local_nrow over
  // all processors
  if (global_nrow == 0)
   {
    if (nproc == 1)
     {
      Nrow = local_nrow;
     }
    else
     {
      Nrow = 0;
      for (int p = 0; p < nproc; p++)
       {
        Nrow += Nrow_local[p];
       }
     }
   }
  else
   {
    Nrow = global_nrow;
   }

  // the distribution is true, unless we only have 1 processor
  if(nproc!=1) {Distributed = true;}
  else {Distributed = false;}

#ifdef OOMPH_HAS_MPI  
#ifdef PARANOID
  // paranoid check that the distribution works


  // check that none of the processors partition overlap
  for (int p = 0; p < nproc; p++)
   {
    if (Nrow_local[p] > 0)
     {
      for (int pp = p+1; pp < nproc; pp++)
       {
        if (Nrow_local[pp] > 0)
         {
          if ((First_row[p] >= First_row[pp] &&
               First_row[p] < First_row[pp] + Nrow_local[pp]) ||
              (First_row[p] + Nrow_local[p] -1 >= First_row[pp] &&
               First_row[p] + Nrow_local[p] -1 <
               First_row[pp] + Nrow_local[pp]))
           {
            std::ostringstream error_message;
            error_message << "The distributed rows on processor " << p
                          << " and processor " << pp << " overlap.\n"
                          << "Processor " << p << " : first_row = "
                          << First_row[p] << ", nrow = " 
                          << Nrow_local[p] << ".\n"
                          << "Processor " << pp << " : first_row = "
                          << First_row[pp] << ", nrow = " 
                          << Nrow_local[pp] << ".\n";
            throw OomphLibWarning(error_message.str(),
                                  "LinearAlgebraDistribution::distribute(...)",
                                  OOMPH_EXCEPTION_LOCATION);
           }
         }
       }
     }
   }
  
  // check that no processor has a row with a global row index greater than
  // the number of global rows
  for (int p = 0; p < nproc; p++)
   {
    if (First_row[p] + Nrow_local[p] > Nrow)
     {
      std::ostringstream error_message;
      error_message << "Processor " << p << " contains rows "
                    << First_row[p] << " to " 
                    << First_row[p]+Nrow_local[p]-1
                    << " but there are only " << Nrow 
                    << " to be distributed." << std::endl;
      throw OomphLibWarning(error_message.str(),
                            "LinearAlgebraDistribution::distribute(...)",
                            OOMPH_EXCEPTION_LOCATION);
     }
   }
#endif
#endif
 }
 
 //============================================================================
 /// \short Uniformly distribute global_nrow over all processors where 
 /// processors 0 holds approximately the first 
 /// global_nrow/n_proc, processor 1 holds the next 
 /// global_nrow/n_proc and so on...
 //============================================================================
 void LinearAlgebraDistribution::build(const OomphCommunicator* 
                                       const comm_pt,
                                       const unsigned& global_nrow,
                                       const bool& distribute)
 {
  // copy the communicator
  delete Comm_pt;
  Comm_pt = new OomphCommunicator(*comm_pt);

  // delete existing storage
  First_row.clear();
  Nrow_local.clear();

  // set global nrow
  Nrow = global_nrow;
  
  // store the distributed flag
  Distributed = distribute;

#ifdef OOMPH_HAS_MPI

  // if distributed object then compute uniform distribution
  if (distribute == true)
   {

    // get the number of processors
    int nproc = Comm_pt->nproc();

    // resize the storage
    First_row.resize(nproc);
    Nrow_local.resize(nproc);

    // compute first row
    for (int p=0;p<nproc;p++)
     {                                          
      First_row[p] = unsigned((double(p)/double(nproc))*double(global_nrow)); 
     }             
    
    // compute local nrow
    for (int p=0; p<nproc-1; p++) 
     {                                                  
      Nrow_local[p] = First_row[p+1] - First_row[p];
     }  
    Nrow_local[nproc-1] = global_nrow -  
     First_row[nproc-1];         
   }
#endif
 }

 //============================================================================
 /// helper method for the =assignment operator and copy constructor
 //============================================================================
 void LinearAlgebraDistribution::build(const LinearAlgebraDistribution& 
                                       new_dist)
  {

   // delete the existing storage
   First_row.clear();
   Nrow_local.clear();

   // if new_dist is not setup
   if (new_dist.communicator_pt() == 0)
    {
     delete Comm_pt;
     Comm_pt = 0;
     Distributed = true;
     Nrow = 0;
     }
   else
    {
     // copy the communicator
     delete Comm_pt;
     Comm_pt = new OomphCommunicator(*new_dist.communicator_pt());
         
     // the new distribution is distributed
     if (new_dist.distributed())
      {
       First_row = new_dist.first_row_vector();
       Nrow_local = new_dist.nrow_local_vector();
       
       Distributed = true;
      }
     // else if the new ditribution is not distributed
     else
      {
       Distributed = false;
      }
     Nrow = new_dist.nrow();
    }
  }


 //============================================================================
 /// operator==
 //============================================================================
 bool LinearAlgebraDistribution::operator==
 (const LinearAlgebraDistribution& other_dist) const
  {   
#ifdef OOMPH_HAS_MPI
   // compare the communicators
   if (!((*Comm_pt) == (*other_dist.communicator_pt())))
    {
     return false;
    }

   // compare Distributed
   if (Distributed != other_dist.distributed())
    {
     return false;
    }

   // if not distributed compare nrow
   if (!Distributed)
    {
     if (other_dist.nrow() == Nrow)
      {
       return true;
      }
     return false;
    }
   
   // compare
   bool flag = true;
   int nproc = Comm_pt->nproc();
   for (int i = 0; i < nproc && flag == true; i++)
    {
     if (other_dist.first_row(i) != First_row[i] ||
         other_dist.nrow_local(i) != Nrow_local[i])
      {
       flag = false;
      }
    }
   return flag;
#else
   if (other_dist.nrow() == Nrow)
    {
     return true;
    }
   return false;
#endif
  }

 //=============================================================================
 /// output operator
 //=============================================================================
 std::ostream& operator<<(std::ostream& stream, 
                          LinearAlgebraDistribution& dist)
 {
  stream << "nrow()=" << dist.nrow() 
         << ", first_row()=" << dist.first_row()
         << ", nrow_local()=" << dist.nrow_local()
         << ", distributed()=" << dist.distributed()
         << std::endl;
  return stream;
 }

//=============================================================================
/// Namespace for helper functions for LinearAlgebraDistributions
//=============================================================================
 namespace LinearAlgebraDistributionHelpers
 {

  //===========================================================================
  /// \short Takes a vector of LinearAlgebraDistribution objects and 
  /// concatenates them such that the nrow_local of the out_distribution
  /// is the sum of the nrow_local of all the in_distributions and the number
  /// of global rows of the out_distribution is the sum of the number of global
  /// rows of all the in_distributions. 
  /// This results in a permutation of the rows in the out_distribution. 
  /// Think of this in terms of DoubleVectors, if we have DoubleVectors with 
  /// distributions A and B, distributed across two processors (p0 and p1),
  /// A: [a0] (on p0)    B: [b0] (on p0) 
  ///    [a1] (on p1)       [b1] (on P1), 
  /// 
  /// then the out_distribution is
  /// [a0  (on p0) 
  ///  b0] (on p0) 
  /// [a1  (on p1) 
  ///  b1] (on p1), 
  /// 
  /// as opposed to 
  /// [a0  (on p0) 
  ///  a1] (on p0) 
  /// [b0  (on p1) 
  ///  b1] (on p1). 
  /// 
  /// Note (1): The out_distribution may not be uniformly distributed even
  /// if the in_distributions are uniform distributions. 
  /// Try this out with two distributions of global rows 3 and 5, uniformly
  /// distributed across two processors. Compare this against a distribution
  /// of global row 8 distributed across two processors.
  /// 
  /// Note (2): There is no equivalent function which takes a Vector of
  /// LinearAlgebraDistribution objects (as opposed to pointers), there should
  /// not be one since we do not want to invoke the assignment operator when
  /// creating the Vector of LinearAlgebraDistribution objects.
  //===========================================================================
  void concatenate(const Vector<LinearAlgebraDistribution*>&in_distribution_pt,
                   LinearAlgebraDistribution &out_distribution)
  {
   // How many distributions are in in_distribution?
   unsigned ndistributions=in_distribution_pt.size();

#ifdef PARANOID
   
   // Are any of the in_distribution pointers null? 
   // We do not want null pointers.
   for (unsigned dist_i = 0; dist_i < ndistributions; dist_i++) 
    {
     if(in_distribution_pt[dist_i]==0)
      {
       std::ostringstream error_message;
       error_message
        << "The pointer in_distribution_pt[" << dist_i << "] is null.\n";
       throw OomphLibError(error_message.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
      }
    }

   ////////////////////////////////
    
   // Check that all in distributions are built.
   for (unsigned dist_i = 0; dist_i < ndistributions; dist_i++)
    {
     if(!in_distribution_pt[dist_i]->built())
      {
       std::ostringstream error_message;
       error_message
        << "The in_distribution_pt[" << dist_i << "] is not built.\n";
       throw OomphLibError(error_message.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
      }
    }

   ////////////////////////////////
   
   // Check that all communicators to concatenate are the same 
   // by comparing all communicators against the first one.
   const OomphCommunicator first_comm 
    = *(in_distribution_pt[0]->communicator_pt());
   
   for (unsigned dist_i = 0; dist_i < ndistributions; dist_i++) 
    {
     // Get the Communicator for the current vector.
     const OomphCommunicator another_comm
      = *(in_distribution_pt[dist_i]->communicator_pt());
     
     if(!(first_comm == another_comm))
      {
       std::ostringstream error_message;
       error_message << "The communicator in position "<<dist_i<<" is \n"
                     << "not the same as the first one.\n";
       throw OomphLibError(error_message.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
      }
    }
   
   ////////////////////////////////

   // Ensure that all distributions are either distributed or not. 
   // This is because we use the distributed() function from the first 
   // distribution to indicate if the result distribution should be distributed
   // or not.
   bool first_distributed = in_distribution_pt[0]->distributed();
   for (unsigned dist_i = 0; dist_i < ndistributions; dist_i++) 
    {
     // Is the current distribution distributed?
     bool another_distributed = in_distribution_pt[dist_i]->distributed();
     if(first_distributed != another_distributed)
      {
       std::ostringstream error_message;
       error_message 
        <<"The distribution in position "<<dist_i<<" has a different\n"
        <<"different distributed boolean than the distribution.\n";
       throw OomphLibError(error_message.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
      }
    }      

   ////////////////////////////////
   
   // Check that the out distribution is not built.
   if(out_distribution.built())
    {
     std::ostringstream error_message;
     error_message
      << "The out distribution is built.\n"
      << "Please clear it.\n";
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
     
#endif
 
   // Get the communicator pointer
   const OomphCommunicator* const comm_pt 
    = in_distribution_pt[0]->communicator_pt();
   
   // Number of processors
   unsigned nproc = comm_pt->nproc();
   
   // First determine the out_nrow_local and the out_nrow (the global nrow)
   unsigned out_nrow_local = 0;
   unsigned out_nrow = 0;
   for (unsigned in_dist_i = 0; in_dist_i < ndistributions; in_dist_i++) 
    {
     out_nrow_local += in_distribution_pt[in_dist_i]->nrow_local();
     out_nrow += in_distribution_pt[in_dist_i]->nrow();
    }
   
   // Now calculate the first row for this processor. We need to know the 
   // out_nrow_local for all the other processors, MPI_Allgather must be used.
   unsigned out_first_row = 0;
   
   // Distributed case: We need to gather all the out_nrow_local from all
   // processors, the out_first_row for this processors is the partial sum up 
   // of the out_nrow_local to my_rank.
   bool distributed = in_distribution_pt[0]->distributed();
   if(distributed)
    {
#ifdef OOMPH_HAS_MPI
     // My rank
     unsigned my_rank = comm_pt->my_rank();

     unsigned *out_nrow_local_all = new unsigned[nproc];
     MPI_Allgather(&out_nrow_local,1,MPI_UNSIGNED,
                   out_nrow_local_all,1,MPI_UNSIGNED,
                   comm_pt->mpi_comm());
     for (unsigned proc_i = 0; proc_i < my_rank; proc_i++) 
      {
       out_first_row += out_nrow_local_all[proc_i];
      }
     delete [] out_nrow_local_all;
#endif
    }
   // Non-distributed case
   else
    {
     out_first_row = 0;
    }
   
   // Build the distribution.
   if(nproc == 1 || !distributed)
    {
     // Some ambiguity here -- on a single processor it doesn't really
     // matter if we think of ourselves as distributed or not but this
     // follows the pattern used elsewhere.
     out_distribution.build(comm_pt,out_nrow,false);
    }
   else
    {
     out_distribution.build(comm_pt,out_first_row,out_nrow_local,out_nrow);
    }
  } // function concatenate

 } // end of namespace LinearAlgebraDistributionHelpers

}//end of oomph namespace