### Install Serinv on Alex Environment (GPU)

Source: https://github.com/vincent-maillou/serinv/blob/main/README.md

Installation steps for the Alex environment on Fau cluster, including MPI, GPU support, and Python package installation.

```bash
# --- Alex-env ---
module load python
module load openmpi/4.1.6-nvhpc23.7-cuda
module load cuda/12.6.1

conda create -n alex
conda activate alex

CFLAGS=-noswitcherror MPICC=$(which mpicc) pip install --no-cache-dir mpi4py

salloc --partition=a40 --nodes=1 --gres=gpu:a40:1 --time 01:00:00
conda activate alex

conda install -c conda-forge cupy-core
conda install blas=*=*mkl
conda install libblas=*=*mkl
conda install numpy scipy
conda install -c conda-forge pytest pytest-mpi pytest-cov coverage black isort ruff just pre-commit matplotlib numba -y
cd /path/to/serinv/
python -m pip install -e .

```

--------------------------------

### Install Serinv on Fritz Environment (CPU)

Source: https://github.com/vincent-maillou/serinv/blob/main/README.md

Installation steps for the Fritz environment on Fau cluster, including MPI support and Python package installation.

```bash
# --- Fritz-env ---
module load python
module load openmpi/4.1.2-gcc11.2.0

conda create -n fritz
conda activate fritz

MPICC=$(which mpicc) pip install --no-cache-dir mpi4py

salloc -N 4 --time 01:00:00
conda activate fritz

conda install blas=*=*mkl
conda install libblas=*=*mkl
conda install numpy scipy
conda install -c conda-forge pytest pytest-mpi pytest-cov coverage black isort ruff just pre-commit matplotlib numba -y
cd /path/to/serinv/
python -m pip install -e .

```

--------------------------------

### Install Miniconda

Source: https://github.com/vincent-maillou/serinv/blob/main/README.md

Downloads and makes the Miniconda installer executable. This script is used to install Miniconda on the system.

```bash
wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-aarch64.sh
chmod u+x Miniconda3-latest-Linux-aarch64.sh
./Miniconda3-latest-Linux-aarch64.sh
```

--------------------------------

### Set up Bare Environment

Source: https://github.com/vincent-maillou/serinv/blob/main/README.md

Initializes a Conda environment named 'bare' and installs essential scientific computing libraries.

```bash
module load python
conda create -n bare
salloc -N 4 --time 01:00:00
conda activate bare

conda install blas=*=*mkl
conda install libblas=*=*mkl
conda install numpy scipy
conda install -c conda-forge pytest pytest-mpi pytest-cov coverage black isort ruff just pre-commit matplotlib numba -y
cd /path/to/serinv/
python -m pip install -e .
```

--------------------------------

### Install Serinv and Dependencies

Source: https://context7.com/vincent-maillou/serinv/llms.txt

Commands to set up a conda environment, install the Serinv package, and optionally install mpi4py for distributed computing or cupy for GPU support.

```bash
conda create --name serinv_env python=3.11
conda activate serinv_env

cd /path/to/serinv/
pip install -e .

pip install mpi4py

pip install cupy
```

--------------------------------

### Install Serinv and Run Tests

Source: https://github.com/vincent-maillou/serinv/blob/main/README.md

Installs the Serinv project in editable mode and runs the project's sequential tests using pytest.

```bash
cd /path/to/serinv/
python -m pip install -e .

# Run the sequential tests.
pytest .
```

--------------------------------

### Pull and Start uenv Image

Source: https://github.com/vincent-maillou/serinv/blob/main/README.md

Finds available `uenv` images, pulls a specific version, and starts it with module view enabled.

```bash
uenv image find
uenv image pull prgenv-gnu/24.11:v1
uenv start --view=modules prgenv-gnu/24.11:v1
```

--------------------------------

### Create Conda Environment and Install Libraries

Source: https://github.com/vincent-maillou/serinv/blob/main/README.md

Creates a new Conda environment named 'myenv', activates it, and installs Python 3.12, NumPy, SciPy, mpi4py, Cupy, and development tools. It also includes a test for NCCL/CuPy installation.

```bash
conda create -n myenv
conda activate myenv

conda install python=3.12
conda install numpy scipy

MPICC=$(which mpicc) python -m pip install --no-cache-dir mpi4py
pip install cupy --no-dependencies --no-cache-dir

conda install -c conda-forge pytest pytest-mpi pytest-cov coverage black isort ruff just pre-commit matplotlib numba -y

# Test the NCCL/CuPy installation
python -c "from cupy.cuda.nccl import *"
```

