### Create and Manipulate Matrix Views in Rust

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Shows how to create immutable and mutable views of matrices using `MatRef` and `MatMut`. Includes examples of slicing, element access, transpose, and splitting.

```rust
use faer::{Mat, mat};

let mut m = mat![
    [1.0, 2.0, 3.0],
    [4.0, 5.0, 6.0],
    [7.0, 8.0, 9.0_f64],
];

// immutable slicing: rows 0..2, cols 1..3  → 2×2 view
let sub = m.as_ref().get(0..2, 1..3);
assert_eq!(sub[(0, 0)], 2.0);
assert_eq!(sub[(1, 1)], 6.0);

// mutable element access
m[(2, 2)] = 99.0;
assert_eq!(m[(2, 2)], 99.0);

// transpose (zero-copy)
let t = m.as_ref().transpose();
assert_eq!(t[(2, 0)], 7.0);  // was m[(0,2)]

// split_at: returns (top, bottom) halves
let (top, _bot) = m.as_ref().split_at_row(1);
assert_eq!(top.nrows(), 1);

```

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

### Low-Level Memory Stack Allocation

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Illustrates using `dyn_stack::MemStack` for allocation-free hot paths. Algorithms expose `_scratch` functions to get `StackReq`, which is satisfied by a `MemBuffer` allocation.

```rust
use dyn_stack::{MemBuffer, MemStack, StackReq};
use faer::prelude::*;
use faer::{Accum, Par, linalg, mat};
use faer::mat::AsMatMut;

fn multiply_and_cholesky(
    mut a: faer::MatMut<'_, f64>,
    b: faer::MatRef<'_, f64>,
    stack: &mut MemStack,
) {
    let n = a.nrows();
    // allocate temporary matrix from the stack (no heap)
    let (mut tmp, stack) = unsafe {
        linalg::temp_mat_uninit::<f64, _, _>(n, n, stack)
    };
    let mut tmp = tmp.as_mat_mut();

    linalg::matmul::matmul(&mut tmp, Accum::Replace, &a, &b, 1.0, Par::Seq);
    a.copy_from(&tmp);
    let _ = linalg::cholesky::llt::factor::cholesky_in_place(
        a.rb_mut(), Default::default(), Par::Seq, stack, Default::default(),
    );
}

fn main() {
    let n = 64usize;
    // pre-compute scratch requirements
    let tmp_req = linalg::temp_mat_scratch::<f64>(n, n);
    let llt_req = linalg::cholesky::llt::factor::cholesky_in_place_scratch::<f64>(
        n, Par::Seq, Default::default(),
    );
    let total = tmp_req.and(StackReq::any_of(&[tmp_req.and(llt_req), tmp_req]));

    // single allocation that covers all scratch needs
    let mut buf = MemBuffer::new(total);
    let stack = MemStack::new(&mut buf);

    let mut a = Mat::<f64>::identity(n, n);
    let b = Mat::<f64>::identity(n, n);
    multiply_and_cholesky(a.as_mut(), b.as_ref(), stack);
}
```

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

### Get Global Parallelism

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Retrieves the currently configured global parallelism level.

```APIDOC
## get_global_parallelism

### Description
Retrieves the currently configured global parallelism strategy. This is useful for checking the active parallelism settings.

### Method
`get_global_parallelism`

### Parameters
None

### Returns
The current `Par` strategy being used globally.

### Request Example
```rust
use faer::{
    Par,
    get_global_parallelism
};

let current_par = get_global_parallelism();
println!("Current parallelism: {:?}", current_par);
```

### Response Example
```json
{
  "current_parallelism": "Par::rayon(4)" 
}
```
```

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

### Create a Matrix in Rust

Source: https://github.com/sarah-quinones/faer-rs/wiki/Home

Demonstrates how to create a matrix using the `mat!` macro and access its elements. Ensure the `faer_core::mat` module is imported.

```rust
using faer_core::mat;

let a = mat![
    [1.0, 2.0, 3.0],
    [4.0, 5.0, 6.0],
];

assert_eq!(a.nrows(), 2);
assert_eq!(a.ncols(), 3);

assert_eq!(a[(1, 2)], 6.0);
```

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

### Create and Manipulate Matrices in Rust

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Demonstrates various methods for creating matrices, including literal arrays, function generation, and macros. Also shows basic arithmetic operations and vector creation.

