laspack_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



// Local includes
#include "libmesh/libmesh_config.h"

#ifdef LIBMESH_HAVE_LASPACK

#include "libmesh/laspack_matrix.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dof_map.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/sparsity_pattern.h"


// C++ Includes
#include <memory>


namespace libMesh
{


//-----------------------------------------------------------------------
// LaspackMatrix members
template <typename T>
void LaspackMatrix<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 = sparsity_pattern.size();

  // Initialize the _row_start data structure,
  // allocate storage for the _csr array
  {
    std::size_t size = 0;

    for (numeric_index_type row=0; row<n_rows; row++)
      size += sparsity_pattern[row].size();

    _csr.resize       (size);
    _row_start.reserve(n_rows + 1);
  }


  // Initialize the _csr data structure.
  {
    std::vector<numeric_index_type>::iterator position = _csr.begin();

    _row_start.push_back (position);

    for (numeric_index_type row=0; row<n_rows; row++)
      {
        // insert the row indices
        for (const auto & col : sparsity_pattern[row])
          {
            libmesh_assert (position != _csr.end());
            *position = col;
            ++position;
          }

        _row_start.push_back (position);
      }
  }


  // Initialize the matrix
  libmesh_assert (!this->initialized());
  this->init ();
  libmesh_assert (this->initialized());
  //libMesh::out << "n_rows=" << n_rows << std::endl;
  //libMesh::out << "m()=" << m() << std::endl;
  libmesh_assert_equal_to (n_rows, this->m());

  // Tell the matrix about its structure.  Initialize it
  // to zero.
  for (numeric_index_type i=0; i<n_rows; i++)
    {
      auto rs = _row_start[i];

      const numeric_index_type length = _row_start[i+1] - rs;

      Q_SetLen (&_QMat, i+1, length);

      for (numeric_index_type l=0; l<length; l++)
        {
          const numeric_index_type j = *(rs+l);

          // sanity check
          //libMesh::out << "m()=" << m() << std::endl;
          //libMesh::out << "(i,j,l) = (" << i
          //          << "," << j
          //          << "," << l
          //           << ")" << std::endl;
          //libMesh::out << "pos(i,j)=" << pos(i,j)
          //              << std::endl;
          libmesh_assert_equal_to (this->pos(i,j), l);
          Q_SetEntry (&_QMat, i+1, l, j+1, 0.);
        }
    }

  // That's it!
  //libmesh_here();
}



template <typename T>
void LaspackMatrix<T>::init (const numeric_index_type m_in,
                             const numeric_index_type n_in,
                             const numeric_index_type m_l,
                             const numeric_index_type n_l,
                             const numeric_index_type nnz,
                             const numeric_index_type,
                             const numeric_index_type)
{
  // noz ignored...  only used for multiple processors!
  libmesh_assert_equal_to (m_in, m_l);
  libmesh_assert_equal_to (n_in, n_l);
  libmesh_assert_equal_to (m_in, n_in);
  libmesh_assert_greater (nnz, 0);

  if ((m_in != m_l) ||
      (n_in != n_l))
    libmesh_not_implemented_msg("Laspack does not support distributed matrices");

  if (m_in != n_in)
    libmesh_not_implemented_msg("Laspack does not support rectangular matrices");

  if (nnz < n_in)
    libmesh_warning("Using inefficient LaspackMatrix allocation via this init() method");

  const dof_id_type n_rows = m_in;

  // Laspack doesn't let us assign sparse indices on the fly, so
  // without a sparsity pattern we're stuck allocating a dense matrix
  {
    _csr.resize       (m_in * m_in);
    _row_start.reserve(m_in + 1);
  }

  // Initialize the _csr data structure.
  {
    std::vector<numeric_index_type>::iterator position = _csr.begin();

    _row_start.push_back (position);

    for (numeric_index_type row=0; row<n_rows; row++)
      {
        // insert the row indices
        for (const auto & col : make_range(n_rows))
          {
            libmesh_assert (position != _csr.end());
            *position = col;
            ++position;
          }

        _row_start.push_back (position);
      }
  }

  Q_Constr(&_QMat, const_cast<char *>("Mat"), m_in, _LPFalse, Rowws, Normal, _LPTrue);

  this->_is_initialized = true;

  // Tell the matrix about its structure.  Initialize it
  // to zero.
  for (numeric_index_type i=0; i<n_rows; i++)
    {
      auto rs = _row_start[i];

      const numeric_index_type length = _row_start[i+1] - rs;

      Q_SetLen (&_QMat, i+1, length);

      for (numeric_index_type l=0; l<length; l++)
        {
          const numeric_index_type j = *(rs+l);

          libmesh_assert_equal_to (this->pos(i,j), l);
          Q_SetEntry (&_QMat, i+1, l, j+1, 0.);
        }
    }
}



