From 971ef9d7e4b6150a9931675477f2b79155e947e2 Mon Sep 17 00:00:00 2001 From: dahlend Date: Mon, 10 Jun 2024 14:30:32 -0700 Subject: [PATCH] reuse allocated memory --- src/neospy_core/benches/spice.rs | 17 +++++++++- src/neospy_core/src/propagation/radau.rs | 38 ++++++++++++++--------- src/neospy_core/src/spice/spk_segments.rs | 7 +++-- 3 files changed, 43 insertions(+), 19 deletions(-) diff --git a/src/neospy_core/benches/spice.rs b/src/neospy_core/benches/spice.rs index d886a5a..48afa06 100644 --- a/src/neospy_core/benches/spice.rs +++ b/src/neospy_core/benches/spice.rs @@ -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) { @@ -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 { @@ -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; diff --git a/src/neospy_core/src/propagation/radau.rs b/src/neospy_core/src/propagation/radau.rs index 1e88a8b..6c1d5ba 100644 --- a/src/neospy_core/src/propagation/radau.rs +++ b/src/neospy_core/src/propagation/radau.rs @@ -104,6 +104,11 @@ where cur_b: OMatrix, g_scratch: OMatrix, + + state_scratch: OVector, + state_der_scratch: OVector, + b_scratch: OVector, + eval_scratch: OVector, } impl<'a, MType, D: Dim> RadauIntegrator<'a, MType, D> @@ -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)( @@ -215,15 +224,14 @@ where /// step guess provided. /// fn step(&mut self, step_size: f64) -> Result { - 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. @@ -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(), @@ -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(), @@ -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. @@ -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) diff --git a/src/neospy_core/src/spice/spk_segments.rs b/src/neospy_core/src/spice/spk_segments.rs index c3ee2bf..5a8257d 100644 --- a/src/neospy_core/src/spice/spk_segments.rs +++ b/src/neospy_core/src/spice/spk_segments.rs @@ -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)?;