petsc_matrix.h

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

};




//-----------------------------------------------------------------------
// PetscMatrix inline members
template <typename T>
inline
PetscMatrix<T>::PetscMatrix()
  : _destroy_mat_on_exit(true)
{}




template <typename T>
inline
PetscMatrix<T>::PetscMatrix(Mat m)
  : _destroy_mat_on_exit(false)
{
  this->_mat = m;
  this->_is_initialized = true;
}




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



template <typename T>
inline
void PetscMatrix<T>::close () const
{
  // BSK - 1/19/2004
  // strictly this check should be OK, but it seems to
  // fail on matrix-free matrices.  Do they falsely
  // state they are assembled?  Check with the developers...
//   if (this->closed())
//     return;
  
  int ierr=0;
 
  ierr = MatAssemblyBegin (_mat, MAT_FINAL_ASSEMBLY);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
  ierr = MatAssemblyEnd   (_mat, MAT_FINAL_ASSEMBLY);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
}



template <typename T>
inline
unsigned int PetscMatrix<T>::m () const
{
  libmesh_assert (this->initialized());
  
  int petsc_m=0, petsc_n=0, ierr=0;

  ierr = MatGetSize (_mat, &petsc_m, &petsc_n);

  return static_cast<unsigned int>(petsc_m);
}



template <typename T>
inline
unsigned int PetscMatrix<T>::n () const
{
  libmesh_assert (this->initialized());
  
  int petsc_m=0, petsc_n=0, ierr=0;

  ierr = MatGetSize (_mat, &petsc_m, &petsc_n);

  return static_cast<unsigned int>(petsc_n);
}



template <typename T>
inline
unsigned int PetscMatrix<T>::row_start () const
{
  libmesh_assert (this->initialized());
  
  int start=0, stop=0, ierr=0;

  ierr = MatGetOwnershipRange(_mat, &start, &stop);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

  return static_cast<unsigned int>(start);
}



template <typename T>
inline
unsigned int PetscMatrix<T>::row_stop () const
{
  libmesh_assert (this->initialized());
  
  int start=0, stop=0, ierr=0;

  ierr = MatGetOwnershipRange(_mat, &start, &stop);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

  return static_cast<unsigned int>(stop);
}



template <typename T>
inline
void PetscMatrix<T>::set (const unsigned int i,
			  const unsigned int j,
			  const T value)
{  
  libmesh_assert (this->initialized());
  
  int ierr=0, i_val=i, j_val=j;

  PetscScalar petsc_value = static_cast<PetscScalar>(value);
  ierr = MatSetValues(_mat, 1, &i_val, 1, &j_val,
		      &petsc_value, INSERT_VALUES);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
}



template <typename T>
inline
void PetscMatrix<T>::add (const unsigned int i,
			  const unsigned int j,
			  const T value)
{
  libmesh_assert (this->initialized());
  
  int ierr=0, i_val=i, j_val=j;

  PetscScalar petsc_value = static_cast<PetscScalar>(value);
  ierr = MatSetValues(_mat, 1, &i_val, 1, &j_val,
		      &petsc_value, ADD_VALUES);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
}



template <typename T>
inline
void PetscMatrix<T>::add_matrix(const DenseMatrix<T>& dm,
				const std::vector<unsigned int>& dof_indices)
{
  this->add_matrix (dm, dof_indices, dof_indices);
}







template <typename T>
inline
void PetscMatrix<T>::add (const T a_in, SparseMatrix<T> &X_in)
{
  libmesh_assert (this->initialized());

  // sanity check. but this cannot avoid 
  // crash due to incompatible sparsity structure...
  libmesh_assert (this->m() == X_in.m());
  libmesh_assert (this->n() == X_in.n());

  PetscScalar     a = static_cast<PetscScalar>      (a_in);
  PetscMatrix<T>* X = dynamic_cast<PetscMatrix<T>*> (&X_in);

  libmesh_assert (X != NULL);
  
  int ierr=0;

  // the matrix from which we copy the values has to be assembled/closed
  X->close ();

// 2.2.x & earlier style
#if PETSC_VERSION_LESS_THAN(2,3,0)  
  
  ierr = MatAXPY(&a,  X->_mat, _mat, SAME_NONZERO_PATTERN);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
	 
// 2.3.x & newer
#else
  
  ierr = MatAXPY(_mat, a, X->_mat, DIFFERENT_NONZERO_PATTERN);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
	 
#endif
}




template <typename T>
inline
T PetscMatrix<T>::operator () (const unsigned int i,
			       const unsigned int j) const
{
  libmesh_assert (this->initialized());

#if PETSC_VERSION_LESS_THAN(2,2,1)

  // PETSc 2.2.0 & older
  PetscScalar *petsc_row;
  int* petsc_cols;
  
#else
  
  // PETSc 2.2.1 & newer
  const PetscScalar *petsc_row;
  const PetscInt    *petsc_cols;

#endif
  
  
  T value=0.;  
  
  int
    ierr=0,
    ncols=0,
    i_val=static_cast<int>(i),
    j_val=static_cast<int>(j);
  

  // the matrix needs to be closed for this to work
  this->close();

  ierr = MatGetRow(_mat, i_val, &ncols, &petsc_cols, &petsc_row);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
	 
  // Perform a binary search to find the contiguous index in
  // petsc_cols (resp. petsc_row) corresponding to global index j_val
  std::pair<const int*, const int*> p =
    std::equal_range (&petsc_cols[0], &petsc_cols[0] + ncols, j_val);

  // Found an entry for j_val
  if (p.first != p.second)
    {
      // The entry in the contiguous row corresponding
      // to the j_val column of interest
      const int j = std::distance (const_cast<int*>(&petsc_cols[0]),
				   const_cast<int*>(p.first));
      
      libmesh_assert (j < ncols);
      libmesh_assert (petsc_cols[j] == j_val);
      
      value = static_cast<T> (petsc_row[j]);
      
      ierr  = MatRestoreRow(_mat, i_val,
			    &ncols, &petsc_cols, &petsc_row);
              CHKERRABORT(libMesh::COMM_WORLD,ierr);
	  
      return value;
    }
  
  // Otherwise the entry is not in the sparse matrix,
  // i.e. it is 0.  
  return 0.;
}




template <typename T>
inline
bool PetscMatrix<T>::closed() const
{
  libmesh_assert (this->initialized());
  
  int ierr=0;
  PetscTruth assembled;

  ierr = MatAssembled(_mat, &assembled);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

  return (assembled == PETSC_TRUE);
}


template <typename T>
inline
void PetscMatrix<T>::swap(PetscMatrix<T> &m) 
{
  std::swap(_mat, m._mat);
  std::swap(_destroy_mat_on_exit, m._destroy_mat_on_exit);
}





template <typename T>
inline
void PetscMatrix<T>::print_personal(std::ostream& os) const
{
  libmesh_assert (this->initialized());

#ifndef NDEBUG
  if (os != std::cout)
    std::cerr << "Warning! PETSc can only print to std::cout!" << std::endl;
#endif
  
  int ierr=0;

  ierr = MatView(_mat, PETSC_VIEWER_STDOUT_SELF);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
}

#endif // #ifdef HAVE_PETSC
#endif // #ifdef __petsc_matrix_h__