SDP preprocessing - facial reduction

external/Spectra/include/Spectra/GenEigsComplexShiftSolver.h

@@ -1,11 +1,11 @@
-// Copyright (C) 2016-2019 Yixuan Qiu <[email protected]>
+// Copyright (C) 2016-2025 Yixuan Qiu <[email protected]>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
 // with this file, You can obtain one at https://mozilla.org/MPL/2.0/.
 
-#ifndef GEN_EIGS_COMPLEX_SHIFT_SOLVER_H
-#define GEN_EIGS_COMPLEX_SHIFT_SOLVER_H
+#ifndef SPECTRA_GEN_EIGS_COMPLEX_SHIFT_SOLVER_H
+#define SPECTRA_GEN_EIGS_COMPLEX_SHIFT_SOLVER_H
 
 #include <Eigen/Core>
 
@@ -15,7 +15,6 @@
 
 namespace Spectra {
 
-
 ///
 /// \ingroup EigenSolver
 ///
@@ -24,33 +23,35 @@ namespace Spectra {
 /// knowledge of the shift-and-invert mode can be found in the documentation
 /// of the SymEigsShiftSolver class.
 ///
-/// \tparam Scalar        The element type of the matrix.
-///                       Currently supported types are `float`, `double` and `long double`.
-/// \tparam SelectionRule An enumeration value indicating the selection rule of
-///                       the shifted-and-inverted eigenvalues.
-///                       The full list of enumeration values can be found in
-///                       \ref Enumerations.
-/// \tparam OpType        The name of the matrix operation class. Users could either
-///                       use the DenseGenComplexShiftSolve wrapper class, or define their
-///                       own that implements all the public member functions as in
-///                       DenseGenComplexShiftSolve.
+/// \tparam OpType  The name of the matrix operation class. Users could either
+///                 use the wrapper classes such as DenseGenComplexShiftSolve and
+///                 SparseGenComplexShiftSolve, or define their own that implements the type
+///                 definition `Scalar` and all the public member functions as in
+///                 DenseGenComplexShiftSolve.
 ///
-template <typename Scalar = double,
-          int SelectionRule = LARGEST_MAGN,
-          typename OpType = DenseGenComplexShiftSolve<double> >
-class GenEigsComplexShiftSolver: public GenEigsBase<Scalar, SelectionRule, OpType, IdentityBOp>
+template <typename OpType = DenseGenComplexShiftSolve<double>>
+class GenEigsComplexShiftSolver : public GenEigsBase<OpType, IdentityBOp>
 {
 private:
-    typedef Eigen::Index Index;
-    typedef std::complex<Scalar> Complex;
-    typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vector;
-    typedef Eigen::Matrix<Complex, Eigen::Dynamic, 1> ComplexVector;
+    using Scalar = typename OpType::Scalar;
+    using Index = Eigen::Index;
+    using Complex = std::complex<Scalar>;
+    using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
+    using ComplexArray = Eigen::Array<Complex, Eigen::Dynamic, 1>;
+
+    using Base = GenEigsBase<OpType, IdentityBOp>;
+    using Base::m_op;
+    using Base::m_n;
+    using Base::m_nev;
+    using Base::m_fac;
+    using Base::m_ritz_val;
+    using Base::m_ritz_vec;
 
     const Scalar m_sigmar;
     const Scalar m_sigmai;
 
     // First transform back the Ritz values, and then sort
-    void sort_ritzpair(int sort_rule)
+    void sort_ritzpair(SortRule sort_rule) override
     {
         using std::abs;
         using std::sqrt;
@@ -78,59 +79,62 @@ class GenEigsComplexShiftSolver: public GenEigsBase<Scalar, SelectionRule, OpTyp
         SimpleRandom<Scalar> rng(0);
         const Scalar shiftr = rng.random() * m_sigmar + rng.random();
         const Complex shift = Complex(shiftr, Scalar(0));
-        this->m_op->set_shift(shiftr, Scalar(0));
+        m_op.set_shift(shiftr, Scalar(0));
 
         // Calculate inv(A - r * I) * vj
-        Vector v_real(this->m_n), v_imag(this->m_n), OPv_real(this->m_n), OPv_imag(this->m_n);
-        const Scalar eps = Eigen::NumTraits<Scalar>::epsilon();
-        for(Index i = 0; i < this->m_nev; i++)
+        Vector v_real(m_n), v_imag(m_n), OPv_real(m_n), OPv_imag(m_n);
+        const Scalar eps = TypeTraits<Scalar>::epsilon();
+        for (Index i = 0; i < m_nev; i++)
         {
-            v_real.noalias() = this->m_fac.matrix_V() * this->m_ritz_vec.col(i).real();
-            v_imag.noalias() = this->m_fac.matrix_V() * this->m_ritz_vec.col(i).imag();
-            this->m_op->perform_op(v_real.data(), OPv_real.data());
-            this->m_op->perform_op(v_imag.data(), OPv_imag.data());
+            v_real.noalias() = m_fac.matrix_V() * m_ritz_vec.col(i).real();
+            v_imag.noalias() = m_fac.matrix_V() * m_ritz_vec.col(i).imag();
+            m_op.perform_op(v_real.data(), OPv_real.data());
+            m_op.perform_op(v_imag.data(), OPv_imag.data());
 
             // Two roots computed from the quadratic equation
-            const Complex nu = this->m_ritz_val[i];
+            const Complex nu = m_ritz_val[i];
             const Complex root_part1 = m_sigmar + Scalar(0.5) / nu;
             const Complex root_part2 = Scalar(0.5) * sqrt(Scalar(1) - Scalar(4) * m_sigmai * m_sigmai * (nu * nu)) / nu;
             const Complex root1 = root_part1 + root_part2;
             const Complex root2 = root_part1 - root_part2;
 
             // Test roots
             Scalar err1 = Scalar(0), err2 = Scalar(0);
-            for(int k = 0; k < this->m_n; k++)
+            for (int k = 0; k < m_n; k++)
             {
                 const Complex rhs1 = Complex(v_real[k], v_imag[k]) / (root1 - shift);
                 const Complex rhs2 = Complex(v_real[k], v_imag[k]) / (root2 - shift);
-                const Complex OPv  = Complex(OPv_real[k], OPv_imag[k]);
+                const Complex OPv = Complex(OPv_real[k], OPv_imag[k]);
                 err1 += norm(OPv - rhs1);
                 err2 += norm(OPv - rhs2);
             }
 
             const Complex lambdaj = (err1 < err2) ? root1 : root2;
-            this->m_ritz_val[i] = lambdaj;
+            m_ritz_val[i] = lambdaj;
 
-            if(abs(Eigen::numext::imag(lambdaj)) > eps)
+            if (abs(Eigen::numext::imag(lambdaj)) > eps)
             {
-                this->m_ritz_val[i + 1] = Eigen::numext::conj(lambdaj);
+                m_ritz_val[i + 1] = Eigen::numext::conj(lambdaj);
                 i++;
-            } else {
-                this->m_ritz_val[i] = Complex(Eigen::numext::real(lambdaj), Scalar(0));
+            }
+            else
+            {
+                m_ritz_val[i] = Complex(Eigen::numext::real(lambdaj), Scalar(0));
             }
         }
 
-        GenEigsBase<Scalar, SelectionRule, OpType, IdentityBOp>::sort_ritzpair(sort_rule);
+        Base::sort_ritzpair(sort_rule);
     }
+
 public:
     ///
     /// Constructor to create a eigen solver object using the shift-and-invert mode.
     ///
-    /// \param op      Pointer to the matrix operation object. This class should implement
+    /// \param op      The matrix operation object that implements
     ///                the complex shift-solve operation of \f$A\f$: calculating
     ///                \f$\mathrm{Re}\{(A-\sigma I)^{-1}v\}\f$ for any vector \f$v\f$. Users could either
-    ///                create the object from the DenseGenComplexShiftSolve wrapper class, or
-    ///                define their own that implements all the public member functions
+    ///                create the object from the wrapper class such as DenseGenComplexShiftSolve, or
+    ///                define their own that implements all the public members
     ///                as in DenseGenComplexShiftSolve.
     /// \param nev     Number of eigenvalues requested. This should satisfy \f$1\le nev \le n-2\f$,
     ///                where \f$n\f$ is the size of matrix.
@@ -142,15 +146,14 @@ class GenEigsComplexShiftSolver: public GenEigsBase<Scalar, SelectionRule, OpTyp
     /// \param sigmar  The real part of the shift.
     /// \param sigmai  The imaginary part of the shift.
     ///
-    GenEigsComplexShiftSolver(OpType* op, Index nev, Index ncv, const Scalar& sigmar, const Scalar& sigmai) :
-        GenEigsBase<Scalar, SelectionRule, OpType, IdentityBOp>(op, NULL, nev, ncv),
+    GenEigsComplexShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigmar, const Scalar& sigmai) :
+        Base(op, IdentityBOp(), nev, ncv),
         m_sigmar(sigmar), m_sigmai(sigmai)
     {
-        this->m_op->set_shift(m_sigmar, m_sigmai);
+        op.set_shift(m_sigmar, m_sigmai);
     }
 };
 
+}  // namespace Spectra
 
