---
title: "Multi-threading Support"
description: >
  Advanced setting to set up the threads for faster processing
  in Sinkhorn, Barycenter, WDL, and WIG models.
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Multi-threading Support}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
library(rwig) |> suppressPackageStartupMessages()
```

The log methods for `sinkhorn()` and `barycenter()`
both require row-by-row and column-by-column *soft-minimum* operations
for each iteration of the algorithm,
and therefore suffer from slow computation time compared to the 
vanilla/parallel version.

If the dimensions of $M$ and $N$ are large 
(recall that cost matrix **C** is of size $M \times N$),
we can use multi-threading to process the rows and columns simultaneously.

This can be done by setting `n_threads` to an integer bigger than 0.
By default it is 0, and it means that threading is disabled.
Also, setting `n_threads` for "vanilla" or "parallel" methods
will be ignored automatically.

But you might ask: if multi-threading is so wonderful,
why don't you set threading as the default?
This is because threading comes with an overhead,
and sometimes for small problems, it can even be slower
than serial processing.
So be sure to benchmark your code and see if threading actually helps.

## BLAS threads

`rwig` calls the BLAS library that R is linked against for its matrix
products. Some BLAS libraries (OpenBLAS, MKL) run multi-threaded by default.
Because `rwig` already parallelizes the expensive log-domain kernels through
`n_threads`, and because most of its BLAS calls are small, letting the BLAS
spawn its own threads on top usually oversubscribes the cores and slows
things down. We therefore recommend a single BLAS thread while using `rwig`.

`rwig` does not change this setting for you (it is a session-wide setting that
also affects every other package). If the optional package `RhpcBLASctl` is
installed, the startup message reports the current BLAS thread count, and you
can set it for the session with

```{r, eval = FALSE}
RhpcBLASctl::blas_set_num_threads(1)
```

The exception is `wdl()`/`wig()` on the CPU with an optimized BLAS: their cost
is dominated by large matrix products, so a multi-threaded BLAS can help
there. Benchmark both settings on your own problem sizes.

## See Also

See also `vignette("sinkhorn")`, `vignette("barycenter")`.

## References

Peyré, G., & Cuturi, M. (2019). 
Computational Optimal Transport: With Applications to Data Science.
*Foundations and Trends® in Machine Learning*, 11(5–6), 355–607.
https://doi.org/10.1561/2200000073

Xie, F. (2025). 
Deriving the Gradients of Some Popular Optimal Transport Algorithms (No. arXiv:2504.08722). *arXiv*. 
https://doi.org/10.48550/arXiv.2504.08722