cl_SpMatrix.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_SPMATRIX_HPP
#define BELFEM_CL_SPMATRIX_HPP
 
#include <cstring>
 
#include "typedefs.hpp"
#include "cl_Cell.hpp"
#include "cl_Graph_Vertex.hpp"
#include "cl_Vector.hpp"
#include "filetools.hpp"
#include "hdf5_tools.hpp"
 
namespace belfem
{
//------------------------------------------------------------------------------
 
    enum class SpMatrixType
    {
        CSC        = 0,         // compressed sparse column
        CSR        = 1,         // compressed sparse row
        UNDEFINED  = 2
    };
 
//------------------------------------------------------------------------------
 
    enum class SpMatrixIndexingBase
    {
        Cpp        = 0,
        Fortran    = 1
    };
 
//------------------------------------------------------------------------------

    class SpMatrix
    {
//------------------------------------------------------------------------------
    private:
//------------------------------------------------------------------------------
 
        SpMatrix * mParent = nullptr ;
        SpMatrix * mChild  = nullptr ;
 
        // type of matrix, CSC or CSR
        SpMatrixType mType;
 
        // size of matrix
        int_t mNumRows = 0;
        int_t mNumCols = 0;
 
        // number of nonzeros
        int_t mNumNonZeros = 0;
 
        // size of pointer array
        int_t mPointerSize = 0;
 
        // container for pointers
        int_t * mPointers = nullptr;
 
        int_t * mRows = nullptr;
 
        int_t * mColumns = nullptr;
 
        // values array
        real * mValues = nullptr;
 
        bool mHaveCooIndices = false;
 
#ifdef BELFEM_NETLIB
        Vector< real > mSwap ;
        int_t mSwapSize = 0;
#endif
 
        // pointer with zero value
        real mZero = 0.0;
 
        // Note: the lookup used to dispatch through a member-function
        // pointer over four near-identical search functions, one per
        // ( type, base ) pair. It is now the inline position() below, which
        // reads the base from mPointers[ 0 ] on every call. That is not only
        // faster - it removes a stale-state bug: a parent and its child share
        // the index arrays, and rebasing the parent did not update a child's
        // cached function pointer.
 
//------------------------------------------------------------------------------
    public:
//------------------------------------------------------------------------------
 
        // empty constructor
        SpMatrix() = default;
 
        // raw-owning class: copy and move construction are not supported,
        // use the assignment operators on an empty matrix instead
        SpMatrix( const SpMatrix & ) = delete;
        SpMatrix( SpMatrix && ) = delete;
 
//------------------------------------------------------------------------------
 
        // standard constructor using a graph
        SpMatrix( Graph & aGraph,
                  const enum SpMatrixType aType = SpMatrixType::CSC,
                          const index_t   aNumRows = 0,
                          const index_t   aNumCols = 0,
                          const bool aSortGraph = true );
 
//------------------------------------------------------------------------------
 
        SpMatrix( SpMatrixType aType,
                  const index_t aNumRows,
                  const index_t aNumCols,
                  const index_t aNumNonZeros,
                  const int_t  * aIndices,
                  const int_t  * aPointers );
 
//------------------------------------------------------------------------------
 
        // file constructor using a HDF5 file
        SpMatrix( const string & aHDF5Path, const string   aLabel="Matrix" );
 
//------------------------------------------------------------------------------
 
        SpMatrix( const hid_t aParent, const string aLabel="Matrix" );
 
//------------------------------------------------------------------------------
 
        // constructor for testing purposes using dense Matrix
        SpMatrix( const Matrix< real > & aMatrix, const SpMatrixType aType = SpMatrixType::CSC );
 
//------------------------------------------------------------------------------
 
        // child constructor: shares the sparsity structure of the parent,
        // owns only its value array
        explicit SpMatrix( SpMatrix * aParent );
 
//------------------------------------------------------------------------------
 
        ~SpMatrix();
 
//------------------------------------------------------------------------------
// Access to data containers
//------------------------------------------------------------------------------

        const SpMatrixType &
        type() const;
 
//------------------------------------------------------------------------------

         void
         fill( const real aValue );
 
//------------------------------------------------------------------------------

        index_t
        n_rows() const;
 
//------------------------------------------------------------------------------

        index_t
        n_cols() const;
 
//------------------------------------------------------------------------------
 
        index_t
        n_pointers() const ;
 
//------------------------------------------------------------------------------
 
        size_t
        memory() const ;
 
//------------------------------------------------------------------------------
 
        bool
        have_coo_indices() const ;
 
//------------------------------------------------------------------------------

