cuddl Namespace Reference

Batch operations over contiguous device rows. More...

Namespaces
namespace	a48
namespace	detail

Classes
struct	all_to_all_shard
	Part cuDDL of count disjoint parts of an all-to-all search. More...
class	all_to_all_tile_queue
	Thread-safe source of the all-to-all tiles one process searches. More...
struct	batch_result_tile
	Non-owning view of one batch-search tile's results. More...
struct	batch_search_requirements
	Caller-owned storage requirements for one bounded query tile. More...
struct	batch_search_result
	Stable query/reference IDs and their exact query-relative bucket comparison counts. More...
struct	CudaError
	CUDA failure payload: status code, call site, and formatted message. More...
class	device_blacklist
	Reusable, move-only device lookup. More...
struct	device_span
	A span that is assumed to point to device-accessible memory. More...
struct	Error
	Error payload carried in Result on failure. More...
struct	genome_slice
	One rank's share of a genome collection split across an MPI communicator. More...
class	host_batch_results
	Host copy of a batch_result_tile holding only the passing results, packed in query-major, reference order, in pinned memory; tile views it. More...
struct	hybrid_cardinality_estimates
	Experimental DDL cardinality estimates computed from one register scan. More...
struct	InvalidArgumentError
	Invalid-argument failure payload. More...
class	kmer_blacklist
	Immutable canonical k-mers excluded from construction, independent of sketch size. More...
struct	packed_pairwise_counts
	One pair's exact bucket comparison counts in 64 bits, the storage and transfer form of batch-search results. More...
struct	pairwise_counts
	Raw per-bucket classification counts for one sketch pair comparison. More...
struct	pairwise_summary
	Pairwise counts with an optional compile-time cardinality field. More...
struct	path_build_options
	Knobs for a build the caller feeds with FASTA/FASTQ paths. More...
class	query_sketch_batch
	Query sketches held on the device, in the shape a search consumes them. More...
struct	reference_build_statistics
	How a build moved sequence bytes to the device. More...
class	reference_database
	Move-only owner of one immutable contiguous reference database. More...
class	reference_database_file
	Host-owned reference sketches and labels, ready for binary storage or GPU upload. More...
struct	reference_database_metadata
	Recorded metadata for an immutable reference database. More...
class	reference_index
	Move-only acceleration data for an immutable reference database. More...
class	reference_index_file
	Index-file persistence and the validation boundary for reference_index. More...
struct	reference_search_result
	Stable reference ID and exact query-relative bucket comparison counts. More...
struct	register_layout
	Compile-time split of a 16-bit DDL score into exponent and mantissa fields. More...
struct	replica_search
	How a replicated_database search splits and delivers its tiles. More...
struct	replica_tile
	One all-to-all tile as a replicated_database search delivers it. More...
class	replicated_database
	One reference database replicated on explicitly chosen devices, searched all-to-all by all of them at once. More...
struct	ResourceError
	Resource failure payload: allocation, I/O, or capacity errors. More...
struct	Result
	Fallible API result: cuda::std::expected<T, Error> with cuDDL factories. More...
struct	Result< void >
class	result_file
	Writes search results to files while the GPUs search on. More...
struct	result_file_options
struct	result_file_statistics
	What a result_file wrote and how long its callers waited for it. More...
struct	score_compatibility
	Construction parameters of compatible score rows. More...
struct	search_options
	Matching-bucket threshold, independent of whether an index accelerates the search. More...
struct	sequence_build_options
	Knobs for a build the caller feeds with bases it already holds. More...
class	shared_tile_counter
	A tile counter shared by every rank of a communicator, for replica_search::next_tile. More...
struct	similarity_record
	One pair of a thresholded search in result_format::similarity. More...
class	sketch
	Move-only owning DDL sketch backed by one contiguous device allocation. More...
struct	SourceLocation
	File, line, and column for a C++ call site or FASTX input position. More...
class	top_k_search
	Keeps each genome's k most similar genomes from search tiles on the device. More...

Typedefs
using	default_register_layout = register_layout< 6, 10 >
using	fasta_parse_result = detail::fasta_parse_result
	Result of parsing one FASTA/FASTQ file into packed canonical k-mers.
using	transfer_mode = detail::transfer_mode
	How a build gets its sequence bytes to the device.
using	decompression_backend = detail::decompression_backend
	Selects decompression independently of host transfers.
using	sequence_record = detail::sequence_record
	One record of a genome the caller already holds: bases only.
using	sequence_genome = detail::sequence_genome
	One genome as its records, in order. No records keeps the genome's ID with a zero row.

Enumerations
enum class	ErrorCategory : uint8_t { cuda , invalid_argument , resource }
	Error category for Result failures. More...
enum class	index_storage { automatic , dense , sparse }
	Dense offsets favor query latency; sparse sorted keys reduce index memory and build time. More...
enum class	result_format : uint8_t { counts , similarity }
	How a result_file stores each pair. More...

Functions
template<size_t BucketCount>
Result< void >	compare_batch_async (device_span< uint32_t const > left_rows, device_span< uint32_t const > right_rows, device_span< pairwise_summary > outputs, cuda::stream_ref stream)
	Compares corresponding rows from two contiguous register batches.
template<size_t BucketCount, typename Layout = default_register_layout>
Result< void >	cardinality_batch_async (device_span< uint32_t const > sketches, device_span< uint64_t > empty_out, device_span< double > estimates_out, cuda::stream_ref stream)
	Computes the cardinality of every stored sketch.
template<size_t BucketCount>
Result< void >	extract_scores_batch_async (device_span< uint32_t const > sketches, device_span< uint16_t > scores_out, cuda::stream_ref stream)
	Copies the winning score of every register into compact row-major scores.
Result< genome_slice >	slice_genomes (std::filesystem::path const &input, MPI_Comm comm, unsigned workers=default_parser_workers())
	Lists a genome collection once and splits it into one contiguous slice per rank of comm, of about equal file bytes.
Result< reference_database_file >	gather_sketches (reference_database_file const &mine, genome_slice const &slice, MPI_Comm comm)
	Gathers every rank's sketches of its genome_slice into one file on every rank.
template<uint32_t K, size_t BucketCount, typename Layout , typename OnTile >
Result< void >	search_all_to_all (replicated_database< K, BucketCount, Layout > &group, OnTile &&on_tile, MPI_Comm comm)
	Compares every pair of group's references, sharing the tiles among the ranks of comm; collective.
