diff --git a/include/fastfields/impl/cuda/distance_mesh.h b/include/fastfields/impl/cuda/distance_mesh.h index 0b75028..ea81742 100755 --- a/include/fastfields/impl/cuda/distance_mesh.h +++ b/include/fastfields/impl/cuda/distance_mesh.h @@ -736,6 +736,7 @@ template < CUGLOB inline void sdt_naive_kernel( offset_t nbatch , // Number of batch dimensions in coord scalar_t * dist , // (*batch) tensor -> Output placeholder for distance + index_t * nearest_vertex , // (*batch) tensor -> Output placeholder for index of nearest vertex const scalar_t * coord , // (*batch, D) tensor -> Coordinates at which to evaluate distance const scalar_t * _vertices , // (N, D) tensor -> All vertices const index_t * _faces , // (M, D) tensor -> All faces (face = D vertex indices) @@ -745,6 +746,7 @@ CUGLOB inline void sdt_naive_kernel( const offset_t * size , // [*batch] list -> Size of `dist` offset_t nb_faces , const offset_t * stride_dist , // [*batch] list -> Strides of `dist` + const offset_t * stride_nearest , // [*batch] list -> Strides of `nearest_vertex` const offset_t * stride_coord , // [*batch, D] list -> Strides of `coord` const offset_t * stride_vertices , // [N, D] list -> Strides of `vertices` const offset_t * stride_faces , // [M, D] list -> Strides of `faces` @@ -756,7 +758,6 @@ CUGLOB inline void sdt_naive_kernel( offset_t stride = blockDim.x * gridDim.x; using Klass = MeshDist; - using Node = typename Klass::Node; using FaceList = ConstStridedPointListSized; using VertexList = ConstStridedPointList; using NormalList = ConstStridedPointList; @@ -779,6 +780,12 @@ CUGLOB inline void sdt_naive_kernel( { offset_t offset_coord = index2offset(i, nbatch, size, stride_coord); offset_t offset_dist = index2offset(i, nbatch, size, stride_dist); + // `nearest_vertex` is optional, exactly as in `sdt_kernel` and in + // cpu-impl's `build_sdt_naive`: when it is null the stride array is + // null too and must not be read. + offset_t offset_nearest = 0; + if (nearest_vertex) + offset_nearest = index2offset(i, nbatch, size, stride_nearest); StaticPoint point( ConstStridedPoint(coord + offset_coord, stride_coord[nbatch])); @@ -789,7 +796,8 @@ CUGLOB inline void sdt_naive_kernel( faces, normfaces, normedges, - normvertices + normvertices, + nearest_vertex + offset_nearest ); } } @@ -1051,6 +1059,19 @@ sdt( normfaces_host = allocHost(nb_faces * ndim); normverts_host = allocHost(nb_vertices * ndim); normedges_host = allocHost(nb_faces * ndim * ndim); + // The vertex normals are ACCUMULATED into, not assigned: + // `MeshDistUtil::build_normals` does `normvertices[v].add_(normal)` + // once per incident face and normalises at the end. cpu-impl allocates + // them with `new scalar_t[...]()` for exactly this reason; `allocHost` + // is `cudaMallocHost`, which does not zero, so without this loop the + // pseudonormals accumulate on top of whatever the pinned allocation + // happened to contain. That does not fail loudly -- it perturbs the + // vertex/edge pseudonormals, i.e. the *sign* of the returned distance + // near vertices and edges. Face normals (`copy_`) and edge normals + // (built in a local map, then `copy_`) are assigned, so only this one + // buffer needs zeroing. + for (offset_t i = 0; i < nb_vertices * ndim; ++i) + normverts_host[i] = static_cast(0); // Build normals build_normals( @@ -1200,6 +1221,211 @@ sdt( ); } +// Naive signed mesh distance transform: no acceleration structure. Structural +// mirror of cpu-impl's `sdt_naive` -- it builds the vertex/face/edge normals on +// the host, uploads them, and launches `sdt_naive_kernel`, which brute-forces +// every face for every point. +// +// This is the reference the tree-accelerated `sdt` above is meant to be +// validated against on real hardware (fastfields-lib#5), so it deliberately +// shares as little machinery with it as possible: no BVH, no POD mirror, no +// per-lane traversal trace. What the two do share is the host-side normal +// construction, which is the same `fastfields-kernels` builder both backends +// use. +// +// Two differences from `sdt` above, both inherited from cpu-impl: +// * no `build_tree`, so the faces are never reordered and the contiguous +// device copy uploaded once at the top stays valid for the launch (`sdt` +// has to re-upload `faces_host` after the in-place BVH sort); +// * `nb_levels` / `treesize` / the trace buffer do not exist here at all. +template < + int _ndim, // Number of spatial dimensions + typename scalar_t, // Value data type + typename index_t, // Index (faces) data type + typename offset_t // Index/Stride data type +> +CUHOST inline void +sdt_naive( + offset_t nbatch, // Number of batch dimensions in coord + scalar_t * dist, // (*batch) tensor -> Output placeholder for distance + index_t * nearest_vertex, // (*batch) tensor -> Output placeholder for index of nearest vertex + const scalar_t * coord, // (*batch, D) tensor -> Coordinates at which to evaluate distance + const scalar_t * vertices, // (N, D) tensor -> All vertices + const index_t * faces, // (M, D) tensor -> All faces (face = D vertex indices) + const offset_t * size, // [*batch] list -> Size of `dist` + offset_t nb_faces, // M -> Number of faces + offset_t nb_vertices, // N -> Number of vertices + const offset_t * stride_dist, // [*batch] list -> Strides of `dist` + const offset_t * stride_nearest, // [*batch] list -> Strides of `nearest_vertex` + const offset_t * stride_coord, // [*batch, D] list -> Strides of `coord` + const offset_t * stride_vertices, // [N, D] list -> Strides of `vertices` + const offset_t * stride_faces , // [M, D] list -> Strides of `faces` + intptr_t stream = 0 // CUDA stream (0 == default stream) +) +{ + static const offset_t ndim = static_cast(_ndim); + const cudaStream_t s = (cudaStream_t)(std::intptr_t)stream; + + index_t * faces_device = nullptr; + scalar_t * verts_device = nullptr; + index_t * faces_host = nullptr; + scalar_t * verts_host = nullptr; + scalar_t * normfaces_host = nullptr; + scalar_t * normverts_host = nullptr; + scalar_t * normedges_host = nullptr; + scalar_t * normfaces_device = nullptr; + scalar_t * normverts_device = nullptr; + scalar_t * normedges_device = nullptr; + offset_t * stride_vec_device = nullptr; + offset_t * stride_mat_device = nullptr; + offset_t * stride_dist_device = nullptr; + offset_t * stride_nearest_device = nullptr; + offset_t * stride_coord_device = nullptr; + offset_t * size_device = nullptr; + + try + { + offset_t size_faces [2] = {nb_faces, ndim}; + offset_t size_verts [2] = {nb_vertices, ndim}; + offset_t stride_vec [2] = {ndim, 1}; + // Strides of the (M, D, D) edge-normal tensor: three values. + offset_t stride_mat [3] = {ndim*ndim, ndim, 1}; + // `faces` and `vertices` are 2-D tensors -- (M, D) and (N, D) -- so the + // rank handed to `copyTensorToContiguous` is 2, never `ndim`. + static const offset_t rank = static_cast(2); + + // The normals are built by a host function, so the mesh has to be on + // the host; the kernel then reads the contiguous device copies with + // `stride_vec`. Same route as `sdt` above: device-side gather into a + // contiguous buffer (which honours the caller's real input strides), + // then one D2H copy. + faces_device = + copyTensorToContiguous(rank, faces, size_faces, stride_faces, s); + verts_device = copyTensorToContiguous(rank, vertices, size_verts, + stride_vertices, s); + + faces_host = copyToHost(faces_device, nb_faces * ndim); + verts_host = copyToHost(verts_device, nb_vertices * ndim); + + // Allocate normals. `normverts_host` must be zeroed: the builder + // accumulates into it (see the same comment in `sdt`). + normfaces_host = allocHost(nb_faces * ndim); + normverts_host = allocHost(nb_vertices * ndim); + normedges_host = allocHost(nb_faces * ndim * ndim); + for (offset_t i = 0; i < nb_vertices * ndim; ++i) + normverts_host[i] = static_cast(0); + + // Build normals + build_normals<_ndim>( + normfaces_host, + normverts_host, + normedges_host, + faces_host, + verts_host, + nb_faces, + nb_vertices, + stride_vec, + stride_vec, + stride_mat, + stride_vec, + stride_vec + ); + + // Copy to device + normfaces_device = copyToDeviceAsync(normfaces_host, nb_faces * ndim, s); + normverts_device = copyToDeviceAsync(normverts_host, nb_vertices * ndim, s); + normedges_device = copyToDeviceAsync(normedges_host, nb_faces * ndim * ndim, s); + stride_vec_device = copyToDeviceAsync(stride_vec, 2, s); + stride_mat_device = copyToDeviceAsync(stride_mat, 3, s); + stride_dist_device = copyToDeviceAsync(stride_dist, nbatch, s); + // Guarded exactly as in `sdt`: `nearest_vertex` is optional and its + // stride array is null when it is, which `copyToDeviceAsync` cannot be + // handed (a null `cudaMemcpyAsync` source fails with + // `cudaErrorInvalidValue`). That is the default call pattern. + stride_nearest_device = nullptr; + if (nearest_vertex) + stride_nearest_device = + copyToDeviceAsync(stride_nearest, nbatch, s); + stride_coord_device = copyToDeviceAsync(stride_coord, nbatch + 1, s); + size_device = copyToDeviceAsync(size, nbatch, s); + + // The grid covers the number of points evaluated -- the batch element + // count, not `nbatch`, which is the batch *rank*. + offset_t numel = prod(size, nbatch); + int num_blocks = GET_BLOCKS(numel); + + // Compute SDT + sdt_naive_kernel + <<>> + ( + nbatch, + dist, + nearest_vertex, + coord, + verts_device, + faces_device, + normfaces_device, + normverts_device, + normedges_device, + size_device, + nb_faces, + stride_dist_device, + stride_nearest_device, + stride_coord_device, + stride_vec_device, + stride_vec_device, + stride_vec_device, + stride_vec_device, + stride_mat_device + ); + } + catch (const std::exception & e) + { + freeDevice( + faces_device, + verts_device, + normfaces_device, + normverts_device, + normedges_device, + stride_vec_device, + stride_mat_device, + stride_dist_device, + stride_nearest_device, + stride_coord_device, + size_device + ); + freeHost( + faces_host, + verts_host, + normfaces_host, + normverts_host, + normedges_host + ); + throw e; + } + + freeDevice( + faces_device, + verts_device, + normfaces_device, + normverts_device, + normedges_device, + stride_vec_device, + stride_mat_device, + stride_dist_device, + stride_nearest_device, + stride_coord_device, + size_device + ); + freeHost( + faces_host, + verts_host, + normfaces_host, + normverts_host, + normedges_host + ); +} + template CUHOST inline void sdt( offset_t nbatch , // Number of batch dimensions in coord @@ -1225,42 +1451,88 @@ CUHOST inline void sdt( const offset_t * stride_normedges , // [M, D, D] list intptr_t stream = 0 ) // CUDA stream (0 == default stream) { - // TODO(host-launcher): this precomputed-tree/normals `sdt` launcher is not - // implemented yet. The previous body was an erroneous copy/paste of a - // Euclidean distance-transform launcher (it referenced a non-existent + // NOT IMPLEMENTED, and deliberately still a throw. Read this before + // "finishing" it -- the missing piece is a contract, not code. + // + // History: the previous body was an erroneous copy/paste of a Euclidean + // distance-transform launcher (it referenced a non-existent // `allocBuffer`/`freeBuffers` callback API and identifiers `f`/`w`/`stride` // that are not parameters here, and launched `sdt_kernel` with the wrong - // signature). Implementing it correctly against the mesh `sdt_kernel` - // (line ~219) is tracked as the separate host-launcher task; see the - // complete `sdt` overload above (which builds the tree/normals itself). + // signature). It was replaced by this throw. + // + // What is actually missing. Mechanically the body is short -- upload the + // stride arrays, size the grid from `prod(size, nbatch)`, launch + // `sdt_kernel` -- but three things this signature does not state have to be + // decided before any of that can be *correct*, and none of them is + // discoverable from a call site: + // + // 1. `const void * tree`. `sdt_kernel` takes + // `const DeviceNode *`: a POD mirror, in + // *device* memory. A `void *` cannot distinguish that from the + // polymorphic host `MeshDist::Node[]` that `build_tree` produces -- + // and passing the latter is precisely the defect + // fastfields-cuda-impl#44 fixed in `sdt_kernel` by giving the + // parameter a real type. Nothing in this repository hands a caller a + // `DeviceNode` array: `flatten_tree` + the upload are private to the + // `sdt` overload above. So there is no producer, and no caller. + // 2. `void * treetrace` / `treesize`. The trace is per *lane*, not per + // point, and is interleaved across the launched lanes, so the buffer + // must be at least `GET_BLOCKS(numel) * CUDA_NUM_THREADS * treesize` + // bytes of device memory -- a size that depends on the launch + // configuration this function picks. A caller cannot size it without + // replicating that choice. + // 3. `faces` must be the BVH-sorted face list. `MeshDist::build_tree` + // reorders faces in place and its leaves store indices into the + // *sorted* order (see the comment on cpu-impl's `sdt`), so the faces + // and pseudonormals passed here have to be the ones the tree was + // built from, not the caller's originals. + // + // The options, for whoever picks this up (fastfields-lib#5): + // + // a. Give it a real signature and an in-repo caller: type `tree` as + // `const DeviceNode *` and `treetrace` as + // `char *`, then make the `sdt` overload above delegate to it -- which + // is exactly how cpu-impl is layered (`sdt` builds, `build_sdt` + // queries). The contract then has a compiled user, and the name should + // follow cpu-impl too (`build_sdt`). + // b. Drop this overload. It has no caller, no producer for its tree + // argument, and no test. + // + // Either is a design change with a reviewable blast radius, so it is not + // being made silently here. Until then this throws rather than pretending: + // a launcher that compiles but cannot be given a valid `tree` would be + // worse than none. throw std::logic_error("distance_mesh::sdt (precomputed tree) not implemented"); } // Top-level mesh distance dispatcher (mirrors cpu-impl distance_mesh::dt). // -// Coverage is PARTIAL, and deliberately so -- of the four branches cpu-impl -// dispatches, only one has a CUDA host launcher: +// Coverage is PARTIAL -- of the four branches cpu-impl dispatches, two have a +// CUDA host launcher: // // _signed naive CPU CUDA // ------- ----- ---------- ------------------------------------------- // true false sdt sdt() -> dispatched below -// true true sdt_naive sdt_naive_kernel -> no launcher, throws +// true true sdt_naive sdt_naive() -> dispatched below // false false udt udt_kernel -> no launcher, throws // false true udt_naive udt_naive_kernel -> no launcher, throws // -// The three device kernels exist and are type-checked by -// `tests/compile_probe_mesh.cu`, but none of them has a `CUHOST` launcher to -// build/upload the BVH and normals and size the grid. Writing those is tracked -// separately; until then those branches throw rather than silently returning -// garbage. See fastfields-lib#5. +// The two unsigned device kernels exist and are type-checked by +// `tests/impl-cuda/compile_probe_mesh.cu`, but neither has a `CUHOST` launcher +// to upload the mesh and size the grid. Writing those is tracked separately; +// until then those branches throw rather than silently returning garbage. See +// fastfields-lib#5. // -// Verification status of the branch that *is* dispatched: the per-element math -// (`MeshDist`, `signed_dist`, ...) is the shared `fastfields-kernels` source, -// compiled here and on the CPU alike and exhaustively covered by cpu-lib's -// `test_distance_mesh.cpp`. What is CUDA-only is the launcher glue -- the host -// BVH build, the POD flatten/upload, the grid sizing and the trace buffer -- -// and that is backed by nvcc compile+link evidence, since there is no GPU in -// CI. Same bar as every other module in cuda-lib's MODULES. +// Verification status of the branches that *are* dispatched: the per-element +// math (`MeshDist`, `signed_dist`, `signed_dist_naive`, ...) is the shared +// kernels source under `impl/kernels/`, compiled here and on the CPU alike and +// exhaustively covered by `tests/lib-cpu/test_distance_mesh.cpp`. What is +// CUDA-only is the launcher glue -- the host BVH build and normal build, the +// POD flatten/upload, the grid sizing and the trace buffer -- and that is +// backed by nvcc compile+link evidence only, since there is no GPU in CI. +// Nothing here has ever been executed: `sdt` vs `sdt_naive` agreement on real +// hardware is still the open acceptance bar of fastfields-lib#5. Same bar as +// every other module in src/lib-cuda's MODULES. template CUHOST inline void dt( @@ -1283,7 +1555,7 @@ dt( intptr_t stream = 0 ) { - // Signed, tree-accelerated: the one branch with a host launcher. + // Signed, tree-accelerated. // // `sdt` builds the BVH and the vertex/face/edge normals on the host, // uploads a POD mirror of the tree, sizes the grid from the batch element @@ -1298,15 +1570,20 @@ dt( nb_faces, nb_vertices, stride_dist, stride_nearest, stride_coord, stride_vertices, stride_faces, stream); - // The remaining three branches have device kernels but no host launcher -- - // see the table above this function. Throw a message that names the missing + // Signed, brute-force: the reference the branch above is meant to be + // validated against. Builds and uploads the normals, then launches + // `sdt_naive_kernel`; no BVH, no traversal trace. Same argument order. + if (_signed && naive) + return sdt_naive( + nbatch, dist, nearest_vertex, coord, vertices, faces, size, + nb_faces, nb_vertices, stride_dist, stride_nearest, stride_coord, + stride_vertices, stride_faces, stream); + + // The two unsigned branches have device kernels but no host launcher -- see + // the table above this function. Throw a message that names the missing // piece rather than the generic "not implemented", so a caller who hits one // knows which variant to ask for. (No `(void)` casts needed: every - // parameter is used by the dispatched branch above.) - if (_signed && naive) - throw std::logic_error( - "distance_mesh::dt (CUDA): the naive signed mesh distance " - "(sdt_naive) has no host launcher yet; use naive=false"); + // parameter is used by the dispatched branches above.) if (!_signed && !naive) throw std::logic_error( "distance_mesh::dt (CUDA): the unsigned mesh distance (udt) has " diff --git a/tests/impl-cuda/compile_probe_mesh.cu b/tests/impl-cuda/compile_probe_mesh.cu index aec1ecc..483add2 100644 --- a/tests/impl-cuda/compile_probe_mesh.cu +++ b/tests/impl-cuda/compile_probe_mesh.cu @@ -13,9 +13,9 @@ // * it covers `dt` for all four (signed x naive) flag combinations and both // the with- and without-`nearest_vertex` shapes, whereas the real build // only instantiates whatever `distance.cpp` happens to reference; -// * it covers `sdt_naive_kernel` / `udt_kernel` / `udt_naive_kernel` and -// `copy_faces`, none of which has a host launcher and none of which any -// build would otherwise instantiate; +// * it covers `udt_kernel` / `udt_naive_kernel` and `copy_faces`, none of +// which has a host launcher and none of which any build would otherwise +// instantiate; // * it applies `-Werror=shadow=local`, which the library build does not. // // It is compiled, never linked and never run (there is no GPU in CI). @@ -74,6 +74,29 @@ probe_sdt(offset_t nbatch, scalar_t * dist, index_t * nearest_vertex, stride_faces, 0); } +// The brute-force signed launcher. It has a caller now (`dt`'s signed+naive +// branch), but probe it directly for the same reason `sdt` is probed directly: +// `distance.cpp` only instantiates whatever dtype/dim combinations it happens +// to reference, and an argument-order slip between two `offset_t` parameters +// (`nb_faces` / `nb_vertices`) is invisible to the type system either way, so +// the launch itself has to be compiled for every combination the dispatcher +// can select. +template +static void +probe_sdt_naive(offset_t nbatch, scalar_t * dist, index_t * nearest_vertex, + const scalar_t * coord, const scalar_t * vertices, + const index_t * faces, const offset_t * size, + offset_t nb_faces, offset_t nb_vertices, + const offset_t * stride_dist, const offset_t * stride_nearest, + const offset_t * stride_coord, const offset_t * stride_vertices, + const offset_t * stride_faces) +{ + M::sdt_naive( + nbatch, dist, nearest_vertex, coord, vertices, faces, size, nb_faces, + nb_vertices, stride_dist, stride_nearest, stride_coord, stride_vertices, + stride_faces, 0); +} + // The launcher above is not what real callers reach: cuda-lib's `distance.cpp` // dispatches into `distance_mesh::dt`, which then selects a launcher. Probing // only `sdt` left `dt` itself uninstantiated, so a mistake in the forwarding @@ -124,9 +147,12 @@ const void * ff_compile_probe_mesh_udt_kernel(int which) } } -// Same reasoning for the two *naive* kernels (the brute-force references the -// signed/unsigned BVH paths are meant to be validated against) -- neither has a -// host launcher, so nothing instantiated them either. +// The two *naive* kernels -- the brute-force references the signed/unsigned BVH +// paths are meant to be validated against. `udt_naive_kernel` still has no host +// launcher, so nothing else instantiates it. `sdt_naive_kernel` does now +// (`M::sdt_naive`, probed above); taking its address as well costs nothing and +// keeps the pair symmetric, so the day `udt_naive` gains a launcher this block +// does not need rethinking. template static const void * probe_naive_kernels(int which) { @@ -182,6 +208,19 @@ void ff_compile_probe_mesh_sdt( probe_sdt<3, double, long, long>(nb64, d64, nv64, c64, v64, f64, sz64, nb64, nb64, st64, st64, st64, st64, st64); + probe_sdt_naive<2, float, int, int>(nb32, d32, nv32, c32, v32, f32, sz32, + nb32, nb32, st32, st32, st32, st32, + st32); + probe_sdt_naive<3, float, int, int>(nb32, d32, nv32, c32, v32, f32, sz32, + nb32, nb32, st32, st32, st32, st32, + st32); + probe_sdt_naive<2, double, long, long>(nb64, d64, nv64, c64, v64, f64, sz64, + nb64, nb64, st64, st64, st64, st64, + st64); + probe_sdt_naive<3, double, long, long>(nb64, d64, nv64, c64, v64, f64, sz64, + nb64, nb64, st64, st64, st64, st64, + st64); + // The real entry point, same four type combinations. probe_dt<2, float, int, int>(nb32, d32, nv32, c32, v32, f32, sz32, nb32, nb32, st32, st32, st32, st32, st32); @@ -199,6 +238,9 @@ void ff_compile_probe_mesh_sdt( probe_dt<3, float, int, int>(nb32, d32, static_cast(nullptr), c32, v32, f32, sz32, nb32, nb32, st32, nullptr, st32, st32, st32); + probe_sdt_naive<3, float, int, int>(nb32, d32, static_cast(nullptr), + c32, v32, f32, sz32, nb32, nb32, st32, + nullptr, st32, st32, st32); probe_copy_faces<2, int, int>(nb32, f32, st32); probe_copy_faces<3, int, int>(nb32, f32, st32);