distributed_vector.h

一个用来实现偏微分方程中网格的计算库

   */
  void localize_to_one (std::vector<T>& v_local,
			const unsigned int proc_id=0) const;
    
private:

  /**
   * Actual vector datatype
   * to hold vector entries
   */
  std::vector<T> _values;

  /**
   * The global vector size
   */
  unsigned int _global_size;

  /**
   * The local vector size
   */
  unsigned int _local_size;

  /**
   * The first component stored locally
   */
  unsigned int _first_local_index;

  /**
   * The last component (+1) stored locally
   */
  unsigned int _last_local_index;
};


//--------------------------------------------------------------------------
// DistributedVector inline methods
template <typename T>
inline
DistributedVector<T>::DistributedVector () :
  _global_size      (0),
  _local_size       (0),
  _first_local_index(0),
  _last_local_index (0)
{
}



template <typename T>
inline
DistributedVector<T>::DistributedVector (const unsigned int n)
{
  this->init(n, n, false);
}



template <typename T>
inline
DistributedVector<T>::DistributedVector (const unsigned int n,
					 const unsigned int n_local)
{
  this->init(n, n_local, false);
}



template <typename T>
inline
DistributedVector<T>::~DistributedVector ()
{
  this->clear ();
}



template <typename T>
inline
void DistributedVector<T>::init (const unsigned int n,
				 const unsigned int n_local,
				 const bool fast)
{
  // This function must be run on all processors at once
  parallel_only();

  libmesh_assert (n_local <= n);

  // Clear the data structures if already initialized
  if (this->initialized())
    this->clear();
    
  // Initialize data structures
  _values.resize(n_local);
  _local_size  = n_local;
  _global_size = n;

  _first_local_index = 0;
  
#ifdef HAVE_MPI

  int n_proc=0, proc_id=0;
  
  MPI_Comm_rank (libMesh::COMM_WORLD, &proc_id);
  MPI_Comm_size (libMesh::COMM_WORLD, &n_proc);
  
  std::vector<int> local_sizes     (n_proc, 0);
  
  local_sizes[proc_id] = n_local;

  Parallel::sum(local_sizes);

  // _first_local_index is the sum of _local_size
  // for all processor ids less than ours
  for (int p=0; p<proc_id; p++)
    _first_local_index += local_sizes[p];


#  ifdef DEBUG
  // Make sure all the local sizes sum up to the global
  // size, otherwise there is big trouble!
  int sum=0;

  for (int p=0; p<n_proc; p++)
    sum += local_sizes[p];

  libmesh_assert (sum == static_cast<int>(n));
  
#  endif
  
#else
  
  // No other options without MPI!
  if (n != n_local)
    {
      std::cerr << "ERROR:  MPI is required for n != n_local!"
		<< std::endl;
      libmesh_error();
    }
  
#endif

  _last_local_index = _first_local_index + n_local;
  
  // Set the initialized flag
  this->_is_initialized = true;

  // Zero the components unless directed otherwise
  if (!fast)
    this->zero();
}



template <typename T>
inline
void DistributedVector<T>::init (const unsigned int n,
				 const bool fast)
{
  this->init(n,n,fast);
}



template <typename T>
inline
void DistributedVector<T>::close ()
{
  libmesh_assert (this->initialized());
  
  this->_is_closed = true;
}



template <typename T>
inline
void DistributedVector<T>::clear ()
{
  _values.clear();
  
  _global_size =
    _local_size =
    _first_local_index =
    _last_local_index = 0;
  
  
  this->_is_closed = this->_is_initialized = false;
}



template <typename T>
inline
void DistributedVector<T>::zero ()
{
  libmesh_assert (this->initialized());
  libmesh_assert (_values.size() == _local_size);
  libmesh_assert ((_last_local_index - _first_local_index) == _local_size);
  
  std::fill (_values.begin(),
	     _values.end(),
	     0.);
}



