Skip to content

SIMD acceleration for spatial interaction kernels #590

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

Merged
bhaller merged 7 commits into from
Dec 17, 2025
Merged

SIMD acceleration for spatial interaction kernels #590

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
30fadda
Add SIMD acceleration for spatial interaction kernels (Gaussian/Expon…
andrewkern Dec 16, 2025
f62d7f3
Enable SIMD path for 2D spatial interaction kernels (except Fixed)
andrewkern Dec 17, 2025
9bbd5fa
Add SIMD acceleration for Student's T spatial interaction kernel
andrewkern Dec 17, 2025
e1eefd1
Add SIMD acceleration for Cauchy spatial interaction kernel
andrewkern Dec 17, 2025
7e1c3f0
Add SIMD acceleration for Linear spatial interaction kernel
andrewkern Dec 17, 2025
e4962d7
Replace benchmark with comprehensive all-kernels version
andrewkern Dec 17, 2025
5baabe7
Update README with spatial interaction kernel benchmark documentation
andrewkern Dec 17, 2025
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
67 changes: 31 additions & 36 deletions core/interaction_type.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@
#include "species.h"
#include "slim_globals.h"
#include "sparse_vector.h"
#include "eidos_simd.h"

#include <utility>
#include <algorithm>
Expand Down Expand Up @@ -2846,12 +2847,26 @@ void InteractionType::FillSparseVectorForReceiverStrengths(SparseVector *sv, Ind
// Figure out what index in the exerter subpopulation, if any, needs to be excluded so self-interaction is zero
slim_popsize_t excluded_index = (exerter_subpop == receiver->subpopulation_) ? receiver->index_ : -1;

// We special-case some builds directly to strength values here, for efficiency,
// with no callbacks and spatiality "xy".
// We special-case Fixed kernel builds directly to strength values here, for efficiency,
// with no callbacks and spatiality "xy". Other kernels use the two-pass path below
// which enables SIMD optimizations for Exponential and Normal kernels.
// ADK 12/16/2025: changed to only use special-case path for Fixed kernel
if ((interaction_callbacks.size() == 0) && (spatiality_ == 2) && (if_type_ == SpatialKernelType::kFixed))
{
sv->SetDataType(SparseVectorDataType::kStrengths);
BuildSV_Strengths_f_2(kd_root, receiver_position, excluded_index, sv, 0);
sv->Finished();
return;
}

// ADK 12/16/2025: The original switch below handled all kernel types in the special-case
// single-pass path. Now only Fixed uses the special path above; other kernels fall
// through to the two-pass distance-then-transform path, enabling SIMD optimizations.
#if 0
if ((interaction_callbacks.size() == 0) && (spatiality_ == 2))
{
sv->SetDataType(SparseVectorDataType::kStrengths);

switch (if_type_)
{
case SpatialKernelType::kFixed: BuildSV_Strengths_f_2(kd_root, receiver_position, excluded_index, sv, 0); break;
Expand All @@ -2863,11 +2878,12 @@ void InteractionType::FillSparseVectorForReceiverStrengths(SparseVector *sv, Ind
default:
EIDOS_TERMINATION << "ERROR (InteractionType::FillSparseVectorForReceiverStrengths): (internal error) unoptimized SpatialKernelType value." << EidosTerminate();
}

sv->Finished();
return;
}

#endif

// Set up to build distances first; this is an internal implementation detail, so we require the sparse vector set up for strengths above
sv->SetDataType(SparseVectorDataType::kDistances);

Expand Down Expand Up @@ -2901,53 +2917,32 @@ void InteractionType::FillSparseVectorForReceiverStrengths(SparseVector *sv, Ind
}
case SpatialKernelType::kLinear:
{
for (uint32_t col_iter = 0; col_iter < nnz; ++col_iter)
{
sv_value_t distance = values[col_iter];

values[col_iter] = (sv_value_t)(if_param1_ * (1.0 - distance / max_distance_));
}
// Use SIMD-optimized kernel: fmax = if_param1_, max_distance = max_distance_
Eidos_SIMD::linear_kernel_float32(values, nnz, (float)if_param1_, (float)max_distance_);
break;
}
case SpatialKernelType::kExponential:
{
for (uint32_t col_iter = 0; col_iter < nnz; ++col_iter)
{
sv_value_t distance = values[col_iter];

values[col_iter] = (sv_value_t)(if_param1_ * exp(-if_param2_ * distance));
}
// SIMD-accelerated exponential kernel: strength = fmax * exp(-lambda * distance)
Eidos_SIMD::exp_kernel_float32(values, nnz, (float)if_param1_, (float)if_param2_);
break;
}
case SpatialKernelType::kNormal:
{
for (uint32_t col_iter = 0; col_iter < nnz; ++col_iter)
{
sv_value_t distance = values[col_iter];

values[col_iter] = (sv_value_t)(if_param1_ * exp(-(distance * distance) / n_2param2sq_));
}
// SIMD-accelerated Gaussian kernel: strength = fmax * exp(-d^2 / 2sigma^2)
Eidos_SIMD::normal_kernel_float32(values, nnz, (float)if_param1_, (float)n_2param2sq_);
break;
}
case SpatialKernelType::kCauchy:
{
for (uint32_t col_iter = 0; col_iter < nnz; ++col_iter)
{
sv_value_t distance = values[col_iter];
double temp = distance / if_param2_;

values[col_iter] = (sv_value_t)(if_param1_ / (1.0 + temp * temp));
}
// Use SIMD-optimized kernel: fmax = if_param1_, lambda = if_param2_
Eidos_SIMD::cauchy_kernel_float32(values, nnz, (float)if_param1_, (float)if_param2_);
break;
}
case SpatialKernelType::kStudentsT:
{
for (uint32_t col_iter = 0; col_iter < nnz; ++col_iter)
{
sv_value_t distance = values[col_iter];

values[col_iter] = (sv_value_t)SpatialKernel::tdist(distance, if_param1_, if_param2_, if_param3_);
}
// Use SIMD-optimized kernel: fmax = if_param1_, nu = if_param2_, tau = if_param3_
Eidos_SIMD::tdist_kernel_float32(values, nnz, (float)if_param1_, (float)if_param2_, (float)if_param3_);
break;
}
default:
Expand Down
Loading
Loading