cl_AR_Matrix.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_AR_MATRIX_HPP
#define BELFEM_CL_AR_MATRIX_HPP
 
#include "typedefs.hpp"
#include "assert.hpp"
 
namespace belfem
{
    template< typename T >
    class Matrix
    {
//------------------------------------------------------------------------------
    public :
//------------------------------------------------------------------------------
 
        typedef arma::Mat<T> MatrixType;
 
//------------------------------------------------------------------------------
    private :
//------------------------------------------------------------------------------
 
        // member class of underlying vector implementation
        MatrixType mMatrix;
 
//------------------------------------------------------------------------------
    public :
//------------------------------------------------------------------------------

        Matrix() = default;
 
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

        Matrix( const size_t aNumRows,
                   const size_t aNumCols ) :
                mMatrix( aNumRows, aNumCols ) {}
 
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

        Matrix( const size_t aNumRows,
                   const size_t aNumCols,
                   const real aValue ) :
                mMatrix( aNumRows, aNumCols )
        {
            this->fill( aValue );
        }
 
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

        Matrix( const std::initializer_list<std::initializer_list<T> > & aInitList )
                : mMatrix( aInitList )
        {
 
        }
 
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

        Matrix( const MatrixType & aExpression )
                : mMatrix( aExpression ) {}
 
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

        template< typename ET, typename OP >
        Matrix( const arma::Op<ET, OP> & aExpression )
                : mMatrix( aExpression ) {}
 
//------------------------------------------------------------------------------

        Matrix( const Matrix< T > & aMatrix ) :
            mMatrix( aMatrix.mMatrix )
        {}

        Matrix( Matrix< T > && aMatrix ) noexcept :
                mMatrix( std::move( aMatrix.mMatrix ) )
        {}
 
 
//------------------------------------------------------------------------------

        virtual ~Matrix() = default;
 
//------------------------------------------------------------------------------
// MEMORY
//------------------------------------------------------------------------------
 
        /*
         * expose the underlying raw pointer
         */
 
        inline T *
        data()
        {
            return mMatrix.memptr();
        }
 
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
        /*
         * expose the underlying raw pointer ( const version )
         */
        inline const T *
        data() const
        {
            return mMatrix.memptr();
        }
 
//------------------------------------------------------------------------------
 
        /*
         * expose the underlying matrix implementation ( writable version )
         */
        inline MatrixType &
        matrix_data()
        {
            return mMatrix;
        }
 
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
        /*
         * expose the underlying matrix implementation ( const version )
         */
        inline const MatrixType &
        matrix_data() const
        {
            return mMatrix;
        }
 
//------------------------------------------------------------------------------
// UTILITIES
//------------------------------------------------------------------------------
 
        inline void
        fill( const T & aValue )
        {
            mMatrix.fill( aValue );
        }
 
//------------------------------------------------------------------------------
 
        inline void
        set_size( const size_t aNumRows, const size_t aNumCols )
        {
            mMatrix.set_size( aNumRows, aNumCols );
        }
 
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
        inline void
        set_size( const size_t aNumRows, const size_t aNumCols, const T & aValue )
        {
            this->set_size( aNumRows, aNumCols );
            this->fill( aValue );
        }
 
//------------------------------------------------------------------------------
// SIZE CHECKS
//------------------------------------------------------------------------------
 
        inline size_t
        n_rows() const
        {
            return mMatrix.n_rows;
        }
 
//------------------------------------------------------------------------------
 
        inline size_t
        n_cols() const
        {
            return mMatrix.n_cols;
        }
 
//------------------------------------------------------------------------------

        inline size_t
        capacity() const
        {
            return mMatrix.n_elem;
        }
 
//------------------------------------------------------------------------------

        inline size_t
        spacing() const
        {
            return mMatrix.n_rows;
        }
 
//------------------------------------------------------------------------------
// OPERATORS
//------------------------------------------------------------------------------

        Matrix< T > &
        operator=( Matrix< T > && aMatrix ) noexcept
        {
            if ( this != & aMatrix )
            {
                mMatrix = std::move( aMatrix.mMatrix );
            }
            return *this;
        }