--------------------------------

### Block Tridiagonal Factorization and Solvers

Source: https://context7.com/vincent-maillou/serinv/llms.txt

Demonstrates Cholesky factorization (pobtf) and solving linear systems Lx=b (pobts) and L^T x = y using forward and backward substitution for block tridiagonal matrices.

```python
pobtf(A_diagonal_blocks, A_lower_diagonal_blocks)

pobts(A_diagonal_blocks, A_lower_diagonal_blocks, B, trans="N")

pobts(A_diagonal_blocks, A_lower_diagonal_blocks, B, trans="T")

print("Solution x:", B[:8])
```

--------------------------------

### Backend Availability Check

Source: https://context7.com/vincent-maillou/serinv/llms.txt

Check which computational backends (CuPy, NCCL, MPI) are available at runtime. This allows for writing portable code that adapts to the available hardware.

```python
from serinv import backend_flags

# Check available backends
print("CuPy available:", backend_flags["cupy_avail"])
print("NCCL available:", backend_flags["nccl_avail"])
print("MPI available:", backend_flags["mpi_avail"])
print("MPI CUDA-aware:", backend_flags["mpi_cuda_aware"])

# Conditional execution based on backend
if backend_flags["cupy_avail"]:
    import cupy as cp
    # Use GPU arrays
    A = cp.random.rand(100, 100)
else:
    import numpy as np
    # Fall back to CPU
    A = np.random.rand(100, 100)
```

--------------------------------

### Load Necessary Modules

Source: https://github.com/vincent-maillou/serinv/blob/main/README.md

Loads essential modules required for the project, including CUDA, compilers, build tools, and MPI libraries.

```bash
module load cuda
module load gcc
module load meson
module load ninja
module load nccl
module load cray-mpich
module load cmake
module load openblas
module load aws-ofi-nccl
```

--------------------------------

### Export Library Paths

Source: https://github.com/vincent-maillou/serinv/blob/main/README.md

Configures environment variables to point to the root, library, and include directories of NCCL and CUDA. This is crucial for the system to locate these libraries during runtime and compilation.

```bash
export NCCL_ROOT=/user-environment/linux-sles15-neoverse_v2/gcc-13.3.0/nccl-2.22.3-1-4j6h3ffzysukqpqbvriorrzk2lm762dd
export NCCL_LIB_DIR=$NCCL_ROOT/lib
export NCCL_INCLUDE_DIR=$NCCL_ROOT/include
export CUDA_DIR=$CUDA_HOME
export CUDA_PATH=$CUDA_HOME
export CPATH=$CUDA_HOME/include:$CPATH
export LIBRARY_PATH=$CUDA_HOME/lib64:$LIBRARY_PATH
export LD_LIBRARY_PATH=$CUDA_HOME/lib64:$LD_LIBRARY_PATH
export CPATH=$NCCL_ROOT/include:$CPATH
export LIBRARY_PATH=$NCCL_ROOT/lib:$LIBRARY_PATH
export LD_LIBRARY_PATH=$NCCL_ROOT/lib:$LD_LIBRARY_PATH
```

--------------------------------

### Direct GPU Computation with CuPy

Source: https://context7.com/vincent-maillou/serinv/llms.txt

Perform matrix factorization and selected inversion entirely on the GPU by passing CuPy arrays directly to the algorithms. Ensure CuPy is available before execution.

```python
import numpy as np
from serinv import backend_flags
from serinv.algs import pobtf, pobtsi

if backend_flags["cupy_avail"]:
    import cupy as cp

    diagonal_blocksize = 4
    n_diag_blocks = 5

    # Create SPD matrix directly on GPU
    A_diagonal_blocks = cp.zeros((n_diag_blocks, diagonal_blocksize, diagonal_blocksize))
    A_lower_diagonal_blocks = cp.zeros((n_diag_blocks - 1, diagonal_blocksize, diagonal_blocksize))

    for i in range(n_diag_blocks):
        temp = cp.random.rand(diagonal_blocksize, diagonal_blocksize)
        A_diagonal_blocks[i] = temp @ temp.T + 10 * cp.eye(diagonal_blocksize)
        if i < n_diag_blocks - 1:
            A_lower_diagonal_blocks[i] = 0.1 * cp.random.rand(diagonal_blocksize, diagonal_blocksize)

    # Factorization runs entirely on GPU
    pobtf(A_diagonal_blocks, A_lower_diagonal_blocks)

    # Selected inversion on GPU
    pobtsi(A_diagonal_blocks, A_lower_diagonal_blocks)

    # Transfer results back to host if needed
    result_host = cp.asnumpy(A_diagonal_blocks)
    print("Inverse diagonal block 0:\n", result_host[0])
```

