trilinos_epetra_vector.c

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

    
    for(unsigned int i=0;i<nl;i++)
      values[i]=v[fli+i];
  }

  /**
   * Case 2: The vector is the same size as our local
   * piece.  Insert directly to the local piece.
   */
  else
  {
    libmesh_assert(v.size()==this->local_size());
    
    const unsigned int nl=this->local_size();
    
    for(unsigned int i=0;i<nl;i++)
      values[i]=v[i];
  }

  return *this;
}



template <typename T>
void EpetraVector<T>::localize (NumericVector<T>& v_local_in) const
{
  EpetraVector<T>* v_local = dynamic_cast<EpetraVector<T>*>(&v_local_in);

  libmesh_assert(v_local);

  Epetra_Map rootMap = Epetra_Util::Create_Root_Map( *_map, libMesh::processor_id() );
  Epetra_Import importer(rootMap, *_map);
  v_local->_vec->ReplaceMap(rootMap);
  
  v_local->_vec->Import(*_vec, importer, Insert);
}



template <typename T>
void EpetraVector<T>::localize (NumericVector<T>& v_local_in,
				const std::vector<unsigned int>& /* send_list */) const
{
  // TODO: optimize to sync only the send list values
  this->localize(v_local_in);
  
//   EpetraVector<T>* v_local =
//   dynamic_cast<EpetraVector<T>*>(&v_local_in);

//   libmesh_assert (v_local != NULL);
//   libmesh_assert (this->_map.get() != NULL);
//   libmesh_assert (v_local->_map.get() != NULL);
//   libmesh_assert (v_local->local_size() == this->size());
//   libmesh_assert (send_list.size() <= v_local->size());

//   Epetra_Import importer (*v_local->_map, *this->_map);
  
//   v_local->_vec->Import (*this->_vec, importer, Insert);
}


template <typename T>
void EpetraVector<T>::localize (const unsigned int /* first_local_idx */,
				const unsigned int /* last_local_idx */,
				const std::vector<unsigned int>& /* send_list */)
{
  LIBMESH_THROW(libMesh::NotImplemented());

//   // Only good for serial vectors.
//   libmesh_assert (this->size() == this->local_size());
//   libmesh_assert (last_local_idx > first_local_idx);
//   libmesh_assert (send_list.size() <= this->size());
//   libmesh_assert (last_local_idx < this->size());
  
//   const unsigned int size       = this->size();
//   const unsigned int local_size = (last_local_idx - first_local_idx + 1);
//   int ierr=0;  
  
//   // Don't bother for serial cases
//   if ((first_local_idx == 0) &&
//       (local_size == size))
//     return;
  
  
//   // Build a parallel vector, initialize it with the local
//   // parts of (*this)
//   EpetraVector<T> parallel_vec;

//   parallel_vec.init (size, local_size);


//   // Copy part of *this into the parallel_vec
//   {
//     IS is;
//     VecScatter scatter;

//     // Create idx, idx[i] = i+first_local_idx;
//     std::vector<int> idx(local_size);
//     Utility::iota (idx.begin(), idx.end(), first_local_idx);

//     // Create the index set & scatter object
//     ierr = ISCreateGeneral(libMesh::COMM_WORLD, local_size, &idx[0], &is); 
//            CHKERRABORT(libMesh::COMM_WORLD,ierr);

//     ierr = VecScatterCreate(_vec,              is,
// 			    parallel_vec._vec, is,
// 			    &scatter);
//            CHKERRABORT(libMesh::COMM_WORLD,ierr);

//     // Perform the scatter
// #if EPETRA_VERSION_LESS_THAN(2,3,3)

//     ierr = VecScatterBegin(_vec, parallel_vec._vec, INSERT_VALUES,
// 			   SCATTER_FORWARD, scatter);
//            CHKERRABORT(libMesh::COMM_WORLD,ierr);
  
//     ierr = VecScatterEnd  (_vec, parallel_vec._vec, INSERT_VALUES,
// 			   SCATTER_FORWARD, scatter);
//            CHKERRABORT(libMesh::COMM_WORLD,ierr);

// #else
	   
//       // API argument order change in Epetra 2.3.3
//     ierr = VecScatterBegin(scatter, _vec, parallel_vec._vec,
// 			   INSERT_VALUES, SCATTER_FORWARD);
//            CHKERRABORT(libMesh::COMM_WORLD,ierr);
  
//     ierr = VecScatterEnd  (scatter, _vec, parallel_vec._vec,
// 			   INSERT_VALUES, SCATTER_FORWARD);
//            CHKERRABORT(libMesh::COMM_WORLD,ierr);
	   
// #endif

//     // Clean up
//     ierr = ISDestroy (is);
//            CHKERRABORT(libMesh::COMM_WORLD,ierr);
  
//     ierr = VecScatterDestroy(scatter);
//            CHKERRABORT(libMesh::COMM_WORLD,ierr);
//   }

//   // localize like normal
//   parallel_vec.close();
//   parallel_vec.localize (*this, send_list);
//   this->close();
}



