Skip to content

Commit

Permalink
Update helixAtZ to use small calculation loops
Browse files Browse the repository at this point in the history
  • Loading branch information
@gartung
gartung committed May 13, 2022
1 parent 3d704ba commit d6ae2c9
Show file tree
Hide file tree
Showing 2 changed files with 238 additions and 69 deletions.
294 changes: 231 additions & 63 deletions RecoTracker/MkFitCore/src/PropagationMPlex.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -708,21 +708,58 @@ namespace mkfit {
errorProp(n, 3, 3) = 1.f;
errorProp(n, 4, 4) = 1.f;
errorProp(n, 5, 5) = 1.f;

const float zout = msZ.constAt(n, 0, 0);

const float zin = inPar.constAt(n, 2, 0);
const float ipt = inPar.constAt(n, 3, 0);
const float phiin = inPar.constAt(n, 4, 0);
const float theta = inPar.constAt(n, 5, 0);

const float k =
inChg.constAt(n, 0, 0) * 100.f /
(-Const::sol * (pflags.use_param_b_field
? Config::bFieldFromZR(zin, hipo(inPar.constAt(n, 0, 0), inPar.constAt(n, 1, 0)))
: Config::Bfield));
const float kinv = 1.f / k;

}
float zout[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
//initialize erroProp to identity matrix, except element 2,2 which is zero
zout[n] = msZ.constAt(n, 0, 0);
}
float zin[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
//initialize erroProp to identity matrix, except element 2,2 which is zero
zin[n] = inPar.constAt(n, 2, 0);
}
float ipt[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
//initialize erroProp to identity matrix, except element 2,2 which is zero
ipt[n] = inPar.constAt(n, 3, 0);
}
float phiin[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
//initialize erroProp to identity matrix, except element 2,2 which is zero
phiin[n] = inPar.constAt(n, 4, 0);
}
float theta[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
//initialize erroProp to identity matrix, except element 2,2 which is zero
theta[n] = inPar.constAt(n, 5, 0);
}
float k[NN];
if (pflags.use_param_b_field) {
#pragma omp simd
for (int n = 0; n < NN; ++n) {
k[n] = inChg.constAt(n, 0, 0) * 100.f / (-Const::sol *
Config::bFieldFromZR(zin[n], hipo(inPar.constAt(n, 0, 0), inPar.constAt(n, 1, 0))));
}
} else {
#pragma omp simd
for (int n = 0; n < NN; ++n) {
k[n] = inChg.constAt(n, 0, 0) * 100.f / (-Const::sol *
Config::Bfield);
}
}
float kinv[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
kinv[n] = 1.f / k[n];
}
#pragma omp simd
for (int n = 0; n < NN; ++n) {
dprint_np(n,
std::endl
<< "input parameters"
Expand All @@ -732,66 +769,194 @@ namespace mkfit {
<< " inPar.constAt(n, 3, 0)=" << std::setprecision(9) << inPar.constAt(n, 3, 0)
<< " inPar.constAt(n, 4, 0)=" << std::setprecision(9) << inPar.constAt(n, 4, 0)
<< " inPar.constAt(n, 5, 0)=" << std::setprecision(9) << inPar.constAt(n, 5, 0));

const float pt = 1.f / ipt;

float cosahTmp = 0., sinahTmp = 0.;
}
float pt[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
pt[n] = 1.f / ipt[n];
}

//no trig approx here, phi can be large
const float cosP = std::cos(phiin), sinP = std::sin(phiin);
const float cosT = std::cos(theta), sinT = std::sin(theta);
const float tanT = sinT / cosT;
const float icos2T = 1.f / (cosT * cosT);
const float pxin = cosP * pt;
const float pyin = sinP * pt;
float cosP[NN];
float sinP[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
cosP[n] = std::cos(phiin[n]);
}
#pragma omp simd
for (int n = 0; n < NN; ++n) {
sinP[n] = std::sin(phiin[n]);
}
float cosT[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
cosT[n] = std::cos(theta[n]);
}
float sinT[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
sinT[n] = std::sin(theta[n]);
}
float tanT[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
tanT[n] = sinT[n] / cosT[n];
}
float icos2T[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
icos2T[n] = 1.f / (cosT[n] * cosT[n]);
}
float pxin[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
pxin[n] = cosP[n] * pt[n];
}
float pyin[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
pyin[n] = sinP[n] * pt[n];
}
#pragma omp simd
for (int n = 0; n < NN; ++n) {

//fixme, make this printout useful for propagation to z
dprint_np(n,
std::endl
<< "k=" << std::setprecision(9) << k << " pxin=" << std::setprecision(9) << pxin
<< " pyin=" << std::setprecision(9) << pyin << " cosP=" << std::setprecision(9) << cosP
<< " sinP=" << std::setprecision(9) << sinP << " pt=" << std::setprecision(9) << pt);

const float deltaZ = zout - zin;
const float alpha = deltaZ * tanT * ipt * kinv;

if (Config::useTrigApprox) {
sincos4(alpha * 0.5f, sinahTmp, cosahTmp);
} else {
cosahTmp = std::cos(alpha * 0.5f);
sinahTmp = std::sin(alpha * 0.5f);
<< "k=" << std::setprecision(9) << k[n] << " pxin=" << std::setprecision(9) << pxin[n]
<< " pyin=" << std::setprecision(9) << pyin[n] << " cosP=" << std::setprecision(9) << cosP[n]
<< " sinP=" << std::setprecision(9) << sinP[n] << " pt=" << std::setprecision(9) << pt[n]);
}
float deltaZ[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
deltaZ[n] = zout[n] - zin[n];
}
float alpha[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
alpha[n] = deltaZ[n] * tanT[n] * ipt[n] * kinv[n];
}
float cosahTmp[NN] = {0.};
float sinahTmp[NN] = {0.};
if (Config::useTrigApprox) {
#pragma omp simd
for (int n = 0; n < NN; ++n) {
sincos4(alpha[n] * 0.5f, sinahTmp[n], cosahTmp[n]);
}
const float cosah = cosahTmp;
const float sinah = sinahTmp;
const float cosa = 1.f - 2.f * sinah * sinah;
const float sina = 2.f * sinah * cosah;

} else {
#pragma omp simd
for (int n = 0; n < NN; ++n) {
cosahTmp[n] = std::cos(alpha[n] * 0.5f);
}
#pragma omp simd
for (int n = 0; n < NN; ++n) {
sinahTmp[n] = std::sin(alpha[n] * 0.5f);
}
}
float cosah[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
cosah[n] = cosahTmp[n];
}
float sinah[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
sinah[n] = sinahTmp[n];
}
float cosa[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
cosa[n] = 1.f - 2.f * sinah[n] * sinah[n];
}
float sina[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
sina[n] = 2.f * sinah[n] * cosah[n];
}
//update parameters
outPar.At(n, 0, 0) = outPar.At(n, 0, 0) + 2.f * k * sinah * (pxin * cosah - pyin * sinah);
outPar.At(n, 1, 0) = outPar.At(n, 1, 0) + 2.f * k * sinah * (pyin * cosah + pxin * sinah);
outPar.At(n, 2, 0) = zout;
outPar.At(n, 4, 0) = phiin + alpha;
#pragma omp simd
for (int n = 0; n < NN; ++n) {
outPar.At(n, 0, 0) = outPar.At(n, 0, 0) + 2.f * k[n] * sinah[n] * (pxin[n] * cosah[n] - pyin[n] * sinah[n]);
}
#pragma omp simd
for (int n = 0; n < NN; ++n) {
outPar.At(n, 1, 0) = outPar.At(n, 1, 0) + 2.f * k[n] * sinah[n] * (pyin[n] * cosah[n] + pxin[n] * sinah[n]);
}
#pragma omp simd
for (int n = 0; n < NN; ++n) {
outPar.At(n, 2, 0) = zout[n];
}
#pragma omp simd
for (int n = 0; n < NN; ++n) {
outPar.At(n, 4, 0) = phiin[n] + alpha[n];
}

#pragma omp simd
for (int n = 0; n < NN; ++n) {
dprint_np(n,
std::endl
<< "outPar.At(n, 0, 0)=" << outPar.At(n, 0, 0) << " outPar.At(n, 1, 0)=" << outPar.At(n, 1, 0)
<< " pxin=" << pxin << " pyin=" << pyin);

const float pxcaMpysa = pxin * cosa - pyin * sina;
errorProp(n, 0, 2) = -tanT * ipt * pxcaMpysa;
errorProp(n, 0, 3) = k * pt * pt * (cosP * (alpha * cosa - sina) + sinP * 2.f * sinah * (sinah - alpha * cosah));
errorProp(n, 0, 4) = -2 * k * pt * sinah * (sinP * cosah + cosP * sinah);
errorProp(n, 0, 5) = deltaZ * ipt * pxcaMpysa * icos2T;

const float pycaPpxsa = pyin * cosa + pxin * sina;
errorProp(n, 1, 2) = -tanT * ipt * pycaPpxsa;
errorProp(n, 1, 3) = k * pt * pt * (sinP * (alpha * cosa - sina) - cosP * 2.f * sinah * (sinah - alpha * cosah));
errorProp(n, 1, 4) = 2 * k * pt * sinah * (cosP * cosah - sinP * sinah);
errorProp(n, 1, 5) = deltaZ * ipt * pycaPpxsa * icos2T;
<< " pxin=" << pxin[n] << " pyin=" << pyin[n]);
}
float pxcaMpysa[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
pxcaMpysa[n] = pxin[n] * cosa[n] - pyin[n] * sina[n];
}
#pragma omp simd
for (int n = 0; n < NN; ++n) {
errorProp(n, 0, 2) = -tanT[n] * ipt[n] * pxcaMpysa[n];
}
#pragma omp simd
for (int n = 0; n < NN; ++n) {
errorProp(n, 0, 3) = k[n] * pt[n] * pt[n] * (cosP[n] * (alpha[n] * cosa[n] - sina[n]) + sinP[n] * 2.f * sinah[n] * (sinah[n] - alpha[n] * cosah[n]));
}
#pragma omp simd
for (int n = 0; n < NN; ++n) {
errorProp(n, 0, 4) = -2.f * k[n] * pt[n] * sinah[n] * (sinP[n] * cosah[n] + cosP[n] * sinah[n]);
}
#pragma omp simd
for (int n = 0; n < NN; ++n) {
errorProp(n, 0, 5) = deltaZ[n] * ipt[n] * pxcaMpysa[n] * icos2T[n];
}
float pycaPpxsa[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
pycaPpxsa[n] = pyin[n] * cosa[n] + pxin[n] * sina[n];
}
#pragma omp simd
for (int n = 0; n < NN; ++n) {
errorProp(n, 1, 2) = -tanT[n] * ipt[n] * pycaPpxsa[n];
}
#pragma omp simd
for (int n = 0; n < NN; ++n) {
errorProp(n, 1, 3) = k[n] * pt[n] * pt[n] * (sinP[n] * (alpha[n] * cosa[n] - sina[n]) - cosP[n] * 2.f * sinah[n] * (sinah[n] - alpha[n] * cosah[n]));
}
#pragma omp simd
for (int n = 0; n < NN; ++n) {
errorProp(n, 1, 4) = 2.f * k[n] * pt[n] * sinah[n] * (cosP[n] * cosah[n] - sinP[n] * sinah[n]);
}
#pragma omp simd
for (int n = 0; n < NN; ++n) {
errorProp(n, 1, 5) = deltaZ[n] * ipt[n] * pycaPpxsa[n] * icos2T[n];
}