        void
        set_type( const SpMatrixType aType );
 
//------------------------------------------------------------------------------

        index_t
        number_of_nonzeros() const;
 
//------------------------------------------------------------------------------

        int_t *
        pointers();
 
//------------------------------------------------------------------------------

        const int_t *
        pointers() const;
 
//------------------------------------------------------------------------------

         int_t *
         indices();
 
//------------------------------------------------------------------------------

        const int_t *
        indices() const;
 
//------------------------------------------------------------------------------

        int_t *
        rows();
 
//------------------------------------------------------------------------------

        const int_t *
        rows() const;
 
//------------------------------------------------------------------------------

        int_t *
        cols();
 
//------------------------------------------------------------------------------

        const int_t *
        cols() const;
 
//------------------------------------------------------------------------------

        real *
        data();
 
//------------------------------------------------------------------------------

        const real *
        data() const;
 
//------------------------------------------------------------------------------

        real &
        data( const index_t aIndex );
 
//------------------------------------------------------------------------------

        const real &
        data( const index_t aIndex ) const;
 
 
//------------------------------------------------------------------------------
 
        /*
         * get the index of a specific row and col
         *
         * Returns the position into the value array, or number_of_nonzeros()
         * when the entry is not part of the sparsity pattern. Thin wrapper
         * around position().
         * */
        int_t
        index( const index_t aRowIndex,
               const index_t aColIndex ) const;
 
//------------------------------------------------------------------------------

        int_t
        position( const index_t aRowIndex,
                  const index_t aColIndex ) const;
 
//------------------------------------------------------------------------------

        void
        positions_in_slice( const index_t   aSlice,
                            const int_t   * aCols,
                            const uint      aNumCols,
                            int_t         * aPos ) const;
 
//------------------------------------------------------------------------------
// Utilities
//------------------------------------------------------------------------------

        void
        sort_entries();
 
//------------------------------------------------------------------------------

        void
        set_indexing_base( const enum SpMatrixIndexingBase & aBasis );
 
//------------------------------------------------------------------------------

        void
        create_coo_indices();
 
//------------------------------------------------------------------------------

        void
        free_coo_indices();
 
//------------------------------------------------------------------------------

        void
        print( const string aLabel="SparseMatrix" );
 
//------------------------------------------------------------------------------

        void
        print2( const string aLabel="SparseMatrix" );
 
//------------------------------------------------------------------------------

        int_t
        indexing_base() const;
 
//------------------------------------------------------------------------------

        void
        multiply( const Vector< real > & aX,
                        Vector< real > & aY,
                  const real aAlpha,
                  const real aBeta,
                  const bool aTransposedFlag=false );
 
//------------------------------------------------------------------------------

        void
        multiply( const Vector< real > & aX,
                        Vector< real > & aY );
 
//------------------------------------------------------------------------------
 
        void
        transpose();
 
//------------------------------------------------------------------------------
// Saving and Loading
//------------------------------------------------------------------------------

        void
        save(   const string & aPath,
                const string   aLabel="Matrix",
                const enum FileMode aMode=FileMode::NEW );
 
//------------------------------------------------------------------------------

        void
        save(   hid_t        & aGroup,
                herr_t       & aStatus );
 
//------------------------------------------------------------------------------

        void
        load(   const string & aPath,
                const string   aLabel="Matrix"
                );
 
//------------------------------------------------------------------------------

        void
        load(   hid_t        & aGroup,
                herr_t       & aStatus );
 
 
//------------------------------------------------------------------------------
// Operators
//------------------------------------------------------------------------------

        real &
        operator()( const index_t & aRowIndex,
                    const index_t & aColIndex );
 
//------------------------------------------------------------------------------

        const real &
        operator()( const index_t & aRowIndex,
                    const index_t & aColIndex ) const;
//------------------------------------------------------------------------------

        SpMatrix &
        operator=( const SpMatrix & aMatrix );
 
//------------------------------------------------------------------------------

         SpMatrix &
         operator=( SpMatrix && aMatrix );
 
//------------------------------------------------------------------------------
    private:
//------------------------------------------------------------------------------
// Initialization
//------------------------------------------------------------------------------

         void
         check_graph( Graph & aGraph );
 
//------------------------------------------------------------------------------

        index_t
        tidy_graph( Graph & aGraph );
 
//------------------------------------------------------------------------------
 
        void
        create_csr_indices( Graph & aGraph );
 
//------------------------------------------------------------------------------
 
        void
        create_csc_indices( Graph & aGraph );
 
//------------------------------------------------------------------------------

