Mat.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
//  Disable the annoying warning of VC6
#   if _MSC_VER < 1300
#       pragma warning( disable : 4786 )
#   endif   // VC6
#   pragma once
#   pragma warning( push )
#   pragma warning( disable : 4100 4127 4239 4244 4512 4786 )
#endif  // _MSC_VER


#ifndef _EXMAT_MAT_H
#define _EXMAT_MAT_H

#include <iterator>
#include <sstream>
#include <iostream>
#include <algorithm>
#include "PlatformSpec.h"
#include "Config.h"
#include "ErrorCheck.h"
#include "TypeTraits.h"
#include "Commainit.h"
#include "Assignment.h"
#include "SIMD/SIMD.h"

// Forward declaration of Mat and Vec
namespace exmat {
    template<class Rep, class ErrChecker_> class Mat;
    template<class Rep, class ErrChecker_> class Vec;
}   // namespace exmat

/************************************************/
namespace exmat {

namespace PNS {
    template<class Mat, typename IdxT> EXMAT_INLINE2 static
    void IOResize(Mat& m, IdxT r, IdxT c, Int2Type<true>, Int2Type<true>) {
        m.resize(r, c);
    }
    template<class Mat, typename IdxT> EXMAT_INLINE2 static
    void IOResize(Mat& m, IdxT r, IdxT c, Int2Type<true>, Int2Type<false>) {
        m.resizeRow(r);
    }
    template<class Mat, typename IdxT> EXMAT_INLINE2 static
    void IOResize(Mat& m, IdxT r, IdxT c, Int2Type<false>, Int2Type<true>) {
        m.resizeCol(c);
    }
    template<class Mat, typename IdxT> EXMAT_INLINE2 static
    void IOResize(Mat& m, IdxT r, IdxT c, Int2Type<false>, Int2Type<false>) {
    }
};};

template<class Rep, class ErrChecker> EXMAT_INLINE0
std::istream& operator>>(std::istream &is, exmat::Mat<Rep, ErrChecker>& mat) {
    typedef exmat::Mat<Rep, ErrChecker> mat_type;
    typedef typename mat_type::index_type index_type;
    index_type rs=0, cs=0;
    // Read the number of rows
    if(mat_type::IsDR) is >> rs;
    // Read the number of columns
    if(mat_type::IsDC) is >> cs;
    // Resize if needed
    exmat::PNS::IOResize(mat, rs, cs,
        exmat::Int2Type<mat_type::IsDR>(), exmat::Int2Type<mat_type::IsDC>());
    for(index_type i=0; i<mat.rows(); ++i) {
        for(index_type j=0; j<mat.cols(); ++j)
            is >> mat.setAt(i, j);
    }
    return is;
}

template<class Rep, class ErrChecker> EXMAT_INLINE0
std::ostream& operator<<(std::ostream &os, const exmat::Mat<Rep, ErrChecker>& mat) {
    typedef typename exmat::Mat<Rep, ErrChecker>::index_type index_type;
    for(index_type i=0; i<mat.rows(); ++i) {
        for(index_type j=0; j<mat.cols(); ++j)
            os << mat.getAt(i, j) << " ";
        os << std::endl;
    }
    return os;
}

//template<class Rep, class ErrChecker> EXMAT_INLINE0
//std::ostream& operator<<(std::ostream &os, const exmat::Vec<Rep, ErrChecker>& vec) {
//  typedef typename exmat::Vec<Rep, ErrChecker>::index_type index_type;
//  for(index_type i=0; i<vec.size(); ++i)
//      os << vec.getAt(i) << " ";
//
//  return os;
//}


/************************************************/
// Functional operators

// VC6 and 7 have problems on using ExpMat<>::rep_type
#if EXMAT_VC6 || EXMAT_VC7

//  Macro to create operators for 2 Mat
#define Mak_Op(LT, RT, OT, op, OpName)  \
template<class L, class R, class LC, class RC> EXMAT_INLINE2    \
    const exmat::OT<exmat::ExpMat<exmat::OpName<L::rep_type, R::rep_type> >, exmat::ResultantChecker<LC,RC>::RET>   \
op (const exmat::LT<L,LC>& A, const exmat::RT<R,RC>& B) {   \
    using namespace exmat;  \
    typedef ExpMat<OpName<L::rep_type, R::rep_type> > RETType;  \
    return RETType(A, B);   \
}
//  Macro to create operators for Mat and scalar
#define Mak_ScalarOp1(OT, op, OpName)   \
template<class L, class EC> EXMAT_INLINE2   \
    const exmat::OT<exmat::ExpMat<exmat::OpName<L::rep_type,exmat::ScalarCon<L::value_type> > >, EC>    \