```rust
use faer::{Mat, Scale, col, mat, row};

// literal matrix (4 rows × 3 cols)
let a = mat![
    [1.0, 5.0,  9.0],
    [2.0, 6.0, 10.0],
    [3.0, 7.0, 11.0],
    [4.0, 8.0, 12.0f64],
];
assert_eq!(a[(0, 0)], 1.0);
assert_eq!(a[(3, 2)], 12.0);

// function-generated matrix
let b = Mat::from_fn(4, 3, |i, j| (i + j) as f64);

// arithmetic operators
let sum   = &a + &b;           // element-wise add
let diff  = &a - &b;           // element-wise subtract
let scaled = Scale(2.0) * &a;  // scalar multiply
let product = &a * b.transpose(); // (4×3) * (3×4) → 4×4 Mat

// column and row vectors
let v = col![1.0_f64, 2.0, 3.0];
let r = row![10.0_f64, 20.0, 30.0];
assert_eq!(v[2], 3.0);
assert_eq!(r[1], 20.0);

// zeros / identity
let zero  = Mat::<f64>::zeros(3, 3);
let eye   = Mat::<f64>::identity(4, 4);

```

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

### Sparse Matrix Construction and Conversion

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Demonstrates how to build sparse matrices in CSC or CSR format from triplets and convert between formats or to dense matrices.

```APIDOC
## Sparse Matrices — `SparseColMat`, `SparseRowMat`, Triplet construction

Sparse matrices are stored in CSC (`SparseColMat`) or CSR (`SparseRowMat`) format. The most convenient construction path is via `try_new_from_triplets`, which accepts a slice of `(row, col, value)` triplets. Conversions between formats and to dense are provided.

```rust
use faer::sparse::{SparseColMat, SparseRowMat, Triplet};

// build a 4×4 sparse matrix from COO triplets
let triplets: Vec<Triplet<usize, f64>> = vec![
    Triplet::new(0, 0, 4.0),
    Triplet::new(1, 1, 3.0),
    Triplet::new(2, 2, 2.0),
    Triplet::new(3, 3, 1.0),
    Triplet::new(0, 1, 1.0),
    Triplet::new(1, 0, 1.0),
];

let csc = SparseColMat::<usize, f64>::try_new_from_triplets(4, 4, &triplets)
    .expect("invalid triplets");

println!("nnz = {}", csc.compute_nnz());

// convert to CSR
let csr = csc.to_row_major().expect("conversion failed");
println!("CSR rows = {}", csr.nrows());

// convert to dense
let dense = csc.to_dense();
assert_eq!(dense[(0, 0)], 4.0);
assert_eq!(dense[(0, 1)], 1.0);
```
```

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

### Configure KaTeX Rendering

Source: https://github.com/sarah-quinones/faer-rs/blob/main/faer/katex-header.html

Use this JavaScript snippet to initialize KaTeX rendering on the entire document body. It defines custom delimiters for inline and display math.

```javascript
document.addEventListener("DOMContentLoaded", function() { renderMathInElement(document.body, { delimiters: [
  {left: "$$", right: "$$", display: true},
  {left: "\\(", right: "\\)", display: false},
  {left: "$", right: "$", display: false},
  {left: " H\\\", right: " H\\\", display: true}
] }); });
```

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

### Perform Coefficient-wise Operations with `zip!` and `unzip!` in Rust

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Illustrates using the `zip!` and `unzip!` macros for element-wise operations on matrices. The traversal order is unspecified and may be parallelized.

```rust
use faer::{Mat, mat, unzip, zip};

let a = mat![[1.0_f64, 3.0, 5.0], [2.0, 4.0, 6.0]];
let b = mat![[7.0_f64, 9.0, 11.0], [8.0, 10.0, 12.0]];
let mut out = Mat::<f64>::zeros(2, 3);

// element-wise add
zip!(&mut out, &a, &b).for_each(|unzip!(o, a, b)| *o = a + b);

assert_eq!(out[(0, 0)], 8.0);   // 1 + 7
assert_eq!(out[(1, 2)], 18.0);  // 6 + 12

// element-wise conditional (in-place clamp to [0, 5])
zip!(&mut out).for_each(|unzip!(x)| {
    *x = x.clamp(0.0, 5.0);
});
assert_eq!(out[(0, 0)], 5.0);   // clamped from 8

```

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

### QR and Column-Pivoted QR Decompositions

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Use `qr` for standard least-squares problems and `col_piv_qr` for improved stability with rank-deficient matrices. Both implement `SolveLstsq` for solving linear systems.

```rust
use faer::{Mat, mat};
use faer::prelude::*;

// overdetermined system (3 equations, 2 unknowns) — least squares
let a = mat![
    [1.0_f64, 1.0],
    [1.0,     2.0],
    [1.0,     3.0],
];
let b = mat![[6.0_f64], [5.0], [7.0]];

// QR least-squares solve: minimises ||A x - b||
let qr = a.as_ref().qr();
let x = qr.solve_lstsq(&b);  // returns 2×1 solution
println!("LS solution: x0={:.4}, x1={:.4}", x[(0,0)], x[(1,0)]);

// column-pivoted QR — more stable for near-rank-deficient A
let cpqr = a.as_ref().col_piv_qr();
let x2 = cpqr.solve_lstsq(&b);
let res = (&a * &x2 - &b).norm_l2();
println!("residual ||Ax-b|| = {res:.6}");
```

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

