portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-12 07:00:45 +01:00
Code Releases Activity

Tested smeared RHMC Wilson1p1, accepting

Browse Source
This commit is contained in:
Guido Cossu 2016-07-07 11:49:36 +01:00
parent e87182cf98
commit ffb8b3116c
4 changed files with 754 additions and 694 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

565
lib/lattice/Lattice_ET.h
View File

@ -1,73 +1,74 @@
/************************************************************************************* /*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/lattice/Lattice_ET.h Source file: ./lib/lattice/Lattice_ET.h
Copyright (C) 2015 Copyright (C) 2015
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk> Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
Author: Peter Boyle <paboyle@ph.ed.ac.uk> Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: neo <cossu@post.kek.jp> Author: neo <cossu@post.kek.jp>
This program is free software; you can redistribute it and/or modify This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or the Free Software Foundation; either version 2 of the License, or
(at your option) any later version. (at your option) any later version.
This program is distributed in the hope that it will be useful, This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details. GNU General Public License for more details.
You should have received a copy of the GNU General Public License along You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc., with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory See the full license in the file "LICENSE" in the top level distribution
*************************************************************************************/ directory
/* END LEGAL */ *************************************************************************************/
/* END LEGAL */
#ifndef GRID_LATTICE_ET_H #ifndef GRID_LATTICE_ET_H
#define GRID_LATTICE_ET_H #define GRID_LATTICE_ET_H
#include <iostream> #include <iostream>
#include <vector>
#include <tuple> #include <tuple>
#include <typeinfo> #include <typeinfo>
#include <vector>
namespace Grid { namespace Grid {
//////////////////////////////////////////////////// ////////////////////////////////////////////////////
// Predicated where support // Predicated where support
//////////////////////////////////////////////////// ////////////////////////////////////////////////////
template<class iobj,class vobj,class robj> template <class iobj, class vobj, class robj>
inline vobj predicatedWhere(const iobj &predicate,const vobj &iftrue,const robj &iffalse) { inline vobj predicatedWhere(const iobj &predicate, const vobj &iftrue,
const robj &iffalse) {
typename std::remove_const<vobj>::type ret;
typename std::remove_const<vobj>::type ret; typedef typename vobj::scalar_object scalar_object;
typedef typename vobj::scalar_type scalar_type;
typedef typename vobj::vector_type vector_type;
typedef typename vobj::scalar_object scalar_object; const int Nsimd = vobj::vector_type::Nsimd();
typedef typename vobj::scalar_type scalar_type; const int words = sizeof(vobj) / sizeof(vector_type);
typedef typename vobj::vector_type vector_type;
const int Nsimd = vobj::vector_type::Nsimd(); std::vector<Integer> mask(Nsimd);
const int words = sizeof(vobj)/sizeof(vector_type); std::vector<scalar_object> truevals(Nsimd);
std::vector<scalar_object> falsevals(Nsimd);
std::vector<Integer> mask(Nsimd); extract(iftrue, truevals);
std::vector<scalar_object> truevals (Nsimd); extract(iffalse, falsevals);
std::vector<scalar_object> falsevals(Nsimd); extract<vInteger, Integer>(TensorRemove(predicate), mask);
extract(iftrue ,truevals); for (int s = 0; s < Nsimd; s++) {
extract(iffalse ,falsevals); if (mask[s]) falsevals[s] = truevals[s];
extract<vInteger,Integer>(TensorRemove(predicate),mask);
for(int s=0;s<Nsimd;s++){
if (mask[s]) falsevals[s]=truevals[s];
}
merge(ret,falsevals);
return ret;
} }
merge(ret, falsevals);
return ret;
}
//////////////////////////////////////////// ////////////////////////////////////////////
// recursive evaluation of expressions; Could // recursive evaluation of expressions; Could
// switch to generic approach with variadics, a la // switch to generic approach with variadics, a la
@ -75,311 +76,342 @@ namespace Grid {
// from tuple is hideous; C++14 introduces std::make_index_sequence for this // from tuple is hideous; C++14 introduces std::make_index_sequence for this
//////////////////////////////////////////// ////////////////////////////////////////////
// leaf eval of lattice ; should enable if protect using traits
//leaf eval of lattice ; should enable if protect using traits template <typename T>
using is_lattice = std::is_base_of<LatticeBase, T>;
template <typename T> using is_lattice = std::is_base_of<LatticeBase,T >; template <typename T>
using is_lattice_expr = std::is_base_of<LatticeExpressionBase, T>;
template <typename T> using is_lattice_expr = std::is_base_of<LatticeExpressionBase,T >; template <class sobj>
inline sobj eval(const unsigned int ss, const sobj &arg) {
template<class sobj>
inline sobj eval(const unsigned int ss, const sobj &arg)
{
return arg; return arg;
} }
template<class lobj> template <class lobj>
inline const lobj &eval(const unsigned int ss, const Lattice<lobj> &arg) inline const lobj &eval(const unsigned int ss, const Lattice<lobj> &arg) {
{ return arg._odata[ss];
return arg._odata[ss];
} }
// handle nodes in syntax tree // handle nodes in syntax tree
template <typename Op, typename T1> template <typename Op, typename T1>
auto inline eval(const unsigned int ss, const LatticeUnaryExpression<Op,T1 > &expr) // eval one operand auto inline eval(
-> decltype(expr.first.func(eval(ss,std::get<0>(expr.second)))) const unsigned int ss,
{ const LatticeUnaryExpression<Op, T1> &expr) // eval one operand
return expr.first.func(eval(ss,std::get<0>(expr.second))); -> decltype(expr.first.func(eval(ss, std::get<0>(expr.second)))) {
return expr.first.func(eval(ss, std::get<0>(expr.second)));
} }
template <typename Op, typename T1, typename T2> template <typename Op, typename T1, typename T2>
auto inline eval(const unsigned int ss, const LatticeBinaryExpression<Op,T1,T2> &expr) // eval two operands auto inline eval(
-> decltype(expr.first.func(eval(ss,std::get<0>(expr.second)),eval(ss,std::get<1>(expr.second)))) const unsigned int ss,
{ const LatticeBinaryExpression<Op, T1, T2> &expr) // eval two operands
return expr.first.func(eval(ss,std::get<0>(expr.second)),eval(ss,std::get<1>(expr.second))); -> decltype(expr.first.func(eval(ss, std::get<0>(expr.second)),
eval(ss, std::get<1>(expr.second)))) {
return expr.first.func(eval(ss, std::get<0>(expr.second)),
eval(ss, std::get<1>(expr.second)));
} }
template <typename Op, typename T1, typename T2, typename T3> template <typename Op, typename T1, typename T2, typename T3>
auto inline eval(const unsigned int ss, const LatticeTrinaryExpression<Op,T1,T2,T3 > &expr) // eval three operands auto inline eval(const unsigned int ss,
-> decltype(expr.first.func(eval(ss,std::get<0>(expr.second)),eval(ss,std::get<1>(expr.second)),eval(ss,std::get<2>(expr.second)))) const LatticeTrinaryExpression<Op, T1, T2, T3>
{ &expr) // eval three operands
return expr.first.func(eval(ss,std::get<0>(expr.second)),eval(ss,std::get<1>(expr.second)),eval(ss,std::get<2>(expr.second)) ); -> decltype(expr.first.func(eval(ss, std::get<0>(expr.second)),
eval(ss, std::get<1>(expr.second)),
eval(ss, std::get<2>(expr.second)))) {
return expr.first.func(eval(ss, std::get<0>(expr.second)),
eval(ss, std::get<1>(expr.second)),
eval(ss, std::get<2>(expr.second)));
} }
////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////
// Obtain the grid from an expression, ensuring conformable. This must follow a tree recursion // Obtain the grid from an expression, ensuring conformable. This must follow a
// tree recursion
////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////
template<class T1, typename std::enable_if<is_lattice<T1>::value, T1>::type * =nullptr > template <class T1,
inline void GridFromExpression(GridBase * &grid,const T1& lat) // Lattice leaf typename std::enable_if<is_lattice<T1>::value, T1>::type * = nullptr>
{ inline void GridFromExpression(GridBase *&grid, const T1 &lat) // Lattice leaf
if ( grid ) {
conformable(grid,lat._grid);
}
grid=lat._grid;
}
template<class T1,typename std::enable_if<!is_lattice<T1>::value, T1>::type * = nullptr >
inline void GridFromExpression(GridBase * &grid,const T1& notlat) // non-lattice leaf
{ {
if (grid) {
conformable(grid, lat._grid);
}
grid = lat._grid;
} }
template <class T1,
typename std::enable_if<!is_lattice<T1>::value, T1>::type * = nullptr>
inline void GridFromExpression(GridBase *&grid,
const T1 &notlat) // non-lattice leaf
{}
template <typename Op, typename T1> template <typename Op, typename T1>
inline void GridFromExpression(GridBase * &grid,const LatticeUnaryExpression<Op,T1 > &expr) inline void GridFromExpression(GridBase *&grid,
{ const LatticeUnaryExpression<Op, T1> &expr) {
GridFromExpression(grid,std::get<0>(expr.second));// recurse GridFromExpression(grid, std::get<0>(expr.second)); // recurse
} }
template <typename Op, typename T1, typename T2> template <typename Op, typename T1, typename T2>
inline void GridFromExpression(GridBase * &grid,const LatticeBinaryExpression<Op,T1,T2> &expr) inline void GridFromExpression(
{ GridBase *&grid, const LatticeBinaryExpression<Op, T1, T2> &expr) {
GridFromExpression(grid,std::get<0>(expr.second));// recurse GridFromExpression(grid, std::get<0>(expr.second)); // recurse
GridFromExpression(grid,std::get<1>(expr.second)); GridFromExpression(grid, std::get<1>(expr.second));
} }
template <typename Op, typename T1, typename T2, typename T3> template <typename Op, typename T1, typename T2, typename T3>
inline void GridFromExpression( GridBase * &grid,const LatticeTrinaryExpression<Op,T1,T2,T3 > &expr) inline void GridFromExpression(
{ GridBase *&grid, const LatticeTrinaryExpression<Op, T1, T2, T3> &expr) {
GridFromExpression(grid,std::get<0>(expr.second));// recurse GridFromExpression(grid, std::get<0>(expr.second)); // recurse
GridFromExpression(grid,std::get<1>(expr.second)); GridFromExpression(grid, std::get<1>(expr.second));
GridFromExpression(grid,std::get<2>(expr.second)); GridFromExpression(grid, std::get<2>(expr.second));
} }
////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////
// Obtain the CB from an expression, ensuring conformable. This must follow a tree recursion // Obtain the CB from an expression, ensuring conformable. This must follow a
// tree recursion
////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////
template<class T1, typename std::enable_if<is_lattice<T1>::value, T1>::type * =nullptr > template <class T1,
inline void CBFromExpression(int &cb,const T1& lat) // Lattice leaf typename std::enable_if<is_lattice<T1>::value, T1>::type * = nullptr>
inline void CBFromExpression(int &cb, const T1 &lat) // Lattice leaf
{ {
if ( (cb==Odd) || (cb==Even) ) { if ((cb == Odd) || (cb == Even)) {
assert(cb==lat.checkerboard); assert(cb == lat.checkerboard);
} }
cb=lat.checkerboard; cb = lat.checkerboard;
// std::cout<<GridLogMessage<<"Lattice leaf cb "<<cb<<std::endl; // std::cout<<GridLogMessage<<"Lattice leaf cb "<<cb<<std::endl;
} }
template<class T1,typename std::enable_if<!is_lattice<T1>::value, T1>::type * = nullptr > template <class T1,
inline void CBFromExpression(int &cb,const T1& notlat) // non-lattice leaf typename std::enable_if<!is_lattice<T1>::value, T1>::type * = nullptr>
inline void CBFromExpression(int &cb, const T1 &notlat) // non-lattice leaf
{ {
// std::cout<<GridLogMessage<<"Non lattice leaf cb"<<cb<<std::endl; // std::cout<<GridLogMessage<<"Non lattice leaf cb"<<cb<<std::endl;
} }
template <typename Op, typename T1> template <typename Op, typename T1>
inline void CBFromExpression(int &cb,const LatticeUnaryExpression<Op,T1 > &expr) inline void CBFromExpression(int &cb,
{ const LatticeUnaryExpression<Op, T1> &expr) {
CBFromExpression(cb,std::get<0>(expr.second));// recurse CBFromExpression(cb, std::get<0>(expr.second)); // recurse
// std::cout<<GridLogMessage<<"Unary node cb "<<cb<<std::endl; // std::cout<<GridLogMessage<<"Unary node cb "<<cb<<std::endl;
} }
template <typename Op, typename T1, typename T2> template <typename Op, typename T1, typename T2>
inline void CBFromExpression(int &cb,const LatticeBinaryExpression<Op,T1,T2> &expr) inline void CBFromExpression(int &cb,
{ const LatticeBinaryExpression<Op, T1, T2> &expr) {
CBFromExpression(cb,std::get<0>(expr.second));// recurse CBFromExpression(cb, std::get<0>(expr.second)); // recurse
CBFromExpression(cb,std::get<1>(expr.second)); CBFromExpression(cb, std::get<1>(expr.second));
// std::cout<<GridLogMessage<<"Binary node cb "<<cb<<std::endl; // std::cout<<GridLogMessage<<"Binary node cb "<<cb<<std::endl;
} }
template <typename Op, typename T1, typename T2, typename T3> template <typename Op, typename T1, typename T2, typename T3>
inline void CBFromExpression( int &cb,const LatticeTrinaryExpression<Op,T1,T2,T3 > &expr) inline void CBFromExpression(
{ int &cb, const LatticeTrinaryExpression<Op, T1, T2, T3> &expr) {
CBFromExpression(cb,std::get<0>(expr.second));// recurse CBFromExpression(cb, std::get<0>(expr.second)); // recurse
CBFromExpression(cb,std::get<1>(expr.second)); CBFromExpression(cb, std::get<1>(expr.second));
CBFromExpression(cb,std::get<2>(expr.second)); CBFromExpression(cb, std::get<2>(expr.second));
// std::cout<<GridLogMessage<<"Trinary node cb "<<cb<<std::endl; // std::cout<<GridLogMessage<<"Trinary node cb "<<cb<<std::endl;
} }
//////////////////////////////////////////// ////////////////////////////////////////////
// Unary operators and funcs // Unary operators and funcs
//////////////////////////////////////////// ////////////////////////////////////////////
#define GridUnopClass(name,ret)\ #define GridUnopClass(name, ret) \
template <class arg> struct name \ template <class arg> \
{ \ struct name { \
static auto inline func(const arg a)-> decltype(ret) { return ret; } \ static auto inline func(const arg a) -> decltype(ret) { return ret; } \
}; };
GridUnopClass(UnarySub,-a); GridUnopClass(UnarySub, -a);
GridUnopClass(UnaryNot,Not(a)); GridUnopClass(UnaryNot, Not(a));
GridUnopClass(UnaryAdj,adj(a)); GridUnopClass(UnaryAdj, adj(a));
GridUnopClass(UnaryConj,conjugate(a)); GridUnopClass(UnaryConj, conjugate(a));
GridUnopClass(UnaryTrace,trace(a)); GridUnopClass(UnaryTrace, trace(a));
GridUnopClass(UnaryTranspose,transpose(a)); GridUnopClass(UnaryTranspose, transpose(a));
GridUnopClass(UnaryTa,Ta(a)); GridUnopClass(UnaryTa, Ta(a));
GridUnopClass(UnaryProjectOnGroup,ProjectOnGroup(a)); GridUnopClass(UnaryProjectOnGroup, ProjectOnGroup(a));
GridUnopClass(UnaryReal,real(a)); GridUnopClass(UnaryReal, real(a));
GridUnopClass(UnaryImag,imag(a)); GridUnopClass(UnaryImag, imag(a));
GridUnopClass(UnaryToReal,toReal(a)); GridUnopClass(UnaryToReal, toReal(a));
GridUnopClass(UnaryToComplex,toComplex(a)); GridUnopClass(UnaryToComplex, toComplex(a));
GridUnopClass(UnaryTimesI,timesI(a)); GridUnopClass(UnaryTimesI, timesI(a));
GridUnopClass(UnaryTimesMinusI,timesMinusI(a)); GridUnopClass(UnaryTimesMinusI, timesMinusI(a));
GridUnopClass(UnaryAbs,abs(a)); GridUnopClass(UnaryAbs, abs(a));
GridUnopClass(UnarySqrt,sqrt(a)); GridUnopClass(UnarySqrt, sqrt(a));
GridUnopClass(UnaryRsqrt,rsqrt(a)); GridUnopClass(UnaryRsqrt, rsqrt(a));
GridUnopClass(UnarySin,sin(a)); GridUnopClass(UnarySin, sin(a));
GridUnopClass(UnaryCos,cos(a)); GridUnopClass(UnaryCos, cos(a));
GridUnopClass(UnaryAsin,asin(a)); GridUnopClass(UnaryAsin, asin(a));
GridUnopClass(UnaryAcos,acos(a)); GridUnopClass(UnaryAcos, acos(a));
GridUnopClass(UnaryLog,log(a)); GridUnopClass(UnaryLog, log(a));
GridUnopClass(UnaryExp,exp(a)); GridUnopClass(UnaryExp, exp(a));
//////////////////////////////////////////// ////////////////////////////////////////////
// Binary operators // Binary operators
//////////////////////////////////////////// ////////////////////////////////////////////
#define GridBinOpClass(name,combination)\ #define GridBinOpClass(name, combination) \
template <class left,class right>\ template <class left, class right> \
struct name\ struct name { \
{\ static auto inline func(const left &lhs, const right &rhs) \
static auto inline func(const left &lhs,const right &rhs)-> decltype(combination) const \ -> decltype(combination) const { \
{\ return combination; \
return combination;\ } \
}\ }
} GridBinOpClass(BinaryAdd, lhs + rhs);
GridBinOpClass(BinaryAdd,lhs+rhs); GridBinOpClass(BinarySub, lhs - rhs);
GridBinOpClass(BinarySub,lhs-rhs); GridBinOpClass(BinaryMul, lhs *rhs);
GridBinOpClass(BinaryMul,lhs*rhs);
GridBinOpClass(BinaryAnd ,lhs&rhs); GridBinOpClass(BinaryAnd, lhs &rhs);
GridBinOpClass(BinaryOr ,lhs|rhs); GridBinOpClass(BinaryOr, lhs | rhs);
GridBinOpClass(BinaryAndAnd,lhs&&rhs); GridBinOpClass(BinaryAndAnd, lhs &&rhs);
GridBinOpClass(BinaryOrOr ,lhs||rhs); GridBinOpClass(BinaryOrOr, lhs || rhs);
//////////////////////////////////////////////////// ////////////////////////////////////////////////////
// Trinary conditional op // Trinary conditional op
//////////////////////////////////////////////////// ////////////////////////////////////////////////////
#define GridTrinOpClass(name,combination)\ #define GridTrinOpClass(name, combination) \
template <class predicate,class left, class right> \ template <class predicate, class left, class right> \
struct name\ struct name { \
{\ static auto inline func(const predicate &pred, const left &lhs, \
static auto inline func(const predicate &pred,const left &lhs,const right &rhs)-> decltype(combination) const \ const right &rhs) -> decltype(combination) const { \
{\ return combination; \
return combination;\ } \
}\ }
}
GridTrinOpClass(TrinaryWhere,(predicatedWhere<predicate, \ GridTrinOpClass(
typename std::remove_reference<left>::type, \ TrinaryWhere,
typename std::remove_reference<right>::type> (pred,lhs,rhs))); (predicatedWhere<predicate, typename std::remove_reference<left>::type,
typename std::remove_reference<right>::type>(pred, lhs,
rhs)));
//////////////////////////////////////////// ////////////////////////////////////////////
// Operator syntactical glue // Operator syntactical glue
//////////////////////////////////////////// ////////////////////////////////////////////
#define GRID_UNOP(name) name<decltype(eval(0, arg))>
#define GRID_BINOP(name) name<decltype(eval(0, lhs)), decltype(eval(0, rhs))>
#define GRID_TRINOP(name) name<decltype(eval(0, pred)), decltype(eval(0, lhs)), decltype(eval(0, rhs))>
#define GRID_DEF_UNOP(op, name)\ #define GRID_UNOP(name) name<decltype(eval(0, arg))>
template <typename T1,\ #define GRID_BINOP(name) name<decltype(eval(0, lhs)), decltype(eval(0, rhs))>
typename std::enable_if<is_lattice<T1>::value||is_lattice_expr<T1>::value, T1>::type* = nullptr> inline auto op(const T1 &arg) \ #define GRID_TRINOP(name) \
-> decltype(LatticeUnaryExpression<GRID_UNOP(name),const T1&>(std::make_pair(GRID_UNOP(name)(),std::forward_as_tuple(arg)))) \ name<decltype(eval(0, pred)), decltype(eval(0, lhs)), decltype(eval(0, rhs))>
{ return LatticeUnaryExpression<GRID_UNOP(name), const T1 &>(std::make_pair(GRID_UNOP(name)(),std::forward_as_tuple(arg))); }
#define GRID_BINOP_LEFT(op, name)\ #define GRID_DEF_UNOP(op, name) \
template <typename T1,typename T2,\ template <typename T1, \
typename std::enable_if<is_lattice<T1>::value||is_lattice_expr<T1>::value, T1>::type* = nullptr>\ typename std::enable_if<is_lattice<T1>::value || \
inline auto op(const T1 &lhs,const T2&rhs) \ is_lattice_expr<T1>::value, \
-> decltype(LatticeBinaryExpression<GRID_BINOP(name),const T1&,const T2 &>(std::make_pair(GRID_BINOP(name)(),\ T1>::type * = nullptr> \
std::forward_as_tuple(lhs, rhs)))) \ inline auto op(const T1 &arg) \
{\ ->decltype(LatticeUnaryExpression<GRID_UNOP(name), const T1 &>( \
return LatticeBinaryExpression<GRID_BINOP(name), const T1 &, const T2 &>(std::make_pair(GRID_BINOP(name)(),\ std::make_pair(GRID_UNOP(name)(), std::forward_as_tuple(arg)))) { \
std::forward_as_tuple(lhs, rhs))); \ return LatticeUnaryExpression<GRID_UNOP(name), const T1 &>( \
} std::make_pair(GRID_UNOP(name)(), std::forward_as_tuple(arg))); \
}
#define GRID_BINOP_RIGHT(op, name)\ #define GRID_BINOP_LEFT(op, name) \
template <typename T1,typename T2,\ template <typename T1, typename T2, \
typename std::enable_if<!is_lattice<T1>::value && !is_lattice_expr<T1>::value, T1>::type* = nullptr,\ typename std::enable_if<is_lattice<T1>::value || \
typename std::enable_if< is_lattice<T2>::value || is_lattice_expr<T2>::value, T2>::type* = nullptr> \ is_lattice_expr<T1>::value, \
inline auto op(const T1 &lhs,const T2&rhs) \ T1>::type * = nullptr> \
-> decltype(LatticeBinaryExpression<GRID_BINOP(name),const T1&,const T2 &>(std::make_pair(GRID_BINOP(name)(),\ inline auto op(const T1 &lhs, const T2 &rhs) \
std::forward_as_tuple(lhs, rhs)))) \ ->decltype( \
{\ LatticeBinaryExpression<GRID_BINOP(name), const T1 &, const T2 &>( \
return LatticeBinaryExpression<GRID_BINOP(name), const T1 &, const T2 &>(std::make_pair(GRID_BINOP(name)(),\ std::make_pair(GRID_BINOP(name)(), \
std::forward_as_tuple(lhs, rhs))); \ std::forward_as_tuple(lhs, rhs)))) { \
} return LatticeBinaryExpression<GRID_BINOP(name), const T1 &, const T2 &>( \
std::make_pair(GRID_BINOP(name)(), std::forward_as_tuple(lhs, rhs))); \
}
#define GRID_DEF_BINOP(op, name)\ #define GRID_BINOP_RIGHT(op, name) \
GRID_BINOP_LEFT(op,name);\ template <typename T1, typename T2, \
GRID_BINOP_RIGHT(op,name); typename std::enable_if<!is_lattice<T1>::value && \
!is_lattice_expr<T1>::value, \
T1>::type * = nullptr, \
typename std::enable_if<is_lattice<T2>::value || \
is_lattice_expr<T2>::value, \
T2>::type * = nullptr> \
inline auto op(const T1 &lhs, const T2 &rhs) \
->decltype( \
LatticeBinaryExpression<GRID_BINOP(name), const T1 &, const T2 &>( \
std::make_pair(GRID_BINOP(name)(), \
std::forward_as_tuple(lhs, rhs)))) { \
return LatticeBinaryExpression<GRID_BINOP(name), const T1 &, const T2 &>( \
std::make_pair(GRID_BINOP(name)(), std::forward_as_tuple(lhs, rhs))); \
}
#define GRID_DEF_BINOP(op, name) \
GRID_BINOP_LEFT(op, name); \
GRID_BINOP_RIGHT(op, name);
#define GRID_DEF_TRINOP(op, name)\ #define GRID_DEF_TRINOP(op, name) \
template <typename T1,typename T2,typename T3> inline auto op(const T1 &pred,const T2&lhs,const T3 &rhs) \ template <typename T1, typename T2, typename T3> \
-> decltype(LatticeTrinaryExpression<GRID_TRINOP(name),const T1&,const T2 &,const T3&>(std::make_pair(GRID_TRINOP(name)(),\ inline auto op(const T1 &pred, const T2 &lhs, const T3 &rhs) \
std::forward_as_tuple(pred,lhs,rhs)))) \ ->decltype( \
{\ LatticeTrinaryExpression<GRID_TRINOP(name), const T1 &, const T2 &, \
return LatticeTrinaryExpression<GRID_TRINOP(name), const T1 &, const T2 &,const T3&>(std::make_pair(GRID_TRINOP(name)(), \ const T3 &>(std::make_pair( \
std::forward_as_tuple(pred,lhs, rhs))); \ GRID_TRINOP(name)(), std::forward_as_tuple(pred, lhs, rhs)))) { \
} return LatticeTrinaryExpression<GRID_TRINOP(name), const T1 &, const T2 &, \
const T3 &>(std::make_pair( \
GRID_TRINOP(name)(), std::forward_as_tuple(pred, lhs, rhs))); \
}
//////////////////////// ////////////////////////
//Operator definitions // Operator definitions
//////////////////////// ////////////////////////
GRID_DEF_UNOP(operator -,UnarySub); GRID_DEF_UNOP(operator-, UnarySub);
GRID_DEF_UNOP(Not,UnaryNot); GRID_DEF_UNOP(Not, UnaryNot);
GRID_DEF_UNOP(operator !,UnaryNot); GRID_DEF_UNOP(operator!, UnaryNot);
GRID_DEF_UNOP(adj,UnaryAdj); GRID_DEF_UNOP(adj, UnaryAdj);
GRID_DEF_UNOP(conjugate,UnaryConj); GRID_DEF_UNOP(conjugate, UnaryConj);
GRID_DEF_UNOP(trace,UnaryTrace); GRID_DEF_UNOP(trace, UnaryTrace);
GRID_DEF_UNOP(transpose,UnaryTranspose); GRID_DEF_UNOP(transpose, UnaryTranspose);
GRID_DEF_UNOP(Ta,UnaryTa); GRID_DEF_UNOP(Ta, UnaryTa);
GRID_DEF_UNOP(ProjectOnGroup,UnaryProjectOnGroup); GRID_DEF_UNOP(ProjectOnGroup, UnaryProjectOnGroup);
GRID_DEF_UNOP(real,UnaryReal); GRID_DEF_UNOP(real, UnaryReal);
GRID_DEF_UNOP(imag,UnaryImag); GRID_DEF_UNOP(imag, UnaryImag);
GRID_DEF_UNOP(toReal,UnaryToReal); GRID_DEF_UNOP(toReal, UnaryToReal);
GRID_DEF_UNOP(toComplex,UnaryToComplex); GRID_DEF_UNOP(toComplex, UnaryToComplex);
GRID_DEF_UNOP(timesI,UnaryTimesI); GRID_DEF_UNOP(timesI, UnaryTimesI);
GRID_DEF_UNOP(timesMinusI,UnaryTimesMinusI); GRID_DEF_UNOP(timesMinusI, UnaryTimesMinusI);
GRID_DEF_UNOP(abs ,UnaryAbs); //abs overloaded in cmath C++98; DON'T do the abs-fabs-dabs-labs thing GRID_DEF_UNOP(abs, UnaryAbs); // abs overloaded in cmath C++98; DON'T do the
GRID_DEF_UNOP(sqrt ,UnarySqrt); // abs-fabs-dabs-labs thing
GRID_DEF_UNOP(rsqrt,UnaryRsqrt); GRID_DEF_UNOP(sqrt, UnarySqrt);
GRID_DEF_UNOP(sin ,UnarySin); GRID_DEF_UNOP(rsqrt, UnaryRsqrt);
GRID_DEF_UNOP(cos ,UnaryCos); GRID_DEF_UNOP(sin, UnarySin);
GRID_DEF_UNOP(asin ,UnaryAsin); GRID_DEF_UNOP(cos, UnaryCos);
GRID_DEF_UNOP(acos ,UnaryAcos); GRID_DEF_UNOP(asin, UnaryAsin);
GRID_DEF_UNOP(log ,UnaryLog); GRID_DEF_UNOP(acos, UnaryAcos);
GRID_DEF_UNOP(exp ,UnaryExp); GRID_DEF_UNOP(log, UnaryLog);
GRID_DEF_UNOP(exp, UnaryExp);
GRID_DEF_BINOP(operator+,BinaryAdd); GRID_DEF_BINOP(operator+, BinaryAdd);
GRID_DEF_BINOP(operator-,BinarySub); GRID_DEF_BINOP(operator-, BinarySub);
GRID_DEF_BINOP(operator*,BinaryMul); GRID_DEF_BINOP(operator*, BinaryMul);
GRID_DEF_BINOP(operator&,BinaryAnd); GRID_DEF_BINOP(operator&, BinaryAnd);
GRID_DEF_BINOP(operator|,BinaryOr); GRID_DEF_BINOP(operator|, BinaryOr);
GRID_DEF_BINOP(operator&&,BinaryAndAnd); GRID_DEF_BINOP(operator&&, BinaryAndAnd);
GRID_DEF_BINOP(operator||,BinaryOrOr); GRID_DEF_BINOP(operator||, BinaryOrOr);
GRID_DEF_TRINOP(where,TrinaryWhere); GRID_DEF_TRINOP(where, TrinaryWhere);
///////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////
// Closure convenience to force expression to evaluate // Closure convenience to force expression to evaluate
///////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////
template<class Op,class T1> template <class Op, class T1>
auto closure(const LatticeUnaryExpression<Op,T1> & expr) auto closure(const LatticeUnaryExpression<Op, T1> &expr)
-> Lattice<decltype(expr.first.func(eval(0,std::get<0>(expr.second))))> -> Lattice<decltype(expr.first.func(eval(0, std::get<0>(expr.second))))> {
{ Lattice<decltype(expr.first.func(eval(0, std::get<0>(expr.second))))> ret(
Lattice<decltype(expr.first.func(eval(0,std::get<0>(expr.second))))> ret(expr); expr);
return ret; return ret;
} }
template<class Op,class T1, class T2> template <class Op, class T1, class T2>
auto closure(const LatticeBinaryExpression<Op,T1,T2> & expr) auto closure(const LatticeBinaryExpression<Op, T1, T2> &expr)
-> Lattice<decltype(expr.first.func(eval(0,std::get<0>(expr.second)), -> Lattice<decltype(expr.first.func(eval(0, std::get<0>(expr.second)),
eval(0,std::get<1>(expr.second))))> eval(0, std::get<1>(expr.second))))> {
{ Lattice<decltype(expr.first.func(eval(0, std::get<0>(expr.second)),
Lattice<decltype(expr.first.func(eval(0,std::get<0>(expr.second)), eval(0, std::get<1>(expr.second))))>
eval(0,std::get<1>(expr.second))))> ret(expr); ret(expr);
return ret; return ret;
} }
template<class Op,class T1, class T2, class T3> template <class Op, class T1, class T2, class T3>
auto closure(const LatticeTrinaryExpression<Op,T1,T2,T3> & expr) auto closure(const LatticeTrinaryExpression<Op, T1, T2, T3> &expr)
-> Lattice<decltype(expr.first.func(eval(0,std::get<0>(expr.second)), -> Lattice<decltype(expr.first.func(eval(0, std::get<0>(expr.second)),
eval(0,std::get<1>(expr.second)), eval(0, std::get<1>(expr.second)),
eval(0,std::get<2>(expr.second))))> eval(0, std::get<2>(expr.second))))> {
{ Lattice<decltype(expr.first.func(eval(0, std::get<0>(expr.second)),
Lattice<decltype(expr.first.func(eval(0,std::get<0>(expr.second)), eval(0, std::get<1>(expr.second)),
eval(0,std::get<1>(expr.second)), eval(0, std::get<2>(expr.second))))>
eval(0,std::get<2>(expr.second))))> ret(expr); ret(expr);
return ret; return ret;
} }
@ -390,12 +422,11 @@ template<class Op,class T1, class T2, class T3>
#undef GRID_DEF_UNOP #undef GRID_DEF_UNOP
#undef GRID_DEF_BINOP #undef GRID_DEF_BINOP
#undef GRID_DEF_TRINOP #undef GRID_DEF_TRINOP
} }
#if 0 #if 0
using namespace Grid; using namespace Grid;
int main(int argc,char **argv){ int main(int argc,char **argv){
Lattice<double> v1(16); Lattice<double> v1(16);
@ -405,7 +436,7 @@ using namespace Grid;
BinaryAdd<double,double> tmp; BinaryAdd<double,double> tmp;
LatticeBinaryExpression<BinaryAdd<double,double>,Lattice<double> &,Lattice<double> &> LatticeBinaryExpression<BinaryAdd<double,double>,Lattice<double> &,Lattice<double> &>
expr(std::make_pair(tmp, expr(std::make_pair(tmp,
std::forward_as_tuple(v1,v2))); std::forward_as_tuple(v1,v2)));
tmp.func(eval(0,v1),eval(0,v2)); tmp.func(eval(0,v1),eval(0,v2));
auto var = v1+v2; auto var = v1+v2;

334
lib/qcd/action/pseudofermion/OneFlavourEvenOddRational.h
View File

@ -1,212 +1,214 @@
/************************************************************************************* /*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/action/pseudofermion/OneFlavourEvenOddRational.h Source file: ./lib/qcd/action/pseudofermion/OneFlavourEvenOddRational.h
Copyright (C) 2015 Copyright (C) 2015
Author: Peter Boyle <paboyle@ph.ed.ac.uk> Author: Peter Boyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or the Free Software Foundation; either version 2 of the License, or
(at your option) any later version. (at your option) any later version.
This program is distributed in the hope that it will be useful, This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details. GNU General Public License for more details.
You should have received a copy of the GNU General Public License along You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc., with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory See the full license in the file "LICENSE" in the top level distribution
*************************************************************************************/ directory
/* END LEGAL */ *************************************************************************************/
/* END LEGAL */
#ifndef QCD_PSEUDOFERMION_ONE_FLAVOUR_EVEN_ODD_RATIONAL_H #ifndef QCD_PSEUDOFERMION_ONE_FLAVOUR_EVEN_ODD_RATIONAL_H
#define QCD_PSEUDOFERMION_ONE_FLAVOUR_EVEN_ODD_RATIONAL_H #define QCD_PSEUDOFERMION_ONE_FLAVOUR_EVEN_ODD_RATIONAL_H
namespace Grid{ namespace Grid {
namespace QCD{ namespace QCD {
/////////////////////////////////////// ///////////////////////////////////////
// One flavour rational // One flavour rational
/////////////////////////////////////// ///////////////////////////////////////
// S_f = chi^dag * N(Mpc^dag*Mpc)/D(Mpc^dag*Mpc) * chi // S_f = chi^dag * N(Mpc^dag*Mpc)/D(Mpc^dag*Mpc) * chi
//
// Here, M is some operator
// N and D makeup the rat. poly
//
template <class Impl>
class OneFlavourEvenOddRationalPseudoFermionAction
: public Action<typename Impl::GaugeField> {
public:
INHERIT_IMPL_TYPES(Impl);
typedef OneFlavourRationalParams Params;
Params param;
MultiShiftFunction PowerHalf;
MultiShiftFunction PowerNegHalf;
MultiShiftFunction PowerQuarter;
MultiShiftFunction PowerNegQuarter;
private:
FermionOperator<Impl> &FermOp; // the basic operator
// NOT using "Nroots"; IroIro is -- perhaps later, but this wasn't good for us
// historically
// and hasenbusch works better
FermionField PhiEven; // the pseudo fermion field for this trajectory
FermionField PhiOdd; // the pseudo fermion field for this trajectory
public:
OneFlavourEvenOddRationalPseudoFermionAction(FermionOperator<Impl> &Op,
Params &p)
: FermOp(Op),
PhiEven(Op.FermionRedBlackGrid()),
PhiOdd(Op.FermionRedBlackGrid()),
param(p) {
AlgRemez remez(param.lo, param.hi, param.precision);
// MdagM^(+- 1/2)
std::cout << GridLogMessage << "Generating degree " << param.degree
<< " for x^(1/2)" << std::endl;
remez.generateApprox(param.degree, 1, 2);
PowerHalf.Init(remez, param.tolerance, false);
PowerNegHalf.Init(remez, param.tolerance, true);
// MdagM^(+- 1/4)
std::cout << GridLogMessage << "Generating degree " << param.degree
<< " for x^(1/4)" << std::endl;
remez.generateApprox(param.degree, 1, 4);
PowerQuarter.Init(remez, param.tolerance, false);
PowerNegQuarter.Init(remez, param.tolerance, true);
};
virtual void refresh(const GaugeField &U, GridParallelRNG &pRNG) {
// P(phi) = e^{- phi^dag (MpcdagMpc)^-1/2 phi}
// = e^{- phi^dag (MpcdagMpc)^-1/4 (MpcdagMpc)^-1/4 phi}
// Phi = MpcdagMpc^{1/4} eta
// //
// Here, M is some operator // P(eta) = e^{- eta^dag eta}
// N and D makeup the rat. poly
// //
// e^{x^2/2 sig^2} => sig^2 = 0.5.
template<class Impl> //
class OneFlavourEvenOddRationalPseudoFermionAction : public Action<typename Impl::GaugeField> { // So eta should be of width sig = 1/sqrt(2).
public:
INHERIT_IMPL_TYPES(Impl);
typedef OneFlavourRationalParams Params; RealD scale = std::sqrt(0.5);
Params param;
MultiShiftFunction PowerHalf ; FermionField eta(FermOp.FermionGrid());
MultiShiftFunction PowerNegHalf; FermionField etaOdd(FermOp.FermionRedBlackGrid());
MultiShiftFunction PowerQuarter; FermionField etaEven(FermOp.FermionRedBlackGrid());
MultiShiftFunction PowerNegQuarter;
private: gaussian(pRNG, eta);
eta = eta * scale;
FermionOperator<Impl> & FermOp;// the basic operator
// NOT using "Nroots"; IroIro is -- perhaps later, but this wasn't good for us historically pickCheckerboard(Even, etaEven, eta);
// and hasenbusch works better pickCheckerboard(Odd, etaOdd, eta);
FermionField PhiEven; // the pseudo fermion field for this trajectory FermOp.ImportGauge(U);
FermionField PhiOdd; // the pseudo fermion field for this trajectory
public: // mutishift CG
SchurDifferentiableOperator<Impl> Mpc(FermOp);
ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter, PowerQuarter);
msCG(Mpc, etaOdd, PhiOdd);
OneFlavourEvenOddRationalPseudoFermionAction(FermionOperator<Impl> &Op, //////////////////////////////////////////////////////
Params & p ) : FermOp(Op), // FIXME : Clover term not yet..
PhiEven(Op.FermionRedBlackGrid()), //////////////////////////////////////////////////////
PhiOdd (Op.FermionRedBlackGrid()),
param(p)
{
AlgRemez remez(param.lo,param.hi,param.precision);
// MdagM^(+- 1/2) assert(FermOp.ConstEE() == 1);
std::cout<<GridLogMessage << "Generating degree "<<param.degree<<" for x^(1/2)"<<std::endl; PhiEven = zero;
remez.generateApprox(param.degree,1,2); };
PowerHalf.Init(remez,param.tolerance,false);
PowerNegHalf.Init(remez,param.tolerance,true);
// MdagM^(+- 1/4) //////////////////////////////////////////////////////
std::cout<<GridLogMessage << "Generating degree "<<param.degree<<" for x^(1/4)"<<std::endl; // S = phi^dag (Mdag M)^-1/2 phi
remez.generateApprox(param.degree,1,4); //////////////////////////////////////////////////////
PowerQuarter.Init(remez,param.tolerance,false); virtual RealD S(const GaugeField &U) {
PowerNegQuarter.Init(remez,param.tolerance,true); FermOp.ImportGauge(U);
};
virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) {
// P(phi) = e^{- phi^dag (MpcdagMpc)^-1/2 phi} FermionField Y(FermOp.FermionRedBlackGrid());
// = e^{- phi^dag (MpcdagMpc)^-1/4 (MpcdagMpc)^-1/4 phi}
// Phi = MpcdagMpc^{1/4} eta
//
// P(eta) = e^{- eta^dag eta}
//
// e^{x^2/2 sig^2} => sig^2 = 0.5.
//
// So eta should be of width sig = 1/sqrt(2).
RealD scale = std::sqrt(0.5); SchurDifferentiableOperator<Impl> Mpc(FermOp);
FermionField eta (FermOp.FermionGrid()); ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter,
FermionField etaOdd (FermOp.FermionRedBlackGrid()); PowerNegQuarter);
FermionField etaEven(FermOp.FermionRedBlackGrid());
gaussian(pRNG,eta); eta=eta*scale; msCG(Mpc, PhiOdd, Y);
pickCheckerboard(Even,etaEven,eta); RealD action = norm2(Y);
pickCheckerboard(Odd,etaOdd,eta); std::cout << GridLogMessage << "Pseudofermion action FIXME -- is -1/4 "
"solve or -1/2 solve faster??? "
<< action << std::endl;
FermOp.ImportGauge(U); return action;
};
// mutishift CG //////////////////////////////////////////////////////
SchurDifferentiableOperator<Impl> Mpc(FermOp); // Need
ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter,PowerQuarter); // dS_f/dU = chi^dag d[N/D] chi
msCG(Mpc,etaOdd,PhiOdd); //
// N/D is expressed as partial fraction expansion:
//
// a0 + \sum_k ak/(M^dagM + bk)
//
// d[N/D] is then
//
// \sum_k -ak [M^dagM+bk]^{-1} [ dM^dag M + M^dag dM ] [M^dag M +
// bk]^{-1}
//
// Need
// Mf Phi_k = [MdagM+bk]^{-1} Phi
// Mf Phi = \sum_k ak [MdagM+bk]^{-1} Phi
//
// With these building blocks
//
// dS/dU = \sum_k -ak Mf Phi_k^dag [ dM^dag M + M^dag dM ] Mf
// Phi_k
// S = innerprodReal(Phi,Mf Phi);
//////////////////////////////////////////////////////
virtual void deriv(const GaugeField &U, GaugeField &dSdU) {
const int Npole = PowerNegHalf.poles.size();
////////////////////////////////////////////////////// std::vector<FermionField> MPhi_k(Npole, FermOp.FermionRedBlackGrid());
// FIXME : Clover term not yet..
//////////////////////////////////////////////////////
assert(FermOp.ConstEE() == 1); FermionField X(FermOp.FermionRedBlackGrid());
PhiEven = zero; FermionField Y(FermOp.FermionRedBlackGrid());
};
////////////////////////////////////////////////////// GaugeField tmp(FermOp.GaugeGrid());
// S = phi^dag (Mdag M)^-1/2 phi
//////////////////////////////////////////////////////
virtual RealD S(const GaugeField &U) {
FermOp.ImportGauge(U); FermOp.ImportGauge(U);
FermionField Y(FermOp.FermionRedBlackGrid()); SchurDifferentiableOperator<Impl> Mpc(FermOp);
SchurDifferentiableOperator<Impl> Mpc(FermOp);
ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter,PowerNegQuarter); ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter, PowerNegHalf);
msCG(Mpc,PhiOdd,Y); msCG(Mpc, PhiOdd, MPhi_k);
RealD action = norm2(Y); dSdU = zero;
std::cout << GridLogMessage << "Pseudofermion action FIXME -- is -1/4 solve or -1/2 solve faster??? "<<action<<std::endl; for (int k = 0; k < Npole; k++) {
RealD ak = PowerNegHalf.residues[k];
return action; X = MPhi_k[k];
};
////////////////////////////////////////////////////// Mpc.Mpc(X, Y);
// Need Mpc.MpcDeriv(tmp, Y, X);
// dS_f/dU = chi^dag d[N/D] chi dSdU = dSdU + ak * tmp;
// Mpc.MpcDagDeriv(tmp, X, Y);
// N/D is expressed as partial fraction expansion: dSdU = dSdU + ak * tmp;
// }
// a0 + \sum_k ak/(M^dagM + bk)
//
// d[N/D] is then
//
// \sum_k -ak [M^dagM+bk]^{-1} [ dM^dag M + M^dag dM ] [M^dag M + bk]^{-1}
//
// Need
// Mf Phi_k = [MdagM+bk]^{-1} Phi
// Mf Phi = \sum_k ak [MdagM+bk]^{-1} Phi
//
// With these building blocks
//
// dS/dU = \sum_k -ak Mf Phi_k^dag [ dM^dag M + M^dag dM ] Mf Phi_k
// S = innerprodReal(Phi,Mf Phi);
//////////////////////////////////////////////////////
virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
const int Npole = PowerNegHalf.poles.size(); // dSdU = Ta(dSdU);
};
std::vector<FermionField> MPhi_k (Npole,FermOp.FermionRedBlackGrid()); };
}
FermionField X(FermOp.FermionRedBlackGrid());
FermionField Y(FermOp.FermionRedBlackGrid());
GaugeField tmp(FermOp.GaugeGrid());
FermOp.ImportGauge(U);
SchurDifferentiableOperator<Impl> Mpc(FermOp);
ConjugateGradientMultiShift<FermionField> msCG(param.MaxIter,PowerNegHalf);
msCG(Mpc,PhiOdd,MPhi_k);
dSdU = zero;
for(int k=0;k<Npole;k++){
RealD ak = PowerNegHalf.residues[k];
X = MPhi_k[k];
Mpc.Mpc(X,Y);
Mpc.MpcDeriv (tmp , Y, X ); dSdU=dSdU+ak*tmp;
Mpc.MpcDagDeriv(tmp , X, Y ); dSdU=dSdU+ak*tmp;
}
//dSdU = Ta(dSdU);
};
};
}
} }
#endif #endif

