Update Faust::power_iteration by defaultly using random vector initialization...

Update Faust::power_iteration by defaultly using random vector initialization and always the same seed for reproducibility.

src/faust_linear_operator/CPU/faust_linear_algebra.h

@@ -161,7 +161,7 @@ namespace Faust
 	 * \param rand_init: true for random initialization (else the vector is full of ones).
 	 */
 	template<typename FPP, typename FPP2 = double>
-		FPP power_iteration(const LinearOperator<FPP,Cpu> & A, const faust_unsigned_int nbr_iter_max,FPP2 threshold,int & flag, FPP* out_vec_data=nullptr, bool rand_init=false);
+		FPP power_iteration(const LinearOperator<FPP,Cpu> & A, const faust_unsigned_int nbr_iter_max,FPP2 threshold,int & flag, FPP* out_vec_data=nullptr, bool rand_init=true);
 
 
 	//! surcharge d'operateur * pour multiplier des matrices et des vecteurs

src/faust_linear_operator/CPU/faust_linear_algebra.hpp

@@ -44,6 +44,7 @@
 #include <iostream>
 #include <stdexcept>
 #include <complex>
+#include <cstdlib>
 #include "faust_LinearOperator.h"
 #include "faust_MatGeneric.h"
 #include "faust_MatDense.h"
@@ -788,7 +789,7 @@ A.t_add_ext.stop();
 
 // compute the largest eigenvalue of A, A must be semi-definite positive
 template<typename FPP, typename FPP2>
-FPP Faust::power_iteration(const  Faust::LinearOperator<FPP,Cpu> & A, const faust_unsigned_int nbr_iter_max, FPP2 threshold, int & flag, FPP* out_vec/*=nullptr*/, const bool rand_init/*=false*/)
+FPP Faust::power_iteration(const  Faust::LinearOperator<FPP,Cpu> & A, const faust_unsigned_int nbr_iter_max, FPP2 threshold, int & flag, FPP* out_vec/*=nullptr*/, const bool rand_init/*=true*/)
 {
 	#ifdef __COMPILE_TIMERS__
 		A.t_power_iteration.start();
@@ -815,7 +816,10 @@ FPP Faust::power_iteration(const  Faust::LinearOperator<FPP,Cpu> & A, const faus
 	}
 	Faust::Vect<FPP,Cpu> xk(nb_col);
 	if(rand_init)
-		xk.setRand();
+	{
+		srand(0xF4+'u'+57); // always the same seed, not standard but allows reproducibility
+		xk.setRand(); // the main objective is to make very unlikely to be orthogonal to the main eigenvector
+	}
 	else
 		xk.setOnes();
 	Faust::Vect<FPP,Cpu> xk_norm(nb_col);