errorProp(n, 4, 2) = -ipt * tanT * kinv;
errorProp(n, 4, 3) = tanT * deltaZ * kinv;
errorProp(n, 4, 5) = ipt * deltaZ * kinv * icos2T;
#pragma omp simd
for (int n = 0; n < NN; ++n) {
errorProp(n, 4, 2) = -ipt[n] * tanT[n] * kinv[n];
}
#pragma omp simd
for (int n = 0; n < NN; ++n) {
errorProp(n, 4, 3) = tanT[n] * deltaZ[n] * kinv[n];
}
#pragma omp simd
for (int n = 0; n < NN; ++n) {
errorProp(n, 4, 5) = ipt[n] * deltaZ[n] * kinv[n] * icos2T[n];
}

#pragma omp simd
for (int n = 0; n < NN; ++n) {
dprint_np(
n,
"propagation end, dump parameters"
Expand All @@ -802,8 +967,11 @@ namespace mkfit {
<< 1. / (outPar.At(n, 3, 0) * tan(outPar.At(n, 5, 0)))
<< " r=" << std::sqrt(outPar.At(n, 0, 0) * outPar.At(n, 0, 0) + outPar.At(n, 1, 0) * outPar.At(n, 1, 0))
<< " pT=" << 1. / std::abs(outPar.At(n, 3, 0)) << std::endl);
}

#ifdef DEBUG
#pragma omp simd
for (int n = 0; n < NN; ++n) {
if (n < N_proc) {
dmutex_guard;
std::cout << n << ": jacobian" << std::endl;
Expand Down Expand Up @@ -850,8 +1018,8 @@ namespace mkfit {
errorProp(n, 5, 4),
errorProp(n, 5, 5));
}
#endif
}
#endif
}