template <typename T>
void EpetraVector<T>::localize (std::vector<T>& v_local) const
{
  // This function must be run on all processors at once
  parallel_only();

  const unsigned int n  = this->size();
  const unsigned int nl = this->local_size();

  libmesh_assert (this->_vec.get() != NULL);

  v_local.clear();
  v_local.reserve(n);

  // build up my local part
  for (unsigned int i=0; i<nl; i++)
    v_local.push_back((*this->_vec)[0][i]);

  Parallel::allgather (v_local);
}



template <typename T>
void EpetraVector<T>::localize_to_one (std::vector<T>&  v_local,
				       const unsigned int  pid) const
{
  // This function must be run on all processors at once
  parallel_only();

  const unsigned int n  = this->size();
  const unsigned int nl = this->local_size();

  libmesh_assert (pid < libMesh::n_processors());
  libmesh_assert (this->_vec.get() != NULL);
  
  v_local.clear();
  v_local.reserve(n);
  
  
  // build up my local part
  for (unsigned int i=0; i<nl; i++)
    v_local.push_back((*this->_vec)[0][i]);

  Parallel::gather (pid, v_local);
}



template <typename T>
void EpetraVector<T>::print_matlab (const std::string /* name */) const
{
  LIBMESH_THROW(libMesh::NotImplemented());

//   libmesh_assert (this->initialized());
//   libmesh_assert (this->closed());
  
//   int ierr=0; 
//   EpetraViewer epetra_viewer;


//   ierr = EpetraViewerCreate (libMesh::COMM_WORLD,
// 			    &epetra_viewer);
//          CHKERRABORT(libMesh::COMM_WORLD,ierr);

//   /**
//    * Create an ASCII file containing the matrix
//    * if a filename was provided.  
//    */
//   if (name != "NULL")
//     {
//       ierr = EpetraViewerASCIIOpen( libMesh::COMM_WORLD,
// 				   name.c_str(),
// 				   &epetra_viewer);
//              CHKERRABORT(libMesh::COMM_WORLD,ierr);
      
//       ierr = EpetraViewerSetFormat (epetra_viewer,
// 				   EPETRA_VIEWER_ASCII_MATLAB);
//              CHKERRABORT(libMesh::COMM_WORLD,ierr);
  
//       ierr = VecView (_vec, epetra_viewer);
//              CHKERRABORT(libMesh::COMM_WORLD,ierr);
//     }

//   /**
//    * Otherwise the matrix will be dumped to the screen.
//    */
//   else
//     {
//       ierr = EpetraViewerSetFormat (EPETRA_VIEWER_STDOUT_WORLD,
// 				   EPETRA_VIEWER_ASCII_MATLAB);
//              CHKERRABORT(libMesh::COMM_WORLD,ierr);
  
//       ierr = VecView (_vec, EPETRA_VIEWER_STDOUT_WORLD);
//              CHKERRABORT(libMesh::COMM_WORLD,ierr);
//     }


//   /**
//    * Destroy the viewer.
//    */
//   ierr = EpetraViewerDestroy (epetra_viewer);
//          CHKERRABORT(libMesh::COMM_WORLD,ierr);
}





