Update Faust::power_iteration to workaround Faust::Transform<float> NaN...

Update Faust::power_iteration to workaround Faust::Transform<float> NaN eigenvalue (due to zero eigenvector approximate) by recomputing the power_iteration in double precision.

Issue #236.
The PALM4MSA 2020 patch of float to double spectralNorm recomputation is kept for the moment even if it is made useless by this one (particularly because it is not yet handled by the GPU2 power_iteration impl.).

src/faust_linear_operator/CPU/faust_MatDense.h

@@ -391,6 +391,11 @@ namespace Faust
 			//	If the matrix is real, it does nothing.
 			void real();
 			void real(MatDense<Real<FPP>, Cpu> & real_mat) const;
+			/**
+			 * Contrary to real() functions it returns a new MatDense of type Real<FPP2> but left this one untouched.
+			 */
+			template<typename FPP2>
+				MatDense<FPP2, Cpu> to_real() const;
 			Real<FPP> normL1(const bool transpose=false) const;
 			Real<FPP> normL1(faust_unsigned_int&, const bool transpose=false) const;
 			Real<FPP> normInf(const bool transpose=false) const;

src/faust_linear_operator/CPU/faust_MatDense.hpp

@@ -438,6 +438,16 @@ namespace Faust
 			real_mat.mat = mat.real().eval().template cast<Real<FPP>>();
 		}
 
+	template<typename FPP>
+		template<typename FPP2>
+		MatDense<FPP2, Cpu> MatDense<FPP, Cpu>::to_real() const
+		{
+			Faust::MatDense<FPP2, Cpu> ddmat;
+			ddmat.resize(this->getNbRow(), this->getNbCol());
+			ddmat.mat = mat.real().eval().template cast<Real<FPP2>>();
+			return ddmat;
+		}
+
 	template<typename FPP>
 		void MatDense<FPP,Cpu>::multiply(MatSparse<FPP,Cpu> & M, char opThis) const
 		{
@@ -1131,7 +1141,7 @@ void MatDense<FPP,Cpu>::normalize(int norm_type/*=-2*/)
 		default:
 			throw std::runtime_error("Unknown kind of norm asked for normalization.");
 	}
-	if(n != FPP(0))
+	if(n != Real<FPP>(0))
 		scalarMultiply(FPP(1.0/n));
 	else
 		throw std::domain_error("the norm is zero, can't normalize.");

src/faust_linear_operator/CPU/faust_Vect.h

@@ -147,7 +147,7 @@ namespace Faust
 	void scalarMultiply(FPP1 const scalar){vec *= scalar;}
 		void conjugate();
 		FPP normL1() const;
-		void normalize(){scalarMultiply(1/norm());}
+		void normalize();
 
 
         // multiply (*this) =  A * (*this)

src/faust_linear_operator/CPU/faust_Vect.hpp

@@ -378,4 +378,13 @@ FPP Faust::Vect<FPP, Cpu>::max_coeff(int *index) const
 	return max;
 }
 
+template<typename FPP>
+void Faust::Vect<FPP, Cpu>::normalize()
+{
+	auto n = norm();
+	if(n == Real<FPP>(0))
+		throw std::domain_error("Can't normalize a zero-norm vector.");
+	scalarMultiply(1/norm());
+}
+
 #endif

src/faust_linear_operator/CPU/faust_linear_algebra.hpp

@@ -832,7 +832,35 @@ FPP Faust::power_iteration(const  Faust::LinearOperator<FPP,Cpu> & A, const faus
 			i++;
       		lambda_old = lambda;
       		xk_norm = xk;
-      		xk_norm.normalize();
+			try
+			{
+				xk_norm.normalize();
+			}
+			catch(std::domain_error de)
+			{
+				// this block catches division by zero error that might occur in xk_norm normalization
+				// known cases of PALM4MSA 2020 running 2-norm computation (through power_iteration) on float precision produces sometimes this error while in double precision the power iteration goes well
+				// so here we process these known cases that occur on Transform<float> 2-norm computations
+				// other cases where A is a MatDense might occur but are not handled until it makes sense (see the 2nd message below)
+				const Faust::Transform<FPP, Cpu>* t = nullptr;
+				if(std::is_same<FPP, float>::value)
+				{
+					if(t = dynamic_cast<const Faust::Transform<FPP, Cpu>*>(&A))
+					{
+						// TODO: ideally each factor should be converted in double precision
+						// here we compute the product and then convert it to float precision, it's simpler
+						auto dmat = t->get_product();
+						auto ddmat = dmat.template to_real<double>();
+						std::cerr << "WARNING: power_iteration detected a zero norm vector, trying to recompute the Transform power_iteration using double precision instead of float precision to contravene this error." << std::endl;
+						auto res = FPP(power_iteration(ddmat, nbr_iter_max, threshold, flag, (double*) nullptr, rand_init));
+						std::cerr << "The error was fixed with double precision!" << std::endl;
+						return res;
+					}
+					else
+						std::cerr << "power_iteration error: vector normalization failed (because of zero norm) in float precision context but A is a: " << typeid(A).name() << ". The workaround of recomputation in double precision is only enabled for Transform<float, Cpu> so far (this error might suggest a new need in the workaround, please contact the developers)." << std::endl;
+				}
+				throw de;
+			}
       		xk = A.multiply(xk_norm);
       		lambda = xk_norm.dot(xk);
      		//std::cout << "i = " << i << " ; lambda=" << lambda << std::endl;