//==============================================================================
Expand Down
13 changes: 7 additions & 6 deletions RecoTracker/MkFitCore/standalone/Makefile.config
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,16 +51,16 @@ OPT := -g -O3

# 4. Vectorization settings
ifdef AVX_512
VEC_GCC := -march=native # -fopt-info-vec -mavx512f -mavx512cd -fopt-info-vec-all
VEC_GCC := -march=native # -fopt-info-vec -mavx512f -mavx512cd
VEC_ICC := -xHost -qopt-zmm-usage=high # -xcore-avx512
else ifdef AVX2
VEC_GCC := -mavx2 -mfma #-fopt-info-vec-all
VEC_GCC := -mavx2 -mfma
VEC_ICC := -mavx2 -mfma
else ifdef SSE3
VEC_GCC := -msse3 # -fopt-info-vec-all
VEC_GCC := -msse3
VEC_ICC := -msse3
else
VEC_GCC := -mavx # -fopt-info-vec-all
VEC_GCC := -mavx
VEC_ICC := -mavx
endif

Expand Down Expand Up @@ -129,8 +129,9 @@ endif

ifeq ($(CXX), g++)
CXXFLAGS += -std=c++1z -ftree-vectorize -Werror=main -Werror=pointer-arith -Werror=overlength-strings -Wno-vla -Werror=overflow -Wstrict-overflow -Werror=array-bounds -Werror=format-contains-nul -Werror=type-limits -fvisibility-inlines-hidden -fno-math-errno --param vect-max-version-for-alias-checks=50 -Xassembler --compress-debug-sections -felide-constructors -fmessage-length=0 -Wall -Wno-non-template-friend -Wno-long-long -Wreturn-type -Wunused -Wparentheses -Wno-deprecated -Werror=return-type -Werror=missing-braces -Werror=unused-value -Werror=address -Werror=format -Werror=sign-compare -Werror=write-strings -Werror=delete-non-virtual-dtor -Wstrict-aliasing -Werror=narrowing -Werror=unused-but-set-variable -Werror=reorder -Werror=unused-variable -Werror=conversion-null -Werror=return-local-addr -Wnon-virtual-dtor -Werror=switch -fdiagnostics-show-option -Wno-unused-local-typedefs -Wno-attributes -Wno-psabi
CXXFLAGS += -fdiagnostics-color=auto -fdiagnostics-show-option -pthread -pipe -fopenmp-simd -ffast-math
#CXXFLAGS += -mveclibabi=svml -lsvml -L/cvmfs/projects.cern.ch/intelsw/oneAPI/linux/x86_64/2022/compiler/latest/linux/compiler/lib/intel64 -Wl,-rpath=/cvmfs/projects.cern.ch/intelsw/oneAPI/linux/x86_64/2022/compiler/latest/linux/compiler/lib/intel64 -funsafe-math-optimizations
CXXFLAGS += -fdiagnostics-color=auto -fdiagnostics-show-option -pthread -pipe -fopenmp-simd
CXXFLAGS += -ffast-math
CXXFLAGS += -mveclibabi=svml -lsvml -L/cvmfs/projects.cern.ch/intelsw/oneAPI/linux/x86_64/2022/compiler/latest/linux/compiler/lib/intel64 -Wl,-rpath=/cvmfs/projects.cern.ch/intelsw/oneAPI/linux/x86_64/2022/compiler/latest/linux/compiler/lib/intel64 -funsafe-math-optimizations
endif

# Try to find a new enough TBB
Expand Down

0 comments on commit d6ae2c9

Please sign in to comment.