template <typename T>
void LaspackMatrix<T>::init (const ParallelType)
{
  // Ignore calls on initialized objects
  if (this->initialized())
    return;

  // We need a sparsity pattern for this!
  libmesh_assert(this->_sp);

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

  const numeric_index_type n_rows   = this->_dof_map->n_dofs();
#ifndef NDEBUG
  // The following variables are only used for assertions,
  // so avoid declaring them when asserts are inactive.
  const numeric_index_type n_cols   = n_rows;
  const numeric_index_type n_l = this->_dof_map->n_dofs_on_processor(0);
  const numeric_index_type m_l = n_l;
#endif

  // Laspack Matrices only work for uniprocessor cases
  libmesh_assert_equal_to (n_rows, n_cols);
  libmesh_assert_equal_to (m_l, n_rows);
  libmesh_assert_equal_to (n_l, n_cols);

#ifndef NDEBUG
  // The following variables are only used for assertions,
  // so avoid declaring them when asserts are inactive.
  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();
#endif

  // 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);

  if (n_rows==0)
    return;

  Q_Constr(&_QMat, const_cast<char *>("Mat"), n_rows, _LPFalse, Rowws, Normal, _LPTrue);

  this->_is_initialized = true;

  libmesh_assert_equal_to (n_rows, this->m());
}



template <typename T>
void LaspackMatrix<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());
  unsigned int n_rows = cast_int<unsigned int>(rows.size());
  unsigned int n_cols = cast_int<unsigned int>(cols.size());
  libmesh_assert_equal_to (dm.m(), n_rows);
  libmesh_assert_equal_to (dm.n(), n_cols);


  for (unsigned int i=0; i<n_rows; i++)
    for (unsigned int j=0; j<n_cols; j++)
      this->add(rows[i],cols[j],dm(i,j));
}



template <typename T>
Real LaspackMatrix<T>::l1_norm () const
{
  // There does not seem to be a straightforward way to iterate over
  // the columns of a LaspackMatrix.  So we use some extra storage and
  // keep track of the column sums while going over the row entries...
  std::vector<Real> abs_col_sums(this->n());

  const numeric_index_type n_rows = this->m();

  for (numeric_index_type row : make_range(n_rows))
    {
      auto r_start = _row_start[row];

      const numeric_index_type len = (_row_start[row+1] - _row_start[row]);

      // Make sure we agree on the row length
      libmesh_assert_equal_to (len, Q_GetLen(&_QMat, row+1));

      for (numeric_index_type l=0; l<len; l++)
        {
          const numeric_index_type j = *(r_start + l);

          // Make sure the data structures are working
          libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, row+1, l));

          abs_col_sums[j] += std::abs(Q_GetVal (&_QMat, row+1, l));
        }
    }

  return *(std::max_element(abs_col_sums.begin(), abs_col_sums.end()));
}



template <typename T>
Real LaspackMatrix<T>::linfty_norm () const
{
  const numeric_index_type n_rows = this->m();

  Real max_row_sum = 0;

  for (numeric_index_type row : make_range(n_rows))
    {
      Real row_sum = 0;

      const numeric_index_type len = (_row_start[row+1] - _row_start[row]);

      // Make sure we agree on the row length
      libmesh_assert_equal_to (len, Q_GetLen(&_QMat, row+1));

      for (numeric_index_type l=0; l<len; l++)
        {
          // Make sure the data structures are working
          libmesh_assert_equal_to ((*(_row_start[row] + l)+1), Q_GetPos
                                   (&_QMat, row+1, l));

          row_sum += std::abs(Q_GetVal (&_QMat, row+1, l));
        }

      max_row_sum = std::max(max_row_sum, row_sum);
    }

  return max_row_sum;
}



template <typename T>
void LaspackMatrix<T>::get_diagonal (NumericVector<T> & dest) const
{
  const numeric_index_type n_rows = this->m();

  dest.init(n_rows, n_rows);

  for (numeric_index_type i : make_range(n_rows))
    dest.set(i, Q_GetEl (const_cast<QMatrix*>(&_QMat), i+1, i+1));
}