template<uint32_t K, size_t BucketCount, typename Layout , typename OnTile >
Result< void >	search_all_to_all (replicated_database< K, BucketCount, Layout > &group, OnTile &&on_tile, search_options options, MPI_Comm comm)
	The threshold form of search_all_to_all: only pairs passing options reach on_tile.
void	cuda_abort_on_error (cudaError_t error, std::source_location location=std::source_location::current())
	Terminates the process when error is not cudaSuccess (for RAII teardown only).
cuda::std::unexpected< Error >	Err (Error error)
	Failure return; converts to any Result<T> via cuda::std::unexpected.
Result< void >	Ok () noexcept
	Success return for Result<void>; same as return {}.
Result< void >	cuda_try (cudaError_t error, std::source_location location=std::source_location::current())
	Checks a CUDA runtime call and returns an error on failure.
template<std::invocable Operation>
auto	cuda_try (Operation &&operation, std::source_location location=std::source_location::current())
	Converts CCCL exceptions to the same Result errors as CUDA status-returning APIs.
Result< fasta_parse_result >	parse_fasta_file (std::string const &path, uint32_t k, unsigned threads=std::thread::hardware_concurrency())
	Parses a FASTA/FASTQ file's sequence bases into packed canonical k-mers of length k.
__host__ __device__ constexpr float	weighted_kmer_identity (pairwise_counts const &counts) noexcept
	Weighted k-mer identity of a pair, equal / max(equal + min(lower, higher), 6), as sketch::wkid computes it, or 0 when no bucket is comparable.
Result< void >	download (batch_result_tile const &tile, cuda::stream_ref stream, host_batch_results &host)
	Copies a device tile's passing results into host's pinned buffers, growing them only when needed, and synchronizes stream.
Result< host_batch_results >	download (batch_result_tile const &tile, cuda::stream_ref stream)
	Copies a device tile's passing results to new pinned host memory; see the overload that downloads into reusable host_batch_results.
template<typename F >
Result< void >	for_each_passing (batch_result_tile const &tile, cuda::stream_ref stream, F &&f)
	Calls f on the host with every passing result of a device tile, query-major and in reference order; synchronizes stream.
unsigned	default_parser_workers ()
	Loaders a path build runs unless the caller asks for another count.
Result< std::vector< std::filesystem::path > >	genome_paths (std::filesystem::path const &directory)
	The FASTA/FASTQ files under directory, searched recursively and sorted by path, which is the order that numbers reference IDs.
Result< std::vector< std::filesystem::path > >	list_genomes (std::filesystem::path const &input)
	Paths from a directory (see genome_paths) or from a text file with one path per line, in the file's order.
template<uint32_t K, size_t BucketCount, typename Layout = default_register_layout>
Result< cuda::device_buffer< uint32_t > >	build_sketch_store (std::span< std::filesystem::path const > paths, cuda::stream_ref stream, path_build_options options={})
	Sketches one genome per plain or gzip/BGZF FASTA/FASTQ file into a device store.
score_compatibility	decoded_compatibility (uint8_t k, uint16_t buckets, uint8_t exponent_bits, uint64_t seed)
	Builds a cuDDL score-compatibility object describing a decoded BBTools A48 database.
std::vector< batch_search_result >	select_top_k (std::vector< batch_search_result > rows, uint32_t k)
	The k best rows of each query among rows, as top_k_search ranks them, best first and query by query.

Detailed Description

Batch operations over contiguous device rows.

Every function takes flat row-major spans of BucketCount uint32_t registers per sketch, the layout of a single sketch allocation and of the stores the file builders write, and runs the whole batch in one kernel launch, so a collection of N rows costs one launch instead of N.

Inputs and outputs must remain valid until stream completes.

Typedef Documentation

decompression_backend

using cuddl::decompression_backend = typedef detail::decompression_backend

Selects decompression independently of host transfers.

Definition at line 155 of file reference_database_file.cuh.

default_register_layout

using cuddl::default_register_layout = typedef register_layout<6, 10>

Definition at line 22 of file register.cuh.

fasta_parse_result

using cuddl::fasta_parse_result = typedef detail::fasta_parse_result

Result of parsing one FASTA/FASTQ file into packed canonical k-mers.

Definition at line 16 of file fastx.hpp.

sequence_genome

using cuddl::sequence_genome = typedef detail::sequence_genome

One genome as its records, in order. No records keeps the genome's ID with a zero row.

Definition at line 161 of file reference_database_file.cuh.

sequence_record

using cuddl::sequence_record = typedef detail::sequence_record

One record of a genome the caller already holds: bases only.

Definition at line 158 of file reference_database_file.cuh.

transfer_mode

using cuddl::transfer_mode = typedef detail::transfer_mode

How a build gets its sequence bytes to the device.

Definition at line 152 of file reference_database_file.cuh.

Enumeration Type Documentation

ErrorCategory

enum class cuddl::ErrorCategory : uint8_t

strong

Error category for Result failures.

Enumerator
cuda
invalid_argument
resource

Definition at line 57 of file error.hpp.

                         : uint8_t {
    cuda,
    invalid_argument,
    resource,
};

index_storage

enum class cuddl::index_storage

strong

Dense offsets favor query latency; sparse sorted keys reduce index memory and build time.

Automatic picks whichever needs fewer device bytes for the database's reference count.

Enumerator
automatic
dense
sparse

Definition at line 60 of file reference_database.cuh.

{ automatic, dense, sparse };

result_format

enum class cuddl::result_format : uint8_t

strong

How a result_file stores each pair.

Enumerator
counts	One 24-byte batch_search_result per pair: query ID, reference ID and the lower, equal, higher and both-empty bucket counts, six little-endian uint32.
similarity	The pair's weighted_kmer_identity as a 4-byte float, computed on the GPU. A tile that holds every pair is written as its first query ID and query count (two uint32) followed by one float per pair in the tile's order: query-major, each query against every later reference in an all-to-all search and against every reference otherwise. Every other pair is a 12-byte similarity_record.

Definition at line 35 of file result_file.cuh.

                         : uint8_t {
    counts,
    similarity,
};

Function Documentation

build_sketch_store()

template<uint32_t K, size_t BucketCount, typename Layout = default_register_layout>