--------------------------------

### GPU Acceleration with Device Streaming

Source: https://context7.com/vincent-maillou/serinv/llms.txt

Utilizes pinned memory and device streaming to perform factorization and inversion on matrices that exceed GPU memory capacity.

```python
import numpy as np
from serinv import backend_flags
from serinv.algs import pobtf, pobtsi

if backend_flags["cupy_avail"]:
    import cupyx as cpx

    diagonal_blocksize = 4
    n_diag_blocks = 5

    # Create SPD matrix on host
    A_diagonal_blocks = np.zeros((n_diag_blocks, diagonal_blocksize, diagonal_blocksize))
    A_lower_diagonal_blocks = np.zeros((n_diag_blocks - 1, diagonal_blocksize, diagonal_blocksize))

    for i in range(n_diag_blocks):
        temp = np.random.rand(diagonal_blocksize, diagonal_blocksize)
        A_diagonal_blocks[i] = temp @ temp.T + 10 * np.eye(diagonal_blocksize)
        if i < n_diag_blocks - 1:
            A_lower_diagonal_blocks[i] = 0.1 * np.random.rand(diagonal_blocksize, diagonal_blocksize)

    # Allocate pinned memory for efficient host-device transfer
    A_diagonal_blocks_pinned = cpx.zeros_like_pinned(A_diagonal_blocks)
    A_diagonal_blocks_pinned[:] = A_diagonal_blocks
    A_lower_diagonal_blocks_pinned = cpx.zeros_like_pinned(A_lower_diagonal_blocks)
    A_lower_diagonal_blocks_pinned[:] = A_lower_diagonal_blocks

    # Factorization with GPU streaming
    pobtf(
        A_diagonal_blocks_pinned,
        A_lower_diagonal_blocks_pinned,
        device_streaming=True
    )

    # Selected inversion with GPU streaming
    pobtsi(
        A_diagonal_blocks_pinned,
        A_lower_diagonal_blocks_pinned,
        device_streaming=True
    )

    print("Result computed with GPU streaming")
```

--------------------------------

### Block Tridiagonal Factorization (pobtf)

Source: https://context7.com/vincent-maillou/serinv/llms.txt

Performs Cholesky factorization of a symmetric positive definite block tridiagonal matrix in-place. Ensure input blocks are initialized to guarantee positive definiteness.

```python
import numpy as np
from serinv.algs import pobtf

diagonal_blocksize = 4
n_diag_blocks = 5

A_diagonal_blocks = np.zeros((n_diag_blocks, diagonal_blocksize, diagonal_blocksize))
A_lower_diagonal_blocks = np.zeros((n_diag_blocks - 1, diagonal_blocksize, diagonal_blocksize))

for i in range(n_diag_blocks):
    temp = np.random.rand(diagonal_blocksize, diagonal_blocksize)
    A_diagonal_blocks[i] = temp @ temp.T + 10 * np.eye(diagonal_blocksize)  # Ensure SPD
    if i < n_diag_blocks - 1:
        A_lower_diagonal_blocks[i] = 0.1 * np.random.rand(diagonal_blocksize, diagonal_blocksize)

pobtf(A_diagonal_blocks, A_lower_diagonal_blocks)

print("L diagonal block 0:\n", A_diagonal_blocks[0])
```

--------------------------------

### Block Tridiagonal Arrowhead Selected Inversion (pobtasi)

Source: https://context7.com/vincent-maillou/serinv/llms.txt

Computes the selected inversion of a block tridiagonal arrowhead matrix after its Cholesky factorization. The input arrays are modified in-place to store the inverse blocks.