### Add Two Matrices in Rust

Source: https://github.com/sarah-quinones/faer-rs/wiki/Home

Shows how to perform element-wise addition of two matrices using `cwise().zip().for_each()`. Requires importing `faer_core::mat`.

```rust
using faer_core::mat;

let a = mat![
    [1.0, 2.0, 3.0],
    [4.0, 5.0, 6.0],
];
let b = mat![
    [7.0, 8.0, 9.0],
    [10.0, 11.0, 12.0],
];

let mut sum = mat![
    [0.0, 0.0, 0.0],
    [0.0, 0.0, 0.0],
];

sum
    .as_mut()
    .cwise()
    .zip(a.as_ref())
    .zip(b.as_ref())
    .for_each(|sum, a, b| *sum = a + b);
```

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

### Sparse Matrix Solvers (LLT, LU, QR)

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Performs sparse Cholesky (LLT), LU, and QR factorizations and solves linear systems. Useful for fixed sparsity patterns with changing values.

```rust
use faer::prelude::*;
use faer::sparse::{SparseColMat, Triplet};
use faer::Mat;

// symmetric positive-definite sparse matrix
let triplets = vec![
    Triplet::new(0usize, 0usize, 4.0_f64),
    Triplet::new(1, 1, 5.0),
    Triplet::new(2, 2, 6.0),
    Triplet::new(0, 1, 1.0),
    Triplet::new(1, 0, 1.0),
    Triplet::new(1, 2, 2.0),
    Triplet::new(2, 1, 2.0),
];
let a = SparseColMat::<usize, f64>::try_new_from_triplets(3, 3, &triplets).unwrap();
let rhs = Mat::from_fn(3, 1, |i, _| (i + 1) as f64);

// sparse LLT (Cholesky)
let llt = a.rb().sp_llt(faer::Side::Lower).expect("not positive definite");
let x_llt = llt.solve(&rhs);
println!("sparse LLT solution: {:?}", x_llt.col(0).iter().collect::<Vec<_>>());

// sparse LU (general square)
let lu = a.rb().sp_lu().expect("LU failed");
let x_lu = lu.solve(&rhs);

// verify residual
let res = (a.to_dense() * &x_lu - &rhs).norm_l2();
assert!(res < 1e-10);

// sparse QR (least squares)
let qr = a.rb().sp_qr().expect("QR failed");
let x_qr = qr.solve_lstsq(&rhs);

```

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

### Implement and Use a Matrix-Free Scaling Operator

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Shows how to implement the `LinOp` trait for a custom matrix-free operator, specifically a scaled identity operator. Demonstrates its usage with `partial_eigen` for computing dominant eigenvalues and eigenvectors.

```rust
use faer::{Accum, Col, Conj, Mat, Par, linalg};
use faer::mat::{AsMatMut, MatMut, MatRef};
use faer::matrix_free::{LinOp, eigen::partial_eigen, eigen::partial_eigen_scratch};
use faer::dyn_stack::{MemBuffer, MemStack, StackReq};

// simple scaling operator: applies 2*I
struct ScaledIdentity { n: usize }
impl LinOp<f64> for ScaledIdentity {
    fn nrows(&self) -> usize { self.n }
    fn ncols(&self) -> usize { self.n }
    fn apply_scratch(&self, rhs_ncols: usize, _par: Par) -> StackReq {
        StackReq::empty()
    }
    fn apply(&self, mut out: MatMut<'_, f64>, rhs: MatRef<'_, f64>,
             _par: Par, _stack: &mut MemStack) {
        use faer::{zip, unzip};
        zip!(&mut out, &rhs).for_each(|unzip!(o, r)| *o = 2.0 * r);
    }
    fn conj_apply(&self, out: MatMut<'_, f64>, rhs: MatRef<'_, f64>,
                  par: Par, stack: &mut MemStack) {
        self.apply(out, rhs, par, stack);
    }
}

let op = ScaledIdentity { n: 8 };
let n_eig = 2usize;
let par = Par::Seq;

let scratch_req = partial_eigen_scratch(&op as &dyn LinOp<f64>, n_eig, par, Default::default());
let mut buf = MemBuffer::new(scratch_req);
let stack = MemStack::new(&mut buf);

let mut eigvecs = Mat::<f64>::zeros(8, n_eig);
let mut eigvals = vec![faer::complex::Complex::new(0.0_f64, 0.0); n_eig];
let v0 = Col::<f64>::ones(8);

partial_eigen(
    eigvecs.as_mut(), &mut eigvals, &op, v0.as_ref(),
    f64::EPSILON, par, stack, Default::default(),
);

// 2*I has all eigenvalues = 2
for ev in &eigvals {
    assert!((ev.re - 2.0).abs() < 1e-8, "eigenvalue = {}", ev.re);
}
```

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