Result< cuda::device_buffer< uint32_t > > cuddl::build_sketch_store	(	std::span< std::filesystem::path const >	paths,
		cuda::stream_ref	stream,
		path_build_options	options = {}
	)

Sketches one genome per plain or gzip/BGZF FASTA/FASTQ file into a device store.

Each genome gets BucketCount registers, the row layout the batch operations in batch.cuh consume, in the order of paths. The rows never leave the device. The build is synchronous.

Definition at line 685 of file reference_database_file.cuh.

                                 {}
) try {
    return [&]() -> Result<cuda::device_buffer<uint32_t>> {
        if (options.blacklist) CUDDL_TRY(options.blacklist->get().validate(K, BucketCount, stream));
 
        auto store = cuda::make_device_buffer<uint32_t>(
            stream, stream.device(), paths.size() * BucketCount, cuda::no_init
        );
 
        CUDDL_TRY((detail::stage_paths<K, BucketCount, Layout>(
            paths,
            stream,
            options.staging_bytes,
            options.parser_workers,
            options.transfer,
            options.decompression,
            options.direct_reads,
            options.statistics,
            [&](device_span<uint32_t const> group, size_t base, size_t count) -> Result<void> {
                CUDDL_CUDA_TRY(
                    cuda::copy_bytes(
                        stream,
                        group,
                        device_span<uint32_t>{
                            store.data() + base * BucketCount, count * BucketCount
                        }
                    )
                );
                return Ok();
            },
            options.blacklist ? options.blacklist->get().view() : cuda::std::nullopt
        )));
        return store;
    }();
} catch (cuda::cuda_error const& error) {
    return Err(Error::cuda(static_cast<cudaError_t>(error.status())));
} catch (std::system_error const& error) {
    return Err(Error::resource(error.what()));
} catch (std::bad_alloc const& error) {
    return Err(Error::resource(error.what()));
}

cardinality_batch_async()

template<size_t BucketCount, typename Layout = default_register_layout>

Result< void > cuddl::cardinality_batch_async ( device_span< uint32_t const >  sketches, device_span< uint64_t >  empty_out, device_span< double >  estimates_out, cuda::stream_ref  stream  )

inline

Computes the cardinality of every stored sketch.

empty_out receives the empty-register count per row and estimates_out the cardinality estimate per row.

Definition at line 80 of file batch.cuh.

  {
    static_assert(BucketCount >= (size_t{1} << 11) && BucketCount <= (size_t{1} << 13));
    static_assert((BucketCount & (BucketCount - 1U)) == 0U);
 
    if (sketches.size() % BucketCount != 0U) {
        return Err(Error::invalid_argument("sketch store must contain whole sketches"));
    }
    auto const row_count = sketches.size() / BucketCount;
    if (empty_out.size() < row_count || estimates_out.size() < row_count) {
        return Err(Error::invalid_argument("sketch store outputs are too small"));
    }
    if (row_count == 0U) {
        return Ok();
    }
    if (sketches.data() == nullptr || empty_out.data() == nullptr ||
        estimates_out.data() == nullptr) {
        return Err(Error::invalid_argument("nonempty sketch store buffers must not be null"));
    }
    if (row_count > detail::maximum_batch_rows) {
        return Err(Error::invalid_argument("sketch store has too many rows for one launch"));
    }
 
    detail::batch_cardinality_kernel<BucketCount, Layout>
        <<<static_cast<uint32_t>(row_count), detail::block_size, 0, stream.get()>>>(
            sketches.data(),
            static_cast<uint32_t>(row_count),
            empty_out.data(),
            estimates_out.data()
        );
    return cuda_try(cudaGetLastError());
}

Here is the call graph for this function:

compare_batch_async()

template<size_t BucketCount>

Result< void > cuddl::compare_batch_async	(	device_span< uint32_t const >	left_rows,
		device_span< uint32_t const >	right_rows,
		device_span< pairwise_summary >	outputs,
		cuda::stream_ref	stream
	)

Compares corresponding rows from two contiguous register batches.

The input spans are row-major and must contain the same whole number of BucketCount rows. outputs must hold at least one summary per row pair. Inputs and outputs must remain valid until stream completes. The whole batch runs in one kernel launch.

Definition at line 35 of file batch.cuh.

  {
    static_assert(BucketCount >= (size_t{1} << 11) && BucketCount <= (size_t{1} << 13));
    static_assert((BucketCount & (BucketCount - 1U)) == 0U);
 
    if (left_rows.size() != right_rows.size()) {
        return Err(Error::invalid_argument("batch register spans must have equal sizes"));
    }
    if (left_rows.size() % BucketCount != 0U) {
        return Err(Error::invalid_argument("batch register spans must contain whole rows"));
    }
    auto const pair_count = left_rows.size() / BucketCount;
    if (outputs.size() < pair_count) {
        return Err(Error::invalid_argument("batch output span is too small"));
    }
    if (pair_count != 0U && (left_rows.data() == nullptr || right_rows.data() == nullptr ||
                             outputs.data() == nullptr)) {
        return Err(Error::invalid_argument("nonempty batch buffers must not be null"));
    }
    if (pair_count == 0U) {
        return Ok();
    }
 
    constexpr size_t warps_per_block = detail::block_size / 32U;
    constexpr size_t maximum_blocks = 65535U;
    auto const blocks =
        std::min((pair_count + warps_per_block - 1U) / warps_per_block, maximum_blocks);
    detail::batch_summary_kernel<BucketCount>
        <<<static_cast<uint32_t>(blocks), detail::block_size, 0, stream.get()>>>(
            left_rows.data(), right_rows.data(), pair_count, outputs.data()
        );
    return cuda_try(cudaGetLastError());
}

Here is the call graph for this function:

cuda_abort_on_error()

void cuddl::cuda_abort_on_error ( cudaError_t  error, std::source_location  location = std::source_location::current()  )

inline

Terminates the process when error is not cudaSuccess (for RAII teardown only).

Definition at line 144 of file error.hpp.

  {
    if (error == cudaSuccess) {
        return;
    }
    const Error err = Error::cuda(error, location);
    std::fputs(err.message().c_str(), stderr);
    std::fputc('\n', stderr);
    std::abort();
}

Here is the call graph for this function:

cuda_try() [1/2]

Result< void > cuddl::cuda_try ( cudaError_t  error, std::source_location  location = std::source_location::current()  )

inline

Checks a CUDA runtime call and returns an error on failure.