        void
        set_sizes(
                const index_t aNumRows,
                const index_t aNumCols );
 
        void
        set_nnz( const index_t aNumberOfNonzeros );
 
//------------------------------------------------------------------------------

        void
        deallocate();
 
//------------------------------------------------------------------------------

        void
        allocate_values();
 
//------------------------------------------------------------------------------
 
#ifdef BELFEM_NETLIB
        void
        allocate_swap();
#endif
 
//------------------------------------------------------------------------------
// Parent/Child logic
//------------------------------------------------------------------------------
 
        void
        set_child( SpMatrix * aChild );
 
//------------------------------------------------------------------------------
 
        void
        update_from_parent();
 
//------------------------------------------------------------------------------
    };
 
 
//------------------------------------------------------------------------------
// external operators
//------------------------------------------------------------------------------

    Vector< real >
    inline operator * ( SpMatrix     & aA,
                        const Vector<real> & aX )
    {
        Vector<real> aY( aA.n_rows(), 0.0 );
        aA.multiply( aX, aY );
        return aY ;
    }
 
//------------------------------------------------------------------------------
 
    inline const SpMatrixType &
    SpMatrix::type() const
    {
        return mType;
    }
 
//------------------------------------------------------------------------------
 
    inline index_t
    SpMatrix::n_rows() const
    {
        return ( index_t ) mNumRows;
    }
 
//------------------------------------------------------------------------------
 
    inline index_t
    SpMatrix::n_cols() const
    {
        return ( index_t ) mNumCols;
    }
 
//------------------------------------------------------------------------------
 
    inline index_t
    SpMatrix::n_pointers() const
    {
        return ( index_t ) mPointerSize;
    }
 
//------------------------------------------------------------------------------
 
    inline bool
    SpMatrix::have_coo_indices() const
    {
        return mHaveCooIndices;
    }
 
//------------------------------------------------------------------------------
 
    inline index_t
    SpMatrix::number_of_nonzeros() const
    {
        return ( index_t ) mNumNonZeros;
    }
 
//------------------------------------------------------------------------------
 
    inline int_t *
    SpMatrix::indices()
    {
        if( mType == SpMatrixType::CSC )
        {
            return mRows;
        }
        else if ( mType == SpMatrixType::CSR )
        {
            return mColumns;
        }
        else
        {
            return nullptr;
        }
    }
 
//------------------------------------------------------------------------------
 
    inline const int_t *
    SpMatrix::indices() const
    {
        if( mType == SpMatrixType::CSC )
        {
            return mRows;
        }
        else if ( mType == SpMatrixType::CSR )
        {
            return mColumns;
        }
        else
        {
            return nullptr;
        }
    }
 
//------------------------------------------------------------------------------
 
    inline int_t *
    SpMatrix::rows()
    {
        return mRows;
    }
 
//------------------------------------------------------------------------------
 
    inline const int_t *
    SpMatrix::rows() const
    {
        return mRows;
    }
 
//------------------------------------------------------------------------------
 
    inline int_t *
    SpMatrix::cols()
    {
        return mColumns;
    }
 
//------------------------------------------------------------------------------
 
    inline const int_t *
    SpMatrix::cols() const
    {
        return mColumns;
    }
 
//------------------------------------------------------------------------------
 
    inline int_t *
    SpMatrix::pointers()
    {
        return mPointers;
    }
 
//------------------------------------------------------------------------------
 
    const inline int_t *
    SpMatrix::pointers() const
    {
        return mPointers;
    }
 
//------------------------------------------------------------------------------
 
    inline real *
    SpMatrix::data()
    {
        return mValues;
    }
 
//------------------------------------------------------------------------------
 
    inline const real *
    SpMatrix::data() const
    {
        return mValues;
    }
 
//------------------------------------------------------------------------------
 
    inline real &
    SpMatrix::data( const index_t aIndex )
    {
        BELFEM_ASSERT( aIndex < ( index_t ) mNumNonZeros,
            "Index %lu for sparse matix out of bounds ( must be less than %lu )",
                      ( long unsigned int ) aIndex,
                      ( long unsigned int ) mNumNonZeros );
 
        return mValues[ aIndex ];
    }
 
//------------------------------------------------------------------------------
 
    inline const real &
    SpMatrix::data( const index_t aIndex ) const
    {
        BELFEM_ASSERT( aIndex < ( index_t ) mNumNonZeros,
                      "Index %lu for sparse matix out of bounds ( must be less than %lu )",
                      ( long unsigned int ) aIndex,
                      ( long unsigned int ) mNumNonZeros );
 
        return mValues[ aIndex ];
    }
 
//------------------------------------------------------------------------------
 
