SimpleRandom.h
Go to the documentation of this file.
1 // Copyright (C) 2016-2019 Yixuan Qiu <yixuan.qiu@cos.name>
2 //
3 // This Source Code Form is subject to the terms of the Mozilla
4 // Public License v. 2.0. If a copy of the MPL was not distributed
5 // with this file, You can obtain one at https://mozilla.org/MPL/2.0/.
6 
7 #ifndef SIMPLE_RANDOM_H
8 #define SIMPLE_RANDOM_H
9 
10 #include <Eigen/Core>
11 
13 
14 namespace Spectra {
15 
16 // We need a simple pseudo random number generator here:
17 // 1. It is used to generate initial and restarted residual vector.
18 // 2. It is not necessary to be so "random" and advanced. All we hope
19 // is that the residual vector is not in the space spanned by the
20 // current Krylov space. This should be met almost surely.
21 // 3. We don't want to call RNG in C++, since we actually want the
22 // algorithm to be deterministic. Also, calling RNG in C/C++ is not
23 // allowed in R packages submitted to CRAN.
24 // 4. The method should be as simple as possible, so an LCG is enough.
25 // 5. Based on public domain code by Ray Gardner
26 // http://stjarnhimlen.se/snippets/rg_rand.c
27 
28 template <typename Scalar = double>
29 class SimpleRandom
30 {
31 private:
32  typedef Eigen::Index Index;
34 
35  const unsigned int m_a; // multiplier
36  const unsigned long m_max; // 2^31 - 1
37  long m_rand;
38 
39  inline long next_long_rand(long seed)
40  {
41  unsigned long lo, hi;
42 
43  lo = m_a * (long) (seed & 0xFFFF);
44  hi = m_a * (long) ((unsigned long) seed >> 16);
45  lo += (hi & 0x7FFF) << 16;
46  if (lo > m_max)
47  {
48  lo &= m_max;
49  ++lo;
50  }
51  lo += hi >> 15;
52  if (lo > m_max)
53  {
54  lo &= m_max;
55  ++lo;
56  }
57  return (long) lo;
58  }
59 
60 public:
61  SimpleRandom(unsigned long init_seed) :
62  m_a(16807),
63  m_max(2147483647L),
64  m_rand(init_seed ? (init_seed & m_max) : 1)
65  {}
66 
67  Scalar random()
68  {
69  m_rand = next_long_rand(m_rand);
70  return Scalar(m_rand) / Scalar(m_max) - Scalar(0.5);
71  }
72 
73  // Vector of random numbers of type Scalar
74  // Ranging from -0.5 to 0.5
75  Vector random_vec(const Index len)
76  {
77  Vector res(len);
78  for (Index i = 0; i < len; i++)
79  {
80  m_rand = next_long_rand(m_rand);
81  res[i] = Scalar(m_rand) / Scalar(m_max) - Scalar(0.5);
82  }
83  return res;
84  }
85 };
86 
87 } // namespace Spectra
88 
90 
91 #endif // SIMPLE_RANDOM_H
SCALAR Scalar
Definition: bench_gemm.cpp:46
MatrixXd L
Definition: LLT_example.cpp:6
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
size_t len(handle h)
Get the length of a Python object.
Definition: pytypes.h:2244
Eigen::Matrix< double, Eigen::Dynamic, 1 > Vector


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:35:43