op (const exmat::OT<L,EC>& A, L::const_param_type scalar) { \
    using namespace exmat;  \
    typedef ExpMat<OpName<L::rep_type,ScalarCon<L::value_type> > > RETType; \
    return RETType(A, ScalarCon<L::value_type>(scalar));    \
}
//  Macro to create operators for scalar and Mat
#define Mak_ScalarOp2(OT, op, OpName)   \
template<class L, class EC> EXMAT_INLINE2   \
    const exmat::OT<exmat::ExpMat<exmat::OpName<exmat::ScalarCon<L::value_type>, L::rep_type> >, EC>    \
op (typename L::const_param_type scalar, const exmat::OT<L,EC>& A) {    \
    using namespace exmat;  \
    typedef ExpMat<OpName<ScalarCon<L::value_type>, L::rep_type> > RETType; \
    return RETType(ScalarCon<L::value_type>(scalar), A);    \
}

#else

#define Mak_Op(LT, RT, OT, op, OpName)  \
template<class L, class R, class LC, class RC> EXMAT_INLINE2 static \
    const exmat::OT<exmat::ExpMat<typename exmat::OpName<typename L::rep_type, typename R::rep_type>::rep_type>, typename exmat::ResultantChecker<LC,RC>::RET>  \
op (const exmat::LT<L,LC>& A, const exmat::RT<R,RC>& B) {   \
    using namespace exmat;  \
    typedef ExpMat<TTYPENAME OpName<TTYPENAME L::rep_type,TTYPENAME R::rep_type>::rep_type> RETType;    \
    return RETType(A, B);   \
}

#define Mak_ScalarOp1(OT, op, OpName)   \
template<class L, class EC> EXMAT_INLINE2 static    \
    const exmat::OT<exmat::ExpMat<typename exmat::OpName<typename L::rep_type,exmat::ScalarCon<typename L::value_type> >::rep_type>, EC>    \
op (const exmat::OT<L,EC>& A, typename L::const_param_type scalar) {    \
    using namespace exmat;  \
    typedef ExpMat<TTYPENAME OpName<TTYPENAME L::rep_type,ScalarCon<TTYPENAME L::value_type> >::rep_type> RETType;  \
    return RETType(A, ScalarCon<TTYPENAME L::value_type>(scalar));  \
}

#define Mak_ScalarOp2(OT, op, OpName)   \
template<class L, class EC> EXMAT_INLINE2 static    \
    const exmat::OT<exmat::ExpMat<typename exmat::OpName<exmat::ScalarCon<typename L::value_type>, typename L::rep_type>::rep_type>, EC>    \
op (typename L::const_param_type scalar, const exmat::OT<L,EC>& A) {    \
    using namespace exmat;  \
    typedef ExpMat<TTYPENAME OpName<ScalarCon<TTYPENAME L::value_type>, TTYPENAME L::rep_type>::rep_type> RETType;  \
    return RETType(ScalarCon<TTYPENAME L::value_type>(scalar), A);  \
}

#endif  // EXMAT_VC6 || EXMAT_VC7



Mak_Op(Mat, Mat, Mat, operator+, AddExp);
Mak_Op(Vec, Vec, Vec, operator+, AddExp);

// Addition of Mat with scalar
Mak_ScalarOp1(Mat, operator+, AddExp);
Mak_ScalarOp2(Mat, operator+, AddExp);

// Addition of Vec with scalar
Mak_ScalarOp1(Vec, operator+, AddExp);
Mak_ScalarOp2(Vec, operator+, AddExp);



Mak_Op(Mat, Mat, Mat, operator-, SubExp);
Mak_Op(Vec, Vec, Vec, operator-, SubExp);

// Subtraction of Mat with scalar
Mak_ScalarOp1(Mat, operator-, SubExp);
Mak_ScalarOp2(Mat, operator-, SubExp);

// Subtraction of Vec with scalar
Mak_ScalarOp1(Vec, operator-, SubExp);
Mak_ScalarOp2(Vec, operator-, SubExp);



Mak_Op(Mat, Mat, Mat, operator*, MulExp);
Mak_ScalarOp1(Mat, operator*, ScalarMulExp);
Mak_ScalarOp2(Mat, operator*, ScalarMulExp);
Mak_Op(Vec, Vec, Vec, operator*, MulExp);
Mak_Op(Mat, Vec, Vec, operator*, MulExp);
Mak_Op(Vec, Mat, Vec, operator*, MulExp);
Mak_ScalarOp1(Vec, operator*, ScalarMulExp);
Mak_ScalarOp2(Vec, operator*, ScalarMulExp);


Mak_ScalarOp1(Mat, operator/, ScalarDivExp);
Mak_ScalarOp1(Vec, operator/, ScalarDivExp);

// Clean up the Macro
#undef Mak_Op
#undef Mak_ScalarOp1
#undef Mak_ScalarOp2



