From a1d843c01be29004d03ccd1935c83595b06dfe0d Mon Sep 17 00:00:00 2001 From: Christian Zimmermann Date: Tue, 15 Jan 2019 14:34:59 +0100 Subject: [PATCH] change operation arithmetics: distinguish between assignment and contraction (= vs +=) --- src/include/expressions.h | 5 ++- src/include/mbase_def.h | 2 +- src/include/multi_array_operation.cc.h | 48 ++++++++++++--------- src/include/multi_array_operation.h | 22 +++++++--- src/tests/op_perf_test.cc | 59 ++++++++++++++------------ src/tests/op_unit_test.cc | 14 +++--- 6 files changed, 89 insertions(+), 61 deletions(-) diff --git a/src/include/expressions.h b/src/include/expressions.h index 2d40b38..2f5a6b5 100644 --- a/src/include/expressions.h +++ b/src/include/expressions.h @@ -30,7 +30,10 @@ namespace MultiArrayTools using OX = Operation,oo...>; template - using AEXT = typename OperationMaster>::AssignmentExpr; + using AEXT = typename OperationMaster,Second,DynamicRange>::AssignmentExpr; + + template + using AEXT_P = typename OperationMaster,Second,DynamicRange>::AssignmentExpr; template class OpF, class... MAs> using AEX = AEXT>; diff --git a/src/include/mbase_def.h b/src/include/mbase_def.h index 90d432e..1f91e6a 100644 --- a/src/include/mbase_def.h +++ b/src/include/mbase_def.h @@ -33,7 +33,7 @@ namespace MultiArrayTools class OperationTemplate; // multi_array_operation.h - template + template class OperationMaster; // multi_array_operation.h diff --git a/src/include/multi_array_operation.cc.h b/src/include/multi_array_operation.cc.h index d75e6ff..25b523a 100644 --- a/src/include/multi_array_operation.cc.h +++ b/src/include/multi_array_operation.cc.h @@ -110,22 +110,22 @@ namespace MultiArrayTools * OperationMaster::AssignmentExpr * *****************************************/ - template - OperationMaster::AssignmentExpr:: + template + OperationMaster::AssignmentExpr:: AssignmentExpr(OperationMaster& m, const OpClass& sec) : mM(m), mSec(sec) {} - template - inline void OperationMaster::AssignmentExpr:: + template + inline void OperationMaster::AssignmentExpr:: operator()(size_t start, ExtType last) const { //VCHECK(mSec.template get(last)); - mM.add(start, mSec.template get(last) ); + mM.set(start, mSec.template get(last) ); } - template - typename OperationMaster::AssignmentExpr::ExtType - OperationMaster::AssignmentExpr:: + template + typename OperationMaster::AssignmentExpr::ExtType + OperationMaster::AssignmentExpr:: rootSteps(std::intptr_t iPtrNum) const { return mSec.rootSteps(iPtrNum); @@ -136,8 +136,8 @@ namespace MultiArrayTools * OperationMaster * *************************/ - template - OperationMaster:: + template + OperationMaster:: OperationMaster(MutableMultiArrayBase& ma, const OpClass& second, IndexType& index) : mSecond(second), mDataPtr(ma.data()), @@ -146,8 +146,8 @@ namespace MultiArrayTools performAssignment(0); } - template - OperationMaster:: + template + OperationMaster:: OperationMaster(T* data, const OpClass& second, IndexType& index) : mSecond(second), mDataPtr(data), @@ -156,17 +156,16 @@ namespace MultiArrayTools performAssignment(0); } - template - void OperationMaster::performAssignment(std::intptr_t blockIndexNum) + template + void OperationMaster::performAssignment(std::intptr_t blockIndexNum) { AssignmentExpr ae(*this, mSecond); // Expression to be executed within loop - const auto loop = mSecond.template loop( mIndex.ifor(1,ae) ); - // hidden Loops outside ! -> auto vectorizable + const auto loop = mIndex.ifor( 1, mSecond.loop(ae) ); loop(); // execute overall loop(s) and so internal hidden loops and so the inherited expressions } - template - inline T OperationMaster::get(size_t pos) const + template + inline T OperationMaster::get(size_t pos) const { return mDataPtr[pos]; } @@ -354,13 +353,20 @@ namespace MultiArrayTools template template - OperationMaster OperationRoot::operator=(const OpClass& in) + OperationMaster,OpClass,Ranges...> OperationRoot::operator=(const OpClass& in) { - return OperationMaster(mDataPtr, in, mIndex); + return OperationMaster,OpClass,Ranges...>(mDataPtr, in, mIndex); } template - OperationMaster,Ranges...> + template + OperationMaster,OpClass,Ranges...> OperationRoot::operator+=(const OpClass& in) + { + return OperationMaster,OpClass,Ranges...>(mDataPtr, in, mIndex); + } + + template + OperationMaster,OperationRoot,Ranges...> OperationRoot::operator=(const OperationRoot& in) { return operator= >(in); diff --git a/src/include/multi_array_operation.h b/src/include/multi_array_operation.h index 6869454..89773bb 100644 --- a/src/include/multi_array_operation.h +++ b/src/include/multi_array_operation.h @@ -88,8 +88,16 @@ namespace MultiArrayTools friend OperationClass; }; + template + struct SelfIdentity + { + static inline T apply(const T& a, const T& b) + { + return b; + } + }; - template + template class OperationMaster { public: @@ -132,8 +140,9 @@ namespace MultiArrayTools OperationMaster(T* data, const OpClass& second, IndexType& index); - inline void set(size_t pos, T val) { mDataPtr[pos] = val; } - inline void add(size_t pos, T val) { mDataPtr[pos] += val; } + inline void set(size_t pos, T val) { mDataPtr[pos] = AOp::apply(mDataPtr[pos],val); } + + //inline void add(size_t pos, T val) { mDataPtr[pos] += val; } inline T get(size_t pos) const; private: @@ -276,9 +285,12 @@ namespace MultiArrayTools OperationRoot(T* data, const IndexType& ind); template - OperationMaster operator=(const OpClass& in); + OperationMaster,OpClass,Ranges...> operator=(const OpClass& in); - OperationMaster operator=(const OperationRoot& in); + template + OperationMaster,OpClass,Ranges...> operator+=(const OpClass& in); + + OperationMaster,OperationRoot,Ranges...> operator=(const OperationRoot& in); template inline T get(ET pos) const; diff --git a/src/tests/op_perf_test.cc b/src/tests/op_perf_test.cc index f42ecfd..43495e7 100644 --- a/src/tests/op_perf_test.cc +++ b/src/tests/op_perf_test.cc @@ -126,12 +126,15 @@ namespace { { public: + typedef ClassicRF CRF; + typedef ClassicRange CR; + typedef SpinRF SRF; typedef SpinRange SR; typedef MultiRangeFactory SR8F; typedef SR8F::oType SR8; - static const size_t s = 65536; + static const size_t s = 65536*1000; OpTest_Spin() { @@ -142,6 +145,8 @@ namespace { } SRF f; sr = std::dynamic_pointer_cast(f.create()); + CRF cf(1000); + cr = std::dynamic_pointer_cast(cf.create()); } void contract(); @@ -150,13 +155,15 @@ namespace { std::vector data; std::shared_ptr sr; + std::shared_ptr cr; }; void OpTest_Spin::contract() { - MultiArray ma(sr, sr, sr, sr, sr, sr, sr, sr, data); - MultiArray res1( sr, sr ); - + MultiArray ma( cr, sr, sr, sr, sr, sr, sr, sr, sr, data); + MultiArray res1( cr, sr, sr ); + + auto ii = MAT::getIndex(cr); auto alpha = MAT::getIndex(); auto beta = MAT::getIndex(); auto gamma = MAT::getIndex(); @@ -166,14 +173,14 @@ namespace { auto mix = MAT::mkMIndex( alpha, beta, gamma ); std::clock_t begin = std::clock(); - for(size_t i = 0; i != 1000; ++i){ - res1(delta, deltap) = ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix); - } + //for(size_t i = 0; i != 1000; ++i){ + res1(ii ,delta, deltap) += ma(ii, delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix); + //} std::clock_t end = std::clock(); std::cout << "MultiArray time: " << static_cast( end - begin ) / CLOCKS_PER_SEC << std::endl; - std::vector vres(4*4); + std::vector vres(4*4*1000); for(size_t d = 0; d != 4; ++d){ for(size_t p = 0; p != 4; ++p){ const size_t tidx = d*4 + p; @@ -187,8 +194,8 @@ namespace { for(size_t c = 0; c != 4; ++c){ for(size_t d = 0; d != 4; ++d){ for(size_t p = 0; p != 4; ++p){ - const size_t tidx = d*4 + p; - const size_t sidx = d*4*4*4*4*4*4*4 + a*5*4*4*4*4*4 + b*5*4*4*4 + + c*5*4 + p; + const size_t tidx = i*4*4 + d*4 + p; + const size_t sidx = i*65536 + d*4*4*4*4*4*4*4 + a*5*4*4*4*4*4 + b*5*4*4*4 + c*5*4 + p; vres[tidx] += data[sidx]; } } @@ -198,25 +205,25 @@ namespace { } std::clock_t end2 = std::clock(); - assert( xround(res1.at(mkts(0,0))) == xround(vres[0]) ); - assert( xround(res1.at(mkts(0,1))) == xround(vres[1]) ); - assert( xround(res1.at(mkts(0,2))) == xround(vres[2]) ); - assert( xround(res1.at(mkts(0,3))) == xround(vres[3]) ); + assert( xround(res1.at(mkts(0,0,0))) == xround(vres[0]) ); + assert( xround(res1.at(mkts(0,0,1))) == xround(vres[1]) ); + assert( xround(res1.at(mkts(0,0,2))) == xround(vres[2]) ); + assert( xround(res1.at(mkts(0,0,3))) == xround(vres[3]) ); - assert( xround(res1.at(mkts(1,0))) == xround(vres[4]) ); - assert( xround(res1.at(mkts(1,1))) == xround(vres[5]) ); - assert( xround(res1.at(mkts(1,2))) == xround(vres[6]) ); - assert( xround(res1.at(mkts(1,3))) == xround(vres[7]) ); + assert( xround(res1.at(mkts(0,1,0))) == xround(vres[4]) ); + assert( xround(res1.at(mkts(0,1,1))) == xround(vres[5]) ); + assert( xround(res1.at(mkts(0,1,2))) == xround(vres[6]) ); + assert( xround(res1.at(mkts(0,1,3))) == xround(vres[7]) ); - assert( xround(res1.at(mkts(2,0))) == xround(vres[8]) ); - assert( xround(res1.at(mkts(2,1))) == xround(vres[9]) ); - assert( xround(res1.at(mkts(2,2))) == xround(vres[10]) ); - assert( xround(res1.at(mkts(2,3))) == xround(vres[11]) ); + assert( xround(res1.at(mkts(0,2,0))) == xround(vres[8]) ); + assert( xround(res1.at(mkts(0,2,1))) == xround(vres[9]) ); + assert( xround(res1.at(mkts(0,2,2))) == xround(vres[10]) ); + assert( xround(res1.at(mkts(0,2,3))) == xround(vres[11]) ); - assert( xround(res1.at(mkts(3,0))) == xround(vres[12]) ); - assert( xround(res1.at(mkts(3,1))) == xround(vres[13]) ); - assert( xround(res1.at(mkts(3,2))) == xround(vres[14]) ); - assert( xround(res1.at(mkts(3,3))) == xround(vres[15]) ); + assert( xround(res1.at(mkts(0,3,0))) == xround(vres[12]) ); + assert( xround(res1.at(mkts(0,3,1))) == xround(vres[13]) ); + assert( xround(res1.at(mkts(0,3,2))) == xround(vres[14]) ); + assert( xround(res1.at(mkts(0,3,3))) == xround(vres[15]) ); std::cout << "std::vector - for loop time: " << static_cast( end2 - begin2 ) / CLOCKS_PER_SEC << std::endl; diff --git a/src/tests/op_unit_test.cc b/src/tests/op_unit_test.cc index 6d1d66d..a0a8300 100644 --- a/src/tests/op_unit_test.cc +++ b/src/tests/op_unit_test.cc @@ -367,10 +367,10 @@ namespace { (*jj)(ii1,ii2); ///auto jj = mkMapI( std::make_shared >(), ii1, ii1 ); - res(jj) = ma1(ii1,ii2); + res(jj) += ma1(ii1,ii2); auto mult = mr->mapMultiplicity(); auto jjx = jj->outIndex(); - res2(jj) = ma1(ii1,ii2) / staticcast( mult(jjx) ); + res2(jj) += ma1(ii1,ii2) / staticcast( mult(jjx) ); MultiArray form = res.format(mpr1ptr->outRange()); MultiArray form2 = res2.format(mpr1ptr->outRange()); @@ -419,7 +419,7 @@ namespace { auto i1 = MAT::getIndex(sr1ptr); auto i2 = MAT::getIndex(sr2ptr); - res(i1) = fma(i1,i2).c(i2); + res(i1) += fma(i1,i2).c(i2); auto i = res.begin(); @@ -452,12 +452,12 @@ namespace { auto mix = MAT::mkMIndex( alpha, beta, gamma ); std::clock_t begin = std::clock(); - res1(delta, deltap) = ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix); + res1(delta, deltap) += ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix); std::clock_t end = std::clock(); std::cout << "MultiArray time: " << static_cast( end - begin ) / CLOCKS_PER_SEC << std::endl; - res2(delta, deltap) = ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(alpha).c(beta).c(gamma); + res2(delta, deltap) += ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(alpha).c(beta).c(gamma); std::vector vres(4*4); @@ -556,7 +556,7 @@ namespace { auto si = MAT::getIndex( subptr ); (*si)(i2); - res(i3,i1) = (ma2(i3,i2) - ma1(i1,i2,i3)).c(si); + res(i3,i1) += (ma2(i3,i2) - ma1(i1,i2,i3)).c(si); res2(i3,si,i1) = ma2(i3,i2) - ma1(i1,i2,i3); EXPECT_EQ( res2.size(), static_cast(8) ); @@ -628,7 +628,7 @@ namespace { auto i1 = MAT::getIndex( sr2ptr ); auto i2 = MAT::getIndex( sr4ptr ); - res(i1) = (ma1(i1) * ma2(i2)).c(i2); + res(i1) += (ma1(i1) * ma2(i2)).c(i2); EXPECT_EQ( xround( res.at('1') ), xround(2.917 * 8.870 + 2.917 * 4.790) ); EXPECT_EQ( xround( res.at('2') ), xround(9.436 * 8.870 + 9.436 * 4.790) );