Metadata-Version: 2.2
Name: selinv
Version: 1.0.0
Summary: Fast selected inversion of sparse symmetric matrices (Python wrapper for SelInv)
Keywords: sparse,linear algebra,selected inversion,LDL factorization
Author: Mathieu Bernard
License: BSD-3-Clause
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: BSD License
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Fortran
Classifier: Topic :: Scientific/Engineering :: Mathematics
Project-URL: Repository, https://gitlab.in2p3.fr/lemaitre/selinv
Requires-Python: >=3.10
Requires-Dist: numpy>=1.21
Requires-Dist: scipy>=1.7
Provides-Extra: test
Requires-Dist: pytest>=7; extra == "test"
Provides-Extra: docs
Requires-Dist: sphinx>=7; extra == "docs"
Requires-Dist: sphinx-autodoc-typehints; extra == "docs"
Requires-Dist: myst-parser; extra == "docs"
Requires-Dist: furo; extra == "docs"
Description-Content-Type: text/markdown

# selinv

> [!WARNING]
> The Python bindings and build infrastructure of this package have been
> largely generated with the help of AI (LLM). While the code is tested,
> it may contain inaccuracies or rough edges. Review and contributions are
> welcome.

**selinv** is a Python package for fast *selected inversion* of sparse
symmetric matrices. It wraps the Fortran **SelInv** library by Lin Lin,
Chao Yang, Juan Meza, Jianfeng Lu, Lexing Ying, and Weinan E.

Selected inversion computes only the entries of A⁻¹ that correspond to the
non-zero pattern of the sparse LDL' factor of A (or a superset thereof),
making it far cheaper than a full inversion for large sparse systems.

---

## Features

- Diagonal of A⁻¹ via `selinv()`, or full selected inverse via `selinvfull()`
- Convenience wrapper `fullselinv` for sub-block extraction and diagonal access
- Input via any `scipy.sparse` format (CSC, CSR, COO, …)
- Fortran core with BLAS/LAPACK acceleration
- Optional METIS fill-reducing reordering

---

## Installation

The `selinv` package works on `python>=3.11` for Linux and macos. It should work on Windows too, but it has not been tested.


### Install from Pypi

