Fix PRNG seed for random vectors and matrices in tests

ndarray-linalg/Cargo.toml

@@ -51,6 +51,7 @@ paste = "1.0.5"
 criterion = "0.3.4"
 # Keep the same version as ndarray's dependency!
 approx = { version = "0.4.0", features = ["num-complex"] }
+rand_pcg = "0.3.1"
 
 [[bench]]
 name = "truncated_eig"

ndarray-linalg/src/lobpcg/lobpcg.rs

@@ -460,7 +460,8 @@ mod tests {
     /// Test the `sorted_eigen` function
     #[test]
     fn test_sorted_eigen() {
-        let matrix: Array2<f64> = generate::random((10, 10)) * 10.0;
+        let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+        let matrix: Array2<f64> = generate::random_using((10, 10), &mut rng) * 10.0;
         let matrix = matrix.t().dot(&matrix);
 
         // return all eigenvectors with largest first
@@ -476,7 +477,8 @@ mod tests {
     /// Test the masking function
     #[test]
     fn test_masking() {
-        let matrix: Array2<f64> = generate::random((10, 5)) * 10.0;
+        let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+        let matrix: Array2<f64> = generate::random_using((10, 5), &mut rng) * 10.0;
         let masked_matrix = ndarray_mask(matrix.view(), &[true, true, false, true, false]);
         close_l2(
             &masked_matrix.slice(s![.., 2]),
@@ -488,7 +490,8 @@ mod tests {
     /// Test orthonormalization of a random matrix
     #[test]
     fn test_orthonormalize() {
-        let matrix: Array2<f64> = generate::random((10, 10)) * 10.0;
+        let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+        let matrix: Array2<f64> = generate::random_using((10, 10), &mut rng) * 10.0;
 
         let (n, l) = orthonormalize(matrix.clone()).unwrap();
 
@@ -509,7 +512,8 @@ mod tests {
         assert_symmetric(a);
 
         let n = a.len_of(Axis(0));
-        let x: Array2<f64> = generate::random((n, num));
+        let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+        let x: Array2<f64> = generate::random_using((n, num), &mut rng);
 
         let result = lobpcg(|y| a.dot(&y), x, |_| {}, None, 1e-5, n * 2, order);
         match result {
@@ -553,7 +557,8 @@ mod tests {
     #[test]
     fn test_eigsolver_constructed() {
         let n = 50;
-        let tmp = generate::random((n, n));
+        let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+        let tmp = generate::random_using((n, n), &mut rng);
         //let (v, _) = tmp.qr_square().unwrap();
         let (v, _) = orthonormalize(tmp).unwrap();
 
@@ -570,7 +575,8 @@ mod tests {
     fn test_eigsolver_constrained() {
         let diag = arr1(&[1., 2., 3., 4., 5., 6., 7., 8., 9., 10.]);
         let a = Array2::from_diag(&diag);
-        let x: Array2<f64> = generate::random((10, 1));
+        let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+        let x: Array2<f64> = generate::random_using((10, 1), &mut rng);
         let y: Array2<f64> = arr2(&[
             [1.0, 0., 0., 0., 0., 0., 0., 0., 0., 0.],
             [0., 1.0, 0., 0., 0., 0., 0., 0., 0., 0.],

ndarray-linalg/src/lobpcg/svd.rs

@@ -214,7 +214,8 @@ mod tests {
 
     #[test]
     fn test_truncated_svd_random() {
-        let a: Array2<f64> = generate::random((50, 10));
+        let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+        let a: Array2<f64> = generate::random_using((50, 10), &mut rng);
 
         let res = TruncatedSvd::new(a.clone(), Order::Largest)
             .precision(1e-5)

ndarray-linalg/src/solve.rs

@@ -5,20 +5,13 @@
 //! Solve `A * x = b`:
 //!
 //! ```
-//! #[macro_use]
-//! extern crate ndarray;
-//! extern crate ndarray_linalg;
-//!
 //! use ndarray::prelude::*;
 //! use ndarray_linalg::Solve;
-//! # fn main() {
 //!
 //! let a: Array2<f64> = array![[3., 2., -1.], [2., -2., 4.], [-2., 1., -2.]];
 //! let b: Array1<f64> = array![1., -2., 0.];
 //! let x = a.solve_into(b).unwrap();
 //! assert!(x.abs_diff_eq(&array![1., -2., -2.], 1e-9));
-//!
-//! # }
 //! ```
 //!
 //! There are also special functions for solving `A^T * x = b` and
@@ -29,21 +22,18 @@
 //! the beginning than solving directly using `A`:
 //!
 //! ```
-//! # extern crate ndarray;
-//! # extern crate ndarray_linalg;
-//!
 //! use ndarray::prelude::*;
 //! use ndarray_linalg::*;
-//! # fn main() {
 //!
-//! let a: Array2<f64> = random((3, 3));
+//! /// Use fixed algorithm and seed of PRNG for reproducible test
+//! let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+//!
+//! let a: Array2<f64> = random_using((3, 3), &mut rng);
 //! let f = a.factorize_into().unwrap(); // LU factorize A (A is consumed)
 //! for _ in 0..10 {
-//!     let b: Array1<f64> = random(3);
+//!     let b: Array1<f64> = random_using(3, &mut  rng);
 //!     let x = f.solve_into(b).unwrap(); // Solve A * x = b using factorized L, U
 //! }
-//!
-//! # }
 //! ```
 
 use ndarray::*;

ndarray-linalg/src/solveh.rs

@@ -8,13 +8,8 @@
 //! Solve `A * x = b`, where `A` is a Hermitian (or real symmetric) matrix:
 //!
 //! ```
-//! #[macro_use]
-//! extern crate ndarray;
-//! extern crate ndarray_linalg;
-//!
 //! use ndarray::prelude::*;
 //! use ndarray_linalg::SolveH;
-//! # fn main() {
 //!
 //! let a: Array2<f64> = array![
 //!     [3., 2., -1.],
@@ -24,29 +19,25 @@
 //! let b: Array1<f64> = array![11., -12., 1.];
 //! let x = a.solveh_into(b).unwrap();
 //! assert!(x.abs_diff_eq(&array![1., 3., -2.], 1e-9));
-//!
-//! # }
 //! ```
 //!
 //! If you are solving multiple systems of linear equations with the same
 //! Hermitian or real symmetric coefficient matrix `A`, it's faster to compute
 //! the factorization once at the beginning than solving directly using `A`:
 //!
 //! ```
-//! # extern crate ndarray;
-//! # extern crate ndarray_linalg;
 //! use ndarray::prelude::*;
 //! use ndarray_linalg::*;
-//! # fn main() {
 //!
-//! let a: Array2<f64> = random((3, 3));
+//! /// Use fixed algorithm and seed of PRNG for reproducible test
+//! let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+//!
+//! let a: Array2<f64> = random_using((3, 3), &mut rng);
 //! let f = a.factorizeh_into().unwrap(); // Factorize A (A is consumed)
 //! for _ in 0..10 {
-//!     let b: Array1<f64> = random(3);
+//!     let b: Array1<f64> = random_using(3, &mut rng);
 //!     let x = f.solveh_into(b).unwrap(); // Solve A * x = b using the factorization
 //! }
-//!
-//! # }
 //! ```
 
 use ndarray::*;

ndarray-linalg/tests/arnoldi.rs

@@ -3,8 +3,9 @@ use ndarray_linalg::{krylov::*, *};
 
 #[test]
 fn aq_qh_mgs() {
-    let a: Array2<f64> = random((5, 5));
-    let v: Array1<f64> = random(5);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<f64> = random_using((5, 5), &mut rng);
+    let v: Array1<f64> = random_using(5, &mut rng);
     let (q, h) = arnoldi_mgs(a.clone(), v, 1e-9);
     println!("A = \n{:?}", &a);
     println!("Q = \n{:?}", &q);
@@ -18,8 +19,9 @@ fn aq_qh_mgs() {
 
 #[test]
 fn aq_qh_householder() {
-    let a: Array2<f64> = random((5, 5));
-    let v: Array1<f64> = random(5);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<f64> = random_using((5, 5), &mut rng);
+    let v: Array1<f64> = random_using(5, &mut rng);
     let (q, h) = arnoldi_mgs(a.clone(), v, 1e-9);
     println!("A = \n{:?}", &a);
     println!("Q = \n{:?}", &q);
@@ -33,8 +35,9 @@ fn aq_qh_householder() {
 
 #[test]
 fn aq_qh_mgs_complex() {
-    let a: Array2<c64> = random((5, 5));
-    let v: Array1<c64> = random(5);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<c64> = random_using((5, 5), &mut rng);
+    let v: Array1<c64> = random_using(5, &mut rng);
     let (q, h) = arnoldi_mgs(a.clone(), v, 1e-9);
     println!("A = \n{:?}", &a);
     println!("Q = \n{:?}", &q);
@@ -48,8 +51,9 @@ fn aq_qh_mgs_complex() {
 
 #[test]
 fn aq_qh_householder_complex() {
-    let a: Array2<c64> = random((5, 5));
-    let v: Array1<c64> = random(5);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<c64> = random_using((5, 5), &mut rng);
+    let v: Array1<c64> = random_using(5, &mut rng);
     let (q, h) = arnoldi_mgs(a.clone(), v, 1e-9);
     println!("A = \n{:?}", &a);
     println!("Q = \n{:?}", &q);

ndarray-linalg/tests/cholesky.rs

@@ -6,7 +6,8 @@ macro_rules! cholesky {
         paste::item! {
             #[test]
             fn [<cholesky_ $elem>]() {
-                let a_orig: Array2<$elem> = random_hpd(3);
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a_orig: Array2<$elem> = random_hpd_using(3, &mut rng);
                 println!("a = \n{:?}", a_orig);
 
                 let upper = a_orig.cholesky(UPLO::Upper).unwrap();
@@ -79,7 +80,8 @@ macro_rules! cholesky_into_lower_upper {
         paste::item! {
             #[test]
             fn [<cholesky_into_lower_upper_ $elem>]() {
-                let a: Array2<$elem> = random_hpd(3);
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a: Array2<$elem> = random_hpd_using(3, &mut rng);
                 println!("a = \n{:?}", a);
                 let upper = a.cholesky(UPLO::Upper).unwrap();
                 let fac_upper = a.factorizec(UPLO::Upper).unwrap();
@@ -106,7 +108,8 @@ macro_rules! cholesky_into_inverse {
         paste::item! {
             #[test]
             fn [<cholesky_inverse_ $elem>]() {
-                let a: Array2<$elem> = random_hpd(3);
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a: Array2<$elem> = random_hpd_using(3, &mut rng);
                 println!("a = \n{:?}", a);
                 let inv = a.invc().unwrap();
                 assert_close_l2!(&a.dot(&inv), &Array2::eye(3), $rtol);
@@ -134,7 +137,8 @@ macro_rules! cholesky_det {
         paste::item! {
             #[test]
             fn [<cholesky_det_ $elem>]() {
-                let a: Array2<$elem> = random_hpd(3);
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a: Array2<$elem> = random_hpd_using(3, &mut rng);
                 println!("a = \n{:?}", a);
                 let ln_det = a
                     .eigvalsh(UPLO::Upper)
@@ -168,8 +172,9 @@ macro_rules! cholesky_solve {
         paste::item! {
             #[test]
             fn [<cholesky_solve_ $elem>]() {
-                let a: Array2<$elem> = random_hpd(3);
-                let x: Array1<$elem> = random(3);
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a: Array2<$elem> = random_hpd_using(3, &mut rng);
+                let x: Array1<$elem> = random_using(3, &mut rng);
                 let b = a.dot(&x);
                 println!("a = \n{:?}", a);
                 println!("x = \n{:?}", x);

ndarray-linalg/tests/convert.rs

@@ -3,7 +3,8 @@ use ndarray_linalg::*;
 
 #[test]
 fn generalize() {
-    let a: Array3<f64> = random((3, 2, 4).f());
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array3<f64> = random_using((3, 2, 4).f(), &mut rng);
     let ans = a.clone();
     let a: Array3<f64> = convert::generalize(a);
     assert_eq!(a, ans);

ndarray-linalg/tests/det.rs

@@ -136,23 +136,45 @@ fn det() {
             assert_rclose!(result.1, ln_det, rtol);
         }
     }
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
     for rows in 1..5 {
-        det_impl(random_regular::<f64>(rows), 1e-9);
-        det_impl(random_regular::<f32>(rows), 1e-4);
-        det_impl(random_regular::<c64>(rows), 1e-9);
-        det_impl(random_regular::<c32>(rows), 1e-4);
-        det_impl(random_regular::<f64>(rows).t().to_owned(), 1e-9);
-        det_impl(random_regular::<f32>(rows).t().to_owned(), 1e-4);
-        det_impl(random_regular::<c64>(rows).t().to_owned(), 1e-9);
-        det_impl(random_regular::<c32>(rows).t().to_owned(), 1e-4);
+        det_impl(random_regular_using::<f64, _>(rows, &mut rng), 1e-9);
+        det_impl(random_regular_using::<f32, _>(rows, &mut rng), 1e-4);
+        det_impl(random_regular_using::<c64, _>(rows, &mut rng), 1e-9);
+        det_impl(random_regular_using::<c32, _>(rows, &mut rng), 1e-4);
+        det_impl(
+            random_regular_using::<f64, _>(rows, &mut rng)
+                .t()
+                .to_owned(),
+            1e-9,
+        );
+        det_impl(
+            random_regular_using::<f32, _>(rows, &mut rng)
+                .t()
+                .to_owned(),
+            1e-4,
+        );
+        det_impl(
+            random_regular_using::<c64, _>(rows, &mut rng)
+                .t()
+                .to_owned(),
+            1e-9,
+        );
+        det_impl(
+            random_regular_using::<c32, _>(rows, &mut rng)
+                .t()
+                .to_owned(),
+            1e-4,
+        );
     }
 }
 
 #[test]
 fn det_nonsquare() {
     macro_rules! det_nonsquare {
         ($elem:ty, $shape:expr) => {
-            let a: Array2<$elem> = random($shape);
+            let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+            let a: Array2<$elem> = random_using($shape, &mut rng);
             assert!(a.factorize().unwrap().det().is_err());
             assert!(a.factorize().unwrap().sln_det().is_err());
             assert!(a.factorize().unwrap().det_into().is_err());

ndarray-linalg/tests/deth.rs

@@ -72,7 +72,8 @@ fn deth_zero_nonsquare() {
 fn deth() {
     macro_rules! deth {
         ($elem:ty, $rows:expr, $atol:expr) => {
-            let a: Array2<$elem> = random_hermite($rows);
+            let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+            let a: Array2<$elem> = random_hermite_using($rows, &mut rng);
             println!("a = \n{:?}", a);
 
             // Compute determinant from eigenvalues.

ndarray-linalg/tests/eigh.rs

@@ -79,7 +79,8 @@ fn fixed_t_lower() {
 
 #[test]
 fn ssqrt() {
-    let a: Array2<f64> = random_hpd(3);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<f64> = random_hpd_using(3, &mut rng);
     let ans = a.clone();
     let s = a.ssqrt(UPLO::Upper).unwrap();
     println!("a = {:?}", &ans);
@@ -92,7 +93,8 @@ fn ssqrt() {
 
 #[test]
 fn ssqrt_t() {
-    let a: Array2<f64> = random_hpd(3).reversed_axes();
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<f64> = random_hpd_using(3, &mut rng).reversed_axes();
     let ans = a.clone();
     let s = a.ssqrt(UPLO::Upper).unwrap();
     println!("a = {:?}", &ans);
@@ -105,7 +107,8 @@ fn ssqrt_t() {
 
 #[test]
 fn ssqrt_lower() {
-    let a: Array2<f64> = random_hpd(3);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<f64> = random_hpd_using(3, &mut rng);
     let ans = a.clone();
     let s = a.ssqrt(UPLO::Lower).unwrap();
     println!("a = {:?}", &ans);
@@ -118,7 +121,8 @@ fn ssqrt_lower() {
 
 #[test]
 fn ssqrt_t_lower() {
-    let a: Array2<f64> = random_hpd(3).reversed_axes();
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<f64> = random_hpd_using(3, &mut rng).reversed_axes();
     let ans = a.clone();
     let s = a.ssqrt(UPLO::Lower).unwrap();
     println!("a = {:?}", &ans);