template <typename T>
void LaspackMatrix<T>::get_transpose (SparseMatrix<T> & dest) const
{
  LaspackMatrix<T> & target = cast_ref<LaspackMatrix<T> &>(dest);
  target.clear();

  const numeric_index_type N = this->n();
  libmesh_assert_equal_to(N, this->m());

  std::vector<numeric_index_type> col_sizes(N);
  for (numeric_index_type row : make_range(N))
    {
      const numeric_index_type len = (_row_start[row+1] - _row_start[row]);
      auto r_start = _row_start[row];
      for (numeric_index_type l=0; l<len; l++)
        {
          const numeric_index_type col = *(r_start + l);
          ++col_sizes[col];
        }
    }

  target._csr.resize(_csr.size());
  target._row_start.resize(N+1);
  std::vector<std::vector<numeric_index_type>::iterator> writeable_row_start(N+1);
  writeable_row_start[0] = target._csr.begin();
  target._row_start[0] = target._csr.begin();
  for (numeric_index_type col: make_range(N))
    {
      writeable_row_start[col+1] = writeable_row_start[col] + col_sizes[col];
      target._row_start[col+1] = target._row_start[col] + col_sizes[col];
    }

  std::vector<numeric_index_type> col_offsets;
  // Reuse memory
  col_offsets.swap(col_sizes);
  std::fill(col_offsets.begin(), col_offsets.end(), 0);

  Q_Constr(&target._QMat, const_cast<char *>("Mat"), N, _LPFalse, Rowws, Normal, _LPTrue);
  target._is_initialized = true;

  for (numeric_index_type row : make_range(N))
    {
      const numeric_index_type len = (_row_start[row+1] - _row_start[row]);
      auto r_start = _row_start[row];
      for (numeric_index_type l=0; l<len; l++)
        {
          const numeric_index_type col = *(r_start + l);
          *(writeable_row_start[col]+col_offsets[col]) = row;
          ++col_offsets[col];
        }
    }

  for (numeric_index_type col : make_range(N))
    {
      auto len = col_offsets[col];
      libmesh_assert(target._row_start[col] + len ==
                     target._row_start[col+1]);
      Q_SetLen(&target._QMat, col+1, len);

      auto rs = target._row_start[col];
      for (numeric_index_type l : make_range(len))
        {
          const numeric_index_type j = *(rs+l);
          const auto value = Q_GetEl(const_cast<QMatrix*>(&_QMat), j+1, col+1);
          Q_SetEntry (&target._QMat, col+1, l, j+1, value);
        }
    }
}



template <typename T>
void LaspackMatrix<T>::get_row(numeric_index_type i,
                               std::vector<numeric_index_type> & indices,
                               std::vector<T> & values) const
{
  indices.clear();
  values.clear();

  const numeric_index_type len = (_row_start[i+1] - _row_start[i]);

  // Make sure we agree on the row length
  libmesh_assert_equal_to (len, Q_GetLen(&_QMat, i+1));

  auto r_start = _row_start[i];

  for (numeric_index_type l=0; l<len; l++)
    {
      const numeric_index_type j = *(r_start + l);

      // Make sure the data structures are working
      libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, i+1, l));

      indices.push_back(j);

      values.push_back(Q_GetVal (&_QMat, i+1, l));
    }
}



template <typename T>
LaspackMatrix<T>::LaspackMatrix (const Parallel::Communicator & comm) :
  SparseMatrix<T>(comm),
  _closed (false)
{
}



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



template <typename T>
void LaspackMatrix<T>::clear ()
{
  if (this->initialized())
    {
      Q_Destr(&_QMat);
    }

  _csr.clear();
  _row_start.clear();
  _closed = false;
  this->_is_initialized = false;
}



template <typename T>
void LaspackMatrix<T>::zero ()
{
  const numeric_index_type n_rows = this->m();

  for (numeric_index_type row=0; row<n_rows; row++)
    {
      auto r_start = _row_start[row];

      const numeric_index_type len = (_row_start[row+1] - _row_start[row]);

      // Make sure we agree on the row length
      libmesh_assert_equal_to (len, Q_GetLen(&_QMat, row+1));

      for (numeric_index_type l=0; l<len; l++)
        {
          const numeric_index_type j = *(r_start + l);

          // Make sure the data structures are working
          libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, row+1, l));

          Q_SetEntry (&_QMat, row+1, l, j+1, 0.);
        }
    }

  this->close();
}



template <typename T>
std::unique_ptr<SparseMatrix<T>> LaspackMatrix<T>::zero_clone () const
{
  // Make empty copy with matching comm, initialize, zero, and return.
  auto mat_copy = std::make_unique<LaspackMatrix<T>>(this->comm());
  if (this->_dof_map)
    mat_copy->attach_dof_map(*this->_dof_map);
  else
    mat_copy->init(this->m(), this->n(), this->local_m(),
                   this->local_n(), this->local_n());

  mat_copy->zero();

  return mat_copy;
}



