Skip to content

Commit

Permalink
Squashed commit of the following:
Browse files Browse the repository at this point in the history 
commit a565f52
Author: Patrick Gartung <[email protected]>
Date:   Thu May 19 15:23:44 2022 +0200

    Remove fast-math flag and attribute needed to prevent segfault when fast-math is used

commit c3fbce9
Author: Patrick Gartung <[email protected]>
Date:   Mon May 16 17:28:41 2022 +0200

    Code format

commit c388ab3
Author: Patrick Gartung <[email protected]>
Date:   Mon May 16 17:27:00 2022 +0200

    Update location of #endif after enabling debug. #pragma omp simd not needed for debug statements. Move debug statements to own loop.

commit 468b476
Author: Patrick Gartung <[email protected]>
Date:   Sat May 14 19:51:00 2022 +0200

    No need to link to libmvec. It is called from libm on el8. Always use -ffastmath

commit 0ef9036
Author: Patrick Gartung <[email protected]>
Date:   Fri May 13 20:01:43 2022 +0200

    Only works on el8_amd64 arch

commit 8866ca6
Author: Patrick Gartung <[email protected]>
Date:   Fri May 13 16:10:23 2022 +0200

    Update for macro name change

commit fe4a38f
Author: Patrick Gartung <[email protected]>
Date:   Thu May 12 21:43:02 2022 +0200

    code format?

commit 8326129
Author: Patrick Gartung <[email protected]>
Date:   Thu May 12 21:18:44 2022 +0200

	    Add CMS_UNROLL_LOOP defines for Intel compilers and use with pragme omp simd instructions that cause internal icc error

commit 056f7ef
Author: Patrick Gartung <[email protected]>
Date:   Thu May 12 03:42:01 2022 +0200

    Add CMS_UNROLL_LOOP_COUNT(Config::Niter) and initialize arrays with += and -= assignment statements to 0. first

commit 31e2c83
Author: Patrick Gartung <[email protected]>
Date:   Wed May 11 05:16:27 2022 +0200

    Add CMS_UNROLL_LOOP_COUNT(Config::Niter)
    and #include "FWCore/Utilities/interface/CMSUnrollLoop.h"

commit 8eab5b1
Author: Patrick Gartung <[email protected]>
Date:   Wed May 11 04:57:47 2022 +0200

    Change atrribute from no-fast-math to match-errno

commit 688e905
Author: Patrick Gartung <[email protected]>
Date:   Tue May 10 17:37:25 2022 +0200

    Code format

commit ed58228
Author: Patrick Gartung <[email protected]>
Date:   Tue May 10 17:23:08 2022 +0200

    Change attribute from no-inline to no-fast-math. Back out changes to applyMaterialEffects since it does not vectorize with the conditional on radL.

commit c19f9c8
Author: Patrick Gartung <[email protected]>
Date:   Mon May 9 19:49:15 2022 +0200

    code format

commit 373ea0d
Author: Patrick Gartung <[email protected]>
Date:   Mon May 9 19:02:48 2022 +0200

    missed a few initialization of variable length arrays

commit 753934c
Author: Patrick Gartung <[email protected]>
Date:   Mon May 9 18:51:13 2022 +0200

    for all el8 archs

commit 2eb25ef
Author: Patrick Gartung <[email protected]>
Date:   Mon May 9 18:31:35 2022 +0200

    clang does not allow initialization of variable length arrays

commit 7f23762
Author: Patrick Gartung <[email protected]>
Date:   Mon May 9 17:41:29 2022 +0200

    Add fast-math and mvec when compiling under el8

commit 5f6efe2
Author: Patrick Gartung <[email protected]>
Date:   Mon May 9 16:07:08 2022 +0200

    Comment out changes that are CERN specific

commit 1ba6e06
Author: Patrick Gartung <[email protected]>
Date:   Thu May 5 21:25:33 2022 +0200

    Try to break loop in applyMaterialEffects into smaller loops
    Combine arithmetic ops into one loop.
    Add -lmvec which is avaiable on RHEL8

commit dbc748c
Author: Patrick Gartung <[email protected]>
Date:   Wed May 4 20:02:05 2022 +0200

    Incorporate Dan Riley's suggested changes

commit d6ae2c9
Author: Patrick Gartung <[email protected]>
Date:   Tue May 3 20:20:47 2022 +0200

    Update helixAtZ to use small calculation loops

commit 3d704ba
Author: Patrick Gartung <[email protected]>
Date:   Fri Apr 29 21:43:00 2022 +0200

    Work in progress

commit 31550b4
Author: Patrick Gartung <[email protected]>
Date:   Fri Apr 29 20:03:09 2022 +0200

    Fixups

commit f10ef52
Author: Patrick Gartung <[email protected]>
Date:   Fri Apr 29 00:57:47 2022 +0200

    WIP
  • Loading branch information
@gartung
gartung committed Jun 6, 2022
1 parent 2a721c8 commit 45977ea
Show file tree
Hide file tree
Showing 2 changed files with 531 additions and 188 deletions.
239 changes: 181 additions & 58 deletions RecoTracker/MkFitCore/src/PropagationMPlex.cc
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
#include "FWCore/Utilities/interface/CMSUnrollLoop.h"

