@@ -1921,7 +1921,9 @@ sp_auxlib::run_aupd_plain
19211921 n = X.n_rows ; // The size of the matrix (should already be set outside).
19221922 blas_int nev = n_eigvals;
19231923
1924- resid.zeros (n);
1924+ // resid.zeros(n);
1925+ eigs_randu_filler<T> randu_filler;
1926+ randu_filler.fill (resid, n); // use deterministic starting point
19251927
19261928 // Two contraints on NCV: (NCV > NEV) for sym problems or
19271929 // (NCV > NEV + 2) for gen problems and (NCV <= N)
@@ -1956,7 +1958,8 @@ sp_auxlib::run_aupd_plain
19561958 // Real work array of length lworkl.
19571959 workl.zeros (lworkl);
19581960
1959- info = 0 ; // Set to 0 initially to use random initial vector.
1961+ // info = 0; // resid to be filled with random values by ARPACK (non-deterministic)
1962+ info = 1 ; // resid is already filled with random values (deterministic)
19601963
19611964 // All the parameters have been set or created. Time to loop a lot.
19621965 while (ido != 99 )
@@ -2109,7 +2112,9 @@ sp_auxlib::run_aupd_shiftinvert
21092112 n = X.n_rows ; // The size of the matrix (should already be set outside).
21102113 blas_int nev = n_eigvals;
21112114
2112- resid.zeros (n);
2115+ // resid.zeros(n);
2116+ eigs_randu_filler<T> randu_filler;
2117+ randu_filler.fill (resid, n); // use deterministic starting point
21132118
21142119 // Two contraints on NCV: (NCV > NEV) for sym problems or
21152120 // (NCV > NEV + 2) for gen problems and (NCV <= N)
@@ -2147,7 +2152,8 @@ sp_auxlib::run_aupd_shiftinvert
21472152 // Real work array of length lworkl.
21482153 workl.zeros (lworkl);
21492154
2150- info = 0 ; // Set to 0 initially to use random initial vector.
2155+ // info = 0; // resid to be filled with random values by ARPACK (non-deterministic)
2156+ info = 1 ; // resid is already filled with random values (deterministic)
21512157
21522158 superlu_opts superlu_opts_default;
21532159 superlu::superlu_options_t options;
@@ -2439,6 +2445,85 @@ sp_auxlib::rudimentary_sym_check(const SpMat< std::complex<T> >& X)
24392445
24402446
24412447
2448+ //
2449+
2450+
2451+
2452+ template <typename eT>
2453+ inline
2454+ eigs_randu_filler<eT>::eigs_randu_filler()
2455+ {
2456+ arma_extra_debug_sigprint ();
2457+
2458+ typedef typename std::mt19937_64::result_type local_seed_type;
2459+
2460+ local_engine.seed (local_seed_type (123 ));
2461+
2462+ typedef typename std::uniform_real_distribution<eT>::param_type local_param_type;
2463+
2464+ local_u_distr.param (local_param_type (-1.0 , +1.0 ));
2465+ }
2466+
2467+
2468+ template <typename eT>
2469+ inline
2470+ void
2471+ eigs_randu_filler<eT>::fill(podarray<eT>& X, const uword N)
2472+ {
2473+ arma_extra_debug_sigprint ();
2474+
2475+ X.set_size (N);
2476+
2477+ eT* X_mem = X.memptr ();
2478+
2479+ for (uword i=0 ; i<N; ++i) { X_mem[i] = eT ( local_u_distr (local_engine) ); }
2480+ }
2481+
2482+
2483+ template <typename T>
2484+ inline
2485+ eigs_randu_filler< std::complex <T> >::eigs_randu_filler()
2486+ {
2487+ arma_extra_debug_sigprint ();
2488+
2489+ typedef typename std::mt19937_64::result_type local_seed_type;
2490+
2491+ local_engine.seed (local_seed_type (123 ));
2492+
2493+ typedef typename std::uniform_real_distribution<T>::param_type local_param_type;
2494+
2495+ local_u_distr.param (local_param_type (-1.0 , +1.0 ));
2496+ }
2497+
2498+
2499+ template <typename T>
2500+ inline
2501+ void
2502+ eigs_randu_filler< std::complex <T> >::fill(podarray< std::complex <T> >& X, const uword N)
2503+ {
2504+ arma_extra_debug_sigprint ();
2505+
2506+ typedef typename std::complex <T> eT;
2507+
2508+ X.set_size (N);
2509+
2510+ eT* X_mem = X.memptr ();
2511+
2512+ for (uword i=0 ; i<N; ++i)
2513+ {
2514+ eT& X_mem_i = X_mem[i];
2515+
2516+ X_mem_i.real ( T (local_u_distr (local_engine)) );
2517+ X_mem_i.imag ( T (local_u_distr (local_engine)) );
2518+ }
2519+ }
2520+
2521+
2522+
2523+ //
2524+
2525+
2526+
24422527#if defined(ARMA_USE_SUPERLU)
24432528
24442529inline
0 commit comments