From b1fe48f9be91cea1520f5b4605958640cf88f1b8 Mon Sep 17 00:00:00 2001
From: Toshiki Teramura <toshiki.teramura@gmail.com>
Date: Sat, 17 Sep 2022 00:50:34 +0900
Subject: [PATCH 1/4] Add rand_pcg as dev-dependency

---
 ndarray-linalg/Cargo.toml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/ndarray-linalg/Cargo.toml b/ndarray-linalg/Cargo.toml
index bfab605f..c7006cb2 100644
--- a/ndarray-linalg/Cargo.toml
+++ b/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"

From 0470a82fc8dc85b0d5a4e088ac8079d2b9fe4502 Mon Sep 17 00:00:00 2001
From: Toshiki Teramura <toshiki.teramura@gmail.com>
Date: Sat, 17 Sep 2022 00:56:58 +0900
Subject: [PATCH 2/4] Use random_*_using API in ndarray_linalg/tests

---
 ndarray-linalg/tests/arnoldi.rs            | 20 +++++---
 ndarray-linalg/tests/cholesky.rs           | 17 +++---
 ndarray-linalg/tests/convert.rs            |  3 +-
 ndarray-linalg/tests/det.rs                | 40 +++++++++++----
 ndarray-linalg/tests/deth.rs               |  3 +-
 ndarray-linalg/tests/eigh.rs               | 12 +++--
 ndarray-linalg/tests/householder.rs        |  9 ++--
 ndarray-linalg/tests/inner.rs              |  3 +-
 ndarray-linalg/tests/inv.rs                |  6 ++-
 ndarray-linalg/tests/least_squares.rs      | 27 ++++++----
 ndarray-linalg/tests/least_squares_nrhs.rs | 60 +++++++++++++---------
 ndarray-linalg/tests/mgs.rs                |  9 ++--
 ndarray-linalg/tests/normalize.rs          |  6 ++-
 ndarray-linalg/tests/qr.rs                 | 24 ++++++---
 ndarray-linalg/tests/solve.rs              | 56 ++++++++++++--------
 ndarray-linalg/tests/solveh.rs             | 30 ++++++-----
 ndarray-linalg/tests/svd.rs                |  6 ++-
 ndarray-linalg/tests/svddc.rs              | 18 ++++---
 ndarray-linalg/tests/trace.rs              |  3 +-
 ndarray-linalg/tests/triangular.rs         | 60 +++++++++++++---------
 ndarray-linalg/tests/tridiagonal.rs        | 21 +++++---
 21 files changed, 277 insertions(+), 156 deletions(-)

diff --git a/ndarray-linalg/tests/arnoldi.rs b/ndarray-linalg/tests/arnoldi.rs
index bbc15553..dd56e0a0 100644
--- a/ndarray-linalg/tests/arnoldi.rs
+++ b/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);
diff --git a/ndarray-linalg/tests/cholesky.rs b/ndarray-linalg/tests/cholesky.rs
index a498afc3..d3e9942b 100644
--- a/ndarray-linalg/tests/cholesky.rs
+++ b/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);
diff --git a/ndarray-linalg/tests/convert.rs b/ndarray-linalg/tests/convert.rs
index 3a6155d8..1e20d916 100644
--- a/ndarray-linalg/tests/convert.rs
+++ b/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);
diff --git a/ndarray-linalg/tests/det.rs b/ndarray-linalg/tests/det.rs
index c3986528..d14fc8d0 100644
--- a/ndarray-linalg/tests/det.rs
+++ b/ndarray-linalg/tests/det.rs
@@ -136,15 +136,36 @@ 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,
+        );
     }
 }
 
@@ -152,7 +173,8 @@ fn det() {
 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());
diff --git a/ndarray-linalg/tests/deth.rs b/ndarray-linalg/tests/deth.rs
index abd54105..4785aadc 100644
--- a/ndarray-linalg/tests/deth.rs
+++ b/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.
diff --git a/ndarray-linalg/tests/eigh.rs b/ndarray-linalg/tests/eigh.rs
index dd1445be..8d8ce385 100644
--- a/ndarray-linalg/tests/eigh.rs
+++ b/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);
diff --git a/ndarray-linalg/tests/householder.rs b/ndarray-linalg/tests/householder.rs
index adc4d9b2..83b500f7 100644
--- a/ndarray-linalg/tests/householder.rs
+++ b/ndarray-linalg/tests/householder.rs
@@ -3,7 +3,8 @@ use ndarray_linalg::{krylov::*, *};
 
 fn over<A: Scalar + Lapack>(rtol: A::Real) {
     const N: usize = 4;
-    let a: Array2<A> = random((N, N * 2));
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<A> = random_using((N, N * 2), &mut rng);
 
     // Terminate
     let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate);
@@ -45,7 +46,8 @@ fn over_c64() {
 
 fn full<A: Scalar + Lapack>(rtol: A::Real) {
     const N: usize = 5;
-    let a: Array2<A> = random((N, N));
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<A> = random_using((N, N), &mut rng);
     let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate);
     let qc: Array2<A> = conjugate(&q);
     assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol; "Check Q^H Q = I");