namespace exmat {


template<class MAT>
struct ConstErrChkedRow {
    typedef typename MAT::index_type        index_type;
    typedef typename MAT::reference         reference;
    typedef typename MAT::const_reference   const_reference;
    typedef typename MAT::ErrChecker        ErrChecker;

    EXMAT_INLINE2 ConstErrChkedRow(const MAT& mat_, index_type i)
        : mat(mat_), row(i) {}
    EXMAT_INLINE2 const_reference operator[](const index_type i) const {
        ErrChecker::AssertBounds(index_type(0), index_type(mat.cols()-1), i);
        return mat.getAt(row, i);
    }
    const MAT& mat;
    index_type row;
};


template<class MAT>
struct ErrChkedRow {
    typedef typename MAT::index_type        index_type;
    typedef typename MAT::reference         reference;
    typedef typename MAT::const_reference   const_reference;
    typedef typename MAT::ErrChecker        ErrChecker;

    EXMAT_INLINE2 ErrChkedRow(MAT& mat_, index_type i)
        : mat(mat_), row(i) {}
    EXMAT_INLINE2 reference operator[](const index_type i) {
        ErrChecker::AssertBounds(index_type(0), index_type(mat.cols()-1), i);
        return mat.setAt(row, i);
    }
    EXMAT_INLINE2 const_reference operator[](const index_type i) const {
        ErrChecker::AssertBounds(index_type(0), index_type(mat.cols()-1), i);
        return mat.getAt(row, i);
    }
    MAT& mat;
    index_type row;
};



template<class Rep, class ErrChecker_=DefaultErrorChecker>
class Mat : public Rep {
public:
    typedef Rep                     super_type;
    typedef Mat<Rep, ErrChecker_>   self_type;
    enum {
        IsFullErrorChecker = ISSAMETYPE<ErrChecker_,FullErrorChecker>::RET,
//      IsDense = HasTagOf<Rep, DenseTag>::RET,
        IsDR = (Rep::ROWS == 0),
        IsDC = (Rep::COLS == 0)
    };
    typedef typename
    IF<IsFullErrorChecker, ErrChkedRow<self_type>, typename Rep::row_type
    >::RET row_type;
    typedef typename
    IF<IsFullErrorChecker, ConstErrChkedRow<self_type>, typename Rep::const_row_type
    >::RET const_row_type;
    typedef typename Rep::col_type          col_type;
    typedef typename Rep::value_type        value_type;
    typedef typename Rep::index_type        index_type;
    typedef typename Rep::reference         reference;
    typedef typename Rep::const_reference   const_reference;
    typedef typename Rep::param_type        param_type;
    typedef typename Rep::const_param_type  const_param_type;
    typedef ErrChecker_                     ErrChecker;
    typedef typename Rep::rep_type          rep_type;

private:

    EXMAT_INLINE2 row_type IndexOp(index_type i, Int2Type<true>) {  // Error checked
        ErrChecker::AssertBounds(index_type(0), index_type(Rep::rows()-1), i);
        return row_type(*this, i);
    }
    EXMAT_INLINE2 row_type IndexOp(index_type i, Int2Type<false>) { // No error checking
        return Rep::operator[](i);
    }
    EXMAT_INLINE2 const_row_type IndexOp(index_type i, Int2Type<true>) const {
        ErrChecker::AssertBounds(index_type(0), index_type(Rep::rows()-1), i);
        return const_row_type(*this, i);
    }
    EXMAT_INLINE2 const_row_type IndexOp(index_type i, Int2Type<false>) const {
        return Rep::operator[](i);
    }

private:

    EXMAT_INLINE2 void CTORresizer(size_t r, size_t c, Int2Type<true>, Int2Type<true>)  // Dynamic row, dynamic col
    {   Rep::resize(r, c);  }
    EXMAT_INLINE2 void CTORresizer(size_t r, size_t c, Int2Type<true>, Int2Type<false>) // Dynamic row, static col
    {   Rep::resize(r); }
    EXMAT_INLINE2 void CTORresizer(size_t r, size_t c, Int2Type<false>, Int2Type<true>) // Static row, dynamic col
    {   Rep::resize(c); }
    EXMAT_INLINE2 void CTORresizer(size_t r, size_t c, Int2Type<false>, Int2Type<false>)// Static row, static col
    {}

public:
    /************************************************/
    EXMAT_INLINE2 Mat() {}

    // Helper class for enable_if
    template<typename T>
    struct _AlwaysTrue {    enum {RET=true};    };
    template<class Rep2, class ErrChk2> EXMAT_INLINE2
    Mat(const Mat<Rep2, ErrChk2>& A
#if !EXMAT_NO_SFINAE
        ,typename DisableIf_c<  // The SFINAE principle need Rep2 to be referenced
            IsDerivedFrom<Rep,ViewTag>::RET && _AlwaysTrue<Rep2>::RET
        >::RET* dummy=NULL
#endif  // !EXMAT_NO_SFINAE
    ) {
        // Do resizing of this matrix
        CTORresizer(A.rows(), A.cols(), Int2Type<IsDR>(), Int2Type<IsDC>());
        assign(*this, A);
//      AssignDispatcher(*this, A);
    }