template <typename T>
void EpetraVector<T>::create_subvector(NumericVector<T>& /* subvector */,
				       const std::vector<unsigned int>& /* rows */) const
{
  LIBMESH_THROW(libMesh::NotImplemented());

//   // Epetra data structures
//   IS parent_is, subvector_is;
//   VecScatter scatter;
//   int ierr = 0;
  
//   // Make sure the passed int subvector is really a EpetraVector
//   EpetraVector<T>* epetra_subvector = dynamic_cast<EpetraVector<T>*>(&subvector);
//   libmesh_assert(epetra_subvector != NULL);
  
//   // If the epetra_subvector is already initialized, we assume that the
//   // user has already allocated the *correct* amount of space for it.
//   // If not, we use the appropriate Epetra routines to initialize it.
//   if (!epetra_subvector->initialized())
//     {
//       // Initialize the epetra_subvector to have enough space to hold
//       // the entries which will be scattered into it.  Note: such an
//       // init() function (where we let Epetra decide the number of local
//       // entries) is not currently offered by the EpetraVector
//       // class.  Should we differentiate here between sequential and
//       // parallel vector creation based on libMesh::n_processors() ?
//       ierr = VecCreateMPI(libMesh::COMM_WORLD,
// 			  EPETRA_DECIDE,          // n_local
// 			  rows.size(),           // n_global
// 			  &(epetra_subvector->_vec)); CHKERRABORT(libMesh::COMM_WORLD,ierr);

//       ierr = VecSetFromOptions (epetra_subvector->_vec); CHKERRABORT(libMesh::COMM_WORLD,ierr);

//       // Mark the subvector as initialized
//       epetra_subvector->_is_initialized = true;
//     }
  
//   // Use iota to fill an array with entries [0,1,2,3,4,...rows.size()]
//   std::vector<int> idx(rows.size());
//   Utility::iota (idx.begin(), idx.end(), 0);

//   // Construct index sets
//   ierr = ISCreateGeneral(libMesh::COMM_WORLD,
// 			 rows.size(),
// 			 (int*) &rows[0],
// 			 &parent_is); CHKERRABORT(libMesh::COMM_WORLD,ierr);

//   ierr = ISCreateGeneral(libMesh::COMM_WORLD,
// 			 rows.size(),
// 			 (int*) &idx[0],
// 			 &subvector_is); CHKERRABORT(libMesh::COMM_WORLD,ierr);

//   // Construct the scatter object
//   ierr = VecScatterCreate(this->_vec,
// 			  parent_is,
// 			  epetra_subvector->_vec,
// 			  subvector_is,
// 			  &scatter); CHKERRABORT(libMesh::COMM_WORLD,ierr);

//   // Actually perform the scatter
// #if EPETRA_VERSION_LESS_THAN(2,3,3)
//   ierr = VecScatterBegin(this->_vec,
// 			 epetra_subvector->_vec,
// 			 INSERT_VALUES,
// 			 SCATTER_FORWARD,
// 			 scatter); CHKERRABORT(libMesh::COMM_WORLD,ierr);

//   ierr = VecScatterEnd(this->_vec,
// 		       epetra_subvector->_vec,
// 		       INSERT_VALUES,
// 		       SCATTER_FORWARD,
// 		       scatter); CHKERRABORT(libMesh::COMM_WORLD,ierr);
// #else
//   // API argument order change in Epetra 2.3.3
//   ierr = VecScatterBegin(scatter,
// 			 this->_vec,
// 			 epetra_subvector->_vec,
// 			 INSERT_VALUES,
// 			 SCATTER_FORWARD); CHKERRABORT(libMesh::COMM_WORLD,ierr);

//   ierr = VecScatterEnd(scatter,
// 		       this->_vec,
// 		       epetra_subvector->_vec,
// 		       INSERT_VALUES,
// 		       SCATTER_FORWARD); CHKERRABORT(libMesh::COMM_WORLD,ierr);
// #endif
  
//   // Clean up 
//   ierr = ISDestroy(parent_is);       CHKERRABORT(libMesh::COMM_WORLD,ierr);
//   ierr = ISDestroy(subvector_is);    CHKERRABORT(libMesh::COMM_WORLD,ierr);
//   ierr = VecScatterDestroy(scatter); CHKERRABORT(libMesh::COMM_WORLD,ierr); 
}




//------------------------------------------------------------------
// Explicit instantiations
template class EpetraVector<Number>;

#endif // #ifdef HAVE_EPETRA