444
lib/tensors/Tensor_class.h
View File

@ -1,31 +1,32 @@
/************************************************************************************* /*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/tensors/Tensor_class.h Source file: ./lib/tensors/Tensor_class.h
Copyright (C) 2015 Copyright (C) 2015
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk> Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
Author: Peter Boyle <paboyle@ph.ed.ac.uk> Author: Peter Boyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or the Free Software Foundation; either version 2 of the License, or
(at your option) any later version. (at your option) any later version.
This program is distributed in the hope that it will be useful, This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details. GNU General Public License for more details.
You should have received a copy of the GNU General Public License along You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc., with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory See the full license in the file "LICENSE" in the top level distribution
*************************************************************************************/ directory
/* END LEGAL */ *************************************************************************************/
/* END LEGAL */
#ifndef GRID_MATH_TENSORS_H #ifndef GRID_MATH_TENSORS_H
#define GRID_MATH_TENSORS_H #define GRID_MATH_TENSORS_H
@ -38,17 +39,18 @@ namespace Grid {
// It is useful to NOT have any constructors // It is useful to NOT have any constructors
// so that these classes assert "is_pod<class> == true" // so that these classes assert "is_pod<class> == true"
// because then the standard C++ valarray container eliminates fill overhead on new allocation and // because then the standard C++ valarray container eliminates fill overhead on
// new allocation and
// non-move copying. // non-move copying.
// //
// However note that doing this eliminates some syntactical sugar such as // However note that doing this eliminates some syntactical sugar such as
// calling the constructor explicitly or implicitly // calling the constructor explicitly or implicitly
// //
class GridTensorBase {}; class GridTensorBase {};
template<class vtype> class iScalar template <class vtype>
{ class iScalar {
public: public:
vtype _internal; vtype _internal;
typedef vtype element; typedef vtype element;
@ -60,13 +62,14 @@ public:
typedef iScalar<recurse_scalar_object> scalar_object; typedef iScalar<recurse_scalar_object> scalar_object;
// substitutes a real or complex version with same tensor structure // substitutes a real or complex version with same tensor structure
typedef iScalar<typename GridTypeMapper<vtype>::Complexified > Complexified; typedef iScalar<typename GridTypeMapper<vtype>::Complexified> Complexified;
typedef iScalar<typename GridTypeMapper<vtype>::Realified > Realified; typedef iScalar<typename GridTypeMapper<vtype>::Realified> Realified;
enum { TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1}; enum { TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1 };
// Scalar no action // Scalar no action
// template<int Level> using tensor_reduce_level = typename iScalar<GridTypeMapper<vtype>::tensor_reduce_level<Level> >; // template<int Level> using tensor_reduce_level = typename
// iScalar<GridTypeMapper<vtype>::tensor_reduce_level<Level> >;
iScalar() = default; iScalar() = default;
/* /*
iScalar(const iScalar<vtype> &copyme)=default; iScalar(const iScalar<vtype> &copyme)=default;
@ -74,83 +77,106 @@ public:
iScalar<vtype> & operator= (const iScalar<vtype> &copyme) = default; iScalar<vtype> & operator= (const iScalar<vtype> &copyme) = default;
iScalar<vtype> & operator= (iScalar<vtype> &&copyme) = default; iScalar<vtype> & operator= (iScalar<vtype> &&copyme) = default;
*/ */
iScalar(scalar_type s) : _internal(s) {};// recurse down and hit the constructor for vector_type iScalar(scalar_type s)
iScalar(const Zero &z){ *this = zero; }; : _internal(s){}; // recurse down and hit the constructor for vector_type
iScalar(const Zero &z) { *this = zero; };
iScalar<vtype> & operator= (const Zero &hero){ iScalar<vtype> &operator=(const Zero &hero) {
zeroit(*this); zeroit(*this);
return *this; return *this;
} }
friend strong_inline void vstream(iScalar<vtype> &out,const iScalar<vtype> &in){ friend strong_inline void vstream(iScalar<vtype> &out,
vstream(out._internal,in._internal); const iScalar<vtype> &in) {
vstream(out._internal, in._internal);
} }
friend strong_inline void zeroit(iScalar<vtype> &that){ friend strong_inline void zeroit(iScalar<vtype> &that) {
zeroit(that._internal); zeroit(that._internal);
} }
friend strong_inline void prefetch(iScalar<vtype> &that){ friend strong_inline void prefetch(iScalar<vtype> &that) {
prefetch(that._internal); prefetch(that._internal);
} }
friend strong_inline void permute(iScalar<vtype> &out,const iScalar<vtype> &in,int permutetype){ friend strong_inline void permute(iScalar<vtype> &out,
permute(out._internal,in._internal,permutetype); const iScalar<vtype> &in, int permutetype) {
permute(out._internal, in._internal, permutetype);
} }
// Unary negation // Unary negation
friend strong_inline iScalar<vtype> operator -(const iScalar<vtype> &r) { friend strong_inline iScalar<vtype> operator-(const iScalar<vtype> &r) {
iScalar<vtype> ret; iScalar<vtype> ret;
ret._internal= -r._internal; ret._internal = -r._internal;
return ret; return ret;
} }
// *=,+=,-= operators inherit from corresponding "*,-,+" behaviour // *=,+=,-= operators inherit from corresponding "*,-,+" behaviour
strong_inline iScalar<vtype> &operator *=(const iScalar<vtype> &r) { strong_inline iScalar<vtype> &operator*=(const iScalar<vtype> &r) {
*this = (*this)*r; *this = (*this) * r;
return *this; return *this;
} }
strong_inline iScalar<vtype> &operator -=(const iScalar<vtype> &r) { strong_inline iScalar<vtype> &operator-=(const iScalar<vtype> &r) {
*this = (*this)-r; *this = (*this) - r;
return *this; return *this;
} }
strong_inline iScalar<vtype> &operator +=(const iScalar<vtype> &r) { strong_inline iScalar<vtype> &operator+=(const iScalar<vtype> &r) {
*this = (*this)+r; *this = (*this) + r;
return *this; return *this;
} }
strong_inline vtype & operator ()(void) { strong_inline vtype &operator()(void) { return _internal; }
return _internal; strong_inline const vtype &operator()(void) const { return _internal; }
}
strong_inline const vtype & operator ()(void) const {
return _internal;
}
// Type casts meta programmed, must be pure scalar to match TensorRemove // Type casts meta programmed, must be pure scalar to match TensorRemove
template<class U=vtype,class V=scalar_type,IfComplex<V> = 0,IfNotSimd<U> = 0> operator ComplexF () const { return(TensorRemove(_internal)); }; template <class U = vtype, class V = scalar_type, IfComplex<V> = 0,
template<class U=vtype,class V=scalar_type,IfComplex<V> = 0,IfNotSimd<U> = 0> operator ComplexD () const { return(TensorRemove(_internal)); }; IfNotSimd<U> = 0>
// template<class U=vtype,class V=scalar_type,IfComplex<V> = 0,IfNotSimd<U> = 0> operator RealD () const { return(real(TensorRemove(_internal))); } operator ComplexF() const {
template<class U=vtype,class V=scalar_type,IfReal<V> = 0,IfNotSimd<U> = 0> operator RealD () const { return TensorRemove(_internal); } return (TensorRemove(_internal));
template<class U=vtype,class V=scalar_type,IfInteger<V> = 0,IfNotSimd<U> = 0> operator Integer () const { return Integer(TensorRemove(_internal)); } };
template <class U = vtype, class V = scalar_type, IfComplex<V> = 0,
// convert from a something to a scalar via constructor of something arg IfNotSimd<U> = 0>
template<class T,typename std::enable_if<!isGridTensor<T>::value, T>::type* = nullptr > strong_inline iScalar<vtype> operator = (T arg) operator ComplexD() const {
{ return (TensorRemove(_internal));
_internal = arg; };
return *this; // template<class U=vtype,class V=scalar_type,IfComplex<V> = 0,IfNotSimd<U> =
} // 0> operator RealD () const { return(real(TensorRemove(_internal))); }
template <class U = vtype, class V = scalar_type, IfReal<V> = 0,
IfNotSimd<U> = 0>
operator RealD() const {
return TensorRemove(_internal);
}
template <class U = vtype, class V = scalar_type, IfInteger<V> = 0,
IfNotSimd<U> = 0>
operator Integer() const {
return Integer(TensorRemove(_internal));
}
friend std::ostream& operator<< (std::ostream& stream, const iScalar<vtype> &o){ // convert from a something to a scalar via constructor of something arg
stream<< "S {"<<o._internal<<"}"; template <class T, typename std::enable_if<!isGridTensor<T>::value, T>::type
return stream; * = nullptr>
}; strong_inline iScalar<vtype> operator=(T arg) {
_internal = arg;
return *this;
}
friend std::ostream &operator<<(std::ostream &stream,
const iScalar<vtype> &o) {
stream << "S {" << o._internal << "}";
return stream;
};
}; };
/////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////
// Allows to turn scalar<scalar<scalar<double>>>> back to double. // Allows to turn scalar<scalar<scalar<double>>>> back to double.
/////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////
template<class T> strong_inline typename std::enable_if<!isGridTensor<T>::value, T>::type TensorRemove(T arg) { return arg;} template <class T>
template<class vtype> strong_inline auto TensorRemove(iScalar<vtype> arg) -> decltype(TensorRemove(arg._internal)) strong_inline typename std::enable_if<!isGridTensor<T>::value, T>::type
{ TensorRemove(T arg) {
return arg;
}
template <class vtype>
strong_inline auto TensorRemove(iScalar<vtype> arg)
-> decltype(TensorRemove(arg._internal)) {
return TensorRemove(arg._internal); return TensorRemove(arg._internal);
} }
template<class vtype,int N> class iVector template <class vtype, int N>
{ class iVector {
public: public:
vtype _internal[N]; vtype _internal[N];
typedef vtype element; typedef vtype element;
@ -159,23 +185,23 @@ public:
typedef typename GridTypeMapper<vtype>::tensor_reduced tensor_reduced_v; typedef typename GridTypeMapper<vtype>::tensor_reduced tensor_reduced_v;
typedef typename GridTypeMapper<vtype>::scalar_object recurse_scalar_object; typedef typename GridTypeMapper<vtype>::scalar_object recurse_scalar_object;
typedef iScalar<tensor_reduced_v> tensor_reduced; typedef iScalar<tensor_reduced_v> tensor_reduced;
typedef iVector<recurse_scalar_object,N> scalar_object; typedef iVector<recurse_scalar_object, N> scalar_object;
// substitutes a real or complex version with same tensor structure // substitutes a real or complex version with same tensor structure
typedef iVector<typename GridTypeMapper<vtype>::Complexified,N > Complexified; typedef iVector<typename GridTypeMapper<vtype>::Complexified, N> Complexified;
typedef iVector<typename GridTypeMapper<vtype>::Realified,N > Realified; typedef iVector<typename GridTypeMapper<vtype>::Realified, N> Realified;
template<class T,typename std::enable_if<!isGridTensor<T>::value, T>::type* = nullptr > strong_inline auto operator = (T arg) -> iVector<vtype,N> template <class T, typename std::enable_if<!isGridTensor<T>::value, T>::type
{ * = nullptr>
zeroit(*this); strong_inline auto operator=(T arg) -> iVector<vtype, N> {
for(int i=0;i<N;i++) zeroit(*this);
_internal[i] = arg; for (int i = 0; i < N; i++) _internal[i] = arg;
return *this; return *this;
} }
enum { TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1}; enum { TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1 };
iVector(const Zero &z){ *this = zero; }; iVector(const Zero &z) { *this = zero; };
iVector() =default; iVector() = default;
/* /*
iVector(const iVector<vtype,N> &copyme)=default; iVector(const iVector<vtype,N> &copyme)=default;
iVector(iVector<vtype,N> &&copyme)=default; iVector(iVector<vtype,N> &&copyme)=default;
@ -183,71 +209,71 @@ public:
iVector<vtype,N> & operator= (iVector<vtype,N> &&copyme) = default; iVector<vtype,N> & operator= (iVector<vtype,N> &&copyme) = default;
*/ */
iVector<vtype,N> & operator= (const Zero &hero){ iVector<vtype, N> &operator=(const Zero &hero) {
zeroit(*this); zeroit(*this);
return *this; return *this;
} }
friend strong_inline void zeroit(iVector<vtype,N> &that){ friend strong_inline void zeroit(iVector<vtype, N> &that) {
for(int i=0;i<N;i++){ for (int i = 0; i < N; i++) {
zeroit(that._internal[i]); zeroit(that._internal[i]);
} }
} }
friend strong_inline void prefetch(iVector<vtype,N> &that){ friend strong_inline void prefetch(iVector<vtype, N> &that) {
for(int i=0;i<N;i++) prefetch(that._internal[i]); for (int i = 0; i < N; i++) prefetch(that._internal[i]);
} }
friend strong_inline void vstream(iVector<vtype,N> &out,const iVector<vtype,N> &in){ friend strong_inline void vstream(iVector<vtype, N> &out,
for(int i=0;i<N;i++){ const iVector<vtype, N> &in) {
vstream(out._internal[i],in._internal[i]); for (int i = 0; i < N; i++) {
vstream(out._internal[i], in._internal[i]);
} }
} }
friend strong_inline void permute(iVector<vtype,N> &out,const iVector<vtype,N> &in,int permutetype){ friend strong_inline void permute(iVector<vtype, N> &out,
for(int i=0;i<N;i++){ const iVector<vtype, N> &in,
permute(out._internal[i],in._internal[i],permutetype); int permutetype) {
for (int i = 0; i < N; i++) {
permute(out._internal[i], in._internal[i], permutetype);
} }
} }
// Unary negation // Unary negation
friend strong_inline iVector<vtype,N> operator -(const iVector<vtype,N> &r) { friend strong_inline iVector<vtype, N> operator-(const iVector<vtype, N> &r) {
iVector<vtype,N> ret; iVector<vtype, N> ret;
for(int i=0;i<N;i++) ret._internal[i]= -r._internal[i]; for (int i = 0; i < N; i++) ret._internal[i] = -r._internal[i];
return ret; return ret;
} }
// *=,+=,-= operators inherit from corresponding "*,-,+" behaviour // *=,+=,-= operators inherit from corresponding "*,-,+" behaviour
strong_inline iVector<vtype,N> &operator *=(const iScalar<vtype> &r) { strong_inline iVector<vtype, N> &operator*=(const iScalar<vtype> &r) {
*this = (*this)*r; *this = (*this) * r;
return *this; return *this;
} }
strong_inline iVector<vtype,N> &operator -=(const iVector<vtype,N> &r) { strong_inline iVector<vtype, N> &operator-=(const iVector<vtype, N> &r) {
*this = (*this)-r; *this = (*this) - r;
return *this; return *this;
} }
strong_inline iVector<vtype,N> &operator +=(const iVector<vtype,N> &r) { strong_inline iVector<vtype, N> &operator+=(const iVector<vtype, N> &r) {
*this = (*this)+r; *this = (*this) + r;
return *this; return *this;
} }
strong_inline vtype & operator ()(int i) { strong_inline vtype &operator()(int i) { return _internal[i]; }
return _internal[i]; strong_inline const vtype &operator()(int i) const { return _internal[i]; }
} friend std::ostream &operator<<(std::ostream &stream,
strong_inline const vtype & operator ()(int i) const { const iVector<vtype, N> &o) {
return _internal[i]; stream << "V<" << N << ">{";
} for (int i = 0; i < N; i++) {
friend std::ostream& operator<< (std::ostream& stream, const iVector<vtype,N> &o){ stream << o._internal[i];
stream<< "V<"<<N<<">{"; if (i < N - 1) stream << ",";
for(int i=0;i<N;i++) {
stream<<o._internal[i];
if (i<N-1) stream<<",";
} }
stream<<"}"; stream << "}";
return stream; return stream;
}; };
// strong_inline vtype && operator ()(int i) { // strong_inline vtype && operator ()(int i) {
// return _internal[i]; // return _internal[i];
// } // }
}; };
template<class vtype,int N> class iMatrix template <class vtype, int N>
{ class iMatrix {
public: public:
vtype _internal[N][N]; vtype _internal[N][N];
typedef vtype element; typedef vtype element;
@ -257,29 +283,27 @@ public:
typedef typename GridTypeMapper<vtype>::scalar_object recurse_scalar_object; typedef typename GridTypeMapper<vtype>::scalar_object recurse_scalar_object;
// substitutes a real or complex version with same tensor structure // substitutes a real or complex version with same tensor structure
typedef iMatrix<typename GridTypeMapper<vtype>::Complexified,N > Complexified; typedef iMatrix<typename GridTypeMapper<vtype>::Complexified, N> Complexified;
typedef iMatrix<typename GridTypeMapper<vtype>::Realified,N > Realified; typedef iMatrix<typename GridTypeMapper<vtype>::Realified, N> Realified;
// Tensure removal // Tensure removal
typedef iScalar<tensor_reduced_v> tensor_reduced; typedef iScalar<tensor_reduced_v> tensor_reduced;
typedef iMatrix<recurse_scalar_object,N> scalar_object; typedef iMatrix<recurse_scalar_object, N> scalar_object;
enum { TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1}; enum { TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1 };
iMatrix(const Zero &z) { *this = zero; };
iMatrix() = default;
iMatrix(const Zero &z){ *this = zero; }; iMatrix &operator=(const iMatrix &rhs) {
iMatrix() =default; for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++) vstream(_internal[i][j], rhs._internal[i][j]);
iMatrix& operator=(const iMatrix& rhs){
for(int i=0;i<N;i++)
for(int j=0;j<N;j++)
vstream(_internal[i][j],rhs._internal[i][j]);
return *this; return *this;
}; };
iMatrix(scalar_type s) { (*this) = s ;};// recurse down and hit the constructor for vector_type iMatrix(scalar_type s) {
(*this) = s;
}; // recurse down and hit the constructor for vector_type
/* /*
iMatrix(const iMatrix<vtype,N> &copyme)=default; iMatrix(const iMatrix<vtype,N> &copyme)=default;
@ -288,118 +312,118 @@ public:
iMatrix<vtype,N> & operator= (iMatrix<vtype,N> &&copyme) = default; iMatrix<vtype,N> & operator= (iMatrix<vtype,N> &&copyme) = default;
*/ */
iMatrix<vtype, N> &operator=(const Zero &hero) {
iMatrix<vtype,N> & operator= (const Zero &hero){
zeroit(*this); zeroit(*this);
return *this; return *this;
} }
template<class T,typename std::enable_if<!isGridTensor<T>::value, T>::type* = nullptr > strong_inline auto operator = (T arg) -> iMatrix<vtype,N> template <class T, typename std::enable_if<!isGridTensor<T>::value, T>::type
{ * = nullptr>
zeroit(*this); strong_inline auto operator=(T arg) -> iMatrix<vtype, N> {
for(int i=0;i<N;i++) zeroit(*this);
_internal[i][i] = arg; for (int i = 0; i < N; i++) _internal[i][i] = arg;
return *this; return *this;
}
friend strong_inline void zeroit(iMatrix<vtype, N> &that) {
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
zeroit(that._internal[i][j]);
}
} }
friend strong_inline void zeroit(iMatrix<vtype,N> &that){
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
zeroit(that._internal[i][j]);
}}
} }
friend strong_inline void prefetch(iMatrix<vtype,N> &that){ friend strong_inline void prefetch(iMatrix<vtype, N> &that) {
for(int i=0;i<N;i++) for (int i = 0; i < N; i++)
for(int j=0;j<N;j++) for (int j = 0; j < N; j++) prefetch(that._internal[i][j]);
prefetch(that._internal[i][j]);
} }
friend strong_inline void vstream(iMatrix<vtype,N> &out,const iMatrix<vtype,N> &in){ friend strong_inline void vstream(iMatrix<vtype, N> &out,
for(int i=0;i<N;i++){ const iMatrix<vtype, N> &in) {
for(int j=0;j<N;j++){ for (int i = 0; i < N; i++) {
vstream(out._internal[i][j],in._internal[i][j]); for (int j = 0; j < N; j++) {
}} vstream(out._internal[i][j], in._internal[i][j]);
}
} }
friend strong_inline void permute(iMatrix<vtype,N> &out,const iMatrix<vtype,N> &in,int permutetype){
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
permute(out._internal[i][j],in._internal[i][j],permutetype);
}}
} }
friend strong_inline void permute(iMatrix<vtype, N> &out,
const iMatrix<vtype, N> &in,
int permutetype) {
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
permute(out._internal[i][j], in._internal[i][j], permutetype);
}
}
}
// Unary negation // Unary negation
friend strong_inline iMatrix<vtype,N> operator -(const iMatrix<vtype,N> &r) { friend strong_inline iMatrix<vtype, N> operator-(const iMatrix<vtype, N> &r) {
iMatrix<vtype,N> ret; iMatrix<vtype, N> ret;
for(int i=0;i<N;i++){ for (int i = 0; i < N; i++) {
for(int j=0;j<N;j++){ for (int j = 0; j < N; j++) {
ret._internal[i][j]= -r._internal[i][j]; ret._internal[i][j] = -r._internal[i][j];
}} }
}
return ret; return ret;
} }
// *=,+=,-= operators inherit from corresponding "*,-,+" behaviour // *=,+=,-= operators inherit from corresponding "*,-,+" behaviour
template<class T> template <class T>
strong_inline iMatrix<vtype,N> &operator *=(const T &r) { strong_inline iMatrix<vtype, N> &operator*=(const T &r) {
*this = (*this)*r; *this = (*this) * r;
return *this; return *this;
} }
template<class T> template <class T>
strong_inline iMatrix<vtype,N> &operator -=(const T &r) { strong_inline iMatrix<vtype, N> &operator-=(const T &r) {
*this = (*this)-r; *this = (*this) - r;
return *this; return *this;
} }
template<class T> template <class T>
strong_inline iMatrix<vtype,N> &operator +=(const T &r) { strong_inline iMatrix<vtype, N> &operator+=(const T &r) {
*this = (*this)+r; *this = (*this) + r;
return *this; return *this;
} }
// returns an lvalue reference // returns an lvalue reference
strong_inline vtype & operator ()(int i,int j) { strong_inline vtype &operator()(int i, int j) { return _internal[i][j]; }
strong_inline const vtype &operator()(int i, int j) const {
return _internal[i][j]; return _internal[i][j];
} }
strong_inline const vtype & operator ()(int i,int j) const { friend std::ostream &operator<<(std::ostream &stream,
return _internal[i][j]; const iMatrix<vtype, N> &o) {
} stream << "M<" << N << ">{";
friend std::ostream& operator<< (std::ostream& stream, const iMatrix<vtype,N> &o){ for (int i = 0; i < N; i++) {
stream<< "M<"<<N<<">{"; stream << "{";
for(int i=0;i<N;i++) { for (int j = 0; j < N; j++) {
stream<< "{"; stream << o._internal[i][j];
for(int j=0;j<N;j++) { if (i < N - 1) stream << ",";
stream<<o._internal[i][j];
if (i<N-1) stream<<",";
} }
stream<<"}"; stream << "}";
if(i!=N-1) stream<<"\n\t\t"; if (i != N - 1) stream << "\n\t\t";
} }
stream<<"}"; stream << "}";
return stream; return stream;
}; };
// strong_inline vtype && operator ()(int i,int j) { // strong_inline vtype && operator ()(int i,int j) {
// return _internal[i][j]; // return _internal[i][j];
// } // }
}; };
template<class v> void vprefetch(const iScalar<v> &vv) template <class v>
{ void vprefetch(const iScalar<v> &vv) {
vprefetch(vv._internal); vprefetch(vv._internal);
} }
template<class v,int N> void vprefetch(const iVector<v,N> &vv) template <class v, int N>
{ void vprefetch(const iVector<v, N> &vv) {
for(int i=0;i<N;i++){ for (int i = 0; i < N; i++) {
vprefetch(vv._internal[i]); vprefetch(vv._internal[i]);
} }
} }
template<class v,int N> void vprefetch(const iMatrix<v,N> &vv) template <class v, int N>
{ void vprefetch(const iMatrix<v, N> &vv) {
for(int i=0;i<N;i++){ for (int i = 0; i < N; i++) {
for(int j=0;j<N;j++){ for (int j = 0; j < N; j++) {
vprefetch(vv._internal[i][j]); vprefetch(vv._internal[i][j]);
}} }
}
} }
} }
#endif #endif

