russell
russell copied to clipboard
cpmech
Rust Scientific Libary. ODE and DAE (Runge-Kutta) solvers. Special functions (Bessel, Elliptic, Beta, Gamma, Erf). Linear algebra. Sparse solvers (MUMPS, UMFPACK). Probability distributions. Tensor ca...
Russell - Rust Scientific Library
(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:
chk Functions to check vectors and other data in tests
lab Matrix-vector laboratory including linear algebra tools
openblas Thin wrapper to some OpenBLAS routines
sparse Sparse matrix tools and solvers
stat Statistics calculations, probability distributions, and pseudo random numbers
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
cpmech