Line data    Source code

       1             : // The libMesh Finite Element Library.
       2             : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
       3             : 
       4             : // This library is free software; you can redistribute it and/or
       5             : // modify it under the terms of the GNU Lesser General Public
       6             : // License as published by the Free Software Foundation; either
       7             : // version 2.1 of the License, or (at your option) any later version.
       8             : 
       9             : // This library is distributed in the hope that it will be useful,
      10             : // but WITHOUT ANY WARRANTY; without even the implied warranty of
      11             : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      12             : // Lesser General Public License for more details.
      13             : 
      14             : // You should have received a copy of the GNU Lesser General Public
      15             : // License along with this library; if not, write to the Free Software
      16             : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
      17             : 
      18             : #include "libmesh/reduced_sparse_matrix.h"
      19             : 
      20             : #ifdef LIBMESH_HAVE_PETSC
      21             : 
      22             : #include "libmesh/reduced_dof_map.h"
      23             : #include "libmesh/petsc_matrix.h"
      24             : #include "libmesh/dense_matrix.h"
      25             : #include "libmesh/dof_object.h"
      26             : #include "libmesh/int_range.h"
      27             : #include "libmesh/numeric_vector.h"
      28             : 
      29             : namespace libMesh
      30             : {
      31             : 
      32             : template <typename T>
      33          36 : ReducedSparseMatrix<T>::ReducedSparseMatrix(const Parallel::Communicator & comm_in,
      34             :                                             const ReducedDofMap & reduced_dof_map)
      35          36 :   : PetscMatrixBase<T>(comm_in), _reduced_dof_map(reduced_dof_map)
      36             : {
      37             :   // The wrapped matrix owns the Mat
      38          36 :   this->_destroy_mat_on_exit = false;
      39          36 : }
      40             : 
      41             : template <typename T>
      42          36 : ReducedSparseMatrix<T>::~ReducedSparseMatrix() = default;
      43             : 
      44             : template <typename T>
      45          24 : void ReducedSparseMatrix<T>::adopt_reduced_matrix()
      46             : {
      47          36 :   this->_mat = _reduced->mat();
      48          36 :   this->_is_initialized = true;
      49             :   // Solvers recover this object from the Mat, so it, not the wrapped matrix, is the context
      50          36 :   this->set_context();
      51          24 : }
      52             : 
      53             : template <typename T>
      54          36 : void ReducedSparseMatrix<T>::init(ParallelType)
      55             : {
      56          12 :   libmesh_assert(_reduced_dof_map.initialized());
      57             : 
      58          36 :   const numeric_index_type n = _reduced_dof_map.n_dofs();
      59          12 :   const numeric_index_type n_local = _reduced_dof_map.n_local_dofs();
      60             : 
      61          60 :   _reduced = std::make_unique<PetscMatrix<T>>(this->comm());
      62          48 :   _reduced->init(
      63          36 :       n, n, n_local, n_local, _reduced_dof_map.reduced_n_nz(), _reduced_dof_map.reduced_n_oz());
      64          24 :   this->adopt_reduced_matrix();
      65          36 : }
      66             : 
      67             : template <typename T>
      68           0 : void ReducedSparseMatrix<T>::init(const numeric_index_type m,
      69             :                                   const numeric_index_type n,
      70             :                                   const numeric_index_type m_l,
      71             :                                   const numeric_index_type n_l,
      72             :                                   const numeric_index_type nnz,
      73             :                                   const numeric_index_type noz,
      74             :                                   const numeric_index_type blocksize)
      75             : {
      76           0 :   _reduced = std::make_unique<PetscMatrix<T>>(this->comm());
      77           0 :   _reduced->init(m, n, m_l, n_l, nnz, noz, blocksize);
      78           0 :   this->adopt_reduced_matrix();
      79           0 : }
      80             : 
      81             : template <typename T>
      82          36 : void ReducedSparseMatrix<T>::clear() noexcept
      83             : {
      84          12 :   _reduced.reset();
      85          36 :   this->_mat = nullptr;
      86          36 :   this->_is_initialized = false;
      87          36 : }
      88             : 
      89             : template <typename T>
      90           0 : PetscMatrix<T> & ReducedSparseMatrix<T>::reduced_matrix()
      91             : {
      92           0 :   libmesh_assert(_reduced);
      93           0 :   return *_reduced;
      94             : }
      95             : 
      96             : template <typename T>
      97           0 : SparseMatrix<T> & ReducedSparseMatrix<T>::operator=(const SparseMatrix<T> & v)
      98             : {
      99           0 :   const auto & other = cast_ref<const ReducedSparseMatrix<T> &>(v);
     100           0 :   libmesh_assert(&other._reduced_dof_map == &_reduced_dof_map);
     101           0 :   *_reduced = *other._reduced;
     102           0 :   return *this;
     103             : }
     104             : 
     105             : template <typename T>
     106        3990 : void ReducedSparseMatrix<T>::to_reduced(const std::vector<numeric_index_type> & full,
     107             :                                         std::vector<PetscInt> & reduced) const
     108             : {
     109        5320 :   reduced.resize(full.size());
     110       54402 :   for (const auto i : index_range(full))
     111             :     {
     112       67216 :       const dof_id_type r = _reduced_dof_map.find_reduced_dof(full[i]);
     113       50412 :       reduced[i] = (r == DofObject::invalid_id) ? PetscInt(-1) : static_cast<PetscInt>(r);
     114             :     }
     115        3990 : }
     116             : 
     117             : template <typename T>
     118        1995 : void ReducedSparseMatrix<T>::add_matrix(const DenseMatrix<T> & dm,
     119             :                                         const std::vector<numeric_index_type> & rows,
     120             :                                         const std::vector<numeric_index_type> & cols)
     121             : {
     122         665 :   libmesh_assert(this->initialized());
     123         665 :   libmesh_assert_equal_to(rows.size(), dm.m());
     124         665 :   libmesh_assert_equal_to(cols.size(), dm.n());
     125             : 
     126        1995 :   if (rows.empty() || cols.empty())
     127           0 :     return;
     128             : 
     129        1995 :   this->to_reduced(rows, _rows);
     130        1995 :   this->to_reduced(cols, _cols);
     131        3325 :   LibmeshPetscCall(MatSetValues(this->_mat,
     132             :                                 cast_int<PetscInt>(_rows.size()),
     133             :                                 _rows.data(),
     134             :                                 cast_int<PetscInt>(_cols.size()),
     135             :                                 _cols.data(),
     136             :                                 pPS(const_cast<T *>(dm.get_values().data())),
     137             :                                 ADD_VALUES));
     138             : }
     139             : 
     140             : template <typename T>
     141        1995 : void ReducedSparseMatrix<T>::add_matrix(const DenseMatrix<T> & dm,
     142             :                                         const std::vector<numeric_index_type> & dof_indices)
     143             : {
     144        1995 :   this->add_matrix(dm, dof_indices, dof_indices);
     145        1995 : }
     146             : 
     147             : template <typename T>
     148           0 : void ReducedSparseMatrix<T>::set(const numeric_index_type i,
     149             :                                  const numeric_index_type j,
     150             :                                  const T value)
     151             : {
     152           0 :   const dof_id_type ri = _reduced_dof_map.find_reduced_dof(i);
     153           0 :   const dof_id_type rj = _reduced_dof_map.find_reduced_dof(j);
     154           0 :   if (ri != DofObject::invalid_id && rj != DofObject::invalid_id)
     155           0 :     _reduced->set(ri, rj, value);
     156           0 : }
     157             : 
     158             : template <typename T>
     159           0 : void ReducedSparseMatrix<T>::add(const numeric_index_type i,
     160             :                                  const numeric_index_type j,
     161             :                                  const T value)
     162             : {
     163           0 :   const dof_id_type ri = _reduced_dof_map.find_reduced_dof(i);
     164           0 :   const dof_id_type rj = _reduced_dof_map.find_reduced_dof(j);
     165           0 :   if (ri != DofObject::invalid_id && rj != DofObject::invalid_id)
     166           0 :     _reduced->add(ri, rj, value);
     167           0 : }
     168             : 
     169             : template <typename T>
     170         240 : void ReducedSparseMatrix<T>::zero()
     171             : {
     172         240 :   _reduced->zero();
     173         240 : }
     174             : 
     175             : template <typename T>
     176         210 : void ReducedSparseMatrix<T>::close()
     177             : {
     178         210 :   _reduced->close();
     179         210 : }
     180             : 
     181             : template <typename T>
     182           0 : void ReducedSparseMatrix<T>::flush()
     183             : {
     184           0 :   _reduced->flush();
     185           0 : }
     186             : 
     187             : template <typename T>
     188           0 : std::unique_ptr<SparseMatrix<T>> ReducedSparseMatrix<T>::zero_clone() const
     189             : {
     190           0 :   auto ret = std::make_unique<ReducedSparseMatrix<T>>(this->comm(), _reduced_dof_map);
     191           0 :   ret->_reduced.reset(cast_ptr<PetscMatrix<T> *>(_reduced->zero_clone().release()));
     192           0 :   ret->adopt_reduced_matrix();
     193           0 :   return ret;
     194           0 : }
     195             : 
     196             : template <typename T>
     197           0 : std::unique_ptr<SparseMatrix<T>> ReducedSparseMatrix<T>::clone() const
     198             : {
     199           0 :   auto ret = std::make_unique<ReducedSparseMatrix<T>>(this->comm(), _reduced_dof_map);
     200           0 :   ret->_reduced.reset(cast_ptr<PetscMatrix<T> *>(_reduced->clone().release()));
     201           0 :   ret->adopt_reduced_matrix();
     202           0 :   return ret;
     203           0 : }
     204             : 
     205             : template <typename T>
     206           0 : void ReducedSparseMatrix<T>::add(const T a, const SparseMatrix<T> & X)
     207             : {
     208           0 :   const auto & other = cast_ref<const ReducedSparseMatrix<T> &>(X);
     209           0 :   _reduced->add(a, *other._reduced);
     210           0 : }
     211             : 
     212             : template <typename T>
     213        3108 : T ReducedSparseMatrix<T>::operator()(const numeric_index_type i, const numeric_index_type j) const
     214             : {
     215        3108 :   return (*_reduced)(i, j);
     216             : }
     217             : 
     218             : template <typename T>
     219           6 : Real ReducedSparseMatrix<T>::l1_norm() const
     220             : {
     221           6 :   return _reduced->l1_norm();
     222             : }
     223             : 
     224             : template <typename T>
     225           0 : Real ReducedSparseMatrix<T>::linfty_norm() const
     226             : {
     227           0 :   return _reduced->linfty_norm();
     228             : }
     229             : 
     230             : template <typename T>
     231           0 : void ReducedSparseMatrix<T>::print_personal(std::ostream & os) const
     232             : {
     233           0 :   _reduced->print_personal(os);
     234           0 : }
     235             : 
     236             : template <typename T>
     237           0 : void ReducedSparseMatrix<T>::get_diagonal(NumericVector<T> & dest) const
     238             : {
     239           0 :   _reduced->get_diagonal(dest);
     240           0 : }
     241             : 
     242             : template <typename T>
     243           0 : void ReducedSparseMatrix<T>::get_transpose(SparseMatrix<T> & dest) const
     244             : {
     245           0 :   auto & other = cast_ref<ReducedSparseMatrix<T> &>(dest);
     246           0 :   libmesh_assert(&other._reduced_dof_map == &_reduced_dof_map);
     247           0 :   if (!other._reduced)
     248           0 :     other._reduced = std::make_unique<PetscMatrix<T>>(this->comm());
     249             :   // The transpose replaces the destination's Mat, which the destination must then adopt
     250           0 :   _reduced->get_transpose(*other._reduced);
     251           0 :   other.adopt_reduced_matrix();
     252           0 : }
     253             : 
     254             : template <typename T>
     255           0 : void ReducedSparseMatrix<T>::get_row(numeric_index_type i,
     256             :                                      std::vector<numeric_index_type> & indices,
     257             :                                      std::vector<T> & values) const
     258             : {
     259           0 :   _reduced->get_row(i, indices, values);
     260           0 : }
     261             : 
     262             : template <typename T>
     263           0 : void ReducedSparseMatrix<T>::scale(const T scale)
     264             : {
     265           0 :   _reduced->scale(scale);
     266           0 : }
     267             : 
     268             : template class LIBMESH_EXPORT ReducedSparseMatrix<Number>;
     269             : 
     270             : } // namespace libMesh
     271             : 
     272             : #endif // LIBMESH_HAVE_PETSC