russell icon indicating copy to clipboard operation
russell copied to clipboard

Rust Scientific Libary. Matrix-vector laboratory, OpenBLAS, sparse direct solvers.

Russell - Rust Scientific Library

codecov

Bertrand Russell

(CC0. Photo: Bertrand Russell)

Russell assists in the development of scientific computations using the Rust language. We focus on numerical methods and solvers for differential equations; however, anything is possible ๐Ÿ˜‰.

An essential goal of this library is to bring the best (fastest) solutions while maintaining a very clean (and idiomatic) code, thoroughly tested (min coverage of 95%), and yet simple to use. The best solutions are brought by wrapping powerful libraries such as OpenBLAS, MUMPS, and SuiteSparse (UMFPACK).

Available crates:

  • Crates.io chk Functions to check vectors and other data in tests
  • Crates.io lab Matrix-vector laboratory including linear algebra tools
  • Crates.io openblas Thin wrapper to some OpenBLAS routines
  • Crates.io sparse Sparse matrix tools and solvers
  • Crates.io stat Statistics calculations, probability distributions, and pseudo random numbers
  • Crates.io tensor Tensor analysis structures and functions for continuum mechanics

External recommended crate:

  • plotpy Plotting tools using Python3/Matplotlib as an engine

Installation

Install the following Debian packages:

sudo apt-get install \
    liblapacke-dev \
    libmumps-seq-dev \
    libopenblas-dev \
    libsuitesparse-dev

Add this to your Cargo.toml (select only the crates you want and replace the right version):

[dependencies]
russell_chk = "*"
russell_lab = "*"
russell_openblas = "*"
russell_sparse = "*"
russell_stat = "*"
russell_tensor = "*"

Number of threads

By default OpenBLAS will use all available threads, including Hyper-Threads that make the performance worse. Thus, it is best to set the following environment variable:

export OPENBLAS_NUM_THREADS=<real-core-count>

Furthermore, if working on a multi-threaded application, it is recommended to set:

export OPENBLAS_NUM_THREADS=1

Examples

Compute a singular value decomposition

use russell_lab::{sv_decomp, Matrix, Vector, StrError};

fn main() -> Result<(), StrError> {
    // set matrix
    let mut a = Matrix::from(&[
        [2.0, 4.0],
        [1.0, 3.0],
        [0.0, 0.0],
        [0.0, 0.0],
    ]);

    // allocate output structures
    let (m, n) = a.dims();
    let min_mn = if m < n { m } else { n };
    let mut s = Vector::new(min_mn);
    let mut u = Matrix::new(m, m);
    let mut vt = Matrix::new(n, n);

    // perform SVD
    sv_decomp(&mut s, &mut u, &mut vt, &mut a)?;

    // define correct data
    let s_correct = "โ”Œ      โ”\n\
                     โ”‚ 5.46 โ”‚\n\
                     โ”‚ 0.37 โ”‚\n\
                     โ””      โ”˜";
    let u_correct = "โ”Œ                         โ”\n\
                     โ”‚ -0.82 -0.58  0.00  0.00 โ”‚\n\
                     โ”‚ -0.58  0.82  0.00  0.00 โ”‚\n\
                     โ”‚  0.00  0.00  1.00  0.00 โ”‚\n\
                     โ”‚  0.00  0.00  0.00  1.00 โ”‚\n\
                     โ””                         โ”˜";
    let vt_correct = "โ”Œ             โ”\n\
                      โ”‚ -0.40 -0.91 โ”‚\n\
                      โ”‚ -0.91  0.40 โ”‚\n\
                      โ””             โ”˜";

    // check solution
    assert_eq!(format!("{:.2}", s), s_correct);
    assert_eq!(format!("{:.2}", u), u_correct);
    assert_eq!(format!("{:.2}", vt), vt_correct);

    // check SVD: a == u * s * vt
    let mut usv = Matrix::new(m, n);
    for i in 0..m {
        for j in 0..n {
            for k in 0..min_mn {
                usv[i][j] += u[i][k] * s[k] * vt[k][j];
            }
        }
    }
    let usv_correct = "โ”Œ     โ”\n\
                       โ”‚ 2 4 โ”‚\n\
                       โ”‚ 1 3 โ”‚\n\
                       โ”‚ 0 0 โ”‚\n\
                       โ”‚ 0 0 โ”‚\n\
                       โ””     โ”˜";
    assert_eq!(format!("{}", usv), usv_correct);
    Ok(())
}

