Expression.h

Go to the documentation of this file.

/*
 * Expression Template Matrix Library
 *
 * Copyright (C) 2004 - 2006 Ricky Lung <mtlung@users.sourceforge.net>
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */


#ifdef _MSC_VER
#   pragma once
#   pragma warning( push )
#   pragma warning( disable : 4100 4127 4512 4786 )
#endif  // _MSC_VER

#ifndef _EXMAT_EXPRESSION_H
#define _EXMAT_EXPRESSION_H

#include "PlatformSpec.h"
#include "TypeTraits.h"
#include "Container/Container.h"
#include "loki/TypeTraits.h"

namespace exmat {



template<class Rep>
struct ExpMat : public Rep {
    enum {
        IsLinear = Rep::IsLinear,
        IsLinearRecursive = Rep::IsLinearRecursive
    };

    typedef ExpMat<Rep>                 self_type;
    typedef typename Rep::rep_type      rep_type;
    // Should use Rep::super_type:
    //typedef typename Rep::super_type  super_type;


    template<typename T1> EXMAT_INLINE2
    ExpMat(T1& t1) : Rep(t1) {}
    template<typename T1, typename T2> EXMAT_INLINE2
    ExpMat(T1& t1, T2& t2) : Rep(t1, t2) {}
    template<typename T1, typename T2, typename T3> EXMAT_INLINE2
    ExpMat(T1& t1, T2& t2, T3& t3) : Rep(t1, t2, t3) {}
    template<typename T1, typename T2, typename T3, typename T4> EXMAT_INLINE2
    ExpMat(T1& t1, T2& t2, T3& t3, T4& t4) : Rep(t1, t2, t3, t4) {}
    template<typename T1, typename T2, typename T3, typename T4, typename T5> EXMAT_INLINE2
    ExpMat(T1& t1, T2& t2, T3& t3, T4& t4, T5& t5) : Rep(t1, t2, t3, t4, t5) {}
#if EXMAT_GNUC || EXMAT_VC2003 || EXMAT_VC2005
    // VC7.1, VC8 and GCC need this to allow force inlining
    EXMAT_INLINE2 ExpMat(const self_type& self) : Rep(self) {}
#endif
#if !EXMAT_VC6 && !EXMAT_VC7
    template<typename T1, typename T2> EXMAT_INLINE2
    ExpMat(const T1& t1, ScalarCon<T2> t2) : Rep(t1, t2) {}
    template<typename T1, typename T2> EXMAT_INLINE2
    ExpMat(ScalarCon<T1> t1, const T2& t2) : Rep(t1, t2) {}
#endif  // HAVE_PARTIAL_SPECIALIZATION


    template<typename T1> EXMAT_INLINE2
    ExpMat& operator=(const T1& t1) {
        (*static_cast<Rep*>(this)) = t1;
        return *this;
    }
};



template<class EXP>
struct ExpRow {
    typedef typename EXP::value_type value_type;
    typedef typename EXP::index_type index_type;