    template<class Rep2> EXMAT_INLINE2
    Mat(ExpMat<Rep2>& A) : Rep(A) {}
    template<class Rep2> EXMAT_INLINE2
    Mat(const ExpMat<Rep2>& A) : Rep(A) {}


//  Mat(typename Rep::super_type& A) : Rep(A) {}
//  Mat(const typename Rep::super_type& A) : Rep(A) {}

    EXMAT_INLINE2 Mat(const std::string& s) {
        // It's only for dense matrix!
//      STATIC_CHECK(IsDense, Mat_StrCtrOnlyForDenseMat);
        std::istringstream ins(s);
        ::operator>>(ins, *this);
    }


    template<class T1, class T2>
    EXMAT_INLINE2 Mat(T1& t1, T2& t2) : Rep(t1, t2) {}
    template<class T1, class T2, class T3>
    EXMAT_INLINE2 Mat(T1& t1, T2& t2, T3& t3) : Rep(t1, t2, t3) {}
    template<class T1, class T2, class T3, class T4>
    EXMAT_INLINE2 Mat(T1& t1, T2& t2, T3& t3, T4& t4) : Rep(t1, t2, t3, t4) {}
    template<class T1, class T2, class T3, class T4, class T5>
    EXMAT_INLINE2 Mat(T1& t1, T2& t2, T3& t3, T4& t4, T5& t5) : Rep(t1, t2, t3, t4, t5) {}
    template<class T1, class T2, class T3, class T4, class T5, class T6, class T7, class T8, class T9>
    EXMAT_INLINE2 Mat(T1& t1, T2& t2, T3& t3, T4& t4, T5& t5, T6& t6, T7& t7, T8& t8, T9& t9) : Rep(t1, t2, t3, t4, t5, t6, t7, t8, t9) {}
    template<class T1, class T2, class T3, class T4, class T5, class T6, class T7, class T8, class T9, class T10, class T11, class T12, class T13, class T14, class T15, class T16>
    EXMAT_INLINE2 Mat(T1& t1, T2& t2, T3& t3, T4& t4, T5& t5, T6& t6, T7& t7, T8& t8, T9& t9, T10 t10, T11& t11, T12& t12, T13& t13, T14& t14, T15& t15, T16& t16) : Rep(t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, t12, t13, t14, t15, t16) {}
    // Const version of generic constructor, fail in VC6
#if HAVE_PARTIAL_SPECIALIZATION
    template<class T1, class T2>
    EXMAT_INLINE2 Mat(const T1& t1, const T2& t2) : Rep(t1, t2) {}
    template<class T1, class T2, class T3>
    EXMAT_INLINE2 Mat(const T1& t1, const T2& t2, const T3& t3) : Rep(t1, t2, t3) {}
    template<class T1, class T2, class T3, class T4>
    EXMAT_INLINE2 Mat(const T1& t1, const T2& t2, const T3& t3, const T4& t4) : Rep(t1, t2, t3, t4) {}
    template<class T1, class T2, class T3, class T4, class T5>
    EXMAT_INLINE2 Mat(const T1& t1, const T2& t2, const T3& t3, const T4& t4, const T5& t5) : Rep(t1, t2, t3, t4, t5) {}
#endif  // HAVE_PARTIAL_SPECIALIZATION



    EXMAT_INLINE2 explicit Mat(const Dim& d) : Rep(d) {}
    EXMAT_INLINE2 explicit Mat(const Dim& d1, const Dim& d2, const_param_type s) : Rep(d1, d2, s) {}

    EXMAT_INLINE2 Mat(const_param_type s) : Rep(s) {}

    /************************************************/

    // The scalar assignment operator is replaced by the CommaInitializer,
    // The same effect get be achieved by using ScalarInit
    //  EXMAT_INLINE2 self_type& operator=(const value_type& v_) {
    //      (*static_cast<Rep*>(this)) = v_;
    //      return *this;
    //  }

    EXMAT_INLINE2 self_type& ScalarInit(const_param_type s) {
        (*static_cast<Rep*>(this)) = s;
        return *this;
    }
    EXMAT_INLINE2 CommaInitializer<self_type> operator=(const_param_type s) {
        return CommaInitializer<self_type>(*this, s);
    }