Solve a linear system

use russell_lab::{solve_lin_sys, Matrix, Vector, StrError};

fn main() -> Result<(), StrError> {
    // set matrix and right-hand side
    let mut a = Matrix::from(&[
        [1.0,  3.0, -2.0],
        [3.0,  5.0,  6.0],
        [2.0,  4.0,  3.0],
    ]);
    let mut b = Vector::from(&[5.0, 7.0, 8.0]);

    // solve linear system b := aโปยนโ‹…b
    solve_lin_sys(&mut b, &mut a)?;

    // check
    let x_correct = "โ”Œ         โ”\n\
                     โ”‚ -15.000 โ”‚\n\
                     โ”‚   8.000 โ”‚\n\
                     โ”‚   2.000 โ”‚\n\
                     โ””         โ”˜";
    assert_eq!(format!("{:.3}", b), x_correct);
    Ok(())
}

Solve a sparse linear system

use russell_lab::{Matrix, Vector, StrError};
use russell_sparse::{ConfigSolver, Solver, SparseTriplet, Symmetry};

fn main() -> Result<(), StrError> {

    // allocate a square matrix
    let mut trip = SparseTriplet::new(5, 5, 13, Symmetry::No)?;
    trip.put(0, 0,  1.0)?; // << (0, 0, a00/2)
    trip.put(0, 0,  1.0)?; // << (0, 0, a00/2)
    trip.put(1, 0,  3.0)?;
    trip.put(0, 1,  3.0)?;
    trip.put(2, 1, -1.0)?;
    trip.put(4, 1,  4.0)?;
    trip.put(1, 2,  4.0)?;
    trip.put(2, 2, -3.0)?;
    trip.put(3, 2,  1.0)?;
    trip.put(4, 2,  2.0)?;
    trip.put(2, 3,  2.0)?;
    trip.put(1, 4,  6.0)?;
    trip.put(4, 4,  1.0)?;

    // print matrix
    let (m, n) = trip.dims();
    let mut a = Matrix::new(m, n);
    trip.to_matrix(&mut a)?;
    let correct = "โ”Œ                โ”\n\
                   โ”‚  2  3  0  0  0 โ”‚\n\
                   โ”‚  3  0  4  0  6 โ”‚\n\
                   โ”‚  0 -1 -3  2  0 โ”‚\n\
                   โ”‚  0  0  1  0  0 โ”‚\n\
                   โ”‚  0  4  2  0  1 โ”‚\n\
                   โ””                โ”˜";
    assert_eq!(format!("{}", a), correct);

    // allocate x and rhs
    let mut x = Vector::new(5);
    let rhs = Vector::from(&[8.0, 45.0, -3.0, 3.0, 19.0]);

    // initialize, factorize, and solve
    let config = ConfigSolver::new();
    let mut solver = Solver::new(config)?;
    solver.initialize(&trip)?;
    solver.factorize()?;
    solver.solve(&mut x, &rhs)?;
    let correct = "โ”Œ          โ”\n\
                   โ”‚ 1.000000 โ”‚\n\
                   โ”‚ 2.000000 โ”‚\n\
                   โ”‚ 3.000000 โ”‚\n\
                   โ”‚ 4.000000 โ”‚\n\
                   โ”‚ 5.000000 โ”‚\n\
                   โ””          โ”˜";
    assert_eq!(format!("{:.6}", x), correct);
    Ok(())
}

Todo list

  • [x] Add complex numbers functions to russell_openblas
  • [ ] Add more complex numbers functions to russell_lab
  • [ ] Add fundamental functions to russell_lab
    • [ ] Implement the modified Bessel functions
  • [ ] Implement some numerical methods in russell_lab
    • [ ] Implement Brent's solver
    • [ ] Implement solver for the cubic equation
    • [ ] Implement numerical derivation
    • [ ] Implement numerical Jacobian function
    • [ ] Implement Newton's method for nonlinear systems
    • [ ] Implement numerical quadrature
  • [ ] Add interpolation and polynomials to russell_lab
    • [ ] Implement Chebyshev interpolation and polynomials
    • [ ] Implement Orthogonal polynomials
    • [ ] Implement Lagrange interpolation
  • [x] Add probability distribution functions to russell_stat
  • [x] Finalize drawing of ASCII histogram in russell_stat
  • [ ] Implement standard continuum mechanics tensors in russell_tensor
  • [ ] Implement more integration tests for linear algebra
  • [ ] Implement more examples