    EXMAT_INLINE2 ExpRow(const EXP& exp_, index_type i)
        : exp(exp_), row(i) {}
    EXMAT_INLINE2 value_type operator[](const index_type i) const {
        return exp.getAt(row, i);
    }
    const EXP& exp;
    index_type row;
};


template<class OperandType>
struct OperandTypeChooser {
    typedef typename
    IF<
        HasTagOf<OperandType,ScalarTag>::RET ||
        HasTagOf<OperandType,UnaExpTag>::RET ||
        HasTagOf<OperandType,BinExpTag>::RET,
        const OperandType,
        const OperandType&
    >::RET RET;
};


template<class OpT1, class OpT2, typename VT>
struct BinViewValueTypeChooser {
    enum {
        CanLRef = !(
            HasTagOf<OpT1,ScalarTag>::RET ||
            HasTagOf<OpT1,UnaExpTag>::RET ||
            HasTagOf<OpT1,BinExpTag>::RET),
        CanRRef = !(
            HasTagOf<OpT2,ScalarTag>::RET ||
            HasTagOf<OpT2,UnaExpTag>::RET ||
            HasTagOf<OpT2,BinExpTag>::RET)
    };
    typedef typename
    IF<
        // Only both left and right can use reference then the final type can be reference
        CanLRef && CanRRef,
        VT&,
        VT
    >::RET RET;
};


template<typename T, int ROWS, int COLS, class Left, class Right>
struct GetBinExpTmpType {
    typedef Mat<DenseMatCon<T,ROWS,COLS>, EmptyErrorChecker> RET;
};



template<class Rep>
struct UnaExp : public UnaExpTag {
    typedef UnaExp<Rep> self_type;
    // -> row_type should be override by child expression
    typedef NullTag                         row_type;
    typedef NullTag                         const_row_type;
    typedef typename Rep::col_type          col_type;
    typedef typename Rep::value_type        value_type;
    typedef NullTag                         iterator;
    typedef NullTag                         const_iterator;
    // -> Expression should return a new value, so reference will be treated as value
    typedef value_type                      reference;
    typedef const value_type                const_reference;
    typedef typename Rep::index_type        index_type;
    typedef typename Rep::param_type        param_type;
    typedef typename Rep::const_param_type  const_param_type;

    typedef typename OperandTypeChooser<Rep>::RET OpType;

    EXMAT_INLINE2 UnaExp(const self_type& self) : rep(self.rep) {}
    EXMAT_INLINE2 UnaExp(const Rep& r) : rep(r) {}

    // Constant reference/value to the operand
    OpType rep;
};



template<class Left, class Right>
struct BinExp : public BinExpTag {
    typedef BinExp<Left, Right>             self_type;
    typedef NullTag                         rep_type;
    // -> row_type should be override by child expression
    typedef NullTag                         row_type;
    typedef NullTag                         const_row_type;
    typedef typename Left::col_type         col_type;
    typedef HEleType(typename Left::value_type, typename Right::value_type) value_type;
    typedef NullTag                         iterator;
    typedef NullTag                         const_iterator;
    // -> Expression should return a new value, so reference will be treated as value
    typedef value_type                      reference;
    typedef const value_type                const_reference;
    // --> Assume both operands' index_type are the same
    typedef typename Left::index_type       index_type;
    typedef typename Loki::TypeTraits<value_type>::ParameterType param_type;
    typedef typename Loki::TypeTraits<const value_type>::ParameterType const_param_type;

    typedef Left left_type;
    typedef Right right_type;

    EXMAT_INLINE2 BinExp(const self_type& self)
        : leftOp(self.leftOp), rightOp(self.rightOp) {}
    EXMAT_INLINE2 BinExp(const Left& l, const Right& r)
        : leftOp(l), rightOp(r) {}

    typedef typename OperandTypeChooser<Left>::RET lOpType;
    typedef typename OperandTypeChooser<Right>::RET rOpType;

    // Constant reference/value to the operands
    lOpType leftOp;
    rOpType rightOp;
};



template<class Left, class Right>
struct LinearBinExpSize {
    enum {
        ROWS =
            (int(Left::ROWS) == 0 && int(Right::ROWS) == 0) ?
            0 : (
            int(Left::ROWS) == 0 ? int(Right::ROWS) :
                int(Left::ROWS)
            ),
        COLS =
            (int(Left::COLS) == 0 && int(Right::COLS) == 0) ?
            0 : (
            int(Left::COLS) == 0 ? int(Right::COLS) :
                int(Left::COLS)
            )
    };
};



template<class Left, class Right>
struct LinBinExp : public BinExp<Left, Right>
{
    typedef BinExp<Left, Right>             super_type;
    typedef LinBinExp<Left, Right>          self_type;
    typedef typename super_type::value_type value_type;

