Revision 353ee6cd90b8e10ec132ce3d430319a7b29795f4 authored by Jonas Rembser on 17 April 2024, 17:24:06 UTC, committed by Jonas Rembser on 19 April 2024, 11:35:25 UTC
After commit a27e60a6d4f, it is not important anymore that only the variables used by the expression are passed to RooFormula. Removing the corresponding warnings helps to get rid of useless warnings in the case where you want to try out variations of the formula that omit certain terms, and in particular it helps in `RooAbsData::reduce()`, where the formula is always passed all the varaiables in the dataset, whether the reduction uses them or not.
1 parent 311b78e
stressMathCore.cxx
// @(#)root/test:$Id$
// Author: Lorenzo Moneta 06/2005
///////////////////////////////////////////////////////////////////////////////////
//
// MathCore Benchmark test suite
// ==============================
//
// This program performs tests :
// - mathematical functions in particular the statistical functions by estimating
// pdf, cdf and quantiles. cdf are estimated directly and compared with calculated integral from pdf
// - physics vectors (2D, 3D and 4D) including I/O for every type and for both double and Double32_t
// - SMatrix and SVectors including I/O for double and Double32_t types
// - I/O of complex objects which dictionary has been generated using CINT (default) or Reflex
// TrackD and TrackD32 which contain physics vectors of double and Double32_t
// TrackErrD and TrackErrD32 which contain physics vectors and an SMatrix of double and Double32_t
// VecTrackD which contains an std::vector<TrackD>
//
//
// the program cun run only in compiled mode.
// To run outside ROOT do:
//
// > cd $ROOTSYS/test
// > make stressMathMore
// > ./stressMathMore
//
// to run using REflex set before compiling the environment variable useReflex.
//
// > export useReflex=1
// > make stressMathMore
// > ./stressMathMore
//
// to run inside ROOT using ACliC
// for using CINT you need first to have the library libTrackMathCoreDict.so
// (type: make libTrackMathCoreDict.so to make it)
//
// root> gSystem->Load("libMathCore");
// root> gSystem->Load("libTree");
// root> gSystem->Load("libHist");
// root> .x stressMathCore.cxx+
//
// for using Reflex dictionaries you need first to have the library libTrackMathCoreRflx.so
// (type: make libTrackMathCoreRflx.so to make it)
//
// root> gSystem->Load("libMathCore");
// root> gSystem->Load("libTree");
// root> gSystem->Load("libHist");
// root> gSystem->Load("libReflex");
// root> gSystem->SetIncludePath("-DUSE_REFLEX");
// root> .x stressMathCore.cxx+
//
//
#ifndef __CINT__
#include "Math/DistFuncMathCore.h"
//#define USE_MATHMORE
#ifdef USE_MATHMORE
#include "Math/DistFuncMathMore.h"
#endif
#include "Math/IParamFunction.h"
#include "Math/Integrator.h"
#include "Math/IntegratorOptions.h"
#include <iostream>
#include <iomanip>
#include <limits>
#include <cmath>
#include "TBenchmark.h"
#include "TROOT.h"
#include "TRandom3.h"
#include "TSystem.h"
#include "TTree.h"
#include "TFile.h"
#include "TF1.h"
#include "TMath.h"
#include "Math/Vector2D.h"
#include "Math/Vector3D.h"
#include "Math/Vector4D.h"
#include "Math/VectorUtil.h"
#include "Math/SVector.h"
#include "Math/SMatrix.h"
R__ADD_INCLUDE_PATH($ROOTSYS/test)
#include "TrackMathCore.h"
#include "Math/GenVector/RotationZ.h" // Workaround to autoload libGenVector ROOT-7056
using namespace ROOT::Math;
#endif
//#define DEBUG
bool debug = true; // print out reason of test failures
bool debugTime = false; // print out separate timings for vectors
bool removeFiles = true; // remove Output root files
void PrintTest(std::string name) {
std::cout << std::left << std::setw(40) << name;
}
void PrintStatus(int iret) {
if (iret == 0)
std::cout <<"\t\t................ OK" << std::endl;
else
std::cout <<"\t\t............ FAILED " << std::endl;
}
int compare( std::string name, double v1, double v2, double scale = 2.0) {
// ntest = ntest + 1;
//std::cout << std::setw(50) << std::left << name << ":\t";
// numerical double limit for epsilon
double eps = scale* std::numeric_limits<double>::epsilon();
int iret = 0;
double delta = v2 - v1;
double d = 0;
if (delta < 0 ) delta = - delta;
if (v1 == 0 || v2 == 0) {
if (delta > eps ) {
iret = 1;
}
}
// skip case v1 or v2 is infinity
else {
d = v1;
if ( v1 < 0) d = -d;
// add also case when delta is small by default
if ( delta/d > eps && delta > eps )
iret = 1;
}
if (iret) {
if (debug) {
int pr = std::cout.precision (18);
std::cout << "\nDiscrepancy in " << name.c_str() << "() :\n " << v1 << " != " << v2 << " discr = " << int(delta/d/eps)
<< " (Allowed discrepancy is " << eps << ")\n\n";
std::cout.precision (pr);
//nfail = nfail + 1;
}
}
//else
// std::cout <<".";
return iret;
}
#ifndef __CINT__
// trait class for distinguishing the number of parameters for the various functions
template<class Func, unsigned int NPAR>
struct Evaluator {
static double F(Func f, double x, const double * ) {
return f(x);
}
};
template<class Func>
struct Evaluator<Func, 1> {
static double F(Func f, double x, const double * p) {
return f(x,p[0]);
}
};
template<class Func>
struct Evaluator<Func, 2> {
static double F(Func f, double x, const double * p) {
return f(x,p[0],p[1]);
}
};
template<class Func>
struct Evaluator<Func, 3> {
static double F(Func f, double x, const double * p) {
return f(x,p[0],p[1],p[2]);
}
};
// global test variable
int NFuncTest = 100;
// statistical function class
// template on the number of parameters
template<class Func, class FuncQ, int NPAR, int NPARQ=NPAR-1>
class StatFunction : public ROOT::Math::IParamFunction {
public:
StatFunction(Func pdf, Func cdf, FuncQ quant) : fPdf(pdf), fCdf(cdf), fQuant(quant)
{
fScale1 = 1.0E6; //scale for cdf test (integral)
fScale2 = 10; //scale for quantile test
for(int i = 0; i< NPAR; ++i) fParams[i]=0;
}
unsigned int NPar() const override { return NPAR; }
const double * Parameters() const override { return fParams; }
ROOT::Math::IGenFunction * Clone() const override { return new StatFunction(fPdf,fCdf,fQuant); }
void SetParameters(const double * p) override { std::copy(p,p+NPAR,fParams); }
void SetParameters(double p0) { *fParams = p0; }
void SetParameters(double p0, double p1) { *fParams = p0; *(fParams+1) = p1; }
void SetParameters(double p0, double p1, double p2) { *fParams = p0; *(fParams+1) = p1; *(fParams+2) = p2; }
static void SetNTest(int n) { NFuncTest = n; }
double Cdf(double x) const {
return Evaluator<Func,NPAR>::F(fCdf,x, fParams);
}
double Quantile(double x) const {
double z = Evaluator<FuncQ,NPARQ>::F(fQuant,x, fParams);
if ((NPAR - NPARQ) == 1)
z += fParams[NPAR-1]; // adjust the offset
return z;
}
// test cumulative function
int Test(double x1, double x2, double xl = 1, double xu = 0, bool cumul = false);
void ScaleTol1(double s) { fScale1 *= s; }
void ScaleTol2(double s) { fScale2 *= s; }
private:
double DoEvalPar(double x, const double * ) const override {
// implement explicitly using cached parameter values
return Evaluator<Func,NPAR>::F(fPdf,x, fParams);
}
Func fPdf;
Func fCdf;
FuncQ fQuant;
double fParams[NPAR];
double fScale1;
double fScale2;
};
// test cdf at value f
template<class F1, class F2, int N1, int N2>
int StatFunction<F1,F2,N1,N2>::Test(double xmin, double xmax, double xlow, double xup, bool c) {
int iret = 0;
// scan all values from xmin to xmax
double dx = (xmax-xmin)/NFuncTest;
#ifdef USE_MATHMORE
for (int i = 0; i < NFuncTest; ++i) {
double v1 = xmin + dx*i; // value used for testing
double q1 = Cdf(v1);
//std::cout << "v1 " << v1 << " pdf " << (*this)(v1) << " cdf " << q1 << " quantile " << Quantile(q1) << std::endl;
// calculate integral of pdf
Integrator ig(IntegrationOneDim::kADAPTIVESINGULAR, 1.E-12,1.E-12,100000);
ig.SetFunction(*this);
double q2 = 0;
if (!c) {
// lower intergal (cdf)
if (xlow >= xup || xlow > xmin)
q2 = ig.IntegralLow(v1);
else
q2 = ig.Integral(xlow,v1);
// use a larger scale (integral error is 10-9)
iret |= compare("test _cdf", q1, q2, fScale1 );
// test the quantile
double v2 = Quantile(q1);
iret |= compare("test _quantile", v1, v2, fScale2 );
}
else {
// upper integral (cdf_c)
if (xlow >= xup || xup < xmax)
q2 = ig.IntegralUp(v1);
else
q2 = ig.Integral(v1,xup);
iret |= compare("test _cdf_c", q1, q2, fScale1);
double v2 = Quantile(q1);
iret |= compare("test _quantile_c", v1, v2, fScale2 );
}
if (iret) {
std::cout << "Failed test for x = " << v1 << " p = ";
for (int j = 0; j < N1; ++j) std::cout << fParams[j] << "\t";
std::cout << std::endl;
break;
}
}
#else
// use TF1 for the integral
// std::cout << "xlow-xuo " << xlow << " " << xup << std::endl;
// std::cout << "xmin-xmax " << xmin << " " << xmax << std::endl;
double x1,x2 = 0;
if (xlow >= xup) {
x1 = -100; x2 = 100;
}
else if (xup < xmax) {
x1 = xlow; x2 = 100;
}
else {
x1=xlow; x2 = xup;
}
//std::cout << "x1-x2 " << x1 << " " << x2 << std::endl;
TF1 * f = new TF1("ftemp",ParamFunctor(*this),x1,x2,0);
//ROOT::Math::IntegratorOneDimOptions::SetDefaultIntegrator("Gauss");
for (int i = 0; i < NFuncTest; ++i) {
double v1 = xmin + dx*i; // value used for testing
double q1 = Cdf(v1);
double q2 = 0;
if (!c) {
q2 = f->Integral(x1,v1,1.E-12);
// use a larger scale (integral error is 10-9)
iret |= compare("test _cdf", q1, q2, fScale1 );
// test the quantile
double v2 = Quantile(q1);
iret |= compare("test _quantile", v1, v2, fScale2 );
}
else {
// upper integral (cdf_c)
q2 = f->Integral(v1,x2,1.E-12);
iret |= compare("test _cdf_c", q1, q2, fScale1);
double v2 = Quantile(q1);
iret |= compare("test _quantile_c", v1, v2, fScale2 );
}
// if (!c) std::cout << "i = " << i << " x1 = " << x1 << " v = " << v1 << " pdf " << (*this)(v1) << " cdf " << q1 << " Int(x1,v)= " << q2 << std::endl;
// else
// std::cout << "i = " << i << " x2 = " << x2 << " v = " << v1 << " pdf " << (*this)(v1) << " cdf_c " << q1 << " Int(v,x2)= " << q2 << std::endl;
}
delete f;
#endif
if (c || iret != 0) PrintStatus(iret);
return iret;
}
// typedef defining the functions
typedef double ( * F0) ( double);
typedef double ( * F1) ( double, double);
typedef double ( * F2) ( double, double, double);
typedef double ( * F3) ( double, double, double, double);
typedef StatFunction<F2,F2,2,2> Dist_beta;
typedef StatFunction<F2,F1,2> Dist_breitwigner;
typedef StatFunction<F2,F1,2> Dist_chisquared;
typedef StatFunction<F3,F2,3> Dist_fdistribution;
typedef StatFunction<F3,F2,3> Dist_gamma;
typedef StatFunction<F2,F1,2> Dist_gaussian;
typedef StatFunction<F3,F2,3> Dist_lognormal;
typedef StatFunction<F2,F1,2> Dist_tdistribution;
typedef StatFunction<F2,F1,2> Dist_exponential;
typedef StatFunction<F2,F1,2> Dist_landau;
typedef StatFunction<F3,F2,3> Dist_uniform;
#define CREATE_DIST(name) Dist_ ##name dist( name ## _pdf, name ## _cdf, name ##_quantile );
#define CREATE_DIST_C(name) Dist_ ##name distc( name ## _pdf, name ## _cdf_c, name ##_quantile_c );
// #define CREATE_DIST_C(name) Dist_ ##name distc( name ## _pdf, name ## _cdf_c, 0 );
// #else
// #endif
template<class Distribution>
int TestDist(Distribution & d, double x1, double x2) {
int ir = 0;
ir |= d.Test(x1,x2);
return ir;
}
int testStatFunctions(int nfunc = 100 ) {
// test statistical functions
int iret = 0;
NFuncTest = nfunc;
{
PrintTest("Beta distribution");
CREATE_DIST(beta);
dist.SetParameters( 2, 2);
iret |= dist.Test(0.01,0.99,0.,1.);
CREATE_DIST_C(beta);
distc.SetParameters( 2, 2);
iret |= distc.Test(0.01,0.99,0.,1.,true);
}
{
PrintTest("Gamma distribution");
#ifdef USE_MATHMORE // gamma_quantile is in mathmore
CREATE_DIST(gamma);
dist.SetParameters( 2, 1);
iret |= dist.Test(0.05,5, 0.,TMath::Infinity());
#endif
CREATE_DIST_C(gamma);
distc.SetParameters( 2, 1);
iret |= distc.Test(0.05,5, 0.,TMath::Infinity(),true);
}
{
PrintTest("Chisquare distribution");
#ifdef USE_MATHMORE
CREATE_DIST(chisquared);
dist.SetParameters( 10, 0);
dist.ScaleTol2(10);
iret |= dist.Test(0.05,30, 0.,1.);
#endif
CREATE_DIST_C(chisquared);
distc.SetParameters( 10, 0);
distc.ScaleTol2(10000000); // t.b.c.
iret |= distc.Test(0.05,30, 0.,TMath::Infinity(),true);
}
{
PrintTest("Normal distribution ");
CREATE_DIST(gaussian);
dist.SetParameters( 2, 1);
dist.ScaleTol2(100);
iret |= dist.Test(-3,5);
CREATE_DIST_C(gaussian);
distc.SetParameters( 1, 0);
distc.ScaleTol2(100);
iret |= distc.Test(-3,5,1,0,true);
}
{
PrintTest("BreitWigner distribution ");
CREATE_DIST(breitwigner);
dist.SetParameters( 1);
dist.ScaleTol1(1E8);
dist.ScaleTol2(10);
iret |= dist.Test(-5,5);
CREATE_DIST_C(breitwigner);
distc.SetParameters( 1);
distc.ScaleTol1(1E8);
distc.ScaleTol2(10);
iret |= distc.Test(-5,5,1,0,true);
}
{
PrintTest("F distribution ");
CREATE_DIST(fdistribution);
dist.SetParameters( 5, 4);
dist.ScaleTol1(1000000);
dist.ScaleTol2(10);
// if enlarge scale test fails
iret |= dist.Test(0.05,5,0,1);
CREATE_DIST_C(fdistribution);
distc.SetParameters( 5, 4);
distc.ScaleTol1(100000000);
distc.ScaleTol2(10);
// if enlarge scale test fails
iret |= distc.Test(0.05,5,0,TMath::Infinity(),true);
}
#ifdef USE_MATHMORE // wait t quantile is in mathcore
{
PrintTest("t distribution ");
CREATE_DIST(tdistribution);
dist.SetParameters( 10 );
// dist.ScaleTol1(1000);
dist.ScaleTol2(5000);
iret |= dist.Test(-10,10);
CREATE_DIST_C(tdistribution);
distc.SetParameters( 10 );
distc.ScaleTol2(10000); // t.b.c.
iret |= distc.Test(-10,10,1,0,true);
}
#endif
{
PrintTest("lognormal distribution");
CREATE_DIST(lognormal);
dist.SetParameters(1,1 );
dist.ScaleTol1(1000);
iret |= dist.Test(0.01,5,0,1);
CREATE_DIST_C(lognormal);
distc.SetParameters(1,1 );
distc.ScaleTol1(1000);
distc.ScaleTol2(1000000); // t.b.c.
iret |= distc.Test(0.01,5,0,TMath::Infinity(),true);
}
{
PrintTest("Exponential distribution");
CREATE_DIST(exponential);
dist.SetParameters( 2);
dist.ScaleTol2(100);
iret |= dist.Test(0.,5.,0.,1.);
CREATE_DIST_C(exponential);
distc.SetParameters( 2);
distc.ScaleTol2(100);
iret |= distc.Test(0.,5.,0.,1.,true);
}
{
PrintTest("Landau distribution");
CREATE_DIST(landau);
dist.SetParameters( 2);
// Landau is not very precise (put prec at 10-6)
// as indicated in Landau paper (
dist.ScaleTol1(10000);
dist.ScaleTol2(1.E10);
iret |= dist.Test(-1,10,-TMath::Infinity(),TMath::Infinity());
CREATE_DIST_C(landau);
distc.SetParameters( 2);
distc.ScaleTol1(1000);
distc.ScaleTol2(1E10);
iret |= distc.Test(-1,10,-TMath::Infinity(), TMath::Infinity(),true);
}
{
PrintTest("Uniform distribution");
CREATE_DIST(uniform);
dist.SetParameters( 1, 2);
iret |= dist.Test(1.,2.,1.,2.);
CREATE_DIST_C(uniform);
distc.SetParameters( 1, 2);
iret |= distc.Test(1.,2.,1.,2.,true);
}
return iret;
}
//*******************************************************************************************************************
// GenVector tests
//*******************************************************************************************************************
// trait for getting vector name
template<class V>
struct VecType {
static std::string name() { return "MathCoreVector";}
};
template<>
struct VecType<XYVector> {
static std::string name() { return "XYVector";}
static std::string name32() { return "ROOT::Math::DisplacementVector2D<ROOT::Math::Cartesian2D<Double32_t> >";}
};
template<>
struct VecType<Polar2DVector> {
static std::string name() { return "Polar2DVector";}
};
template<>
struct VecType<XYZVector> {
static std::string name() { return "XYZVector";}
static std::string name32() { return "ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<Double32_t> >";}
};
template<>
struct VecType<Polar3DVector> {
static std::string name() { return "Polar3DVector";}
};
template<>
struct VecType<RhoEtaPhiVector> {
static std::string name() { return "RhoEtaPhiVector";}
};
template<>
struct VecType<RhoZPhiVector> {
static std::string name() { return "RhoZPhiVector";}
};
template<>
struct VecType<PxPyPzEVector> {
static std::string name() { return "PxPyPzEVector";}
static std::string name32() { return "ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<Double32_t> >";}
};
template<>
struct VecType<PtEtaPhiEVector> {
static std::string name() { return "PtEtaPhiEVector";}
};
template<>
struct VecType<PtEtaPhiMVector> {
static std::string name() { return "PtEtaPhiMVector";}
};
template<>
struct VecType<PxPyPzMVector> {
static std::string name() { return "PxPyPzMVector";}
};
template<>
struct VecType<SVector<double,3> > {
static std::string name() { return "SVector3";}
static std::string name32() { return "ROOT::Math::SVector<Double32_t,3>";}
};
template<>
struct VecType<SVector<double,4> > {
static std::string name() { return "SVector4";}
static std::string name32() { return "ROOT::Math::SVector<Double32_t,4>";}
};
template<>
struct VecType<TrackD> {
static std::string name() { return "TrackD";}
};
template<>
struct VecType<TrackD32> {
static std::string name() { return "TrackD32";}
};
template<>
struct VecType<TrackErrD> {
static std::string name() { return "TrackErrD";}
};
template<>
struct VecType<TrackErrD32> {
static std::string name() { return "TrackErrD32";}
};
template<>
struct VecType<VecTrack<TrackD> > {
static std::string name() { return "VecTrackD";}
};
template<>
struct VecType<VecTrack<TrackErrD> > {
static std::string name() { return "VecTrackErrD";}
};
// generic (2 dim)
template<class V, int Dim>
struct VecOp {
template<class It>
static V Create(It &x, It &y, It & , It& ) { return V(*x++,*y++); }
template<class It>
static void Set(V & v, It &x, It &y, It &, It&) { v.SetXY(*x++,*y++); }
static double Add(const V & v) { return v.X() + v.Y(); }
static double Delta(const V & v1, const V & v2) { double d = ROOT::Math::VectorUtil::DeltaPhi(v1,v2); return d*d; } // is v2-v1
};
// specialized for 3D
template<class V>
struct VecOp<V,3> {
template<class It>
static V Create(It &x, It& y, It& z , It& ) { return V(*x++,*y++,*z++); }
template<class It>
static void Set(V & v, It & x, It &y, It &z, It&) { v.SetXYZ(*x++,*y++,*z++); }
static V Create(double x, double y, double z , double ) { return V(x,y,z); }
static void Set(V & v, double x, double y, double z, double) { v.SetXYZ(x,y,z); }
static double Add(const V & v) { return v.X() + v.Y() + v.Z(); }
static double Delta(const V & v1, const V & v2) { return ROOT::Math::VectorUtil::DeltaR2(v1,v2); }
};
// specialized for 4D
template<class V>
struct VecOp<V,4> {
template<class It>
static V Create(It &x, It &y, It &z , It &t ) { return V(*x++,*y++,*z++,*t++);}
template<class It>
static void Set(V & v, It & x, It &y, It &z, It &t) { v.SetXYZT(*x++,*y++,*z++,*t++); }
static double Add(const V & v) { return v.X() + v.Y() + v.Z() + v.E(); }
static double Delta(const V & v1, const V & v2) {
return ROOT::Math::VectorUtil::DeltaR2(v1,v2) + ROOT::Math::VectorUtil::InvariantMass(v1,v2); }
};
// specialized for SVector<3>
template<>
struct VecOp<SVector<double,3>,3> {
typedef SVector<double,3> V;
template<class It>
static V Create(It &x, It &y, It &z , It & ) { return V(*x++,*y++,*z++);}
static double Add(const V & v) { return v(0) + v(1) + v(2); }
};
// specialized for SVector<4>
template<>
struct VecOp<SVector<double,4>,4> {
typedef SVector<double,4> V;
template<class It>
static V Create(It &x, It &y, It &z , It &t ) { return V(*x++,*y++,*z++,*t++);}
static double Add(const V & v) { return v(0) + v(1) + v(2) + v(3); }
};
// internal structure to measure the time
TStopwatch gTimer;
double gTotTime;
struct Timer {
Timer() {
gTimer.Start();
}
~Timer() {
gTimer.Stop();
gTotTime += Time();
if (debugTime) printTime();
}
void printTime( std::string s = "") {
int pr = std::cout.precision(8);
std::cout << s << "\t" << " time = " << Time() << "\t(sec)\t"
// << time.CpuTime()
<< std::endl;
std::cout.precision(pr);
}
double Time() { return gTimer.RealTime(); } // use real time
TStopwatch gTimer;
double gTotTime;
};
template<int Dim>
class VectorTest {
private:
// global data variables
std::vector<double> dataX;
std::vector<double> dataY;
std::vector<double> dataZ;
std::vector<double> dataE;
int nGen;
int n2Loop ;
double fSum; // total sum of x,y,z,t (for testing first addition)
public:
VectorTest(int n1, int n2=0) :
nGen(n1),
n2Loop(n2)
{
gTotTime = 0;
}
double TotalTime() const { return gTotTime; } // use real time
double Sum() const { return fSum; }
int check(std::string name, double s1, double s2, double scale=1) {
int iret = 0;
PrintTest(name);
iret |= compare(name,s1,s2,scale);
PrintStatus(iret);
return iret;
}
void print(std::string name) {
PrintTest(name);
std::cout <<"\t\t..............\n";
}
void genData() {
// generate for all 4 d data
TRandom3 r(111); // use a fixed seed to be able to reproduce tests
fSum = 0;
for (int i = 0; i < nGen ; ++ i) {
// generate a 4D vector and stores only the interested dimensions
double phi = r.Rndm()*3.1415926535897931;
double eta = r.Uniform(-5.,5.);
double pt = r.Exp(10.);
double m = r.Uniform(0,10.);
if ( i%50 == 0 )
m = r.BreitWigner(1.,0.01);
double E = sqrt( m*m + pt*pt*std::cosh(eta)*std::cosh(eta) );
// fill vectors
PtEtaPhiEVector q( pt, eta, phi, E);
dataX.push_back( q.x() );
dataY.push_back( q.y() );
fSum += q.x() + q.y();
if (Dim >= 3) {
dataZ.push_back( q.z() );
fSum += q.z();
}
if (Dim >=4 ) {
dataE.push_back( q.t() );
fSum += q.t();
}
}
assert( int(dataX.size()) == nGen);
assert( int(dataY.size()) == nGen);
if (Dim >= 3) assert( int(dataZ.size()) == nGen);
if (Dim >=4 ) assert( int(dataE.size()) == nGen);
// // // dataZ.resize(nGen);
// // // dataE.resize(nGen);
}
// gen data for a Ndim matrix or vector
void genDataN() {
// generate for all 4 d data
TRandom3 r(111); // use a fixed seed to be able to reproduce tests
fSum = 0;
dataX.reserve(nGen*Dim);
for (int i = 0; i < nGen*Dim ; ++ i) {
// generate random data between [0,1]
double x = r.Rndm();
fSum += x;
dataX.push_back( x );
}
}
typedef std::vector<double>::const_iterator DataIt;
// test methods
template <class V>
void testCreate( std::vector<V > & dataV) {
Timer tim;
DataIt x = dataX.begin();
DataIt y = dataY.begin();
DataIt z = dataZ.begin();
DataIt t = dataE.begin();
while (x != dataX.end() ) {
dataV.push_back(VecOp<V,Dim>::Create(x,y,z,t) );
assert(int(dataV.size()) <= nGen);
}
}
template <class V>
void testCreateAndSet( std::vector<V > & dataV) {
Timer tim;
DataIt x = dataX.begin();
DataIt y = dataY.begin();
DataIt z = dataZ.begin();
DataIt t = dataE.begin();
while (x != dataX.end() ) {
V v;
VecOp<V,Dim>::Set( v, x,y,z,t);
dataV.push_back(v);
assert(int(dataV.size()) <= nGen);
}
}
template <class V>
double testAddition( const std::vector<V > & dataV) {
V v0;
Timer t;
for (int i = 0; i < nGen; ++i) {
v0 += dataV[i];
}
return VecOp<V,Dim>::Add(v0);
}
template <class V>
double testOperations( const std::vector<V > & dataV) {
double tot = 0;
Timer t;
for (int i = 0; i < nGen-1; ++i) {
const V & v1 = dataV[i];
const V & v2 = dataV[i+1];
double a = v1.R();
double b = v2.mag2(); // mag2 is defined for all dimensions;
double c = 1./v1.Dot(v2);
V v3 = c * ( v1/a + v2/b );
tot += VecOp<V,Dim>::Add(v3);
}
return tot;
}
// mantain loop in gen otherwise is proportional to N**@
template <class V>
double testDelta( const std::vector<V > & dataV) {
double tot = 0;
Timer t;
for (int i = 0; i < nGen-1; ++i) {
const V & v1 = dataV[i];
const V & v2 = dataV[i+1];
tot += VecOp<V,Dim>::Delta(v1,v2);
}
return tot;
}
// template <class V>
// double testDotProduct( const std::vector<V *> & dataV) {
// //unsigned int n = std::min(n2Loop, dataV.size() );
// double tot = 0;
// V v0 = *(dataV[0]);
// Timer t;
// for (unsigned int i = 0; i < nGen-1; ++i) {
// V & v1 = *(dataV[i]);
// tot += v0.Dot(v1);
// }
// return tot;
// }
template <class V1, class V2>
void testConversion( std::vector<V1 > & dataV1, std::vector<V2 > & dataV2) {
Timer t;
for (int i = 0; i < nGen; ++i) {
dataV2.push_back( V2( dataV1[i] ) );
}
}
// rotation
template <class V, class R>
double testRotation( std::vector<V > & dataV ) {
double sum = 0;
double rotAngle = 1;
Timer t;
for (unsigned int i = 0; i < nGen; ++i) {
V & v1 = dataV[i];
V v2 = v1;
v2.Rotate(rotAngle);
sum += VecOp<V,Dim>::Add(v2);
}
return sum;
}
template<class V>
double testWrite(const std::vector<V> & dataV, std::string typeName="", bool compress = false) {
std::string fname = VecType<V>::name() + ".root";
// replace < character with _
TFile file(fname.c_str(),"RECREATE","",compress);
// create tree
std::string tree_name="Tree with" + VecType<V>::name();
TTree tree("VectorTree",tree_name.c_str());
V *v1 = new V();
//std::cout << "typeID written : " << typeid(*v1).name() << std::endl;
// need to add namespace to full type name
if (typeName == "") {
typeName = "ROOT::Math::" + VecType<V>::name();
}
//std::cout << typeName << std::endl;
TBranch * br = tree.Branch("Vector_branch",typeName.c_str(),&v1);
if (br == 0) {
std::cout << "Error creating branch for" << typeName << "\n\t typeid is "
<< typeid(*v1).name() << std::endl;
return -1;
}
Timer timer;
for (int i = 0; i < nGen; ++i) {
*v1 = dataV[i];
tree.Fill();
}
#ifdef DEBUG
tree.Print(); // debug
#endif
file.Write();
file.Close();
return file.GetSize();
}
template<class V>
int testRead(std::vector<V> & dataV) {
dataV.clear();
dataV.reserve(nGen);
std::string fname = VecType<V>::name() + ".root";
TFile f1(fname.c_str());
if (f1.IsZombie() ) {
std::cout << " Error opening file " << fname << std::endl;
return -1;
}
// create tree
TTree *tree = dynamic_cast<TTree*>(f1.Get("VectorTree"));
if (tree == 0) {
std::cout << " Error reading file " << fname << std::endl;
return -1;
}
V *v1 = 0;
//std::cout << "reading typeID : " << typeid(*v1).name() << std::endl;
// cast to void * to avoid a warning
tree->SetBranchAddress("Vector_branch",(void *) &v1);
Timer timer;
int n = (int) tree->GetEntries();
if (n != nGen) {
std::cout << "wrong tree entries from file" << fname << std::endl;
return -1;
}
for (int i = 0; i < n; ++i) {
tree->GetEntry(i);
dataV.push_back(*v1);
}
if (removeFiles) gSystem->Unlink(fname.c_str());
return 0;
}
// test of SVEctor's or SMatrix
template<class V>
void testCreateSV( std::vector<V> & dataV) {
Timer tim;
//DataIt x = dataX.begin();
double * x = &dataX.front();
double * end = x + dataX.size();
//SVector cannot be created from a generic iterator (should be fixed)
while (x != end) {
dataV.push_back( V(x,x+Dim) );
assert(int(dataV.size()) <= nGen);
x += Dim;
}
}
template <class V>
double testAdditionSV( const std::vector<V > & dataV) {
V v0;
Timer t;
for (int i = 0; i < nGen; ++i) {
v0 += dataV[i];
}
double tot = 0;
typedef typename V::const_iterator It;
for (It itr = v0.begin(); itr != v0.end(); ++itr)
tot += *itr;
return tot;
}
template <class V>
double testAdditionTR( const std::vector<V > & dataV) {
V v0;
Timer t;
for (int i = 0; i < nGen; ++i) {
v0 += dataV[i];
}
#ifdef DEBUG
v0.Print();
#endif
return v0.Sum();
}
};
//--------------------------------------------------------------------------------------
// test of all physics vector (GenVector's)
//--------------------------------------------------------------------------------------
template<class V1, class V2, int Dim>
int testVector(int ngen, bool testio=false) {
int iret = 0;
VectorTest<Dim> a(ngen);
a.genData();
std::vector<V1> v1;
std::vector<V2> v2;
v1.reserve(ngen);
v2.reserve(ngen);
double s1, s2 = 0;
double scale = 1;
double sref1, sref2 = 0;
a.testCreate(v1); iret |= a.check(VecType<V1>::name()+" creation",v1.size(),ngen);
s1 = a.testAddition(v1); iret |= a.check(VecType<V1>::name()+" addition",s1,a.Sum(),Dim*20);
sref1 = s1;
v1.clear();
assert(v1.size() == 0);
a.testCreateAndSet(v1); iret |= a.check(VecType<V1>::name()+" creation",v1.size(),ngen);
s2 = a.testAddition(v1); iret |= a.check(VecType<V1>::name()+" setting",s2,s1);
a.testConversion(v1,v2); iret |= a.check(VecType<V1>::name()+" -> " + VecType<V2>::name(),v1.size(),v2.size() );
scale = 1000;
if (Dim == 2) scale = 1.E12; // to be understood
if (Dim == 3) scale = 1.E4; // problem with RhoEtaPhiVector
s2 = a.testAddition(v2); iret |= a.check("Vector conversion",s2,s1,scale);
sref2 = s2;
s1 = a.testOperations(v1); a.print(VecType<V1>::name()+" operations");
scale = Dim*20;
if (Dim==3 && VecType<V2>::name() == "RhoEtaPhiVector") scale *= 12; // for problem with RhoEtaPhi
if (Dim==4 && ( VecType<V2>::name() == "PtEtaPhiMVector" || VecType<V2>::name() == "PxPyPzMVector")) {
#if (defined(__arm__) || defined(__arm64__) || defined(__aarch64__) || defined(__s390x__))
scale *= 1.E7;
#else
scale *= 10;
#endif
#if defined(__FAST_MATH__) && defined(__clang__)
scale *= 1.E6;
#endif
}
// for problem with PtEtaPhiE
#if defined (R__LINUX) && !defined(R__B64)
// problem of precision on linux 32
if (Dim ==4) scale = 1000000000;
#endif
if (Dim==4 && VecType<V2>::name() == "PtEtaPhiEVector") scale = 0.01/(std::numeric_limits<double>::epsilon());
s2 = a.testOperations(v2); iret |= a.check(VecType<V2>::name()+" operations",s2,s1,scale);
s1 = a.testDelta(v1); a.print(VecType<V1>::name()+" delta values");
scale = Dim*16;
if (Dim==4) scale *= 100; // for problem with PtEtaPhiE
s2 = a.testDelta(v2); iret |= a.check(VecType<V2>::name()+" delta values",s2,s1,scale);
double fsize = 0;
int ir = 0;
if (!testio) return iret;
double estSize = ngen*8 * Dim + 10000; //add extra bytes
scale = 0.1 / std::numeric_limits<double>::epsilon();
fsize = a.testWrite(v1); iret |= a.check(VecType<V1>::name()+" write",fsize,estSize,scale);
ir = a.testRead(v1); iret |= a.check(VecType<V1>::name()+" read",ir,0);
s1 = a.testAddition(v1); iret |= a.check(VecType<V1>::name()+" after read",s1,sref1);
// test io vector 2
fsize = a.testWrite(v2); iret |= a.check(VecType<V2>::name()+" write",fsize,estSize,scale);
ir = a.testRead(v2); iret |= a.check(VecType<V2>::name()+" read",ir,0);
scale = 100; // gcc 4.3.2 gives and error for RhoEtaPhiVector for 32 bits
s2 = a.testAddition(v2); iret |= a.check(VecType<V2>::name()+" after read",s2,sref2,scale);
// test io of double 32 for vector 1
if (Dim==2) return iret;
estSize = ngen*4 * Dim + 10000; //add extra bytes
scale = 0.1 / std::numeric_limits<double>::epsilon();
std::string typeName = VecType<V1>::name32();
fsize = a.testWrite(v1,typeName); iret |= a.check(VecType<V1>::name() +"_D32 write",fsize,estSize,scale);
ir = a.testRead(v1); iret |= a.check(VecType<V1>::name()+"_D32 read",ir,0);
s1 = a.testAddition(v1); iret |= a.check(VecType<V1>::name()+"_D32 after read",s1,sref1,1.E9);
return iret;
}
//--------------------------------------------------------------------------------------
// test of Svector of dim 3 or 4
//--------------------------------------------------------------------------------------
template<int Dim>
int testVector34(int ngen, bool testio=false) {
int iret = 0;
VectorTest<Dim> a(ngen);
a.genData();
typedef SVector<double, Dim> SV;
std::vector<SV> v1;
v1.reserve(ngen);
double s1 = 0;
//double scale = 1;
double sref1 = 0;
std::string name = "SVector<double," + Util::ToString(Dim) + ">";
a.testCreate(v1); iret |= a.check(name+" creation",v1.size(),ngen);
s1 = a.testAddition(v1); iret |= a.check(name+" addition",s1,a.Sum(),Dim*20);
sref1 = s1;
// test the io
double fsize = 0;
int ir = 0;
if (!testio) return iret;
std::string typeName = "ROOT::Math::"+name;
double estSize = ngen*8 * Dim + 47000; // add extra bytes (why so much ? ) this is for ngen=10000
double scale = 0.1 / std::numeric_limits<double>::epsilon();
fsize = a.testWrite(v1,typeName); iret |= a.check(name+" write",fsize,estSize,scale);
ir = a.testRead(v1); iret |= a.check(name+" read",ir,0);
s1 = a.testAddition(v1); iret |= a.check(name+" after read",s1,sref1,10);
//std::cout << "File size = " << fsize << " estimated " << 8 * Dim * ngen << std::endl;
// test Double32
estSize = ngen*4 * Dim + 47000; //add extra bytes
scale = 0.1 / std::numeric_limits<double>::epsilon();
typeName = VecType<SV>::name32();
fsize = a.testWrite(v1,typeName); iret |= a.check(VecType<SV>::name() +"_D32 write",fsize,estSize,scale);
ir = a.testRead(v1); iret |= a.check(VecType<SV>::name()+"_D32 read",ir,0);
s1 = a.testAddition(v1); iret |= a.check(VecType<SV>::name()+"_D32 after read",s1,sref1,1.E9);
return iret;
}
//--------------------------------------------------------------------------------------
// test of generic Svector
//--------------------------------------------------------------------------------------
template<int Dim>
int testSVector(int ngen, bool testio=false) {
// test the matrix if D2 is not equal to 1
int iret = 0;
VectorTest<Dim> a(ngen);
a.genDataN();
typedef SVector<double, Dim> SV;
std::vector<SV> v1;
v1.reserve(ngen);
double s1 = 0;
//double scale = 1;
double sref1 = 0;
std::string name = "SVector<double," + Util::ToString(Dim) + ">";
a.testCreateSV(v1); iret |= a.check(name+" creation",v1.size(),ngen);
s1 = a.testAdditionSV(v1); iret |= a.check(name+" addition",s1,a.Sum(),Dim*4);
sref1 = s1;
// test the io
double fsize = 0;
int ir = 0;
if (!testio) return iret;
std::string typeName = "ROOT::Math::"+name;
double estSize = ngen*8 * Dim + 10000;
double scale = 0.1 / std::numeric_limits<double>::epsilon();
fsize = a.testWrite(v1,typeName); iret |= a.check(name+" write",fsize,estSize,scale);
ir = a.testRead(v1); iret |= a.check(name+" read",ir,0);
s1 = a.testAdditionSV(v1);
#if(defined __FAST_MATH__ && defined __clang__)
iret |= a.check(name+" after read",s1,sref1, 10);
#else
iret |= a.check(name+" after read",s1,sref1);
#endif
return iret;
}
template<int D1, int D2>
struct RepStd {
typedef typename ROOT::Math::MatRepStd<double,D1,D2> R;
static std::string name() {
return "ROOT::Math::MatRepStd<double," + Util::ToString(D1) + "," + Util::ToString(D2) + "> ";
}
static std::string name32() {
return "ROOT::Math::MatRepStd<Double32_t," + Util::ToString(D1) + "," + Util::ToString(D2) + "> ";
}
static std::string sname() { return ""; }
};
template<int D1>
struct RepSym {
typedef typename ROOT::Math::MatRepSym<double,D1> R;
static std::string name() {
return "ROOT::Math::MatRepSym<double," + Util::ToString(D1) + "> ";
}
static std::string name32() {
return "ROOT::Math::MatRepSym<Double32_t," + Util::ToString(D1) + "> ";
}
static std::string sname() { return "MatRepSym"; }
};
//--------------------------------------------------------------------------------------
// test of generic SMatrix
//--------------------------------------------------------------------------------------
template<int D1, int D2,class Rep >
int testSMatrix(int ngen, bool testio=false) {
// test the matrix if D2 is not equal to 1
int iret = 0;
typedef typename Rep::R R;
typedef SMatrix<double, D1, D2,R> SM;
// in the case of sym matrices SM::kSize is different than R::kSize
// need to use the R::kSize for Dim
const int Dim = R::kSize;
VectorTest<Dim> a(ngen);
a.genDataN();
std::vector<SM> v1;
v1.reserve(ngen);
double s1 = 0;
//double scale = 1;
double sref1 = 0;
std::string name0 = "SMatrix<double," + Util::ToString(D1) + "," + Util::ToString(D2);
std::string name = name0 + "," + Rep::sname() + ">";
a.testCreateSV(v1); iret |= a.check(name+" creation",v1.size(),ngen);
s1 = a.testAdditionSV(v1); iret |= a.check(name+" addition",s1,a.Sum(),Dim*4);
sref1 = s1;
// test the io
if (!testio) return iret;
double fsize = 0;
int ir = 0;
// the full name is needed for sym matrices
std::string typeName = "ROOT::Math::"+name0 + "," + Rep::name() + ">";
double estSize = ngen*8 * Dim + 10000;
double scale = 0.1 / std::numeric_limits<double>::epsilon();
fsize = a.testWrite(v1,typeName); iret |= a.check(name+" write",fsize,estSize,scale);
ir = a.testRead(v1); iret |= a.check(name+" read",ir,0);
s1 = a.testAdditionSV(v1); iret |= a.check(name+" after read",s1,sref1);
//std::cout << "File size = " << fsize << " estimated " << 8 * Dim * ngen << std::endl;
// test storing as Double32_t
// dictionay exist only for square matrices between 3 and 6
if (D1 != D2) return iret;
if (D1 < 3 || D1 > 6) return iret;
double fsize32 = 0;
name0 = "SMatrix<Double32_t," + Util::ToString(D1) + "," + Util::ToString(D2);
name = name0 + "," + Rep::sname() + ">";
typeName = "ROOT::Math::"+name0+ "," + Rep::name32() + ">";
estSize = ngen* 4 * Dim + 60158;
scale = 0.1 / std::numeric_limits<double>::epsilon();
fsize32 = a.testWrite(v1,typeName); iret |= a.check(name+" write",fsize32,estSize,scale);
ir = a.testRead(v1); iret |= a.check(name+" read",ir,0);
//we read back float (scale errors then)
s1 = a.testAdditionSV(v1); iret |= a.check(name+" after read",s1,sref1,1.E9);
return iret;
}
//--------------------------------------------------------------------------------------
// test of a track an object containing vector and matrices
////////////////////////////////////////////////////////////////////////////////
template<class T>
int testTrack(int ngen) {
int iret = 0;
const int Dim = T::kSize;
std::string name = T::Type();
//std::cout << "Dim is " << Dim << std::endl;
VectorTest<Dim> a(ngen);
a.genDataN();
std::vector<T> v1;
v1.reserve(ngen);
double s1 = 0;
//double scale = 1;
double sref1 = 0;
a.testCreateSV(v1); iret |= a.check(name+" creation",v1.size(),ngen);
s1 = a.testAdditionTR(v1); iret |= a.check(name+" addition",s1,a.Sum(),Dim*4);
sref1 = s1;
//std::cout << " reference sum is " << sref1 << std::endl;
double fsize = 0;
int ir = 0;
// the full name is needed for sym matrices
std::string typeName = name;
int wsize = 8;
if (T::IsD32() ) wsize = 4;
double estSize = ngen*wsize * Dim + 10000;
double scale = 0.2 / std::numeric_limits<double>::epsilon();
fsize = a.testWrite(v1,typeName); iret |= a.check(name+" write",fsize,estSize,scale);
ir = a.testRead(v1); iret |= a.check(name+" read",ir,0);
scale = 1;
if (T::IsD32() ) scale = 1.E9; // use float numbers
s1 = a.testAdditionTR(v1); iret |= a.check(name+" after read",s1,sref1,scale);
if (T::IsD32() ) {
// check at double precision type must fail otherwise Double's are stored
bool pdebug = debug; debug = false;
int ret = compare(" ",s1,sref1);
iret |= a.check(name+" Double32 test",ret,1);
debug = pdebug;
}
return iret;
}
int testGenVectors(int ngen,bool io) {
int iret = 0;
std::cout <<"******************************************************************************\n";
std::cout << "\tTest of Physics Vector (GenVector package)\n";
std::cout <<"******************************************************************************\n";
iret |= testVector<XYVector, Polar2DVector, 2>(ngen,io);
iret |= testVector<XYZVector, Polar3DVector, 3>(ngen,io);
iret |= testVector<XYZVector, RhoEtaPhiVector, 3>(ngen,io);
iret |= testVector<XYZVector, RhoZPhiVector, 3>(ngen,io);
iret |= testVector<XYZTVector, PtEtaPhiEVector, 4>(ngen,io);
iret |= testVector<XYZTVector, PtEtaPhiMVector, 4>(ngen,io);
iret |= testVector<XYZTVector, PxPyPzMVector, 4>(ngen,io);
return iret;
}
int testSMatrix(int ngen,bool io) {
int iret = 0;
std::cout <<"******************************************************************************\n";
std::cout << "\tTest of SMatrix package\n";
std::cout <<"******************************************************************************\n";
iret |= testVector34<3>(ngen,io);
iret |= testVector34<4>(ngen,io);
// test now matrices and vectors
iret |= testSVector<6>(ngen,io);
// asymetric matrices we have only 4x3 and 3x4
iret |= testSMatrix<3,4,RepStd<3,4> >(ngen,io);
iret |= testSMatrix<4,3,RepStd<4,3> >(ngen,io);
iret |= testSMatrix<3,3,RepStd<3,3> >(ngen,io);
// sym matrix
iret |= testSMatrix<5,5,RepSym<5> >(ngen,io);
return iret;
}
int testCompositeObj(int ngen) {
std::cout <<"******************************************************************************\n";
std::cout << "\tTest of a Composite Object (containing Vector's and Matrices)\n";
std::cout <<"******************************************************************************\n";
std::cout << "Test Using CINT library\n\n";
// put path relative to LD_LIBRARY_PATH
std::unique_ptr<const char> dynPath(gSystem->DynamicPathName("../test/libTrackMathCoreDict",
/*quiet*/ true));
int iret = -1;
if (dynPath)
iret = gSystem->Load("../test/libTrackMathCoreDict");
if (iret < 0) {
// if not assume running from top ROOT dir (case of roottest)
iret = gSystem->Load("test/libTrackMathCoreDict");
if (iret < 0) {
std::cerr <<"Error Loading libTrackMathCoreDict" << std::endl;
return iret;
}
}
iret = 0;
iret |= testTrack<TrackD>(ngen);
iret |= testTrack<TrackD32>(ngen);
iret |= testTrack<TrackErrD>(ngen);
iret |= testTrack<TrackErrD32>(ngen);
iret |= testTrack<VecTrack<TrackD> >(ngen);
iret |= testTrack<VecTrack<TrackErrD> >(ngen);
return iret;
}
#endif // endif ifndef __CINT__
int stressMathCore(double nscale = 1) {
int iret = 0;
#ifdef __CINT__
std::cout << "Test must be run in compile mode - use ACLIC to compile!!" << std::endl;
gSystem->Load("libMathCore");
gSystem->Load("libTree");
gROOT->ProcessLine(".L stressMathCore.cxx++");
return stressMathCore();
#endif
// iret |= gSystem->Load("libMathCore");
// iret |= gSystem->Load("libMathMore");
// if (iret !=0) return iret;
TBenchmark bm;
bm.Start("stressMathCore");
const int ntest = 10000;
int n = int(nscale*ntest);
//std::cout << "StressMathCore: test number n = " << n << std::endl;
iret |= testStatFunctions(n/10);
bool io = true;
iret |= ( gSystem->Load("libSmatrix") < 0 ); // iret = 0 or = 1 is fine
if (iret != 0) {
std::cerr <<"Error Loading libSmatrix" << std::endl;
io = false;
}
iret |= testGenVectors(n,io);
iret |= testSMatrix(n,io);
if (io) iret |= testCompositeObj(n);
bm.Stop("stressMathCore");
std::cout <<"******************************************************************************\n";
bm.Print("stressMathCore");
const double reftime = 7.1; // needs to be updated // ref time on pcbrun4
double rootmarks = 860 * reftime / bm.GetCpuTime("stressMathCore");
std::cout << " ROOTMARKS = " << rootmarks << " ROOT version: " << gROOT->GetVersion() << "\t"
<< gROOT->GetGitBranch() << "@" << gROOT->GetGitCommit() << std::endl;
std::cout <<"*******************************************************************************\n";
if (iret !=0) std::cerr << "stressMathCore Test Failed !!" << std::endl;
return iret;
}
int main(int argc,const char *argv[]) {
std::string inclRootSys = ("-I" + TROOT::GetRootSys() + "/test").Data();
TROOT::AddExtraInterpreterArgs({inclRootSys});
double nscale = 1;
if (argc > 1) {
nscale = atof(argv[1]);
//nscale = std::pow(10.0,double(scale));
}
return stressMathCore(nscale);
}
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...