### Sparse Solvers

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Provides high-level extension methods for sparse matrix factorization and solving linear systems, including LLT, LU, and QR decompositions.

```APIDOC
## Sparse Solvers — `sp_llt`, `sp_lu`, `sp_qr`

High-level extension methods on sparse matrix references directly mirror the dense API. The two-phase factorization (symbolic + numeric) is exposed for cases where the sparsity pattern is fixed but values change.

```rust
use faer::prelude::*;
use faer::sparse::{SparseColMat, Triplet};
use faer::Mat;

// symmetric positive-definite sparse matrix
let triplets = vec![
    Triplet::new(0usize, 0usize, 4.0_f64),
    Triplet::new(1, 1, 5.0),
    Triplet::new(2, 2, 6.0),
    Triplet::new(0, 1, 1.0),
    Triplet::new(1, 0, 1.0),
    Triplet::new(1, 2, 2.0),
    Triplet::new(2, 1, 2.0),
];
let a = SparseColMat::<usize, f64>::try_new_from_triplets(3, 3, &triplets).unwrap();
let rhs = Mat::from_fn(3, 1, |i, _| (i + 1) as f64);

// sparse LLT (Cholesky)
let llt = a.rb().sp_llt(faer::Side::Lower).expect("not positive definite");
let x_llt = llt.solve(&rhs);
println!("sparse LLT solution: {:?}", x_llt.col(0).iter().collect::<Vec<_>>());

// sparse LU (general square)
let lu = a.rb().sp_lu().expect("LU failed");
let x_lu = lu.solve(&rhs);

// verify residual
let res = (a.to_dense() * &x_lu - &rhs).norm_l2();
assert!(res < 1e-10);

// sparse QR (least squares)
let qr = a.rb().sp_qr().expect("QR failed");
let x_qr = qr.solve_lstsq(&rhs);
```
```

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

### Singular Value Decomposition (SVD)

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Use `svd` for the full SVD ($A = USV^H$), `thin_svd` for a reduced $U$ and $V$, and `singular_values` to obtain only the singular values. The singular values are returned in descending order.

```rust
use faer::{Mat, mat};
use faer::prelude::*;

let a = mat![
    [1.0_f64, 0.0, 0.0],
    [0.0,     2.0, 0.0],
    [0.0,     0.0, 3.0],
    [0.0,     0.0, 0.0],
];

// full SVD
let svd = a.as_ref().svd().expect("SVD failed");
let u = svd.U();   // 4×4 unitary
let s = svd.S();   // diagonal, singular values descending
let v = svd.V();   // 3×3 unitary

println!("singular values: {:.2}, {:.2}, {:.2}",
    s.column_vector()[0], s.column_vector()[1], s.column_vector()[2]);
// → 3.00, 2.00, 1.00

// thin SVD (only first min(m,n)=3 columns of U)
let tsvd = a.as_ref().thin_svd().unwrap();
println!("thin U shape: {}×{}", tsvd.U().nrows(), tsvd.U().ncols());

// singular values only
let sv = a.as_ref().singular_values().unwrap();
assert!((sv[0] - 3.0).abs() < 1e-10);

// reconstruct A from factors: U * diag(S) * V^H
use faer::Scale;
let a_approx = tsvd.U()
    * faer::mat::Mat::from_fn(3, 3, |i, j| {
        if i == j { tsvd.S().column_vector()[i] } else { 0.0 }
    })
    * tsvd.V().adjoint();
let err = (&a.get(0..3, ..) - &a_approx).norm_l2();
assert!(err < 1e-10);
```

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

### Random Matrix Generation (Distributions)

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Generates random matrices using distributions like StandardNormal, ComplexDistribution, and UnitaryMat, compatible with `rand::Rng`.

```rust
use faer::stats::prelude::*;
use faer::c64;

let rng = &mut StdRng::seed_from_u64(42);

// real N(0,1) matrix
let dist = CwiseMatDistribution {
    nrows: 4usize,
    ncols: 3usize,
    dist: StandardNormal,
};
let m: faer::Mat<f64> = dist.sample(rng);
println!("random 4×3 matrix, element (0,0) = {:.4}", m[(0, 0)]);

// complex Gaussian matrix
let c_dist = CwiseMatDistribution {
    nrows: 3usize,
    ncols: 3usize,
    dist: ComplexDistribution::new(StandardNormal, StandardNormal),
};
let cm: faer::Mat<c64> = c_dist.sample(rng);
println!("complex (0,0) = {:.4}+{:.4}i", cm[(0,0)].re, cm[(0,0)].im);

// Haar-distributed unitary matrix
let u_dist = UnitaryMat { dim: 4usize, standard_normal: StandardNormal };
let u: faer::Mat<f64> = u_dist.sample(rng);
// verify U^T * U ≈ I
let utu = u.transpose() * &u;
let err = (&utu - faer::Mat::<f64>::identity(4, 4)).norm_l2();
assert!(err < 1e-10, "not unitary: {err}");

```

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