```python
import numpy as np
from serinv.algs import pobtaf, pobtasi

diagonal_blocksize = 4
arrowhead_blocksize = 2
n_diag_blocks = 5

# Initialize arrays
A_diagonal_blocks = np.zeros((n_diag_blocks, diagonal_blocksize, diagonal_blocksize))
A_lower_diagonal_blocks = np.zeros((n_diag_blocks - 1, diagonal_blocksize, diagonal_blocksize))
A_lower_arrow_blocks = np.zeros((n_diag_blocks, arrowhead_blocksize, diagonal_blocksize))
A_arrow_tip_block = np.zeros((arrowhead_blocksize, arrowhead_blocksize))

# Create SPD block tridiagonal arrowhead matrix
for i in range(n_diag_blocks):
    temp = np.random.rand(diagonal_blocksize, diagonal_blocksize)
    A_diagonal_blocks[i] = temp @ temp.T + 15 * np.eye(diagonal_blocksize)
    if i < n_diag_blocks - 1:
        A_lower_diagonal_blocks[i] = 0.1 * np.random.rand(diagonal_blocksize, diagonal_blocksize)
    A_lower_arrow_blocks[i] = 0.05 * np.random.rand(arrowhead_blocksize, diagonal_blocksize)

temp_tip = np.random.rand(arrowhead_blocksize, arrowhead_blocksize)
A_arrow_tip_block[:] = temp_tip @ temp_tip.T + 20 * np.eye(arrowhead_blocksize)

# Step 1: Cholesky factorization
pobtaf(A_diagonal_blocks, A_lower_diagonal_blocks, A_lower_arrow_blocks, A_arrow_tip_block)

# Step 2: Selected inversion
pobtasi(A_diagonal_blocks, A_lower_diagonal_blocks, A_lower_arrow_blocks, A_arrow_tip_block)

# Arrays now contain selected inverse blocks
print("Inverse diagonal block 0:\n", A_diagonal_blocks[0])
print("Inverse arrow tip block:\n", A_arrow_tip_block)
```

--------------------------------

### Distributed Parallel Factorization with MPI

Source: https://context7.com/vincent-maillou/serinv/llms.txt

Execute matrix factorization in parallel across multiple processes using MPI for distributed memory systems. Requires MPI to be available and configured.

```python
from mpi4py import MPI
import numpy as np
from serinv import backend_flags
from serinv.wrappers import ppobtf, allocate_pobtrs

if backend_flags["mpi_avail"]:
    comm = MPI.COMM_WORLD
    rank = comm.Get_rank()
    size = comm.Get_size()

    diagonal_blocksize = 4
    n_diag_blocks_per_process = 3

    # Each process owns a portion of the matrix
    A_diagonal_blocks_local = np.zeros(
        (n_diag_blocks_per_process, diagonal_blocksize, diagonal_blocksize)
    )
    A_lower_diagonal_blocks_local = np.zeros(
        (n_diag_blocks_per_process - 1, diagonal_blocksize, diagonal_blocksize)
    )

    # Initialize local blocks with SPD data
    for i in range(n_diag_blocks_per_process):
        temp = np.random.rand(diagonal_blocksize, diagonal_blocksize)
        A_diagonal_blocks_local[i] = temp @ temp.T + 10 * np.eye(diagonal_blocksize)
        if i < n_diag_blocks_per_process - 1:
            A_lower_diagonal_blocks_local[i] = 0.1 * np.random.rand(diagonal_blocksize, diagonal_blocksize)

    # Allocate reduced system arrays for parallel algorithm
    reduced_system = allocate_pobtrs(A_diagonal_blocks_local, comm)

    # Perform parallel factorization
    ppobtf(
        A_diagonal_blocks_local,
        A_lower_diagonal_blocks_local,
        reduced_system=reduced_system,
        comm=comm
    )

    if rank == 0:
        print("Parallel factorization complete")
```

--------------------------------

### Block Tridiagonal Selected Inversion (pobtsi)

Source: https://context7.com/vincent-maillou/serinv/llms.txt

Computes the selected inversion of a block tridiagonal matrix after its Cholesky factorization. The operation is performed in-place, overwriting the factorization results with the inverse blocks.

