diff --git a/MIGRATION.md b/MIGRATION.md new file mode 100644 index 0000000..0840873 --- /dev/null +++ b/MIGRATION.md @@ -0,0 +1,133 @@ +# fastfields migration status & plan + +This document tracks the port of the jitfields kernels into the layered, +DLPack-based `fastfields` C++/CUDA libraries. See `.github/profile/README.md` +for the repository hierarchy. + +## Layer recap + +``` +kernels (voxelwise, header-only, templated) + └─ {cpu,cuda}-impl (loops over elements; CPU thread pool / CUDA kernels; templated, dynamic sizes) + └─ {cpu,cuda}-lib (exported symbols; unsafe pointers -> template dispatch on dtype/dim) + └─ lib (exported symbols; DLTensor in; dispatch on device -> cpu/cuda lib) +``` + +Each ported operation appears at three "library" levels: +- `{cpu,cuda}-lib/.{h,cpp}` — `ff::cpu::` / `ff::cuda::` functions taking + `DLTensor&`, dispatching on dtype (and dim/spline/bound where relevant) to the + templated impl. Pointer/stride arrays are narrowed to 32-bit when + `canUse32BitIndexMath` allows (see `autocast.h` / `copy_if_needed`). +- `lib/.{h,cpp}` — public `ff::` functions taking `DLTensor&`, dispatching + on `device.device_type` to the cpu or cuda lib. + +## Status matrix + +| module | kernels | cpu-impl | cuda-impl | cpu-lib | cuda-lib | lib | CPU tested | +|----------------|:------:|:--------:|:---------:|:-------:|:--------:|:---:|:----------:| +| distance | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | **yes** | +| posdef | ✓ | ✓ | ✓ | ✗ | ✗ | ✗ | no | +| pushpull | ✓ | ✓ | ✓ | ✗ | ✗ | ✗ | no | +| reg_field | ✓ | ✓ | ✓ | ✗ | ✗ | ✗ | no | +| reg_flow | ✓ | ✓ | ✓ | ✗ | ✗ | ✗ | no | +| resize | ✓ | ✓ | ✓ | ✗ | ✗ | ✗ | no | +| restrict | ✓ | ✓ | ✓ | ✗ | ✗ | ✗ | no | +| splinc | ✓ | ✓ | ✓ | ✗ | ✗ | ✗ | no | +| tetrahedron | ✓ | ✓ | — | ✗ | ✗ | ✗ | no | + +"✓" for a lib column means the dispatch layer exists; only `distance` is currently +compiled and tested (CPU). `tetrahedron` has no cuda-impl header yet. + +## Build & test (CPU) + +No `nvcc` is assumed available in CI; the CPU path is the primary automated gate. +Submodules must be checked out (or symlinked) so the include nesting resolves: +`cpu-lib/impl -> cpu-impl`, `cpu-impl/kernels -> kernels`, `lib/cpu -> cpu-lib`. + +``` +make -C fastfields-cpu-lib test CXX=clang++ # builds libfastfields-cpu.so + runs tests/test_*.cpp +make -C fastfields-lib all CXX=clang++ # builds libfastfields.so (links cpu lib) +``` + +`make test` compiles each `tests/test_.cpp` together with the module sources +and runs it. New modules should ship a `tests/test_.cpp` that validates the +CPU path against a brute-force / reference implementation, as `test_distance.cpp` does. + +## Bugs found & fixed (CPU, verified) + +1. **kernels/parallel_impl.h** — missing ``/``; would not compile. +2. **kernels/distance/mesh.h** — missing `` (`std::sort`). +3. **cpu-impl/distance_euclidean.h** — `z`/`d` scratch buffers declared `offset_t*` + but assigned `new scalar_t[n]` and passed where `scalar_t*` is expected; failed + to compile when `scalar_t != offset_t`. +4. **lib/distance.cpp** — `dt_l1` dispatched to `dt_euclidean` in both branches; + the L1 transform was never called. +5. **cpu-lib/Makefile, lib/Makefile** — object compile rule missing `-fPIC`, so the + shared-library link failed (`relocation R_X86_64_PC32 … recompile with -fPIC`). + +## Bugs found, NOT yet fixed (need a CUDA build to verify) + +- **cuda-impl/distance_euclidean.h** — `dt()` allocates a scratch buffer of + `vector_size * (sizeof(offset_t) + 2*sizeof(scalar_t))` bytes, but `dt_kernel` + addresses `z`/`d` at `buf + stride_buf * n * …` and indexes each of the + `stride_buf` (= total launched threads) lanes with its own length-`n` `v/z/d` + region. The allocation therefore appears short by a factor of `stride_buf`, + which would be an out-of-bounds device write. Confirm and fix once a CUDA + toolchain is available; add a CUDA test mirroring `test_distance.cpp`. +- **cpu-lib/distance.cpp** `_dt_spline_*` — the `copy_if_needed(…, ndim)` + calls pass the spatial dim `ndim` (1/2/3) as the array length for `size`/stride + arrays whose true length is `nbatch (+1/+2)`. In the 32-bit index path this + under-copies and the kernel may read past the narrowed arrays. Needs a spline + test (CPU) to confirm before changing; euclidean/l1 are unaffected (their arrays + really are length `ndim`). + +## Porting pattern (per module) + +Use `distance.{h,cpp}` at each level as the template. + +1. **cpu-lib/``.h** — declare `ff::cpu::(DLTensor&, …)` for each op. +2. **cpu-lib/``.cpp** — for each op: + - a `CHECK_*` block (shape/dtype/batch), reusing the macros in `distance.cpp` + (consider hoisting these into a shared `checks.h`); + - a `DISPATCH_*` macro selecting `scalar_t`/`offset_t` (+ `ndim`, spline, bound); + - an anonymous-namespace `_<…>` template that narrows pointers via + `copy_if_needed`, casts `void*`, calls the impl, then `free_if_needed`. +3. **cuda-lib/``.{h,cpp}** — same as cpu-lib but `ff::cuda::`, forwarding + `stream`, including `impl/…` from cuda-impl. (Cannot be compiled here.) +4. **lib/``.{h,cpp}** — public `ff::(DLTensor&, …)` dispatching on + `device.device_type` to `FF_CPU::` / `FF_CUDA::` (guard cuda with `FF_WITH_CUDA`). +5. **Makefiles** — add `` to `MODULES` in cpu-lib, cuda-lib and lib. +6. **tests/test_``.cpp** in cpu-lib — CPU correctness vs. reference. + +### Avoiding merge conflicts when porting in parallel +Module `.cpp/.h` and `tests/test_.cpp` are disjoint per module. The only +shared files are the three `MODULES` lists in the Makefiles. When several modules +are ported concurrently, add each module's `.cpp` compiled *directly* with +its test (`clang++ -I. tests/test_.cpp .cpp -o …`) to validate +without editing the Makefile, and record the required `MODULES +=` line in the PR; +integrate the `MODULES` lists in a single follow-up commit. + +## Task breakdown + +- **T1 (done)** distance: CPU green + tested; fixes above; pushed. +- **T2 (this doc)** migration plan & status matrix. +- **T3** posdef: `sym_matvec[_backward]`, `sym_addmatvec_`, `sym_submatvec_`, + `sym_solve[_]`, `sym_invert[_]` over `DLTensor`. Matrix layouts (full/sym/diag/ + estatics/eye) are selected by a runtime enum — dispatch on it. CPU test: compare + `sym_matvec` then `sym_solve` round-trips a known SPD system. +- **T4** resize + restrict + splinc: smaller loops (`loop`/`loopnd`, `loop2`). + resize/restrict take spline order + bound + a scale/anchor; splinc is the spline + prefilter. CPU tests: resize by integer factor vs. manual interpolation; splinc + followed by evaluation reproduces samples. +- **T5** pushpull: `pull/push/count/grad/hess` + `_backward`, dispatched on + dim × spline × bound × dtype. Large; port `pull`/`push` first with a small + gather/scatter CPU test, then the rest. +- **T6** regularisers (reg_field, reg_flow): `matvec/kernel/diag/relax` for + absolute/membrane/bending (+ RLS). Dispatch on dim × dtype; energies parametrised + by `absolute/membrane/bending` weights + voxel size. CPU test: `matvec` equals the + finite-difference Laplacian for the membrane term. +- **T7** de-templating audit + Makefile/CI hardening: confirm which impl entry + points still template runtime sizes that the migration intends to de-template + (cross-check jitfields); wire `USE_OPENMP` (currently defined but unused) into + `CXXFLAGS`; add cuda-lib `-Xcompiler -fPIC`; add a SessionStart/CI job that runs + the CPU build + tests. diff --git a/Makefile b/Makefile index e43392d..1027cc3 100644 --- a/Makefile +++ b/Makefile @@ -107,7 +107,8 @@ clean-cpu: ######################################################################## MODULES = \ - distance + distance \ + posdef OBJECTS = $(addprefix $(BUILDDIR)/,$(addsuffix .$(MOSUF),$(MODULES))) CPPFILES = $(addsuffix .cpp,$(MODULES)) @@ -136,7 +137,7 @@ $(BUILDDIR)/libfastfields.$(SOSUF): $(OBJECTS) $(BUILDDIR)/%.$(MOSUF): %.cpp | $(BUILDDIR) $(CXX) $(CXXFLAGS) $(INCLUDES) \ - -c -o $@ $< + -fPIC -c -o $@ $< ######################################################################## # Messages diff --git a/distance.cpp b/distance.cpp index 54a2bca..ebb05cc 100644 --- a/distance.cpp +++ b/distance.cpp @@ -33,10 +33,10 @@ void dt_l1( { #ifdef FF_WITH_CUDA if (IS_CUDA(inp_out)) - return FF_CUDA::dt_euclidean(inp_out, voxel_spacing, stream); + return FF_CUDA::dt_l1(inp_out, voxel_spacing, stream); #endif if (IS_CPU(inp_out)) - return FF_CPU::dt_euclidean(inp_out, voxel_spacing, stream); + return FF_CPU::dt_l1(inp_out, voxel_spacing, stream); throw std::invalid_argument("unsupported device"); } diff --git a/posdef.cpp b/posdef.cpp new file mode 100644 index 0000000..9c64708 --- /dev/null +++ b/posdef.cpp @@ -0,0 +1,138 @@ +#include +#include "posdef.h" +#include "cpu/posdef.h" +#ifdef FF_WITH_CUDA +#include "cuda/posdef.h" +#endif + +using namespace FF; + +#define IS_CUDA(tensor) (tensor.device.device_type == DLDeviceType::kDLCUDA) +#define IS_CPU(tensor) (tensor.device.device_type == DLDeviceType::kDLCPU || \ + tensor.device.device_type == DLDeviceType::kDLCUDAHost) + +void sym_matvec( + DLTensor & out , + const DLTensor & hessian , + const DLTensor & inp , + int stream ) +{ +#ifdef FF_WITH_CUDA + if (IS_CUDA(out)) + return FF_CUDA::sym_matvec(out, hessian, inp, stream); +#endif + if (IS_CPU(out)) + return FF_CPU::sym_matvec(out, hessian, inp, stream); + + throw std::invalid_argument("unsupported device"); +} + +void sym_matvec_backward( + DLTensor & out , + const DLTensor & grd , + const DLTensor & inp , + int stream ) +{ +#ifdef FF_WITH_CUDA + if (IS_CUDA(out)) + return FF_CUDA::sym_matvec_backward(out, grd, inp, stream); +#endif + if (IS_CPU(out)) + return FF_CPU::sym_matvec_backward(out, grd, inp, stream); + + throw std::invalid_argument("unsupported device"); +} + +void sym_addmatvec_( + DLTensor & out , + const DLTensor & hessian , + const DLTensor & inp , + int stream ) +{ +#ifdef FF_WITH_CUDA + if (IS_CUDA(out)) + return FF_CUDA::sym_addmatvec_(out, hessian, inp, stream); +#endif + if (IS_CPU(out)) + return FF_CPU::sym_addmatvec_(out, hessian, inp, stream); + + throw std::invalid_argument("unsupported device"); +} + +void sym_submatvec_( + DLTensor & out , + const DLTensor & hessian , + const DLTensor & inp , + int stream ) +{ +#ifdef FF_WITH_CUDA + if (IS_CUDA(out)) + return FF_CUDA::sym_submatvec_(out, hessian, inp, stream); +#endif + if (IS_CPU(out)) + return FF_CPU::sym_submatvec_(out, hessian, inp, stream); + + throw std::invalid_argument("unsupported device"); +} + +void sym_solve( + DLTensor & out , + const DLTensor & hessian , + const DLTensor & inp , + const DLTensor & weight , + int stream ) +{ +#ifdef FF_WITH_CUDA + if (IS_CUDA(out)) + return FF_CUDA::sym_solve(out, hessian, inp, weight, stream); +#endif + if (IS_CPU(out)) + return FF_CPU::sym_solve(out, hessian, inp, weight, stream); + + throw std::invalid_argument("unsupported device"); +} + +void sym_solve_( + DLTensor & inp_out , + const DLTensor & hessian , + const DLTensor & weight , + int stream ) +{ +#ifdef FF_WITH_CUDA + if (IS_CUDA(inp_out)) + return FF_CUDA::sym_solve_(inp_out, hessian, weight, stream); +#endif + if (IS_CPU(inp_out)) + return FF_CPU::sym_solve_(inp_out, hessian, weight, stream); + + throw std::invalid_argument("unsupported device"); +} + +void sym_invert( + DLTensor & out , + const DLTensor & hessian , + int stream ) +{ +#ifdef FF_WITH_CUDA + if (IS_CUDA(out)) + return FF_CUDA::sym_invert(out, hessian, stream); +#endif + if (IS_CPU(out)) + return FF_CPU::sym_invert(out, hessian, stream); + + throw std::invalid_argument("unsupported device"); +} + +void sym_invert_( + DLTensor & hessian , + int stream ) +{ +#ifdef FF_WITH_CUDA + if (IS_CUDA(hessian)) + return FF_CUDA::sym_invert_(hessian, stream); +#endif + if (IS_CPU(hessian)) + return FF_CPU::sym_invert_(hessian, stream); + + throw std::invalid_argument("unsupported device"); +} diff --git a/posdef.h b/posdef.h new file mode 100644 index 0000000..b588af2 --- /dev/null +++ b/posdef.h @@ -0,0 +1,83 @@ +#ifndef FF_LIB_POSDEF +#define FF_LIB_POSDEF +#include "dlpack.h" +#include "defines.h" + +FF_NAMESPACE_BEGIN(FF) + +/** + * Compact symmetric ("Sym") positive-definite matrix operations. + * + * The trailing dimension of a vector tensor holds the C channels; the + * trailing dimension of a matrix tensor holds the C*(C+1)/2 unique entries + * of a symmetric CxC matrix in a "diagonal-then-rows" layout: + * [ h00 h11 ... h(C-1)(C-1) | h01 h02 ... h0(C-1) | h12 ... ] + * e.g. C=2 -> [h00, h11, h01]; C=3 -> [h00, h11, h22, h01, h02, h12]. + * Every leading dimension is treated as batch. + */ + +/** out = H * inp */ +void sym_matvec( + DLTensor & out , + const DLTensor & hessian , + const DLTensor & inp , + int stream = 0 +); + +/** Backward of matvec wrt the matrix: out (*batch, C*(C+1)/2) from grd, inp (*batch, C) */ +void sym_matvec_backward( + DLTensor & out , + const DLTensor & grd , + const DLTensor & inp , + int stream = 0 +); + +/** out += H * inp */ +void sym_addmatvec_( + DLTensor & out , + const DLTensor & hessian , + const DLTensor & inp , + int stream = 0 +); + +/** out -= H * inp */ +void sym_submatvec_( + DLTensor & out , + const DLTensor & hessian , + const DLTensor & inp , + int stream = 0 +); + +/** out = (H + diag(weight)) \ inp (weight optional: pass a null-data DLTensor) */ +void sym_solve( + DLTensor & out , + const DLTensor & hessian , + const DLTensor & inp , + const DLTensor & weight , + int stream = 0 +); + +/** inp_out = (H + diag(weight)) \ inp_out (in place) */ +void sym_solve_( + DLTensor & inp_out , + const DLTensor & hessian , + const DLTensor & weight , + int stream = 0 +); + +/** out = inv(H) (both compact-symmetric) */ +void sym_invert( + DLTensor & out , + const DLTensor & hessian , + int stream = 0 +); + +/** hessian = inv(hessian) (in place, compact-symmetric) */ +void sym_invert_( + DLTensor & hessian , + int stream = 0 +); + +FF_NAMESPACE_END(FF) + +#endif // FF_LIB_POSDEF