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             : 
      19             : 
      20             : // C++ includes
      21             : #include <algorithm> // for std::min
      22             : #include <limits>
      23             : 
      24             : // Local Includes
      25             : #include "libmesh/dense_subvector.h"
      26             : #include "libmesh/dense_vector.h"
      27             : #include "libmesh/eigen_sparse_vector.h"
      28             : #include "libmesh/eigen_sparse_matrix.h"
      29             : #include "libmesh/int_range.h"
      30             : 
      31             : #ifdef LIBMESH_HAVE_EIGEN
      32             : 
      33             : namespace libMesh
      34             : {
      35             : 
      36             : template <typename T>
      37         144 : T EigenSparseVector<T>::sum () const
      38             : {
      39           4 :   libmesh_assert (this->closed());
      40           4 :   libmesh_assert (this->initialized());
      41             : 
      42         144 :   return _vec.sum();
      43             : }
      44             : 
      45             : 
      46             : 
      47             : template <typename T>
      48         716 : Real EigenSparseVector<T>::l1_norm () const
      49             : {
      50          20 :   libmesh_assert (this->closed());
      51          20 :   libmesh_assert (this->initialized());
      52             : 
      53         716 :   return _vec.lpNorm<1>();
      54             : }
      55             : 
      56             : 
      57             : 
      58             : template <typename T>
      59        2950 : Real EigenSparseVector<T>::l2_norm () const
      60             : {
      61           8 :   libmesh_assert (this->closed());
      62           8 :   libmesh_assert (this->initialized());
      63             : 
      64        2950 :   return _vec.lpNorm<2>();
      65             : }
      66             : 
      67             : 
      68             : 
      69             : template <typename T>
      70         143 : Real EigenSparseVector<T>::linfty_norm () const
      71             : {
      72           4 :   libmesh_assert (this->closed());
      73           4 :   libmesh_assert (this->initialized());
      74             : 
      75         143 :   return _vec.lpNorm<Eigen::Infinity>();
      76             : }
      77             : 
      78             : 
      79             : 
      80             : template <typename T>
      81         145 : NumericVector<T> & EigenSparseVector<T>::operator += (const NumericVector<T> & v_in)
      82             : {
      83           4 :   libmesh_assert (this->closed());
      84             : 
      85           4 :   const EigenSparseVector<T> & v = cast_ref<const EigenSparseVector<T> &>(v_in);
      86             : 
      87         145 :   _vec += v._vec;
      88             : 
      89         145 :   return *this;
      90             : }
      91             : 
      92             : 
      93             : 
      94             : 
      95             : template <typename T>
      96         356 : NumericVector<T> & EigenSparseVector<T>::operator -= (const NumericVector<T> & v_in)
      97             : {
      98           8 :   libmesh_assert (this->closed());
      99             : 
     100           8 :   const EigenSparseVector<T> & v = cast_ref<const EigenSparseVector<T> &>(v_in);
     101             : 
     102         356 :   _vec -= v._vec;
     103             : 
     104         356 :   return *this;
     105             : }
     106             : 
     107             : 
     108             : 
     109             : template <typename T>
     110         142 : NumericVector<T> & EigenSparseVector<T>::operator *= (const NumericVector<T> & v_in)
     111             : {
     112           4 :   libmesh_assert (this->closed());
     113           4 :   libmesh_assert_equal_to(size(), v_in.size());
     114             : 
     115           4 :   const EigenSparseVector<T> & v = cast_ref<const EigenSparseVector<T> &>(v_in);
     116             : 
     117         142 :   _vec = _vec.cwiseProduct(v._vec);
     118             : 
     119         142 :   return *this;
     120             : }
     121             : 
     122             : 
     123             : 
     124             : template <typename T>
     125         515 : NumericVector<T> & EigenSparseVector<T>::operator /= (const NumericVector<T> & v_in)
     126             : {
     127           4 :   libmesh_assert (this->closed());
     128           4 :   libmesh_assert_equal_to(size(), v_in.size());
     129             : 
     130           4 :   const EigenSparseVector<T> & v = cast_ref<const EigenSparseVector<T> &>(v_in);
     131             : 
     132         515 :   _vec = _vec.cwiseQuotient(v._vec);
     133             : 
     134         515 :   return *this;
     135             : }
     136             : 
     137             : 
     138             : 
     139             : 
     140             : template <typename T>
     141         142 : void EigenSparseVector<T>::reciprocal()
     142             : {
     143             : #ifndef NDEBUG
     144           4 :   const numeric_index_type n = this->size();
     145             : 
     146          44 :   for (numeric_index_type i=0; i<n; i++)
     147             :     // Don't divide by zero!
     148          40 :     libmesh_assert_not_equal_to ((*this)(i), T(0));
     149             : #endif
     150             : 
     151         142 :   _vec = _vec.cwiseInverse();
     152         142 : }
     153             : 
     154             : 
     155             : 
     156             : template <typename T>
     157           0 : void EigenSparseVector<T>::conjugate()
     158             : {
     159           0 :   _vec = _vec.conjugate();
     160           0 : }
     161             : 
     162             : 
     163             : 
     164             : template <typename T>
     165         143 : void EigenSparseVector<T>::add (const T v)
     166             : {
     167         143 :   _vec += EigenSV::Constant(this->size(), v);
     168         143 : }
     169             : 
     170             : 
     171             : 
     172             : 
     173             : template <typename T>
     174         596 : void EigenSparseVector<T>::add (const NumericVector<T> & v_in)
     175             : {
     176           0 :   libmesh_assert (this->initialized());
     177             : 
     178           0 :   const EigenSparseVector<T> & v = cast_ref<const EigenSparseVector<T> &>(v_in);
     179             : 
     180         596 :   _vec += v._vec;
     181         596 : }
     182             : 
     183             : 
     184             : 
     185             : template <typename T>
     186       26575 : void EigenSparseVector<T>::add (const T a, const NumericVector<T> & v_in)
     187             : {
     188           4 :   libmesh_assert (this->initialized());
     189             : 
     190           4 :   const EigenSparseVector<T> & v = cast_ref<const EigenSparseVector<T> &>(v_in);
     191             : 
     192       26579 :   _vec += v._vec*a;
     193       26575 : }
     194             : 
     195             : 
     196             : 
     197             : template <typename T>
     198       22913 : void EigenSparseVector<T>::add_vector (const NumericVector<T> & vec_in,
     199             :                                        const SparseMatrix<T>  & mat_in)
     200             : {
     201             :   // Make sure the data passed in are really in Eigen types
     202           0 :   const EigenSparseVector<T> * e_vec = cast_ptr<const EigenSparseVector<T> *>(&vec_in);
     203           0 :   const EigenSparseMatrix<T> * mat = cast_ptr<const EigenSparseMatrix<T> *>(&mat_in);
     204             : 
     205           0 :   libmesh_assert(e_vec);
     206           0 :   libmesh_assert(mat);
     207             : 
     208       22913 :   _vec += mat->_mat*e_vec->_vec;
     209       22913 : }
     210             : 
     211             : 
     212             : 
     213             : template <typename T>
     214           0 : void EigenSparseVector<T>::add_vector_transpose (const NumericVector<T> & vec_in,
     215             :                                                  const SparseMatrix<T>  & mat_in)
     216             : {
     217             :   // Make sure the data passed in are really in Eigen types
     218           0 :   const EigenSparseVector<T> * e_vec = cast_ptr<const EigenSparseVector<T> *>(&vec_in);
     219           0 :   const EigenSparseMatrix<T> * mat = cast_ptr<const EigenSparseMatrix<T> *>(&mat_in);
     220             : 
     221           0 :   libmesh_assert(e_vec);
     222           0 :   libmesh_assert(mat);
     223             : 
     224           0 :   _vec += mat->_mat.transpose()*e_vec->_vec;
     225           0 : }
     226             : 
     227             : 
     228             : 
     229             : template <typename T>
     230        1003 : void EigenSparseVector<T>::scale (const T factor)
     231             : {
     232           4 :   libmesh_assert (this->initialized());
     233             : 
     234        1003 :   _vec *= factor;
     235        1003 : }
     236             : 
     237             : 
     238             : 
     239             : template <typename T>
     240           0 : void EigenSparseVector<T>::abs()
     241             : {
     242           0 :   libmesh_assert (this->initialized());
     243             : 
     244           0 :   const numeric_index_type n = this->size();
     245             : 
     246           0 :   for (numeric_index_type i=0; i!=n; ++i)
     247           0 :     this->set(i,std::abs((*this)(i)));
     248           0 : }
     249             : 
     250             : 
     251             : 
     252             : template <typename T>
     253        9102 : T EigenSparseVector<T>::dot (const NumericVector<T> & v_in) const
     254             : {
     255           4 :   libmesh_assert (this->initialized());
     256             : 
     257             :   // Make sure the NumericVector passed in is really a EigenSparseVector
     258           4 :   const EigenSparseVector<T> * v = cast_ptr<const EigenSparseVector<T> *>(&v_in);
     259           4 :   libmesh_assert(v);
     260             : 
     261        9106 :   return _vec.dot(v->_vec);
     262             : }
     263             : 
     264             : 
     265             : 
     266             : template <typename T>
     267             : NumericVector<T> &
     268         142 : EigenSparseVector<T>::operator = (const T s)
     269             : {
     270           4 :   libmesh_assert (this->initialized());
     271           4 :   libmesh_assert (this->closed());
     272             : 
     273         138 :   _vec.fill(s);
     274             : 
     275         142 :   return *this;
     276             : }
     277             : 
     278             : 
     279             : 
     280             : template <typename T>
     281             : NumericVector<T> &
     282        7606 : EigenSparseVector<T>::operator = (const NumericVector<T> & v_in)
     283             : {
     284             :   // Make sure the NumericVector passed in is really a EigenSparseVector
     285          10 :   const EigenSparseVector<T> * v =
     286          10 :     cast_ptr<const EigenSparseVector<T> *>(&v_in);
     287             : 
     288          10 :   libmesh_assert(v);
     289             : 
     290          20 :   *this = *v;
     291             : 
     292        7606 :   return *this;
     293             : }
     294             : 
     295             : 
     296             : 
     297             : template <typename T>
     298             : EigenSparseVector<T> &
     299          91 : EigenSparseVector<T>::operator = (const EigenSparseVector<T> & v)
     300             : {
     301          12 :   libmesh_assert (this->initialized());
     302          12 :   libmesh_assert (v.closed());
     303          12 :   libmesh_assert_equal_to (this->size(), v.size());
     304             : 
     305        7677 :   _vec = v._vec;
     306             : 
     307       20304 :   this->_is_closed = true;
     308             : 
     309          91 :   return *this;
     310             : }
     311             : 
     312             : 
     313             : 
     314             : template <typename T>
     315             : NumericVector<T> &
     316           1 : EigenSparseVector<T>::operator = (const std::vector<T> & v)
     317             : {
     318             :   /**
     319             :    * Case 1:  The vector is the same size of
     320             :    * The global vector.  Only add the local components.
     321             :    */
     322           1 :   if (this->size() == v.size())
     323        4298 :     for (auto i : index_range(v))
     324        4297 :       this->set (i, v[i]);
     325             : 
     326             :   else
     327           0 :     libmesh_error_msg("this->size() = " << this->size() << " must be equal to v.size() = " << v.size());
     328             : 
     329           1 :   return *this;
     330             : }
     331             : 
     332             : 
     333             : template <typename T>
     334        1329 : void EigenSparseVector<T>::localize (NumericVector<T> & v_local_in) const
     335             : {
     336             :   // Make sure the NumericVector passed in is really a EigenSparseVector
     337           0 :   EigenSparseVector<T> * v_local =
     338           0 :     cast_ptr<EigenSparseVector<T> *>(&v_local_in);
     339             : 
     340           0 :   libmesh_assert(v_local);
     341             : 
     342           0 :   *v_local = *this;
     343        1329 : }
     344             : 
     345             : 
     346             : 
     347             : template <typename T>
     348       11298 : void EigenSparseVector<T>::localize (NumericVector<T> & v_local_in,
     349             :                                      const std::vector<numeric_index_type> & libmesh_dbg_var(send_list)) const
     350             : {
     351             :   // Make sure the NumericVector passed in is really a EigenSparseVector
     352           0 :   EigenSparseVector<T> * v_local =
     353           0 :     cast_ptr<EigenSparseVector<T> *>(&v_local_in);
     354             : 
     355           0 :   libmesh_assert(v_local);
     356           0 :   libmesh_assert_less_equal (send_list.size(), v_local->size());
     357             : 
     358           0 :   *v_local = *this;
     359       11298 : }
     360             : 
     361             : 
     362             : 
     363             : template <typename T>
     364         142 : void EigenSparseVector<T>::localize (std::vector<T> & v_local,
     365             :                                      const std::vector<numeric_index_type> & indices) const
     366             : {
     367             :   // EigenSparseVectors are serial, so we can just copy values
     368         146 :   v_local.resize(indices.size());
     369             : 
     370        1562 :   for (auto i : index_range(v_local))
     371        1500 :     v_local[i] = (*this)(indices[i]);
     372         142 : }
     373             : 
     374             : 
     375             : 
     376             : template <typename T>
     377         246 : void EigenSparseVector<T>::localize (const numeric_index_type libmesh_dbg_var(first_local_idx),
     378             :                                      const numeric_index_type libmesh_dbg_var(last_local_idx),
     379             :                                      const std::vector<numeric_index_type> & libmesh_dbg_var(send_list))
     380             : {
     381           0 :   libmesh_assert_equal_to (first_local_idx, 0);
     382           0 :   libmesh_assert_equal_to (last_local_idx+1, this->size());
     383             : 
     384           0 :   libmesh_assert_less_equal (send_list.size(), this->size());
     385             : 
     386         246 :   this->_is_closed = true;
     387         246 : }
     388             : 
     389             : 
     390             : 
     391             : template <typename T>
     392        1310 : void EigenSparseVector<T>::localize (std::vector<T> & v_local) const
     393             : 
     394             : {
     395        1310 :   v_local.resize(this->size());
     396             : 
     397     3039053 :   for (auto i : index_range(v_local))
     398     3037743 :     v_local[i] = (*this)(i);
     399        1310 : }
     400             : 
     401             : 
     402             : 
     403             : template <typename T>
     404         725 : void EigenSparseVector<T>::localize_to_one (std::vector<T> & v_local,
     405             :                                             const processor_id_type libmesh_dbg_var(pid)) const
     406             : {
     407           8 :   libmesh_assert_equal_to (pid, 0);
     408             : 
     409         725 :   this->localize (v_local);
     410         725 : }
     411             : 
     412             : 
     413             : 
     414             : template <typename T>
     415           0 : void EigenSparseVector<T>::pointwise_mult (const NumericVector<T> & /*vec1*/,
     416             :                                            const NumericVector<T> & /*vec2*/)
     417             : {
     418           0 :   libmesh_not_implemented();
     419             : }
     420             : 
     421             : template <typename T>
     422           0 : void EigenSparseVector<T>::pointwise_divide (const NumericVector<T> & /*vec1*/,
     423             :                                              const NumericVector<T> & /*vec2*/)
     424             : {
     425           0 :   libmesh_not_implemented();
     426             : }
     427             : 
     428             : 
     429             : 
     430             : template <typename T>
     431         426 : void EigenSparseVector<T>::create_subvector(NumericVector<T> & subvector,
     432             :                                             const std::vector<numeric_index_type> & rows,
     433             :                                             const bool /* supplying_global_rows */) const
     434             : {
     435             :   // Make sure the passed in subvector is really an EigenSparseVector
     436          12 :   EigenSparseVector<T> * eigen_subvector = cast_ptr<EigenSparseVector<T> *>(&subvector);
     437             : 
     438             :   // If the eigen_subvector is already initialized, we assume that the
     439             :   // user has already allocated the *correct* amount of space for it.
     440             :   // If not, we do so.
     441         426 :   if (!eigen_subvector->initialized())
     442           0 :     eigen_subvector->init(rows.size());
     443             : 
     444         426 :   eigen_subvector->vec() = this->vec()(rows);
     445             : 
     446         426 :   eigen_subvector->_is_closed = true;
     447         426 : }
     448             : 
     449             : 
     450             : 
     451             : template <typename T>
     452             : std::unique_ptr<NumericVector<T>>
     453         142 : EigenSparseVector<T>::get_subvector(const std::vector<numeric_index_type> & rows)
     454             : {
     455         150 :   auto returnval = std::make_unique<EigenSparseVector<T>>(this->comm(), rows.size());
     456             : 
     457             :   // We're just going to make a copy here, since we'll be calling a
     458             :   // restore later anyway.  Someone with more familiarity with Eigen
     459             :   // can see if there's a way to improve this later.
     460         142 :   this->create_subvector(*returnval, rows);
     461             : 
     462         142 :   this->_is_closed = false;
     463             : 
     464         146 :   return returnval;
     465             : }
     466             : 
     467             : template <typename T>
     468             : void
     469         142 : EigenSparseVector<T>::restore_subvector(std::unique_ptr<NumericVector<T>> subvector,
     470             :                                         const std::vector<numeric_index_type> & rows)
     471             : {
     472           4 :   auto * const eigen_subvector = cast_ptr<EigenSparseVector<T> *>(subvector.get());
     473             : 
     474         142 :   this->vec()(rows) = eigen_subvector->vec();
     475         142 :   this->_is_closed = true;
     476         142 : }
     477             : 
     478             : 
     479             : 
     480             : template <typename T>
     481           0 : Real EigenSparseVector<T>::max() const
     482             : {
     483           0 :   libmesh_assert (this->initialized());
     484           0 :   if (!this->size())
     485           0 :     return -std::numeric_limits<Real>::max();
     486             : 
     487             : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
     488           0 :   Real the_max = libmesh_real((*this)(0));
     489             : 
     490             :   const numeric_index_type n = this->size();
     491             : 
     492           0 :   for (numeric_index_type i=1; i<n; i++)
     493           0 :     the_max = std::max (the_max, libmesh_real((*this)(i)));
     494             : 
     495           0 :   return the_max;
     496             : #else
     497           0 :   return libmesh_real(_vec.maxCoeff());
     498             : #endif
     499             : }
     500             : 
     501             : 
     502             : 
     503             : template <typename T>
     504           0 : Real EigenSparseVector<T>::min () const
     505             : {
     506           0 :   libmesh_assert (this->initialized());
     507           0 :   if (!this->size())
     508           0 :     return std::numeric_limits<Real>::max();
     509             : 
     510             : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
     511           0 :   Real the_min = libmesh_real((*this)(0));
     512             : 
     513             :   const numeric_index_type n = this->size();
     514             : 
     515           0 :   for (numeric_index_type i=1; i<n; i++)
     516           0 :     the_min = std::min (the_min, libmesh_real((*this)(i)));
     517             : 
     518           0 :   return the_min;
     519             : #else
     520           0 :   return libmesh_real(_vec.minCoeff());
     521             : #endif
     522             : }
     523             : 
     524             : 
     525             : //------------------------------------------------------------------
     526             : // Explicit instantiations
     527             : template class LIBMESH_EXPORT EigenSparseVector<Number>;
     528             : 
     529             : } // namespace libMesh
     530             : 
     531             : 
     532             : #endif // #ifdef LIBMESH_HAVE_EIGEN