If using [pixi](https://pixi.sh/latest/), add `selinv` to your `pixi.toml` file:

```bash
pixi add --pypi selinv
```

Or with `pip` (similarly with `pipx` ou `uv`):

```bash
pip install selinv
```


### Install from source

#### Using pixi (recommanded)

```bash
# As of pixi-0.66 you must add `preview = ["pixi-build"]`
# to the `workspace` or `project` table of your pixi.toml.
pixi add --git https://gitlab.in2p3.fr/lemaitre/selinv.git selinv --branch main
```

#### Using pip

You need a C, C++, and Fortran compiler, CMake ≥ 3.18, and BLAS/LAPACK
development libraries installed on your system.

**System dependencies (Linux):**

```sh
# Debian / Ubuntu
sudo apt-get install gfortran g++ cmake liblapack-dev libblas-dev
# RHEL / Fedora
sudo dnf install gcc-gfortran gcc-c++ cmake lapack-devel blas-devel
```

**System dependencies (macOS):**

```sh
xcode-select --install
brew install gcc cmake    # provides gfortran via the gcc formula
```

On macOS the Accelerate framework is used automatically for BLAS/LAPACK.

**Install the package:**

```sh
git clone https://gitlab.in2p3.fr/lemaitre/selinv.git
cd selinv
pip install .
```

---

## Quick Start

```python
import numpy as np
import scipy.sparse as sp
from selinv import selinv, selinvfull, fullselinv

# Build a small sparse SPD tridiagonal matrix (CSC format)
diag = np.array([4.0, 4.0, 4.0, 4.0])
off = np.array([-1.0, -1.0, -1.0])
A = sp.diags([off, diag, off], [-1, 0, 1], format="csc")

# 1. Compute the diagonal of A^{-1}
d = selinv(A)
print("Diagonal of A^{-1}:", d)

# 2. Compute the full selected inverse with an explicit permutation
P = np.arange(A.shape[0], dtype=np.int32)   # identity permutation
L, Ainv, perm = selinvfull(A, P)

# 3. Use the fullselinv wrapper for convenient access
result = fullselinv(L, Ainv, perm)
print("Diagonal:", result.diag())           # same as selinv(A)
print("Sub-block [[0,1]]:\n", result[[0, 1]])  # 2x2 dense block of A^{-1}
```

---

## API Overview

### `selinv(mat, P=None) -> numpy.ndarray`

Compute the diagonal of A⁻¹ via the SelInv algorithm.

| Parameter | Type | Description |
|-----------|------|-------------|
| `mat` | `scipy.sparse` matrix | Sparse symmetric matrix (lower triangle stored) |
| `P` | `numpy.ndarray` or `None` | Optional 0-based permutation vector. When `None`, Multiple Minimum Degree (MMD) reordering is used. |

**Returns** a 1-D `numpy.ndarray` of shape `(n,)` with dtype `float64`
containing the diagonal elements of A⁻¹ in the original ordering.

---

### `selinvfull(mat, P) -> tuple[csc_array, csc_array, ndarray]`

Compute the full selected inverse: all entries of A⁻¹ at positions where
the LDL' factor L has a non-zero.

| Parameter | Type | Description |
|-----------|------|-------------|
| `mat` | `scipy.sparse` matrix | Sparse symmetric matrix (lower triangle stored) |
| `P` | `numpy.ndarray` | 0-based permutation vector (required) |

**Returns** a 3-tuple:

| Element | Type | Description |
|---------|------|-------------|
| `L` | `scipy.sparse.csc_array` | Lower-triangular factor (diagonal zeroed) |
| `Ainv` | `scipy.sparse.csc_array` | Selected inverse in the permuted basis |
| `perm` | `numpy.ndarray` (int32) | Inverse permutation vector: `perm[i]` gives the permuted position of original index `i`. May differ from the simple inverse of `P` due to elimination-tree postordering. |

---

### `fullselinv(L, Ainv, perm)`

Convenience class wrapping the output of `selinvfull()`.

| Method / Accessor | Description |
|-------------------|-------------|
| `result.diag()` | Diagonal of A⁻¹ in the original ordering |
| `result[[i, j, ...]]` | Dense sub-block of A⁻¹ for the given row/column indices |
| `result.L` | The lower-triangular factor |
| `result.Ai` | The selected inverse (sparse) |
| `result.perm` | The inverse permutation |

> **Note:** `result[[i, j]]` only returns accurate values for index pairs
> that fall within the fill pattern of L. Entries outside the fill pattern
> will be zero, not the true inverse value. See the docstring for details.

---

## Mathematical Background

Given a sparse symmetric matrix **A**, selected inversion computes the
entries of **A⁻¹** that correspond to the non-zero pattern of the LDL'
factor of **A**. This is significantly cheaper than computing the full
inverse and is exact (not approximate).

The algorithm exploits a recurrence first described by Takahashi, Fagan,
and Chin (1973) and extended to a supernodal setting by Lin et al. (2011).
The computational cost is O(|L|), the same as the factorization itself.

Applications include:

- Green's function calculations in computational physics
- Density matrix evaluation in linear-scaling DFT
- Marginal variances in Gaussian Markov Random Fields (Bayesian inference)
- Uncertainty quantification in large-scale inverse problems

See [`docs/math.md`](docs/math.md) for a detailed derivation.

---

## Building the Documentation

The project includes Sphinx documentation under `docs/`. To build it
locally:

```sh
pixi run -e dev docs
```

Then open `docs/_build/html/index.html` in your browser.

To clean the build:

```sh
pixi run -e dev docs-clean
```

---

## Running the Tests

```sh
pixi run test
```

Or for a specific Python version:

```sh
pixi run -e py314 test
```

---

## Releasing a New Version

> This section is intended for maintainers.

Releases are driven by git tags: pushing a tag matching `vX.Y.Z` triggers
the CI/CD pipeline which builds wheels for all supported Python versions and
publishes them to PyPI automatically.

Version bumping is automated via [`bump-my-version`](https://github.com/callowayproject/bump-my-version),
available in the `dev` environment.

**Step-by-step:**

1. **Update `CHANGELOG.md`** — move the items under `[Unreleased]` to a new
   section named after the upcoming version (e.g. `[0.2.0] - 2025-01-31`).
   Commit this change manually if you want it separate from the version bump.

2. **Bump the version** — choose the appropriate increment:

   ```sh
   pixi run -e dev bump-patch   # 0.1.0 → 0.1.1  (bug fixes)
   pixi run -e dev bump-minor   # 0.1.0 → 0.2.0  (new features)
   pixi run -e dev bump-major   # 0.1.0 → 1.0.0  (breaking changes)
   ```

   This updates both `version` fields in `pyproject.toml`, creates a commit
   (`chore: release X.Y.Z`), and tags it `vX.Y.Z`.

3. **Push the commit and the tag:**

   ```sh
   git push --follow-tags
   ```

The CI/CD pipeline then takes over: it builds the source distribution and
binary wheels, and uploads everything to PyPI via `twine`.

---

## License

This Python wrapper is released under the BSD 3-Clause License.
Copyright © 2025 Mathieu Bernard.

The original Fortran SelInv library (`extern/selinv/`) is distributed under
its own BSD-style license.
Copyright © 2011 The Regents of the University of California.

See [LICENSE](LICENSE) for the full text of both licenses.

---

## Citation

If you use **selinv** in academic work, please cite the original SelInv paper:

> Lin Lin, Chao Yang, Juan C. Meza, Jianfeng Lu, Lexing Ying, and Weinan E.
> **SelInv — An Algorithm for Selected Inversion of a Sparse Symmetric Matrix.**
> *ACM Transactions on Mathematical Software*, **37**(4), Article 40, 2011.
> https://doi.org/10.1145/1916461.1916464

<details>
<summary>BibTeX</summary>

```bibtex
@article{lin2011selinv,
  author  = {Lin, Lin and Yang, Chao and Meza, Juan C. and Lu, Jianfeng
             and Ying, Lexing and E, Weinan},
  title   = {{SelInv} --- An Algorithm for Selected Inversion of a
             Sparse Symmetric Matrix},
  journal = {ACM Transactions on Mathematical Software},
  volume  = {37},
  number  = {4},
  pages   = {40:1--40:19},
  year    = {2011},
  doi     = {10.1145/1916461.1916464},
}
```

</details>
