cl_FEM_DofMgr_EigenValues.hpp

Go to the documentation of this file.

/*
 * BELFEM -- The Berkeley Lab Finite Element Framework
 * Copyright (c) 2026, The Regents of the University of California,
 * through Lawrence Berkeley National Laboratory (subject to receipt of any required
 * approvals from the U.S. Dept. of Energy).  All rights reserved.
 *
 * Developers: Christian Messe, Gregory Giard
 *
 * See the top-level LICENSE file for the complete license and disclaimer.
 */
#ifndef CL_FEM_DOFMGR_EIGENVALUES_HPP
#define CL_FEM_DOFMGR_EIGENVALUES_HPP
 
#include "typedefs.hpp"
#include "cl_Vector.hpp"
#include "cl_Bitset.hpp"
#include "cl_Solver.hpp"
#include "arpacktools.hpp"
 
namespace belfem
{
    class SpMatrix ;
 
#ifdef BELFEM_PARPACK
    namespace sparse
    {
        template< typename T > class DistMatrixCSR ;
    }
#endif

    enum class EigenOutcome
    {
        Ok                  = 0,   
        NotConverged        = 1,   
        NotPositiveDefinite = 2,   
        Complex             = 3,   
        Inconsistent        = 4,   
        Infeasible          = 5,   
        SolverFailed        = 6    
    };
 
    enum class ArpackWhat
    {
        LambaMinMagn = 0,
        LambaMaxMagn = 1,
        LambaMinReal = 2,
        LambaMaxReal = 3,
        LambaMinImag = 4,
        LambaMaxImag = 5
    };
 
    namespace fem
    {
        class DofManager ;
 
        namespace dofmgr
        {
            class EigenValues
            {
                const proc_t mCommRank ;
                const proc_t mCommSize ;
                DofManager * mParent ;
 
                // matrix data
                SpMatrix   * mK = nullptr ;
                SpMatrix   * mM = nullptr ;
                int_t  mN   = 0 ;
                int_t  mNNZ = 0 ;
 
 
                // options
                int_t          mNumMinVals   = 1 ;
                int_t          mNumMaxVals   = 1 ;
                int_t          mNumMaxIter   = 300 ;
 
                // the three tuning knobs used to be hardwired here and were
                // reachable from nowhere -- no caller in the tree ever invoked
                // the setters. They are sized automatically now ( configure() ),
                // and these flags are what keeps an explicit caller in charge:
                // once a setter has been used, the heuristic leaves that knob
                // alone for the rest of the run
                bool           mSubspaceSizeExplicit = false ;
                bool           mToleranceExplicit    = false ;
                bool           mNumMaxIterExplicit   = false ;
 
                // ARPACK accepts a Ritz value when its error bound falls below
                // mEpsilon * |lambda|, and that bound cannot go below the
                // backward error of the matvec, roughly eps_mach * |lambda_max|.
                // The small end is therefore only reachable while
                //
                //     mEpsilon  >  eps_mach * kappa      ( ~ 1e-16 * kappa )
                //
                // so 1e-7 caps the usable condition number at about 1e9 -- below
                // what an h-phi Jacobian reaches on a refined mesh, where SM then
                // grinds to mNumMaxIter and converges nothing. 1e-4 lifts that
                // ceiling to ~1e12 and still pins lambda to four digits, which is
                // far more than a conditioning DIAGNOSTIC needs.
                // Tighten only with that inequality in hand
                real           mEpsilon      = 1e-4 ;
 
                // floor for the Krylov subspace handed to ( p )dnaupd. Remark 4
                // of dnaupd: raising ncv at fixed nev usually REDUCES the total
                // OP*x count, at the price of an n x ncv basis
                int_t          mSubspaceSize = 20 ;
 