105
tests/Test_rhmc_EOWilson1p1.cc
View File

@ -1,97 +1,100 @@
/************************************************************************************* /*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_rhmc_EOWilson1p1.cc Source file: ./tests/Test_rhmc_EOWilson1p1.cc
Copyright (C) 2015 Copyright (C) 2015
Author: Peter Boyle <paboyle@ph.ed.ac.uk> Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: paboyle <paboyle@ph.ed.ac.uk> Author: paboyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or the Free Software Foundation; either version 2 of the License, or
(at your option) any later version. (at your option) any later version.
This program is distributed in the hope that it will be useful, This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details. GNU General Public License for more details.
You should have received a copy of the GNU General Public License along You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc., with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory See the full license in the file "LICENSE" in the top level distribution
*************************************************************************************/ directory
/* END LEGAL */ *************************************************************************************/
/* END LEGAL */
#include "Grid.h" #include "Grid.h"
using namespace std; using namespace std;
using namespace Grid; using namespace Grid;
using namespace Grid::QCD; using namespace Grid::QCD;
namespace Grid { namespace Grid {
namespace QCD { namespace QCD {
class HmcRunner : public NerscHmcRunner { class HmcRunner : public NerscHmcRunner {
public: public:
void BuildTheAction(int argc, char **argv)
void BuildTheAction (int argc, char **argv)
{ {
typedef WilsonImplR ImplPolicy; typedef WilsonImplR ImplPolicy;
typedef WilsonFermionR FermionAction; typedef WilsonFermionR FermionAction;
typedef typename FermionAction::FermionField FermionField; typedef typename FermionAction::FermionField FermionField;
UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); UGrid = SpaceTimeGrid::makeFourDimGrid(
GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()),
GridDefaultMpi());
UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
FGrid = UGrid; FGrid = UGrid;
FrbGrid = UrbGrid; FrbGrid = UrbGrid;
// temporarily need a gauge field // temporarily need a gauge field
LatticeGaugeField U(UGrid); LatticeGaugeField U(UGrid);
// Gauge action // Gauge action
WilsonGaugeActionR Waction(5.6); WilsonGaugeActionR Waction(5.6);
Real mass=-0.77; Real mass = -0.77;
FermionAction FermOp(U,*FGrid,*FrbGrid,mass); FermionAction FermOp(U, *FGrid, *FrbGrid, mass);
// 1+1 flavour // 1+1 flavour
OneFlavourRationalParams Params(1.0e-4,64.0,1000,1.0e-6); OneFlavourRationalParams Params(1.0e-4, 64.0, 2000, 1.0e-6);
OneFlavourEvenOddRationalPseudoFermionAction<WilsonImplR> WilsonNf1a(FermOp,Params); OneFlavourEvenOddRationalPseudoFermionAction<WilsonImplR> WilsonNf1a(
OneFlavourEvenOddRationalPseudoFermionAction<WilsonImplR> WilsonNf1b(FermOp,Params); FermOp, Params);
OneFlavourEvenOddRationalPseudoFermionAction<WilsonImplR> WilsonNf1b(
FermOp, Params);
//Collect actions //Smearing on/off
WilsonNf1a.is_smeared = true;
WilsonNf1b.is_smeared = true;
// Collect actions
ActionLevel<LatticeGaugeField> Level1; ActionLevel<LatticeGaugeField> Level1;
Level1.push_back(&WilsonNf1a); Level1.push_back(&WilsonNf1a);
Level1.push_back(&WilsonNf1b); Level1.push_back(&WilsonNf1b);
Level1.push_back(&Waction); Level1.push_back(&Waction);
TheAction.push_back(Level1); TheAction.push_back(Level1);
Run(argc,argv); Run(argc, argv);
}; };
}; };
}
}}
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
HmcRunner TheHMC;
TheHMC.BuildTheAction(argc,argv);
} }
int main(int argc, char **argv) {
Grid_init(&argc, &argv);
int threads = GridThread::GetThreads();
std::cout << GridLogMessage << "Grid is setup to use " << threads
<< " threads" << std::endl;
HmcRunner TheHMC;
TheHMC.BuildTheAction(argc, argv);
}