dlib C++ Library - matrix_assign.h

    // Copyright (C) 2008  Davis E. King (davisking@users.sourceforge.net)
    // License: Boost Software License   See LICENSE.txt for the full license.
    #ifndef DLIB_MATRIx_ASSIGn_
    #define DLIB_MATRIx_ASSIGn_

    #include "../geometry.h"
    #include "matrix.h"
    #include "matrix_utilities.h"
    #include "../enable_if.h"

    namespace dlib
    {

        namespace ma
        {

        // ------------------------------------------------------------------------------------

            template <
                typename EXP
                >
            const matrix_exp<EXP> make_exp (
                const EXP& exp
            )
            /*!
                The only point of this function is to make it easy to cause the overloads
                of matrix_assign to not trigger for a matrix expression.
            !*/
            {
                return matrix_exp<EXP>(exp);
            }

        // ------------------------------------------------------------------------------------

            template < typename EXP, typename enable = void >
            struct matrix_is_vector { static const bool value = false; };
            template < typename EXP >
            struct matrix_is_vector<EXP, typename enable_if_c<EXP::NR==1 || EXP::NC==1>::type > { static const bool value = true; };

            template < typename EXP, typename enable = void >
            struct is_small_matrix { static const bool value = false; };
            template < typename EXP >
            struct is_small_matrix<EXP, typename enable_if_c<EXP::NR>=1 && EXP::NC>=1 &&
            EXP::NR<=100 && EXP::NC<=100>::type > { static const bool value = true; };

        }

    // ----------------------------------------------------------------------------------------

        template <
            typename matrix_dest_type,
            typename EXP1,
            typename EXP2,
            unsigned long count
            >
        inline typename disable_if_c<ma::matrix_is_vector<EXP1>::value || ma::matrix_is_vector<EXP2>::value ||
                                     ma::is_small_matrix<EXP1>::value || ma::is_small_matrix<EXP2>::value >::type matrix_assign (
            matrix_dest_type& dest,
            const matrix_exp<matrix_multiply_exp<EXP1,EXP2,count> >& src,
            const long row_offset = 0,
            const long col_offset = 0
        )
        /*!
            This overload catches assignments like:
                dest = lhs*rhs
                where lhs and rhs are both not vectors
        !*/
        {
            using namespace ma;
            const matrix_exp<EXP1> lhs(src.ref().lhs);
            const matrix_exp<EXP2> rhs(src.ref().rhs);
            const long bs = 100;

            // if the matrices are small enough then just use the simple multiply algorithm
            if (lhs.nc() <= 2 || rhs.nc() <= 2 || lhs.nr() <= 2 || rhs.nr() <= 2 || (lhs.size() <= bs*10 && rhs.size() <= bs*10) )
            {
                for (long r = 0; r < src.nr(); ++r)
                {
                    for (long c = 0; c < src.nc(); ++c)
                    {
                        dest(r+row_offset,c+col_offset) = make_exp(lhs*rhs)(r,c);
                    }
                }
            }
            else
            {
                // if the lhs and rhs matrices are big enough we should use a cache friendly
                // algorithm that computes the matrix multiply in blocks.

                // Loop over all the blocks in the lhs matrix
                for (long r = 0; r < lhs.nr(); r+=bs)
                {
                    for (long c = 0; c < lhs.nc(); c+=bs)
                    {
                        // make a rect for the block from lhs
                        rectangle lhs_block(c, r, c+bs-1, r+bs-1);
                        lhs_block = lhs_block.intersect(get_rect(lhs));

                        // now loop over all the rhs blocks we have to multiply with the current lhs block
                        for (long i = 0; i < rhs.nc(); i += bs)
                        {
                            // make a rect for the block from rhs
                            rectangle rhs_block(i, c, i+bs-1, c+bs-1);
                            rhs_block = rhs_block.intersect(get_rect(rhs));

                            // make a target rect in res
                            rectangle res_block(rhs_block.left(),lhs_block.top(), rhs_block.right(), lhs_block.bottom());
                            res_block = translate_rect(res_block,col_offset, row_offset);
                            if (c != 0)
                                set_subm(dest, res_block) = subm(dest,res_block) + subm(lhs,lhs_block)*subm(rhs, rhs_block);
                            else
                                set_subm(dest, res_block) = make_exp(subm(lhs,lhs_block)*subm(rhs, rhs_block));
                        }
                    }
                }
            }

        }

    // ----------------------------------------------------------------------------------------

    }

    #endif // DLIB_MATRIx_ASSIGn_