Definition at line 187 of file error.hpp.

                                                                                       {
    if (error == cudaSuccess) {
        return Ok();
    }
    return Err(Error::cuda(error, location));
}

Here is the call graph for this function:

Here is the caller graph for this function:

cuda_try() [2/2]

template<std::invocable Operation>

auto cuddl::cuda_try	(	Operation &&	operation,
		std::source_location	location = std::source_location::current()
	)

Converts CCCL exceptions to the same Result errors as CUDA status-returning APIs.

Definition at line 197 of file error.hpp.

                                                                                           {
    using Value = std::invoke_result_t<Operation>;
    using Return = Result<std::conditional_t<std::is_same_v<Value, cudaError_t>, void, Value>>;
    try {
        if constexpr (std::is_void_v<Value>) {
            std::forward<Operation>(operation)();
            return Return::ok();
        } else if constexpr (std::is_same_v<Value, cudaError_t>) {
            return cuda_try(std::forward<Operation>(operation)(), location);
        } else {
            return Return::ok(std::forward<Operation>(operation)());
        }
    } catch (cuda::cuda_error const& error) {
        return Return::err(Error::cuda(static_cast<cudaError_t>(error.status()), location));
    } catch (std::bad_alloc const& error) {
        return Return::err(Error::resource(error.what()));
    }
}

Here is the call graph for this function:

decoded_compatibility()

score_compatibility cuddl::decoded_compatibility ( uint8_t  k, uint16_t  buckets, uint8_t  exponent_bits, uint64_t  seed  )

inline

Builds a cuDDL score-compatibility object describing a decoded BBTools A48 database.

Decoded rows come from the BBTools encoder/hash pipeline, not the cuDDL one. The encoder and hash identity fields therefore carry nonzero sentinels (2) distinct from cuDDL's own construction constants, while the mantissa/exponent split and bucket/key masks are taken from the asset's header so decoded rows pass cuDDL's normal validation paths.

Definition at line 15 of file refseq_parity.hpp.

                                                                                         {
    score_compatibility compat;
    compat.kmer_length = k;
    compat.bucket_count = buckets;
    compat.indexed_bucket_count = buckets;
    compat.score_encoder_identity = 2U;
    compat.exponent_bits = exponent_bits;
    compat.mantissa_bits = 16U - exponent_bits;
    compat.hash_identity = 2U;
    compat.hash_seed = seed;
    compat.canonicalisation_policy = 1U;
    compat.blacklist_identity = 0U;
    compat.blacklist_version = 0U;
    compat.key_mask = 0xffffU;
    return compat;
}

default_parser_workers()

unsigned cuddl::default_parser_workers ( )

inline

Loaders a path build runs unless the caller asks for another count.

Machine thread count, falling back to 8 when the count is unavailable.

Definition at line 166 of file reference_database_file.cuh.

                                                       {
    auto const hardware = std::thread::hardware_concurrency();
    return hardware == 0U ? 8U : hardware;
}

download() [1/2]

Result< host_batch_results > cuddl::download ( batch_result_tile const &  tile, cuda::stream_ref  stream  )

inline

Copies a device tile's passing results to new pinned host memory; see the overload that downloads into reusable host_batch_results.

Definition at line 583 of file reference_database.cuh.

                                                               {
    host_batch_results host(stream);
    CUDDL_TRY(download(tile, stream, host));
    return Result<host_batch_results>::ok(std::move(host));
}

Here is the call graph for this function:

download() [2/2]

Result< void > cuddl::download ( batch_result_tile const &  tile, cuda::stream_ref  stream, host_batch_results &  host  )

inline

Copies a device tile's passing results into host's pinned buffers, growing them only when needed, and synchronizes stream.

Exhaustive tiles are copied as they are, since every slot passes in order. Threshold tiles are packed on the device: a scan of the pass-bit words gives each word's first packed position, and a gather writes the passing records straight into the host buffer.

