1use crate::cubic_cell_kernel::{DenestedCubicCell, DenestedPartitionCell, LocalSpanCubic};
43use gam_gpu::gpu_error::GpuError;
44
45pub mod kernel_src {
48 pub const DENESTED_PARTITION_CELLS_KERNEL_SRC: &str = r#"
56// f64 throughout (no --use_fast_math).
57
58extern "C" {
59
60__device__ __forceinline__ double pos_inf_f64() {
61 // IEEE-754 +inf bit pattern: 0x7ff0000000000000.
62 return __longlong_as_double((long long)0x7ff0000000000000LL);
63}
64__device__ __forceinline__ double neg_inf_f64() {
65 // IEEE-754 -inf bit pattern: 0xfff0000000000000.
66 return __longlong_as_double((long long)0xfff0000000000000LL);
67}
68
69__global__ void denested_partition_cells_kernel(
70 int n_rows,
71 double scale,
72 const double *a_per_row,
73 const double *b_per_row,
74 double *out_cells_flat, // 18 doubles per row (single cell)
75 unsigned int *out_row_offsets, // length n_rows + 1
76 unsigned char *out_status // length n_rows
77) {
78 int i = blockIdx.x * blockDim.x + threadIdx.x;
79 if (i >= n_rows) return;
80 double a = a_per_row[i];
81 double b = b_per_row[i];
82 double *cell = out_cells_flat + (long long)i * 18;
83 // ── cell: (-inf, +inf, c0=a*scale, c1=b*scale, c2=0, c3=0) ──
84 cell[0] = neg_inf_f64();
85 cell[1] = pos_inf_f64();
86 cell[2] = a * scale;
87 cell[3] = b * scale;
88 cell[4] = 0.0;
89 cell[5] = 0.0;
90 // ── score_span (zero cubic, left=0,right=1) ──
91 cell[6] = 0.0; cell[7] = 1.0;
92 cell[8] = 0.0; cell[9] = 0.0; cell[10] = 0.0; cell[11] = 0.0;
93 // ── link_span (zero cubic, left=0,right=1) ──
94 cell[12] = 0.0; cell[13] = 1.0;
95 cell[14] = 0.0; cell[15] = 0.0; cell[16] = 0.0; cell[17] = 0.0;
96 // ── row offset: one cell per row ──
97 out_row_offsets[i] = (unsigned int)i;
98 if (i == n_rows - 1) {
99 out_row_offsets[n_rows] = (unsigned int)n_rows;
100 }
101 out_status[i] = 0;
102}
103
104} // extern "C"
105"#;
106
107 pub const DENESTED_CELL_PRIMARY_FIXED_PARTIALS_KERNEL_SRC: &str = r#"
125// f64 throughout (no --use_fast_math).
126
127extern "C" {
128
129__global__ void denested_cell_primary_fixed_partials_kernel(
130 int n_cells_total,
131 unsigned int r,
132 unsigned int g_slot,
133 double scale,
134 double *out_partials_flat, // (12 + 40·r) doubles per cell
135 unsigned char *out_status
136) {
137 int cell = blockIdx.x * blockDim.x + threadIdx.x;
138 if (cell >= n_cells_total) return;
139 unsigned int per_cell = 12u + 40u * r;
140 double *base = out_partials_flat + (long long)cell * (long long)per_cell;
141 // Zero the whole block (cheap; r is small).
142 for (unsigned int s = 0; s < per_cell; ++s) {
143 base[s] = 0.0;
144 }
145 // dc_da = [1, 0, 0, 0] · scale
146 base[0] = scale;
147 // dc_daa, dc_daaa already zero.
148 // g-slot fills (offset = 12 + 4·g_slot within each per-cell run).
149 // coeff_u [g] = dc_db = [0, 1, 0, 0] · scale
150 // coeff_au [g] = dc_dab = [0, 0, 0, 0]
151 // coeff_bu [g] = dc_dbb = [0, 0, 0, 0]
152 // coeff_aau [g] = dc_daab = [0, 0, 0, 0]
153 // coeff_abu [g] = dc_dabb = [0, 0, 0, 0]
154 // coeff_bbu [g] = dc_dbbb = [0, 0, 0, 0]
155 // (third partials all zero in the no-runtime case)
156 unsigned int g_off = 12u + 4u * g_slot;
157 base[g_off + 1] = scale; // coeff_u[g][1] = scale
158 out_status[cell] = 0;
159}
160
161} // extern "C"
162"#;
163}
164
165#[derive(Clone, Copy, Debug)]
167pub struct PartitionCellsRowInputs<'a> {
168 pub a: f64,
169 pub b: f64,
170 pub beta_h: Option<&'a [f64]>,
171 pub beta_w: Option<&'a [f64]>,
172}
173
174pub type PartitionCellsOutput = Vec<Vec<DenestedPartitionCell>>;
177
178pub fn try_device_partition_cells(
189 rows: &[PartitionCellsRowInputs<'_>],
190) -> Result<Option<PartitionCellsOutput>, GpuError> {
191 if rows.is_empty() {
192 return Ok(Some(Vec::new()));
193 }
194 let trivial = rows
198 .iter()
199 .all(|r| r.beta_h.is_none() && r.beta_w.is_none());
200 if !trivial {
201 return Ok(None);
202 }
203 device_dispatch::partition_cells_baseline(rows, 1.0)
204}
205
206#[derive(Clone, Copy, Debug)]
208pub struct CellPrimaryFixedPartialsCellInputs {
209 pub score_span: LocalSpanCubic,
210 pub link_span: LocalSpanCubic,
211}
212
213#[derive(Clone, Copy, Debug)]
218pub struct CellPrimaryFixedPartialsRowInputs<'a> {
219 pub cells: &'a [CellPrimaryFixedPartialsCellInputs],
220 pub layout: FlexPrimaryLayout,
221}
222
223#[derive(Clone, Debug, Default)]
229pub struct CellPrimaryFixedPartialsOutput {
230 pub partials: Vec<Vec<Vec<f64>>>,
231}
232
233#[derive(Clone, Copy, Debug)]
239pub struct FlexPrimaryLayout {
240 pub r: u32,
241 pub g_slot: u32,
242}
243
244pub fn try_device_cell_primary_fixed_partials(
255 rows: &[CellPrimaryFixedPartialsRowInputs<'_>],
256) -> Result<Option<CellPrimaryFixedPartialsOutput>, GpuError> {
257 if rows.is_empty() {
258 return Ok(Some(CellPrimaryFixedPartialsOutput::default()));
259 }
260 let trivial_spans = rows.iter().all(|row| {
264 row.cells
265 .iter()
266 .all(|cell| span_is_zero(cell.score_span) && span_is_zero(cell.link_span))
267 });
268 if !trivial_spans {
269 return Ok(None);
270 }
271 let layout0 = rows[0].layout;
275 if !rows
276 .iter()
277 .all(|r| r.layout.r == layout0.r && r.layout.g_slot == layout0.g_slot)
278 {
279 return Ok(None);
280 }
281 let mut row_cell_counts: Vec<usize> = rows.iter().map(|r| r.cells.len()).collect();
284 let total_cells: usize = row_cell_counts.iter().copied().sum();
285 if total_cells == 0 {
286 let mut partials: Vec<Vec<Vec<f64>>> = Vec::with_capacity(rows.len());
287 for _ in 0..rows.len() {
288 partials.push(Vec::new());
289 }
290 return Ok(Some(CellPrimaryFixedPartialsOutput { partials }));
291 }
292 let flat = match device_dispatch::cell_primary_fixed_partials_baseline(layout0, total_cells) {
293 Ok(flat) => flat,
294 Err(_) => return Ok(None),
295 };
296 let per_cell = 12usize + 40usize * (layout0.r as usize);
297 let mut partials: Vec<Vec<Vec<f64>>> = Vec::with_capacity(rows.len());
298 let mut cursor = 0usize;
299 for n_cells in row_cell_counts.drain(..) {
300 let mut row_cells: Vec<Vec<f64>> = Vec::with_capacity(n_cells);
301 for _ in 0..n_cells {
302 row_cells.push(flat[cursor..cursor + per_cell].to_vec());
303 cursor += per_cell;
304 }
305 partials.push(row_cells);
306 }
307 assert_eq!(cursor, flat.len());
308 Ok(Some(CellPrimaryFixedPartialsOutput { partials }))
309}
310
311#[inline]
312fn span_is_zero(span: LocalSpanCubic) -> bool {
313 span.c0 == 0.0 && span.c1 == 0.0 && span.c2 == 0.0 && span.c3 == 0.0
314}
315
316pub fn trivial_partition_cell(a: f64, b: f64, scale: f64) -> DenestedPartitionCell {
320 DenestedPartitionCell {
321 cell: DenestedCubicCell {
322 left: f64::NEG_INFINITY,
323 right: f64::INFINITY,
324 c0: a * scale,
325 c1: b * scale,
326 c2: 0.0,
327 c3: 0.0,
328 },
329 score_span: LocalSpanCubic {
330 left: 0.0,
331 right: 1.0,
332 c0: 0.0,
333 c1: 0.0,
334 c2: 0.0,
335 c3: 0.0,
336 },
337 link_span: LocalSpanCubic {
338 left: 0.0,
339 right: 1.0,
340 c0: 0.0,
341 c1: 0.0,
342 c2: 0.0,
343 c3: 0.0,
344 },
345 left_edge: crate::cubic_cell_kernel::PartitionEdge::Fixed(f64::NEG_INFINITY),
346 right_edge: crate::cubic_cell_kernel::PartitionEdge::Fixed(f64::INFINITY),
347 }
348}
349
350#[cfg(target_os = "linux")]
351mod device_dispatch {
352 use super::kernel_src::DENESTED_PARTITION_CELLS_KERNEL_SRC;
353 use super::{PartitionCellsOutput, PartitionCellsRowInputs, trivial_partition_cell};
354 use cudarc::driver::{LaunchConfig, PushKernelArg};
355 use gam_gpu::device_cache::PtxModuleCache;
356 use gam_gpu::gpu_err as gam_gpu_err;
357 use gam_gpu::gpu_error::{GpuError, GpuResultExt};
358 use gam_gpu::solver::context_and_stream;
359
360 static PARTITION_PTX_CACHE: PtxModuleCache = PtxModuleCache::new();
361
362 const THREADS_PER_BLOCK: u32 = 128;
363
364 pub(super) fn partition_cells_baseline(
366 rows: &[PartitionCellsRowInputs<'_>],
367 scale: f64,
368 ) -> Result<Option<PartitionCellsOutput>, GpuError> {
369 let n = rows.len();
370 let n_u32 = u32::try_from(n)
371 .map_err(|_| gam_gpu_err!("partition_cells_baseline: n_rows={n} exceeds u32"))?;
372 let n_i32 = i32::try_from(n)
373 .map_err(|_| gam_gpu_err!("partition_cells_baseline: n_rows={n} exceeds i32"))?;
374 let (ctx, stream) = match context_and_stream() {
375 Ok(pair) => pair,
376 Err(_) => return Ok(None),
377 };
378 let module = PARTITION_PTX_CACHE.get_or_compile(
379 &ctx,
380 "survival_flex_prep::partition_cells",
381 DENESTED_PARTITION_CELLS_KERNEL_SRC,
382 )?;
383 let func = module
384 .load_function("denested_partition_cells_kernel")
385 .gpu_ctx("survival_flex_prep: load_function partition_cells")?;
386
387 let a_host: Vec<f64> = rows.iter().map(|r| r.a).collect();
388 let b_host: Vec<f64> = rows.iter().map(|r| r.b).collect();
389 let a_dev = stream
390 .clone_htod(&a_host)
391 .gpu_ctx("survival_flex_prep: upload a_per_row")?;
392 let b_dev = stream
393 .clone_htod(&b_host)
394 .gpu_ctx("survival_flex_prep: upload b_per_row")?;
395 let mut cells_dev = stream
396 .alloc_zeros::<f64>(n * 18)
397 .gpu_ctx("survival_flex_prep: alloc cells_flat")?;
398 let mut offsets_dev = stream
399 .alloc_zeros::<u32>(n + 1)
400 .gpu_ctx("survival_flex_prep: alloc row_offsets")?;
401 let mut status_dev = stream
402 .alloc_zeros::<u8>(n)
403 .gpu_ctx("survival_flex_prep: alloc status")?;
404
405 let cfg = LaunchConfig {
406 grid_dim: (n_u32.div_ceil(THREADS_PER_BLOCK).max(1), 1, 1),
407 block_dim: (THREADS_PER_BLOCK, 1, 1),
408 shared_mem_bytes: 0,
409 };
410 unsafe {
415 let mut builder = stream.launch_builder(&func);
416 builder.arg(&n_i32);
417 builder.arg(&scale);
418 builder.arg(&a_dev);
419 builder.arg(&b_dev);
420 builder.arg(&mut cells_dev);
421 builder.arg(&mut offsets_dev);
422 builder.arg(&mut status_dev);
423 builder.launch(cfg)
424 }
425 .map(|_event_pair| ())
426 .gpu_ctx("survival_flex_prep: launch partition_cells")?;
427
428 let cells_host = stream
429 .clone_dtoh(&cells_dev)
430 .gpu_ctx("survival_flex_prep: download cells_flat")?;
431 let status_host = stream
432 .clone_dtoh(&status_dev)
433 .gpu_ctx("survival_flex_prep: download status")?;
434 for (i, st) in status_host.iter().enumerate() {
435 if *st != 0 {
436 return Err(gam_gpu_err!(
437 "survival_flex_prep: row {i} status={st} from device kernel"
438 ));
439 }
440 }
441 assert_eq!(cells_host.len(), n * 18);
442 let mut out: PartitionCellsOutput = Vec::with_capacity(n);
449 for i in 0..n {
450 let base = i * 18;
451 let c0 = cells_host[base + 2];
452 let c1 = cells_host[base + 3];
453 let mut cell = trivial_partition_cell(rows[i].a, rows[i].b, scale);
454 cell.cell.c0 = c0;
457 cell.cell.c1 = c1;
458 out.push(vec![cell]);
459 }
460 Ok(Some(out))
461 }
462
463 pub(super) fn cell_primary_fixed_partials_baseline(
472 layout: super::FlexPrimaryLayout,
473 n_cells_total: usize,
474 ) -> Result<Vec<f64>, GpuError> {
475 use super::kernel_src::DENESTED_CELL_PRIMARY_FIXED_PARTIALS_KERNEL_SRC;
476 static FP_PTX_CACHE: PtxModuleCache = PtxModuleCache::new();
477
478 let n_i32 = i32::try_from(n_cells_total).map_err(|_| {
479 gam_gpu_err!(
480 "cell_primary_fixed_partials_baseline: n_cells={n_cells_total} exceeds i32"
481 )
482 })?;
483 let n_u32 = u32::try_from(n_cells_total).map_err(|_| {
484 gam_gpu_err!(
485 "cell_primary_fixed_partials_baseline: n_cells={n_cells_total} exceeds u32"
486 )
487 })?;
488 let (ctx, stream) = context_and_stream()
489 .map_err(|reason| gam_gpu::gpu_error::GpuError::DriverCallFailed { reason })?;
490 let module = FP_PTX_CACHE.get_or_compile(
491 &ctx,
492 "survival_flex_prep::cell_primary_fixed_partials",
493 DENESTED_CELL_PRIMARY_FIXED_PARTIALS_KERNEL_SRC,
494 )?;
495 let func = module
496 .load_function("denested_cell_primary_fixed_partials_kernel")
497 .gpu_ctx("survival_flex_prep: load_function fixed_partials")?;
498
499 let per_cell = 12usize + 40usize * (layout.r as usize);
500 let scale = 1.0f64;
501 let mut out_dev = stream
502 .alloc_zeros::<f64>(n_cells_total * per_cell)
503 .gpu_ctx("survival_flex_prep: alloc fixed_partials")?;
504 let mut status_dev = stream
505 .alloc_zeros::<u8>(n_cells_total)
506 .gpu_ctx("survival_flex_prep: alloc fixed_partials status")?;
507 let cfg = LaunchConfig {
508 grid_dim: (n_u32.div_ceil(THREADS_PER_BLOCK).max(1), 1, 1),
509 block_dim: (THREADS_PER_BLOCK, 1, 1),
510 shared_mem_bytes: 0,
511 };
512 unsafe {
515 let mut builder = stream.launch_builder(&func);
516 builder.arg(&n_i32);
517 builder.arg(&layout.r);
518 builder.arg(&layout.g_slot);
519 builder.arg(&scale);
520 builder.arg(&mut out_dev);
521 builder.arg(&mut status_dev);
522 builder.launch(cfg)
523 }
524 .map(|_event_pair| ())
525 .gpu_ctx("survival_flex_prep: launch fixed_partials")?;
526 let out_host = stream
527 .clone_dtoh(&out_dev)
528 .gpu_ctx("survival_flex_prep: download fixed_partials")?;
529 let status_host = stream
530 .clone_dtoh(&status_dev)
531 .gpu_ctx("survival_flex_prep: download fixed_partials status")?;
532 for (i, st) in status_host.iter().enumerate() {
533 if *st != 0 {
534 return Err(gam_gpu_err!(
535 "survival_flex_prep: fixed_partials cell {i} status={st}"
536 ));
537 }
538 }
539 Ok(out_host)
540 }
541}
542
543#[cfg(not(target_os = "linux"))]
544mod device_dispatch {
545 use super::{PartitionCellsOutput, PartitionCellsRowInputs};
546 use gam_gpu::gpu_err as gam_gpu_err;
547 use gam_gpu::gpu_error::GpuError;
548
549 pub(super) fn partition_cells_baseline(
550 rows: &[PartitionCellsRowInputs<'_>],
551 scale: f64,
552 ) -> Result<Option<PartitionCellsOutput>, GpuError> {
553 let first = rows.first().map(|row| (row.a, row.b));
557 log::trace!(
558 "survival_flex_prep::partition_cells_baseline declined on non-linux \
559 (n_rows={}, scale={scale}, first_ab={first:?})",
560 rows.len(),
561 );
562 Ok(None)
563 }
564
565 pub(super) fn cell_primary_fixed_partials_baseline(
566 layout: super::FlexPrimaryLayout,
567 n_cells_total: usize,
568 ) -> Result<Vec<f64>, GpuError> {
569 Err(gam_gpu_err!(
570 "survival_flex_prep::cell_primary_fixed_partials_baseline: CUDA only supported on linux \
571 (would have launched n_cells={n_cells_total}, r={}, g_slot={})",
572 layout.r,
573 layout.g_slot
574 ))
575 }
576}
577
578#[cfg(test)]
579mod tests {
580 use super::*;
581
582 #[test]
583 fn empty_partition_inputs_short_circuit() {
584 let out = try_device_partition_cells(&[]).expect("ok");
585 assert!(out.is_some());
586 assert!(out.unwrap().is_empty());
587 }
588
589 #[test]
590 fn nonempty_partition_with_betas_declines() {
591 let beta = [0.0_f64];
592 let inputs = [PartitionCellsRowInputs {
593 a: 0.0,
594 b: 1.0,
595 beta_h: Some(&beta),
596 beta_w: None,
597 }];
598 let out = try_device_partition_cells(&inputs).expect("ok");
599 assert!(out.is_none());
602 }
603
604 #[test]
605 fn empty_fixed_partials_inputs_short_circuit() {
606 let out = try_device_cell_primary_fixed_partials(&[]).expect("ok");
607 assert!(out.is_some());
608 assert!(out.unwrap().partials.is_empty());
609 }
610
611 #[test]
612 fn empty_cells_per_row_returns_empty_partials() {
613 let inputs = [CellPrimaryFixedPartialsRowInputs {
614 cells: &[],
615 layout: FlexPrimaryLayout { r: 4, g_slot: 3 },
616 }];
617 let out = try_device_cell_primary_fixed_partials(&inputs).expect("ok");
618 let some = out.expect("Some when all rows have zero cells");
619 assert_eq!(some.partials.len(), 1);
620 assert!(some.partials[0].is_empty());
621 }
622
623 #[test]
624 fn kernel_src_strings_are_nonempty() {
625 assert!(!kernel_src::DENESTED_PARTITION_CELLS_KERNEL_SRC.is_empty());
626 assert!(!kernel_src::DENESTED_CELL_PRIMARY_FIXED_PARTIALS_KERNEL_SRC.is_empty());
627 }
628
629 #[test]
630 fn trivial_partition_cell_matches_cpu_empty_split_branch() {
631 let cell = trivial_partition_cell(2.5, -1.25, 1.0);
635 assert_eq!(cell.cell.c0, 2.5);
636 assert_eq!(cell.cell.c1, -1.25);
637 assert_eq!(cell.cell.c2, 0.0);
638 assert_eq!(cell.cell.c3, 0.0);
639 assert!(cell.cell.left.is_infinite() && cell.cell.left.is_sign_negative());
640 assert!(cell.cell.right.is_infinite() && cell.cell.right.is_sign_positive());
641 }
642}