    template<class Res, class Exp> EXMAT_INLINE2 static
    void AssignDispatcher(Res& res, const Exp& exp) {
        enum {
            IsSIMD =
                // Both left and right hand side should have SIMDTag
                IsDerivedFrom<Res,SIMDTag>::RET &&
                IsDerivedFrom<Exp,SIMDTag>::RET &&
                // If EXP is not linear then it needs ENABLE_TMP_ANALYSIS
                (ENABLE_TMP_ANALYSIS || Exp::IsLinearRecursive)
        };
#if ENABLE_UNROLLING
        typedef LinearBinExpSize<Res, Exp> BINEXPSize;
        enum {
            // Prefer unroll row before col for row major storage
            UnrollRow = 
                BINEXPSize::ROWS > 0 &&
                (unsigned int)(Exp::ENOP * BINEXPSize::COLS) < (unsigned int)MAX_UROLLROW_NOP
        };
        enum {
            UnrollCol =
                BINEXPSize::COLS > 0 &&
                (unsigned int)Exp::TNOP < (unsigned int)MAX_FULLUROLL_NOP
        };

        enum {
            AssignDispatchID =
                IsSIMD ? exmat::PNS::ASSIGN_TYPE::Ass_SIMD : 
                UnrollRow + 2*UnrollCol
        };
#else
        enum {
            AssignDispatchID =
                IsSIMD ? exmat::PNS::ASSIGN_TYPE::Ass_SIMD : 
                PNS::ASSIGN_TYPE::Ass_NoUnRoll
        };
#endif

        // Do compatibility check here, dimension should be the same
        ErrChecker::AssertCompat(res.rows() == exp.rows() && res.cols() == exp.cols());

        exmat::PNS::Assign(res, exp, Int2Type<AssignDispatchID>());
//      exmat::PNS::Assign(res, exp, Int2Type<false>(), Int2Type<false>());

        // Original assignment loops
//      for(index_type i=res.rows(); i--; ) {
//          for(index_type j=res.cols(); j--; )
//              res.setAt(i, j) = exp.getAt(i, j);
//      }
    }

    template<class Rep2, class ErrChk2> EXMAT_INLINE2
    Mat& AssignNoAna(const Mat<Rep2, ErrChk2>& A) {
        AssignDispatcher(*this, A);
        return *this;
    }
    template<class Rep2, class ErrChk2> EXMAT_INLINE2
    Mat& operator=(const Mat<Rep2, ErrChk2>& A) {
        // Using the global assign function with the default analysis option
        assign(*this, A, EXMAT_DEFAULT_ANALYSIS);
        return *this;
    }
    EXMAT_INLINE2 Mat& operator=(const self_type& A) {
        // Forward the assignment operation to the rep_type
        *((Rep*)(this)) = A;
        return *this;
    }

    /************************************************/
    EXMAT_INLINE2 reference operator() (index_type r, index_type c) {
        ErrChecker::AssertBounds(index_type(0), index_type(Rep::rows()-1), r);
        ErrChecker::AssertBounds(index_type(0), index_type(Rep::cols()-1), c);
        return setAt(r, c);
    }
    EXMAT_INLINE2 const_reference operator() (index_type r, index_type c) const {
        ErrChecker::AssertBounds(index_type(0), index_type(Rep::rows()-1), r);
        ErrChecker::AssertBounds(index_type(0), index_type(Rep::cols()-1), c);
        return getAt(r, c);
    }

    EXMAT_INLINE2 row_type operator[] (index_type r) {
        return IndexOp(r, Int2Type<IsFullErrorChecker>());
    }
    EXMAT_INLINE2 const_row_type operator[] (index_type r) const {
        return IndexOp(r, Int2Type<IsFullErrorChecker>());
    }

    /************************************************/

    EXMAT_INLINE2 self_type& operator+=(const_param_type s) {
        *this = ::operator+(*this, s);
        return *this;
    }
    EXMAT_INLINE2 self_type& operator-=(const_param_type s) {
        *this = ::operator-(*this, s);
        return *this;
    }
    EXMAT_INLINE2 self_type& operator*=(const_param_type s) {
        *this = ::operator*(*this, s);
        return *this;
    }
    EXMAT_INLINE2 self_type& operator/=(const_param_type s) {
        *this = ::operator/(*this, s);
        return *this;
    }
    template<class Rep2, class ErrChk2> EXMAT_INLINE2
    self_type& operator+=(const Mat<Rep2,ErrChk2>& A) {
        *this = ::operator+(*this, A);
        return *this;
    }
    template<class Rep2, class ErrChk2> EXMAT_INLINE2
    self_type& operator-=(const Mat<Rep2,ErrChk2>& A) {
        *this = ::operator-(*this, A);
        return *this;
    }
    template<class Rep2, class ErrChk2> EXMAT_INLINE2
    self_type& operator*=(const Mat<Rep2,ErrChk2>& A) {
        // Dynamic sized matrix is handled by using direct initialization
        self_type tmp( ::operator*(*this, A) );
        *this = tmp;
        return *this;
    }