### Cholesky, LDLT, and Bunch-Kaufman Decompositions

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Use `llt` for positive-definite matrices, `ldlt` for positive-semi-definite, and `lblt` for indefinite self-adjoint matrices. These methods implement `Solve` and `DenseSolve` traits for solving linear systems and provide `reconstruct()` and `inverse()` methods.

```rust
use faer::{Mat, Side, mat};
use faer::prelude::*;

// 3×3 symmetric positive-definite matrix
let a = mat![
    [4.0_f64, 2.0, 1.0],
    [2.0,     5.0, 3.0],
    [1.0,     3.0, 6.0],
];
let rhs = mat![[1.0_f64], [2.0], [3.0]];

// LL^T decomposition (lower triangle is read)
let llt = a.as_ref().llt(Side::Lower).expect("matrix not positive definite");
let x = llt.solve(&rhs);   // solves A x = rhs

// verify: A * x ≈ rhs
let residual = &a * &x - &rhs;
let norm: f64 = residual.as_ref().norm_l2();
assert!(norm < 1e-10, "residual too large: {norm}");

// extract L factor
let l = llt.L();
println!("L factor shape: {}×{}", l.nrows(), l.ncols());

// LDL^T for indefinite case
let ldlt = a.as_ref().ldlt(Side::Lower).unwrap();
let x2 = ldlt.solve(&rhs);
let res2 = (&a * &x2 - &rhs).norm_l2();
assert!(res2 < 1e-10);
```

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

### Concatenate Matrices using `concat!` Macro in Rust

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Demonstrates how to assemble a block matrix from submatrices using the `concat!` macro, following NumPy's `block` convention. Sub-lists are concatenated horizontally, then stacked vertically.

```rust
use faer::{Mat, concat, mat};

let a = mat![[1.0_f64, 2.0], [3.0, 4.0]];        // 2×2
let b = mat![[5.0_f64], [6.0]];                   // 2×1
let c = mat![[7.0_f64, 8.0, 9.0]];               // 1×3

// block layout:
//   [ a | b ]   →  2×3
//   [   c   ]   →  1×3
let block = concat![[&a, &b], [&c]];

assert_eq!(block.nrows(), 3);
assert_eq!(block.ncols(), 3);
assert_eq!(block[(0, 0)], 1.0);
assert_eq!(block[(0, 2)], 5.0);
assert_eq!(block[(2, 1)], 8.0);

```

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

### Build Sparse Matrix from Triplets

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Constructs a sparse matrix in CSC or CSR format using triplets (row, column, value). Supports conversion to other formats and dense matrices.

```rust
use faer::sparse::{SparseColMat, SparseRowMat, Triplet};

// build a 4×4 sparse matrix from COO triplets
let triplets: Vec<Triplet<usize, f64>> = vec![
    Triplet::new(0, 0, 4.0),
    Triplet::new(1, 1, 3.0),
    Triplet::new(2, 2, 2.0),
    Triplet::new(3, 3, 1.0),
    Triplet::new(0, 1, 1.0),
    Triplet::new(1, 0, 1.0),
];

let csc = SparseColMat::<usize, f64>::try_new_from_triplets(4, 4, &triplets)
    .expect("invalid triplets");

println!("nnz = {}", csc.compute_nnz());

// convert to CSR
let csr = csc.to_row_major().expect("conversion failed");
println!("CSR rows = {}", csr.nrows());

// convert to dense
let dense = csc.to_dense();
assert_eq!(dense[(0, 0)], 4.0);
assert_eq!(dense[(0, 1)], 1.0);

```

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

### LU Decomposition (Partial and Full Pivoting)

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Use `partial_piv_lu` for general square invertible matrices and `full_piv_lu` for rectangular or rank-deficient matrices. Both implement `Solve` and `DenseSolve` traits. `full_piv_lu` provides a rank estimate.

