Skip to content

Commit

Permalink
[Trifocal+P2Pt] microopt: using tight float for last element t
Browse files Browse the repository at this point in the history
  • Loading branch information
rfabbri committed Nov 13, 2023
1 parent 9ea934f commit 09513b5
Show file tree
Hide file tree
Showing 3 changed files with 16 additions and 13 deletions.
2 changes: 1 addition & 1 deletion minus/chicago14a-HxH.hxx
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ HxH(const C<F>* __restrict ux /*x and t*/, const C<F> * __restrict uparams, C<F>
const C<F> &X11 = x[11]; // transl
const C<F> &X12 = x[12]; // transl
const C<F> &X13 = x[13]; // transl
const F &t = (F &) ux[14];
const F &t = *((F *)ux + 28);

const C<F> &X15 = params[0];
const C<F> &X16 = params[1];
Expand Down
2 changes: 1 addition & 1 deletion minus/chicago14a-Hxt.hxx
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@ Hxt(const C<F> * __restrict ux, const C<F> * __restrict uparams, C<F> * __restri
const C<F> &X11 = x[11]; // transl
const C<F> &X12 = x[12]; // transl
const C<F> &X13 = x[13]; // transl
const F &t = (F &) ux[14];
const F &t = *((F *)ux + 28);

const C<F> &X15 = params[0];
const C<F> &X16 = params[1];
Expand Down
25 changes: 14 additions & 11 deletions minus/minus.hxx
Original file line number Diff line number Diff line change
Expand Up @@ -94,10 +94,13 @@ track(const track_settings &s, const C<F> s_sols_u[f::nve*f::nsols], const C<F>
solution *raw_solutions = reinterpret_cast<solution *> (__builtin_assume_aligned(raw_solutions_u,64));
assert(sol_min <= sol_max && sol_max <= f::nsols);
alignas(64) C<F> Hxt[NVEPLUS1 * f::nve];
alignas(64) C<F> x0t0[NVEPLUS1];
alignas(64) C<F> xt[NVEPLUS1];
alignas(64) C<F> dxdt[NVEPLUS1];
alignas(64) F x0t0f[f::nve*2+1];
alignas(64) F xtf[f::nve*2+1];
alignas(64) F dxdtf[f::nve*2+1];
alignas(64) C<F> dxi[f::nve];
C<F> *const x0t0 = (C<F> *) x0t0f;
C<F> *const xt = (C<F> *) xtf;
C<F> *const dxdt = (C<F> *) dxdtf;
C<F> *const x0 = x0t0;
F *const t0 = (F *) (x0t0 + f::nve);
C<F> *const x1t1 = xt;
Expand All @@ -107,20 +110,20 @@ track(const track_settings &s, const C<F> s_sols_u[f::nve*f::nsols], const C<F>
C<F> *const HxH=Hxt;
Map<Matrix<C<F>, f::nve, NVEPLUS1>,Aligned> AA((C<F> *)Hxt,f::nve,NVEPLUS1);
static constexpr F the_smallest_number = 1e-13; // XXX BENCHMARK THIS
typedef minus_array<f::nve,F> v; typedef minus_array<NVEPLUS1,F> vp;
typedef minus_array<f::nve,F> v;

alignas(64) C<F> ycHxt[16];
alignas(64) C<F> ycHxH[13];
// memoization_init() :

C<F> previous[13];
// C<F> previous[13];
const F &t_step = s.init_dt_; // initial step
solution *t_s = raw_solutions + sol_min; // current target solution
const C<F>* __restrict s_s = s_sols + sol_min*f::nve; // current start solution
for (unsigned sol_n = sol_min; sol_n < sol_max; ++sol_n) { // solution loop
t_s->status = PROCESSING;
bool end_zone = false;
v::copy(s_s, x0);
v::fcopy(s_s, x0);
*t0 = 0; *dt = t_step;
char predictor_successes = 0;

Expand All @@ -144,7 +147,7 @@ track(const track_settings &s, const C<F> s_sols_u[f::nve*f::nsols], const C<F>
dx3 := solveHxTimesDXequalsminusHt(x0+(1/2)*dx2*dt,t0+(1/2)*dt);
dx4 := solveHxTimesDXequalsminusHt(x0+dx3*dt,t0+dt);
(1/6)*dt*(dx1+2*dx2+2*dx3+dx4) */
vp::copy(x0t0, xt);
v::fcopy(x0t0, xt);

// dx1
evaluate_Hxt_constants(xt, params, ycHxt);
Expand Down Expand Up @@ -177,7 +180,7 @@ track(const track_settings &s, const C<F> s_sols_u[f::nve*f::nsols], const C<F>

// dx4
v::multiply_scalar_to_self(dxi, *dt);
vp::copy(x0t0, xt);
v::fcopy(x0t0, xt);
v::add_to_self(xt, dxi);
v::multiply_scalar_to_self(dxi, 2);
v::add_to_self(dx4, dxi);
Expand All @@ -193,8 +196,8 @@ track(const track_settings &s, const C<F> s_sols_u[f::nve*f::nsols], const C<F>
// dx4_eigen = (dx4_eigen* *dt + dx1_eigen*2 + dx2_eigen*4 + dx3_eigen*2)*(1./6.);

// make prediction
vp::copy(x0t0, x1t1);
vp::add_to_self(x1t1, dxdt);
v::fcopy(x0t0, x1t1);
v::fadd_to_self(x1t1, dxdt);

/// CORRECTOR ///
char n_corr_steps = 0;
Expand Down Expand Up @@ -237,7 +240,7 @@ track(const track_settings &s, const C<F> s_sols_u[f::nve*f::nsols], const C<F>
++predictor_successes;
// std::swap(x1t1,x0t0);
// x0 = x0t0; t0 = (F *) (x0t0 + f::nve); xt = x1t1;
vp::copy(x1t1, x0t0);
v::fcopy(x1t1, x0t0);
if (predictor_successes >= s.num_successes_before_increase_) {
predictor_successes = 0;
*dt *= s.dt_increase_factor_;
Expand Down

0 comments on commit 09513b5

Please sign in to comment.