                // Tolerance for the FOLDED run, which is a different job from
                // the plain one. lambda_min falls out of a cancellation,
                // lambda_min = sigma - mu, and mu is of order sigma, so an
                // error of tol*sigma in mu lands undivided on lambda_min:
                //
                //     rel_err( lambda_min )  ~  max( tol, eps_mach ) * kappa
                //
                // 1e-4 would therefore be worthless at the small end ( it caps
                // kappa at ~1e4 before the answer is all noise ), while the
                // FOLDED problem is an exterior one and can afford a tight
                // tolerance -- it converges in a handful of restarts either way.
                // The floor is eps_mach*kappa and cannot be bought off
                real           mFoldEpsilon  = 1e-10 ;
 
                // NOTE: the spectral FOLD that this class used to run for the
                // small end has been REPLACED by shift-invert
                // ( run_shift_invert ), because an affine shift cannot improve
                // eigenvalue SEPARATION and separation is what sets the
                // convergence rate -- measured, see compute_conditioning().
                // The sigma argument on the Fortran drivers is deliberately
                // KEPT: it costs nothing, the folded operator is still correct
                // where it converges, and the shift-invert shims share those
                // drivers' argument shape
 
                // ceiling on the Arnoldi basis, which is n x ncv doubles. When
                // not even the smallest legal ncv fits, the diagnostic reports
                // Infeasible rather than allocating past it
                size_t         mBasisBudget  = 512UL * 1024UL * 1024UL ;
 
                // output
                Vector< int_t > mInfo ;
 
                // the SIGNED real and imaginary part of the Ritz value that won
                // the last reduction. The function value is a magnitude, which
                // is all the caller used to need ; the fold needs the sign to
                // tell a positive definite spectrum from one straddling zero,
                // and the imaginary part to notice a complex pair. Copied from
                // the winning SLOT, so the two always describe the same value
                real            mExtremalReal = BELFEM_QUIET_NAN ;
                real            mExtremalImag = BELFEM_QUIET_NAN ;
 
                // set once an end of the spectrum is proven out of reach, and
                // never cleared: retrying it every timestep cost 4.7 s per step
                // on tapestack3d and produced nothing. Assigned ONLY from a
                // broadcast outcome -- see EigenOutcome
                bool            mDiagnosticUnavailable = false ;
 
                // what ended the last compute_conditioning(), for the caller
                // and for the one-time message
                EigenOutcome    mOutcome = EigenOutcome::Ok ;
 
                // which of the two runs is in flight, for the messages. A
                // failure that does not say whether the plain end or the fold
                // died cannot be acted on -- measured 2026-08-28, when exactly
                // that ambiguity cost a diagnostic cycle
                const char *    mRunLabel = "eigenvalue" ;
 
                // which ARPACK family drives this field. TOLD, not detected:
                // the thermal Jacobian is symmetric by construction and the
                // magnetic h-phi one is not ( Christian, 2026-08-28 ). Detecting
                // it would mean comparing A against A^T across a distribution
                // that may be transposed column blocks already ( see the CSC
                // note in run_parpack ), which is a bigger job than the caller
                // simply knowing. Defaults to the SAFE answer: the nonsymmetric
                // driver is correct for a symmetric matrix, merely slower
                bool            mSymmetric = false ;
 
                // ---- shift-invert ( mode 3 ) state -------------------------
                // ARPACK's reverse-communication arrays. Members, not locals:
                // the loop hands them to Fortran across many calls and dseupd
                // requires them untouched between the last step and the
                // extraction, so nothing transient may own them
                Vector< real >  mSiResid ;      // [ n ]
                Vector< real >  mSiBasis ;      // [ n * ncv ], the Lanczos basis
                Vector< real >  mSiWorkD ;      // [ 3 * n ]
                Vector< real >  mSiWorkL ;      // [ ncv * ( ncv + 8 ) ]
                Vector< int_t > mSiIparam ;     // [ 11 ]
                Vector< int_t > mSiIpntr ;      // [ 11 ]
                Vector< real >  mSiLambda ;     // [ nev ], eigenvalues of A
 