```rust
use faer::{Mat, mat};
use faer::prelude::*;

let a = mat![
    [2.0_f64, 1.0, -1.0],
    [-3.0,   -1.0,  2.0],
    [-2.0,    1.0,  2.0],
];
let b = mat![[8.0_f64], [-11.0], [-3.0]];

// partial-pivot LU — standard general solver
let lu = a.as_ref().partial_piv_lu();
let x = lu.solve(&b);
let res = (&a * &x - &b).norm_l2();
assert!(res < 1e-10);

// determinant via LU
let det = lu.determinant();
println!("det(A) = {det}");

// full-pivot LU for rectangular / rank-deficient matrices
let a_rect = mat![
    [1.0_f64, 2.0, 3.0],
    [4.0,     5.0, 6.0],
];
let fpu = a_rect.as_ref().full_piv_lu();
println!("rank estimate: {}", fpu.rank());

// inverse
let a_inv = lu.inverse();
let eye_approx = &a * &a_inv;
assert!((eye_approx[(0, 0)] - 1.0).abs() < 1e-10);
```

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

### Random Matrix Generation

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Utilities for generating random matrices using distributions that implement `rand::distributions::Distribution`, including Gaussian and Haar-distributed unitary matrices.

```APIDOC
## Random Matrix Generation — `CwiseMatDistribution`, `UnitaryMat`

With the `rand` feature (enabled by default), `faer::stats` exposes distribution wrappers that implement `rand::distributions::Distribution<Mat<T>>`, making random matrix generation composable with any Rng.

```rust
use faer::stats::prelude::*;
use faer::c64;

let rng = &mut StdRng::seed_from_u64(42);

// real N(0,1) matrix
let dist = CwiseMatDistribution {
    nrows: 4usize,
    ncols: 3usize,
    dist: StandardNormal,
};
let m: faer::Mat<f64> = dist.sample(rng);
println!("random 4×3 matrix, element (0,0) = {:.4}", m[(0, 0)]);

// complex Gaussian matrix
let c_dist = CwiseMatDistribution {
    nrows: 3usize,
    ncols: 3usize,
    dist: ComplexDistribution::new(StandardNormal, StandardNormal),
};
let cm: faer::Mat<c64> = c_dist.sample(rng);
println!("complex (0,0) = {:.4}+{:.4}i", cm[(0,0)].re, cm[(0,0)].im);

// Haar-distributed unitary matrix
let u_dist = UnitaryMat { dim: 4usize, standard_normal: StandardNormal };
let u: faer::Mat<f64> = u_dist.sample(rng);
// verify U^T * U ≈ I
let utu = u.transpose() * &u;
let err = (&utu - faer::Mat::<f64>::identity(4, 4)).norm_l2();
assert!(err < 1e-10, "not unitary: {err}");
```
```

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

### Load and Save NumPy .npy Files with faer

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Demonstrates loading a matrix from a .npy byte buffer and saving a matrix to .npy bytes. Handles zero-copy loading for aligned buffers and copies for unaligned buffers. Supports f32, f64, c32, and c64 dtypes.

```rust
#[cfg(feature = "npy")]
{
    use faer::io::npy::{FromNpy, Npy};
    use faer::c64;

    // assume `bytes` is a 64-byte aligned buffer containing an npy file
    // (e.g. memory-mapped file or aligned Vec<u8>)
    fn load_matrix(bytes: &[u8]) -> faer::Mat<f64> {
        let view = Npy::new(bytes).expect("invalid npy format");
        assert_eq!(view.dtype(), faer::io::npy::NpyDType::F64);

        if view.is_aligned() {
            // zero-copy view
            let mat_ref = view.as_aligned_ref::<f64>();
            mat_ref.to_owned()
        } else {
            // copy into owned matrix
            view.read::<f64>().expect("type mismatch")
        }
    }

    // write a matrix to npy bytes
    fn save_matrix(m: faer::MatRef<'_, f64>) -> Vec<u8> {
        let mut buf = Vec::new();
        faer::io::npy::write_npy(&mut buf, m).expect("write failed");
        buf
    }
}
```

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

### Low-Level Memory Stack Operations

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Provides abstractions for memory management using a bump-allocator (`MemStack`) for allocation-free hot paths, including scratch space calculation and allocation.

```APIDOC
## dyn_stack::MemStack, StackReq, temp_mat_scratch

### Description
These components facilitate efficient memory management for performance-critical code paths by using a bump-allocator (`MemStack`). Algorithms can expose their scratch space requirements via `StackReq`, which can then be satisfied by a single `MemBuffer` allocation. Multiple requirements can be combined using `all_of` or `any_of`.

### Methods
- `dyn_stack::MemStack::new`: Creates a new memory stack from a `MemBuffer`.
- `faer::linalg::temp_mat_scratch`: Computes the scratch space requirement for a temporary matrix.
- `faer::linalg::cholesky::llt::factor::cholesky_in_place_scratch`: Computes the scratch space requirement for in-place Cholesky factorization.
- `StackReq::and`: Combines two stack requirements sequentially.
- `StackReq::any_of`: Combines multiple stack requirements, allowing for reuse.

### Parameters
- `stack` (&mut MemStack): The memory stack to allocate from.
- `size` (usize): The dimensions for matrix scratch space.
- `par` (Par): Parallelism strategy.

### Request Example
```rust
use dyn_stack::{MemBuffer, MemStack, StackReq};
use faer::prelude::*;
use faer::{Accum, Par, linalg, mat};
use faer::mat::AsMatMut;

