trilinos_epetra_matrix.C

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA


// C++ includes
#include "libmesh/libmesh_config.h"

#ifdef LIBMESH_TRILINOS_HAVE_EPETRA

// Local includes
#include "libmesh/trilinos_epetra_matrix.h"
#include "libmesh/trilinos_epetra_vector.h"
#include "libmesh/dof_map.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/parallel.h"
#include "libmesh/sparsity_pattern.h"
#include "libmesh/int_range.h"

namespace libMesh
{



//-----------------------------------------------------------------------
//EpetraMatrix members

template <typename T>
void EpetraMatrix<T>::update_sparsity_pattern (const SparsityPattern::Graph & sparsity_pattern)
{
  // clear data, start over
  this->clear ();

  // big trouble if this fails!
  libmesh_assert(this->_dof_map);

  const numeric_index_type n_rows = cast_int<numeric_index_type>
    (sparsity_pattern.size());

  const numeric_index_type m   = this->_dof_map->n_dofs();
  const numeric_index_type n   = m;
  const numeric_index_type n_l = this->_dof_map->n_local_dofs();
  const numeric_index_type m_l = n_l;

  // error checking
#ifndef NDEBUG
  {
    libmesh_assert_equal_to (n, m);
    libmesh_assert_equal_to (n_l, m_l);

    numeric_index_type
      summed_m_l = m_l,
      summed_n_l = n_l;

    this->comm().sum (summed_m_l);
    this->comm().sum (summed_n_l);

    libmesh_assert_equal_to (m, summed_m_l);
    libmesh_assert_equal_to (n, summed_n_l);
  }
#endif

  // build a map defining the data distribution
  _map = new Epetra_Map (static_cast<int>(m),
                         m_l,
                         0,
                         Epetra_MpiComm (this->comm().get()));

  libmesh_assert_equal_to (static_cast<numeric_index_type>(_map->NumGlobalPoints()), m);
  libmesh_assert_equal_to (static_cast<numeric_index_type>(_map->MaxAllGID()+1), m);

  const std::vector<numeric_index_type> & n_nz = this->_sp->get_n_nz();
  const std::vector<numeric_index_type> & n_oz = this->_sp->get_n_oz();

  // Make sure the sparsity pattern isn't empty
  libmesh_assert_equal_to (n_nz.size(), n_l);
  libmesh_assert_equal_to (n_oz.size(), n_l);

  // Epetra wants the total number of nonzeros, both local and remote.
  std::vector<int> n_nz_tot;  n_nz_tot.reserve(n_nz.size());

  for (auto i : index_range(n_nz))
    n_nz_tot.push_back(std::min(n_nz[i] + n_oz[i], n));

  if (m==0)
    return;

  _graph = new Epetra_CrsGraph(Copy, *_map, n_nz_tot.data());

  // Tell the matrix about its structure.  Initialize it
  // to zero.
  for (numeric_index_type i=0; i<n_rows; i++)
    _graph->InsertGlobalIndices(_graph->GRID(i),
                                cast_int<numeric_index_type>(sparsity_pattern[i].size()),
                                const_cast<int *>(reinterpret_cast<const int *>(sparsity_pattern[i].data())));

  _graph->FillComplete();

  //Initialize the matrix
  libmesh_assert (!this->initialized());
  this->init ();
  libmesh_assert (this->initialized());
}



template <typename T>
void EpetraMatrix<T>::init (const numeric_index_type m,
                            const numeric_index_type n,
                            const numeric_index_type m_l,
                            const numeric_index_type libmesh_dbg_var(n_l),
                            const numeric_index_type nnz,
                            const numeric_index_type noz,
                            const numeric_index_type /* blocksize */)
{
  if ((m==0) || (n==0))
    return;

  {
    // Clear initialized matrices
    if (this->initialized())
      this->clear();

    libmesh_assert (!this->_mat);
    libmesh_assert (!this->_map);

    this->_is_initialized = true;
  }

  // error checking
#ifndef NDEBUG
  {
    libmesh_assert_equal_to (n, m);
    libmesh_assert_equal_to (n_l, m_l);

    numeric_index_type
      summed_m_l = m_l,
      summed_n_l = n_l;

    this->comm().sum (summed_m_l);
    this->comm().sum (summed_n_l);

    libmesh_assert_equal_to (m, summed_m_l);
    libmesh_assert_equal_to (n, summed_n_l);
  }
#endif

  // build a map defining the data distribution
  _map = new Epetra_Map (static_cast<int>(m),
                         m_l,
                         0,
                         Epetra_MpiComm (this->comm().get()));

  libmesh_assert_equal_to (static_cast<numeric_index_type>(_map->NumGlobalPoints()), m);
  libmesh_assert_equal_to (static_cast<numeric_index_type>(_map->MaxAllGID()+1), m);

  _mat = new Epetra_FECrsMatrix (Copy, *_map, nnz + noz);
}




template <typename T>
void EpetraMatrix<T>::init (const ParallelType)
{
  libmesh_assert(this->_dof_map);

  {
    // Clear initialized matrices
    if (this->initialized())
      this->clear();

    this->_is_initialized = true;
  }


  _mat = new Epetra_FECrsMatrix (Copy, *_graph);
}



template <typename T>
void EpetraMatrix<T>::zero ()
{
  libmesh_assert (this->initialized());

  _mat->Scale(0.0);
}



template <typename T>
std::unique_ptr<SparseMatrix<T>> EpetraMatrix<T>::zero_clone () const
{
  // This function is marked as "not implemented" since it hasn't been
  // tested, the code below might serve as a possible implementation.
  libmesh_not_implemented();

  // Make empty copy with matching comm, initialize, and return.
  auto mat_copy = std::make_unique<EpetraMatrix<T>>(this->comm());
  mat_copy->init();
  mat_copy->zero();

  return mat_copy;
}



template <typename T>
std::unique_ptr<SparseMatrix<T>> EpetraMatrix<T>::clone () const
{
  // We don't currently have a faster implementation than making a
  // zero clone and then filling in the values.
  auto mat_copy = this->zero_clone();
  mat_copy->add(1., *this);

  return mat_copy;
}



template <typename T>
void EpetraMatrix<T>::clear () noexcept
{
  // FIXME: clear() doesn't actually free the memory managed by this
  // class, so it probably leaks memory.
  // delete _mat;
  // delete _map;

  this->_is_initialized = false;
}



template <typename T>
Real EpetraMatrix<T>::l1_norm () const
{
  libmesh_assert (this->initialized());

  libmesh_assert(_mat);

  return static_cast<Real>(_mat->NormOne());
}



template <typename T>
Real EpetraMatrix<T>::linfty_norm () const
{
  libmesh_assert (this->initialized());


  libmesh_assert(_mat);

  return static_cast<Real>(_mat->NormInf());
}



template <typename T>
void EpetraMatrix<T>::add_matrix(const DenseMatrix<T> & dm,
                                 const std::vector<numeric_index_type> & rows,
                                 const std::vector<numeric_index_type> & cols)
{
  libmesh_assert (this->initialized());

  const numeric_index_type m = dm.m();
  const numeric_index_type n = dm.n();

  libmesh_assert_equal_to (rows.size(), m);
  libmesh_assert_equal_to (cols.size(), n);

  _mat->SumIntoGlobalValues(m, numeric_trilinos_cast(rows.data()),
                            n, numeric_trilinos_cast(cols.data()),
                            dm.get_values().data());
}






template <typename T>
void EpetraMatrix<T>::get_diagonal (NumericVector<T> & dest) const
{
  // Convert vector to EpetraVector.
  EpetraVector<T> * epetra_dest = cast_ptr<EpetraVector<T> *>(&dest);

  // Call Epetra function.
  _mat->ExtractDiagonalCopy(*(epetra_dest->vec()));
}



template <typename T>
void EpetraMatrix<T>::get_transpose (SparseMatrix<T> & dest) const
{
  // Make sure the SparseMatrix passed in is really a EpetraMatrix
  EpetraMatrix<T> & epetra_dest = cast_ref<EpetraMatrix<T> &>(dest);

  // We currently only support calling get_transpose() with ourself
  // as the destination. Previously, this called the default copy
  // constructor which was not safe because this class manually
  // manages memory.
  if (&epetra_dest != this)
    libmesh_not_implemented();

  epetra_dest._use_transpose = !epetra_dest._use_transpose;
  epetra_dest._mat->SetUseTranspose(epetra_dest._use_transpose);
}



template <typename T>
void EpetraMatrix<T>::get_row(numeric_index_type i,
                              std::vector<numeric_index_type> & indices,
                              std::vector<T> & values) const
{
  libmesh_assert (this->initialized());
  libmesh_assert(this->_mat);
  libmesh_assert (this->_mat->MyGlobalRow(static_cast<int>(i)));
  libmesh_assert_greater_equal (i, this->row_start());
  libmesh_assert_less (i, this->row_stop());

  int row_length;
  int * row_indices;
  double * row_values;

  _mat->ExtractMyRowView (i-this->row_start(),
                          row_length,
                          row_values,
                          row_indices);

  indices.resize(row_length);
  values.resize(row_length);

  for (auto i : make_range(row_length))
    {
      indices[i] = row_indices[i];
      values[i] = row_values[i];
    }
}



template <typename T>
EpetraMatrix<T>::EpetraMatrix(const Parallel::Communicator & comm) :
  SparseMatrix<T>(comm),
  _destroy_mat_on_exit(true),
  _use_transpose(false)
{}




template <typename T>
EpetraMatrix<T>::EpetraMatrix(Epetra_FECrsMatrix * m,
                              const Parallel::Communicator & comm) :
  SparseMatrix<T>(comm),
  _destroy_mat_on_exit(false),
  _use_transpose(false) // dumb guess is the best we can do...
{
  this->_mat = m;
  this->_is_initialized = true;
}




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



template <typename T>
void EpetraMatrix<T>::close ()
{
  libmesh_assert(_mat);

  _mat->GlobalAssemble();
}



template <typename T>
numeric_index_type EpetraMatrix<T>::m () const
{
  libmesh_assert (this->initialized());

  return static_cast<numeric_index_type>(_mat->NumGlobalRows());
}



template <typename T>
numeric_index_type EpetraMatrix<T>::n () const
{
  libmesh_assert (this->initialized());

  return static_cast<numeric_index_type>(_mat->NumGlobalCols());
}



template <typename T>
numeric_index_type EpetraMatrix<T>::row_start () const
{
  libmesh_assert (this->initialized());
  libmesh_assert(_map);

  return static_cast<numeric_index_type>(_map->MinMyGID());
}



template <typename T>
numeric_index_type EpetraMatrix<T>::row_stop () const
{
  libmesh_assert (this->initialized());
  libmesh_assert(_map);

  return static_cast<numeric_index_type>(_map->MaxMyGID())+1;
}



template <typename T>
numeric_index_type EpetraMatrix<T>::col_start () const
{
  libmesh_assert (this->initialized());
  libmesh_assert(_map);

  return static_cast<numeric_index_type>(_map->MinMyGID());
}



template <typename T>
numeric_index_type EpetraMatrix<T>::col_stop () const
{
  libmesh_assert (this->initialized());
  libmesh_assert(_map);

  return static_cast<numeric_index_type>(_map->MaxMyGID())+1;
}



template <typename T>
void EpetraMatrix<T>::set (const numeric_index_type i,
                           const numeric_index_type j,
                           const T value)
{
  libmesh_assert (this->initialized());

  int
    epetra_i = static_cast<int>(i),
    epetra_j = static_cast<int>(j);

  T epetra_value = value;

  if (_mat->Filled())
    _mat->ReplaceGlobalValues (epetra_i, 1, &epetra_value, &epetra_j);
  else
    _mat->InsertGlobalValues (epetra_i, 1, &epetra_value, &epetra_j);
}



template <typename T>
void EpetraMatrix<T>::add (const numeric_index_type i,
                           const numeric_index_type j,
                           const T value)
{
  libmesh_assert (this->initialized());

  int
    epetra_i = static_cast<int>(i),
    epetra_j = static_cast<int>(j);

  T epetra_value = value;

  _mat->SumIntoGlobalValues (epetra_i, 1, &epetra_value, &epetra_j);
}



template <typename T>
void EpetraMatrix<T>::add_matrix(const DenseMatrix<T> & dm,
                                 const std::vector<numeric_index_type> & dof_indices)
{
  this->add_matrix (dm, dof_indices, dof_indices);
}



template <typename T>
void EpetraMatrix<T>::add (const T a_in, const SparseMatrix<T> & X_in)
{
#ifdef LIBMESH_TRILINOS_HAVE_EPETRAEXT
  libmesh_assert (this->initialized());

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

  const EpetraMatrix<T> * X =
    cast_ptr<const EpetraMatrix<T> *> (&X_in);

  EpetraExt::MatrixMatrix::Add (*X->_mat, false, a_in, *_mat, 1.);
#else
  libmesh_error_msg("ERROR: EpetraExt is required for EpetraMatrix::add()!");
#endif
}




template <typename T>
T EpetraMatrix<T>::operator () (const numeric_index_type i,
                                const numeric_index_type j) const
{
  libmesh_assert (this->initialized());
  libmesh_assert(this->_mat);
  libmesh_assert (this->_mat->MyGlobalRow(static_cast<int>(i)));
  libmesh_assert_greater_equal (i, this->row_start());
  libmesh_assert_less (i, this->row_stop());


  int row_length;
  int * row_indices;
  double * values;

  _mat->ExtractMyRowView (i-this->row_start(),
                          row_length,
                          values,
                          row_indices);

  //libMesh::out << "row_length=" << row_length << std::endl;

  int * index = std::lower_bound (row_indices, row_indices+row_length, j);

  libmesh_assert_less (*index, row_length);
  libmesh_assert_equal_to (static_cast<numeric_index_type>(row_indices[*index]), j);

  //libMesh::out << "val=" << values[*index] << std::endl;

  return values[*index];
}




template <typename T>
bool EpetraMatrix<T>::closed() const
{
  libmesh_assert (this->initialized());
  libmesh_assert(this->_mat);

  return this->_mat->Filled();
}


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





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

  os << *_mat;
}



//------------------------------------------------------------------
// Explicit instantiations
template class LIBMESH_EXPORT EpetraMatrix<Number>;

} // namespace libMesh


#endif // LIBMESH_TRILINOS_HAVE_EPETRA