cl_SolverMUMPS.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 BELFEM_CL_SOLVERMUMPS_HPP
#define BELFEM_CL_SOLVERMUMPS_HPP
 
#include "cl_Cell.hpp"
#include "cl_SolverWrapper.hpp"
#include "cl_SolverParameters.hpp"
 
namespace belfem
{
    namespace solver
    {
 
        class MUMPS : public Wrapper
        {
            // parameter object
            const SolverParameters * mParams ;
 
            // Rank of HOST
            const proc_t mMasterRank ;
 
            // id of solver
            int_t mSolverID = 0 ;
 
            // vector containing user parameters
            Vector< int_t >  mIParameters ;
            Vector< real >   mRParameters ;
 
            // vector containing debug information ( rank-local INFO: a
            // failing rank holds the true code, the others the propagated -1;
            // still read by check_warnings(), which is rank-local by design )
            Vector< int_t >  mInfo ;
 
            // rank-uniform INFOG: same code and supplementary value on every
            // rank, so the -9 retry ladder and the error arms key on it --
            // a branch taken from mInfo would strand ranks outside the
            // collective call. Vector< int_t > to match its sibling mInfo,
            // not a linear-algebra object
            Vector< int_t >  mInfoG ;
 
            Vector< real > mRInfoG ;
 
            // pointer to matrix
            SpMatrix * mMatrix = nullptr ;

            struct FrozenFactorization
            {
                bool             mArmed    = false ;

                bool             mHaveIdentity = false ;
 
                const SpMatrix * mSource   = nullptr ;
                index_t          mNumRows  = 0 ;
                index_t          mNumCols  = 0 ;
                index_t          mNumNonZeros = 0 ;
                const real     * mValues   = nullptr ;
                const int_t    * mPointers = nullptr ;
            };
 
            FrozenFactorization mFrozen ;

            int_t mMeasuredBudgetMB = 0 ;

            bool mInitialized = false ;

            SymmetryMode mSymmetryMode = SymmetryMode::Unsymmetric ;
 
//------------------------------------------------------------------------------
        public:
//------------------------------------------------------------------------------
 
            MUMPS( const SolverParameters * aParams, const proc_t aMasterRank = 0 );
 
//------------------------------------------------------------------------------
 
            ~MUMPS() override ;
 
//------------------------------------------------------------------------------

            bool
            supports_factorization_reuse() const override ;
 
            void
            freeze_factorization( const SpMatrix & aMatrix ) override ;
 
            void
            unfreeze_factorization() override ;
 
            bool
            factorization_is_frozen() const override ;
 
//------------------------------------------------------------------------------
        private:
//------------------------------------------------------------------------------

            int_t
            select_job( SpMatrix & aMatrix );

            void
            invalidate_factorization();
 
//------------------------------------------------------------------------------

            bool
            escalate_workspace( int_t & aJob );

            int_t
            memory_budget_mb();
 
//------------------------------------------------------------------------------
        public:
//------------------------------------------------------------------------------
 
//------------------------------------------------------------------------------
 
            void
            solve(
                    SpMatrix & aMatrix,
                    Vector< real > & aLHS,
                    Vector< real > & aRHS ) override ;
 
//------------------------------------------------------------------------------
 
            void
            solve(
                    SpMatrix & aMatrix,
                    Matrix< real > & aLHS,
                    Matrix< real > & aRHS ) override ;
 
//------------------------------------------------------------------------------
 
            void
            free() override ;
 
//------------------------------------------------------------------------------
 
            void
            set_reordering(
                    const MumpsSerialReodrdering   aSerial,
                    const MumpsParallelReodrdering aParallel
            );
 
//------------------------------------------------------------------------------
 
            void
            set_block_low_ranking(
                    const MumpsBlockLowRanking aBLK,
                    const real            aEpsilon
            );
 
//------------------------------------------------------------------------------
 
            void
            set_error_analysis(
                    const MumpsErrorAnalysis aSetting
            );
 
//------------------------------------------------------------------------------

            real
            get_determinant() const override ;
 
//------------------------------------------------------------------------------

            real
            get_cond1() const override ;
 
//------------------------------------------------------------------------------

            real
            get_cond2() const override ;
 
//------------------------------------------------------------------------------

            real
            get_forward_error() const override ;
 
//------------------------------------------------------------------------------

            real
            get_backward_error() const override ;
 
//------------------------------------------------------------------------------

            real
            get_omega2() const override ;
 
//------------------------------------------------------------------------------
        protected :
//------------------------------------------------------------------------------
 
            void
            initialize(
                    SpMatrix & aMatrix,
                    const SymmetryMode aSymmetryMode = SymmetryMode::Unsymmetric,
                    const int_t aNumRhsColumns=1 ) override ;
 
 
//------------------------------------------------------------------------------
 
            string
            error_message(
                    const int_t   * aInfo,
                    const int_t   & aN,
                    const int_t   & aNNZ ) ;
 
//------------------------------------------------------------------------------
 
            // decodes a positive INFO(1) warning ( a usable solution is still
            // returned ) into one human-readable line per active warning bit
            void
            warning_message(
                    const int_t    * aInfo,
                    Cell< string > & aWarnings ) ;
 
//------------------------------------------------------------------------------
        private :
//------------------------------------------------------------------------------
 
            void
            init_defaults();
 
//------------------------------------------------------------------------------
 
            // handles a positive INFO(1) after a solve: raises a hard error on
            // the out-of-range-index bit ( +1 ) and reports the remaining
            // warnings per rank. NOT used by free(), which is teardown
            void
            check_warnings();
 
//------------------------------------------------------------------------------
 
            // rank-0 log line for the soft-fail path, box-styled to match the
            // controller's timestep box, with a compact INFO decode ( the full
            // decode lives in error_message() on the hard path )
            void
            print_soft_fail();
 
//------------------------------------------------------------------------------
        };
    }
}
#endif //BELFEM_CL_SOLVERMUMPS_HPP