    /************************************************/
    template<class Rep2, class ErrChk2> EXMAT_INLINE2
    bool operator!=(const Mat<Rep2,ErrChk2>& A) const {
        return !(*this == A);
    }
    template<class Rep2, class ErrChk2> EXMAT_INLINE2
    bool operator==(const Mat<Rep2,ErrChk2>& A) const
    {
        typedef typename Rep2::value_type T2;
        // At least the dimension should be the same
        if(Rep::rows() != A.rows() || Rep::cols() != A.cols())
            return false;
        typedef LEleType(value_type,T2) lower_type;
        for(index_type i=Rep::rows(); i--; ) {
            for(index_type j=Rep::cols(); j--; )
                if( lower_type(getAt(i,j)) != lower_type(A.getAt(i,j)) )
                    return false;
        }
        return true;
    }
    EXMAT_INLINE2 bool operator!=(const value_type& s) const {
        return !(*this == s);
    }
    friend EXMAT_INLINE2 bool operator!=(const value_type& lhs, const self_type& rhs) {
        return rhs != lhs;
    }
    EXMAT_INLINE2 bool operator==(const value_type& s) const
    {
        for(index_type i=Rep::rows(); i--; ) {
            for(index_type j=Rep::cols(); j--; )
                if( getAt(i,j) != s )
                    return false;
        }
        return true;
    }
    friend EXMAT_INLINE2 bool operator==(const value_type& lhs, const self_type& rhs) {
        return rhs == lhs;
    }

//protected:
    /************************************************/
    EXMAT_INLINE2 reference setAt(index_type i)
    {   return Rep::setAt(i);   }
    EXMAT_INLINE2 const_reference getAt(index_type i) const
    {   return Rep::getAt(i);   }
    EXMAT_INLINE2 reference setAt(index_type r, index_type c)
    {   return Rep::setAt(r, c);    }
    EXMAT_INLINE2 const_reference getAt(index_type r, index_type c) const
    {   return Rep::getAt(r, c);    }

    /************************************************/
    // Make some friends
//  friend std::ostream& operator<<(std::ostream &os, const self_type& mat);
};


/************************************************/
template<class Rep, class ErrChecker> EXMAT_INLINE2 static
const Mat<ExpMat<SubExp<ScalarCon<typename Rep::value_type>, Rep> >, ErrChecker>
operator-(const Mat<Rep, ErrChecker>& mat) {
    typedef typename Rep::value_type value_type;
    // --! Only simply return 0 - mat, which is not efficient
    return Mat<ExpMat<SubExp<ScalarCon<value_type>, Rep> >, ErrChecker>(ScalarCon<value_type>(0), mat);
}



template<class Rep, class ErrChecker_=DefaultErrorChecker>
class Vec : public Mat<Rep, ErrChecker_> {
public:
    typedef Mat<Rep, ErrChecker_>           super_type;
    typedef Vec<Rep, ErrChecker_>           self_type;
    typedef typename Rep::row_type          row_type;
    typedef typename Rep::const_row_type    const_row_type;
    typedef typename Rep::col_type          col_type;
    typedef typename Rep::value_type        value_type;
    typedef typename Rep::index_type        index_type;
    typedef typename Rep::reference         reference;
    typedef typename Rep::const_reference   const_reference;
    typedef typename Rep::param_type        param_type;
    typedef typename Rep::const_param_type  const_param_type;
    typedef ErrChecker_                     ErrChecker;
    typedef typename Rep::rep_type          rep_type;

    enum {
        IsRowVec = (Rep::ROWS == 1),
        IsColVec = (Rep::COLS == 1)
    };

public:
    /************************************************/
    // Constructors
    EXMAT_INLINE2 Vec() {}


#if defined( EXMAT_VC6 ) || defined( EXMAT_VC7 )
    //  VC6 have problems on distinguishing const and non-const template function
    //  therefore extra functions are need for single parameter
    template<class Rep2, class ErrChk2> EXMAT_INLINE2
    Vec(Vec<Rep2, ErrChk2>& vec) : super_type(vec) {}
    template<class Rep2, class ErrChk2> EXMAT_INLINE2
    Vec(const Vec<Rep2, ErrChk2>& vec) : super_type(vec) {}
    template<class Rep2, class ErrChk2> EXMAT_INLINE2
    Vec(Mat<Rep2, ErrChk2>& vec) : super_type(vec) {}
    template<class Rep2, class ErrChk2> EXMAT_INLINE2
    Vec(const Mat<Rep2, ErrChk2>& vec) : super_type(vec) {}
    template<class Rep2> EXMAT_INLINE2
    Vec(ExpMat<Rep2>& A) : super_type(A) {}
    template<class Rep2> EXMAT_INLINE2
    Vec(const ExpMat<Rep2>& A) : super_type(A) {}
    EXMAT_INLINE2 Vec(const std::string& s) : super_type(s) {}
    EXMAT_INLINE2 Vec(const Dim& d) : super_type(d) {}
#else
    template<typename T1> EXMAT_INLINE2
    Vec(T1& t1) : super_type(t1) {}
    template<typename T1> EXMAT_INLINE2
    Vec(const T1& t1) : super_type(t1) {}
//  template<typename T1, typename T2, typename T3> EXMAT_INLINE2
//  Vec(T1 t1, T2 t2, T3 t3) : super_type(t1, t2, t3) {}
//  template<typename T1, typename T2, typename T3, typename T4> EXMAT_INLINE2
//  Vec(T1 t1, T2 t2, T3 t3, T4 t4) : super_type(t1, t2, t3, t4) {}
#endif
    template<typename T1, typename T2> EXMAT_INLINE2
    Vec(const T1& t1, const T2& t2) : super_type(t1, t2) {}
    template<typename T1, typename T2, typename T3> EXMAT_INLINE2
    Vec(T1& t1, T2& t2, T3& t3) : super_type(t1, t2, t3) {}
    template<typename T1, typename T2, typename T3, typename T4> EXMAT_INLINE2
    Vec(T1& t1, T2& t2, T3& t3, T4& t4) : super_type(t1, t2, t3, t4) {}