Definition at line 427 of file reference_database.cuh.

                                                                                         {
    detail::device_scope const scope{stream};
    host.layout_ = tile;
    host.layout_.word_offsets = nullptr;
    host.has_match_counts_ = tile.match_counts != nullptr;
    if (tile.pass_bits == nullptr) {
        auto const slots = static_cast<size_t>(tile.slot_count());
        host.count_ = slots;
        host.reserve(host.results_, slots, stream);
        if (slots != 0U) {
            CUDDL_CUDA_TRY(
                cuda::copy_bytes(
                    stream,
                    cuda::std::span{tile.results, slots},
                    cuda::std::span{host.results_.data(), slots}
                )
            );
        }
        if (host.has_match_counts_) {
            host.reserve(host.match_counts_, slots, stream);
            if (slots != 0U) {
                CUDDL_CUDA_TRY(
                    cuda::copy_bytes(
                        stream,
                        cuda::std::span{tile.match_counts, slots},
                        cuda::std::span{host.match_counts_.data(), slots}
                    )
                );
            }
        }
        CUDDL_CUDA_TRY(stream.sync());
        return Ok();
    }
 
    auto const words = static_cast<size_t>(tile.query_count) * tile.words_per_query();
    host.reserve(host.pass_bits_, words, stream);
    host.reserve(host.word_offsets_, words, stream);
    if (words == 0U) {
        host.count_ = 0U;
        CUDDL_CUDA_TRY(stream.sync());
        return Ok();
    }
    auto const popcounts = cuda::make_transform_iterator(tile.pass_bits, detail::word_popcount{});
    size_t scan_bytes = 0;
    CUDDL_CUDA_TRY(
        cub::DeviceScan::ExclusiveSum(
            nullptr,
            scan_bytes,
            popcounts,
            host.word_offsets_.data(),
            static_cast<int64_t>(words),
            stream.get()
        )
    );
    auto scan_storage = CUDDL_CUDA_TRY(
        cuda::make_device_buffer<uint8_t>(stream, stream.device(), scan_bytes, cuda::no_init)
    );
    CUDDL_CUDA_TRY(
        cub::DeviceScan::ExclusiveSum(
            scan_storage.data(),
            scan_bytes,
            popcounts,
            host.word_offsets_.data(),
            static_cast<int64_t>(words),
            stream.get()
        )
    );
    CUDDL_CUDA_TRY(
        cuda::copy_bytes(
            stream,
            cuda::std::span{tile.pass_bits, words},
            cuda::std::span{host.pass_bits_.data(), words}
        )
    );
    CUDDL_CUDA_TRY(stream.sync());
    host.count_ = static_cast<uint64_t>(host.word_offsets_.data()[words - 1U]) +
                  static_cast<uint64_t>(cuda::std::popcount(host.pass_bits_.data()[words - 1U]));
    auto const count = static_cast<size_t>(host.count_);
    host.reserve(host.results_, count, stream);
    if (host.has_match_counts_) host.reserve(host.match_counts_, count, stream);
    if (count == 0U) {
        CUDDL_CUDA_TRY(stream.sync());
        return Ok();
    }
    // Passing records are gathered into bounded device staging, so each chunk crosses the bus as
    // one bulk copy rather than as scattered small writes.
    constexpr size_t staging_records = size_t{1} << 23U;
    auto const staged = std::min(count, staging_records);
    auto staging = CUDDL_CUDA_TRY(
        cuda::make_device_buffer<packed_pairwise_counts>(
            stream, stream.device(), staged, cuda::no_init
        )
    );
    auto staging_matches = CUDDL_CUDA_TRY(
        cuda::make_device_buffer<uint32_t>(
            stream, stream.device(), tile.match_counts == nullptr ? 0U : staged, cuda::no_init
        )
    );
    auto const kernel = tile.upper_triangle
                            ? detail::gather_passing_kernel<true, packed_pairwise_counts>
                            : detail::gather_passing_kernel<false, packed_pairwise_counts>;
    auto const* offsets = host.word_offsets_.data();
    for (size_t first_word = 0U; first_word < words;) {
        auto const base = offsets[first_word];
        // The last word whose results all fit the staging buffer ends the chunk.
        auto const last_word = static_cast<size_t>(
            std::upper_bound(offsets + first_word + 1U, offsets + words, base + staged - 32U) -
            offsets
        );
        auto const end = last_word == words ? count : static_cast<size_t>(offsets[last_word]);
        kernel<<<
            detail::warp_grid_blocks(last_word - first_word),
            detail::block_size,
            0,
            stream.get()>>>(
            tile.results,
            tile.match_counts,
            tile.pass_bits,
            host.word_offsets_.data(),
            tile.first_query_id,
            tile.reference_count,
            first_word,
            last_word,
            base,
            staging.data(),
            tile.match_counts == nullptr ? nullptr : staging_matches.data()
        );
        CUDDL_CUDA_TRY(cudaGetLastError());
        auto const records = end - base;
        if (records != 0U) {
            CUDDL_CUDA_TRY(
                cuda::copy_bytes(
                    stream,
                    cuda::std::span{staging.data(), records},
                    cuda::std::span{host.results_.data() + base, records}
                )
            );
            if (tile.match_counts != nullptr) {
                CUDDL_CUDA_TRY(
                    cuda::copy_bytes(
                        stream,
                        cuda::std::span{staging_matches.data(), records},
                        cuda::std::span{host.match_counts_.data() + base, records}
                    )
                );
            }
        }
        first_word = last_word;
    }
    CUDDL_CUDA_TRY(stream.sync());
    return Ok();
}

Here is the caller graph for this function:

Err()

cuda::std::unexpected< Error > cuddl::Err ( Error  error )

inline

Failure return; converts to any Result<T> via cuda::std::unexpected.

Definition at line 176 of file error.hpp.

                                                               {
    return cuda::std::unexpected<Error>(std::move(error));
}

extract_scores_batch_async()

template<size_t BucketCount>

Result< void > cuddl::extract_scores_batch_async ( device_span< uint32_t const >  sketches, device_span< uint16_t >  scores_out, cuda::stream_ref  stream  )

inline

Copies the winning score of every register into compact row-major scores.

scores_out holds one uint16_t-per-register row per stored sketch.

Definition at line 123 of file batch.cuh.

  {
    static_assert(BucketCount >= (size_t{1} << 11) && BucketCount <= (size_t{1} << 13));
    static_assert((BucketCount & (BucketCount - 1U)) == 0U);
 
    if (sketches.size() % BucketCount != 0U) {
        return Err(Error::invalid_argument("sketch store must contain whole sketches"));
    }
    auto const row_count = sketches.size() / BucketCount;
    if (scores_out.size() < row_count * BucketCount) {
        return Err(Error::invalid_argument("sketch store outputs are too small"));
    }
    if (row_count == 0U) {
        return Ok();
    }
    if (sketches.data() == nullptr || scores_out.data() == nullptr) {
        return Err(Error::invalid_argument("nonempty sketch store buffers must not be null"));
    }
    if (row_count > detail::maximum_batch_rows) {
        return Err(Error::invalid_argument("sketch store has too many rows for one launch"));
    }
 
    detail::batch_scores_kernel<BucketCount>
        <<<static_cast<uint32_t>(row_count), detail::block_size, 0, stream.get()>>>(
            sketches.data(), static_cast<uint32_t>(row_count), scores_out.data()
        );
    return cuda_try(cudaGetLastError());
}

Here is the call graph for this function:

for_each_passing()

template<typename F >

Result< void > cuddl::for_each_passing	(	batch_result_tile const &	tile,
		cuda::stream_ref	stream,
		F &&	f
	)

Calls f on the host with every passing result of a device tile, query-major and in reference order; synchronizes stream.

Results cross the bus in chunks through two pinned buffers, so the copy of chunk k + 1 runs while f visits chunk k. Threshold tiles are packed on the device first, so only passing results are copied. This streams where download materializes the whole tile.