    enum {
        ROWS = LinearBinExpSize<Left, Right>::ROWS,
        COLS = LinearBinExpSize<Left, Right>::COLS,
        IsLinear = true,
        IsLinearRecursive = IsLinear && Left::IsLinearRecursive && Right::IsLinearRecursive
    };

    using super_type::leftOp;
    using super_type::rightOp;

    typedef typename GetBinExpTmpType<value_type,ROWS,COLS,Left,Right>::RET temp_type;

    enum {
        TNOP = (ROWS != 0 && COLS != 0) && (Left::TNOP >= 0 && Left::TNOP >= 0) ?
            Left::TNOP + Right::TNOP + ROWS*COLS : -1,
        ENOP = (ROWS != 0 && COLS != 0) && (Left::TNOP >= 0 && Left::TNOP >= 0) ?
            Left::ENOP + Right::ENOP + 1 : -1,
        EXPLevel = Left::EXPLevel + Right::EXPLevel + 2
    };

    EXMAT_INLINE2 LinBinExp(const self_type& self)
        : super_type(self.leftOp, self.rightOp) {}

    EXMAT_INLINE2 LinBinExp(const Left& l, const Right& r)
        : super_type(l, r)
    {
        // Do Incompatibility check here
#if EXMAT_STATIC_CHECK
        // Static check
        StaticAssert1(Left::ROWS == Right::ROWS, RowNotMatch);
        StaticAssert1(Left::COLS == Right::COLS, ColNotMatch);
#endif  // EXMAT_STATIC_CHECK

        // Dynamic check
        if(!EXMAT_STATIC_CHECK || ROWS == 0 || COLS == 0) {
            // Check for scalar, which have zero dimension
            if(!(   leftOp.rows() == 0 || rightOp.rows() == 0 ||
                    leftOp.cols() == 0 || rightOp.cols() == 0))
            DefaultErrorChecker::AssertCompat(leftOp.rows() == rightOp.rows() && leftOp.cols() == rightOp.cols());
        }
    }

    EXMAT_INLINE2 size_t size() const {
        return rows() * cols();
    }
    EXMAT_INLINE2 size_t rows() const {
        return ROWS != 0 ? ROWS : (
            leftOp.rows() > rightOp.rows() ?
            leftOp.rows() : rightOp.rows());
    }
    EXMAT_INLINE2 size_t cols() const {
        return COLS != 0 ? COLS : (
            leftOp.cols() > rightOp.cols() ?
            leftOp.cols() : rightOp.cols());
    }
};



template<class Left, class Right>
struct MulBinExpSize {
    enum {
        ROWS = Left::ROWS,
        COLS = Right::COLS,
        MulDotSize =
            (int(Left::COLS) == 0 && int(Right::ROWS) == 0) ?
            0 : (
            int(Left::COLS) == 0 ? int(Right::ROWS) :
                int(Left::COLS)
            )
    };
};



template<class Left, class Right>
struct MulBinExp : public BinExp<Left, Right>
{
    typedef BinExp<Left, Right>             super_type;
    typedef MulBinExp<Left, Right>          self_type;
    typedef self_type                       rep_type;
    typedef typename super_type::value_type value_type;

    enum {
        ROWS = MulBinExpSize<Left, Right>::ROWS,
        COLS = MulBinExpSize<Left, Right>::COLS,
        MulDotSize = MulBinExpSize<Left, Right>::MulDotSize,
        IsLinear = false,
        IsLinearRecursive = IsLinear
    };

    using super_type::leftOp;
    using super_type::rightOp;

    typedef typename GetBinExpTmpType<value_type,ROWS,COLS,Left,Right>::RET temp_type;

    enum {
        TNOP = (ROWS != 0 && COLS != 0) && (Left::TNOP >= 0 && Left::TNOP >= 0) ?
            Left::TNOP + Right::TNOP + (ROWS*COLS)*(2*MulDotSize-1) : -1,
        ENOP = (ROWS != 0 && COLS != 0) && (Left::ENOP >= 0 && Left::ENOP >= 0) ?
            (Left::ENOP + Right::ENOP + 2)*MulDotSize-1 : -1,
        EXPLevel = Left::EXPLevel + Right::EXPLevel + 2
    };

