From 85f6344b4e8f36c327308dc2c2c988e0b9d9d211 Mon Sep 17 00:00:00 2001
From: Toshiki Teramura <toshiki.teramura@gmail.com>
Date: Fri, 16 Sep 2022 22:23:04 +0900
Subject: [PATCH] ndarray_linalg::generate::random_*_using API for using given
 RNG

---
 ndarray-linalg/src/generate.rs | 100 ++++++++++++++++++++++++++++++---
 1 file changed, 92 insertions(+), 8 deletions(-)

diff --git a/ndarray-linalg/src/generate.rs b/ndarray-linalg/src/generate.rs
index 67def14a..2efd15b2 100644
--- a/ndarray-linalg/src/generate.rs
+++ b/ndarray-linalg/src/generate.rs
@@ -22,7 +22,10 @@ where
     a
 }
 
-/// Generate random array
+/// Generate random array with given shape
+///
+/// - This function uses [rand::thread_rng].
+///   See [random_using] for using another RNG
 pub fn random<A, S, Sh, D>(sh: Sh) -> ArrayBase<S, D>
 where
     A: Scalar,
@@ -31,29 +34,77 @@ where
     Sh: ShapeBuilder<Dim = D>,
 {
     let mut rng = thread_rng();
-    ArrayBase::from_shape_fn(sh, |_| A::rand(&mut rng))
+    random_using(sh, &mut rng)
+}
+
+/// Generate random array with given RNG
+///
+/// - See [random] for using default RNG
+pub fn random_using<A, S, Sh, D, R>(sh: Sh, rng: &mut R) -> ArrayBase<S, D>
+where
+    A: Scalar,
+    S: DataOwned<Elem = A>,
+    D: Dimension,
+    Sh: ShapeBuilder<Dim = D>,
+    R: Rng,
+{
+    ArrayBase::from_shape_fn(sh, |_| A::rand(rng))
 }
 
 /// Generate random unitary matrix using QR decomposition
 ///
-/// Be sure that this it **NOT** a uniform distribution. Use it only for test purpose.
+/// - Be sure that this it **NOT** a uniform distribution.
+///   Use it only for test purpose.
+/// - This function uses [rand::thread_rng].
+///   See [random_unitary_using] for using another RNG.
 pub fn random_unitary<A>(n: usize) -> Array2<A>
 where
     A: Scalar + Lapack,
 {
-    let a: Array2<A> = random((n, n));
+    let mut rng = thread_rng();
+    random_unitary_using(n, &mut rng)
+}
+
+/// Generate random unitary matrix using QR decomposition with given RNG
+///
+/// - Be sure that this it **NOT** a uniform distribution.
+///   Use it only for test purpose.
+/// - See [random_unitary] for using default RNG.
+pub fn random_unitary_using<A, R>(n: usize, rng: &mut R) -> Array2<A>
+where
+    A: Scalar + Lapack,
+    R: Rng,
+{
+    let a: Array2<A> = random_using((n, n), rng);
     let (q, _r) = a.qr_into().unwrap();
     q
 }
 
 /// Generate random regular matrix
 ///
-/// Be sure that this it **NOT** a uniform distribution. Use it only for test purpose.
+/// - Be sure that this it **NOT** a uniform distribution.
+///   Use it only for test purpose.
+/// - This function uses [rand::thread_rng].
+///   See [random_regular_using] for using another RNG.
 pub fn random_regular<A>(n: usize) -> Array2<A>
 where
     A: Scalar + Lapack,
 {
-    let a: Array2<A> = random((n, n));
+    let mut rng = rand::thread_rng();
+    random_regular_using(n, &mut rng)
+}
+
+/// Generate random regular matrix with given RNG
+///
+/// - Be sure that this it **NOT** a uniform distribution.
+///   Use it only for test purpose.
+/// - See [random_regular] for using default RNG.
+pub fn random_regular_using<A, R>(n: usize, rng: &mut R) -> Array2<A>
+where
+    A: Scalar + Lapack,
+    R: Rng,
+{
+    let a: Array2<A> = random_using((n, n), rng);
     let (q, mut r) = a.qr_into().unwrap();
     for i in 0..n {
         r[(i, i)] = A::one() + A::from_real(r[(i, i)].abs());
@@ -62,12 +113,28 @@ where
 }
 
 /// Random Hermite matrix
+///
+/// - This function uses [rand::thread_rng].
+///   See [random_hermite_using] for using another RNG.
 pub fn random_hermite<A, S>(n: usize) -> ArrayBase<S, Ix2>
 where
     A: Scalar,
     S: DataOwned<Elem = A> + DataMut,
 {
-    let mut a: ArrayBase<S, Ix2> = random((n, n));
+    let mut rng = rand::thread_rng();
+    random_hermite_using(n, &mut rng)
+}
+
+/// Random Hermite matrix
+///
+/// - See [random_hermite] for using default RNG.
+pub fn random_hermite_using<A, S, R>(n: usize, rng: &mut R) -> ArrayBase<S, Ix2>
+where
+    A: Scalar,
+    S: DataOwned<Elem = A> + DataMut,
+    R: Rng,
+{
+    let mut a: ArrayBase<S, Ix2> = random_using((n, n), rng);
     for i in 0..n {
         a[(i, i)] = a[(i, i)] + a[(i, i)].conj();
         for j in (i + 1)..n {
@@ -80,13 +147,30 @@ where
 /// Random Hermite Positive-definite matrix
 ///
 /// - Eigenvalue of matrix must be larger than 1 (thus non-singular)
+/// - This function uses [rand::thread_rng].
+///   See [random_hpd_using] for using another RNG.
 ///
 pub fn random_hpd<A, S>(n: usize) -> ArrayBase<S, Ix2>
 where
     A: Scalar,
     S: DataOwned<Elem = A> + DataMut,
 {
-    let a: Array2<A> = random((n, n));
+    let mut rng = rand::thread_rng();
+    random_hpd_using(n, &mut rng)
+}
+
+/// Random Hermite Positive-definite matrix
+///
+/// - Eigenvalue of matrix must be larger than 1 (thus non-singular)
+/// - See [random_hpd] for using default RNG.
+///
+pub fn random_hpd_using<A, S, R>(n: usize, rng: &mut R) -> ArrayBase<S, Ix2>
+where
+    A: Scalar,
+    S: DataOwned<Elem = A> + DataMut,
+    R: Rng,
+{
+    let a: Array2<A> = random_using((n, n), rng);
     let ah: Array2<A> = conjugate(&a);
     ArrayBase::eye(n) + &ah.dot(&a)
 }