Definition at line 615 of file reference_database.cuh.

                                                                              {
    detail::device_scope const scope{stream};
    using pinned_counts = cuda::
        buffer<packed_pairwise_counts, cuda::mr::host_accessible, cuda::mr::device_accessible>;
    auto const memory =
        detail::device_pinned_pool(stream.device(), detail::pinned_placement::device_node);
    constexpr auto chunk_records = detail::visit_chunk_records;
    auto const threshold = tile.pass_bits != nullptr;
    auto const words = static_cast<size_t>(tile.query_count) * tile.words_per_query();
 
    // Chunk plan, and for threshold tiles the host copies of the pass bits and word offsets.
    std::vector<detail::visit_chunk> chunks;
    std::vector<uint32_t> pass_bits;
    std::vector<uint32_t> offsets;
    cuda::device_buffer<uint32_t> device_offsets(
        stream, cuda::device_default_memory_pool(stream.device())
    );
    if (!threshold) {
        size_t record = 0U;
        detail::visit_chunk chunk{};
        for (uint32_t query = 0; query < tile.query_count; ++query) {
            auto const query_id = tile.first_query_id + query;
            auto const row = static_cast<size_t>(
                tile.upper_triangle
                    ? (tile.reference_count > query_id + 1U ? tile.reference_count - query_id - 1U
                                                            : 0U)
                    : tile.reference_count
            );
            if (record + row - chunk.first_record > chunk_records && query != chunk.first) {
                chunk.last = query;
                chunk.end_record = record;
                chunks.push_back(chunk);
                // NVCC 13.4 mislowers assignments that skip aggregate members.
                chunk.first = query;
                chunk.first_record = record;
            }
            record += row;
        }
        chunk.last = tile.query_count;
        chunk.end_record = record;
        if (record != 0U) chunks.push_back(chunk);
    } else if (words != 0U) {
        device_offsets = CUDDL_CUDA_TRY(
            cuda::make_device_buffer<uint32_t>(stream, stream.device(), words, cuda::no_init)
        );
        auto const popcounts =
            cuda::make_transform_iterator(tile.pass_bits, detail::word_popcount{});
        size_t scan_bytes = 0;
        CUDDL_CUDA_TRY(
            cub::DeviceScan::ExclusiveSum(
                nullptr,
                scan_bytes,
                popcounts,
                device_offsets.data(),
                static_cast<int64_t>(words),
                stream.get()
            )
        );
        auto scan_storage = CUDDL_CUDA_TRY(
            cuda::make_device_buffer<uint8_t>(stream, stream.device(), scan_bytes, cuda::no_init)
        );
        CUDDL_CUDA_TRY(
            cub::DeviceScan::ExclusiveSum(
                scan_storage.data(),
                scan_bytes,
                popcounts,
                device_offsets.data(),
                static_cast<int64_t>(words),
                stream.get()
            )
        );
        pass_bits.resize(words);
        offsets.resize(words);
        CUDDL_CUDA_TRY(
            cuda::copy_bytes(
                stream, cuda::std::span{tile.pass_bits, words}, cuda::std::span{pass_bits}
            )
        );
        CUDDL_CUDA_TRY(cuda::copy_bytes(stream, device_offsets, cuda::std::span{offsets}));
        CUDDL_CUDA_TRY(stream.sync());
        auto const total = static_cast<size_t>(offsets.back()) +
                           static_cast<size_t>(cuda::std::popcount(pass_bits.back()));
        for (size_t first = 0U; first < words;) {
            auto const base = static_cast<size_t>(offsets[first]);
            // Whole words only: a word holds at most 32 results.
            auto const last = static_cast<size_t>(
                std::upper_bound(
                    offsets.begin() + static_cast<std::ptrdiff_t>(first) + 1,
                    offsets.end(),
                    static_cast<uint32_t>(base + chunk_records - 32U)
                ) -
                offsets.begin()
            );
            auto const end = last == words ? total : static_cast<size_t>(offsets[last]);
            if (end != base) chunks.push_back({first, last, base, end});
            first = last;
        }
    }
    if (chunks.empty()) {
        CUDDL_CUDA_TRY(stream.sync());
        return Ok();
    }
 
    auto staging = CUDDL_CUDA_TRY(
        cuda::make_device_buffer<packed_pairwise_counts>(
            stream, stream.device(), threshold ? chunk_records : 0U, cuda::no_init
        )
    );
    std::array<pinned_counts, 2> buffers{
        pinned_counts(stream, memory, chunk_records, cuda::no_init),
        pinned_counts(stream, memory, chunk_records, cuda::no_init),
    };
    std::array<cuda::event, 2> ready{cuda::event(stream.device()), cuda::event(stream.device())};
    auto const enqueue = [&](size_t k) -> Result<void> {
        auto const& chunk = chunks[k];
        auto const records = chunk.end_record - chunk.first_record;
        packed_pairwise_counts const* source = tile.results + chunk.first_record;
        if (threshold) {
            auto const kernel = tile.upper_triangle
                                    ? detail::gather_passing_kernel<true, packed_pairwise_counts>
                                    : detail::gather_passing_kernel<false, packed_pairwise_counts>;
            kernel<<<
                detail::warp_grid_blocks(chunk.last - chunk.first),
                detail::block_size,
                0,
                stream.get()>>>(
                tile.results,
                nullptr,
                tile.pass_bits,
                device_offsets.data(),
                tile.first_query_id,
                tile.reference_count,
                chunk.first,
                chunk.last,
                static_cast<uint32_t>(chunk.first_record),
                staging.data(),
                nullptr
            );
            CUDDL_CUDA_TRY(cudaGetLastError());
            source = staging.data();
        }
        CUDDL_CUDA_TRY(
            cuda::copy_bytes(
                stream,
                cuda::std::span{source, records},
                cuda::std::span{buffers[k % 2U].data(), records}
            )
        );
        ready[k % 2U].record(stream);
        return Ok();
    };
 
    CUDDL_TRY(enqueue(0U));
    for (size_t k = 0; k < chunks.size(); ++k) {
        // Buffer (k + 1) % 2 last held chunk k - 1, which the host has finished visiting.
        if (k + 1U < chunks.size()) CUDDL_TRY(enqueue(k + 1U));
        CUDDL_CUDA_TRY(ready[k % 2U].sync());
        auto const* packed = buffers[k % 2U].data();
        auto const& chunk = chunks[k];
        size_t next = 0U;
        auto const emit = [&](uint32_t query_id, uint32_t reference_id) {
            f(batch_search_result{
                .query_id = query_id,
                .reference_id = reference_id,
                .counts = packed[next++].unpack(tile.bucket_count),
            });
        };
        if (!threshold) {
            for (auto query = static_cast<uint32_t>(chunk.first); query < chunk.last; ++query) {
                auto const query_id = tile.first_query_id + query;
                for (auto reference = tile.upper_triangle ? query_id + 1U : 0U;
                     reference < tile.reference_count;
                     ++reference) {
                    emit(query_id, reference);
                }
            }
            continue;
        }
        auto const words_per_query = tile.words_per_query();
        auto query = static_cast<uint32_t>(chunk.first / words_per_query);
        auto column = static_cast<uint32_t>(chunk.first % words_per_query);
        auto w = chunk.first;
        while (w < chunk.last) {
            // A transfer chunk may start or end inside a query row.
            auto const row_end = std::min(chunk.last, w + words_per_query - column);
            auto const query_id = tile.first_query_id + query;
            for (auto base = column * 32U; w < row_end; ++w, base += 32U) {
                for (auto bits = pass_bits[w]; bits != 0U; bits &= bits - 1U) {
                    emit(query_id, base + static_cast<uint32_t>(cuda::std::countr_zero(bits)));
                }
            }
            ++query;
            column = 0U;
        }
    }
    return Ok();
}

