matrix_factory.hpp

Go to the documentation of this file.

/*
 *   Copyright 2023 Huawei Technologies Co., Ltd.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#ifndef _H_GRB_MATRIX_FACTORY
#define _H_GRB_MATRIX_FACTORY

#include <random>
#include <iostream>
#include <algorithm>

#include <graphblas/utils/iterators/adapter.hpp>
#include <graphblas/utils/iterators/regular.hpp>

#include <graphblas.hpp>


namespace grb::algorithms {

    template<
        typename D,
        grb::IOMode mode = grb::PARALLEL,
        grb::Backend backend = grb::config::default_backend,
        typename RIT = grb::config::RowIndexType,
        typename CIT = grb::config::ColIndexType,
        typename NIT = grb::config::NonzeroIndexType
    >
    class matrices {

        friend class matrices< void, mode, backend >;

        private:

            typedef Matrix< D, backend, RIT, CIT, NIT > MatrixType;

            static size_t getP() {
                return mode == grb::SEQUENTIAL ? 1 : grb::spmd<backend>::nprocs();
            }

            static size_t getPID() {
                return mode == grb::SEQUENTIAL ? 0 : grb::spmd<backend>::pid();
            }

            static size_t compute_diag_length(
                const size_t m, const size_t n,
                const long k
            ) {
                constexpr const long zero = static_cast< long >( 0 );
                const size_t k_abs = static_cast< size_t >(
                    (k < zero) ? -k : k );
                // catch out-of-bounds offsets
                if( k < zero && k_abs >= m ) {
                    return 0;
                } else if( k > zero && k_abs >= n ) {
                    return 0;
                }
                return std::min( std::min( m, n ), k < zero
                        ? m - k_abs
                        : n - k_abs
                    );
            }

            template< class IteratorV >
            static MatrixType createIdentity_generic(
                const size_t m, const size_t n, const long k,
                const IteratorV V_iter, const IteratorV V_end
            ) {
#ifdef _DEBUG
                std::cout << "createIdentity_generic called with m = " << m << ", n = "
                    << n << ", k = " << k << " (non-void variant)\n";
#endif
                // some useful scalars
                constexpr const long s_zero = static_cast< long >( 0 );
                constexpr const size_t u_zero = static_cast< size_t >( 0 );
                const size_t diag_length = compute_diag_length( m, n, k );
#ifdef NDEBUG
                (void) V_end;
#endif
#ifdef _DEBUG
                std::cout << "Computed diag_length = " << diag_length << "\n";
#endif

                // declare matrix-to-be-returned
                MatrixType matrix( m, n, diag_length );

                // get row- and column-wise starting positions
                const RIT k_i_incr = static_cast< RIT >(
                    (k < s_zero) ? std::abs( k ) : u_zero );
                const CIT k_j_incr = static_cast< CIT >(
                    (k < s_zero) ? u_zero : std::abs( k ) );

                // translate it to a range so we can get iterators
                grb::utils::containers::Range< RIT > I( k_i_incr, diag_length + k_i_incr );
                grb::utils::containers::Range< CIT > J( k_j_incr, diag_length + k_j_incr );

                // construct the matrix from the given iterators
                const size_t s = getPID();
                const size_t P = getP();
                assert( static_cast< size_t >(std::distance( V_iter, V_end )) ==
                    static_cast< size_t >(std::distance( I.begin( s, P ), I.end( s, P ) )) );
                assert( static_cast< size_t >(std::distance( V_iter, V_end )) ==
                    static_cast< size_t >(std::distance( J.begin( s, P ), J.end( s, P ) )) );
                const RC rc = buildMatrixUnique(
                    matrix,
                    I.begin( s, P ), I.end( s, P ),
                    J.begin( s, P ), J.end( s, P ),
                    V_iter, V_end,
                    mode
                );

                if( rc != SUCCESS ) {
                    throw std::runtime_error(
                        "Error: createIdentity_generic failed: rc = " + grb::toString( rc )
                    );
                }
                return matrix;
            }


        public:

            static MatrixType empty(
                const size_t m, const size_t n
            ) {
                return MatrixType( m, n, 0 );
            }

            template< class ValueIterator >
            static MatrixType diag(
                const size_t m, const size_t n,
                const ValueIterator values, const ValueIterator valEnd,
                const long k = static_cast< long >( 0 )
            ) {
                // static sanity checks
                static_assert( mode == grb::SEQUENTIAL,
                    "matrices<>::diag is currently only supported with sequential I/O" );
                static_assert(
                    std::is_convertible<
                        typename std::iterator_traits< ValueIterator >::value_type,
                        D
                    >::value,
                    "grb::factory::diag called with value types that are not convertible to the "
                    "matrix nonzero type."
                );

                // check trivial dispatch
                if( m == 0 || n == 0 ) {
                    return empty( n, n );
                }

                // call generic implementation with given iterator
                return createIdentity_generic( m, n, k, values, valEnd );
            }

            static MatrixType eye(
                const size_t m,
                const size_t n,
                const D value = static_cast< D >( 1 ),
                const long k = static_cast< long >( 0 )
            ) {
                // check for possible trivial dispatch
                if( m == 0 || n == 0 ) {
                    return empty( m, n );
                }

                // we are in the non-pattern case, so we need a repeated iterator for the
                // values-- determine worst-case length (cheaper than computing the actual
                // diagonal length)
                const size_t diag_length = compute_diag_length( m, n, k );
                const grb::utils::containers::ConstantVector< D > V( value,
                    diag_length );

                // call generic implementation
                const size_t s = getPID();
                const size_t P = getP();
                return createIdentity_generic( m, n, k, V.cbegin( s, P ), V.cend( s, P ) );
            }

            static MatrixType eye( const size_t n ) {
                return eye( n, n );
            }

            template<
                typename Semiring = grb::Semiring<
                    grb::operators::add< D >, grb::operators::mul< D >,
                    grb::identities::zero, grb::identities::one
                >
            >
            static MatrixType identity(
                const size_t n,
                const Semiring &ring = Semiring()
            ) {
                const D identity_value = ring.template getOne< D >();
                return eye( n, n, identity_value, 0 );
            }

            static MatrixType full(
                const size_t m, const size_t n,
                const D value
            ) {
                if( m == 0 || n == 0 ) {
                    return empty( m, n );
                }

                const size_t nz = m * n;
                if( nz / m != n ) {
                    throw std::runtime_error( "Requested dense matrix overflows in number of "
                        "nonzeroes." );
                }

                const size_t s = getPID();
                const size_t P = getP();

                // Initialise rows indices container with a range from 0 to nrows,
                // each value repeated ncols times.
                grb::utils::containers::Range< RIT > I( 0, m, 1, n );

                // Initialise columns values container with a range from 0 to ncols
                // repeated nrows times. There are two ways of doing this:
                //  1) using InterleavedIterators, or
                //  2) using the iterator::Adapter.
                // We select here way #2, and disable way #1:
#if 0
                grb::utils::containers::InterleavedIteratorsVector<
                        typename grb::utils::containers::Range< CIT >::const_iterator
                    > J( m );
                for( size_t i = 0; i < m; ++i ) {
                    J.emplace_back( grb::utils::containers::Range< CIT >( 0, n ) );
                }
#endif
                const grb::utils::containers::Range< size_t > J_raw( 0, nz );
                const auto entryInd2colInd = [&m] (const CIT k) -> CIT {
                        return k / m;
                    };
                auto J_begin = grb::utils::iterators::make_adapter_iterator(
                    J_raw.cbegin( s, P ), J_raw.cend( s, P ), entryInd2colInd );
                const auto J_end = grb::utils::iterators::make_adapter_iterator(
                    J_raw.cend( s, P ), J_raw.cend( s, P ), entryInd2colInd );

                // Initialise values container with the given value.
                const size_t local_nz = std::distance( I.begin( s, P ), I.end( s, P ) );
#ifndef NDEBUG
                const size_t Isz = local_nz;
                const size_t Jsz = std::distance( J_begin, J_end );
                assert( Isz == Jsz );
#endif
                grb::utils::containers::ConstantVector< D > V( value, local_nz );
#ifndef NDEBUG
                const size_t Vsz = std::distance( V.begin(), V.end() );
                assert( Isz == Vsz );
#endif

                // allocate and build
                MatrixType matrix( m, n, nz );
                const RC rc = buildMatrixUnique(
                    matrix,
                    I.begin( s, P ), I.end( s, P ),
                    J_begin, J_end,
                    V.begin(), V.end(),
                    mode
                );

                if( rc != SUCCESS ) {
                    throw std::runtime_error(
                        "Error: factory::full<void> failed: rc = " + grb::toString( rc )
                    );
                }

                return matrix;
            }

            static MatrixType dense(
                const size_t m, const size_t n, const D value
            ) {
                return full( m, n, value );
            }

            template<
                typename Semiring = grb::Semiring<
                    grb::operators::add< D >, grb::operators::mul< D >,
                    grb::identities::zero, grb::identities::one
                >
            >
            static MatrixType dense(
                const MatrixType &A,
                const Semiring &ring = Semiring()
            ) {
                const size_t m = grb::nrows( A );
                const size_t n = grb::ncols( A );
                if( n == 0 || m == 0 ) {
                    return empty( m, n );
                }

                if( grb::nnz( A ) == 0 ) {
                    return zeros( m, n, ring );
                }

                const size_t nz = m * n;
                if( nz / m != n ) {
                    throw std::runtime_error( "Requested dense matrix overflows in number of "
                        "nonzeroes." );
                }

                const auto addMon = ring.getAdditiveMonoid();
                const D zero = ring.template getZero< D >();
                MatrixType matrix( m, n, nz );
                grb::RC rc = grb::set( matrix, zero );
                rc = rc ? rc : grb::foldl( matrix, A, addMon );

                if( rc != grb::SUCCESS) {
                    throw std::runtime_error( "Could not promote input matrix to a dense one: "
                        + grb::toString( rc ) );
                }

                return matrix;
            }

            template<
                typename Semiring = grb::Semiring<
                    grb::operators::add< D >, grb::operators::mul< D >,
                    grb::identities::zero, grb::identities::one
                >
            >
            static MatrixType zeros(
                const size_t m, const size_t n,
                const Semiring &ring = Semiring()
            ) {
                const D zero = ring.template getZero< D >();
                return full( m, n, static_cast< D >( zero ) );
            }

            template<
                typename Semiring = grb::Semiring<
                    grb::operators::add< D >, grb::operators::mul< D >,
                    grb::identities::zero, grb::identities::one
                >
            >
            static MatrixType ones(
                const size_t m, const size_t n,
                const Semiring &ring = Semiring()
            ) {
                const D one = ring.template getOne< D >();
                return full( m, n, static_cast< D >( one ) );
            }

    }; // end class matrices

    template<
        grb::IOMode mode, grb::Backend backend,
        typename RIT, typename CIT, typename NIT
    >
    class matrices< void, mode, backend, RIT, CIT, NIT > {

        private:

            typedef Matrix< void, backend, RIT, CIT, NIT > MatrixType;

            typedef matrices< int, mode, backend, RIT, CIT, NIT > BaseType;

            static MatrixType createIdentity_generic(
                const size_t m, const size_t n, const long k
            ) {
                // pattern matrix variant of the above
                constexpr const long s_zero = static_cast< long >( 0 );
                constexpr const size_t u_zero = static_cast< size_t >( 0 );
                const size_t diag_length = BaseType::compute_diag_length( m, n, k );
                MatrixType matrix( m, n, diag_length );
                const RIT k_i_incr = static_cast< RIT >(
                    (k < s_zero) ? std::abs( k ) : u_zero );
                const CIT k_j_incr = static_cast< CIT >(
                    (k < s_zero) ? u_zero : std::abs( k ) );
                grb::utils::containers::Range< RIT > I( k_i_incr, diag_length + k_i_incr );
                grb::utils::containers::Range< CIT > J( k_j_incr, diag_length + k_j_incr );
                const size_t s = BaseType::getPID();
                const size_t P = BaseType::getP();
                const RC rc = buildMatrixUnique(
                    matrix,
                    I.begin( s, P ), I.end( s, P ),
                    J.begin( s, P ), J.end( s, P ),
                    mode
                );
                if( rc != SUCCESS ) {
                    throw std::runtime_error(
                        "Error: createIdentity_generic<void> failed: rc = " + grb::toString( rc )
                    );
                }
                return matrix;
            }


        public:

            static MatrixType empty(
                const size_t m, const size_t n
            ) {
                return MatrixType( m, n, 0 );
            }

            MatrixType diag(
                const size_t m, const size_t n,
                const long k = static_cast< long >( 0 )
            ) {
                // check trivial dispatch
                if( m == 0 || n == 0 ) {
                    return empty( n, n );
                }

                // call generic implementation with given iterator
                return createIdentity_generic( m, n, k );
            }

            static MatrixType eye(
                const size_t m,
                const size_t n,
                const long k = static_cast< long >( 0 )
            ) {
                // check trivial case
                if( m == 0 || n == 0 ) {
                    return empty( m, n );
                }

                // dispatch to generic function
                return createIdentity_generic( m, n, k );
            }

            static MatrixType eye( const size_t n ) {
                return eye( n, n );
            }

            static MatrixType identity(
                const size_t n
            ) {
                return eye( n, n, 0 );
            }

            static MatrixType full( const size_t m, const size_t n ) {
                if( m == 0 || n == 0 ) {
                    return empty( m, n );
                }

                const size_t nz = m * n;
                if( nz / m != n ) {
                    throw std::runtime_error( "Requested dense matrix overflows in number of "
                        "nonzeroes." );
                }

                const size_t s = BaseType::getPID();
                const size_t P = BaseType::getP();

                // Initialise rows indices container with a range from 0 to nrows,
                // each value repeated ncols times.
                grb::utils::containers::Range< RIT > I( 0, m, 1, n );

                // Initialise columns values container with a range from 0 to ncols
                // repeated nrows times. As mentioned above, there are two ways to provide
                // iterators, we disable the first way (via InterleavedIterators) and enable
                // the iterators::adaptor way:
#if 0
                grb::utils::containers::InterleavedIteratorsVector<
                        typename grb::utils::containers::Range< CIT >::const_iterator
                    > J( m );
                for( size_t i = 0; i < m; ++i ) {
                    J.emplace_back( grb::utils::containers::Range< CIT >( 0, n ) );
                }
#endif
                const grb::utils::containers::Range< size_t > J_raw( 0, nz );
                const auto nonzeroInd2colInd = [&m] (const size_t k) -> size_t {
                        return k / m;
                    };
                auto J_begin = grb::utils::iterators::make_adapter_iterator(
                    J_raw.cbegin( s, P ), J_raw.cend( s, P ), nonzeroInd2colInd );
                auto J_end = grb::utils::iterators::make_adapter_iterator(
                    J_raw.cend( s, P ), J_raw.cend( s, P ), nonzeroInd2colInd );

                assert( std::distance( I.begin( s, P ), I.end( s, P ) ) ==
                    std::distance( J_begin, J_end ) );

                // construct and populate matrix
                MatrixType matrix( m, n, nz );
                const RC rc = buildMatrixUnique(
                    matrix,
                    I.begin( s, P ), I.end( s, P ),
                    J_begin, J_end,
                    mode
                );

                if( rc != SUCCESS ) {
                    throw std::runtime_error(
                        "Error: factory::full<void> failed: rc = " + grb::toString( rc )
                    );
                }

                return matrix;
            }

            static MatrixType dense( const size_t m, const size_t n ) {
                return full( m, n );
            }

            static MatrixType dense( const MatrixType &A ) {
                (void) A;
                return full( grb::nrows( A ), grb::ncols( A ) );
            }

            static MatrixType zeros( const size_t m, const size_t n ) {
                return full( m, n );
            }

            static MatrixType ones( const size_t m, const size_t n ) {
                return full( m, n );
            }

    }; // end class matrices (pattern specialisation)

} // namespace grb::algorithms

#endif // end _H_GRB_MATRIX_FACTORY