Skip to content

Commit

Permalink
reuse allocated memory
Browse files Browse the repository at this point in the history
  • Loading branch information
@dahlend
dahlend committed Jun 10, 2024
1 parent e5012be commit 971ef9d
Show file tree
Hide file tree
Showing 3 changed files with 43 additions and 19 deletions.
17 changes: 16 additions & 1 deletion src/neospy_core/benches/spice.rs
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
extern crate criterion;
use criterion::{black_box, criterion_group, criterion_main, Criterion};
use neospy_core::spice::get_spk_singleton;
use neospy_core::{spice::get_spk_singleton, state::State};
use pprof::criterion::{Output, PProfProfiler};

fn spice_get_raw_state(jd: f64) {
Expand All @@ -10,6 +10,14 @@ fn spice_get_raw_state(jd: f64) {
}
}

fn spice_change_center(mut state: State) {
let spice = get_spk_singleton().try_read().unwrap();
for _ in 0..500 {
spice.try_change_center(&mut state, 10).unwrap();
spice.try_change_center(&mut state, 0).unwrap();
}
}

fn spice_get_state(jd: f64) {
let spice = get_spk_singleton().try_read().unwrap();
for _ in 0..1000 {
Expand All @@ -20,12 +28,19 @@ fn spice_get_state(jd: f64) {
}

pub fn spice_benchmark(c: &mut Criterion) {
let spice = get_spk_singleton().try_read().unwrap();
let state = spice
.try_get_state(5, 2451545.0, 10, neospy_core::frames::Frame::Ecliptic)
.unwrap();
c.bench_function("spice_get_raw_state", |b| {
b.iter(|| spice_get_raw_state(black_box(2451545.0)))
});
c.bench_function("spice_get_state", |b| {
b.iter(|| spice_get_state(black_box(2451545.0)))
});
c.bench_function("spice_change_center", |b| {
b.iter(|| spice_change_center(black_box(state.clone())))
});
}

criterion_group!(name=spice;
Expand Down
38 changes: 23 additions & 15 deletions src/neospy_core/src/propagation/radau.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,6 +104,11 @@ where

cur_b: OMatrix<f64, D, U7>,
g_scratch: OMatrix<f64, D, U7>,

state_scratch: OVector<f64, D>,
state_der_scratch: OVector<f64, D>,
b_scratch: OVector<f64, D>,
eval_scratch: OVector<f64, D>,
}

impl<'a, MType, D: Dim> RadauIntegrator<'a, MType, D>
Expand Down Expand Up @@ -134,6 +139,10 @@ where
cur_state_der_der: Matrix::zeros_generic(dim, U1),
cur_b: Matrix::zeros_generic(dim, U7),
g_scratch: Matrix::zeros_generic(dim, U7),
b_scratch: Matrix::zeros_generic(dim, U1),
state_scratch: Matrix::zeros_generic(dim, U1),
state_der_scratch: Matrix::zeros_generic(dim, U1),
eval_scratch: Matrix::zeros_generic(dim, U1),
};

res.cur_state_der_der = (res.func)(
Expand Down Expand Up @@ -215,15 +224,14 @@ where
/// step guess provided.
///
fn step(&mut self, step_size: f64) -> Result<f64, NEOSpyError> {
let (dim, _) = self.cur_state.shape_generic();
let mut func_eval_tmp = self.cur_state_der_der.clone();
self.g_scratch.fill(0.0);

let mut tmp_state = Matrix::zeros_generic(dim, U1);
let mut tmp_state_der = Matrix::zeros_generic(dim, U1);
self.g_scratch.fill(0.0);
self.state_scratch.fill(0.0);
self.state_der_scratch.fill(0.0);
self.eval_scratch.set_column(0, &self.cur_state_der_der);

for _ in 0..10 {
let last_b = self.cur_b.column(6).clone_owned();
self.b_scratch.set_column(0, &self.cur_b.column(6));
// Calculate b and g
for (idj, gauss_radau_frac) in GAUSS_RADAU_SPACINGS.iter().enumerate().skip(1) {
// the sample point at the Guass-Radau spacings.
Expand All @@ -244,7 +252,7 @@ where
);

izip!(
tmp_state.iter_mut(),
self.state_scratch.iter_mut(),
self.cur_state.iter(),
self.cur_state_der.iter(),
self.cur_state_der_der.iter(),
Expand All @@ -259,7 +267,7 @@ where
});

izip!(
tmp_state_der.iter_mut(),
self.state_der_scratch.iter_mut(),
self.cur_state_der.iter(),
self.cur_state_der_der.iter(),
self.cur_b.row_iter(),
Expand All @@ -272,15 +280,15 @@ where
});

// Evaluate the function at this new intermediate state.
func_eval_tmp = (self.func)(
self.eval_scratch.set_column(0, &(self.func)(
self.cur_time + gauss_radau_frac * step_size,
&tmp_state,
&tmp_state_der,
&self.state_scratch,
&self.state_der_scratch,
&mut self.metadata,
false,
)?;
)?);

let diff = &func_eval_tmp - &self.cur_state_der_der;
let diff = &self.eval_scratch - &self.cur_state_der_der;

// Use the result of that evaluation to update the current G and B
// matrices for the next gauss spacing.
Expand All @@ -302,8 +310,8 @@ where
// Compute the largest b value and compare it to the last b value.
// convergence is decided if B stops changing.
self.g_scratch.mul_to(&C_MAT, &mut self.cur_b);
let b_diff = (self.cur_b.column(6) - last_b).abs();
let func_eval_max = func_eval_tmp.abs().add_scalar(1e-6);
let b_diff = (self.cur_b.column(6) - &self.b_scratch).abs();
let func_eval_max = self.eval_scratch.abs().add_scalar(1e-6);

// This is using the convergence criterion as defined in
// https://arxiv.org/pdf/1409.4779.pdf equation (8)
Expand Down
7 changes: 4 additions & 3 deletions src/neospy_core/src/spice/spk_segments.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -356,14 +356,15 @@ impl SpkSegmentType2 {
) -> Result<([f64; 3], [f64; 3]), NEOSpyError> {
let record_index = ((jd - segment.jd_start) / self.jd_step).floor() as usize;
let record = unsafe { self.records.get_unchecked(record_index) };
let t_mid = record[0];
let t_step = record[1];
let t = (jd - t_mid) / t_step;

let x_coef = unsafe { record.get_unchecked(2..(self.n_coef + 2)) };
let y_coef = unsafe { record.get_unchecked((self.n_coef + 2)..(2 * self.n_coef + 2)) };
let z_coef = unsafe { record.get_unchecked((2 * self.n_coef + 2)..(3 * self.n_coef + 2)) };

let t_mid = unsafe{record.get_unchecked(0)};
let t_step = unsafe{record.get_unchecked(1)};
let t = (jd - t_mid) / t_step;

let (x, vx) = chebyshev_evaluate_both(t, x_coef)?;
let (y, vy) = chebyshev_evaluate_both(t, y_coef)?;
let (z, vz) = chebyshev_evaluate_both(t, z_coef)?;
Expand Down

0 comments on commit 971ef9d

Please sign in to comment.