        Matrix< T > &
        operator=( const Matrix< T > & aMatrix )
        {
            if ( this != &aMatrix )
            {
                mMatrix = aMatrix.mMatrix;
            }
            return *this;
        }
 
 
//------------------------------------------------------------------------------

        inline T &
        operator()( const size_t aRowIndex, const size_t aColIndex )
        {
            BELFEM_ASSERT( aRowIndex< this->n_rows(),
                          "Row index %lu out of bounds, which must be smaller than %lu.",
                          ( long unsigned int ) aRowIndex,
                          ( long unsigned int ) this->n_rows() );
 
            BELFEM_ASSERT( aColIndex< this->n_cols(),
                          "Col index %lu out of bounds, which must be smaller than %lu.",
                          ( long unsigned int ) aColIndex,
                          ( long unsigned int ) this->n_cols() );
 
            return mMatrix( aRowIndex, aColIndex );
        }
 
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

        virtual inline const T &
        operator()( const size_t aRowIndex, const size_t aColIndex ) const
        {
            BELFEM_ASSERT( aRowIndex< this->n_rows(),
                          "Row index %lu out of bounds, which must be smaller than %lu.",
                          ( long unsigned int ) aRowIndex,
                          ( long unsigned int ) this->n_rows() );
 
            BELFEM_ASSERT( aColIndex< this->n_cols(),
                          "Col index %lu out of bounds, which must be smaller than %lu.",
                          ( long unsigned int ) aColIndex,
                          ( long unsigned int ) this->n_cols() );
 
            return mMatrix( aRowIndex, aColIndex );
        }
 
//------------------------------------------------------------------------------
// Cols and Rows
//------------------------------------------------------------------------------
 
        inline auto
        row( const size_t aRowIndex )
            ->decltype( mMatrix.row( aRowIndex ) )
        {
            BELFEM_ASSERT( aRowIndex < this->n_rows(),
                          "Row index %lu out of bounds, which must be less than %lu.",
                          ( long unsigned int ) aRowIndex,
                          ( long unsigned int ) this->n_rows() );
 
            return  mMatrix.row( aRowIndex );
        }
 
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
        inline auto
        row( const size_t aRowIndex ) const
            ->decltype( mMatrix.row( aRowIndex ) )
        {
            BELFEM_ASSERT( aRowIndex < this->n_rows(),
                          "Row index %lu out of bounds, which must be less than %lu.",
                          ( long unsigned int ) aRowIndex,
                          ( long unsigned int ) this->n_rows() );
 
            return  mMatrix.row( aRowIndex );
        }
 
//------------------------------------------------------------------------------
 
        inline auto
        col( const size_t aColIndex )
            ->decltype( mMatrix.col( aColIndex ) )
        {
            BELFEM_ASSERT( aColIndex < this->n_cols(),
                          "Col index %lu out of bounds, which must be less than %lu.",
                          ( long unsigned int ) aColIndex,
                          ( long unsigned int ) this->n_cols() );
 
            return  mMatrix.col( aColIndex );
        }
 
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
        inline auto
        col( const size_t aColIndex ) const
            ->decltype( mMatrix.col( aColIndex ) )
        {
            BELFEM_ASSERT( aColIndex < this->n_cols(),
                          "Col index %lu out of bounds, which must be less than %lu.",
                          ( long unsigned int ) aColIndex,
                          ( long unsigned int ) this->n_cols() );
 
            return  mMatrix.col( aColIndex );
        }
 
//------------------------------------------------------------------------------
 