                // right-hand side and solution of one inverse application.
                // Sized once ; the solve writes into them every iteration
                Vector< real >  mSiRhs ;
                Vector< real >  mSiLhs ;
 
                // dedicated solver for the inverse operator. NOT mSolver,
                // which belongs to compute_lambda_max() and is built from the
                // parent's solver type -- this one is MUMPS specifically,
                // because it is the wrapper that can hold a factorization
                // still and solve against it repeatedly
                Solver *        mShiftInvertSolver = nullptr ;
 
                // the tolerance the last run was actually given. configure()
                // picks between mEpsilon and mFoldEpsilon and an explicit
                // set_tolerance() overrides both, so recomputing that choice at
                // the point of use duplicates the rule and lets the two drift.
                // Recorded once, read where it is reported
                real            mEffectiveTolerance = BELFEM_QUIET_NAN ;
 
                // set by a driver when configure() found no legal ncv. It is
                // the difference between "this iterate did not converge", which
                // gets a few strikes, and "no run of this size can ever be
                // attempted", which latches at once -- and a bare NaN cannot
                // tell the two apart. Cleared at the top of every run
                bool            mSubspaceInfeasible = false ;
 
                // consecutive failures. Every outcome but Infeasible gets a few
                // strikes ( a non-positive iterate may be transient: the Jacobian
                // is rebuilt every step ); Infeasible latches at once because it
                // depends on size and budget, not on the numbers. Every rank
                // counts the same broadcast outcome
                int_t           mFailureCount = 0 ;
                int_t           mMaxFailures  = 3 ;
 
                // which driver produced the last result. It depends on the
                // rank count AND on the build ( multi-rank without PARPACK
                // still runs ARPACK on the master ), so a caller cannot
                // derive it and has to be told
                string          mBackendLabel = "ARPACK" ;
 
                // work data
                string          mWhich;
                real mLambdaMax = BELFEM_QUIET_NAN ;
 
                Vector< real >  mLambdaReal ;
                Vector< real >  mLambdaImag ;
 
                Vector< real >  mXreal ;
                Vector< real >  mYreal ;
                Vector< real >  mZreal ;
 
                Vector< real >  mXimag ;
                Vector< real >  mYimag ;
                Vector< real >  mZimag ;
 
                Vector< cplx >  mX ;
                Vector< cplx >  mY ;
                Vector< cplx >  mZ ;
 
                bool mMatrixFlag = false ;
                bool mFirstRun = true ;
 
                Solver * mSolver = nullptr ;
 
                int_t mOriginalBase = 0 ;
 
#ifdef BELFEM_PARPACK
                // row-distributed view of the master's Jacobian. The master
                // owns the full matrix ; every other rank holds only its
                // assembly submatrix, so the row blocks PARPACK needs have to
                // be scattered from rank 0
                sparse::DistMatrixCSR< int_t > * mDistMatrix = nullptr ;
 
                // one-based copies of the distributed pattern. DistMatrixCSR
                // is zero-based by construction ( distribute_values() forces
                // SpMatrixIndexingBase::Cpp on the source ), and its arrays
                // are handed out const, so the shift lives here
                Vector< int_t > mParpackPointers ;
                Vector< int_t > mParpackIndices ;
 
                // structural fingerprint of the matrix the pattern was built
                // from. NOT reset() -- that runs after every Jacobian
                // computation, so rebuilding there would redistribute the
                // pattern once per Newton iteration
                const SpMatrix * mDistSource = nullptr ;
                int_t mDistNumRows = 0 ;
                int_t mDistNumNonzeros = 0 ;
 
#endif
 
//------------------------------------------------------------------------------
            public:
//------------------------------------------------------------------------------
 
                EigenValues( DofManager * aParent ) ;
 
                ~EigenValues();
 
