change operation arithmetics: distinguish between assignment and contraction (= vs +=)
This commit is contained in:
Christian Zimmermann
parent
e04d9aa5bc
commit
a1d843c01b
6 changed files with 89 additions and 61 deletions
|
|
@ -30,7 +30,10 @@ namespace MultiArrayTools
|
|||
using OX = Operation<double,OpF<double>,oo<EC,MAs>...>;
|
||||
|
||||
template <class EC, class Second>
|
||||
using AEXT = typename OperationMaster<double,Second,DynamicRange<EC>>::AssignmentExpr;
|
||||
using AEXT = typename OperationMaster<double,SelfIdentity<double>,Second,DynamicRange<EC>>::AssignmentExpr;
|
||||
|
||||
template <class EC, class Second>
|
||||
using AEXT_P = typename OperationMaster<double,plus<double>,Second,DynamicRange<EC>>::AssignmentExpr;
|
||||
|
||||
template <class EC, template <class> class OpF, class... MAs>
|
||||
using AEX = AEXT<EC,OX<EC,OpF,MAs...>>;
|
||||
|
|
|
|||
|
|
@ -33,7 +33,7 @@ namespace MultiArrayTools
|
|||
class OperationTemplate;
|
||||
|
||||
// multi_array_operation.h
|
||||
template <typename T, class OpClass, class... Ranges>
|
||||
template <typename T, class AOp, class OpClass, class... Ranges>
|
||||
class OperationMaster;
|
||||
|
||||
// multi_array_operation.h
|
||||
|
|
|
|||
|
|
@ -110,22 +110,22 @@ namespace MultiArrayTools
|
|||
* OperationMaster::AssignmentExpr *
|
||||
*****************************************/
|
||||
|
||||
template <typename T, class OpClass, class... Ranges>
|
||||
OperationMaster<T,OpClass,Ranges...>::AssignmentExpr::
|
||||
template <typename T, class AOp, class OpClass, class... Ranges>
|
||||
OperationMaster<T,AOp,OpClass,Ranges...>::AssignmentExpr::
|
||||
AssignmentExpr(OperationMaster& m, const OpClass& sec) :
|
||||
mM(m), mSec(sec) {}
|
||||
|
||||
template <typename T, class OpClass, class... Ranges>
|
||||
inline void OperationMaster<T,OpClass,Ranges...>::AssignmentExpr::
|
||||
template <typename T, class AOp, class OpClass, class... Ranges>
|
||||
inline void OperationMaster<T,AOp,OpClass,Ranges...>::AssignmentExpr::
|
||||
operator()(size_t start, ExtType last) const
|
||||
{
|
||||
//VCHECK(mSec.template get<ExtType>(last));
|
||||
mM.add(start, mSec.template get<ExtType>(last) );
|
||||
mM.set(start, mSec.template get<ExtType>(last) );
|
||||
}
|
||||
|
||||
template <typename T, class OpClass, class... Ranges>
|
||||
typename OperationMaster<T,OpClass,Ranges...>::AssignmentExpr::ExtType
|
||||
OperationMaster<T,OpClass,Ranges...>::AssignmentExpr::
|
||||
template <typename T, class AOp, class OpClass, class... Ranges>
|
||||
typename OperationMaster<T,AOp,OpClass,Ranges...>::AssignmentExpr::ExtType
|
||||
OperationMaster<T,AOp,OpClass,Ranges...>::AssignmentExpr::
|
||||
rootSteps(std::intptr_t iPtrNum) const
|
||||
{
|
||||
return mSec.rootSteps(iPtrNum);
|
||||
|
|
@ -136,8 +136,8 @@ namespace MultiArrayTools
|
|||
* OperationMaster *
|
||||
*************************/
|
||||
|
||||
template <typename T, class OpClass, class... Ranges>
|
||||
OperationMaster<T,OpClass,Ranges...>::
|
||||
template <typename T, class AOp, class OpClass, class... Ranges>
|
||||
OperationMaster<T,AOp,OpClass,Ranges...>::
|
||||
OperationMaster(MutableMultiArrayBase<T,Ranges...>& ma, const OpClass& second,
|
||||
IndexType& index) :
|
||||
mSecond(second), mDataPtr(ma.data()),
|
||||
|
|
@ -146,8 +146,8 @@ namespace MultiArrayTools
|
|||
performAssignment(0);
|
||||
}
|
||||
|
||||
template <typename T, class OpClass, class... Ranges>
|
||||
OperationMaster<T,OpClass,Ranges...>::
|
||||
template <typename T, class AOp, class OpClass, class... Ranges>
|
||||
OperationMaster<T,AOp,OpClass,Ranges...>::
|
||||
OperationMaster(T* data, const OpClass& second,
|
||||
IndexType& index) :
|
||||
mSecond(second), mDataPtr(data),
|
||||
|
|
@ -156,17 +156,16 @@ namespace MultiArrayTools
|
|||
performAssignment(0);
|
||||
}
|
||||
|
||||
template <typename T, class OpClass, class... Ranges>
|
||||
void OperationMaster<T,OpClass,Ranges...>::performAssignment(std::intptr_t blockIndexNum)
|
||||
template <typename T, class AOp, class OpClass, class... Ranges>
|
||||
void OperationMaster<T,AOp,OpClass,Ranges...>::performAssignment(std::intptr_t blockIndexNum)
|
||||
{
|
||||
AssignmentExpr ae(*this, mSecond); // Expression to be executed within loop
|
||||
const auto loop = mSecond.template loop<decltype(mIndex.ifor(1,ae))>( 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 <typename T, class OpClass, class... Ranges>
|
||||
inline T OperationMaster<T,OpClass,Ranges...>::get(size_t pos) const
|
||||
template <typename T, class AOp, class OpClass, class... Ranges>
|
||||
inline T OperationMaster<T,AOp,OpClass,Ranges...>::get(size_t pos) const
|
||||
{
|
||||
return mDataPtr[pos];
|
||||
}
|
||||
|
|
@ -354,13 +353,20 @@ namespace MultiArrayTools
|
|||
|
||||
template <typename T, class... Ranges>
|
||||
template <class OpClass>
|
||||
OperationMaster<T,OpClass,Ranges...> OperationRoot<T,Ranges...>::operator=(const OpClass& in)
|
||||
OperationMaster<T,SelfIdentity<T>,OpClass,Ranges...> OperationRoot<T,Ranges...>::operator=(const OpClass& in)
|
||||
{
|
||||
return OperationMaster<T,OpClass,Ranges...>(mDataPtr, in, mIndex);
|
||||
return OperationMaster<T,SelfIdentity<T>,OpClass,Ranges...>(mDataPtr, in, mIndex);
|
||||
}
|
||||
|
||||
template <typename T, class... Ranges>
|
||||
OperationMaster<T,OperationRoot<T,Ranges...>,Ranges...>
|
||||
template <class OpClass>
|
||||
OperationMaster<T,plus<T>,OpClass,Ranges...> OperationRoot<T,Ranges...>::operator+=(const OpClass& in)
|
||||
{
|
||||
return OperationMaster<T,plus<T>,OpClass,Ranges...>(mDataPtr, in, mIndex);
|
||||
}
|
||||
|
||||
template <typename T, class... Ranges>
|
||||
OperationMaster<T,SelfIdentity<T>,OperationRoot<T,Ranges...>,Ranges...>
|
||||
OperationRoot<T,Ranges...>::operator=(const OperationRoot<T,Ranges...>& in)
|
||||
{
|
||||
return operator=<OperationRoot<T,Ranges...> >(in);
|
||||
|
|
|
|||
|
|
@ -88,8 +88,16 @@ namespace MultiArrayTools
|
|||
friend OperationClass;
|
||||
};
|
||||
|
||||
template <typename T>
|
||||
struct SelfIdentity
|
||||
{
|
||||
static inline T apply(const T& a, const T& b)
|
||||
{
|
||||
return b;
|
||||
}
|
||||
};
|
||||
|
||||
template <typename T, class OpClass, class... Ranges>
|
||||
template <typename T, class AOp, class OpClass, class... Ranges>
|
||||
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 <class OpClass>
|
||||
OperationMaster<T,OpClass,Ranges...> operator=(const OpClass& in);
|
||||
OperationMaster<T,SelfIdentity<T>,OpClass,Ranges...> operator=(const OpClass& in);
|
||||
|
||||
OperationMaster<T,OperationRoot,Ranges...> operator=(const OperationRoot& in);
|
||||
template <class OpClass>
|
||||
OperationMaster<T,plus<T>,OpClass,Ranges...> operator+=(const OpClass& in);
|
||||
|
||||
OperationMaster<T,SelfIdentity<T>,OperationRoot,Ranges...> operator=(const OperationRoot& in);
|
||||
|
||||
template <class ET>
|
||||
inline T get(ET pos) const;
|
||||
|
|
|
|||
|
|
@ -126,12 +126,15 @@ namespace {
|
|||
{
|
||||
public:
|
||||
|
||||
typedef ClassicRF CRF;
|
||||
typedef ClassicRange CR;
|
||||
|
||||
typedef SpinRF SRF;
|
||||
typedef SpinRange SR;
|
||||
typedef MultiRangeFactory<SR,SR,SR,SR,SR,SR,SR,SR> 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<SR>(f.create());
|
||||
CRF cf(1000);
|
||||
cr = std::dynamic_pointer_cast<CR>(cf.create());
|
||||
}
|
||||
|
||||
void contract();
|
||||
|
|
@ -150,13 +155,15 @@ namespace {
|
|||
|
||||
std::vector<double> data;
|
||||
std::shared_ptr<SR> sr;
|
||||
std::shared_ptr<CR> cr;
|
||||
};
|
||||
|
||||
void OpTest_Spin::contract()
|
||||
{
|
||||
MultiArray<double,SR,SR,SR,SR,SR,SR,SR,SR> ma(sr, sr, sr, sr, sr, sr, sr, sr, data);
|
||||
MultiArray<double,SR,SR> res1( sr, sr );
|
||||
MultiArray<double,CR,SR,SR,SR,SR,SR,SR,SR,SR> ma( cr, sr, sr, sr, sr, sr, sr, sr, sr, data);
|
||||
MultiArray<double,CR,SR,SR> res1( cr, sr, sr );
|
||||
|
||||
auto ii = MAT::getIndex<CR>(cr);
|
||||
auto alpha = MAT::getIndex<SR>();
|
||||
auto beta = MAT::getIndex<SR>();
|
||||
auto gamma = MAT::getIndex<SR>();
|
||||
|
|
@ -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<double>( end - begin ) / CLOCKS_PER_SEC
|
||||
<< std::endl;
|
||||
|
||||
std::vector<double> vres(4*4);
|
||||
std::vector<double> 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<double>( end2 - begin2 ) / CLOCKS_PER_SEC
|
||||
<< std::endl;
|
||||
|
|
|
|||
|
|
@ -367,10 +367,10 @@ namespace {
|
|||
(*jj)(ii1,ii2);
|
||||
///auto jj = mkMapI( std::make_shared<plus<size_t> >(), 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<double>( mult(jjx) );
|
||||
res2(jj) += ma1(ii1,ii2) / staticcast<double>( mult(jjx) );
|
||||
|
||||
MultiArray<double,TRange> form = res.format(mpr1ptr->outRange());
|
||||
MultiArray<double,TRange> 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<double>( 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<double> 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<size_t>(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) );
|
||||
|
|
|
|||
Loading…
Reference in a new issue