    inline int_t
    SpMatrix::indexing_base() const
    {
        return mPointers != nullptr ? mPointers[ 0 ] : 0;
    }
 
//------------------------------------------------------------------------------
 
    inline int_t
    SpMatrix::position( const index_t aRowIndex, const index_t aColIndex ) const
    {
        BELFEM_ASSERT( aRowIndex < ( index_t ) mNumRows, "aRowIndex out of bounds ( %lu >= %lu )",
                      ( long unsigned int ) aRowIndex ,
                      ( long unsigned int ) mNumRows );
 
        BELFEM_ASSERT( aColIndex < ( index_t ) mNumCols, "aColIndex out of bounds ( %lu >= %lu )",
                      ( long unsigned int ) aColIndex ,
                      ( long unsigned int ) mNumCols );
 
        // CSR searches a row for a column, CSC a column for a row. Both the
        // type and the base are loop-invariant in every caller, so these
        // branches predict perfectly and the whole function inlines.
        const int_t  tBase  = mPointers[ 0 ];
        const bool   tIsCsr = ( mType == SpMatrixType::CSR );
 
        const index_t tSlice  = tIsCsr ? aRowIndex : aColIndex ;
        const int_t   tTarget = ( int_t )( tIsCsr ? aColIndex : aRowIndex ) + tBase ;
 
        const int_t * tIndices = tIsCsr ? mColumns : mRows ;
 
        const int_t tBegin = mPointers[ tSlice ]     - tBase ;
        const int_t tEnd   = mPointers[ tSlice + 1 ] - tBase ;
 
        // Binary search, always. A linear scan with an early exit was tried
        // for short slices on the theory that FE rows are only tens of
        // entries long. On the machine it was measured on it lost at every
        // row length ( 8, 16, 24, 32, 50 ): a slice sits in one or two cache
        // lines, so lower_bound costs ~6 well-predicted steps while the scan
        // averages k/2 iterations behind a data-dependent branch.
        //
        // A second reviewer running the same experiment on different hardware
        // did see linear win at 8 entries per row, so the crossover is
        // microarchitecture-dependent and this is not a universal result. The
        // threshold knob was removed rather than guessed at; if it is ever
        // revisited, re-measure on the production nodes rather than trusting
        // either of those numbers. Both agree the difference is small next to
        // what the assembly restructure targets.
        // See devlog/dl20260818_spmatrix_accessor.md.
        //
        // ( indices within a slice are sorted, see sort_entries )
        const int_t * tFound = std::lower_bound( tIndices + tBegin,
                                                 tIndices + tEnd,
                                                 tTarget );
 
        return ( tFound < tIndices + tEnd && *tFound == tTarget )
               ? ( int_t )( tFound - tIndices )
               : mNumNonZeros ;
    }
 
//------------------------------------------------------------------------------
 
    inline int_t
    SpMatrix::index( const index_t aRowIndex, const index_t aColIndex ) const
    {
        return this->position( aRowIndex, aColIndex );
    }
 
//------------------------------------------------------------------------------

    inline real &
    SpMatrix::operator()( const index_t & aRowIndex,
                const index_t & aColIndex )
    {
        BELFEM_ASSERT( this->indexing_base() == 0,
            "operator() called while matrix is in Fortran indexing mode. "
            "Call set_indexing_base( SpMatrixIndexingBase::Cpp ) first." );
 
        const int_t tIndex = this->position( aRowIndex, aColIndex );
 
        BELFEM_ASSERT( tIndex < mNumNonZeros, "tried to access zero value in writable mode( %lu, %lu )",
                      ( long unsigned int ) aRowIndex,
                      ( long unsigned int ) aColIndex );
 
        return mValues[ tIndex ];
    }
 
//------------------------------------------------------------------------------
 
    inline const real &
    SpMatrix::operator()( const index_t & aRowIndex,
                const index_t & aColIndex ) const
    {
        BELFEM_ASSERT( this->indexing_base() == 0,
            "operator() called while matrix is in Fortran indexing mode. "
            "Call set_indexing_base( SpMatrixIndexingBase::Cpp ) first." );
 
        const int_t tIndex = this->position( aRowIndex, aColIndex );
 
        if( tIndex < mNumNonZeros )
        {
            return mValues[ tIndex ];
        }
        else
        {
            return mZero;
        }
    }
 
//------------------------------------------------------------------------------
}
 
#endif //BELFEM_CL_SPMATRIX_HPP