                // NON-COPYABLE, NON-MOVABLE. This class owns four raw pointers
                // and deletes all four in its destructor -- mSolver,
                // mShiftInvertSolver, mM and ( under PARPACK ) mDistMatrix --
                // so the implicit copy would be a shallow pointer copy and the
                // second destructor a double free. Nothing copies it today
                // ( DofManager holds it by pointer and news it once ), which is
                // exactly why the hazard was latent rather than loud.
                //
                // Same shape and same remedy as SpMatrix ( cl_SpMatrix.hpp )
                // and Solver ( cl_Solver.hpp ). Deleting them makes any future attempt
                // a compile error instead of a run-time double free
                EigenValues( const EigenValues & ) = delete ;
                EigenValues( EigenValues && ) = delete ;
                EigenValues & operator=( const EigenValues & ) = delete ;
                EigenValues & operator=( EigenValues && ) = delete ;
 
//------------------------------------------------------------------------------
 
                void
                set_num_minvals( const index_t aNumMinVals );
 
//------------------------------------------------------------------------------
 
                void
                set_num_maxvals( const index_t aNumMaxVals );
 
//------------------------------------------------------------------------------

                void
                set_tolerance( const real aTolerance );

                void
                set_subspace_size( const index_t aSubspaceSize );

                void
                set_max_iterations( const index_t aNumMaxIter );

                void
                set_symmetric( const bool aSymmetric );

                bool
                is_symmetric() const ;
 
//------------------------------------------------------------------------------
 
                real
                compute_conditioning() ;
 
//------------------------------------------------------------------------------
 
                real
                compute_lambda_max() ;
 
//------------------------------------------------------------------------------

                real
                compute_smallest_eigenvalues();
 
                real
                compute_largest_eigenvalues();
 
//------------------------------------------------------------------------------

                const Vector< real > &
                lambda_real() const ;
 
                const Vector< real > &
                lambda_imag() const ;

                int_t
                number_of_converged_values() const ;

                const string &
                backend_label() const ;

                EigenOutcome
                outcome() const ;
 
//------------------------------------------------------------------------------
 
                void
                reset();
 
//------------------------------------------------------------------------------
            private:
//------------------------------------------------------------------------------
 
                void
                compute_matrices() ;

                bool
                configure( const int_t aJob,
                           const bool  aFolded,
                           const int_t aNumEigenValues,
                           int_t     & aSubspaceSize,
                           int_t     & aNumMaxIter,
                           real      & aTolerance ) ;

                real
                run( const int_t aJob, const real aSigma );
 
                void
                link_matrix();
 
                real
                run_arpack( const int_t aJob, const real aSigma );
 
                real
                run_parpack( const int_t aJob, const real aSigma );

                real
                reduce_extremum( const int_t aJob, const int_t aNumConverged );

                real
                report_unavailable( const int_t aOutcome );

                real
                run_shift_invert( int_t & aOutcome );
 
                void
                restore_indexing_base();
//------------------------------------------------------------------------------
 
            };
 
//------------------------------------------------------------------------------
 
            inline const Vector< real > &
            EigenValues::lambda_real() const
            {
                return mLambdaReal ;
            }
 
//------------------------------------------------------------------------------
 
            inline const Vector< real > &
            EigenValues::lambda_imag() const
            {
                return mLambdaImag ;
            }
 
//------------------------------------------------------------------------------
 
            inline int_t
            EigenValues::number_of_converged_values() const
            {
                return mInfo( arpack::gInfoNumConverged ) ;
            }
 
//------------------------------------------------------------------------------
 
            inline const string &
            EigenValues::backend_label() const
            {
                return mBackendLabel ;
            }
 
//------------------------------------------------------------------------------
 
            inline bool
            EigenValues::is_symmetric() const
            {
                return mSymmetric ;
            }
 
//------------------------------------------------------------------------------
 
            inline EigenOutcome
            EigenValues::outcome() const
            {
                return mOutcome ;
            }
 
//------------------------------------------------------------------------------
        }
    }
}
#endif //CL_FEM_DOFMGR_EIGENVALUES_HPP