scivex-core 0.1.1

Scivex — Tensors, linear algebra, FFT, and math primitives
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84

//! Linear algebra operations.
//!
//! All routines are implemented from scratch — no external BLAS/LAPACK
//! bindings. The API is split into three BLAS levels plus matrix
//! decompositions:
//!
//! | Level | Operations | Complexity |
//! |-------|-----------|------------|
//! | L1 | `dot`, `axpy`, `nrm2`, `asum`, `scal`, `iamax` | O(n) |
//! | L2 | `gemv` (matrix-vector multiply) | O(n^2) |
//! | L3 | `gemm` (matrix-matrix multiply) | O(n^3) |
//!
//! Decompositions: `LuDecomposition`, `QrDecomposition`,
//! `CholeskyDecomposition`, `SvdDecomposition`, `EigDecomposition`

pub mod blas;
pub mod decomp;
pub mod sparse;
pub mod tensor_decomp;

pub use blas::{asum, axpy, dot, gemm, gemv, iamax, nrm2, scal};
pub use decomp::CholeskyDecomposition;
pub use decomp::EigDecomposition;
pub use decomp::LuDecomposition;
pub use decomp::QrDecomposition;
pub use decomp::SvdDecomposition;
pub use decomp::lstsq;
pub use sparse::{CooMatrix, CscMatrix, CsrMatrix};
pub use tensor_decomp::{CpDecomposition, NtfDecomposition, TuckerDecomposition};

use crate::Float;
use crate::error::Result;
use crate::tensor::Tensor;

/// Solve the linear system `Ax = b` for a square matrix `A`.
///
/// Uses LU decomposition with partial pivoting internally.
///
/// ```
/// # use scivex_core::tensor::Tensor;
/// # use scivex_core::linalg;
/// let a = Tensor::from_vec(vec![2.0_f64, 1.0, 1.0, 4.0], vec![2, 2]).unwrap();
/// let b = Tensor::from_vec(vec![5.0_f64, 6.0], vec![2]).unwrap();
/// let x = linalg::solve(&a, &b).unwrap();
/// assert!((x.as_slice()[0] - 2.0).abs() < 1e-10);
/// assert!((x.as_slice()[1] - 1.0).abs() < 1e-10);
/// ```
pub fn solve<T: Float>(a: &Tensor<T>, b: &Tensor<T>) -> Result<Tensor<T>> {
    LuDecomposition::decompose(a)?.solve(b)
}

/// Compute the inverse of a square matrix.
///
/// Uses LU decomposition with partial pivoting internally.
/// Returns [`CoreError::SingularMatrix`](crate::CoreError::SingularMatrix) if the
/// matrix is singular.
///
/// ```
/// # use scivex_core::tensor::Tensor;
/// # use scivex_core::linalg;
/// let a = Tensor::from_vec(vec![2.0_f64, 1.0, 1.0, 4.0], vec![2, 2]).unwrap();
/// let inv = linalg::inv(&a).unwrap();
/// // A * A^-1 ≈ I
/// let eye = a.matmul(&inv).unwrap();
/// assert!((eye.as_slice()[0] - 1.0).abs() < 1e-10);
/// ```
pub fn inv<T: Float>(a: &Tensor<T>) -> Result<Tensor<T>> {
    LuDecomposition::decompose(a)?.inverse()
}

/// Compute the determinant of a square matrix.
///
/// Uses LU decomposition with partial pivoting internally.
///
/// ```
/// # use scivex_core::tensor::Tensor;
/// # use scivex_core::linalg;
/// let a = Tensor::from_vec(vec![2.0_f64, 1.0, 1.0, 4.0], vec![2, 2]).unwrap();
/// let det = linalg::det(&a).unwrap();
/// assert!((det - 7.0).abs() < 1e-10);
/// ```
pub fn det<T: Float>(a: &Tensor<T>) -> Result<T> {
    Ok(LuDecomposition::decompose(a)?.det())
}