    EXMAT_INLINE2 MulBinExp(const self_type& self)
        : super_type(self) {}
    EXMAT_INLINE2 MulBinExp(const Left& l, const Right& r)
        : super_type(l, r)
    {
        // Do compatability check here
        // The check have been delaied to the child class, to allow customs
        // MulExp (eg. MulExp in G3D) without throwing exception. 
    }

    EXMAT_INLINE2 size_t size() const { return rows() * cols(); }
    EXMAT_INLINE2 size_t rows() const { return leftOp.rows();   }
    EXMAT_INLINE2 size_t cols() const { return rightOp.cols();  }
};



template<class Left, class Right, class SpecTag=NullTag>
class AddExp :
public LinBinExp<Left, Right>,
public AddExpTag
{
    typedef LinBinExp<Left, Right>          super_type;
    typedef AddExp<Left, Right, SpecTag>    self_type;
public:
    typedef self_type                       rep_type;
    typedef const ExpRow<self_type>         const_row_type;
    typedef typename super_type::value_type value_type;
    typedef typename super_type::index_type index_type;
    typedef AddExpTag exp_tag;  // Tag indicating the type of expression

    using super_type::leftOp;
    using super_type::rightOp;

    EXMAT_INLINE2 AddExp(const self_type& self)
        : super_type(self) {}
    EXMAT_INLINE2 AddExp(const Left& l, const Right& r)
        : super_type(l, r) {}

    /************************************************/
    // Index function
    EXMAT_INLINE2 const_row_type operator[] (const index_type i) const {
        return const_row_type(*this, i);
    }
    EXMAT_INLINE2 value_type getAt(const index_type i) const {
        return (value_type)(leftOp.getAt(i)) + (value_type)(rightOp.getAt(i));
    }
    EXMAT_INLINE2 value_type getAt(const index_type i, const index_type j) const {
        return (value_type)(leftOp.getAt(i,j)) + (value_type)(rightOp.getAt(i,j));
    }
};


template<class Left, class Right, class SpecTag=NullTag>
class SubExp :
public LinBinExp<Left, Right>,
public SubExpTag
{
    typedef LinBinExp<Left, Right>          super_type;
    typedef SubExp<Left, Right, SpecTag>    self_type;
public:
    typedef self_type                       rep_type;
    typedef const ExpRow<self_type>         const_row_type;
    typedef typename super_type::value_type value_type;
    typedef typename super_type::index_type index_type;
    typedef SubExpTag exp_tag;  // Tag indicating the type of expression

    using super_type::leftOp;
    using super_type::rightOp;

    EXMAT_INLINE2 SubExp(const self_type& self)
        : super_type(self) {}
    EXMAT_INLINE2 SubExp(const Left& l, const Right& r)
        : super_type(l, r) {}

    /************************************************/
    // Index function
    EXMAT_INLINE2 const_row_type operator[] (const index_type i) const {
        return const_row_type(*this, i);
    }
    EXMAT_INLINE2 value_type getAt(const index_type i) const {
        return (value_type)(leftOp.getAt(i)) - (value_type)(rightOp.getAt(i));
    }
    EXMAT_INLINE2 value_type getAt(const index_type i, const index_type j) const {
        return (value_type)(leftOp.getAt(i,j)) - (value_type)(rightOp.getAt(i,j));
    }
};


template<class Left, class Right, class SpecTag=NullTag>
class MulExp :
public MulBinExp<Left, Right>,
public MulExpTag
{
    typedef MulBinExp<Left, Right>          super_type;
    typedef MulExp<Left, Right, SpecTag>    self_type;
public:
    typedef self_type                       rep_type;
    typedef const ExpRow<self_type>         const_row_type;
    typedef typename super_type::value_type value_type;
    typedef typename super_type::index_type index_type;
    typedef MulExpTag                       exp_tag;

    using super_type::leftOp;
    using super_type::rightOp;

    enum {
        TNOP = super_type::TNOP,
        ENOP = super_type::ENOP,
        MulDotSize = super_type::MulDotSize,
        EXPLevel = super_type::EXPLevel
    };

    EXMAT_INLINE2 MulExp(const self_type& self)
        : super_type(self) {}
    EXMAT_INLINE2 MulExp(const Left& l, const Right& r)
        : super_type(l, r)
    {
        // Do compatability check here
#if EXMAT_STATIC_CHECK
        // Static check
        StaticAssert(Left::ROWS == Right::COLS, LeftRowNotMatchRightCol);
        StaticAssert(Left::COLS == Right::ROWS, LeftColNotMatchRightRow);
#endif  // EXMAT_STATIC_CHECK
        // Dynamic check
        if(!EXMAT_STATIC_CHECK || Left::COLS == 0 || Right::COLS == 0) {
            // Check for scalar, which have zero dimension
            if(!(   leftOp.rows() == 0 || rightOp.rows() == 0 ||
                    leftOp.cols() == 0 || rightOp.cols() == 0))
            DefaultErrorChecker::AssertCompat(leftOp.cols() == rightOp.rows());
        }
    }

    /************************************************/
    // Index function
    EXMAT_INLINE2 const_row_type operator[] (const index_type i) const {
        return const_row_type(*this, i);
    }

protected:
    // Unrolling for inner_product
#if HAVE_PARTIAL_SPECIALIZATION
    // Loop full unrolling template for the dot product
    template<unsigned int I, class A, class B, typename value_type, typename index_type>
    struct MulDot_FullUnRoll {
        enum { go = (I-1) != 0 };
        static EXMAT_INLINE2 value_type         //    A's row               B's col
        exec(const A& a, const B& b, const index_type& ar, const index_type& bc) {
            return a.getAt(ar, I-1) * b.getAt(I-1, bc) +
                MulDot_FullUnRoll< go?(I-1):0, A, B, value_type, index_type >::exec(a, b, ar, bc);
        }
    };
    // Partial specialization to end the recursion
    template<class A, class B, typename value_type, typename index_type>
    struct MulDot_FullUnRoll<0, A, B, value_type, index_type> {
        static EXMAT_INLINE2 value_type
        exec(const A& a, const B& b, const index_type& ar, const index_type& bc) {
            return value_type(0);
        }
    };
#else       // Workaround for MSVC6
    template<unsigned int SIZE, class A, class B, typename value_type, typename index_type>
    struct MulDot_FullUnRoll {
        template<unsigned int I> struct inner {
            enum { go = (I-1) != 0 };
            static EXMAT_INLINE2 value_type
            exec(const A& a, const B& b, const index_type& ar, const index_type& bc) {
                return a.getAt(ar, I-1) * b.getAt(I-1, bc) +
                    inner< go?(I-1):0 >::exec(a, b, ar, bc);
            }
        };
        template<> struct inner<0> {
            static EXMAT_INLINE2 value_type
            exec(const A& a, const B& b, const index_type& ar, const index_type& bc) {
                return value_type(0);
            }
        };
        static EXMAT_INLINE2 value_type
        exec(const A& a, const B& b, const index_type& ar, const index_type& bc) {
            return inner<SIZE>::exec(a, b, ar, bc);
        }
    };
#endif  // HAVE_PARTIAL_SPECIALIZATION

    /************************************************/
    //  Overloaded function to choose the appropriate unroller
    EXMAT_INLINE2 value_type
    MulDot(const index_type& i, const index_type& j, Int2Type<true>) const {
        return MulDot_FullUnRoll<MulDotSize, Left, Right, value_type, index_type>::
            exec(leftOp, rightOp, i, j);
    }
    EXMAT_INLINE2 value_type
    MulDot(const index_type& i, const index_type& j, Int2Type<false>) const {
        value_type sum = value_type(0);
        for(index_type k=leftOp.cols(); k--; ) {
            sum += (value_type)(leftOp.getAt(i,k)) * (value_type)(rightOp.getAt(k,j));
        }
        return sum;
    }

public:
    EXMAT_INLINE2 value_type getAt(const index_type i) const {
        return getAt(i/this->cols(), i%this->cols());
    }
    EXMAT_INLINE2 value_type getAt(const index_type i, const index_type j) const {
#if ENABLE_UNROLLING
        enum {
            CanUnroll = (MulDotSize > 0) && (ENOP < MulDotSize * MAX_MULDOTUROLL_ENOP)
        };
        return MulDot(i, j, Int2Type<CanUnroll>());
#else
        return MulDot(i, j, Int2Type<false>());
#endif
    }
};


template<class Left, class Right, class SpecTag=NullTag>
class ScalarMulExp :
public LinBinExp<Left, Right>,
public ScalarMulExpTag
{
    typedef LinBinExp<Left, Right>          super_type;
    typedef ScalarMulExp<Left, Right, SpecTag>  self_type;
public:
    typedef self_type                       rep_type;
    typedef const ExpRow<self_type>         const_row_type;
    typedef typename super_type::value_type value_type;
    typedef typename super_type::index_type index_type;
    typedef ScalarMulExpTag exp_tag;    // Tag indicating the type of expression

    using super_type::leftOp;
    using super_type::rightOp;

    EXMAT_INLINE2 ScalarMulExp(const self_type& self)
        : super_type(self) {}
    EXMAT_INLINE2 ScalarMulExp(const Left& l, const Right& r)
        : super_type(l, r) {}

    /************************************************/
    // Index function
    EXMAT_INLINE2 const_row_type operator[] (const index_type i) const {
        return const_row_type(*this, i);
    }
    EXMAT_INLINE2 value_type getAt(const index_type i) const {
        return (value_type)(leftOp.getAt(i)) * (value_type)(rightOp.getAt(i));
    }
    EXMAT_INLINE2 value_type getAt(const index_type i, const index_type j) const {
        return (value_type)(leftOp.getAt(i,j)) * (value_type)(rightOp.getAt(i,j));
    }
};


template<class Left, class Right, class SpecTag=NullTag>
class ScalarDivExp :
public LinBinExp<Left, Right>,
public ScalarDivExpTag
{
    typedef LinBinExp<Left, Right>          super_type;
    typedef ScalarDivExp<Left, Right, SpecTag>  self_type;
public:
    typedef self_type                       rep_type;
    typedef const ExpRow<self_type>         const_row_type;
    typedef typename super_type::value_type value_type;
    typedef typename super_type::index_type index_type;
    typedef ScalarDivExpTag exp_tag;    // Tag indicating the type of expression

    using super_type::leftOp;
    using super_type::rightOp;

    EXMAT_INLINE2 ScalarDivExp(const self_type& self)
        : super_type(self) {}
    EXMAT_INLINE2 ScalarDivExp(const Left& l, const Right& r)
        : super_type(l, r) {}

    /************************************************/
    // Index function
    EXMAT_INLINE2 const_row_type operator[] (const index_type i) const {
        return const_row_type(*this, i);
    }
    EXMAT_INLINE2 value_type getAt(const index_type i) const {
        return (value_type)(leftOp.getAt(i)) / (value_type)(rightOp.getAt(i));
    }
    EXMAT_INLINE2 value_type getAt(const index_type i, const index_type j) const {
        return (value_type)(leftOp.getAt(i,j)) / (value_type)(rightOp.getAt(i,j));
    }
};


}   // namespace exmat

#if _MSC_VER
#   pragma warning( pop )
#endif  // EXMAT_VC

#endif  // _EXMAT_EXPRESSION_H

---