Here is the call graph for this function:

gather_sketches()

Result< reference_database_file > cuddl::gather_sketches ( reference_database_file const &  mine, genome_slice const &  slice, MPI_Comm  comm  )

inline

Gathers every rank's sketches of its genome_slice into one file on every rank.

Collective.

mine must hold the sketches of slice.mine(), in order, and every rank must have built with the same k, buckets, layout and blacklist; both are checked. The result equals a single-process build of slice.paths: rows in reference-ID order, labelled with the paths. Rows travel over MPI from host memory; every rank receives all of them.

Definition at line 285 of file distributed.cuh.

                                                                                               {
    auto const counts = CUDDL_TRY(detail::slice_counts(mine, slice, comm));
    std::vector<int> sizes(counts.begin(), counts.end()), offsets(counts.size());
    std::exclusive_scan(sizes.begin(), sizes.end(), offsets.begin(), 0);
    auto const buckets = mine.metadata().compatibility.bucket_count;
    auto row = CUDDL_TRY(detail::row_type(buckets));
    std::vector<uint16_t> rows(slice.paths.size() * buckets);
    auto const gathered = detail::mpi_try(
        MPI_Allgatherv(
            mine.rows().data(),
            static_cast<int>(mine.metadata().reference_count),
            row,
            rows.data(),
            sizes.data(),
            offsets.data(),
            row,
            comm
        ),
        "MPI_Allgatherv"
    );
    MPI_Type_free(&row);
    CUDDL_TRY(gathered);
    return detail::reference_file_assembly::make(mine, slice.paths, std::move(rows));
}

Here is the call graph for this function:

genome_paths()

Result< std::vector< std::filesystem::path > > cuddl::genome_paths ( std::filesystem::path const &  directory )

inline

The FASTA/FASTQ files under directory, searched recursively and sorted by path, which is the order that numbers reference IDs.

Recognised extensions are .fa, .fna, .fasta, .ffn, .frn, .fq and .fastq, in any case, optionally followed by .gz, .bgz or .bgzf. Entries are filtered by the type the directory reports and symlinks are kept without being followed, so listing a directory of symlinked genomes costs no per-file metadata lookup.

Definition at line 214 of file reference_database_file.cuh.

      {
    std::vector<std::filesystem::path> paths;
    for (auto const& entry : std::filesystem::recursive_directory_iterator(
             directory, std::filesystem::directory_options::skip_permission_denied
         )) {
        if (!entry.is_symlink() && !entry.is_regular_file()) continue;
        auto name = entry.path().filename().string();
        std::ranges::transform(name, name.begin(), [](unsigned char c) {
            return static_cast<char>(std::tolower(c));
        });
        auto extension = std::filesystem::path(name).extension().string();
        if (extension == ".gz" || extension == ".bgz" || extension == ".bgzf") {
            extension = std::filesystem::path(name).stem().extension().string();
        }
        if (extension == ".fa" || extension == ".fna" || extension == ".fasta" ||
            extension == ".ffn" || extension == ".frn" || extension == ".fq" ||
            extension == ".fastq") {
            paths.push_back(entry.path());
        }
    }
    std::ranges::sort(paths);
    return paths;
} catch (std::filesystem::filesystem_error const& error) {
    return Err(Error::resource(error.what()));
}

Here is the caller graph for this function:

list_genomes()

Result< std::vector< std::filesystem::path > > cuddl::list_genomes ( std::filesystem::path const &  input )

inline

Paths from a directory (see genome_paths) or from a text file with one path per line, in the file's order.

Definition at line 244 of file reference_database_file.cuh.

  {
    if (std::filesystem::is_directory(input)) return genome_paths(input);
    std::ifstream list(input);
    if (!list) return Err(Error::resource("cannot read genome list " + input.string()));
    std::vector<std::filesystem::path> paths;
    for (std::string line; std::getline(list, line);) {
        if (!line.empty()) paths.emplace_back(line);
    }
    return paths;
}

Here is the call graph for this function:

Here is the caller graph for this function:

Ok()

Result< void > cuddl::Ok ( )

inlinenoexcept

Success return for Result<void>; same as return {}.

Definition at line 181 of file error.hpp.

                                                {
    return {};
}

Here is the caller graph for this function:

parse_fasta_file()

Result< fasta_parse_result > cuddl::parse_fasta_file ( std::string const &  path, uint32_t  k, unsigned  threads = std::thread::hardware_concurrency()  )

inline

Parses a FASTA/FASTQ file's sequence bases into packed canonical k-mers of length k.

Gzip/BGZF is detected by magic bytes and decompressed before parsing. Format is detected from the first record header; FASTQ qualities are validated and ignored. This is the FASTX entry point. It reads every record as one combined genome and emits a packed canonical k-mer per contiguous run of k valid bases; an invalid or ambiguous base breaks the rolling window, and no k-mer spanning it is emitted. Record boundaries also reset the window, without contributing to invalid_windows. See detail::parse_fasta.

Parameters

path	FASTA/FASTQ file to parse.
k	K-mer length.
threads	Worker count for parallel files, defaulting to std::thread::hardware_concurrency(); files at or below 1 MiB parse serially.

Returns

Parsed k-mers and counts, or an error if the file cannot be opened.

Definition at line 34 of file fastx.hpp.

  {
    return detail::parse_fasta(path, k, threads);
}

Here is the call graph for this function:

Here is the caller graph for this function:

search_all_to_all() [1/2]

template<uint32_t K, size_t BucketCount, typename Layout , typename OnTile >

Result< void > cuddl::search_all_to_all	(	replicated_database< K, BucketCount, Layout > &	group,
		OnTile &&	on_tile,
		MPI_Comm	comm
	)

Compares every pair of group's references, sharing the tiles among the ranks of comm; collective.

