For Development
Building wheels
Linux
Start by installing the build dependencies, if not already installed (check Installation section):
Then create a virtual environment and install the build tools:
And if you want to be explicit:
Finally build the wheel:
This will create the wheel file in the dist/ folder.
MacOS
For MacOS we recommend using homebrew and uv. Start by installing homebrew if not already installed (see Installation section). Then install the c++ toolchain and buiild tools if not already installed:
Then install uv if not already installed:
Then build the wheel with uv:
This will create the wheel file in the dist/ folder.
Windows
Start by creating a virtual environment and upgrading pip:
Then install the build tools:
And if you want to be explicit:
Finally build the wheel:
This will create the wheel file in the dist/ folder.
Debugging Python and C++ code
VSCode
To debug both Python and C++ code in VSCode you need to install the Python extension and the C++ extension by Microsoft. After you create your python environment, you also need to set up a launch configuration. Below is an example launch.json configuration that you can use as a starting point:
{
"version": "0.2.0",
"configurations": [
{
"name": "Python: Launch (frontend)",
"type": "debugpy",
"request": "launch",
"program": "${file}", // your current file
"console": "integratedTerminal",
"justMyCode": true
},
{
"name": "C++: Attach to python (lldb)",
"type": "cppdbg",
"request": "attach",
"MIMode": "lldb",
"program": "${command:python.interpreterPath}", // selected env's python
"processId": "${command:pickProcess}"
}
]
}
This will allow for a two-step debugging process that is the most portable and works on most platforms. To install SERVAL in debug mode you need to satisfy all the requirements listed in the Installation section depending on your platfrom. Then run:
or if you are using uv:
Once the code is installed the Debugging process is as follows:
- Start by selecting the correct interpreter that matches your virtual environment:
Ctrl+Shift+P -> Python: Select Interpreter -> select your env
-
Create breakpoints in both Python and C++ code as needed. It is recommended to set breakpoints in the binding functions in the C++ code. If a breakpoint is not set in the C++ code the debugger will not access the C++ source.
-
Start the Python debugger by selecting the "Python: Launch (frontend)" configuration and pressing F5. This will start the Python code execution and stop at any Python breakpoints.
-
When a breakpoint is hit in Python, select from the RUN AND DEBUG pane the "C++: Attach to python (lldb)" configuration and hit F5. This will try to attach the C++ debugger by asking you to select the correct running process. You can search for "python" in the list of processes to find the correct one. It should be the newest one (last one). If uncertain you can print the process id in the Python code with:
-
Wait to see the C++ debugger attach in your CALL STACK pane. If successful you should be able to hit C++ breakpoints now. If your code stops but the VSCode interface does not update, go to the CALL STACK pane and click on the breakpoint on the C++ debugger. This should update the interface and show you the current line of code.
-
In case all works out but your C++ breakpoints are not hit include
SIGTRAPsomewhere in the C++ function:
#include <csignal>
void some_function() {
// some code
std::raise(SIGTRAP); // this will trigger a breakpoint if the debugger is attached
// some more code
}
This will force the debugger to stop at that line if it is attached. If you cannot move forward after this point. Remove the SIGTRAP line and try to set a breakpoint on the next line. This should allow you to continue debugging without the SIGTRAP line.'
Debugging unit tests
If you would like to do debugging during unit tests, from the VScode test explorer, then add the following configuration to your launch.json:
{
"name": "Python: Debug Tests",
"type": "debugpy",
"request": "launch",
"program": "${file}",
"purpose": [
"debug-test"
],
"console": "integratedTerminal",
"justMyCode": true
},
The remainder of the steps are the same as above. Just make sure to start the "Python: Debug Tests" configuration first.
NOTE:
1) The selected Debug flags in CMakeLists.txt can change numerical results for very tight comparison
(1E-9). If you would still want some Debug information with optimization closer to release use the RelWithDebInfo flag instead of Debug.
Configuring the target CPU architecture
The SERVAL_MARCH CMake variable controls the CPU architecture targeted by the C++ compiler through the -march option (GCC/Clang). This determines which processor-specific instructions and optimizations the compiler can use.
By default, SERVAL_MARCH is set to native, meaning that the compiler optimizes the code for the CPU of the machine performing the compilation.
To specify a different architecture, pass the variable during compilation. For example, to target AMD Zen 4:
On Apple Silicon, a more generic ARM architecture can be requested, if supported by the compiler:
To build a wheel with a specific architecture:
The selected architecture must be supported by the compiler and compatible with the machine on which the compiled package will run. Using native may produce binaries that are not compatible with older processors.
Clean rebuild
When changing the target architecture, a clean rebuild is recommended, particularly when comparing performance.
On Linux or macOS:
rm -rf build/
pip install . --no-cache-dir --force-reinstall \
--config-settings=cmake.define.SERVAL_MARCH=znver4
To rebuild a wheel:
To restore the default behavior, remove the build directory and rebuild without specifying SERVAL_MARCH.
Note: Changing -march does not necessarily produce a significant performance difference. Its effect depends on the CPU, compiler, and workload.
Streaming Gaunt dot12: mathematics, implementation, and optimization
Mathematical formulation
The Gaunt contraction combines two sets of spherical-harmonic coefficients through the Gaunt integral, which describes the coupling of three spherical harmonics.
For the dot12 contraction, the output coefficients are
where \(b\) denotes the batch channel and \(m_2=-m_1-m_3\).
The Gaunt coefficient is expressed in terms of Wigner-3j symbols:
Not every combination of angular indices contributes to the contraction. The Wigner-3j symbols impose several selection rules:
- Magnetic quantum numbers: \(m_1+m_2+m_3=0\).
- Triangle condition: \(|l_2-l_3|\leq l_1\leq l_2+l_3\).
- Magnetic index constraint: \(l_1\geq |m_1|\).
- Parity condition: \(l_1+l_2+l_3\) must be even.
Consequently, for fixed \(l_2\) and \(l_3\), the sum over \(l_1\) only includes values satisfying these constraints. In particular, the parity condition allows the \(l_1\) index to advance in steps of two.
The contraction can therefore be rewritten as
This form separates the calculation into an inner weighted sum over \(l_1\), followed by multiplication by the corresponding \(a^{(2)}\) coefficient and accumulation over \(l_2\).
This separation is central to the optimized implementation: the inner sum can be evaluated simultaneously for a small group of batch channels, while the Gaunt coefficients are shared across those channels.
Direct implementation and performance limitations
The Gaunt contraction is evaluated through nested loops over \(l_3\), \(l_2\), and \(l_1\), with OpenMP parallelization over \((m_1,m_3)\).
For each \(l_2\), the implementation computes
then updates the output:
A simplified representation is:
#pragma omp for schedule(dynamic)
for (size_t idx = 0; idx < m1xm3_len; ++idx) {
for (int l3 = abs_m3; l3 <= l3max; ++l3) {
for (int l2 = abs_m2; l2 <= l2max; ++l2) {
// Compute valid l1 range and Gaunt coefficients.
std::fill(inner_sum.begin(), inner_sum.end(), 0.0);
for (int l1 = l1start; l1 <= l1end; l1 += 2) {
for (size_t b = 0; b < batch_size; ++b) {
inner_sum[b] += gaunt * alm1(im1, l1, b);
}
}
for (size_t b = 0; b < batch_size; ++b) {
result(lim1, l3, im3, b) +=
alm2(im2, l2, b) * inner_sum[b];
}
}
}
}
Performance limitations
The direct implementation has two main bottlenecks:
-
Large accumulation buffer:
inner_sumcontains one complex number per channel. For large batches, it cannot remain entirely in CPU registers, requiring repeated memory accesses. -
Limited cache reuse: Neighboring \(l_2\) values access overlapping
alm1coefficients, but processing all channels before advancing to the next \(l_2\) increases the working set and can cause repeated L2 cache misses.
The optimized implementation addresses both problems through channel blocking and improved data locality, without changing the mathematical contraction.
Channel blocking and register-based accumulation
The first optimization divides the batch dimension into blocks of \(B\) channels, instead of processing all channels simultaneously.
For each block, the inner contraction becomes
The corresponding C++ computation is:
double accumulators[2 * kDot12ChannelBlockSize] = {};
for (int j = 0; j < l1_count; ++j) {
const double* packed_row =
packed_alm1 + j * 2 * kDot12ChannelBlockSize;
for (size_t component = 0;
component < 2 * kDot12ChannelBlockSize;
++component) {
accumulators[component] +=
gaunt_coefficients[j] * packed_row[component];
}
}
Each complex channel uses two real accumulators, one for the real part and one for the imaginary part.
Why is this faster?
- Register accumulation: A small, fixed-size accumulator array gives the compiler a better opportunity to keep intermediate sums in CPU registers.
- SIMD vectorization: The contiguous real-valued data and fixed loop bounds facilitate vectorized instructions.
- Reduced memory traffic: Intermediate sums need not be repeatedly loaded from and stored to a large temporary buffer.
The block size is controlled by kDot12ChannelBlockSize, with a default of 16 channels.
The final block may contain fewer than \(B\) valid channels; the implementation handles this through valid_channels.
The main benefit is that the innermost contraction operates on a small, predictable amount of data rather than the entire batch.
Memory layout and parity-based packing
To support channel blocking, alm1 is reorganized into a contiguous buffer with the layout
For each channel block, the coefficients are separated according to the parity of \(l_1\):
This follows directly from the Gaunt selection rule
For fixed \(l_2\) and \(l_3\), only one parity of \(l_1\) contributes.
The packing operation is implemented in pack_alm1_block():
for (int l1 = 0; l1 <= l1max; ++l1) {
double* row = packed_block
+ (l1 & 1) * parity_stride_doubles
+ (l1 / 2) * doubles_per_row;
std::copy_n(
reinterpret_cast<const double*>(&alm1(im1, l1, b0)),
2 * valid_channels,
row);
std::fill(
row + 2 * valid_channels,
row + doubles_per_row,
0.0);
}
Why is this faster?
- Contiguous access: Consecutive valid \(l_1\) values are stored in adjacent rows, eliminating the stride between same-parity coefficients.
- Efficient SIMD processing: Real and imaginary components are stored contiguously, allowing the inner contraction to operate on a regular array of doubles.
- Simpler indexing: Each contraction accesses consecutive packed rows without repeatedly calculating the original array indices.
The final channel block is zero-padded to maintain a fixed row size.
The main benefit is that the data layout matches the access pattern imposed by the mathematical selection rules.
Loop reordering and cache reuse
The next optimization changes the loop order to process all \(l_2\) contractions for one channel block before moving to the next block.
The direct implementation processes the full batch for each \(l_2\):
for (int l2 = abs_m2; l2 <= l2max; ++l2) {
for (size_t block = 0; block < nblocks; ++block) {
// Contract this channel block.
}
}
The optimized implementation reverses these loops:
for (size_t block = 0; block < nblocks; ++block) {
const double* packed_alm1_block = /* packed block */;
for (const auto& contraction : l2_contractions) {
contract_dot12_block(
gaunt_coefficients.data() + contraction.gaunt_offset,
contraction.l1_count,
packed_alm1_block + contraction.packed_alm1_offset,
&alm2(im2, contraction.l2, b0),
valid_channels,
res_ptr + b0);
}
}
The Gaunt coefficients and valid \(l_1\) ranges are computed beforehand and stored in gaunt_coefficients and l2_contractions.
Why is this faster?
- Better cache reuse: Neighboring \(l_2\) contractions access overlapping
alm1coefficients. Processing one channel block across all \(l_2\) values increases the likelihood that these coefficients remain in L2 cache. - Smaller working set: Only the packed data for one channel block needs to be accessed repeatedly, rather than cycling through the entire batch.
- No repeated coefficient generation: Gaunt coefficients are computed once for each \((m_1,m_3,l_3)\) combination and reused across channel blocks.
The mathematical summation is unchanged; only the order in which its contributions are evaluated differs.
The main benefit is improved temporal locality: the same packed alm1 data is reused while it is still likely to be in cache.
Overall implementation and performance benefits
The optimized gaunt_opt_batch_dot12() combines the previous techniques into a single computational strategy:
// Pack alm1 into channel blocks.
#pragma omp for schedule(static)
for (size_t i1 = 0; i1 < m1len; ++i1) {
for (size_t block = 0; block < nblocks; ++block) {
pack_alm1_block(...);
}
}
// Perform the Gaunt contraction.
#pragma omp for schedule(dynamic)
for (size_t idx = 0; idx < m1xm3_len; ++idx) {
for (int l3 = abs_m3; l3 <= l3max; ++l3) {
build_dot12_coefficients(...);
for (size_t block = 0; block < nblocks; ++block) {
for (const auto& contraction : l2_contractions) {
contract_dot12_block(...);
}
}
}
}
OpenMP scheduling: schedule(static) is used for packing because the workload is approximately uniform across m1 indices, minimizing scheduling overhead. In contrast, schedule(dynamic) is used for the Gaunt contraction because different \((m_1,m_3)\) combinations require different amounts of computation, allowing threads to share the workload more evenly.
The implementation separates data preparation, coefficient generation, and numerical contraction. OpenMP distributes independent \((m_1,m_3)\) combinations across threads.
Summary of performance improvements
| Optimization | Main benefit |
|---|---|
| Channel blocking | Keeps small accumulation buffers in CPU registers |
| Parity-based packing | Provides contiguous access to valid \(l_1\) coefficients |
| Loop reordering | Improves reuse of alm1 data in CPU caches |
| Coefficient reuse | Avoids regenerating Gaunt coefficients for each channel block |
| OpenMP parallelization | Distributes independent contractions across CPU cores |
Together, these changes reduce memory traffic and improve cache utilization and SIMD efficiency.
The mathematical contraction and its selection rules remain unchanged. The optimization primarily improves how data is accessed and reused, rather than reducing the number of terms being evaluated.
Configuring the Gaunt contraction channel block size
The SERVAL_CHANNEL_BLOCK_SIZE CMake variable controls the number of batch channels processed together by the optimized gaunt_opt_batch_dot12 contraction. It determines the channel block size used for packing input data and performing the contraction.
The default value is 16, but it can be changed at compilation time to tune performance for different CPU architectures.
For example, to compile with a block size of 32:
Alternatively, when building a wheel:
If the variable is not specified, the default value of 16 is used. The block size must be a positive integer.
This is a compile-time setting, not a runtime parameter. Different values may affect performance depending on the CPU architecture, cache behavior, and batch size.
Clean rebuild
When changing the block size, a clean rebuild is recommended, especially for performance benchmarking. From the project root, remove the existing CMake build directory before compiling again.
On Linux or macOS:
rm -rf build/
pip install . --no-cache-dir --force-reinstall \
--config-settings=cmake.define.SERVAL_CHANNEL_BLOCK_SIZE=32
To rebuild a wheel instead:
To return to the default block size, remove the build directory and rebuild without specifying SERVAL_CHANNEL_BLOCK_SIZE.