diff --git a/include/armadillo_bits/newarp_GenEigsSolver_bones.hpp b/include/armadillo_bits/newarp_GenEigsSolver_bones.hpp
index 3ceace6f..eabaf063 100644
--- a/include/armadillo_bits/newarp_GenEigsSolver_bones.hpp
+++ b/include/armadillo_bits/newarp_GenEigsSolver_bones.hpp
@@ -25,18 +25,18 @@ template<typename eT, int SelectionRule, typename OpType>
 class GenEigsSolver
   {
   protected:
-
+  
   const OpType&           op;        // object to conduct matrix operation, eg. matrix-vector product
   const uword             nev;       // number of eigenvalues requested
   Col< std::complex<eT> > ritz_val;  // ritz values
-
+  
   // Sort the first nev Ritz pairs in decreasing magnitude order
   // This is used to return the final results
   virtual void sort_ritzpair();
-
-
+  
+  
   private:
-
+  
   const uword             dim_n;     // dimension of matrix A
   const uword             ncv;       // number of ritz values
   uword                   nmatop;    // number of matrix operations called
@@ -54,46 +54,50 @@ class GenEigsSolver
                                      // used to test the orthogonality of vectors,
                                      // and in convergence test, tol*approx0 is
                                      // the absolute tolerance
-
+  
+  std::mt19937_64         local_rng; // local random number generator
+  
+  inline void fill_rand(eT* dest, const uword N, const uword seed_val);
+  
   // Arnoldi factorisation starting from step-k
   inline void factorise_from(uword from_k, uword to_m, const Col<eT>& fk);
-
+  
   // Implicitly restarted Arnoldi factorisation
   inline void restart(uword k);
-
+  
   // Calculate the number of converged Ritz values
   inline uword num_converged(eT tol);
-
+  
   // Return the adjusted nev for restarting
   inline uword nev_adjusted(uword nconv);
-
+  
   // Retrieve and sort ritz values and ritz vectors
   inline void retrieve_ritzpair();
-
-
+  
+  
   public:
-
+  
   //! Constructor to create a solver object.
   inline GenEigsSolver(const OpType& op_, uword nev_, uword ncv_);
-
+  
   //! Providing the initial residual vector for the algorithm.
   inline void init(eT* init_resid);
-
+  
   //! Providing a random initial residual vector.
   inline void init();
-
+  
   //! Conducting the major computation procedure.
   inline uword compute(uword maxit = 1000, eT tol = 1e-10);
-
+  
   //! Returning the number of iterations used in the computation.
   inline int num_iterations() { return niter; }
-
+  
   //! Returning the number of matrix operations used in the computation.
   inline int num_operations() { return nmatop; }
-
+  
   //! Returning the converged eigenvalues.
   inline Col< std::complex<eT> > eigenvalues();
-
+  
   //! Returning the eigenvectors associated with the converged eigenvalues.
   inline Mat< std::complex<eT> > eigenvectors(uword nvec);
   
diff --git a/include/armadillo_bits/newarp_GenEigsSolver_meat.hpp b/include/armadillo_bits/newarp_GenEigsSolver_meat.hpp
index 80eee36e..290fa4f2 100644
--- a/include/armadillo_bits/newarp_GenEigsSolver_meat.hpp
+++ b/include/armadillo_bits/newarp_GenEigsSolver_meat.hpp
@@ -20,6 +20,24 @@ namespace newarp
 {
 
 
+template<typename eT, int SelectionRule, typename OpType>
+inline
+void
+GenEigsSolver<eT, SelectionRule, OpType>::fill_rand(eT* dest, const uword N, const uword seed_val)
+  {
+  arma_extra_debug_sigprint();
+  
+  typedef typename std::mt19937_64::result_type seed_type;
+  
+  local_rng.seed( seed_type(seed_val) );
+  
+  std::uniform_real_distribution<double> dist(-1.0, +1.0);
+  
+  for(uword i=0; i < N; ++i)  { dest[i] = eT(dist(local_rng)); }
+  }
+
+
+
 template<typename eT, int SelectionRule, typename OpType>
 inline
 void
@@ -44,12 +62,16 @@ GenEigsSolver<eT, SelectionRule, OpType>::factorise_from(uword from_k, uword to_
     // to the current V, which we call a restart
     if(beta < eps)
       {
+      // // Generate new random vector for fac_f
+      // blas_int idist = 2;
+      // blas_int iseed[4] = {1, 3, 5, 7};
+      // iseed[0] = (i + 100) % 4095;
+      // blas_int n = dim_n;
+      // lapack::larnv(&idist, &iseed[0], &n, fac_f.memptr());
+      
       // Generate new random vector for fac_f
-      blas_int idist = 2;
-      blas_int iseed[4] = {1, 3, 5, 7};
-      iseed[0] = (i + 100) % 4095;
-      blas_int n = dim_n;
-      lapack::larnv(&idist, &iseed[0], &n, fac_f.memptr());
+      fill_rand(fac_f.memptr(), dim_n, i+1);
+      
       // f <- f - V * V' * f, so that f is orthogonal to V
       Mat<eT> Vs(fac_V.memptr(), dim_n, i, false); // First i columns
       Col<eT> Vf = Vs.t() * fac_f;
@@ -362,11 +384,17 @@ GenEigsSolver<eT, SelectionRule, OpType>::init()
   {
   arma_extra_debug_sigprint();
   
+  // podarray<eT> init_resid(dim_n);
+  // blas_int idist = 2;                // Uniform(-1, 1)
+  // blas_int iseed[4] = {1, 3, 5, 7};  // Fixed random seed
+  // blas_int n = dim_n;
+  // lapack::larnv(&idist, &iseed[0], &n, init_resid.memptr());
+  // init(init_resid.memptr());
+  
   podarray<eT> init_resid(dim_n);
-  blas_int idist = 2;                // Uniform(-1, 1)
-  blas_int iseed[4] = {1, 3, 5, 7};  // Fixed random seed
-  blas_int n = dim_n;
-  lapack::larnv(&idist, &iseed[0], &n, init_resid.memptr());
+  
+  fill_rand(init_resid.memptr(), dim_n, 0);
+  
   init(init_resid.memptr());
   }