#include "MaterialEffects.h"
#include "PropagationMPlex.h"

Expand Down Expand Up @@ -331,6 +333,7 @@ namespace mkfit {
float pxin = cosP / ipt;
float pyin = sinP / ipt;

CMS_UNROLL_LOOP_COUNT(Config::Niter)
for (int i = 0; i < Config::Niter; ++i) {
dprint_np(n,
std::endl
Expand All @@ -343,7 +346,7 @@ namespace mkfit {
const float ialpha = (r - r0) * ipt / k;
//alpha+=ialpha;

if (Config::useTrigApprox) {
if constexpr (Config::useTrigApprox) {
sincos4(ialpha * 0.5f, sinah, cosah);
} else {
cosah = std::cos(ialpha * 0.5f);
Expand Down Expand Up @@ -708,21 +711,44 @@ namespace mkfit {
errorProp(n, 3, 3) = 1.f;
errorProp(n, 4, 4) = 1.f;
errorProp(n, 5, 5) = 1.f;
}
float zout[NN];
float zin[NN];
float ipt[NN];
float phiin[NN];
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
zout[n] = msZ.constAt(n, 0, 0);
zin[n] = inPar.constAt(n, 2, 0);
ipt[n] = inPar.constAt(n, 3, 0);
phiin[n] = inPar.constAt(n, 4, 0);
theta[n] = inPar.constAt(n, 5, 0);
}

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);
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);
}
}

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 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 +758,160 @@ 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));
}

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

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

float cosahTmp = 0., sinahTmp = 0.;
//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;
#pragma omp simd
for (int n = 0; n < NN; ++n) {
sinP[n] = std::sin(phiin[n]);
}

float cosT[NN];
float sinT[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
cosT[n] = std::cos(theta[n]);
}

#pragma omp simd
for (int n = 0; n < NN; ++n) {
sinT[n] = std::sin(theta[n]);
}

float tanT[NN];
float icos2T[NN];
float pxin[NN];
float pyin[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
tanT[n] = sinT[n] / cosT[n];
icos2T[n] = 1.f / (cosT[n] * cosT[n]);
pxin[n] = cosP[n] * pt[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];
float alpha[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
deltaZ[n] = zout[n] - zin[n];
alpha[n] = deltaZ[n] * tanT[n] * ipt[n] * kinv[n];
}

float cosahTmp[NN];
float sinahTmp[NN];
if constexpr (Config::useTrigApprox) {
#ifndef __INTEL_COMPILER
#pragma omp simd
#endif
for (int n = 0; n < NN; ++n) {
sincos4(alpha[n] * 0.5f, sinahTmp[n], cosahTmp[n]);
}
} else {
#ifndef __INTEL_COMPILER
#pragma omp simd
#endif
for (int n = 0; n < NN; ++n) {
cosahTmp[n] = std::cos(alpha[n] * 0.5f);
}
#ifndef __INTEL_COMPILER
#pragma omp simd
#endif
for (int n = 0; n < NN; ++n) {
sinahTmp[n] = std::sin(alpha[n] * 0.5f);
}
const float cosah = cosahTmp;
const float sinah = sinahTmp;
const float cosa = 1.f - 2.f * sinah * sinah;
const float sina = 2.f * sinah * cosah;
}

//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;
float cosah[NN];
float sinah[NN];
float cosa[NN];
float sina[NN];
#pragma omp simd
for (int n = 0; n < NN; ++n) {
cosah[n] = cosahTmp[n];
sinah[n] = sinahTmp[n];
cosa[n] = 1.f - 2.f * sinah[n] * sinah[n];
sina[n] = 2.f * sinah[n] * cosah[n];
}
//update parameters
#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]);
outPar.At(n, 1, 0) = outPar.At(n, 1, 0) + 2.f * k[n] * sinah[n] * (pyin[n] * cosah[n] + pxin[n] * sinah[n]);
outPar.At(n, 2, 0) = zout[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);
<< " 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];
}

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;
#pragma omp simd
for (int n = 0; n < NN; ++n) {
errorProp(n, 0, 2) = -tanT[n] * ipt[n] * pxcaMpysa[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]));
errorProp(n, 0, 4) = -2.f * k[n] * pt[n] * sinah[n] * (sinP[n] * cosah[n] + cosP[n] * sinah[n]);
errorProp(n, 0, 5) = deltaZ[n] * ipt[n] * pxcaMpysa[n] * icos2T[n];
}

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

#pragma omp simd
for (int n = 0; n < NN; ++n) {
errorProp(n, 4, 2) = -ipt[n] * tanT[n] * kinv[n];
errorProp(n, 4, 3) = tanT[n] * deltaZ[n] * kinv[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 +922,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 +973,8 @@ namespace mkfit {
errorProp(n, 5, 4),
errorProp(n, 5, 5));
}
#endif
}
#endif
}

//==============================================================================
Expand Down
Loading

0 comments on commit 45977ea

Please sign in to comment.