@@ -71,7 +73,8 @@ fn full_c64() {
 
 fn half<A: Scalar + Lapack>(rtol: A::Real) {
     const N: usize = 4;
-    let a: Array2<A> = random((N, N / 2));
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<A> = random_using((N, N / 2), &mut rng);
     let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate);
     let qc: Array2<A> = conjugate(&q);
     assert_close_l2!(&qc.dot(&q), &Array::eye(N / 2), rtol; "Check Q^H Q = I");
diff --git a/ndarray-linalg/tests/inner.rs b/ndarray-linalg/tests/inner.rs
index 7fc42c83..076b2791 100644
--- a/ndarray-linalg/tests/inner.rs
+++ b/ndarray-linalg/tests/inner.rs
@@ -19,7 +19,8 @@ fn size_longer() {
 
 #[test]
 fn abs() {
-    let a: Array1<c32> = random(1);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array1<c32> = random_using(1, &mut rng);
     let aa = a.inner(&a);
     assert_aclose!(aa.re(), a.norm().powi(2), 1e-5);
     assert_aclose!(aa.im(), 0.0, 1e-5);
diff --git a/ndarray-linalg/tests/inv.rs b/ndarray-linalg/tests/inv.rs
index 71e8973a..030773c1 100644
--- a/ndarray-linalg/tests/inv.rs
+++ b/ndarray-linalg/tests/inv.rs
@@ -5,7 +5,8 @@ fn test_inv_random<A>(n: usize, set_f: bool, rtol: A::Real)
 where
     A: Scalar + Lapack,
 {
-    let a: Array2<A> = random([n; 2].set_f(set_f));
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<A> = random_using([n; 2].set_f(set_f), &mut rng);
     let identity = Array2::eye(n);
     assert_close_l2!(&a.inv().unwrap().dot(&a), &identity, rtol);
     assert_close_l2!(
@@ -24,7 +25,8 @@ fn test_inv_into_random<A>(n: usize, set_f: bool, rtol: A::Real)
 where
     A: Scalar + Lapack,
 {
-    let a: Array2<A> = random([n; 2].set_f(set_f));
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<A> = random_using([n; 2].set_f(set_f), &mut rng);
     let identity = Array2::eye(n);
     assert_close_l2!(&a.clone().inv_into().unwrap().dot(&a), &identity, rtol);
     assert_close_l2!(
diff --git a/ndarray-linalg/tests/least_squares.rs b/ndarray-linalg/tests/least_squares.rs
index 33e20ca7..17f62b64 100644
--- a/ndarray-linalg/tests/least_squares.rs
+++ b/ndarray-linalg/tests/least_squares.rs
@@ -5,7 +5,8 @@ use ndarray_linalg::*;
 
 /// A is square. `x = A^{-1} b`, `|b - Ax| = 0`
 fn test_exact<T: Scalar + Lapack>(a: Array2<T>) {
-    let b: Array1<T> = random(3);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let b: Array1<T> = random_using(3, &mut rng);
     let result = a.least_squares(&b).unwrap();
     // unpack result
     let x = result.solution;
@@ -27,13 +28,15 @@ macro_rules! impl_exact {
         paste::item! {
             #[test]
             fn [<least_squares_ $scalar _exact>]() {
-                let a: Array2<$scalar> = random((3, 3));
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a: Array2<$scalar> = random_using((3, 3), &mut rng);
                 test_exact(a)
             }
 
             #[test]
             fn [<least_squares_ $scalar _exact_t>]() {
-                let a: Array2<$scalar> = random((3, 3).f());
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a: Array2<$scalar> = random_using((3, 3).f(), &mut rng);
                 test_exact(a)
             }
         }
@@ -51,7 +54,8 @@ fn test_overdetermined<T: Scalar + Lapack>(a: Array2<T>)
 where
     T::Real: AbsDiffEq<Epsilon = T::Real>,
 {
-    let b: Array1<T> = random(4);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let b: Array1<T> = random_using(4, &mut rng);
     let result = a.least_squares(&b).unwrap();
     // unpack result
     let x = result.solution;
@@ -73,13 +77,15 @@ macro_rules! impl_overdetermined {
         paste::item! {
             #[test]
             fn [<least_squares_ $scalar _overdetermined>]() {
-                let a: Array2<$scalar> = random((4, 3));
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a: Array2<$scalar> = random_using((4, 3), &mut rng);
                 test_overdetermined(a)
             }
 
             #[test]
             fn [<least_squares_ $scalar _overdetermined_t>]() {
-                let a: Array2<$scalar> = random((4, 3).f());
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a: Array2<$scalar> = random_using((4, 3).f(), &mut rng);
                 test_overdetermined(a)
             }
         }
@@ -94,7 +100,8 @@ impl_overdetermined!(c64);
 /// #column > #row case.
 /// Linear problem is underdetermined, `|b - Ax| = 0` and `x` is not unique
 fn test_underdetermined<T: Scalar + Lapack>(a: Array2<T>) {
-    let b: Array1<T> = random(3);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let b: Array1<T> = random_using(3, &mut rng);
     let result = a.least_squares(&b).unwrap();
     assert_eq!(result.rank, 3);
     assert!(result.residual_sum_of_squares.is_none());
@@ -110,13 +117,15 @@ macro_rules! impl_underdetermined {
         paste::item! {
             #[test]
             fn [<least_squares_ $scalar _underdetermined>]() {
-                let a: Array2<$scalar> = random((3, 4));
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a: Array2<$scalar> = random_using((3, 4), &mut rng);
                 test_underdetermined(a)
             }
 
             #[test]
             fn [<least_squares_ $scalar _underdetermined_t>]() {
-                let a: Array2<$scalar> = random((3, 4).f());
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a: Array2<$scalar> = random_using((3, 4).f(), &mut rng);
                 test_underdetermined(a)
             }
         }
diff --git a/ndarray-linalg/tests/least_squares_nrhs.rs b/ndarray-linalg/tests/least_squares_nrhs.rs
index dd7d283c..bcf6d013 100644
--- a/ndarray-linalg/tests/least_squares_nrhs.rs
+++ b/ndarray-linalg/tests/least_squares_nrhs.rs
@@ -32,29 +32,33 @@ macro_rules! impl_exact {
         paste::item! {
             #[test]
             fn [<least_squares_ $scalar _exact_ac_bc>]() {
-                let a: Array2<$scalar> = random((3, 3));
-                let b: Array2<$scalar> = random((3, 2));
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a: Array2<$scalar> = random_using((3, 3), &mut rng);
+                let b: Array2<$scalar> = random_using((3, 2), &mut rng);
                 test_exact(a, b)
             }
 
             #[test]
             fn [<least_squares_ $scalar _exact_ac_bf>]() {
-                let a: Array2<$scalar> = random((3, 3));
-                let b: Array2<$scalar> = random((3, 2).f());
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a: Array2<$scalar> = random_using((3, 3), &mut rng);
+                let b: Array2<$scalar> = random_using((3, 2).f(), &mut rng);
                 test_exact(a, b)
             }
 
             #[test]
             fn [<least_squares_ $scalar _exact_af_bc>]() {
-                let a: Array2<$scalar> = random((3, 3).f());
-                let b: Array2<$scalar> = random((3, 2));
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a: Array2<$scalar> = random_using((3, 3).f(), &mut rng);
+                let b: Array2<$scalar> = random_using((3, 2), &mut rng);
                 test_exact(a, b)
             }
 
             #[test]
             fn [<least_squares_ $scalar _exact_af_bf>]() {
-                let a: Array2<$scalar> = random((3, 3).f());
-                let b: Array2<$scalar> = random((3, 2).f());
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a: Array2<$scalar> = random_using((3, 3).f(), &mut rng);
+                let b: Array2<$scalar> = random_using((3, 2).f(), &mut rng);
                 test_exact(a, b)
             }
         }
@@ -100,29 +104,33 @@ macro_rules! impl_overdetermined {
         paste::item! {
             #[test]
             fn [<least_squares_ $scalar _overdetermined_ac_bc>]() {
-                let a: Array2<$scalar> = random((4, 3));
-                let b: Array2<$scalar> = random((4, 2));
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a: Array2<$scalar> = random_using((4, 3), &mut rng);
+                let b: Array2<$scalar> = random_using((4, 2), &mut rng);
                 test_overdetermined(a, b)
             }
 
             #[test]
             fn [<least_squares_ $scalar _overdetermined_af_bc>]() {
-                let a: Array2<$scalar> = random((4, 3).f());
-                let b: Array2<$scalar> = random((4, 2));
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a: Array2<$scalar> = random_using((4, 3).f(), &mut rng);
+                let b: Array2<$scalar> = random_using((4, 2), &mut rng);
                 test_overdetermined(a, b)
             }
 
             #[test]
             fn [<least_squares_ $scalar _overdetermined_ac_bf>]() {
-                let a: Array2<$scalar> = random((4, 3));
-                let b: Array2<$scalar> = random((4, 2).f());
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a: Array2<$scalar> = random_using((4, 3), &mut rng);
+                let b: Array2<$scalar> = random_using((4, 2).f(), &mut rng);
                 test_overdetermined(a, b)
             }
 
             #[test]
             fn [<least_squares_ $scalar _overdetermined_af_bf>]() {
-                let a: Array2<$scalar> = random((4, 3).f());
-                let b: Array2<$scalar> = random((4, 2).f());
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a: Array2<$scalar> = random_using((4, 3).f(), &mut rng);
+                let b: Array2<$scalar> = random_using((4, 2).f(), &mut rng);
                 test_overdetermined(a, b)
             }
         }
@@ -155,29 +163,33 @@ macro_rules! impl_underdetermined {
         paste::item! {
             #[test]
             fn [<least_squares_ $scalar _underdetermined_ac_bc>]() {
-                let a: Array2<$scalar> = random((3, 4));
-                let b: Array2<$scalar> = random((3, 2));
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a: Array2<$scalar> = random_using((3, 4), &mut rng);
+                let b: Array2<$scalar> = random_using((3, 2), &mut rng);
                 test_underdetermined(a, b)
             }
 
             #[test]
             fn [<least_squares_ $scalar _underdetermined_af_bc>]() {
-                let a: Array2<$scalar> = random((3, 4).f());
-                let b: Array2<$scalar> = random((3, 2));
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a: Array2<$scalar> = random_using((3, 4).f(), &mut rng);
+                let b: Array2<$scalar> = random_using((3, 2), &mut rng);
                 test_underdetermined(a, b)
             }
 
             #[test]
             fn [<least_squares_ $scalar _underdetermined_ac_bf>]() {
-                let a: Array2<$scalar> = random((3, 4));
-                let b: Array2<$scalar> = random((3, 2).f());
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a: Array2<$scalar> = random_using((3, 4), &mut rng);
+                let b: Array2<$scalar> = random_using((3, 2).f(), &mut rng);
                 test_underdetermined(a, b)
             }
 
             #[test]
             fn [<least_squares_ $scalar _underdetermined_af_bf>]() {
-                let a: Array2<$scalar> = random((3, 4).f());
-                let b: Array2<$scalar> = random((3, 2).f());
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a: Array2<$scalar> = random_using((3, 4).f(), &mut rng);
+                let b: Array2<$scalar> = random_using((3, 2).f(), &mut rng);
                 test_underdetermined(a, b)
             }
         }
diff --git a/ndarray-linalg/tests/mgs.rs b/ndarray-linalg/tests/mgs.rs
index 35c860de..9e9aa29e 100644
--- a/ndarray-linalg/tests/mgs.rs
+++ b/ndarray-linalg/tests/mgs.rs
@@ -5,7 +5,8 @@ fn qr_full<A: Scalar + Lapack>() {
     const N: usize = 5;
     let rtol: A::Real = A::real(1e-9);
 
-    let a: Array2<A> = random((N, N));
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<A> = random_using((N, N), &mut rng);
     let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate);
     assert_close_l2!(&q.dot(&r), &a, rtol);
 
@@ -27,7 +28,8 @@ fn qr<A: Scalar + Lapack>() {
     const N: usize = 4;
     let rtol: A::Real = A::real(1e-9);
 
-    let a: Array2<A> = random((N, N / 2));
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<A> = random_using((N, N / 2), &mut rng);
     let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate);
     assert_close_l2!(&q.dot(&r), &a, rtol);
 
@@ -49,7 +51,8 @@ fn qr_over<A: Scalar + Lapack>() {
     const N: usize = 4;
     let rtol: A::Real = A::real(1e-9);
 
-    let a: Array2<A> = random((N, N * 2));
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<A> = random_using((N, N * 2), &mut rng);
 
     // Terminate
     let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate);
diff --git a/ndarray-linalg/tests/normalize.rs b/ndarray-linalg/tests/normalize.rs
index ca50912e..8d71a009 100644
--- a/ndarray-linalg/tests/normalize.rs
+++ b/ndarray-linalg/tests/normalize.rs
@@ -3,14 +3,16 @@ use ndarray_linalg::*;
 
 #[test]
 fn n_columns() {
-    let a: Array2<f64> = random((3, 2));
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<f64> = random_using((3, 2), &mut rng);
     let (n, v) = normalize(a.clone(), NormalizeAxis::Column);
     assert_close_l2!(&n.dot(&from_diag(&v)), &a, 1e-7);
 }
 
 #[test]
 fn n_rows() {
-    let a: Array2<f64> = random((3, 2));
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<f64> = random_using((3, 2), &mut rng);
     let (n, v) = normalize(a.clone(), NormalizeAxis::Row);
     assert_close_l2!(&from_diag(&v).dot(&n), &a, 1e-7);
 }
diff --git a/ndarray-linalg/tests/qr.rs b/ndarray-linalg/tests/qr.rs
index a69d89e2..702ed060 100644
--- a/ndarray-linalg/tests/qr.rs
+++ b/ndarray-linalg/tests/qr.rs
@@ -26,48 +26,56 @@ fn test_square(a: &Array2<f64>, n: usize, m: usize) {
 
 #[test]
 fn qr_sq() {
-    let a = random((3, 3));
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a = random_using((3, 3), &mut rng);
     test_square(&a, 3, 3);
 }
 
 #[test]
 fn qr_sq_t() {
-    let a = random((3, 3).f());
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a = random_using((3, 3).f(), &mut rng);
     test_square(&a, 3, 3);
 }
 
 #[test]
 fn qr_3x3() {
-    let a = random((3, 3));
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a = random_using((3, 3), &mut rng);
     test(&a, 3, 3);
 }
 
 #[test]
 fn qr_3x3_t() {
-    let a = random((3, 3).f());
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a = random_using((3, 3).f(), &mut rng);
     test(&a, 3, 3);
 }
 
 #[test]
 fn qr_3x4() {
-    let a = random((3, 4));
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a = random_using((3, 4), &mut rng);
     test(&a, 3, 4);
 }
 
 #[test]
 fn qr_3x4_t() {
-    let a = random((3, 4).f());
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a = random_using((3, 4).f(), &mut rng);
     test(&a, 3, 4);
 }
 
 #[test]
 fn qr_4x3() {
-    let a = random((4, 3));
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a = random_using((4, 3), &mut rng);
     test(&a, 4, 3);
 }
 
 #[test]
 fn qr_4x3_t() {
-    let a = random((4, 3).f());
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a = random_using((4, 3).f(), &mut rng);
     test(&a, 4, 3);
 }
diff --git a/ndarray-linalg/tests/solve.rs b/ndarray-linalg/tests/solve.rs
index b3d0cb53..57c29b67 100644
--- a/ndarray-linalg/tests/solve.rs
+++ b/ndarray-linalg/tests/solve.rs
@@ -1,6 +1,7 @@
 use ndarray::prelude::*;
 use ndarray_linalg::{
-    assert_aclose, assert_close_l2, c32, c64, random, random_hpd, solve::*, OperationNorm, Scalar,
+    assert_aclose, assert_close_l2, c32, c64, random_hpd_using, random_using, solve::*,
+    OperationNorm, Scalar,
 };
 
 macro_rules! test_solve {
@@ -97,12 +98,13 @@ macro_rules! test_solve_all {
 
 #[test]
 fn solve_random_float() {
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
     for n in 0..=8 {
         for &set_f in &[false, true] {
             test_solve_all!(
                 [f32 => 1e-3, f64 => 1e-9],
-                a = random([n; 2].set_f(set_f)),
-                x = random(n),
+                a = random_using([n; 2].set_f(set_f), &mut rng),
+                x = random_using(n, &mut rng),
                 b = a.dot(&x),
                 [solve, solve_into, solve_inplace],
             );
@@ -112,12 +114,13 @@ fn solve_random_float() {
 
 #[test]
 fn solve_random_complex() {
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
     for n in 0..=8 {
         for &set_f in &[false, true] {
             test_solve_all!(
                 [c32 => 1e-3, c64 => 1e-9],
-                a = random([n; 2].set_f(set_f)),
-                x = random(n),
+                a = random_using([n; 2].set_f(set_f), &mut rng),
+                x = random_using(n, &mut rng),
                 b = a.dot(&x),
                 [solve, solve_into, solve_inplace],
             );
@@ -128,19 +131,21 @@ fn solve_random_complex() {
 #[should_panic]
 #[test]
 fn solve_shape_mismatch() {
-    let a: Array2<f64> = random((3, 3));
-    let b: Array1<f64> = random(2);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<f64> = random_using((3, 3), &mut rng);
+    let b: Array1<f64> = random_using(2, &mut rng);
     let _ = a.solve_into(b);
 }
 
 #[test]
 fn solve_t_random_float() {
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
     for n in 0..=8 {
         for &set_f in &[false, true] {
             test_solve_all!(
                 [f32 => 1e-3, f64 => 1e-9],
-                a = random([n; 2].set_f(set_f)),
-                x = random(n),
+                a = random_using([n; 2].set_f(set_f), &mut rng),
+                x = random_using(n, &mut rng),
                 b = a.t().dot(&x),
                 [solve_t, solve_t_into, solve_t_inplace],
             );
@@ -151,19 +156,21 @@ fn solve_t_random_float() {
 #[should_panic]
 #[test]
 fn solve_t_shape_mismatch() {
-    let a: Array2<f64> = random((3, 3).f());
-    let b: Array1<f64> = random(4);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<f64> = random_using((3, 3).f(), &mut rng);
+    let b: Array1<f64> = random_using(4, &mut rng);
     let _ = a.solve_into(b);
 }
 
 #[test]
 fn solve_t_random_complex() {
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
     for n in 0..=8 {
         for &set_f in &[false, true] {
             test_solve_all!(
                 [c32 => 1e-3, c64 => 1e-9],
-                a = random([n; 2].set_f(set_f)),
-                x = random(n),
+                a = random_using([n; 2].set_f(set_f), &mut rng),
+                x = random_using(n, &mut rng),
                 b = a.t().dot(&x),
                 [solve_t, solve_t_into, solve_t_inplace],
             );
@@ -174,20 +181,22 @@ fn solve_t_random_complex() {
 #[should_panic]
 #[test]
 fn solve_factorized_shape_mismatch() {
-    let a: Array2<f64> = random((3, 3));
-    let b: Array1<f64> = random(4);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<f64> = random_using((3, 3), &mut rng);
+    let b: Array1<f64> = random_using(4, &mut rng);
     let f = a.factorize_into().unwrap();
     let _ = f.solve_into(b);
 }
 
 #[test]
 fn solve_h_random_float() {
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
     for n in 0..=8 {
         for &set_f in &[false, true] {
             test_solve_all!(
                 [f32 => 1e-3, f64 => 1e-9],
-                a = random([n; 2].set_f(set_f)),
-                x = random(n),
+                a = random_using([n; 2].set_f(set_f), &mut rng),
+                x = random_using(n, &mut rng),
                 b = a.t().mapv(|x| x.conj()).dot(&x),
                 [solve_h, solve_h_into, solve_h_inplace],
             );
@@ -198,20 +207,22 @@ fn solve_h_random_float() {
 #[should_panic]
 #[test]
 fn solve_factorized_t_shape_mismatch() {
-    let a: Array2<f64> = random((3, 3).f());
-    let b: Array1<f64> = random(4);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<f64> = random_using((3, 3).f(), &mut rng);
+    let b: Array1<f64> = random_using(4, &mut rng);
     let f = a.factorize_into().unwrap();
     let _ = f.solve_into(b);
 }
 
 #[test]
 fn solve_h_random_complex() {
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
     for n in 0..=8 {
         for &set_f in &[false, true] {
             test_solve_all!(
                 [c32 => 1e-3, c64 => 1e-9],
-                a = random([n; 2].set_f(set_f)),
-                x = random(n),
+                a = random_using([n; 2].set_f(set_f), &mut rng),
+                x = random_using(n, &mut rng),
                 b = a.t().mapv(|x| x.conj()).dot(&x),
                 [solve_h, solve_h_into, solve_h_inplace],
             );
@@ -223,7 +234,8 @@ fn solve_h_random_complex() {
 fn rcond() {
     macro_rules! rcond {
         ($elem:ty, $rows:expr, $atol:expr) => {
-            let a: Array2<$elem> = random_hpd($rows);
+            let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+            let a: Array2<$elem> = random_hpd_using($rows, &mut rng);
             let rcond = 1. / (a.opnorm_one().unwrap() * a.inv().unwrap().opnorm_one().unwrap());
             assert_aclose!(a.rcond().unwrap(), rcond, $atol);
             assert_aclose!(a.rcond_into().unwrap(), rcond, $atol);
diff --git a/ndarray-linalg/tests/solveh.rs b/ndarray-linalg/tests/solveh.rs
index 56d016c4..25513551 100644
--- a/ndarray-linalg/tests/solveh.rs
+++ b/ndarray-linalg/tests/solveh.rs
@@ -4,24 +4,27 @@ use ndarray_linalg::*;
 #[should_panic]
 #[test]
 fn solveh_shape_mismatch() {
-    let a: Array2<f64> = random_hpd(3);
-    let b: Array1<f64> = random(2);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<f64> = random_hpd_using(3, &mut rng);
+    let b: Array1<f64> = random_using(2, &mut rng);
     let _ = a.solveh_into(b);
 }
 
 #[should_panic]
 #[test]
 fn factorizeh_solveh_shape_mismatch() {
-    let a: Array2<f64> = random_hpd(3);
-    let b: Array1<f64> = random(2);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<f64> = random_hpd_using(3, &mut rng);
+    let b: Array1<f64> = random_using(2, &mut rng);
     let f = a.factorizeh_into().unwrap();
     let _ = f.solveh_into(b);
 }
 
 #[test]
 fn solveh_random() {
-    let a: Array2<f64> = random_hpd(3);
-    let x: Array1<f64> = random(3);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<f64> = random_hpd_using(3, &mut rng);
+    let x: Array1<f64> = random_using(3, &mut rng);
     let b = a.dot(&x);
     let y = a.solveh_into(b).unwrap();
     assert_close_l2!(&x, &y, 1e-7);
@@ -35,24 +38,27 @@ fn solveh_random() {
 #[should_panic]
 #[test]
 fn solveh_t_shape_mismatch() {
-    let a: Array2<f64> = random_hpd(3).reversed_axes();
-    let b: Array1<f64> = random(2);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<f64> = random_hpd_using(3, &mut rng).reversed_axes();
+    let b: Array1<f64> = random_using(2, &mut rng);
     let _ = a.solveh_into(b);
 }
 
 #[should_panic]
 #[test]
 fn factorizeh_solveh_t_shape_mismatch() {
-    let a: Array2<f64> = random_hpd(3).reversed_axes();
-    let b: Array1<f64> = random(2);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<f64> = random_hpd_using(3, &mut rng).reversed_axes();
+    let b: Array1<f64> = random_using(2, &mut rng);
     let f = a.factorizeh_into().unwrap();
     let _ = f.solveh_into(b);
 }
 
 #[test]
 fn solveh_random_t() {
-    let a: Array2<f64> = random_hpd(3).reversed_axes();
-    let x: Array1<f64> = random(3);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<f64> = random_hpd_using(3, &mut rng).reversed_axes();
+    let x: Array1<f64> = random_using(3, &mut rng);
     let b = a.dot(&x);
     let y = a.solveh_into(b).unwrap();
     assert_close_l2!(&x, &y, 1e-7);
diff --git a/ndarray-linalg/tests/svd.rs b/ndarray-linalg/tests/svd.rs
index c83885e1..0eac35ea 100644
--- a/ndarray-linalg/tests/svd.rs
+++ b/ndarray-linalg/tests/svd.rs
@@ -53,13 +53,15 @@ macro_rules! test_svd_impl {
         paste::item! {
             #[test]
             fn [<svd_ $type _ $test _ $n x $m>]() {
-                let a = random(($n, $m));
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a = random_using(($n, $m), &mut rng);
                 $test::<$type>(&a);
             }
 
             #[test]
             fn [<svd_ $type _ $test _ $n x $m _t>]() {
-                let a = random(($n, $m).f());
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a = random_using(($n, $m).f(), &mut rng);
                 $test::<$type>(&a);
             }
         }
diff --git a/ndarray-linalg/tests/svddc.rs b/ndarray-linalg/tests/svddc.rs
index ed28fc30..60c6bd66 100644
--- a/ndarray-linalg/tests/svddc.rs
+++ b/ndarray-linalg/tests/svddc.rs
@@ -32,37 +32,43 @@ macro_rules! test_svd_impl {
         paste::item! {
             #[test]
             fn [<svddc_ $scalar _full_ $n x $m>]() {
-                let a = random(($n, $m));
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a = random_using(($n, $m), &mut rng);
                 test::<$scalar>(&a, JobSvd::All);
             }
 
             #[test]
             fn [<svddc_ $scalar _some_ $n x $m>]() {
-                let a = random(($n, $m));
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a = random_using(($n, $m), &mut rng);
                 test::<$scalar>(&a, JobSvd::Some);
             }
 
             #[test]
             fn [<svddc_ $scalar _none_ $n x $m>]() {
-                let a = random(($n, $m));
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a = random_using(($n, $m), &mut rng);
                 test::<$scalar>(&a, JobSvd::None);
             }
 
             #[test]
             fn [<svddc_ $scalar _full_ $n x $m _t>]() {
-                let a = random(($n, $m).f());
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a = random_using(($n, $m).f(), &mut rng);
                 test::<$scalar>(&a, JobSvd::All);
             }
 
             #[test]
             fn [<svddc_ $scalar _some_ $n x $m _t>]() {
-                let a = random(($n, $m).f());
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a = random_using(($n, $m).f(), &mut rng);
                 test::<$scalar>(&a, JobSvd::Some);
             }
 
             #[test]
             fn [<svddc_ $scalar _none_ $n x $m _t>]() {
-                let a = random(($n, $m).f());
+                let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+                let a = random_using(($n, $m).f(), &mut rng);
                 test::<$scalar>(&a, JobSvd::None);
             }
         }
diff --git a/ndarray-linalg/tests/trace.rs b/ndarray-linalg/tests/trace.rs
index 9127be39..f93bb52b 100644
--- a/ndarray-linalg/tests/trace.rs
+++ b/ndarray-linalg/tests/trace.rs
@@ -3,6 +3,7 @@ use ndarray_linalg::*;
 
 #[test]
 fn trace() {
-    let a: Array2<f64> = random((3, 3));
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let a: Array2<f64> = random_using((3, 3), &mut rng);
     assert_rclose!(a.trace().unwrap(), a[(0, 0)] + a[(1, 1)] + a[(2, 2)], 1e-7);
 }
diff --git a/ndarray-linalg/tests/triangular.rs b/ndarray-linalg/tests/triangular.rs
index ca609c5c..4ddf0a8b 100644
--- a/ndarray-linalg/tests/triangular.rs
+++ b/ndarray-linalg/tests/triangular.rs
@@ -34,87 +34,99 @@ where
 #[test]
 fn triangular_1d_upper() {
     let n = 3;
-    let b: Array1<f64> = random(n);
-    let a: Array2<f64> = random((n, n)).into_triangular(UPLO::Upper);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let b: Array1<f64> = random_using(n, &mut rng);
+    let a: Array2<f64> = random_using((n, n), &mut rng).into_triangular(UPLO::Upper);
     test1d(UPLO::Upper, &a, &b, 1e-7);
 }
 
 #[test]
 fn triangular_1d_lower() {
     let n = 3;
-    let b: Array1<f64> = random(n);
-    let a: Array2<f64> = random((n, n)).into_triangular(UPLO::Lower);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let b: Array1<f64> = random_using(n, &mut rng);
+    let a: Array2<f64> = random_using((n, n), &mut rng).into_triangular(UPLO::Lower);
     test1d(UPLO::Lower, &a, &b, 1e-7);
 }
 
 #[test]
 fn triangular_1d_upper_t() {
     let n = 3;
-    let b: Array1<f64> = random(n);
-    let a: Array2<f64> = random((n, n).f()).into_triangular(UPLO::Upper);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let b: Array1<f64> = random_using(n, &mut rng);
+    let a: Array2<f64> = random_using((n, n).f(), &mut rng).into_triangular(UPLO::Upper);
     test1d(UPLO::Upper, &a, &b, 1e-7);
 }
 
 #[test]
 fn triangular_1d_lower_t() {
     let n = 3;
-    let b: Array1<f64> = random(n);
-    let a: Array2<f64> = random((n, n).f()).into_triangular(UPLO::Lower);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let b: Array1<f64> = random_using(n, &mut rng);
+    let a: Array2<f64> = random_using((n, n).f(), &mut rng).into_triangular(UPLO::Lower);
     test1d(UPLO::Lower, &a, &b, 1e-7);
 }
 
 #[test]
 fn triangular_2d_upper() {
-    let b: Array2<f64> = random((3, 4));
-    let a: Array2<f64> = random((3, 3)).into_triangular(UPLO::Upper);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let b: Array2<f64> = random_using((3, 4), &mut rng);
+    let a: Array2<f64> = random_using((3, 3), &mut rng).into_triangular(UPLO::Upper);
     test2d(UPLO::Upper, &a, &b, 1e-7);
 }
 
 #[test]
 fn triangular_2d_lower() {
-    let b: Array2<f64> = random((3, 4));
-    let a: Array2<f64> = random((3, 3)).into_triangular(UPLO::Lower);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let b: Array2<f64> = random_using((3, 4), &mut rng);
+    let a: Array2<f64> = random_using((3, 3), &mut rng).into_triangular(UPLO::Lower);
     test2d(UPLO::Lower, &a, &b, 1e-7);
 }
 
 #[test]
 fn triangular_2d_lower_t() {
-    let b: Array2<f64> = random((3, 4));
-    let a: Array2<f64> = random((3, 3).f()).into_triangular(UPLO::Lower);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let b: Array2<f64> = random_using((3, 4), &mut rng);
+    let a: Array2<f64> = random_using((3, 3).f(), &mut rng).into_triangular(UPLO::Lower);
     test2d(UPLO::Lower, &a, &b, 1e-7);
 }
 
 #[test]
 fn triangular_2d_upper_t() {
-    let b: Array2<f64> = random((3, 4));
-    let a: Array2<f64> = random((3, 3).f()).into_triangular(UPLO::Upper);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let b: Array2<f64> = random_using((3, 4), &mut rng);
+    let a: Array2<f64> = random_using((3, 3).f(), &mut rng).into_triangular(UPLO::Upper);
     test2d(UPLO::Upper, &a, &b, 1e-7);
 }
 
 #[test]
 fn triangular_2d_upper_bt() {
-    let b: Array2<f64> = random((3, 4).f());
-    let a: Array2<f64> = random((3, 3)).into_triangular(UPLO::Upper);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let b: Array2<f64> = random_using((3, 4).f(), &mut rng);
+    let a: Array2<f64> = random_using((3, 3), &mut rng).into_triangular(UPLO::Upper);
     test2d(UPLO::Upper, &a, &b, 1e-7);
 }
 
 #[test]
 fn triangular_2d_lower_bt() {
-    let b: Array2<f64> = random((3, 4).f());
-    let a: Array2<f64> = random((3, 3)).into_triangular(UPLO::Lower);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let b: Array2<f64> = random_using((3, 4).f(), &mut rng);
+    let a: Array2<f64> = random_using((3, 3), &mut rng).into_triangular(UPLO::Lower);
     test2d(UPLO::Lower, &a, &b, 1e-7);
 }
 
 #[test]
 fn triangular_2d_lower_t_bt() {
-    let b: Array2<f64> = random((3, 4).f());
-    let a: Array2<f64> = random((3, 3).f()).into_triangular(UPLO::Lower);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let b: Array2<f64> = random_using((3, 4).f(), &mut rng);
+    let a: Array2<f64> = random_using((3, 3).f(), &mut rng).into_triangular(UPLO::Lower);
     test2d(UPLO::Lower, &a, &b, 1e-7);
 }
 
 #[test]
 fn triangular_2d_upper_t_bt() {
-    let b: Array2<f64> = random((3, 4).f());
-    let a: Array2<f64> = random((3, 3).f()).into_triangular(UPLO::Upper);
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let b: Array2<f64> = random_using((3, 4).f(), &mut rng);
+    let a: Array2<f64> = random_using((3, 3).f(), &mut rng).into_triangular(UPLO::Upper);
     test2d(UPLO::Upper, &a, &b, 1e-7);
 }
diff --git a/ndarray-linalg/tests/tridiagonal.rs b/ndarray-linalg/tests/tridiagonal.rs
index 38278951..513d625b 100644
--- a/ndarray-linalg/tests/tridiagonal.rs
+++ b/ndarray-linalg/tests/tridiagonal.rs
@@ -28,7 +28,8 @@ fn tridiagonal_index() {
 
 #[test]
 fn opnorm_tridiagonal() {
-    let mut a: Array2<f64> = random((4, 4));
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let mut a: Array2<f64> = random_using((4, 4), &mut rng);
     a[[0, 2]] = 0.0;
     a[[0, 3]] = 0.0;
     a[[1, 3]] = 0.0;
@@ -129,10 +130,11 @@ fn solve_tridiagonal_c64() {
 
 #[test]
 fn solve_tridiagonal_random() {
-    let mut a: Array2<f64> = random((3, 3));
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let mut a: Array2<f64> = random_using((3, 3), &mut rng);
     a[[0, 2]] = 0.0;
     a[[2, 0]] = 0.0;
-    let x: Array1<f64> = random(3);
+    let x: Array1<f64> = random_using(3, &mut rng);
     let b1 = a.dot(&x);
     let b2 = b1.clone();
     let y1 = a.solve_tridiagonal_into(b1).unwrap();
@@ -143,10 +145,11 @@ fn solve_tridiagonal_random() {
 
 #[test]
 fn solve_tridiagonal_random_t() {
-    let mut a: Array2<f64> = random((3, 3));
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let mut a: Array2<f64> = random_using((3, 3), &mut rng);
     a[[0, 2]] = 0.0;
     a[[2, 0]] = 0.0;
-    let x: Array1<f64> = random(3);
+    let x: Array1<f64> = random_using(3, &mut rng);
     let at = a.t();
     let b1 = at.dot(&x);
     let b2 = b1.clone();
@@ -158,11 +161,12 @@ fn solve_tridiagonal_random_t() {
 
 #[test]
 fn extract_tridiagonal_solve_random() {
-    let mut a: Array2<f64> = random((3, 3));
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let mut a: Array2<f64> = random_using((3, 3), &mut rng);
     a[[0, 2]] = 0.0;
     a[[2, 0]] = 0.0;
     let tridiag = a.extract_tridiagonal().unwrap();
-    let x: Array1<f64> = random(3);
+    let x: Array1<f64> = random_using(3, &mut rng);
     let b1 = a.dot(&x);
     let b2 = b1.clone();
     let y1 = tridiag.solve_tridiagonal_into(b1).unwrap();
@@ -180,7 +184,8 @@ fn det_tridiagonal_f64() {
 
 #[test]
 fn det_tridiagonal_random() {
-    let mut a: Array2<f64> = random((3, 3));
+    let mut rng = rand_pcg::Mcg128Xsl64::new(0xcafef00dd15ea5e5);
+    let mut a: Array2<f64> = random_using((3, 3), &mut rng);
     a[[0, 2]] = 0.0;
     a[[2, 0]] = 0.0;
     assert_aclose!(a.det_tridiagonal().unwrap(), a.det().unwrap(), 1e-7);

From 36a0276957d28fba2b6f13c2ff9c65bf19c82ff8 Mon Sep 17 00:00:00 2001
From: Toshiki Teramura <toshiki.teramura@gmail.com>
Date: Sat, 17 Sep 2022 15:33:52 +0900
Subject: [PATCH 3/4] Update doctests to use random_*_using

---
 ndarray-linalg/src/solve.rs  | 20 +++++---------------
 ndarray-linalg/src/solveh.rs | 19 +++++--------------
 2 files changed, 10 insertions(+), 29 deletions(-)

diff --git a/ndarray-linalg/src/solve.rs b/ndarray-linalg/src/solve.rs
index 2e4b568b..79df1f77 100644
--- a/ndarray-linalg/src/solve.rs
+++ b/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::*;
diff --git a/ndarray-linalg/src/solveh.rs b/ndarray-linalg/src/solveh.rs
index c6dccba4..df0b3a6f 100644
--- a/ndarray-linalg/src/solveh.rs
+++ b/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,8 +19,6 @@
 //! 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
@@ -33,20 +26,18 @@
 //! 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::*;

From 4ea4f369d5485e5c7af44046eb3ab5da8c1413c6 Mon Sep 17 00:00:00 2001
From: Toshiki Teramura <toshiki.teramura@gmail.com>
Date: Sat, 17 Sep 2022 15:48:24 +0900
Subject: [PATCH 4/4] Fix RNG for testing LOBPCG

---
 ndarray-linalg/src/lobpcg/lobpcg.rs | 18 ++++++++++++------
 ndarray-linalg/src/lobpcg/svd.rs    |  3 ++-
 2 files changed, 14 insertions(+), 7 deletions(-)

diff --git a/ndarray-linalg/src/lobpcg/lobpcg.rs b/ndarray-linalg/src/lobpcg/lobpcg.rs
index 8b7ab616..10d12da2 100644
--- a/ndarray-linalg/src/lobpcg/lobpcg.rs
+++ b/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.],
diff --git a/ndarray-linalg/src/lobpcg/svd.rs b/ndarray-linalg/src/lobpcg/svd.rs
index f46857cc..62d18b49 100644
--- a/ndarray-linalg/src/lobpcg/svd.rs
+++ b/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)