```python
import numpy as np
from serinv.algs import pobtf, pobtsi

diagonal_blocksize = 4
n_diag_blocks = 5

A_diagonal_blocks = np.zeros((n_diag_blocks, diagonal_blocksize, diagonal_blocksize))
A_lower_diagonal_blocks = np.zeros((n_diag_blocks - 1, diagonal_blocksize, diagonal_blocksize))

for i in range(n_diag_blocks):
    temp = np.random.rand(diagonal_blocksize, diagonal_blocksize)
    A_diagonal_blocks[i] = temp @ temp.T + 10 * np.eye(diagonal_blocksize)
    if i < n_diag_blocks - 1:
        A_lower_diagonal_blocks[i] = 0.1 * np.random.rand(diagonal_blocksize, diagonal_blocksize)

pobtf(A_diagonal_blocks, A_lower_diagonal_blocks)

pobtsi(A_diagonal_blocks, A_lower_diagonal_blocks)

print("Inverse diagonal block 0:\n", A_diagonal_blocks[0])
```

--------------------------------

### Block Tridiagonal Arrowhead Cholesky Factorization (pobtaf)

Source: https://context7.com/vincent-maillou/serinv/llms.txt

Performs Cholesky factorization for a symmetric positive definite block tridiagonal matrix with an arrowhead structure. Ensure matrix blocks are initialized correctly before calling.

```python
import numpy as np
from serinv.algs import pobtaf

diagonal_blocksize = 4
arrowhead_blocksize = 2
n_diag_blocks = 5

# Initialize block arrays
A_diagonal_blocks = np.zeros((n_diag_blocks, diagonal_blocksize, diagonal_blocksize))
A_lower_diagonal_blocks = np.zeros((n_diag_blocks - 1, diagonal_blocksize, diagonal_blocksize))
A_lower_arrow_blocks = np.zeros((n_diag_blocks, arrowhead_blocksize, diagonal_blocksize))
A_arrow_tip_block = np.zeros((arrowhead_blocksize, arrowhead_blocksize))

# Create SPD matrix blocks
for i in range(n_diag_blocks):
    temp = np.random.rand(diagonal_blocksize, diagonal_blocksize)
    A_diagonal_blocks[i] = temp @ temp.T + 15 * np.eye(diagonal_blocksize)
    if i < n_diag_blocks - 1:
        A_lower_diagonal_blocks[i] = 0.1 * np.random.rand(diagonal_blocksize, diagonal_blocksize)
    A_lower_arrow_blocks[i] = 0.05 * np.random.rand(arrowhead_blocksize, diagonal_blocksize)

temp_tip = np.random.rand(arrowhead_blocksize, arrowhead_blocksize)
A_arrow_tip_block[:] = temp_tip @ temp_tip.T + 20 * np.eye(arrowhead_blocksize)

# Perform Cholesky factorization (in-place)
pobtaf(
    A_diagonal_blocks,
    A_lower_diagonal_blocks,
    A_lower_arrow_blocks,
    A_arrow_tip_block
)

print("L diagonal block 0:\n", A_diagonal_blocks[0])
print("L arrow tip block:\n", A_arrow_tip_block)
```

--------------------------------

### Block Tridiagonal Solve (pobts)

Source: https://context7.com/vincent-maillou/serinv/llms.txt

Solves a linear system Ax = b using the Cholesky factorization of A. The `trans` parameter controls forward or backward substitution. The solution overwrites the input vector b.

```python
import numpy as np
from serinv.algs import pobtf, pobts

diagonal_blocksize = 4
n_diag_blocks = 5
total_size = diagonal_blocksize * n_diag_blocks

A_diagonal_blocks = np.zeros((n_diag_blocks, diagonal_blocksize, diagonal_blocksize))
A_lower_diagonal_blocks = np.zeros((n_diag_blocks - 1, diagonal_blocksize, diagonal_blocksize))

for i in range(n_diag_blocks):
    temp = np.random.rand(diagonal_blocksize, diagonal_blocksize)
    A_diagonal_blocks[i] = temp @ temp.T + 10 * np.eye(diagonal_blocksize)
    if i < n_diag_blocks - 1:
        A_lower_diagonal_blocks[i] = 0.1 * np.random.rand(diagonal_blocksize, diagonal_blocksize)

B = np.random.rand(total_size)
B_original = B.copy()
```

--------------------------------

### Perform In-place Schur Complement

Source: https://context7.com/vincent-maillou/serinv/llms.txt

Executes the Schur complement operation directly on the provided diagonal and off-diagonal blocks.

```python
ddbtsc(
    A_diagonal_blocks,
    A_lower_diagonal_blocks,
    A_upper_diagonal_blocks
)

# A_diagonal_blocks now contains the Schur complement reduced blocks
print("Schur complement diagonal block 0:\n", A_diagonal_blocks[0])
```