        inline auto
        submat( const size_t aFirstRow,
                const size_t aFirstCol,
                const size_t aLastRow,
                const size_t aLastCol )
            -> decltype( mMatrix.submat( aFirstRow, aFirstCol, aLastRow, aLastCol ) )
        {
            BELFEM_ASSERT( aFirstRow < this->n_rows(),
                          "First row index %lu out of bounds, which must be less than %lu.",
                          ( long unsigned int ) aFirstRow,
                          ( long unsigned int ) this->n_rows() );
 
            BELFEM_ASSERT( aLastRow < this->n_rows(),
                          "Last row index %lu out of bounds, which must be less than %lu.",
                          ( long unsigned int ) aLastRow,
                          ( long unsigned int ) this->n_rows() );
 
            BELFEM_ASSERT( aFirstCol < this->n_cols(),
                          "First col index %lu out of bounds, which must be less than %lu.",
                          ( long unsigned int ) aFirstCol,
                          ( long unsigned int ) this->n_cols() );
 
            BELFEM_ASSERT( aLastCol < this->n_cols(),
                          "Last col index %lu out of bounds, which must be less than %lu.",
                          ( long unsigned int ) aLastCol,
                          ( long unsigned int ) this->n_cols() );
 
            return mMatrix.submat( aFirstRow, aFirstCol, aLastRow, aLastCol );
        }
 
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
        inline auto
        submat( const size_t aFirstRow,
                const size_t aFirstCol,
                const size_t aLastRow,
                const size_t aLastCol ) const
            -> decltype( mMatrix.submat( aFirstRow, aFirstCol, aLastRow, aLastCol ) )
        {
            BELFEM_ASSERT( aFirstRow < this->n_rows(),
                          "First row index %lu out of bounds, which must be less than %lu.",
                          ( long unsigned int ) aFirstRow,
                          ( long unsigned int ) this->n_rows() );
 
            BELFEM_ASSERT( aLastRow < this->n_rows(),
                          "Last row index %lu out of bounds, which must be less than %lu.",
                          ( long unsigned int ) aLastRow,
                          ( long unsigned int ) this->n_rows() );
 
            BELFEM_ASSERT( aFirstCol < this->n_cols(),
                          "First col index %lu out of bounds, which must be less than %lu.",
                          ( long unsigned int ) aFirstCol,
                          ( long unsigned int ) this->n_cols() );
 
            BELFEM_ASSERT( aLastCol < this->n_cols(),
                          "Last col index %lu out of bounds, which must be less than %lu.",
                          ( long unsigned int ) aLastCol,
                          ( long unsigned int ) this->n_cols() );
 
            return mMatrix.submat( aFirstRow, aFirstCol, aLastRow, aLastCol );
        }
 
//------------------------------------------------------------------------------
 
        inline void
        set_row( const size_t aRowIndex, const Vector <T> & aVector )
        {
            // make sure that length is correct
            BELFEM_ASSERT( aVector.length() == this->n_cols(),
                          "Wrong length of vector : %lu ( expect %lu )",
                          ( long unsigned int ) aVector.length(),
                          ( long unsigned int ) this->n_cols() );
 
            this->row( aRowIndex ) = arma::trans( aVector.vector_data() );
        }
 
        template < typename ET >
        inline void
        set_row( const size_t aColIndex, const ET & aExpression )
        {
            this->row( aColIndex ) = aExpression;
        }
 
//------------------------------------------------------------------------------
 
        inline void
        set_col( const size_t aColIndex, const Vector <T> & aVector )
        {
            // make sure that length is correct
            BELFEM_ASSERT( aVector.length() == this->n_rows(),
                          "Wrong length of vector : %lu ( expect %lu )",
                          ( long unsigned int ) aVector.length(),
                          ( long unsigned int ) this->n_rows() );
 
            this->col( aColIndex ) = aVector.vector_data();
        }
 
        template < typename ET >
        inline void
        set_col( const size_t aColIndex, const ET & aExpression )
        {
            this->col( aColIndex ) = aExpression;
        }
 
 
//------------------------------------------------------------------------------
// internal operators
//------------------------------------------------------------------------------

        Matrix< T > &
        operator=( const T & aValue )
        {
            this->fill( aValue );
            return *this;
        }
 
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

        template < typename ET >
        Matrix< T > &
        operator=( const ET & aExpression )
        {
            mMatrix = aExpression;
            return *this;
        }
 
//------------------------------------------------------------------------------
// ADD OPERATORS
//------------------------------------------------------------------------------
 
        Matrix< T > &
        operator+=( const T & aValue )
        {
            mMatrix += aValue;
            return *this;
        }
 
