mirror of
https://github.com/paboyle/Grid.git
synced 2025-06-15 14:27:06 +01:00
Compare commits
1 Commits
feature/se
...
feature/ad
Author | SHA1 | Date | |
---|---|---|---|
b284d50863 |
1
.gitignore
vendored
1
.gitignore
vendored
@ -88,7 +88,6 @@ Thumbs.db
|
|||||||
# build directory #
|
# build directory #
|
||||||
###################
|
###################
|
||||||
build*/*
|
build*/*
|
||||||
Documentation/_build
|
|
||||||
|
|
||||||
# IDE related files #
|
# IDE related files #
|
||||||
#####################
|
#####################
|
||||||
|
@ -7,6 +7,7 @@ Source file: ./lib/qcd/modules/plaquette.h
|
|||||||
Copyright (C) 2017
|
Copyright (C) 2017
|
||||||
|
|
||||||
Author: Guido Cossu <guido.cossu@ed.ac.uk>
|
Author: Guido Cossu <guido.cossu@ed.ac.uk>
|
||||||
|
Author: Chulwoo Jung <chulwoo@bnl.gov>
|
||||||
|
|
||||||
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
|
||||||
@ -35,7 +36,7 @@ template <class Gimpl>
|
|||||||
class WilsonFlow: public Smear<Gimpl>{
|
class WilsonFlow: public Smear<Gimpl>{
|
||||||
unsigned int Nstep;
|
unsigned int Nstep;
|
||||||
unsigned int measure_interval;
|
unsigned int measure_interval;
|
||||||
mutable RealD epsilon, taus;
|
mutable RealD epsilon, taus,tolerance;
|
||||||
|
|
||||||
|
|
||||||
mutable WilsonGaugeAction<Gimpl> SG;
|
mutable WilsonGaugeAction<Gimpl> SG;
|
||||||
@ -47,13 +48,15 @@ class WilsonFlow: public Smear<Gimpl>{
|
|||||||
public:
|
public:
|
||||||
INHERIT_GIMPL_TYPES(Gimpl)
|
INHERIT_GIMPL_TYPES(Gimpl)
|
||||||
|
|
||||||
explicit WilsonFlow(unsigned int Nstep, RealD epsilon, unsigned int interval = 1):
|
explicit WilsonFlow(unsigned int Nstep, RealD epsilon, unsigned int interval = 1, RealD tol = 1e-3):
|
||||||
Nstep(Nstep),
|
Nstep(Nstep),
|
||||||
epsilon(epsilon),
|
epsilon(epsilon),
|
||||||
|
tolerance(tol),
|
||||||
measure_interval(interval),
|
measure_interval(interval),
|
||||||
SG(WilsonGaugeAction<Gimpl>(3.0)) {
|
SG(WilsonGaugeAction<Gimpl>(3.0)) {
|
||||||
// WilsonGaugeAction with beta 3.0
|
// WilsonGaugeAction with beta 3.0
|
||||||
assert(epsilon > 0.0);
|
assert(epsilon > 0.0);
|
||||||
|
assert(tolerance > 0.0);
|
||||||
LogMessage();
|
LogMessage();
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -64,6 +67,8 @@ public:
|
|||||||
<< "[WilsonFlow] epsilon : " << epsilon << std::endl;
|
<< "[WilsonFlow] epsilon : " << epsilon << std::endl;
|
||||||
std::cout << GridLogMessage
|
std::cout << GridLogMessage
|
||||||
<< "[WilsonFlow] full trajectory : " << Nstep * epsilon << std::endl;
|
<< "[WilsonFlow] full trajectory : " << Nstep * epsilon << std::endl;
|
||||||
|
std::cout << GridLogMessage
|
||||||
|
<< "[WilsonFlow] tolerance : " << tolerance << std::endl;
|
||||||
}
|
}
|
||||||
|
|
||||||
virtual void smear(GaugeField&, const GaugeField&) const;
|
virtual void smear(GaugeField&, const GaugeField&) const;
|
||||||
@ -106,11 +111,14 @@ void WilsonFlow<Gimpl>::evolve_step_adaptive(typename Gimpl::GaugeField &U, Real
|
|||||||
if (maxTau - taus < epsilon){
|
if (maxTau - taus < epsilon){
|
||||||
epsilon = maxTau-taus;
|
epsilon = maxTau-taus;
|
||||||
}
|
}
|
||||||
//std::cout << GridLogMessage << "Integration epsilon : " << epsilon << std::endl;
|
std::cout << GridLogMessage << "Integration epsilon : " << epsilon << std::endl;
|
||||||
GaugeField Z(U.Grid());
|
GaugeField Z(U.Grid());
|
||||||
GaugeField Zprime(U.Grid());
|
GaugeField Zprime(U.Grid());
|
||||||
GaugeField tmp(U.Grid()), Uprime(U.Grid());
|
GaugeField tmp(U.Grid()), Uprime(U.Grid()),Usave(U.Grid());
|
||||||
|
|
||||||
Uprime = U;
|
Uprime = U;
|
||||||
|
Usave = U;
|
||||||
|
|
||||||
SG.deriv(U, Z);
|
SG.deriv(U, Z);
|
||||||
Zprime = -Z;
|
Zprime = -Z;
|
||||||
Z *= 0.25; // Z0 = 1/4 * F(U)
|
Z *= 0.25; // Z0 = 1/4 * F(U)
|
||||||
@ -128,18 +136,33 @@ void WilsonFlow<Gimpl>::evolve_step_adaptive(typename Gimpl::GaugeField &U, Real
|
|||||||
Z *= 3.0/4.0; // Z = 17/36*Z0 -8/9*Z1 +3/4*Z2
|
Z *= 3.0/4.0; // Z = 17/36*Z0 -8/9*Z1 +3/4*Z2
|
||||||
Gimpl::update_field(Z, U, -2.0*epsilon); // V(t+e) = exp(ep*Z)*W2
|
Gimpl::update_field(Z, U, -2.0*epsilon); // V(t+e) = exp(ep*Z)*W2
|
||||||
|
|
||||||
// Ramos
|
// Ramos arXiv:1301.4388
|
||||||
Gimpl::update_field(Zprime, Uprime, -2.0*epsilon); // V'(t+e) = exp(ep*Z')*W0
|
Gimpl::update_field(Zprime, Uprime, -2.0*epsilon); // V'(t+e) = exp(ep*Z')*W0
|
||||||
// Compute distance as norm^2 of the difference
|
// Compute distance as norm^2 of the difference
|
||||||
GaugeField diffU = U - Uprime;
|
GaugeField diffU = U - Uprime;
|
||||||
RealD diff = norm2(diffU);
|
// Wrong
|
||||||
// adjust integration step
|
// RealD diff = norm2(diffU);
|
||||||
|
// std::cout << GridLogMessage << "norm2: " << diff << std::endl;
|
||||||
|
|
||||||
|
// RealD tol=1e-3;
|
||||||
|
|
||||||
|
RealD diff = real(rankInnerMax(diffU,diffU));
|
||||||
|
diff = sqrt(diff)/18.; // distance defined in Ramos
|
||||||
|
|
||||||
|
GridBase *grid = diffU.Grid();
|
||||||
|
std::cout << GridLogMessage << "max: " << diff << std::endl;
|
||||||
|
grid->GlobalMax(diff);
|
||||||
|
std::cout << GridLogMessage << "max: " << diff << std::endl;
|
||||||
|
|
||||||
|
if(diff < tolerance) {
|
||||||
taus += epsilon;
|
taus += epsilon;
|
||||||
//std::cout << GridLogMessage << "Adjusting integration step with distance: " << diff << std::endl;
|
// std::cout << GridLogMessage << "Adjusting integration step with distance: " << diff << std::endl;
|
||||||
|
} else {
|
||||||
|
U = Usave;
|
||||||
|
}
|
||||||
|
|
||||||
epsilon = epsilon*0.95*std::pow(1e-4/diff,1./3.);
|
epsilon = epsilon*0.95*std::pow(tolerance/diff,1./3.);
|
||||||
//std::cout << GridLogMessage << "New epsilon : " << epsilon << std::endl;
|
std::cout << GridLogMessage << "Distance : "<<diff<<"New epsilon : " << epsilon << std::endl;
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -184,8 +207,11 @@ void WilsonFlow<Gimpl>::smear(GaugeField& out, const GaugeField& in) const {
|
|||||||
template <class Gimpl>
|
template <class Gimpl>
|
||||||
void WilsonFlow<Gimpl>::smear_adaptive(GaugeField& out, const GaugeField& in, RealD maxTau){
|
void WilsonFlow<Gimpl>::smear_adaptive(GaugeField& out, const GaugeField& in, RealD maxTau){
|
||||||
out = in;
|
out = in;
|
||||||
taus = epsilon;
|
// taus = epsilon;
|
||||||
|
taus = 0.;
|
||||||
unsigned int step = 0;
|
unsigned int step = 0;
|
||||||
|
double measTau = epsilon*measure_interval;
|
||||||
|
std::cout << GridLogMessage << "measTau :"<< measTau << std::endl;
|
||||||
do{
|
do{
|
||||||
step++;
|
step++;
|
||||||
//std::cout << GridLogMessage << "Evolution time :"<< taus << std::endl;
|
//std::cout << GridLogMessage << "Evolution time :"<< taus << std::endl;
|
||||||
@ -193,10 +219,12 @@ void WilsonFlow<Gimpl>::smear_adaptive(GaugeField& out, const GaugeField& in, Re
|
|||||||
std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : "
|
std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : "
|
||||||
<< step << " " << taus << " "
|
<< step << " " << taus << " "
|
||||||
<< energyDensityPlaquette(out) << std::endl;
|
<< energyDensityPlaquette(out) << std::endl;
|
||||||
if( step % measure_interval == 0){
|
// if( step % measure_interval == 0){
|
||||||
|
if( taus > measTau ) {
|
||||||
std::cout << GridLogMessage << "[WilsonFlow] Top. charge : "
|
std::cout << GridLogMessage << "[WilsonFlow] Top. charge : "
|
||||||
<< step << " "
|
<< step << " "
|
||||||
<< WilsonLoops<PeriodicGimplR>::TopologicalCharge(out) << std::endl;
|
<< WilsonLoops<PeriodicGimplR>::TopologicalCharge(out) << std::endl;
|
||||||
|
measTau += epsilon*measure_interval;
|
||||||
}
|
}
|
||||||
} while (taus < maxTau);
|
} while (taus < maxTau);
|
||||||
|
|
||||||
|
@ -1,35 +0,0 @@
|
|||||||
/*************************************************************************************
|
|
||||||
|
|
||||||
Grid physics library, www.github.com/paboyle/Grid
|
|
||||||
|
|
||||||
Source file: ./lib/serialisation/BaseIO.h
|
|
||||||
|
|
||||||
Copyright (C) 2015
|
|
||||||
|
|
||||||
Author: Michael Marshall <michael.marshall@ed.ac.uk>
|
|
||||||
|
|
||||||
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
|
|
||||||
the Free Software Foundation; either version 2 of the License, or
|
|
||||||
(at your option) any later version.
|
|
||||||
|
|
||||||
This program is distributed in the hope that it will be useful,
|
|
||||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
||||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
||||||
GNU General Public License for more details.
|
|
||||||
|
|
||||||
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.,
|
|
||||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
|
||||||
|
|
||||||
See the full license in the file "LICENSE" in the top level distribution directory
|
|
||||||
*************************************************************************************/
|
|
||||||
/* END LEGAL */
|
|
||||||
|
|
||||||
#include <Grid/GridCore.h>
|
|
||||||
|
|
||||||
NAMESPACE_BEGIN(Grid)
|
|
||||||
|
|
||||||
std::uint64_t EigenIO::EigenResizeCounter(0);
|
|
||||||
|
|
||||||
NAMESPACE_END(Grid)
|
|
@ -9,7 +9,6 @@
|
|||||||
Author: Antonin Portelli <antonin.portelli@me.com>
|
Author: Antonin Portelli <antonin.portelli@me.com>
|
||||||
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||||
Author: Guido Cossu <guido.cossu@ed.ac.uk>
|
Author: Guido Cossu <guido.cossu@ed.ac.uk>
|
||||||
Author: Michael Marshall <michael.marshall@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
|
||||||
@ -31,7 +30,6 @@ Author: Michael Marshall <michael.marshall@ed.ac.uk>
|
|||||||
#ifndef GRID_SERIALISATION_ABSTRACT_READER_H
|
#ifndef GRID_SERIALISATION_ABSTRACT_READER_H
|
||||||
#define GRID_SERIALISATION_ABSTRACT_READER_H
|
#define GRID_SERIALISATION_ABSTRACT_READER_H
|
||||||
|
|
||||||
#include <atomic>
|
|
||||||
#include <type_traits>
|
#include <type_traits>
|
||||||
#include <Grid/tensors/Tensors.h>
|
#include <Grid/tensors/Tensors.h>
|
||||||
#include <Grid/serialisation/VectorUtils.h>
|
#include <Grid/serialisation/VectorUtils.h>
|
||||||
@ -112,10 +110,6 @@ namespace Grid {
|
|||||||
template <typename ET>
|
template <typename ET>
|
||||||
inline typename std::enable_if<is_tensor_of_container<ET>::value, typename Traits<ET>::scalar_type *>::type
|
inline typename std::enable_if<is_tensor_of_container<ET>::value, typename Traits<ET>::scalar_type *>::type
|
||||||
getFirstScalar(ET &eigenTensor) { return eigenTensor.data()->begin(); }
|
getFirstScalar(ET &eigenTensor) { return eigenTensor.data()->begin(); }
|
||||||
|
|
||||||
// Counter for resized EigenTensors (poor man's substitute for allocator)
|
|
||||||
// Defined in BinaryIO.cc
|
|
||||||
extern std::uint64_t EigenResizeCounter;
|
|
||||||
}
|
}
|
||||||
|
|
||||||
// Abstract writer/reader classes ////////////////////////////////////////////
|
// Abstract writer/reader classes ////////////////////////////////////////////
|
||||||
@ -503,14 +497,8 @@ namespace Grid {
|
|||||||
typename std::enable_if<EigenIO::is_tensor_variable<ETensor>::value, void>::type
|
typename std::enable_if<EigenIO::is_tensor_variable<ETensor>::value, void>::type
|
||||||
Reader<T>::Reshape(ETensor &t, const std::array<typename ETensor::Index, ETensor::NumDimensions> &dims )
|
Reader<T>::Reshape(ETensor &t, const std::array<typename ETensor::Index, ETensor::NumDimensions> &dims )
|
||||||
{
|
{
|
||||||
#ifdef GRID_OMP
|
|
||||||
// The memory counter is the reason this must be done from the primary thread
|
|
||||||
assert(omp_in_parallel()==0 && "Deserialisation which resizes Eigen tensor must happen from primary thread");
|
|
||||||
#endif
|
|
||||||
EigenIO::EigenResizeCounter -= static_cast<uint64_t>(t.size()) * sizeof(typename ETensor::Scalar);
|
|
||||||
//t.reshape( dims );
|
//t.reshape( dims );
|
||||||
t.resize( dims );
|
t.resize( dims );
|
||||||
EigenIO::EigenResizeCounter += static_cast<uint64_t>(t.size()) * sizeof(typename ETensor::Scalar);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
template <typename T>
|
template <typename T>
|
||||||
|
@ -1,34 +1,3 @@
|
|||||||
/*************************************************************************************
|
|
||||||
|
|
||||||
Grid physics library, www.github.com/paboyle/Grid
|
|
||||||
|
|
||||||
Source file: ./Grid/serialisation/VectorUtils.h
|
|
||||||
|
|
||||||
Copyright (C) 2015
|
|
||||||
|
|
||||||
Author: Antonin Portelli <antonin.portelli@me.com>
|
|
||||||
Author: Peter Boyle <paboyle@ed.ac.uk>
|
|
||||||
Author: Guido Cossu <guido.cossu@ed.ac.uk>
|
|
||||||
Author: Michael Marshall <michael.marshall@ed.ac.uk>
|
|
||||||
|
|
||||||
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
|
|
||||||
the Free Software Foundation; either version 2 of the License, or
|
|
||||||
(at your option) any later version.
|
|
||||||
|
|
||||||
This program is distributed in the hope that it will be useful,
|
|
||||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
||||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
||||||
GNU General Public License for more details.
|
|
||||||
|
|
||||||
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.,
|
|
||||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
|
||||||
|
|
||||||
See the full license in the file "LICENSE" in the top level distribution directory
|
|
||||||
*************************************************************************************/
|
|
||||||
/* END LEGAL */
|
|
||||||
|
|
||||||
#include <Grid/Grid.h>
|
#include <Grid/Grid.h>
|
||||||
|
|
||||||
using namespace Grid;
|
using namespace Grid;
|
||||||
|
@ -1,34 +1,3 @@
|
|||||||
/*************************************************************************************
|
|
||||||
|
|
||||||
Grid physics library, www.github.com/paboyle/Grid
|
|
||||||
|
|
||||||
Source file: ./Grid/serialisation/VectorUtils.h
|
|
||||||
|
|
||||||
Copyright (C) 2015
|
|
||||||
|
|
||||||
Author: Peter Boyle <paboyle@ed.ac.uk>
|
|
||||||
Author: Antonin Portelli <antonin.portelli@me.com>
|
|
||||||
Author: Guido Cossu <guido.cossu@ed.ac.uk>
|
|
||||||
Author: Michael Marshall <michael.marshall@ed.ac.uk>
|
|
||||||
|
|
||||||
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
|
|
||||||
the Free Software Foundation; either version 2 of the License, or
|
|
||||||
(at your option) any later version.
|
|
||||||
|
|
||||||
This program is distributed in the hope that it will be useful,
|
|
||||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
||||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
||||||
GNU General Public License for more details.
|
|
||||||
|
|
||||||
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.,
|
|
||||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
|
||||||
|
|
||||||
See the full license in the file "LICENSE" in the top level distribution directory
|
|
||||||
*************************************************************************************/
|
|
||||||
/* END LEGAL */
|
|
||||||
|
|
||||||
#ifndef GRID_SERIALISATION_HDF5_H
|
#ifndef GRID_SERIALISATION_HDF5_H
|
||||||
#define GRID_SERIALISATION_HDF5_H
|
#define GRID_SERIALISATION_HDF5_H
|
||||||
|
|
||||||
@ -65,13 +34,11 @@ namespace Grid
|
|||||||
template <typename U>
|
template <typename U>
|
||||||
void writeDefault(const std::string &s, const U &x);
|
void writeDefault(const std::string &s, const U &x);
|
||||||
template <typename U>
|
template <typename U>
|
||||||
void writeRagged(const std::string &s, const std::vector<U> &x);
|
typename std::enable_if<element<std::vector<U>>::is_number, void>::type
|
||||||
template <typename U>
|
|
||||||
typename std::enable_if<is_flattenable<std::vector<U>>::value>::type
|
|
||||||
writeDefault(const std::string &s, const std::vector<U> &x);
|
writeDefault(const std::string &s, const std::vector<U> &x);
|
||||||
template <typename U>
|
template <typename U>
|
||||||
typename std::enable_if<!is_flattenable<std::vector<U>>::value>::type
|
typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
|
||||||
writeDefault(const std::string &s, const std::vector<U> &x) { writeRagged(s, x); }
|
writeDefault(const std::string &s, const std::vector<U> &x);
|
||||||
template <typename U>
|
template <typename U>
|
||||||
void writeMultiDim(const std::string &s, const std::vector<size_t> & Dimensions, const U * pDataRowMajor, size_t NumElements);
|
void writeMultiDim(const std::string &s, const std::vector<size_t> & Dimensions, const U * pDataRowMajor, size_t NumElements);
|
||||||
H5NS::Group & getGroup(void);
|
H5NS::Group & getGroup(void);
|
||||||
@ -97,13 +64,11 @@ namespace Grid
|
|||||||
template <typename U>
|
template <typename U>
|
||||||
void readDefault(const std::string &s, U &output);
|
void readDefault(const std::string &s, U &output);
|
||||||
template <typename U>
|
template <typename U>
|
||||||
void readRagged(const std::string &s, std::vector<U> &x);
|
typename std::enable_if<element<std::vector<U>>::is_number, void>::type
|
||||||
template <typename U>
|
|
||||||
typename std::enable_if<is_flattenable<std::vector<U>>::value>::type
|
|
||||||
readDefault(const std::string &s, std::vector<U> &x);
|
readDefault(const std::string &s, std::vector<U> &x);
|
||||||
template <typename U>
|
template <typename U>
|
||||||
typename std::enable_if<!is_flattenable<std::vector<U>>::value>::type
|
typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
|
||||||
readDefault(const std::string &s, std::vector<U> &x) { readRagged(s, x); }
|
readDefault(const std::string &s, std::vector<U> &x);
|
||||||
template <typename U>
|
template <typename U>
|
||||||
void readMultiDim(const std::string &s, std::vector<U> &buf, std::vector<size_t> &dim);
|
void readMultiDim(const std::string &s, std::vector<U> &buf, std::vector<size_t> &dim);
|
||||||
H5NS::Group & getGroup(void);
|
H5NS::Group & getGroup(void);
|
||||||
@ -211,13 +176,11 @@ namespace Grid
|
|||||||
}
|
}
|
||||||
|
|
||||||
template <typename U>
|
template <typename U>
|
||||||
typename std::enable_if<is_flattenable<std::vector<U>>::value>::type
|
typename std::enable_if<element<std::vector<U>>::is_number, void>::type
|
||||||
Hdf5Writer::writeDefault(const std::string &s, const std::vector<U> &x)
|
Hdf5Writer::writeDefault(const std::string &s, const std::vector<U> &x)
|
||||||
{
|
|
||||||
if (isRegularShape(x))
|
|
||||||
{
|
{
|
||||||
// alias to element type
|
// alias to element type
|
||||||
using Scalar = typename is_flattenable<std::vector<U>>::type;
|
typedef typename element<std::vector<U>>::type Element;
|
||||||
|
|
||||||
// flatten the vector and getting dimensions
|
// flatten the vector and getting dimensions
|
||||||
Flatten<std::vector<U>> flat(x);
|
Flatten<std::vector<U>> flat(x);
|
||||||
@ -225,16 +188,12 @@ namespace Grid
|
|||||||
const auto &flatx = flat.getFlatVector();
|
const auto &flatx = flat.getFlatVector();
|
||||||
for (auto &d: flat.getDim())
|
for (auto &d: flat.getDim())
|
||||||
dim.push_back(d);
|
dim.push_back(d);
|
||||||
writeMultiDim<Scalar>(s, dim, &flatx[0], flatx.size());
|
writeMultiDim<Element>(s, dim, &flatx[0], flatx.size());
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
writeRagged(s, x);
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
|
|
||||||
template <typename U>
|
template <typename U>
|
||||||
void Hdf5Writer::writeRagged(const std::string &s, const std::vector<U> &x)
|
typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
|
||||||
|
Hdf5Writer::writeDefault(const std::string &s, const std::vector<U> &x)
|
||||||
{
|
{
|
||||||
push(s);
|
push(s);
|
||||||
writeSingleAttribute(x.size(), HDF5_GRID_GUARD "vector_size",
|
writeSingleAttribute(x.size(), HDF5_GRID_GUARD "vector_size",
|
||||||
@ -270,7 +229,7 @@ namespace Grid
|
|||||||
void Hdf5Reader::readMultiDim(const std::string &s, std::vector<U> &buf, std::vector<size_t> &dim)
|
void Hdf5Reader::readMultiDim(const std::string &s, std::vector<U> &buf, std::vector<size_t> &dim)
|
||||||
{
|
{
|
||||||
// alias to element type
|
// alias to element type
|
||||||
using Scalar = typename is_flattenable<std::vector<U>>::type;
|
typedef typename element<std::vector<U>>::type Element;
|
||||||
|
|
||||||
// read the dimensions
|
// read the dimensions
|
||||||
H5NS::DataSpace dataSpace;
|
H5NS::DataSpace dataSpace;
|
||||||
@ -301,33 +260,26 @@ namespace Grid
|
|||||||
H5NS::DataSet dataSet;
|
H5NS::DataSet dataSet;
|
||||||
|
|
||||||
dataSet = group_.openDataSet(s);
|
dataSet = group_.openDataSet(s);
|
||||||
dataSet.read(buf.data(), Hdf5Type<Scalar>::type());
|
dataSet.read(buf.data(), Hdf5Type<Element>::type());
|
||||||
}
|
}
|
||||||
else
|
else
|
||||||
{
|
{
|
||||||
H5NS::Attribute attribute;
|
H5NS::Attribute attribute;
|
||||||
|
|
||||||
attribute = group_.openAttribute(s);
|
attribute = group_.openAttribute(s);
|
||||||
attribute.read(Hdf5Type<Scalar>::type(), buf.data());
|
attribute.read(Hdf5Type<Element>::type(), buf.data());
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
template <typename U>
|
template <typename U>
|
||||||
typename std::enable_if<is_flattenable<std::vector<U>>::value>::type
|
typename std::enable_if<element<std::vector<U>>::is_number, void>::type
|
||||||
Hdf5Reader::readDefault(const std::string &s, std::vector<U> &x)
|
Hdf5Reader::readDefault(const std::string &s, std::vector<U> &x)
|
||||||
{
|
|
||||||
if (H5Lexists (group_.getId(), s.c_str(), H5P_DEFAULT) > 0
|
|
||||||
&& H5Aexists_by_name(group_.getId(), s.c_str(), HDF5_GRID_GUARD "vector_size", H5P_DEFAULT ) > 0)
|
|
||||||
{
|
|
||||||
readRagged(s, x);
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
{
|
||||||
// alias to element type
|
// alias to element type
|
||||||
using Scalar = typename is_flattenable<std::vector<U>>::type;
|
typedef typename element<std::vector<U>>::type Element;
|
||||||
|
|
||||||
std::vector<size_t> dim;
|
std::vector<size_t> dim;
|
||||||
std::vector<Scalar> buf;
|
std::vector<Element> buf;
|
||||||
readMultiDim( s, buf, dim );
|
readMultiDim( s, buf, dim );
|
||||||
|
|
||||||
// reconstruct the multidimensional vector
|
// reconstruct the multidimensional vector
|
||||||
@ -335,10 +287,10 @@ namespace Grid
|
|||||||
|
|
||||||
x = r.getVector();
|
x = r.getVector();
|
||||||
}
|
}
|
||||||
}
|
|
||||||
|
|
||||||
template <typename U>
|
template <typename U>
|
||||||
void Hdf5Reader::readRagged(const std::string &s, std::vector<U> &x)
|
typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
|
||||||
|
Hdf5Reader::readDefault(const std::string &s, std::vector<U> &x)
|
||||||
{
|
{
|
||||||
uint64_t size;
|
uint64_t size;
|
||||||
|
|
||||||
|
@ -118,13 +118,13 @@ static inline std::string SerialisableClassName(void) {return std::string(#cname
|
|||||||
static constexpr bool isEnum = false; \
|
static constexpr bool isEnum = false; \
|
||||||
GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_MEMBER,__VA_ARGS__))\
|
GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_MEMBER,__VA_ARGS__))\
|
||||||
template <typename T>\
|
template <typename T>\
|
||||||
static inline void write(::Grid::Writer<T> &WR,const std::string &s, const cname &obj){ \
|
static inline void write(Writer<T> &WR,const std::string &s, const cname &obj){ \
|
||||||
push(WR,s);\
|
push(WR,s);\
|
||||||
GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_WRITE_MEMBER,__VA_ARGS__)) \
|
GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_WRITE_MEMBER,__VA_ARGS__)) \
|
||||||
pop(WR);\
|
pop(WR);\
|
||||||
}\
|
}\
|
||||||
template <typename T>\
|
template <typename T>\
|
||||||
static inline void read(::Grid::Reader<T> &RD,const std::string &s, cname &obj){ \
|
static inline void read(Reader<T> &RD,const std::string &s, cname &obj){ \
|
||||||
if (!push(RD,s))\
|
if (!push(RD,s))\
|
||||||
{\
|
{\
|
||||||
std::cout << ::Grid::GridLogWarning << "IO: Cannot open node '" << s << "'" << std::endl; \
|
std::cout << ::Grid::GridLogWarning << "IO: Cannot open node '" << s << "'" << std::endl; \
|
||||||
|
@ -9,7 +9,6 @@
|
|||||||
Author: Antonin Portelli <antonin.portelli@me.com>
|
Author: Antonin Portelli <antonin.portelli@me.com>
|
||||||
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>
|
||||||
Author: Michael Marshall <michael.marshall@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
|
||||||
@ -237,36 +236,21 @@ namespace Grid {
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
// is_flattenable<T>::value is true if T is a std::vector<> which can be flattened //////////////////////
|
// Vector element trait //////////////////////////////////////////////////////
|
||||||
template <typename T, typename V = void>
|
template <typename T>
|
||||||
struct is_flattenable : std::false_type
|
struct element
|
||||||
{
|
{
|
||||||
using type = T;
|
typedef T type;
|
||||||
using grid_type = T;
|
static constexpr bool is_number = false;
|
||||||
static constexpr int vecRank = 0;
|
|
||||||
static constexpr bool isGridTensor = false;
|
|
||||||
static constexpr bool children_flattenable = std::is_arithmetic<T>::value or is_complex<T>::value;
|
|
||||||
};
|
};
|
||||||
|
|
||||||
template <typename T>
|
template <typename T>
|
||||||
struct is_flattenable<T, typename std::enable_if<isGridTensor<T>::value>::type> : std::false_type
|
struct element<std::vector<T>>
|
||||||
{
|
{
|
||||||
using type = typename GridTypeMapper<T>::scalar_type;
|
typedef typename element<T>::type type;
|
||||||
using grid_type = T;
|
static constexpr bool is_number = std::is_arithmetic<T>::value
|
||||||
static constexpr int vecRank = 0;
|
or is_complex<T>::value
|
||||||
static constexpr bool isGridTensor = true;
|
or element<T>::is_number;
|
||||||
static constexpr bool children_flattenable = true;
|
|
||||||
};
|
|
||||||
|
|
||||||
template <typename T>
|
|
||||||
struct is_flattenable<std::vector<T>, typename std::enable_if<is_flattenable<T>::children_flattenable>::type>
|
|
||||||
: std::true_type
|
|
||||||
{
|
|
||||||
using type = typename is_flattenable<T>::type;
|
|
||||||
using grid_type = typename is_flattenable<T>::grid_type;
|
|
||||||
static constexpr bool isGridTensor = is_flattenable<T>::isGridTensor;
|
|
||||||
static constexpr int vecRank = is_flattenable<T>::vecRank + 1;
|
|
||||||
static constexpr bool children_flattenable = true;
|
|
||||||
};
|
};
|
||||||
|
|
||||||
// Vector flattening utility class ////////////////////////////////////////////
|
// Vector flattening utility class ////////////////////////////////////////////
|
||||||
@ -275,29 +259,22 @@ namespace Grid {
|
|||||||
class Flatten
|
class Flatten
|
||||||
{
|
{
|
||||||
public:
|
public:
|
||||||
using Scalar = typename is_flattenable<V>::type;
|
typedef typename element<V>::type Element;
|
||||||
static constexpr bool isGridTensor = is_flattenable<V>::isGridTensor;
|
|
||||||
public:
|
public:
|
||||||
explicit Flatten(const V &vector);
|
explicit Flatten(const V &vector);
|
||||||
const V & getVector(void) const { return vector_; }
|
const V & getVector(void);
|
||||||
const std::vector<Scalar> & getFlatVector(void) const { return flatVector_; }
|
const std::vector<Element> & getFlatVector(void);
|
||||||
const std::vector<size_t> & getDim(void) const { return dim_; }
|
const std::vector<size_t> & getDim(void);
|
||||||
private:
|
private:
|
||||||
template <typename W> typename std::enable_if<!is_flattenable<W>::value && !is_flattenable<W>::isGridTensor>::type
|
void accumulate(const Element &e);
|
||||||
accumulate(const W &e);
|
template <typename W>
|
||||||
template <typename W> typename std::enable_if<!is_flattenable<W>::value && is_flattenable<W>::isGridTensor>::type
|
void accumulate(const W &v);
|
||||||
accumulate(const W &e);
|
void accumulateDim(const Element &e);
|
||||||
template <typename W> typename std::enable_if< is_flattenable<W>::value>::type
|
template <typename W>
|
||||||
accumulate(const W &v);
|
void accumulateDim(const W &v);
|
||||||
template <typename W> typename std::enable_if<!is_flattenable<W>::value && !is_flattenable<W>::isGridTensor>::type
|
|
||||||
accumulateDim(const W &e) {} // Innermost is a scalar - do nothing
|
|
||||||
template <typename W> typename std::enable_if<!is_flattenable<W>::value && is_flattenable<W>::isGridTensor>::type
|
|
||||||
accumulateDim(const W &e);
|
|
||||||
template <typename W> typename std::enable_if< is_flattenable<W>::value>::type
|
|
||||||
accumulateDim(const W &v);
|
|
||||||
private:
|
private:
|
||||||
const V &vector_;
|
const V &vector_;
|
||||||
std::vector<Scalar> flatVector_;
|
std::vector<Element> flatVector_;
|
||||||
std::vector<size_t> dim_;
|
std::vector<size_t> dim_;
|
||||||
};
|
};
|
||||||
|
|
||||||
@ -306,32 +283,23 @@ namespace Grid {
|
|||||||
class Reconstruct
|
class Reconstruct
|
||||||
{
|
{
|
||||||
public:
|
public:
|
||||||
using Scalar = typename is_flattenable<V>::type;
|
typedef typename element<V>::type Element;
|
||||||
static constexpr bool isGridTensor = is_flattenable<V>::isGridTensor;
|
|
||||||
public:
|
public:
|
||||||
Reconstruct(const std::vector<Scalar> &flatVector,
|
Reconstruct(const std::vector<Element> &flatVector,
|
||||||
const std::vector<size_t> &dim);
|
const std::vector<size_t> &dim);
|
||||||
const V & getVector(void) const { return vector_; }
|
const V & getVector(void);
|
||||||
const std::vector<Scalar> & getFlatVector(void) const { return flatVector_; }
|
const std::vector<Element> & getFlatVector(void);
|
||||||
const std::vector<size_t> & getDim(void) const { return dim_; }
|
const std::vector<size_t> & getDim(void);
|
||||||
private:
|
private:
|
||||||
template <typename W> typename std::enable_if<!is_flattenable<W>::value && !is_flattenable<W>::isGridTensor>::type
|
void fill(std::vector<Element> &v);
|
||||||
fill(W &v);
|
template <typename W>
|
||||||
template <typename W> typename std::enable_if<!is_flattenable<W>::value && is_flattenable<W>::isGridTensor>::type
|
void fill(W &v);
|
||||||
fill(W &v);
|
void resize(std::vector<Element> &v, const unsigned int dim);
|
||||||
template <typename W> typename std::enable_if< is_flattenable<W>::value>::type
|
template <typename W>
|
||||||
fill(W &v);
|
void resize(W &v, const unsigned int dim);
|
||||||
template <typename W> typename std::enable_if< is_flattenable<W>::value && is_flattenable<W>::vecRank==1>::type
|
|
||||||
resize(W &v, const unsigned int dim);
|
|
||||||
template <typename W> typename std::enable_if< is_flattenable<W>::value && (is_flattenable<W>::vecRank>1)>::type
|
|
||||||
resize(W &v, const unsigned int dim);
|
|
||||||
template <typename W> typename std::enable_if<!is_flattenable<W>::isGridTensor>::type
|
|
||||||
checkInnermost(const W &e) {} // Innermost is a scalar - do nothing
|
|
||||||
template <typename W> typename std::enable_if< is_flattenable<W>::isGridTensor>::type
|
|
||||||
checkInnermost(const W &e);
|
|
||||||
private:
|
private:
|
||||||
V vector_;
|
V vector_;
|
||||||
const std::vector<Scalar> &flatVector_;
|
const std::vector<Element> &flatVector_;
|
||||||
std::vector<size_t> dim_;
|
std::vector<size_t> dim_;
|
||||||
size_t ind_{0};
|
size_t ind_{0};
|
||||||
unsigned int dimInd_{0};
|
unsigned int dimInd_{0};
|
||||||
@ -339,24 +307,14 @@ namespace Grid {
|
|||||||
|
|
||||||
// Flatten class template implementation
|
// Flatten class template implementation
|
||||||
template <typename V>
|
template <typename V>
|
||||||
template <typename W> typename std::enable_if<!is_flattenable<W>::value && !is_flattenable<W>::isGridTensor>::type
|
void Flatten<V>::accumulate(const Element &e)
|
||||||
Flatten<V>::accumulate(const W &e)
|
|
||||||
{
|
{
|
||||||
flatVector_.push_back(e);
|
flatVector_.push_back(e);
|
||||||
}
|
}
|
||||||
|
|
||||||
template <typename V>
|
template <typename V>
|
||||||
template <typename W> typename std::enable_if<!is_flattenable<W>::value && is_flattenable<W>::isGridTensor>::type
|
template <typename W>
|
||||||
Flatten<V>::accumulate(const W &e)
|
void Flatten<V>::accumulate(const W &v)
|
||||||
{
|
|
||||||
for (const Scalar &x: e) {
|
|
||||||
flatVector_.push_back(x);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
template <typename V>
|
|
||||||
template <typename W> typename std::enable_if<is_flattenable<W>::value>::type
|
|
||||||
Flatten<V>::accumulate(const W &v)
|
|
||||||
{
|
{
|
||||||
for (auto &e: v)
|
for (auto &e: v)
|
||||||
{
|
{
|
||||||
@ -365,17 +323,11 @@ namespace Grid {
|
|||||||
}
|
}
|
||||||
|
|
||||||
template <typename V>
|
template <typename V>
|
||||||
template <typename W> typename std::enable_if<!is_flattenable<W>::value && is_flattenable<W>::isGridTensor>::type
|
void Flatten<V>::accumulateDim(const Element &e) {};
|
||||||
Flatten<V>::accumulateDim(const W &e)
|
|
||||||
{
|
|
||||||
using Traits = GridTypeMapper<typename is_flattenable<W>::grid_type>;
|
|
||||||
for (int rank=0; rank < Traits::Rank; ++rank)
|
|
||||||
dim_.push_back(Traits::Dimension(rank));
|
|
||||||
}
|
|
||||||
|
|
||||||
template <typename V>
|
template <typename V>
|
||||||
template <typename W> typename std::enable_if<is_flattenable<W>::value>::type
|
template <typename W>
|
||||||
Flatten<V>::accumulateDim(const W &v)
|
void Flatten<V>::accumulateDim(const W &v)
|
||||||
{
|
{
|
||||||
dim_.push_back(v.size());
|
dim_.push_back(v.size());
|
||||||
accumulateDim(v[0]);
|
accumulateDim(v[0]);
|
||||||
@ -385,26 +337,32 @@ namespace Grid {
|
|||||||
Flatten<V>::Flatten(const V &vector)
|
Flatten<V>::Flatten(const V &vector)
|
||||||
: vector_(vector)
|
: vector_(vector)
|
||||||
{
|
{
|
||||||
accumulateDim(vector_);
|
|
||||||
std::size_t TotalSize{ dim_[0] };
|
|
||||||
for (int i = 1; i < dim_.size(); ++i) {
|
|
||||||
TotalSize *= dim_[i];
|
|
||||||
}
|
|
||||||
flatVector_.reserve(TotalSize);
|
|
||||||
accumulate(vector_);
|
accumulate(vector_);
|
||||||
|
accumulateDim(vector_);
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename V>
|
||||||
|
const V & Flatten<V>::getVector(void)
|
||||||
|
{
|
||||||
|
return vector_;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename V>
|
||||||
|
const std::vector<typename Flatten<V>::Element> &
|
||||||
|
Flatten<V>::getFlatVector(void)
|
||||||
|
{
|
||||||
|
return flatVector_;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename V>
|
||||||
|
const std::vector<size_t> & Flatten<V>::getDim(void)
|
||||||
|
{
|
||||||
|
return dim_;
|
||||||
}
|
}
|
||||||
|
|
||||||
// Reconstruct class template implementation
|
// Reconstruct class template implementation
|
||||||
template <typename V>
|
template <typename V>
|
||||||
template <typename W> typename std::enable_if<!is_flattenable<W>::value && !is_flattenable<W>::isGridTensor>::type
|
void Reconstruct<V>::fill(std::vector<Element> &v)
|
||||||
Reconstruct<V>::fill(W &v)
|
|
||||||
{
|
|
||||||
v = flatVector_[ind_++];
|
|
||||||
}
|
|
||||||
|
|
||||||
template <typename V>
|
|
||||||
template <typename W> typename std::enable_if<!is_flattenable<W>::value && is_flattenable<W>::isGridTensor>::type
|
|
||||||
Reconstruct<V>::fill(W &v)
|
|
||||||
{
|
{
|
||||||
for (auto &e: v)
|
for (auto &e: v)
|
||||||
{
|
{
|
||||||
@ -413,8 +371,8 @@ namespace Grid {
|
|||||||
}
|
}
|
||||||
|
|
||||||
template <typename V>
|
template <typename V>
|
||||||
template <typename W> typename std::enable_if<is_flattenable<W>::value>::type
|
template <typename W>
|
||||||
Reconstruct<V>::fill(W &v)
|
void Reconstruct<V>::fill(W &v)
|
||||||
{
|
{
|
||||||
for (auto &e: v)
|
for (auto &e: v)
|
||||||
{
|
{
|
||||||
@ -423,15 +381,14 @@ namespace Grid {
|
|||||||
}
|
}
|
||||||
|
|
||||||
template <typename V>
|
template <typename V>
|
||||||
template <typename W> typename std::enable_if<is_flattenable<W>::value && is_flattenable<W>::vecRank==1>::type
|
void Reconstruct<V>::resize(std::vector<Element> &v, const unsigned int dim)
|
||||||
Reconstruct<V>::resize(W &v, const unsigned int dim)
|
|
||||||
{
|
{
|
||||||
v.resize(dim_[dim]);
|
v.resize(dim_[dim]);
|
||||||
}
|
}
|
||||||
|
|
||||||
template <typename V>
|
template <typename V>
|
||||||
template <typename W> typename std::enable_if<is_flattenable<W>::value && (is_flattenable<W>::vecRank>1)>::type
|
template <typename W>
|
||||||
Reconstruct<V>::resize(W &v, const unsigned int dim)
|
void Reconstruct<V>::resize(W &v, const unsigned int dim)
|
||||||
{
|
{
|
||||||
v.resize(dim_[dim]);
|
v.resize(dim_[dim]);
|
||||||
for (auto &e: v)
|
for (auto &e: v)
|
||||||
@ -441,31 +398,34 @@ namespace Grid {
|
|||||||
}
|
}
|
||||||
|
|
||||||
template <typename V>
|
template <typename V>
|
||||||
template <typename W> typename std::enable_if<is_flattenable<W>::isGridTensor>::type
|
Reconstruct<V>::Reconstruct(const std::vector<Element> &flatVector,
|
||||||
Reconstruct<V>::checkInnermost(const W &)
|
|
||||||
{
|
|
||||||
using Traits = GridTypeMapper<typename is_flattenable<W>::grid_type>;
|
|
||||||
const int gridRank{Traits::Rank};
|
|
||||||
const int dimRank{static_cast<int>(dim_.size())};
|
|
||||||
assert(dimRank >= gridRank && "Tensor rank too low for Grid tensor");
|
|
||||||
for (int i=0; i<gridRank; ++i) {
|
|
||||||
assert(dim_[dimRank - gridRank + i] == Traits::Dimension(i) && "Tensor dimension doesn't match Grid tensor");
|
|
||||||
}
|
|
||||||
dim_.resize(dimRank - gridRank);
|
|
||||||
}
|
|
||||||
|
|
||||||
template <typename V>
|
|
||||||
Reconstruct<V>::Reconstruct(const std::vector<Scalar> &flatVector,
|
|
||||||
const std::vector<size_t> &dim)
|
const std::vector<size_t> &dim)
|
||||||
: flatVector_(flatVector)
|
: flatVector_(flatVector)
|
||||||
, dim_(dim)
|
, dim_(dim)
|
||||||
{
|
{
|
||||||
checkInnermost(vector_);
|
|
||||||
assert(dim_.size() == is_flattenable<V>::vecRank && "Tensor rank doesn't match nested std::vector rank");
|
|
||||||
resize(vector_, 0);
|
resize(vector_, 0);
|
||||||
fill(vector_);
|
fill(vector_);
|
||||||
}
|
}
|
||||||
|
|
||||||
|
template <typename V>
|
||||||
|
const V & Reconstruct<V>::getVector(void)
|
||||||
|
{
|
||||||
|
return vector_;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename V>
|
||||||
|
const std::vector<typename Reconstruct<V>::Element> &
|
||||||
|
Reconstruct<V>::getFlatVector(void)
|
||||||
|
{
|
||||||
|
return flatVector_;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename V>
|
||||||
|
const std::vector<size_t> & Reconstruct<V>::getDim(void)
|
||||||
|
{
|
||||||
|
return dim_;
|
||||||
|
}
|
||||||
|
|
||||||
// Vector IO utilities ///////////////////////////////////////////////////////
|
// Vector IO utilities ///////////////////////////////////////////////////////
|
||||||
// helper function to read space-separated values
|
// helper function to read space-separated values
|
||||||
template <typename T>
|
template <typename T>
|
||||||
@ -499,64 +459,6 @@ namespace Grid {
|
|||||||
|
|
||||||
return os;
|
return os;
|
||||||
}
|
}
|
||||||
|
|
||||||
// In general, scalar types are considered "flattenable" (regularly shaped)
|
|
||||||
template <typename T>
|
|
||||||
bool isRegularShapeHelper(const std::vector<T> &, std::vector<std::size_t> &, int, bool)
|
|
||||||
{
|
|
||||||
return true;
|
|
||||||
}
|
|
||||||
|
|
||||||
template <typename T>
|
|
||||||
bool isRegularShapeHelper(const std::vector<std::vector<T>> &v, std::vector<std::size_t> &Dims, int Depth, bool bFirst)
|
|
||||||
{
|
|
||||||
if( bFirst)
|
|
||||||
{
|
|
||||||
assert( Dims.size() == Depth && "Bug: Delete this message after testing" );
|
|
||||||
Dims.push_back(v[0].size());
|
|
||||||
if (!Dims[Depth])
|
|
||||||
return false;
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
assert( Dims.size() >= Depth + 1 && "Bug: Delete this message after testing" );
|
|
||||||
}
|
|
||||||
for (std::size_t i = 0; i < v.size(); ++i)
|
|
||||||
{
|
|
||||||
if (v[i].size() != Dims[Depth] || !isRegularShapeHelper(v[i], Dims, Depth + 1, bFirst && i==0))
|
|
||||||
{
|
|
||||||
return false;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
return true;
|
|
||||||
}
|
|
||||||
|
|
||||||
template <typename T>
|
|
||||||
bool isRegularShape(const T &t) { return true; }
|
|
||||||
|
|
||||||
template <typename T>
|
|
||||||
bool isRegularShape(const std::vector<T> &v) { return !v.empty(); }
|
|
||||||
|
|
||||||
// Return non-zero if all dimensions of this std::vector<std::vector<T>> are regularly shaped
|
|
||||||
template <typename T>
|
|
||||||
bool isRegularShape(const std::vector<std::vector<T>> &v)
|
|
||||||
{
|
|
||||||
if (v.empty() || v[0].empty())
|
|
||||||
return false;
|
|
||||||
// Make sure all of my rows are the same size
|
|
||||||
std::vector<std::size_t> Dims;
|
|
||||||
Dims.reserve(is_flattenable<T>::vecRank);
|
|
||||||
Dims.push_back(v.size());
|
|
||||||
Dims.push_back(v[0].size());
|
|
||||||
for (std::size_t i = 0; i < Dims[0]; ++i)
|
|
||||||
{
|
|
||||||
if (v[i].size() != Dims[1] || !isRegularShapeHelper(v[i], Dims, 2, i==0))
|
|
||||||
{
|
|
||||||
return false;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
return true;
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
|
|
||||||
// helper function to read space-separated values
|
// helper function to read space-separated values
|
||||||
|
@ -417,7 +417,7 @@ public:
|
|||||||
stream << "{";
|
stream << "{";
|
||||||
for (int j = 0; j < N; j++) {
|
for (int j = 0; j < N; j++) {
|
||||||
stream << o._internal[i][j];
|
stream << o._internal[i][j];
|
||||||
if (j < N - 1) stream << ",";
|
if (i < N - 1) stream << ",";
|
||||||
}
|
}
|
||||||
stream << "}";
|
stream << "}";
|
||||||
if (i != N - 1) stream << "\n\t\t";
|
if (i != N - 1) stream << "\n\t\t";
|
||||||
|
Binary file not shown.
@ -1787,7 +1787,7 @@ Hdf5Writer Hdf5Reader HDF5
|
|||||||
|
|
||||||
Write interfaces, similar to the XML facilities in QDP++ are presented. However,
|
Write interfaces, similar to the XML facilities in QDP++ are presented. However,
|
||||||
the serialisation routines are automatically generated by the macro, and a virtual
|
the serialisation routines are automatically generated by the macro, and a virtual
|
||||||
reader and writer interface enables writing to any of a number of formats.
|
reader adn writer interface enables writing to any of a number of formats.
|
||||||
|
|
||||||
**Example**::
|
**Example**::
|
||||||
|
|
||||||
@ -1814,91 +1814,6 @@ reader and writer interface enables writing to any of a number of formats.
|
|||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
Eigen tensor support -- added 2019H1
|
|
||||||
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
|
|
||||||
|
|
||||||
The Serialisation library was expanded in 2019 to support de/serialisation of
|
|
||||||
Eigen tensors. De/serialisation of existing types was not changed. Data files
|
|
||||||
without Eigen tensors remain compatible with earlier versions of Grid and other readers.
|
|
||||||
Conversely, data files containing serialised Eigen tensors is a breaking change.
|
|
||||||
|
|
||||||
Eigen tensor serialisation support was added to BaseIO, which was modified to provide a Traits class
|
|
||||||
to recognise Eigen tensors with elements that are either: primitive scalars (arithmetic and complex types);
|
|
||||||
or Grid tensors.
|
|
||||||
|
|
||||||
**Traits determining de/serialisable scalars**::
|
|
||||||
|
|
||||||
// Is this an Eigen tensor
|
|
||||||
template<typename T> struct is_tensor : std::integral_constant<bool,
|
|
||||||
std::is_base_of<Eigen::TensorBase<T, Eigen::ReadOnlyAccessors>, T>::value> {};
|
|
||||||
// Is this an Eigen tensor of a supported scalar
|
|
||||||
template<typename T, typename V = void> struct is_tensor_of_scalar : public std::false_type {};
|
|
||||||
template<typename T> struct is_tensor_of_scalar<T, typename std::enable_if<is_tensor<T>::value && is_scalar<typename T::Scalar>::value>::type> : public std::true_type {};
|
|
||||||
// Is this an Eigen tensor of a supported container
|
|
||||||
template<typename T, typename V = void> struct is_tensor_of_container : public std::false_type {};
|
|
||||||
template<typename T> struct is_tensor_of_container<T, typename std::enable_if<is_tensor<T>::value && isGridTensor<typename T::Scalar>::value>::type> : public std::true_type {};
|
|
||||||
|
|
||||||
|
|
||||||
Eigen tensors are regular, multidimensional objects, and each Reader/Writer
|
|
||||||
was extended to support this new datatype. Where the Eigen tensor contains
|
|
||||||
a Grid tensor, the dimensions of the data written are the dimensions of the
|
|
||||||
Eigen tensor plus the dimensions of the underlying Grid scalar. Dimensions
|
|
||||||
of size 1 are preserved.
|
|
||||||
|
|
||||||
**New Reader/Writer methods for multi-dimensional data**::
|
|
||||||
|
|
||||||
template <typename U>
|
|
||||||
void readMultiDim(const std::string &s, std::vector<U> &buf, std::vector<size_t> &dim);
|
|
||||||
template <typename U>
|
|
||||||
void writeMultiDim(const std::string &s, const std::vector<size_t> & Dimensions, const U * pDataRowMajor, size_t NumElements);
|
|
||||||
|
|
||||||
|
|
||||||
On readback, the Eigen tensor rank must match the data being read, but the tensor
|
|
||||||
dimensions will be resized if necessary. Resizing is not possible for Eigen::TensorMap<T>
|
|
||||||
because these tensors use a buffer provided at construction, and this buffer cannot be changed.
|
|
||||||
Deserialisation failures cause Grid to assert.
|
|
||||||
|
|
||||||
|
|
||||||
HDF5 Optimisations -- added June 2021
|
|
||||||
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
|
|
||||||
|
|
||||||
Grid serialisation is intended to be light, deterministic and provide a layer of abstraction over
|
|
||||||
multiple file formats. HDF5 excels at handling multi-dimensional data, and the Grid HDF5Reader/HDF5Writer exploits this.
|
|
||||||
When serialising nested ``std::vector<T>``, where ``T`` is an arithmetic or complex type,
|
|
||||||
the Hdf5Writer writes the data as an Hdf5 DataSet object.
|
|
||||||
|
|
||||||
However, nested ``std::vector<std::vector<...T>>`` might be "ragged", i.e. not necessarily regular. E.g. a 3d nested
|
|
||||||
``std::vector`` might contain 2 rows, the first being a 2x2 block and the second row being a 1 x 2 block.
|
|
||||||
A bug existed whereby this was not checked on write, so nested, ragged vectors
|
|
||||||
were written as a regular dataset, with a buffer under/overrun and jumbled contents.
|
|
||||||
|
|
||||||
Clearly this was not used in production, as the bug went undetected until now. Fixing this bug
|
|
||||||
is an opportunity to further optimise the HDF5 file format.
|
|
||||||
|
|
||||||
The goals of this change are to:
|
|
||||||
|
|
||||||
* Make changes to the Hdf5 file format only -- i.e. do not impact other file formats
|
|
||||||
|
|
||||||
* Implement file format changes in such a way that they are transparent to the Grid reader
|
|
||||||
|
|
||||||
* Correct the bug for ragged vectors of numeric / complex types
|
|
||||||
|
|
||||||
* Extend the support of nested std::vector<T> to arbitrarily nested Grid tensors
|
|
||||||
|
|
||||||
|
|
||||||
The trait class ``element`` has been redefined to ``is_flattenable``, which is a trait class for
|
|
||||||
potentially "flattenable" objects. These are (possibly nested) ``std::vector<T>`` where ``T`` is
|
|
||||||
an arithmetic, complex or Grid tensor type. Flattenable objects are tested on write
|
|
||||||
(with the function ``isRegularShape``) to see whether they actually are regular.
|
|
||||||
|
|
||||||
Flattenable, regular objects are written to a multidimensional HDF5 DataSet.
|
|
||||||
Otherwise, an Hdf5 sub group is created with the object "name", and each element of the outer dimension is
|
|
||||||
recursively written to as object "name_n", where n is a 0-indexed number.
|
|
||||||
|
|
||||||
On readback (by Grid)), the presence of a subgroup containing the attribute ``Grid_vector_size`` triggers a
|
|
||||||
"ragged read", otherwise a read from a DataSet is attempted.
|
|
||||||
|
|
||||||
|
|
||||||
Data parallel field IO
|
Data parallel field IO
|
||||||
-----------------------
|
-----------------------
|
||||||
|
|
||||||
|
@ -48,9 +48,7 @@ public:
|
|||||||
std::vector<double>, array,
|
std::vector<double>, array,
|
||||||
std::vector<std::vector<double> >, twodimarray,
|
std::vector<std::vector<double> >, twodimarray,
|
||||||
std::vector<std::vector<std::vector<std::complex<double>> > >, cmplx3darray,
|
std::vector<std::vector<std::vector<std::complex<double>> > >, cmplx3darray,
|
||||||
SpinColourMatrix, scm,
|
SpinColourMatrix, scm
|
||||||
std::vector<std::vector<std::vector<int> > >, ragged,
|
|
||||||
std::vector<std::vector<SpinColourMatrix> >, vscm
|
|
||||||
);
|
);
|
||||||
myclass() {}
|
myclass() {}
|
||||||
myclass(int i)
|
myclass(int i)
|
||||||
@ -58,10 +56,6 @@ public:
|
|||||||
, twodimarray(3,std::vector<double>(5, 1.23456))
|
, twodimarray(3,std::vector<double>(5, 1.23456))
|
||||||
, cmplx3darray(3,std::vector<std::vector<std::complex<double>>>(5, std::vector<std::complex<double>>(7, std::complex<double>(1.2, 3.4))))
|
, cmplx3darray(3,std::vector<std::vector<std::complex<double>>>(5, std::vector<std::complex<double>>(7, std::complex<double>(1.2, 3.4))))
|
||||||
, ve(2, myenum::blue)
|
, ve(2, myenum::blue)
|
||||||
, ragged( {{{i+1},{i+2,i+3}}, // ragged
|
|
||||||
{{i+4,i+5,i+6,i+7},{i+8,i+9,i+10,i+11},{i+12,i+13,i+14,i+15}}, // block
|
|
||||||
{{i+16,i+17},{i+18,i+19,i+20}}} ) //ragged
|
|
||||||
, vscm(3, std::vector<SpinColourMatrix>(5))
|
|
||||||
{
|
{
|
||||||
e=myenum::red;
|
e=myenum::red;
|
||||||
x=i;
|
x=i;
|
||||||
@ -74,13 +68,6 @@ public:
|
|||||||
scm()(0, 2)(1, 1) = 6.336;
|
scm()(0, 2)(1, 1) = 6.336;
|
||||||
scm()(2, 1)(2, 2) = 7.344;
|
scm()(2, 1)(2, 2) = 7.344;
|
||||||
scm()(1, 1)(2, 0) = 8.3534;
|
scm()(1, 1)(2, 0) = 8.3534;
|
||||||
int Counter = i;
|
|
||||||
for( auto & v : vscm ) {
|
|
||||||
for( auto & j : v ) {
|
|
||||||
j = std::complex<double>(Counter, -Counter);
|
|
||||||
Counter++;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
};
|
};
|
||||||
|
|
||||||
|
@ -33,6 +33,7 @@ namespace Grid{
|
|||||||
GRID_SERIALIZABLE_CLASS_MEMBERS(WFParameters,
|
GRID_SERIALIZABLE_CLASS_MEMBERS(WFParameters,
|
||||||
int, steps,
|
int, steps,
|
||||||
double, step_size,
|
double, step_size,
|
||||||
|
double, tol,
|
||||||
int, meas_interval,
|
int, meas_interval,
|
||||||
double, maxTau); // for the adaptive algorithm
|
double, maxTau); // for the adaptive algorithm
|
||||||
|
|
||||||
@ -82,13 +83,27 @@ int main(int argc, char **argv) {
|
|||||||
SU<Nc>::HotConfiguration(pRNG, Umu);
|
SU<Nc>::HotConfiguration(pRNG, Umu);
|
||||||
|
|
||||||
typedef Grid::XmlReader Serialiser;
|
typedef Grid::XmlReader Serialiser;
|
||||||
Serialiser Reader("input.xml");
|
// Serialiser Reader("input.xml");
|
||||||
WFParameters WFPar(Reader);
|
// WFParameters WFPar(Reader);
|
||||||
ConfParameters CPar(Reader);
|
// ConfParameters CPar(Reader);
|
||||||
CheckpointerParameters CPPar(CPar.conf_prefix, CPar.rng_prefix);
|
// WFParameters WFPar;
|
||||||
|
int steps = 800;
|
||||||
|
double step_size=0.02;
|
||||||
|
double tol=1e-4;
|
||||||
|
int meas_interval=50;
|
||||||
|
double maxTau = 16;
|
||||||
|
// ConfParameters CPar;
|
||||||
|
// CPar. conf_prefix="configurations/ckpoint_lat";
|
||||||
|
// CPar. rng_prefix="rngs/ckpoint_rng";
|
||||||
|
// CPar. StartConfiguration=100,
|
||||||
|
// CPar. EndConfiguration=110,
|
||||||
|
// CPar. Skip=1;
|
||||||
|
// CheckpointerParameters CPPar(CPar.conf_prefix, CPar.rng_prefix);
|
||||||
|
CheckpointerParameters CPPar("configurations/ckpoint_lat","rngs/ckpoint_rng");
|
||||||
BinaryHmcCheckpointer<PeriodicGimplR> CPBin(CPPar);
|
BinaryHmcCheckpointer<PeriodicGimplR> CPBin(CPPar);
|
||||||
|
|
||||||
for (int conf = CPar.StartConfiguration; conf <= CPar.EndConfiguration; conf+= CPar.Skip){
|
// for (int conf = CPar.StartConfiguration; conf <= CPar.EndConfiguration; conf+= CPar.Skip){
|
||||||
|
for (int conf = 100; conf <= 110; conf+= 1){
|
||||||
|
|
||||||
CPBin.CheckpointRestore(conf, Umu, sRNG, pRNG);
|
CPBin.CheckpointRestore(conf, Umu, sRNG, pRNG);
|
||||||
|
|
||||||
@ -96,9 +111,10 @@ int main(int argc, char **argv) {
|
|||||||
std::cout << GridLogMessage << "Initial plaquette: "
|
std::cout << GridLogMessage << "Initial plaquette: "
|
||||||
<< WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu) << std::endl;
|
<< WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu) << std::endl;
|
||||||
|
|
||||||
WilsonFlow<PeriodicGimplR> WF(WFPar.steps, WFPar.step_size, WFPar.meas_interval);
|
WilsonFlow<PeriodicGimplR> WF(steps, step_size, meas_interval);
|
||||||
|
|
||||||
WF.smear_adaptive(Uflow, Umu, WFPar.maxTau);
|
// WF.smear_adaptive(Uflow, Umu, maxTau);
|
||||||
|
WF.smear(Uflow, Umu);
|
||||||
|
|
||||||
RealD WFlow_plaq = WilsonLoops<PeriodicGimplR>::avgPlaquette(Uflow);
|
RealD WFlow_plaq = WilsonLoops<PeriodicGimplR>::avgPlaquette(Uflow);
|
||||||
RealD WFlow_TC = WilsonLoops<PeriodicGimplR>::TopologicalCharge(Uflow);
|
RealD WFlow_TC = WilsonLoops<PeriodicGimplR>::TopologicalCharge(Uflow);
|
||||||
|
Reference in New Issue
Block a user