// Copyright (C) 2006 Davis E. King (davisking@users.sourceforge.net)
// License: Boost Software License See LICENSE.txt for the full license.
#ifndef DLIB_MATRIx_
#define DLIB_MATRIx_
#include "matrix_abstract.h"
#include "../algs.h"
#include "../serialize.h"
#include "../enable_if.h"
#include <sstream>
#include <algorithm>
#include "../memory_manager.h"
#include "../is_kind.h"
#ifdef _MSC_VER
// Disable the following warnings for Visual Studio
// This warning is:
// "warning C4355: 'this' : used in base member initializer list"
// Which we get from this code but it is not an error so I'm turning this
// warning off and then turning it back on at the end of the file.
#pragma warning(disable : 4355)
#endif
namespace dlib
{
// ----------------------------------------------------------------------------------------
template <
typename T,
long num_rows = 0,
long num_cols = 0,
typename mem_manager = memory_manager<char>::kernel_1a
>
class matrix;
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
template <
typename T,
long num_rows,
long num_cols,
typename mem_manager
>
class matrix_ref
{
public:
typedef T type;
typedef matrix_ref ref_type;
typedef mem_manager mem_manager_type;
const static long NR = num_rows;
const static long NC = num_cols;
matrix_ref (
const matrix<T,num_rows,num_cols,mem_manager>& m_
) : m(m_) {}
matrix_ref (
const matrix_ref& i_
) : m(i_.m) {}
const T& operator() (
long r,
long c
) const { return m(r,c); }
long nr (
) const { return m.nr(); }
long nc (
) const { return m.nc(); }
long size (
) const { return m.size(); }
template <typename U, long iNR, long iNC, typename mm >
bool aliases (
const matrix<U,iNR,iNC,mm>& item
) const { return false; }
template <typename U, long iNR, long iNC, typename mm>
bool destructively_aliases (
const matrix<U,iNR,iNC,mm>& item
) const { return false; }
bool aliases (
const matrix<T,num_rows,num_cols,mem_manager>& item
) const { return (&m == &item); }
const matrix_ref ref(
) const { return *this; }
private:
// no assignment operator
matrix_ref& operator=(const matrix_ref&);
const matrix<T,num_rows,num_cols,mem_manager>& m; // This is the item contained by this expression.
};
// ----------------------------------------------------------------------------------------
// this is a hack to avoid a compile time error in visual studio 8. I would just
// use sizeof(T) and be done with it but that won't compile. The idea here
// is to avoid using the stack allocation of the matrix_data object if it
// is going to contain another matrix and also avoid asking for the sizeof()
// the contained matrix.
template <typename T>
struct get_sizeof_helper
{
const static std::size_t val = sizeof(T);
};
template <typename T, long NR, long NC, typename mm>
struct get_sizeof_helper<matrix<T,NR,NC,mm> >
{
const static std::size_t val = 1000000;
};
template <
typename T,
long num_rows,
long num_cols,
typename mem_manager,
int val = static_switch <
// when the sizes are all non zero and small
(num_rows*num_cols*get_sizeof_helper<T>::val <= 64) && (num_rows != 0 && num_cols != 0),
// when the sizes are all non zero and big
(num_rows*num_cols*get_sizeof_helper<T>::val >= 65) && (num_rows != 0 && num_cols != 0),
num_rows == 0 && num_cols != 0,
num_rows != 0 && num_cols == 0,
num_rows == 0 && num_cols == 0
>::value
>
class matrix_data ;
/*!
WHAT THIS OBJECT REPRESENTS
This object represents the actual allocation of space for a matrix.
Small matrices allocate all their data on the stack and bigger ones
use a memory_manager to get their memory.
!*/
// ----------------------------------------------------------------------------------------
template <
typename T,
long num_rows,
long num_cols,
typename mem_manager
>
class matrix_data<T,num_rows,num_cols,mem_manager,1> : noncopyable // when the sizes are all non zero and small
{
public:
const static long NR = num_rows;
const static long NC = num_cols;
matrix_data() {}
T& operator() (
long r,
long c
) { return data[r][c]; }
const T& operator() (
long r,
long c
) const { return data[r][c]; }
T& operator() (
long i
) { return *(*data + i); }
const T& operator() (
long i
) const { return *(*data + i); }
void swap(
matrix_data& item
)
{
for (long r = 0; r < num_rows; ++r)
{
for (long c = 0; c < num_cols; ++c)
{
exchange((*this)(r,c),item(r,c));
}
}
}
long nr (
) const { return num_rows; }
long nc (
) const { return num_cols; }
void set_size (
long nr,
long nc
)
{
}
private:
T data[num_rows][num_cols];
};
// ----------------------------------------------------------------------------------------
template <
typename T,
long num_rows,
long num_cols,
typename mem_manager
>
class matrix_data<T,num_rows,num_cols,mem_manager,2> : noncopyable // when the sizes are all non zero and big
{
public:
const static long NR = num_rows;
const static long NC = num_cols;
matrix_data (
) { data = pool.allocate_array(num_rows*num_cols); }
~matrix_data ()
{ pool.deallocate_array(data); }
T& operator() (
long r,
long c
) { return data[r*num_cols + c]; }
const T& operator() (
long r,
long c
) const { return data[r*num_cols + c]; }
T& operator() (
long i
) { return data[i]; }
const T& operator() (
long i
) const { return data[i]; }
void swap(
matrix_data& item
)
{
std::swap(item.data,data);
pool.swap(item.pool);
}
long nr (
) const { return num_rows; }
long nc (
) const { return num_cols; }
void set_size (
long nr,
long nc
)
{
}
private:
T* data;
typename mem_manager::template rebind<T>::other pool;
};
// ----------------------------------------------------------------------------------------
template <
typename T,
long num_rows,
long num_cols,
typename mem_manager
>
class matrix_data<T,num_rows,num_cols,mem_manager,3> : noncopyable // when num_rows == 0 && num_cols != 0,
{
public:
const static long NR = num_rows;
const static long NC = num_cols;
matrix_data (
):data(0), nr_(0) { }
~matrix_data ()
{
if (data)
pool.deallocate_array(data);
}
T& operator() (
long r,
long c
) { return data[r*num_cols + c]; }
const T& operator() (
long r,
long c
) const { return data[r*num_cols + c]; }
T& operator() (
long i
) { return data[i]; }
const T& operator() (
long i
) const { return data[i]; }
void swap(
matrix_data& item
)
{
std::swap(item.data,data);
std::swap(item.nr_,nr_);
pool.swap(item.pool);
}
long nr (
) const { return nr_; }
long nc (
) const { return num_cols; }
void set_size (
long nr,
long nc
)
{
if (data)
{
pool.deallocate_array(data);
}
data = pool.allocate_array(nr*nc);
nr_ = nr;
}
private:
T* data;
long nr_;
typename mem_manager::template rebind<T>::other pool;
};
// ----------------------------------------------------------------------------------------
template <
typename T,
long num_rows,
long num_cols,
typename mem_manager
>
class matrix_data<T,num_rows,num_cols,mem_manager,4> : noncopyable // when num_rows != 0 && num_cols == 0
{
public:
const static long NR = num_rows;
const static long NC = num_cols;
matrix_data (
):data(0), nc_(0) { }
~matrix_data ()
{
if (data)
{
pool.deallocate_array(data);
}
}
T& operator() (
long r,
long c
) { return data[r*nc_ + c]; }
const T& operator() (
long r,
long c
) const { return data[r*nc_ + c]; }
T& operator() (
long i
) { return data[i]; }
const T& operator() (
long i
) const { return data[i]; }
void swap(
matrix_data& item
)
{
std::swap(item.data,data);
std::swap(item.nc_,nc_);
pool.swap(item.pool);
}
long nr (
) const { return num_rows; }
long nc (
) const { return nc_; }
void set_size (
long nr,
long nc
)
{
if (data)
{
pool.deallocate_array(data);
}
data = pool.allocate_array(nr*nc);
nc_ = nc;
}
private:
T* data;
long nc_;
typename mem_manager::template rebind<T>::other pool;
};
// ----------------------------------------------------------------------------------------
template <
typename T,
long num_rows,
long num_cols,
typename mem_manager
>
class matrix_data<T,num_rows,num_cols,mem_manager,5> : noncopyable // when num_rows == 0 && num_cols == 0
{
public:
const static long NR = num_rows;
const static long NC = num_cols;
matrix_data (
):data(0), nr_(0), nc_(0) { }
~matrix_data ()
{
if (data)
{
pool.deallocate_array(data);
}
}
T& operator() (
long r,
long c
) { return data[r*nc_ + c]; }
const T& operator() (
long r,
long c
) const { return data[r*nc_ + c]; }
T& operator() (
long i
) { return data[i]; }
const T& operator() (
long i
) const { return data[i]; }
void swap(
matrix_data& item
)
{
std::swap(item.data,data);
std::swap(item.nc_,nc_);
std::swap(item.nr_,nr_);
pool.swap(item.pool);
}
long nr (
) const { return nr_; }
long nc (
) const { return nc_; }
void set_size (
long nr,
long nc
)
{
if (data)
{
pool.deallocate_array(data);
}
data = pool.allocate_array(nr*nc);
nr_ = nr;
nc_ = nc;
}
private:
T* data;
long nr_;
long nc_;
typename mem_manager::template rebind<T>::other pool;
};
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
// We want to return the compile time constant if our NR and NC dimensions
// aren't zero but if they are then we want to call ref_.nx() and return
// the correct values.
template < typename ref_type, long NR >
struct get_nr_helper
{
static inline long get(const ref_type&) { return NR; }
};
template < typename ref_type >
struct get_nr_helper<ref_type,0>
{
static inline long get(const ref_type& m) { return m.nr(); }
};
template < typename ref_type, long NC >
struct get_nc_helper
{
static inline long get(const ref_type&) { return NC; }
};
template < typename ref_type >
struct get_nc_helper<ref_type,0>
{
static inline long get(const ref_type& m) { return m.nc(); }
};
// the matrix_exp for statically sized matrices
template <
typename EXP
>
class matrix_exp
{
public:
typedef typename EXP::type type;
typedef typename EXP::ref_type ref_type;
typedef typename EXP::mem_manager_type mem_manager_type;
const static long NR = EXP::NR;
const static long NC = EXP::NC;
typedef matrix<type,NR,NC,mem_manager_type> matrix_type;
matrix_exp (
const EXP& exp
) : ref_(exp.ref()) {}
inline const type operator() (
long r,
long c
) const
{
DLIB_ASSERT(r < nr() && c < nc() && r >= 0 && c >= c,
"\tconst type matrix_exp::operator(r,c)"
<< "\n\tYou must give a valid row and column"
<< "\n\tr: " << r
<< "\n\tc: " << c
<< "\n\tnr(): " << nr()
<< "\n\tnc(): " << nc()
<< "\n\tthis: " << this
);
return ref_(r,c);
}
const type operator() (
long i
) const
{
COMPILE_TIME_ASSERT(NC == 1 || NC == 0 || NR == 1 || NR == 0);
DLIB_ASSERT(nc() == 1 || nr() == 1,
"\tconst type matrix_exp::operator(i)"
<< "\n\tYou can only use this operator on column or row vectors"
<< "\n\ti: " << i
<< "\n\tnr(): " << nr()
<< "\n\tnc(): " << nc()
<< "\n\tthis: " << this
);
DLIB_ASSERT( ((nc() == 1 && i < nr()) || (nr() == 1 && i < nc())) && i >= 0,
"\tconst type matrix_exp::operator(i)"
<< "\n\tYou must give a valid row/column number"
<< "\n\ti: " << i
<< "\n\tnr(): " << nr()
<< "\n\tnc(): " << nc()
<< "\n\tthis: " << this
);
if (nc() == 1)
return ref_(i,0);
else
return ref_(0,i);
}
long size (
) const { return nr()*nc(); }
long nr (
) const { return get_nr_helper<ref_type,NR>::get(ref_); }
long nc (
) const { return get_nc_helper<ref_type,NC>::get(ref_); }
template <typename U, long iNR, long iNC, typename mm >
bool aliases (
const matrix<U,iNR,iNC,mm>& item
) const { return ref_.aliases(item); }
template <typename U, long iNR, long iNC , typename mm>
bool destructively_aliases (
const matrix<U,iNR,iNC,mm>& item
) const { return ref_.destructively_aliases(item); }
const ref_type& ref (
) const { return ref_; }
inline operator const type (
) const
{
COMPILE_TIME_ASSERT(NC == 1 || NC == 0);
COMPILE_TIME_ASSERT(NR == 1 || NR == 0);
DLIB_ASSERT(nr() == 1 && nc() == 1,
"\tmatrix_exp::operator const type&() const"
<< "\n\tYou can only use this operator on a 1x1 matrix"
<< "\n\tnr(): " << nr()
<< "\n\tnc(): " << nc()
<< "\n\tthis: " << this
);
return ref_(0,0);
}
private:
const ref_type ref_;
};
// ----------------------------------------------------------------------------------------
template <typename T>
struct is_matrix<matrix_exp<T> > { static const bool value = true; };
template <typename T, long NR, long NC, typename mm>
struct is_matrix<matrix_ref<T,NR,NC,mm> > { static const bool value = true; };
template <typename T, long NR, long NC, typename mm>
struct is_matrix<matrix<T,NR,NC,mm> > { static const bool value = true; };
template <typename T>
struct is_matrix<T&> { static const bool value = is_matrix<T>::value; };
template <typename T>
struct is_matrix<const T&> { static const bool value = is_matrix<T>::value; };
template <typename T>
struct is_matrix<const T> { static const bool value = is_matrix<T>::value; };
/*
is_matrix<T>::value == 1 if T is a matrix type else 0
*/
// ----------------------------------------------------------------------------------------
// This template will perform the needed loop for element multiplication using whichever
// dimension is provided as a compile time constant (if one is at all).
template <
typename LHS,
typename RHS,
long lhs_nc = LHS::NC,
long rhs_nr = RHS::NR
>
struct matrix_multiply_helper
{
typedef typename LHS::type type;
inline const static type eval (
const RHS& rhs,
const LHS& lhs,
long r,
long c
)
{
type temp = type();
for (long i = 0; i < rhs.nr(); ++i)
{
temp += lhs(r,i)*rhs(i,c);
}
return temp;
}
};
template <
typename LHS,
typename RHS,
long lhs_nc
>
struct matrix_multiply_helper <LHS,RHS,lhs_nc,0>
{
typedef typename LHS::type type;
inline const static type eval (
const RHS& rhs,
const LHS& lhs,
long r,
long c
)
{
type temp = type();
for (long i = 0; i < lhs.nc(); ++i)
{
temp += lhs(r,i)*rhs(i,c);
}
return temp;
}
};
template <
typename LHS,
typename RHS,
unsigned long count = 0
>
class matrix_multiply_exp
{
/*!
REQUIREMENTS ON LHS AND RHS
- they must be matrix_exp or matrix_ref objects (or
objects with a compatible interface).
!*/
public:
typedef typename LHS::type type;
typedef matrix_multiply_exp ref_type;
typedef typename LHS::mem_manager_type mem_manager_type;
const static long NR = LHS::NR;
const static long NC = RHS::NC;
matrix_multiply_exp (
const matrix_multiply_exp& item
) : lhs(item.lhs), rhs(item.rhs) {}
inline matrix_multiply_exp (
const LHS& lhs_,
const RHS& rhs_
) :
lhs(lhs_),
rhs(rhs_)
{
// You are trying to multiply two incompatible matrices together. The number of columns
// in the matrix on the left must match the number of rows in the matrix on the right.
COMPILE_TIME_ASSERT(LHS::NC == RHS::NR || LHS::NC*RHS::NR == 0);
DLIB_ASSERT(lhs.nc() == rhs.nr(),
"\tconst matrix_exp operator*(const matrix_exp& lhs, const matrix_exp& rhs)"
<< "\n\tYou are trying to multiply two incompatible matrices together"
<< "\n\tlhs.nr(): " << lhs.nr()
<< "\n\tlhs.nc(): " << lhs.nc()
<< "\n\trhs.nr(): " << rhs.nr()
<< "\n\trhs.nc(): " << rhs.nc()
<< "\n\t&lhs: " << &lhs
<< "\n\t&rhs: " << &rhs
);
// You can't multiply matrices together if they don't both contain the same type of elements.
COMPILE_TIME_ASSERT((is_same_type<typename LHS::type, typename RHS::type>::value == true));
}
inline const type operator() (
long r,
long c
) const
{
return matrix_multiply_helper<LHS,RHS>::eval(rhs,lhs,r,c);
}
long nr (
) const { return lhs.nr(); }
long nc (
) const { return rhs.nc(); }
template <typename U, long iNR, long iNC, typename mm >
bool aliases (
const matrix<U,iNR,iNC,mm>& item
) const { return lhs.aliases(item) || rhs.aliases(item); }
template <typename U, long iNR, long iNC , typename mm>
bool destructively_aliases (
const matrix<U,iNR,iNC,mm>& item
) const { return aliases(item); }
const ref_type& ref(
) const { return *this; }
const LHS lhs;
const RHS rhs;
};
template <
typename T,
long NR,
long NC,
typename EXP1,
typename EXP2,
typename MM
>
inline const matrix_exp<matrix_multiply_exp<EXP1, matrix_multiply_exp<EXP2,typename matrix<T,NR,NC,MM>::ref_type >,0 > > operator* (
const matrix_exp<matrix_multiply_exp<EXP1,EXP2,1> >& m1,
const matrix<T,NR,NC,MM>& m2
)
{
// We are going to reorder the order of evaluation of the terms here. This way the
// multiplication will go faster.
typedef matrix_multiply_exp<EXP2,typename matrix<T,NR,NC,MM>::ref_type > exp_inner;
typedef matrix_multiply_exp<EXP1, exp_inner,0 > exp_outer;
return matrix_exp<exp_outer>(exp_outer(m1.ref().lhs,exp_inner(m1.ref().rhs,m2)));
}
template <
typename EXP1,
typename EXP2
>
inline const matrix_exp<matrix_multiply_exp<EXP1, EXP2 > > operator* (
const matrix_exp<EXP1>& m1,
const matrix_exp<EXP2>& m2
)
{
typedef matrix_multiply_exp<EXP1, EXP2> exp;
return matrix_exp<exp>(exp(m1.ref(),m2.ref()));
}
template <
typename T,
long NR,
long NC,
typename EXP,
typename MM
>
inline const matrix_exp<matrix_multiply_exp<typename matrix<T,NR,NC,MM>::ref_type, matrix_exp<EXP> > > operator* (
const matrix<T,NR,NC,MM>& m1,
const matrix_exp<EXP>& m2
)
{
typedef matrix_multiply_exp<typename matrix<T,NR,NC,MM>::ref_type, matrix_exp<EXP> > exp;
return matrix_exp<exp>(exp(m1,m2));
}
template <
typename T,
long NR,
long NC,
typename EXP,
typename MM
>
inline const matrix_exp<matrix_multiply_exp< matrix_exp<EXP>, typename matrix<T,NR,NC,MM>::ref_type, 1> > operator* (
const matrix_exp<EXP>& m1,
const matrix<T,NR,NC,MM>& m2
)
{
typedef matrix_multiply_exp< matrix_exp<EXP>, typename matrix<T,NR,NC,MM>::ref_type, 1 > exp;
return matrix_exp<exp>(exp(m1,m2));
}
template <
typename T,
long NR1,
long NC1,
long NR2,
long NC2,
typename MM1,
typename MM2
>
inline const matrix_exp<matrix_multiply_exp<typename matrix<T,NR1,NC1,MM1>::ref_type,typename matrix<T,NR2,NC2,MM2>::ref_type > > operator* (
const matrix<T,NR1,NC1,MM1>& m1,
const matrix<T,NR2,NC2,MM2>& m2
)
{
typedef matrix_multiply_exp<typename matrix<T,NR1,NC1,MM1>::ref_type, typename matrix<T,NR2,NC2,MM2>::ref_type > exp;
return matrix_exp<exp>(exp(m1,m2));
}
// ----------------------------------------------------------------------------------------
template <
typename LHS,
typename RHS
>
class matrix_add_expression
{
/*!
REQUIREMENTS ON LHS AND RHS
- they must be matrix_exp or matrix_ref objects (or
objects with a compatible interface).
!*/
public:
typedef typename LHS::type type;
typedef typename LHS::mem_manager_type mem_manager_type;
typedef matrix_add_expression ref_type;
const static long NR = (RHS::NR > LHS::NR) ? RHS::NR : LHS::NR;
const static long NC = (RHS::NC > LHS::NC) ? RHS::NC : LHS::NC;
matrix_add_expression (
const matrix_add_expression& item
) : lhs(item.lhs), rhs(item.rhs) {}
matrix_add_expression (
const LHS& lhs_,
const RHS& rhs_
) :
lhs(lhs_),
rhs(rhs_)
{
// You can only add matrices together if they both have the same number of rows and columns.
COMPILE_TIME_ASSERT(LHS::NR == RHS::NR || LHS::NR == 0 || RHS::NR == 0);
COMPILE_TIME_ASSERT(LHS::NC == RHS::NC || LHS::NC == 0 || RHS::NC == 0);
DLIB_ASSERT(lhs.nc() == rhs.nc() &&
lhs.nr() == rhs.nr(),
"\tconst matrix_exp operator+(const matrix_exp& lhs, const matrix_exp& rhs)"
<< "\n\tYou are trying to add two incompatible matrices together"
<< "\n\tlhs.nr(): " << lhs.nr()
<< "\n\tlhs.nc(): " << lhs.nc()
<< "\n\trhs.nr(): " << rhs.nr()
<< "\n\trhs.nc(): " << rhs.nc()
<< "\n\t&lhs: " << &lhs
<< "\n\t&rhs: " << &rhs
);
// You can only add matrices together if they both contain the same types of elements.
COMPILE_TIME_ASSERT((is_same_type<typename LHS::type, typename RHS::type>::value == true));
}
const type operator() (
long r,
long c
) const { return lhs(r,c) + rhs(r,c); }
template <typename U, long iNR, long iNC , typename mm>
bool aliases (
const matrix<U,iNR,iNC,mm>& item
) const { return lhs.aliases(item) || rhs.aliases(item); }
template <typename U, long iNR, long iNC, typename mm >
bool destructively_aliases (
const matrix<U,iNR,iNC,mm>& item
) const { return lhs.destructively_aliases(item) || rhs.destructively_aliases(item); }
const ref_type& ref(
) const { return *this; }
long nr (
) const { return lhs.nr(); }
long nc (
) const { return lhs.nc(); }
const LHS lhs;
const RHS rhs;
};
template <
typename EXP1,
typename EXP2
>
inline const matrix_exp<matrix_add_expression<EXP1, EXP2 > > operator+ (
const matrix_exp<EXP1>& m1,
const matrix_exp<EXP2>& m2
)
{
typedef matrix_add_expression<EXP1, EXP2 > exp;
return matrix_exp<exp>(exp(m1.ref(),m2.ref()));
}
// ----------------------------------------------------------------------------------------
template <
typename LHS,
typename RHS
>
class matrix_subtract_exp
{
/*!
REQUIREMENTS ON LHS AND RHS
- they must be matrix_exp or matrix_ref objects (or
objects with a compatible interface).
!*/
public:
typedef typename LHS::type type;
typedef typename LHS::mem_manager_type mem_manager_type;
typedef matrix_subtract_exp ref_type;
const static long NR = (RHS::NR > LHS::NR) ? RHS::NR : LHS::NR;
const static long NC = (RHS::NC > LHS::NC) ? RHS::NC : LHS::NC;
matrix_subtract_exp (
const LHS& lhs_,
const RHS& rhs_
) :
lhs(lhs_),
rhs(rhs_)
{
// You can only subtract one matrix from another if they both have the same number of rows and columns.
COMPILE_TIME_ASSERT(LHS::NR == RHS::NR || LHS::NR == 0 || RHS::NR == 0);
COMPILE_TIME_ASSERT(LHS::NC == RHS::NC || LHS::NC == 0 || RHS::NC == 0);
DLIB_ASSERT(lhs.nc() == rhs.nc() &&
lhs.nr() == rhs.nr(),
"\tconst matrix_exp operator-(const matrix_exp& lhs, const matrix_exp& rhs)"
<< "\n\tYou are trying to add two incompatible matrices together"
<< "\n\tlhs.nr(): " << lhs.nr()
<< "\n\tlhs.nc(): " << lhs.nc()
<< "\n\trhs.nr(): " << rhs.nr()
<< "\n\trhs.nc(): " << rhs.nc()
<< "\n\t&lhs: " << &lhs
<< "\n\t&rhs: " << &rhs
);
// You can only subtract one matrix from another if they both contain elements of the same type.
COMPILE_TIME_ASSERT((is_same_type<typename LHS::type, typename RHS::type>::value == true));
}
const type operator() (
long r,
long c
) const { return lhs(r,c) - rhs(r,c); }
template <typename U, long iNR, long iNC, typename mm >
bool aliases (
const matrix<U,iNR,iNC, mm>& item
) const { return lhs.aliases(item) || rhs.aliases(item); }
template <typename U, long iNR, long iNC , typename mm>
bool destructively_aliases (
const matrix<U,iNR,iNC,mm>& item
) const { return lhs.destructively_aliases(item) || rhs.destructively_aliases(item); }
const ref_type& ref(
) const { return *this; }
long nr (
) const { return lhs.nr(); }
long nc (
) const { return lhs.nc(); }
const LHS lhs;
const RHS rhs;
};
template <
typename EXP1,
typename EXP2
>
inline const matrix_exp<matrix_subtract_exp<EXP1, EXP2 > > operator- (
const matrix_exp<EXP1>& m1,
const matrix_exp<EXP2>& m2
)
{
typedef matrix_subtract_exp<EXP1, EXP2 > exp;
return matrix_exp<exp>(exp(m1.ref(),m2.ref()));
}
// ----------------------------------------------------------------------------------------
template <
typename M,
typename S
>
class matrix_divscal_exp
{
/*!
REQUIREMENTS ON M
- must be a matrix_exp or matrix_ref object (or
an object with a compatible interface).
REQUIREMENTS ON S
- must be some kind of scalar type
!*/
public:
typedef typename M::type type;
typedef typename M::mem_manager_type mem_manager_type;
typedef matrix_divscal_exp ref_type;
const static long NR = M::NR;
const static long NC = M::NC;
matrix_divscal_exp (
const M& m_,
const S& s_
) :
m(m_),
s(s_)
{}
const type operator() (
long r,
long c
) const { return m(r,c)/s; }
template <typename U, long iNR, long iNC, typename mm >
bool aliases (
const matrix<U,iNR,iNC,mm>& item
) const { return m.aliases(item); }
template <typename U, long iNR, long iNC, typename mm >
bool destructively_aliases (
const matrix<U,iNR,iNC,mm>& item
) const { return m.destructively_aliases(item); }
const ref_type& ref(
) const { return *this; }
long nr (
) const { return m.nr(); }
long nc (
) const { return m.nc(); }
const M m;
const S s;
};
template <
typename EXP,
typename S
>
inline const matrix_exp<matrix_divscal_exp<matrix_exp<EXP>, S> > operator/ (
const matrix_exp<EXP>& m,
const S& s
)
{
typedef matrix_divscal_exp<matrix_exp<EXP>,S > exp;
return matrix_exp<exp>(exp(m,s));
}
// ----------------------------------------------------------------------------------------
template <
typename M,
typename S
>
class matrix_mulscal_exp
{
/*!
REQUIREMENTS ON M
- must be a matrix_exp or matrix_ref object (or
an object with a compatible interface).
REQUIREMENTS ON S
- must be some kind of scalar type
!*/
public:
typedef typename M::type type;
typedef typename M::mem_manager_type mem_manager_type;
typedef matrix_mulscal_exp ref_type;
const static long NR = M::NR;
const static long NC = M::NC;
matrix_mulscal_exp (
const M& m_,
const S& s_
) :
m(m_),
s(s_)
{}
const type operator() (
long r,
long c
) const { return m(r,c)*s; }
template <typename U, long iNR, long iNC , typename mm>
bool aliases (
const matrix<U,iNR,iNC,mm>& item
) const { return m.aliases(item); }
template <typename U,