-} // namespace Spectra
-
-#endif // GEN_EIGS_COMPLEX_SHIFT_SOLVER_H
+#endif  // SPECTRA_GEN_EIGS_COMPLEX_SHIFT_SOLVER_H

external/Spectra/include/Spectra/GenEigsRealShiftSolver.h

@@ -1,11 +1,11 @@
-// Copyright (C) 2016-2019 Yixuan Qiu <[email protected]>
+// Copyright (C) 2016-2025 Yixuan Qiu <[email protected]>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
 // with this file, You can obtain one at https://mozilla.org/MPL/2.0/.
 
-#ifndef GEN_EIGS_REAL_SHIFT_SOLVER_H
-#define GEN_EIGS_REAL_SHIFT_SOLVER_H
+#ifndef SPECTRA_GEN_EIGS_REAL_SHIFT_SOLVER_H
+#define SPECTRA_GEN_EIGS_REAL_SHIFT_SOLVER_H
 
 #include <Eigen/Core>
 
@@ -15,7 +15,6 @@
 
 namespace Spectra {
 
-
 ///
 /// \ingroup EigenSolver
 ///
@@ -24,48 +23,45 @@ namespace Spectra {
 /// knowledge of the shift-and-invert mode can be found in the documentation
 /// of the SymEigsShiftSolver class.
 ///
-/// \tparam Scalar        The element type of the matrix.
-///                       Currently supported types are `float`, `double` and `long double`.
-/// \tparam SelectionRule An enumeration value indicating the selection rule of
-///                       the shifted-and-inverted eigenvalues.
-///                       The full list of enumeration values can be found in
-///                       \ref Enumerations.
-/// \tparam OpType        The name of the matrix operation class. Users could either
-///                       use the wrapper classes such as DenseGenRealShiftSolve and
-///                       SparseGenRealShiftSolve, or define their
-///                       own that implements all the public member functions as in
-///                       DenseGenRealShiftSolve.
+/// \tparam OpType  The name of the matrix operation class. Users could either
+///                 use the wrapper classes such as DenseGenRealShiftSolve and
+///                 SparseGenRealShiftSolve, or define their own that implements the type
+///                 definition `Scalar` and all the public member functions as in
+///                 DenseGenRealShiftSolve.
 ///
-template <typename Scalar = double,
-          int SelectionRule = LARGEST_MAGN,
-          typename OpType = DenseGenRealShiftSolve<double> >
-class GenEigsRealShiftSolver: public GenEigsBase<Scalar, SelectionRule, OpType, IdentityBOp>
+template <typename OpType = DenseGenRealShiftSolve<double>>
+class GenEigsRealShiftSolver : public GenEigsBase<OpType, IdentityBOp>
 {
 private:
-    typedef Eigen::Index Index;
-    typedef std::complex<Scalar> Complex;
-    typedef Eigen::Array<Complex, Eigen::Dynamic, 1> ComplexArray;
+    using Scalar = typename OpType::Scalar;
+    using Index = Eigen::Index;
+    using Complex = std::complex<Scalar>;
+    using ComplexArray = Eigen::Array<Complex, Eigen::Dynamic, 1>;
+
+    using Base = GenEigsBase<OpType, IdentityBOp>;
+    using Base::m_nev;
+    using Base::m_ritz_val;
 
     const Scalar m_sigma;
 
     // First transform back the Ritz values, and then sort
-    void sort_ritzpair(int sort_rule)
+    void sort_ritzpair(SortRule sort_rule) override
     {
         // The eigenvalues we get from the iteration is nu = 1 / (lambda - sigma)
         // So the eigenvalues of the original problem is lambda = 1 / nu + sigma
-        ComplexArray ritz_val_org = Scalar(1.0) / this->m_ritz_val.head(this->m_nev).array() + m_sigma;
-        this->m_ritz_val.head(this->m_nev) = ritz_val_org;
-        GenEigsBase<Scalar, SelectionRule, OpType, IdentityBOp>::sort_ritzpair(sort_rule);
+        m_ritz_val.head(m_nev) = Scalar(1) / m_ritz_val.head(m_nev).array() + m_sigma;
+        Base::sort_ritzpair(sort_rule);
     }
+
 public:
     ///
     /// Constructor to create a eigen solver object using the shift-and-invert mode.
     ///
-    /// \param op     Pointer to the matrix operation object. This class should implement
+    /// \param op     The matrix operation object that implements
     ///               the shift-solve operation of \f$A\f$: calculating
     ///               \f$(A-\sigma I)^{-1}v\f$ for any vector \f$v\f$. Users could either
     ///               create the object from the wrapper class such as DenseGenRealShiftSolve, or
-    ///               define their own that implements all the public member functions
+    ///               define their own that implements all the public members
     ///               as in DenseGenRealShiftSolve.
     /// \param nev    Number of eigenvalues requested. This should satisfy \f$1\le nev \le n-2\f$,
     ///               where \f$n\f$ is the size of matrix.
@@ -76,15 +72,14 @@ class GenEigsRealShiftSolver: public GenEigsBase<Scalar, SelectionRule, OpType,
     ///               and is advised to take \f$ncv \ge 2\cdot nev + 1\f$.
     /// \param sigma  The real-valued shift.
     ///
-    GenEigsRealShiftSolver(OpType* op, Index nev, Index ncv, Scalar sigma) :
-        GenEigsBase<Scalar, SelectionRule, OpType, IdentityBOp>(op, NULL, nev, ncv),
+    GenEigsRealShiftSolver(OpType& op, Index nev, Index ncv, const Scalar& sigma) :
+        Base(op, IdentityBOp(), nev, ncv),
         m_sigma(sigma)
     {
-        this->m_op->set_shift(m_sigma);
+        op.set_shift(m_sigma);
     }
 };
 
+}  // namespace Spectra
 
-} // namespace Spectra
-
-#endif // GEN_EIGS_REAL_SHIFT_SOLVER_H
+#endif  // SPECTRA_GEN_EIGS_REAL_SHIFT_SOLVER_H