fn multiply_and_cholesky(
    mut a: faer::MatMut<'_, f64>,
    b: faer::MatRef<'_, f64>,
    stack: &mut MemStack,
) {
    let n = a.nrows();
    // allocate temporary matrix from the stack (no heap)
    let (mut tmp, stack) = unsafe {
        linalg::temp_mat_uninit::<f64, _, _>(n, n, stack)
    };
    let mut tmp = tmp.as_mat_mut();

    linalg::matmul::matmul(&mut tmp, Accum::Replace, &a, &b, 1.0, Par::Seq);
    a.copy_from(&tmp);
    let _ = linalg::cholesky::llt::factor::cholesky_in_place(
        a.rb_mut(), Default::default(), Par::Seq, stack, Default::default(),
    );
}

fn main() {
    let n = 64usize;
    // pre-compute scratch requirements
    let tmp_req = linalg::temp_mat_scratch::<f64>(n, n);
    let llt_req = linalg::cholesky::llt::factor::cholesky_in_place_scratch::<f64>(
        n, Par::Seq, Default::default(),
    );
    let total = tmp_req.and(StackReq::any_of(&[tmp_req.and(llt_req), tmp_req]));

    // single allocation that covers all scratch needs
    let mut buf = MemBuffer::new(total);
    let stack = MemStack::new(&mut buf);

    let mut a = Mat::<f64>::identity(n, n);
    let b = Mat::<f64>::identity(n, n);
    multiply_and_cholesky(a.as_mut(), b.as_ref(), stack);
}
```

### Response Example
```json
{
  "allocation_status": "success",
  "memory_used": "1024 bytes"
}
```
```

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

### LU Decompositions (Partial & Full Pivoting)

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Performs LU decomposition on square invertible matrices using partial pivoting ($PA=LU$) or rectangular/rank-deficient matrices using full pivoting ($PAQ^T=LU$). These decompositions are useful for solving linear systems, calculating determinants, and finding matrix inverses.

```APIDOC
## Dense Decompositions — LU (Partial & Full Pivoting)

`Mat::partial_piv_lu` decomposes a square invertible matrix as $PA = LU$ (recommended default for general square systems and determinants). `Mat::full_piv_lu` decomposes a rectangular matrix as $PAQ^T = LU$ and is more stable for rank-deficient matrices. Both implement `Solve` and `DenseSolve`.

```rust
use faer::{Mat, mat};
use faer::prelude::*;

let a = mat![
    [2.0_f64, 1.0, -1.0],
    [-3.0,   -1.0,  2.0],
    [-2.0,    1.0,  2.0],
];
let b = mat![[8.0_f64], [-11.0], [-3.0]];

// partial-pivot LU — standard general solver
let lu = a.as_ref().partial_piv_lu();
let x = lu.solve(&b);
let res = (&a * &x - &b).norm_l2();
assert!(res < 1e-10);

// determinant via LU
let det = lu.determinant();
println!("det(A) = {det}");

// full-pivot LU for rectangular / rank-deficient matrices
let a_rect = mat![
    [1.0_f64, 2.0, 3.0],
    [4.0,     5.0, 6.0],
];
let fpu = a_rect.as_ref().full_piv_lu();
println!("rank estimate: {}", fpu.rank());

// inverse
let a_inv = lu.inverse();
let eye_approx = &a * &a_inv;
assert!((eye_approx[(0, 0)] - 1.0).abs() < 1e-10);
```
```

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

### Multiply Two Matrices in Rust

Source: https://github.com/sarah-quinones/faer-rs/wiki/Home

Illustrates matrix multiplication using the `matmul` function. Requires importing `faer_core::{mat, mul::matmul, Conj, Parallelism}`.

```rust
using faer_core::{mat, mul::matmul, Conj, Parallelism};

let a = mat![
    [1.0, 2.0, 3.0],
    [4.0, 5.0, 6.0],
];
let b = mat![
    [7.0, 8.0],
    [9.0, 10.0],
    [11.0, 12.0],
];

let mut mul = mat![
    [0.0, 0.0],
    [0.0, 0.0],
];

matmul(
    mul.as_mut(),
    Conj::No,
    a.as_ref(),
    Conj::No,
    b.as_ref(),
    Conj::No,
    None,
    1.0,
    Parallelism::None
);
```

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

### Self-Adjoint Eigendecomposition (Eigenvalues Only)

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Efficiently computes only the eigenvalues for a self-adjoint matrix, useful when eigenvectors are not needed. Requires the `faer::prelude::*` import.

