Skip to content

avoid RK4Vec reallocations with a pool #1523

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
baltzell wants to merge 2 commits into
base: development
Choose a base branch
from
Open

avoid RK4Vec reallocations with a pool #1523

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
d9ada6d
use a pool of RK4Vec
baltzell Oct 16, 2026
b8dd1ca
Merge branch 'development' into rungekutta-objects
baltzell Oct 16, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
79 changes: 63 additions & 16 deletions ...on-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/utilities/RungeKuttaDoca.java
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
package org.jlab.clas.tracking.utilities;

import java.util.concurrent.ConcurrentLinkedQueue;
import org.jlab.clas.pdg.PhysicsConstants;
import org.jlab.jnp.matrix.*;
import org.jlab.clas.swimtools.Swim;
Expand All @@ -23,12 +24,11 @@ public void SwimToZ(int sector, StateVec vec, Swim swim, double z0, float[] bf)
swim.Bfield(sector, vec.x, vec.y, vec.z, bf);

vec.B = Math.sqrt(bf[0]*bf[0] + bf[1]*bf[1] + bf[2]*bf[2]);
double s = vec.B;
final double travelSign = Math.signum(z0 - vec.z);
double BatMeas = 0;

while(travelSign * vec.z < travelSign * z0) {
s = travelSign * stepSize;
double s = travelSign * stepSize;
if (travelSign*(vec.z+s) > travelSign*z0) s = travelSign*Math.abs(z0-vec.z);

this.RK4transport(sector, s, swim, vec);
Expand Down Expand Up @@ -82,8 +82,8 @@ private double[] RK4step(int sector, double z0, double h, Swim swim, double[] sI
public void RK4transport(int sector, double h, Swim swim, Matrix cMat, StateVec vec) {
// Set initial state and Jacobian.
double qv = vec.Q*v;
RK4Vec s0 = new RK4Vec(vec);
RK4Vec sNull = new RK4Vec();
RK4Vec s0 = RK4VEC_POOL.poll(vec);
RK4Vec sNull = RK4VEC_POOL.poll();

// Perform steps.
RK4Vec s1 = RK4step(sector, vec.z, 0, swim, s0, sNull, qv);
Expand All @@ -94,6 +94,11 @@ public void RK4transport(int sector, double h, Swim swim, Matrix cMat, StateVec
// Compute and set final state and covariance matrix.
RK4Vec sF = computeFinalState(h, s0, s1, s2, s3, s4);

RK4VEC_POOL.offer(s1);
RK4VEC_POOL.offer(s2);
RK4VEC_POOL.offer(s3);
RK4VEC_POOL.offer(s4);

vec.x = sF.x;
vec.y = sF.y;
vec.tx = sF.tx;
Expand All @@ -102,12 +107,17 @@ public void RK4transport(int sector, double h, Swim swim, Matrix cMat, StateVec
vec.B = Math.sqrt(_b[0]*_b[0]+_b[1]*_b[1]+_b[2]*_b[2]);
vec.deltaPath += Math.sqrt((s0.x-vec.x)*(s0.x-vec.x) + (s0.y-vec.y)*(s0.y-vec.y) + h*h);
vec.CM.set(computeCovMat(cMat, sF));

RK4VEC_POOL.offer(sF);

RK4VEC_POOL.offer(s0);
RK4VEC_POOL.offer(sNull);
}

/** Perform one RK4 step, updating the covariance matrix. */
private RK4Vec RK4step(int sector, double z0, double h, Swim swim,
RK4Vec sInit, RK4Vec sPrev, double qv) {
RK4Vec sNext = new RK4Vec();
RK4Vec sNext = RK4VEC_POOL.poll();
swim.Bfield(sector, sInit.x + h*sPrev.x, sInit.y + h*sPrev.y, z0 + h, _b);

// State.
Expand Down Expand Up @@ -148,7 +158,7 @@ private RK4Vec RK4step(int sector, double z0, double h, Swim swim,
/** Compute the final state for each entry in the internal state matrix. */
private RK4Vec computeFinalState(double h, RK4Vec s0, RK4Vec s1, RK4Vec s2, RK4Vec s3,
RK4Vec s4) {
RK4Vec sF = new RK4Vec();
RK4Vec sF = RK4VEC_POOL.poll();

sF.x = s0.x + this.RK4(s1.x, s2.x, s3.x, s4.x, h);
sF.dxdtx0 = this.RK4(s1.dxdtx0, s2.dxdtx0, s3.dxdtx0, s4.dxdtx0, h);
Expand Down Expand Up @@ -329,8 +339,34 @@ private double dty_dq0(double qv, double Ay, double dAy_dtx, double dAy_dty,
return v*Ay + qv * (dAy_dtx*dtx_dq0 + dAy_dty*dty_dq0);
}

static class RK4VecPool extends ConcurrentLinkedQueue<RK4Vec> {
public RK4VecPool(int size) {
super();
for (int i=0; i<size; i++) add(new RK4Vec());
}
@Override
public RK4Vec poll() {
RK4Vec v = super.poll();
v.reset();
return v;
}
public RK4Vec poll(StateVec vec) {
RK4Vec v = poll();
v.x = vec.x;
v.y = vec.y;
v.tx = vec.tx;
v.ty = vec.ty;
v.dtxdtx0 = 1;
v.dtydty0 = 1;
return v;
}
}

// At 16 doubles each, 10k is about 10 MB:
static final RK4VecPool RK4VEC_POOL =new RK4VecPool(10000);

/** State vector and its derivatives used internally. */
private class RK4Vec {
static class RK4Vec {
public double x = 0;
public double dxdtx0 = 0;
public double dxdty0 = 0;
Expand All @@ -351,15 +387,26 @@ private class RK4Vec {
public double dtydty0 = 0;
public double dtydq0 = 0;

RK4Vec() {}
RK4Vec(StateVec vec) {
this.x = vec.x;
this.y = vec.y;
this.tx = vec.tx;
this.ty = vec.ty;

this.dtxdtx0 = 1;
this.dtydty0 = 1;
void reset() {
x = 0;
dxdtx0 = 0;
dxdty0 = 0;
dxdq0 = 0;

y = 0;
dydtx0 = 0;
dydty0 = 0;
dydq0 = 0;

tx = 0;
dtxdtx0 = 0;
dtxdty0 = 0;
dtxdq0 = 0;

ty = 0;
dtydtx0 = 0;
dtydty0 = 0;
dtydq0 = 0;
}
}
}
Loading