template <typename T>
inline
AutoPtr<NumericVector<T> > DistributedVector<T>::clone () const
{
  AutoPtr<NumericVector<T> > cloned_vector (new DistributedVector<T>);

  *cloned_vector = *this;

  return cloned_vector;
}



template <typename T>
inline
unsigned int DistributedVector<T>::size () const
{
  libmesh_assert (this->initialized());
  libmesh_assert (_values.size() == _local_size);
  libmesh_assert ((_last_local_index - _first_local_index) == _local_size);

  return _global_size;
}



template <typename T>
inline
unsigned int DistributedVector<T>::local_size () const
{
  libmesh_assert (this->initialized());
  libmesh_assert (_values.size() == _local_size);
  libmesh_assert ((_last_local_index - _first_local_index) == _local_size);

  return _local_size;
}



template <typename T>
inline
unsigned int DistributedVector<T>::first_local_index () const
{
  libmesh_assert (this->initialized());
  libmesh_assert (_values.size() == _local_size);
  libmesh_assert ((_last_local_index - _first_local_index) == _local_size);

  return _first_local_index;
}



template <typename T>
inline
unsigned int DistributedVector<T>::last_local_index () const
{
  libmesh_assert (this->initialized());
  libmesh_assert (_values.size() == _local_size);
  libmesh_assert ((_last_local_index - _first_local_index) == _local_size);

  return _last_local_index;
}



template <typename T>
inline
T DistributedVector<T>::operator() (const unsigned int i) const
{
  libmesh_assert (this->initialized());
  libmesh_assert (_values.size() == _local_size);
  libmesh_assert ((_last_local_index - _first_local_index) == _local_size);
  libmesh_assert ( ((i >= first_local_index()) &&
	    (i <  last_local_index())) );

  return _values[i - _first_local_index];
}



template <typename T>
inline
void DistributedVector<T>::set (const unsigned int i, const T value)
{
  libmesh_assert (this->initialized());
  libmesh_assert (_values.size() == _local_size);
  libmesh_assert ((_last_local_index - _first_local_index) == _local_size);
  libmesh_assert (i<size());
  libmesh_assert (i-first_local_index() < local_size());
  
  _values[i - _first_local_index] = value;
}



template <typename T>
inline
void DistributedVector<T>::add (const unsigned int i, const T value)
{
  libmesh_assert (this->initialized());
  libmesh_assert (_values.size() == _local_size);
  libmesh_assert ((_last_local_index - _first_local_index) == _local_size);
  libmesh_assert (i<size());
  libmesh_assert (i-first_local_index() < local_size());
  
  _values[i - _first_local_index] += value;
}



template <typename T>
inline
Real DistributedVector<T>::min () const
{
  // This function must be run on all processors at once
  parallel_only();

  libmesh_assert (this->initialized());
  libmesh_assert (_values.size() == _local_size);
  libmesh_assert ((_last_local_index - _first_local_index) == _local_size);

  Real local_min = _values.size() ? 
    libmesh_real(_values[0]) : std::numeric_limits<Real>::max();
  for (unsigned int i = 1; i < _values.size(); ++i)
    local_min = std::min(libmesh_real(_values[i]), local_min);
  
  Parallel::min(local_min);

  return local_min;
}



template <typename T>
inline
Real DistributedVector<T>::max() const
{
  // This function must be run on all processors at once
  parallel_only();

  libmesh_assert (this->initialized());
  libmesh_assert (_values.size() == _local_size);
  libmesh_assert ((_last_local_index - _first_local_index) == _local_size);

  Real local_max = _values.size() ? 
    libmesh_real(_values[0]) : -std::numeric_limits<Real>::max();
  for (unsigned int i = 1; i < _values.size(); ++i)
    local_max = std::max(libmesh_real(_values[i]), local_max);
  
  Parallel::max(local_max);

  return local_max;
}

#endif  // #ifdef __distributed_vector_h__