Every rank calls it with a replicated_database of the same rows. Each rank delivers the tiles it claims to on_tile, as replicated_database::search_all_to_all does, and every pair is delivered on exactly one rank. The ranks share one shared_tile_counter for the call.

Definition at line 475 of file distributed.cuh.

  {
    return detail::shared_search(comm, [&](replica_search split) {
        return group.search_all_to_all(on_tile, std::move(split));
    });
}

Here is the call graph for this function:

Here is the caller graph for this function:

search_all_to_all() [2/2]

template<uint32_t K, size_t BucketCount, typename Layout , typename OnTile >

Result< void > cuddl::search_all_to_all	(	replicated_database< K, BucketCount, Layout > &	group,
		OnTile &&	on_tile,
		search_options	options,
		MPI_Comm	comm
	)

The threshold form of search_all_to_all: only pairs passing options reach on_tile.

group needs its index (replicated_database::build_index); collective.

Definition at line 488 of file distributed.cuh.

  {
    return detail::shared_search(comm, [&](replica_search split) {
        return group.search_all_to_all(on_tile, options, std::move(split));
    });
}

Here is the call graph for this function:

select_top_k()

std::vector< batch_search_result > cuddl::select_top_k ( std::vector< batch_search_result >  rows, uint32_t  k  )

inline

The k best rows of each query among rows, as top_k_search ranks them, best first and query by query.

Merges the downloads of collectors that saw disjoint tiles.

Definition at line 533 of file top_k.cuh.

                                                              {
    std::ranges::sort(rows, [](batch_search_result const& a, batch_search_result const& b) {
        if (a.query_id != b.query_id) return a.query_id < b.query_id;
        return detail::top_k_key(a.counts, a.reference_id) >
               detail::top_k_key(b.counts, b.reference_id);
    });
    std::vector<batch_search_result> kept;
    for (size_t first = 0U; first < rows.size();) {
        auto last = first;
        while (last < rows.size() && rows[last].query_id == rows[first].query_id) ++last;
        kept.insert(
            kept.end(),
            rows.begin() + static_cast<std::ptrdiff_t>(first),
            rows.begin() + static_cast<std::ptrdiff_t>(std::min<size_t>(last, first + k))
        );
        first = last;
    }
    return kept;
}

Here is the call graph for this function:

slice_genomes()

Result< genome_slice > cuddl::slice_genomes ( std::filesystem::path const &  input, MPI_Comm  comm, unsigned  workers = default_parser_workers()  )

inline

Lists a genome collection once and splits it into one contiguous slice per rank of comm, of about equal file bytes.

Collective.

Rank 0 lists input, a directory (see genome_paths) or a text file with one path per line, and shares the list. Every rank then looks up the sizes of an equal share of the files with workers threads, because on a parallel file system each lookup waits for a metadata round trip, and the ranks exchange them. Slices follow path order, so concatenating them in rank order gives the reference IDs of a single-process build. Bytes stand in for sketching work; a rank may get no genomes when there are fewer genomes than ranks.

Definition at line 159 of file distributed.cuh.

  {
    if (workers == 0U)
        return Err(Error::invalid_argument("slice_genomes needs at least one worker"));
    int rank = 0, size = 1;
    CUDDL_TRY(detail::mpi_try(MPI_Comm_rank(comm, &rank), "MPI_Comm_rank"));
    CUDDL_TRY(detail::mpi_try(MPI_Comm_size(comm, &size), "MPI_Comm_size"));
    genome_slice slice;
    slice.paths = CUDDL_TRY(detail::share_listing(input, rank, comm));
    auto const& paths = slice.paths;
    if (paths.size() > static_cast<size_t>(std::numeric_limits<int>::max())) {
        return Err(Error::resource("genome collection exceeds MPI counts"));
    }
 
    auto const share = [&](int r) {
        return paths.size() * static_cast<size_t>(r) / static_cast<size_t>(size);
    };
    auto const own = share(rank);
    std::vector<uint64_t> mine(share(rank + 1) - own);
    std::vector<std::error_code> errors(mine.size());
    {
        std::vector<std::jthread> threads;
        for (unsigned t = 0; t < workers && t < mine.size(); ++t) {
            threads.emplace_back([&, t] {
                for (size_t i = t; i < mine.size(); i += workers) {
                    mine[i] = std::filesystem::file_size(paths[own + i], errors[i]);
                }
            });
        }
    }
    auto sized = [&]() -> Result<void> {
        for (size_t i = 0; i < mine.size(); ++i) {
            if (errors[i]) {
                return Err(Error::resource(paths[own + i].string() + ": " + errors[i].message()));
            }
        }
        return Ok();
    }();
    CUDDL_TRY(detail::agree(std::move(sized), comm));
 
    std::vector<int> counts(size), offsets(size);
    for (int r = 0; r < size; ++r) {
        offsets[r] = static_cast<int>(share(r));
        counts[r] = static_cast<int>(share(r + 1) - share(r));
    }
    slice.bytes.resize(paths.size());
    CUDDL_TRY(
        detail::mpi_try(
            MPI_Allgatherv(
                mine.data(),
                static_cast<int>(mine.size()),
                MPI_UINT64_T,
                slice.bytes.data(),
                counts.data(),
                offsets.data(),
                MPI_UINT64_T,
                comm
            ),
            "MPI_Allgatherv"
        )
    );
    auto const total = std::accumulate(slice.bytes.begin(), slice.bytes.end(), uint64_t{0});
    slice.first = detail::slice_begin(slice.bytes, total, rank, size);
    slice.last =
        rank + 1 == size ? paths.size() : detail::slice_begin(slice.bytes, total, rank + 1, size);
    return slice;
}

Here is the call graph for this function:

weighted_kmer_identity()

__host__ __device__ constexpr float cuddl::weighted_kmer_identity ( pairwise_counts const &  counts )

constexprnoexcept

Weighted k-mer identity of a pair, equal / max(equal + min(lower, higher), 6), as sketch::wkid computes it, or 0 when no bucket is comparable.

The divisor's floor of six suppresses similarity from sparse random collisions.

Definition at line 83 of file pairwise_counts.cuh.

           {
    auto const divisor =
        counts.equal + (counts.lower < counts.higher ? counts.lower : counts.higher);
    if (divisor == 0U) return 0.0F;
    return static_cast<float>(counts.equal) / static_cast<float>(divisor < 6U ? 6U : divisor);
}

Here is the caller graph for this function: