Line data    Source code

       1              : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
       2              : // SPDX-License-Identifier: Apache-2.0
       3              : 
       4              : #ifndef PALACE_LINALG_CHEBYSHEV_SMOOTHER_HPP
       5              : #define PALACE_LINALG_CHEBYSHEV_SMOOTHER_HPP
       6              : 
       7              : #include "linalg/operator.hpp"
       8              : #include "linalg/solver.hpp"
       9              : #include "linalg/vector.hpp"
      10              : 
      11              : namespace palace
      12              : {
      13              : 
      14              : //
      15              : // Matrix-free diagonally-scaled Chebyshev smoothing. This is largely the same as
      16              : // mfem::OperatorChebyshevSmoother allows a nonzero initial guess and uses alternative
      17              : // methods to estimate the largest eigenvalue. We use a smoother based on Chebyshev
      18              : // polynomials of the 4th-kind as proposed in recent work.
      19              : // Reference: Phillips and Fischer, Optimal Chebyshev smoothers and one-sided V-cycles,
      20              : //            arXiv:2210.03179v1 (2022).
      21              : //
      22              : template <typename OperType>
      23              : class ChebyshevSmoother : public Solver<OperType>
      24              : {
      25              :   using VecType = typename Solver<OperType>::VecType;
      26              : 
      27              : private:
      28              :   // MPI communicator associated with the solver operator and vectors.
      29              :   MPI_Comm comm;
      30              : 
      31              :   // Number of smoother iterations and polynomial order.
      32              :   const int pc_it, order;
      33              : 
      34              :   // System matrix (not owned).
      35              :   const OperType *A;
      36              : 
      37              :   // Inverse diagonal scaling of the operator (real-valued for now).
      38              :   VecType dinv;
      39              : 
      40              :   // Maximum operator eigenvalue for Chebyshev polynomial smoothing.
      41              :   double lambda_max, sf_max;
      42              : 
      43              :   // Temporary vector for smoother application.
      44              :   mutable VecType d, r;
      45              : 
      46              : public:
      47              :   ChebyshevSmoother(MPI_Comm comm, int smooth_it, int poly_order, double sf_max);
      48              : 
      49              :   void SetOperator(const OperType &op) override;
      50              : 
      51            0 :   void Mult(const VecType &x, VecType &y) const override
      52              :   {
      53            0 :     if (r.Size() != y.Size())
      54              :     {
      55            0 :       r.SetSize(y.Size());
      56            0 :       r.UseDevice(true);
      57              :     }
      58            0 :     Mult2(x, y, r);
      59            0 :   }
      60              : 
      61            0 :   void MultTranspose(const VecType &x, VecType &y) const override
      62              :   {
      63            0 :     if (r.Size() != y.Size())
      64              :     {
      65            0 :       r.SetSize(y.Size());
      66            0 :       r.UseDevice(true);
      67              :     }
      68            0 :     MultTranspose2(x, y, r);
      69            0 :   }
      70              : 
      71              :   void Mult2(const VecType &x, VecType &y, VecType &r) const override;
      72              : 
      73            0 :   void MultTranspose2(const VecType &x, VecType &y, VecType &r) const override
      74              :   {
      75            0 :     Mult2(x, y, r);  // Assumes operator symmetry
      76            0 :   }
      77              : };
      78              : 
      79              : //
      80              : // Matrix-free diagonally-scaled Chebyshev smoothing using standard 1st-kind Chebyshev
      81              : // polynomials.
      82              : // Reference: Adams et al., Parallel multigrid smoothing: polynomial versus Gauss–Seidel,
      83              : //            JCP (2003).
      84              : //
      85              : template <typename OperType>
      86              : class ChebyshevSmoother1stKind : public Solver<OperType>
      87              : {
      88              :   using VecType = typename Solver<OperType>::VecType;
      89              : 
      90              : private:
      91              :   // MPI communicator associated with the solver operator and vectors.
      92              :   MPI_Comm comm;
      93              : 
      94              :   // Number of smoother iterations and polynomial order.
      95              :   const int pc_it, order;
      96              : 
      97              :   // System matrix (not owned).
      98              :   const OperType *A;
      99              : 
     100              :   // Inverse diagonal scaling of the operator (real-valued for now).
     101              :   VecType dinv;
     102              : 
     103              :   // Parameters depending on maximum and minimum operator eigenvalue estimates for Chebyshev
     104              :   // polynomial smoothing.
     105              :   double theta, delta, sf_max, sf_min;
     106              : 
     107              :   // Temporary vector for smoother application.
     108              :   mutable VecType d, r;
     109              : 
     110              : public:
     111              :   ChebyshevSmoother1stKind(MPI_Comm comm, int smooth_it, int poly_order, double sf_max,
     112              :                            double sf_min);
     113              : 
     114              :   void SetOperator(const OperType &op) override;
     115              : 
     116            0 :   void Mult(const VecType &x, VecType &y) const override
     117              :   {
     118            0 :     if (r.Size() != y.Size())
     119              :     {
     120            0 :       r.SetSize(y.Size());
     121            0 :       r.UseDevice(true);
     122              :     }
     123            0 :     Mult2(x, y, r);
     124            0 :   }
     125              : 
     126            0 :   void MultTranspose(const VecType &x, VecType &y) const override
     127              :   {
     128            0 :     if (r.Size() != y.Size())
     129              :     {
     130            0 :       r.SetSize(y.Size());
     131            0 :       r.UseDevice(true);
     132              :     }
     133            0 :     MultTranspose2(x, y, r);
     134            0 :   }
     135              : 
     136              :   void Mult2(const VecType &x, VecType &y, VecType &r) const override;
     137              : 
     138            0 :   void MultTranspose2(const VecType &x, VecType &y, VecType &r) const override
     139              :   {
     140            0 :     Mult2(x, y, r);  // Assumes operator symmetry
     141            0 :   }
     142              : };
     143              : 
     144              : }  // namespace palace
     145              : 
     146              : #endif  // PALACE_LINALG_CHEBYSHEV_SMOOTHER_HPP