    EXMAT_INLINE2 Vec(const_param_type s) : super_type(s) {}

    /************************************************/

    EXMAT_INLINE2 CommaInitializer<self_type> operator=(const_param_type v_) {
        return CommaInitializer<self_type>(*this, v_);
    }
    template<class Rep2, class ErrChk2> EXMAT_INLINE2
    Vec& operator=(const Vec<Rep2, ErrChk2>& A) {
        ((super_type*)(this))->operator=(A);
        return *this;
    }
    template<class Rep2, class ErrChk2> EXMAT_INLINE2
    Vec& operator=(const Mat<Rep2, ErrChk2>& A) {
        ((super_type*)(this))->operator=(A);
        return *this;
    }
    EXMAT_INLINE2 Vec& operator=(const self_type& A) {
        // Forward the assignment operation to the rep_type
        ((super_type*)(this))->operator=(A);
        return *this;
    }

    /************************************************/
    EXMAT_INLINE2 reference operator() (index_type i) {
        ErrChecker::AssertBounds(index_type(0), index_type(this->size()-1), i);
        return setAt(i);
    }
    EXMAT_INLINE2 const_reference operator() (index_type i) const {
        ErrChecker::AssertBounds(index_type(0), index_type(this->size()-1), i);
        return getAt(i);
    }
    EXMAT_INLINE2 reference operator[] (index_type i) {
        ErrChecker::AssertBounds(index_type(0), index_type(this->size()-1), i);
        return setAt(i);
    }
    EXMAT_INLINE2 const_reference operator[] (index_type i) const {
        ErrChecker::AssertBounds(index_type(0), index_type(this->size()-1), i);
        return getAt(i);
    }

};

/************************************************/
template<class Rep, class ErrChecker> EXMAT_INLINE2 static
const Vec<ExpMat<SubExp<ScalarCon<typename Rep::value_type>, typename Rep::rep_type> >, ErrChecker>
operator-(const Vec<Rep, ErrChecker>& vec) {
    typedef typename Rep::value_type value_type;
    // --! Only simply return 0 - mat, which is not efficient
    return Vec<ExpMat<SubExp<ScalarCon<value_type>, TTYPENAME Rep::rep_type> >, ErrChecker>(ScalarCon<value_type>(0), vec);
}

}   // namespace exmat


#define EXMAT_DERIVEDCLASS_IMP(NAME, MatOrVec)                  \
public:                                                         \
typedef typename super_type::rep_type   rep_type;               \
typedef typename super_type::ErrChecker ErrChecker;             \
typedef typename super_type::value_type value_type;             \
EXMAT_INLINE2 NAME(value_type s) : super_type(s) {}             \
EXMAT_INLINE2 NAME(const std::string& s) : super_type(s) {}     \
template<class Rep2, class EC2> EXMAT_INLINE2                   \
NAME(::exmat::MatOrVec<Rep2,EC2>& A) : super_type(A) {}         \
template<class Rep2, class EC2> EXMAT_INLINE2                   \
NAME(const ::exmat::MatOrVec<Rep2,EC2>& A) : super_type(A) {}   \
template<class Rep2> EXMAT_INLINE2                              \
NAME(::exmat::ExpMat<Rep2>& A) : super_type(A) {}               \
template<class Rep2> EXMAT_INLINE2                              \
NAME(const ::exmat::ExpMat<Rep2>& A) : super_type(A) {}         \
template<class T1, class T2> EXMAT_INLINE2                      \
NAME(const T1& t1, const T2& t2) : super_type(t1, t2) {}        \
template<class T1, class T2, class T3> EXMAT_INLINE2            \
NAME(const T1& t1, const T2& t2, const T3& t3)                  \
 : super_type(t1, t2, t3) {}                                    \