template <typename T>
std::unique_ptr<SparseMatrix<T>> LaspackMatrix<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>
numeric_index_type LaspackMatrix<T>::m () const
{
  libmesh_assert (this->initialized());

  return static_cast<numeric_index_type>(Q_GetDim(const_cast<QMatrix*>(&_QMat)));
}



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

  return static_cast<numeric_index_type>(Q_GetDim(const_cast<QMatrix*>(&_QMat)));
}



template <typename T>
numeric_index_type LaspackMatrix<T>::row_start () const
{
  return 0;
}



template <typename T>
numeric_index_type LaspackMatrix<T>::row_stop () const
{
  return this->m();
}



template <typename T>
numeric_index_type LaspackMatrix<T>::col_start () const
{
  return 0;
}



template <typename T>
numeric_index_type LaspackMatrix<T>::col_stop () const
{
  return this->n();
}



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

  const numeric_index_type position = this->pos(i,j);

  // Sanity check
  libmesh_assert_equal_to (*(_row_start[i]+position), j);
  libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, i+1, position));

  Q_SetEntry (&_QMat, i+1, position, j+1, value);
}



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

  const numeric_index_type position = this->pos(i,j);

  // Sanity check
  libmesh_assert_equal_to (*(_row_start[i]+position), j);

  Q_AddVal (&_QMat, i+1, position, value);
}



template <typename T>
void LaspackMatrix<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 LaspackMatrix<T>::add (const T a_in, const SparseMatrix<T> & X_in)
{
  libmesh_assert (this->initialized());
  libmesh_assert_equal_to (this->m(), X_in.m());
  libmesh_assert_equal_to (this->n(), X_in.n());

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

  _LPNumber a = static_cast<_LPNumber> (a_in);

  libmesh_assert(X);

  // loops taken from LaspackMatrix<T>::zero ()

  const numeric_index_type n_rows = this->m();

  for (numeric_index_type row=0; row<n_rows; row++)
    {
      auto r_start = _row_start[row];

      const numeric_index_type len = (_row_start[row+1] - _row_start[row]);

      // Make sure we agree on the row length
      libmesh_assert_equal_to (len, Q_GetLen(&_QMat, row+1));
      // compare matrix sparsity structures
      libmesh_assert_equal_to (len, Q_GetLen(&(X->_QMat), row+1));


      for (numeric_index_type l=0; l<len; l++)
        {
          const numeric_index_type j = *(r_start + l);

          // Make sure the data structures are working
          libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, row+1, l));

          const _LPNumber value = a * Q_GetEl(const_cast<QMatrix*>(&(X->_QMat)), row+1, j+1);
          Q_AddVal   (&_QMat, row+1, l, value);
        }
    }
}




template <typename T>
T LaspackMatrix<T>::operator () (const numeric_index_type i,
                                 const numeric_index_type j) const
{
  libmesh_assert (this->initialized());
  libmesh_assert_less (i, this->m());
  libmesh_assert_less (j, this->n());

  return Q_GetEl (const_cast<QMatrix*>(&_QMat), i+1, j+1);
}



template <typename T>
numeric_index_type LaspackMatrix<T>::pos (const numeric_index_type i,
                                          const numeric_index_type j) const
{
  libmesh_assert_less (i, this->m());
  libmesh_assert_less (j, this->n());
  libmesh_assert_less (i+1, _row_start.size());
  libmesh_assert (_row_start.back() == _csr.end());

  // note this requires the _csr to be sorted
  auto p = std::equal_range (_row_start[i], _row_start[i+1], j);

  // Make sure the row contains the element j
  libmesh_assert (p.first != p.second);

  // Make sure the values match
  libmesh_assert (*p.first == j);

  // Return the position in the compressed row
  return std::distance (_row_start[i], p.first);
}



template <typename T>
void LaspackMatrix<T>::close()
{
  libmesh_assert(this->initialized());

  this->_closed = true;

  // We've probably changed some entries so we need to tell LASPACK
  // that cached data is now invalid.
  *_QMat.DiagElAlloc = _LPFalse;
  *_QMat.ElSorted = _LPFalse;
  if (*_QMat.ILUExists)
    {
      *_QMat.ILUExists = _LPFalse;
      Q_Destr(_QMat.ILU);
    }
}



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

} // namespace libMesh


#endif // #ifdef LIBMESH_HAVE_LASPACK