        Matrix< T > &
        operator+=( T & aValue )
        {
            mMatrix += aValue;
            return *this;
        }
 
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
        Matrix< T > &
        operator+=( const Matrix< T > & aMatrix )
        {
            mMatrix += aMatrix.matrix_data();
            return *this;
        }
 
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
        Matrix< T > &
        operator+=( const MatrixType & aExpression )
        {
            mMatrix += aExpression;
            return *this;
        }
 
//------------------------------------------------------------------------------
// SUBTRACT OPERATORS
//------------------------------------------------------------------------------
 
        Matrix< T > &
        operator-=( const T & aValue )
        {
            mMatrix -= aValue;
            return *this;
        }
 
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
        Matrix< T > &
        operator-=( const Matrix< T > & aMatrix )
        {
            mMatrix -= aMatrix.matrix_data();
            return *this;
        }
 
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
        Matrix< T > &
        operator-=( const MatrixType & aExpression )
        {
            mMatrix -= aExpression;
            return *this;
        }
 
//------------------------------------------------------------------------------
// MULTIPLY OPERATORS
//------------------------------------------------------------------------------
 
        Matrix< T > &
        operator*=( const T & aValue )
        {
            mMatrix *= aValue;
            return *this;
        }
 
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
        Matrix< T > &
        operator*=( const Matrix< T > & aMatrix )
        {
            mMatrix *= aMatrix.matrix_data();
            return *this;
        }
 
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
        Matrix< T > &
        operator*=( const MatrixType & aExpression )
        {
            mMatrix *= aExpression;
            return *this;
        }
 
//------------------------------------------------------------------------------
// DIVIDE OPERATORS
//------------------------------------------------------------------------------
 
        Matrix< T > &
        operator/=( const T & aValue )
        {
            mMatrix /= aValue;
            return *this;
        }
 
//------------------------------------------------------------------------------
// PRINT OPERATION
//------------------------------------------------------------------------------
 
        inline void
        print( const std::string aLabel="Matrix" ) const;
 
//------------------------------------------------------------------------------
    };
 
//------------------------------------------------------------------------------
// turn off annoying waning
#ifdef BELFEM_GCC
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wformat"
#endif
    template < typename T > void
    Matrix< T >::print( const std::string aLabel ) const
    {
        FILE * tOutFile = stdout;
 
        fprintf( tOutFile, "%s = [ ... \n", aLabel.c_str() );
 
        // get a ref of this
        const Matrix< T > & tThis = *this;
 
        uint tRows = tThis.n_rows();
        uint tCols = tThis.n_cols();
 
        for( uint i=0; i< tRows; ++i )
        {
            for( uint j=0; j< tCols; ++j )
            {
                if( j < tCols-1 )
                {
                    fprintf( tOutFile, "%d, ", ( int ) tThis( i, j ) );
                }
                else
                {
                    fprintf( tOutFile, "%d; ", ( int ) tThis( i, j ) );
                }
            }
 
            if( i < tRows-1 )
            {
                fprintf( tOutFile, "...\n" );
            }
            else
            {
                fprintf( tOutFile, "];\n" );
            }
        }
    }
 
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
    template <>
    inline void
    Matrix<real>::print( const std::string aLabel ) const
    {
        FILE * tOutFile = stdout;
 
        fprintf( tOutFile, "%s = [ ... \n", aLabel.c_str() );
 
        const Matrix< real > & tThis = *this;
 
        uint tRows = tThis.n_rows();
        uint tCols = tThis.n_cols();
 
        for( uint i=0; i< tRows; ++i )
        {
            for( uint j=0; j< tCols; ++j )
            {
                if( j < tCols-1 )
                {
                    fprintf( tOutFile, "%+.15e, ", ( double ) tThis( i, j ) );
                }
                else
                {
                    fprintf( tOutFile, "%+.15e; ", ( double ) tThis( i, j ) );
                }
            }
 
            if( i < tRows-1 )
            {
                fprintf( tOutFile, "...\n" );
            }
            else
            {
                fprintf( tOutFile, "];\n" );
            }
        }
    }
 
 
#ifdef BELFEM_GCC
#pragma GCC diagnostic pop
#endif
}
#endif //BELFEM_CL_AR_MATRIX_HPP