template<class T1, class T2, class T3, class T4> EXMAT_INLINE2  \
NAME(const T1& t1, const T2& t2, const T3& t3, const T4& t4)    \
 : super_type(t1, t2, t3, t4) {}                                \
private:

#define EXMAT_DERIVEDCLASS_ASSIGNOP(NAME, MatOrVec)             \
public:                                                         \
template<class Rep2, class EC> EXMAT_INLINE2                    \
self_type& operator=(const ::exmat::MatOrVec<Rep2, EC>& v)      \
{   super_type::operator=(v);   return *this;   }               \
EXMAT_INLINE2 self_type& operator=(value_type s)                \
{   super_type::operator=(s);   return *this;   }


#define EXMAT_DERIVEDCLASS_Mak_Op(LT, RT, OT, op, OpName)   \
template<class L, class R, class LC, class RC> EXMAT_INLINE2    \
    const OT< ::exmat::ExpMat<OpName<typename L::rep_type,typename R::rep_type> >, typename ::exmat::ResultantChecker<LC,RC>::RET> >    \
op (const LT<L,LC> >& A, const RT<R,RC> >& B) { \
    typedef ::exmat::ExpMat<OpName<TTYPENAME L::rep_type,TTYPENAME R::rep_type> > RETType;  \
    return RETType(A, B);                                       \
}
#define EXMAT_DERIVEDCLASS_Mak_ScalarOp1(OT, op, OpName)    \
template<class L, class ErrChk> EXMAT_INLINE2   \
    const OT< ::exmat::ExpMat<OpName<typename L::rep_type,::exmat::ScalarCon<typename L::value_type> > >, ErrChk> > \
op (const OT<L,ErrChk> >& A, typename L::const_param_type scalar) { \
    typedef ::exmat::ExpMat<OpName<TTYPENAME L::rep_type,::exmat::ScalarCon<TTYPENAME L::value_type> > > RETType;   \
    return RETType(A, ::exmat::ScalarCon<TTYPENAME L::value_type>(scalar)); \
}
#define EXMAT_DERIVEDCLASS_Mak_ScalarOp2(OT, op, OpName)    \
template<class L, class ErrChk> EXMAT_INLINE2   \
    const OT< ::exmat::ExpMat<OpName< ::exmat::ScalarCon<typename L::value_type>, typename L::rep_type> >, ErrChk> >    \
op (typename L::const_param_type scalar, const OT<L,ErrChk> >& A) { \
    typedef ::exmat::ExpMat<OpName< ::exmat::ScalarCon<TTYPENAME L::value_type>, TTYPENAME L::rep_type> > RETType;  \
    return RETType(exmat::ScalarCon<TTYPENAME L::value_type>(scalar), A);   \
}

#define EXMAT_DERIVEDCLASS_DEFALUTOPS(ClassName, MatOrVec)  \
EXMAT_DERIVEDCLASS_Mak_Op(ClassName< ::exmat::MatOrVec, ClassName< ::exmat::MatOrVec, ClassName< ::exmat::MatOrVec, operator+, ::exmat::AddExp);    \
EXMAT_DERIVEDCLASS_Mak_ScalarOp1(ClassName< ::exmat::MatOrVec, operator+, ::exmat::AddExp);         \
EXMAT_DERIVEDCLASS_Mak_ScalarOp2(ClassName< ::exmat::MatOrVec, operator+, ::exmat::AddExp);         \
EXMAT_DERIVEDCLASS_Mak_Op(ClassName< ::exmat::MatOrVec, ClassName< ::exmat::MatOrVec, ClassName< ::exmat::MatOrVec, operator-, ::exmat::SubExp);    \
EXMAT_DERIVEDCLASS_Mak_ScalarOp1(ClassName< ::exmat::MatOrVec, operator-, ::exmat::SubExp);         \
EXMAT_DERIVEDCLASS_Mak_ScalarOp2(ClassName< ::exmat::MatOrVec, operator-, ::exmat::SubExp);         \
EXMAT_DERIVEDCLASS_Mak_ScalarOp1(ClassName< ::exmat::MatOrVec, operator*, ::exmat::ScalarMulExp);   \
EXMAT_DERIVEDCLASS_Mak_ScalarOp2(ClassName< ::exmat::MatOrVec, operator*, ::exmat::ScalarMulExp);   \
EXMAT_DERIVEDCLASS_Mak_ScalarOp1(ClassName< ::exmat::MatOrVec, operator/, ::exmat::ScalarDivExp);   \
EXMAT_DERIVEDCLASS_Mak_ScalarOp2(ClassName< ::exmat::MatOrVec, operator/, ::exmat::ScalarDivExp);


#ifdef _MSC_VER
#   pragma warning( pop )
#endif  // _MSC_VER


#endif  // _EXMAT_MAT_H

---