```rust
use faer::{Mat, Side, mat};
use faer::prelude::*;

// eigenvalues only (cheaper)
let a = mat![
    [2.0_f64, 1.0],
    [1.0,     3.0],
];
let ev = a.as_ref().self_adjoint_eigenvalues(Side::Lower).unwrap();
let eigenvalues = ev;
assert!((ev[0] - 1.3820).abs() < 1e-4);
```

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

### General Eigendecomposition

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Computes the eigendecomposition for general square matrices, which may produce complex eigenvalues. Requires the `faer::prelude::*` import.

```rust
use faer::{Mat, Side, mat};
use faer::prelude::*;

// general (non-symmetric) eigendecomposition — produces complex results
let b = mat![[0.0_f64, -1.0], [1.0, 0.0]];  // rotation matrix
let gen = b.as_ref().eigen().unwrap();
let s_vals = gen.S().column_vector();
println!("eigenvalues of rotation: {:.3}+{:.3}i, {:.3}+{:.3}i",
    s_vals[0].re, s_vals[0].im, s_vals[1].re, s_vals[1].im);
// → eigenvalues ≈ +i and -i
```

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

### QR and Column-Pivoted QR Decompositions

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Computes the QR decomposition ($A=QR$) for solving linear least-squares problems. Column-pivoted QR ($AP^T=QR$) is used for improved stability with rank-deficient matrices.

```APIDOC
## Dense Decompositions — QR and Column-Pivoted QR

`Mat::qr` computes $A = QR$ where $Q$ is unitary and $R$ is upper trapezoidal — the standard least-squares solver. `Mat::col_piv_qr` computes $AP^T = QR$ (column pivoting) for improved stability on rank-deficient matrices and implements `SolveLstsq`.

```rust
use faer::{Mat, mat};
use faer::prelude::*;

// overdetermined system (3 equations, 2 unknowns) — least squares
let a = mat![
    [1.0_f64, 1.0],
    [1.0,     2.0],
    [1.0,     3.0],
];
let b = mat![[6.0_f64], [5.0], [7.0]];

// QR least-squares solve: minimises ||A x - b||
let qr = a.as_ref().qr();
let x = qr.solve_lstsq(&b);  // returns 2×1 solution
println!("LS solution: x0={:.4}, x1={:.4}", x[(0,0)], x[(1,0)]);

// column-pivoted QR — more stable for near-rank-deficient A
let cpqr = a.as_ref().col_piv_qr();
let x2 = cpqr.solve_lstsq(&b);
let res = (&a * &x2 - &b).norm_l2();
println!("residual ||Ax-b|| = {res:.6}");
```
```

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

### Column and Row Statistics (Mean, Variance)

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Calculates column-wise and row-wise means and variances. Supports configurable NaN handling (Propagate or Ignore).

```rust
use faer::mat;
use faer::stats::{NanHandling, col_mean, col_varm, row_mean};
use faer::Col;

let data = mat![
    [1.0_f64, 2.0, 3.0],
    [4.0,     5.0, 6.0],
    [7.0,     8.0, 9.0],
];

// column means → [4.0, 5.0, 6.0]
let mut means = Col::<f64>::zeros(3);
col_mean(means.as_mut(), data.as_ref(), NanHandling::Propagate);
assert!((means[0] - 4.0).abs() < 1e-10);
assert!((means[1] - 5.0).abs() < 1e-10);

// row means → [2.0, 5.0, 8.0]
let mut row_means = Col::<f64>::zeros(3);
row_mean(row_means.as_mut(), data.as_ref(), NanHandling::Propagate);
assert!((row_means[0] - 2.0).abs() < 1e-10);

// column variance (ddof=1, i.e. unbiased)
let mut var = Col::<f64>::zeros(3);
col_varm(var.as_mut(), data.as_ref(), means.as_ref(), NanHandling::Propagate);
// each column has values spaced 3 apart, variance = (9+0+9)/2 = 9
assert!((var[0] - 9.0).abs() < 1e-10);

```

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

### Explicit Parallelism for Matrix Multiplication

Source: https://context7.com/sarah-quinones/faer-rs/llms.txt

Demonstrates passing `Par` explicitly to a low-level routine like `matmul` for fine-grained control over parallelism. Requires `faer::linalg::matmul::matmul` and other Faer types.

```rust
use faer::{Par, Mat, Accum};
use faer::linalg::matmul::matmul;

let a = Mat::<f64>::identity(4, 4);
let b = Mat::<f64>::identity(4, 4);
let mut c = Mat::<f64>::zeros(4, 4);
matmul(&mut c, Accum::Replace, &a, &b, 1.0, Par::Seq);
assert_eq!(c[(0, 0)], 1.0);
```