--------------------------------

### General Matrix Schur Complement (ddbtsc)

Source: https://context7.com/vincent-maillou/serinv/llms.txt

Computes the Schur complement for diagonally dominant block tridiagonal matrices. This function can handle standard Schur complements and quadratic forms with optional right-hand side transformations.

```python
import numpy as np
from serinv.algs import ddbtsc

diagonal_blocksize = 4
n_diag_blocks = 5

# Create diagonally dominant block tridiagonal matrix
A_diagonal_blocks = np.zeros((n_diag_blocks, diagonal_blocksize, diagonal_blocksize))
A_lower_diagonal_blocks = np.zeros((n_diag_blocks - 1, diagonal_blocksize, diagonal_blocksize))
A_upper_diagonal_blocks = np.zeros((n_diag_blocks - 1, diagonal_blocksize, diagonal_blocksize))

for i in range(n_diag_blocks):
    A_diagonal_blocks[i] = np.random.rand(diagonal_blocksize, diagonal_blocksize)
    A_diagonal_blocks[i] = A_diagonal_blocks[i] + 10 * np.eye(diagonal_blocksize)  # Make diagonally dominant
    if i < n_diag_blocks - 1:
        A_lower_diagonal_blocks[i] = 0.1 * np.random.rand(diagonal_blocksize, diagonal_blocksize)
        A_upper_diagonal_blocks[i] = 0.1 * np.random.rand(diagonal_blocksize, diagonal_blocksize)
```

--------------------------------

### Schur Complement with Quadratic RHS

Source: https://context7.com/vincent-maillou/serinv/llms.txt

Computes the Schur complement while simultaneously transforming a right-hand side matrix for covariance propagation.

```python
import numpy as np
from serinv.algs import ddbtsc

diagonal_blocksize = 4
n_diag_blocks = 5

# Create diagonally dominant matrix A
A_diagonal_blocks = np.zeros((n_diag_blocks, diagonal_blocksize, diagonal_blocksize))
A_lower_diagonal_blocks = np.zeros((n_diag_blocks - 1, diagonal_blocksize, diagonal_blocksize))
A_upper_diagonal_blocks = np.zeros((n_diag_blocks - 1, diagonal_blocksize, diagonal_blocksize))

# Create RHS matrix B (same structure)
B_diagonal_blocks = np.zeros((n_diag_blocks, diagonal_blocksize, diagonal_blocksize))
B_lower_diagonal_blocks = np.zeros((n_diag_blocks - 1, diagonal_blocksize, diagonal_blocksize))
B_upper_diagonal_blocks = np.zeros((n_diag_blocks - 1, diagonal_blocksize, diagonal_blocksize))

for i in range(n_diag_blocks):
    A_diagonal_blocks[i] = np.eye(diagonal_blocksize) * 10 + np.random.rand(diagonal_blocksize, diagonal_blocksize)
    B_diagonal_blocks[i] = np.eye(diagonal_blocksize) + 0.1 * np.random.rand(diagonal_blocksize, diagonal_blocksize)
    if i < n_diag_blocks - 1:
        A_lower_diagonal_blocks[i] = 0.1 * np.random.rand(diagonal_blocksize, diagonal_blocksize)
        A_upper_diagonal_blocks[i] = 0.1 * np.random.rand(diagonal_blocksize, diagonal_blocksize)
        B_lower_diagonal_blocks[i] = 0.05 * np.random.rand(diagonal_blocksize, diagonal_blocksize)
        B_upper_diagonal_blocks[i] = 0.05 * np.random.rand(diagonal_blocksize, diagonal_blocksize)

# Schur complement with quadratic RHS transformation
rhs = {
    "B_diagonal_blocks": B_diagonal_blocks,
    "B_lower_diagonal_blocks": B_lower_diagonal_blocks,
    "B_upper_diagonal_blocks": B_upper_diagonal_blocks,
}

ddbtsc(
    A_diagonal_blocks,
    A_lower_diagonal_blocks,
    A_upper_diagonal_blocks,
    rhs=rhs,
    quadratic=True
)

print("Transformed B diagonal block 0:\n", B_diagonal_blocks[0])
```

=== COMPLETE CONTENT === This response contains all available snippets from this library. No additional content exists. Do not make further requests.