Bugzilla – Attachment 400 Details for
Bug 99
Expression evaluator
Home
|
New
|
Browse
|
Search
|
[?]
|
Reports
|
Requests
|
Help
|
Log In
[x]
|
Forgot Password
Login:
[x]
This bugzilla service is closed. All entries have been migrated to
https://gitlab.com/libeigen/eigen
Proof of concept V3
proof_of_concept_v3.cpp (text/x-c++src), 18.93 KB, created by
Gael Guennebaud
on 2013-11-25 13:39:10 UTC
(
hide
)
Description:
Proof of concept V3
Filename:
MIME Type:
Creator:
Gael Guennebaud
Created:
2013-11-25 13:39:10 UTC
Size:
18.93 KB
patch
obsolete
>#include <iostream> > >#define TRACK std::cout << __LINE__ << "\n"; > >//----------------------------- >// --- Forward declarations --- >//----------------------------- > >struct Matrix; >struct SparseMatrix; >template<typename A> struct Noalias; >template<typename A, typename B> struct Sum; >template<typename A, typename B> struct Product; >template<typename A, typename B> struct LazyProduct; > >template<typename T, int StorageKind = T::StorageKind> struct evaluator; > >template<int EvalKindDst, int EvalKindSrc> struct assignment_kind; > >template< typename A, typename B, typename Op, > int AssignmentKind = assignment_kind< evaluator<A>::Kind , evaluator<B>::Kind >::Kind> >struct Assignment; > >enum {DenseStorage=1, SparseStorage=2, DenseAssignment=4, SparseAssignment=8}; > >template<typename Dst, typename Src, typename Func> >void assign(Dst& dst, const Src& src, const Func& func); > > >//--------------------------------- >// --- Functors for assignments --- >//--------------------------------- > >template<typename Scalar> struct assign_default { void assign(Scalar* a, const Scalar& b) const { *a = b; } std::string name() const { return " = "; } }; >template<typename Scalar> struct assign_add { void assign(Scalar* a, const Scalar& b) const { *a += b; } std::string name() const { return " += "; } }; > > > > >//---------------------------------------- >// --- The expression template classes --- >//---------------------------------------- > >/// The base class (like DenseBase/MatrixBase) >template<typename Derived> struct Base >{ > const Derived& derived() const { return *static_cast<const Derived*>(this); } > Derived& derived() { return *static_cast<Derived*>(this); } > > template<typename Other> Sum<Derived,Other> operator+(const Base<Other>& other) const > { return Sum<Derived,Other>(derived(), other.derived()); } > > template<typename Other> Product<Derived,Other> operator*(const Base<Other>& other) const > { return Product<Derived,Other>(derived(), other.derived()); } > > template<typename Other> LazyProduct<Derived,Other> operator%(const Base<Other>& other) const > { return LazyProduct<Derived,Other>(derived(), other.derived()); } > > Noalias<Derived> noalias() { return Noalias<Derived>(derived()); } > > template<typename Other> Derived& operator=(const Base<Other>& other) > { > assign(derived(), other.derived(), assign_default<int>()); > return derived(); > } > > template<typename Other> Derived& operator+=(const Base<Other>& other) > { > assign(derived(), other.derived(), assign_add<int>()); > return derived(); > } >}; > >template<typename Derived> struct DenseBase : Base<Derived> >{ > enum { > StorageKind = DenseStorage > }; > using Base<Derived>::operator=; >}; > >template<typename Derived> struct SparseBase : Base<Derived> >{ > enum { > StorageKind = SparseStorage > }; > using Base<Derived>::operator=; >}; > >template<typename Derived, int StorageKind> >struct base_class_dispacher; > >template<typename Derived> >struct base_class_dispacher<Derived,DenseStorage> >{ typedef DenseBase<Derived> type; }; > >template<typename Derived> >struct base_class_dispacher<Derived,SparseStorage> >{ typedef SparseBase<Derived> type; }; > >/// A 2x2 Matrix class with storage >struct Matrix : DenseBase<Matrix> >{ > typedef DenseBase<Matrix> Base; > using Base::operator=; > Matrix(const std::string& name) : m_name(name) { data = new int[4]; std::cout << "- Create matrix " << m_name << "\n"; } > ~Matrix() { delete[] data; std::cout << "- Delete matrix " << m_name << "\n"; } > Matrix(const Matrix& m) { std::cerr << "- COPY matrix\n"; } > > int coeff(int i) const { return data[i]; } > int& coeffRef(int i) { return data[i]; } > > int size() const { return 4; } > void setZero() { > std::cout << "- Set " << m_name << " to zero\n"; > for(int i=0;i<4;++i) > data[i] = 0; > } > > std::string name() const { return m_name; } > > typedef const Matrix& Nested; > int* data; > std::string m_name; >}; > >/// A 2x2 sparse matrix class with storage >struct SparseMatrix : SparseBase<SparseMatrix> >{ > enum { > FOO = SparseStorage > }; > typedef SparseBase<SparseMatrix> Base; > using Base::operator=; > SparseMatrix(const std::string& name) : m_name(name) { > indices = new int[4]; > data = new int[4]; > std::cout << "- Create sparse matrix " << m_name << "\n"; > } > ~SparseMatrix() { delete[] data; delete[] indices; std::cout << "- Delete sparse matrix " << m_name << "\n"; } > SparseMatrix(const SparseMatrix& m) { std::cerr << "- COPY sparse matrix\n"; } > > int size() const { return 4; } > > std::string name() const { return m_name; } > > typedef const SparseMatrix& Nested; > int* indices; > int* data; > std::string m_name; >}; > >/// Expression of A + B >template<typename A, typename B> struct Sum > : base_class_dispacher<Sum<A,B>, A::StorageKind>::type >{ > typedef const Sum Nested; > explicit Sum(const A& a, const B& b) : m_a(a), m_b(b) {} > > std::string name() const { return std::string("(") + m_a.name() + " + " + m_b.name() + ")"; } > > typename A::Nested m_a; > typename B::Nested m_b; >}; > >/// Expression of A * B >template<typename A, typename B> struct Product > : base_class_dispacher<Product<A,B>, A::StorageKind>::type >{ >// enum { StorageKind = A::StorageKind }; > typedef const Product Nested; > explicit Product(const A& a, const B& b) : m_a(a), m_b(b) {} > typename A::Nested m_a; > typename B::Nested m_b; > > std::string name() const { return std::string("(") + m_a.name() + " * " + m_b.name() + ")"; } >}; > >/// Expression of A * B with lazy evaluation >template<typename A, typename B> struct LazyProduct > : base_class_dispacher<LazyProduct<A,B>, A::StorageKind>::type >{ >// enum { StorageKind = A::StorageKind }; > typedef const LazyProduct Nested; > explicit LazyProduct(const A& a, const B& b) : m_a(a), m_b(b) {} > typename A::Nested m_a; > typename B::Nested m_b; > > std::string name() const { return std::string("(") + m_a.name() + " % " + m_b.name() + ")"; } >}; > >/// Expression of A.noalias() >template<typename A> struct Noalias > : base_class_dispacher<Noalias<A>, A::StorageKind >::type >{ > typedef typename base_class_dispacher<Noalias, A::StorageKind >::type BaseType; > using BaseType::operator=; > > typedef Noalias Nested; > explicit Noalias(A& a) : m_a(a) {} > int size() const { return m_a.size(); } > > std::string name() const { return std::string("noalias(") + m_a.name() + ")"; } > > A &m_a; >}; > >/// Expression of explicit evaluation to a temporary >template<typename A> struct EvalToTemp > : base_class_dispacher<EvalToTemp<A>, A::StorageKind >::type >{ > typedef typename base_class_dispacher<EvalToTemp, A::StorageKind >::type BaseType; > using BaseType::operator=; > > typedef EvalToTemp Nested; > EvalToTemp(const A& a) : m_a(a) {} > int size() const { return m_a.size(); } > > std::string name() const { return std::string("eval(") + m_a.name() + ")"; } > > const A &m_a; >}; > > > > > > >//------------------------- >// --- Dense Evaluators --- >//------------------------- > >/// Default evaluator assuming the expression is it's own evaluator > >template<> >struct evaluator<Matrix,DenseStorage> >{ > enum { > MayAlias = 0, > Kind = DenseStorage > }; > typedef evaluator type; > evaluator(const Matrix& x) : xpr(x) {/*std::cout << __LINE__ << "\n";*/} > int coeff(int i) const { return xpr.coeff(i); } > int& coeffRef(int i) { return const_cast<Matrix&>(xpr).coeffRef(i); } > const Matrix& xpr; > > std::string name() const { return xpr.name(); } >}; > >/// evaluator of A+B >template<typename A, typename B> >struct evaluator<Sum<A,B>, DenseStorage> >{ > enum { > MayAlias = evaluator<A>::MayAlias || evaluator<B>::MayAlias, > Kind = DenseStorage > }; > typedef evaluator type; > evaluator(const Sum<A,B>& plus) : m_lhs(plus.m_a), m_rhs(plus.m_b) {} > int coeff(int i) const { return m_lhs.coeff(i) + m_rhs.coeff(i); } > typename evaluator<A>::type m_lhs; > typename evaluator<B>::type m_rhs; > > std::string name() const { return std::string("(") + m_lhs.name() + " + " + m_rhs.name() + ")"; } >}; > >/// evaluator of A.noalias() >template<typename A, int StorageKind> >struct evaluator<Noalias<A>, StorageKind> : evaluator<A, StorageKind> >{ > enum { > MayAlias = 0, > Kind = StorageKind > }; > typedef evaluator type; > evaluator(Noalias<A>& a) : evaluator<A>(a.m_a) {} > > std::string name() const { return std::string("Noalias(") + evaluator<A>::name() + ")"; } >}; > >/// evaluator of A%B >template<typename A, typename B> >struct evaluator<LazyProduct<A,B>, DenseStorage> >{ > enum { > MayAlias = 1, > Kind = LazyProduct<A,B>::StorageKind > }; > typedef evaluator type; > evaluator(const LazyProduct<A,B>& prod) : m_lhs(prod.m_a), m_rhs(prod.m_b) {} > int coeff(int index) const { > int i = index%2; > int j = index/2; > return m_lhs.coeff(i+0*2)*m_rhs.coeff(0+2*j) + m_lhs.coeff(i+1*2)*m_rhs.coeff(1+2*j); > } > typename evaluator<A>::type m_lhs; > typename evaluator<B>::type m_rhs; > > std::string name() const { return std::string("(") + m_lhs.name() + " % " + m_rhs.name() + ")"; } >}; > >/// low level product implementation (c += a*b) - like blas's gemm >void product_impl_low_level(int* c, const int* a, const int* b) >{ > for(int i=0;i<2;++i) > for(int j=0;j<2;++j) > for(int k=0;k<2;++k) > c[i+j*2] += a[i+k*2] * b[k+j*2]; >} > >/// special product evaluator (in Eigen it is currently equivalent to blas_traits) >template<typename T> >struct product_evaluator >{ > product_evaluator(const T& x) > : tmp("product_arg_temporary_of_\"" + x.name() + "\"") > { > tmp = x; > } > const int* data() const { return tmp.data; } > std::string name() const { return tmp.m_name; } > Matrix tmp; >}; > >template<> >struct product_evaluator<Matrix> >{ > product_evaluator(const Matrix& x) > : m_mat(x) > {} > const int* data() const { return m_mat.data; } > std::string name() const { return m_mat.m_name; } > const Matrix& m_mat; >}; > >/// performs dest += a * b; >/// A general implementation must be available outisde the evaluator because it is needed by both >/// evaluator<Product> and Assignment<..., Product<...> > for in-place product evaluation >template<typename Dest, typename A, typename B> >void product_impl(Dest& dest, const Product<A,B>& prod) >{ > product_evaluator<A> lhs(prod.m_a); > product_evaluator<B> rhs(prod.m_b); > std::cout << "- performs: " << dest.m_name << " += " << lhs.name() << " * " << rhs.name() + " in place\n"; > product_impl_low_level(dest.data, lhs.data(), rhs.data()); >}; > >/// Matrix product evaluator, products are complex and must be evaluated into a temporary >/// unless the evaluator is by-passed by specializations of Assignment<> >template<typename A, typename B> >struct evaluator<Product<A,B>, DenseStorage > : evaluator<Matrix, DenseStorage> >{ > enum { > MayAlias = 1, > Kind = DenseStorage > }; > > typedef evaluator<Matrix> Base; > typedef evaluator type; > evaluator(const Product<A,B>& prod) : Base(res), res(std::string("product_res_of_\"") + prod.name() + "\"") > { > res.setZero(); > product_impl(res, prod); > } > Matrix res; >}; > >template<typename A> >struct evaluator<EvalToTemp<A>, DenseStorage> : evaluator<Matrix, DenseStorage> >{ > enum { > MayAlias = 0 > }; > typedef evaluator<Matrix> Base; > typedef evaluator type; > evaluator(const EvalToTemp<A>& a) : Base(res), res(std::string("temp_of_\"") + a.m_a.name() + "\"") > { > res.noalias() = a.m_a; > } > Matrix res; >}; > >// Generic assignment loop for dense evaluators. >// In Eigen it is split in many more variants to deal with unrolling and vectorization >template<typename A, typename B, typename Op> >void dense_assignment_loop(A& a, const B& b, const Op& op) >{ > typename evaluator<A>::type ea(a); > typename evaluator<B>::type eb(b); > std::cout << "- dense assignment: " << a.name() << op.name() << b.name() << "\n"; > std::cout << " as: " << ea.name() << op.name() << eb.name() << "\n"; > for(int i=0; i<a.size(); ++i) > op.assign(&ea.coeffRef(i), eb.coeff(i)); >} > > > >//-------------------------- >// --- Sparse Evaluators --- >//-------------------------- > >/// evaluator of a sparse matrix >template<> >struct evaluator<SparseMatrix,SparseStorage> >{ > enum { > MayAlias = 0, > Kind = SparseStorage > }; > typedef evaluator type; > evaluator(const SparseMatrix& x) : xpr(x) {} > const SparseMatrix& xpr; > > struct iterator { /* TODO */ }; > > std::string name() const { return xpr.name(); } >}; > > >/// evaluator of A+B >template<typename A, typename B> >struct evaluator<Sum<A,B>, SparseStorage> >{ > enum { > MayAlias = 0, // sparse expressions never alias > Kind = SparseStorage > }; > typedef evaluator type; > evaluator(const Sum<A,B>& plus) : m_lhs(plus.m_a), m_rhs(plus.m_b) {} > typename evaluator<A>::type m_lhs; > typename evaluator<B>::type m_rhs; > > struct iterator { /* TODO */ }; > > std::string name() const { return std::string("(") + m_lhs.name() + " + " + m_rhs.name() + ")"; } >}; > >/* ... TODO ... */ > >// Generic assignment loop for sparse evaluators. >template<typename A, typename B> >void sparse_assignment_loop(A& a, const B& b) >{ > typename evaluator<A>::type ea(a); > typename evaluator<B>::type eb(b); > std::cout << "- sparse assignment: " << a.name() << " = " << b.name() << "\n"; > std::cout << " as: " << ea.name() << " <- " << eb.name() << "\n"; > > // TODO >} > > > > > > > > >//------------------- >// --- Assignment --- >//------------------- > >// The only purpose of this assign() function is to deal with noalias/MayAlias >// so that the main Assignment do not has to bother about these annoying details > >template<typename Dst, typename Src, typename Func> >void assign(Dst& dst, const Src& src, const Func& func) >{ > typedef typename std::conditional<evaluator<Src>::MayAlias==1, EvalToTemp<Src>, Src>::type ActualSrc; > Assignment<Dst,ActualSrc,Func>::run(dst, src, func); >} > >// by-pass MayAlias >template<typename Dst, typename Src, typename Func> >void assign(Noalias<Dst>& dst, const Src& src, const Func& func) >{ > Assignment<Dst,Src,Func>::run(dst.m_a, src, func); >} > > > > > > >// The assignment_kind structure permits to select the appropriate assignment method depending on the kind of evaluators. > >template<> >struct assignment_kind<DenseStorage,DenseStorage> { > enum { Kind = DenseAssignment }; >}; > >template<> >struct assignment_kind<SparseStorage,SparseStorage> { > enum { Kind = SparseAssignment }; >}; > > > >// In Assignement, aliasing and noalias have already been resolved, >// so we assume no aliasing: >// - The type A cannot be a Noalias<>. >// - evaluator<B>::MayAlias has to be ignored. > > >template<typename A, typename B, typename Op> >struct Assignment<A,B,Op,DenseAssignment> >{ > static void run(A& a, const B& b, const Op& op) > { > dense_assignment_loop(a, b, op); > } >}; > >template<typename A, typename B> >struct Assignment<A,B,assign_default<int>,SparseAssignment> >{ > static void run(A& a, const B& b, const assign_default<int>& ) > { > sparse_assignment_loop(a, b); > } >}; > >template<typename A, typename B> >struct Assignment<A,B,assign_add<int>,SparseAssignment> >{ > static void run(A& a, const B& b, const assign_add<int>& ) > { > sparse_assignment_loop(a, a+b); > } >}; > > > >// This overload should never happens that is why it is empty (calling it will raise a compilation error) >template<typename A, typename B, typename Op, int AssignmentKind> >struct Assignment<Noalias<A>, B, Op, AssignmentKind> {}; > > >// specialization for a = b * c (dense) >template<typename A, typename ProdLhs, typename ProdRhs> >struct Assignment<A,Product<ProdLhs,ProdRhs>,assign_default<int>, DenseAssignment> >{ > static void run(A& a, const Product<ProdLhs,ProdRhs>& b, const assign_default<int>& ) > { > a.setZero(); > product_impl(a, b); > } >}; > >// specialization for a += b * c (dense) >template<typename A, typename ProdLhs, typename ProdRhs> >struct Assignment<A,Product<ProdLhs,ProdRhs>,assign_add<int>, DenseAssignment> >{ > static void run(A& a, const Product<ProdLhs,ProdRhs>& b, const assign_add<int>& ) > { > product_impl(a, b); > } >}; > >// specialization for a = d + b * c (dense) >template<typename A, typename D, typename ProdLhs, typename ProdRhs, typename Op> >struct Assignment<A,Sum<D,Product<ProdLhs,ProdRhs> >, Op, DenseAssignment> >{ > static void run(A& a, const Sum<D,Product<ProdLhs,ProdRhs> >& b, const Op& op) > { > Assignment<A,D,Op>::run(a,b.m_a,op); > Assignment<A,Product<ProdLhs,ProdRhs>, assign_add<int> >::run(a, b.m_b, assign_add<int>() ); > } >}; > > > > > >//------------------------ >// --- Tree Optimizers --- >//------------------------ > >/// Default tree optimizer that copies the expression by value >template<typename Xpr> struct TreeOpt >{ > typedef Xpr NewXpr; > static NewXpr make(const Xpr& xpr) { return xpr; } >}; > >/// Tree optimizer for Matrix that copies by reference >template<> struct TreeOpt<Matrix> >{ > typedef Matrix NewXpr; > static const Matrix& make(const Matrix& xpr) { return xpr; } >}; > >// Needed to forward the optimizer to the children >template<typename A, typename B> struct TreeOpt<Sum<A,B> > >{ > typedef Sum<A,B> Xpr; > typedef typename TreeOpt<A>::NewXpr NewA; > typedef typename TreeOpt<B>::NewXpr NewB; > typedef Sum<NewA,NewB> NewXpr; > static NewXpr make(const Xpr& xpr) { return TreeOpt<A>::make(xpr.m_a) + TreeOpt<B>::make(xpr.m_b); } >}; > >// catch A * B + Y and builds Y' + A' * B' >template<typename A, typename B, typename Y> >struct TreeOpt<Sum<Product<A,B>,Y> > >{ > typedef Sum<Product<A,B>, Y> Xpr; > typedef typename TreeOpt<Y>::NewXpr NewY; > typedef Sum<NewY,Product<A,B> > NewXpr; > static NewXpr make(const Xpr& xpr) { return TreeOpt<Y>::make(xpr.m_b) + xpr.m_a; } >}; > >// catch X + A * B + Y and builds (X' + Y') + (A' * B') >template<typename A, typename B, typename X, typename Y> >struct TreeOpt<Sum<Sum<X,Product<A,B> >,Y> > >{ > typedef Sum<Sum<X,Product<A,B> >,Y> Xpr; > typedef typename TreeOpt<X>::NewXpr NewX; > typedef typename TreeOpt<Y>::NewXpr NewY; > typedef Sum<Sum<NewX,NewY>,Product<A,B> > NewXpr; > static NewXpr make(const Xpr& xpr) { return (TreeOpt<X>::make(xpr.m_a.m_a) + TreeOpt<Y>::make(xpr.m_b)) + xpr.m_a.m_b; } >}; > > > > >int main() >{ > Matrix a("a"), b("b"), c("c"), d("d"); > for (int i=0;i<4;++i) > { > a.data[i] = i+10; > b.data[i] = i+20; > c.data[i] = i+30; > d.data[i] = i+40; > } > >// auto xpr = a + c*d + (a*(d+d)) + c; >// typedef __typeof(xpr) Xpr; >// >// std::cout << "init version:"; >// std::cout << " " << xpr.name() << "\n"; >// >// >// auto xpr1 = TreeOpt<Xpr>::make(xpr); >// typedef __typeof(xpr1) Xpr1; >// std::cout << "\n optimized version 1:"; >// std::cout << " " << xpr1.name() << "\n"; >// >// auto xpr2 = TreeOpt<Xpr1>::make(xpr1); >// typedef __typeof(xpr2) Xpr2; >// std::cout << "\n optimized version 2:"; >// std::cout << " " << xpr2.name() << "\n"; > > Matrix r("r"); > > std::cout << "\nEval r = a + b*c;\n"; > r = a + b*c; > > std::cout << "\nEval r.noalias() += b*c;\n"; > r.noalias() += b*c; > > std::cout << "\nEval r = b*c;\n"; > r = b*c; > > std::cout << "\nEval r.noalias() = b*c;\n"; > r.noalias() = b*c; > > std::cout << "\nEval r = b%c;\n"; > r = b%c; > > std::cout << "\nEval r.noalias() = b%c;\n"; > r.noalias() = b%c; > > std::cout << "\nEval r.noalias() += (a+b+c+d) + b*c;\n"; > r.noalias() += (a+b+c+d) + b*c; > > std::cout << "\n"; > > > SparseMatrix sa("sa"), sb("sb"), sc("sc"); > > std::cout << "\nEval sa = sb + sc;\n"; > sa = sb + sc; > > std::cout << "\nEval sa += sb;\n"; > sa += sb; > > return 0; >}
You cannot view the attachment while viewing its details because your browser does not support IFRAMEs.
View the attachment on a separate page
.
View Attachment As Raw
Actions:
View
Attachments on
bug 99
:
28
|
30
|
143
|
156
|
165
|
295
| 400