Add Transform<FPP, GPU2>::adjoint and rewrite power_iteration and spectralNorm...

Add Transform<FPP, GPU2>::adjoint and rewrite power_iteration and spectralNorm to be independent from gpu_mod.

src/faust_linear_operator/CPU/faust_linear_algebra.hpp

@@ -814,7 +814,7 @@ FPP Faust::power_iteration(const  Faust::LinearOperator<FPP,Cpu> & A, const faus
 	    #ifdef __COMPILE_TIMERS__
             A.t_power_iteration.stop();
         #endif
-        handleError("linear_algebra "," power_iteration : Faust::Transform<FPP,Cpu> 1 must be a squared matrix");
+        handleError("linear_algebra "," power_iteration : Faust::Transform<FPP,Cpu> must be a squared matrix");
 	}
 	Faust::Vect<FPP,Cpu> xk(nb_col);
 	if(rand_init)

src/faust_linear_operator/GPU2/faust_Transform_gpu.cpp.in

@@ -155,27 +155,6 @@ namespace Faust
 			marr_funcs->scalar_mul(gpu_mat_arr, reinterpret_cast<const @GM_SCALAR@*>(&alpha));
 		}
 
-
-	template<>
-		Real<@FAUST_SCALAR_FOR_GM@> Transform<@FAUST_SCALAR_FOR_GM@,GPU2>::spectralNorm(int32_t nb_iter_max, float threshold, int& flag)
-		{
-			auto marr_funcs = GPUModHandler::get_singleton()->marr_funcs((@FAUST_SCALAR_FOR_GM@)(0));
-			if(gpu_mat_arr == nullptr) throw std::runtime_error("gpu_mat_arr is nullptr");
-			//TODO: update gpu_mod to handle output flag
-			return marr_funcs->spectral_norm(gpu_mat_arr, threshold, nb_iter_max);
-		}
-
-	template<>
-		@FAUST_SCALAR_FOR_GM@ Transform<@FAUST_SCALAR_FOR_GM@,GPU2>::power_iteration(int32_t nb_iter_max, float threshold, int& flag)
-		{
-			auto marr_funcs = GPUModHandler::get_singleton()->marr_funcs((@FAUST_SCALAR_FOR_GM@)(0));
-			if(gpu_mat_arr == nullptr) throw std::runtime_error("gpu_mat_arr is nullptr");
-			//TODO: update gpu_mod to handle output flag
-			@FAUST_SCALAR_FOR_GM@ lambda;
-			marr_funcs->power_iteration(gpu_mat_arr, threshold, nb_iter_max, reinterpret_cast<@GM_SCALAR@*>(&lambda));
-			return lambda;
-		}
-
 	template<>
 		MatDense<@FAUST_SCALAR_FOR_GM@, GPU2> Transform<@FAUST_SCALAR_FOR_GM@,GPU2>::sliceMultiply(const Slice s[2], MatDense<@FAUST_SCALAR_FOR_GM@, GPU2>& gpu_X, const char opThis) const
 		{

src/faust_linear_operator/GPU2/faust_Transform_gpu.h

@@ -55,6 +55,7 @@ namespace Faust
 			bool is_fact_dense(int id) const;
 			bool is_fact_bsr(int id) const;
 			void transpose();
+			void adjoint();
 			faust_unsigned_int getNbRow()const;
 			faust_unsigned_int getNbCol()const;
 			void Display(bool transpose=false) const;
@@ -73,8 +74,10 @@ namespace Faust
 			void multiplyLeft(const Transform<FPP,GPU2> & A);
 			void multiply(const FPP& a);
 			Vect<FPP,GPU2> multiply(const Vect<FPP,GPU2>& x, const char opThis='N');
-			Real<FPP> spectralNorm(int32_t nb_iter_max, float threshold, int& flag);
-			FPP power_iteration(int32_t nb_iter_max, float threshold, int& flag);
+			Real<FPP> spectralNorm(int32_t nb_iter_max, float threshold, int& flag) const;
+			template<typename FPP2 = double>
+				FPP power_iteration(const faust_unsigned_int nbr_iter_max, const FPP2 threshold, int & flag, const bool rand_init=true) const;
+
 			Real<FPP> normL1(const bool transpose = false, const bool full_array=true, const int batch_size=1) const;
 			void tocpu(Transform<FPP, Cpu>& cpu_transf) const;
 			Transform<FPP, Cpu> tocpu() const;
@@ -110,6 +113,7 @@ namespace Faust
 			friend TransformHelper<FPP,GPU2>;
 			friend TransformHelperGen<FPP,GPU2>;
 		};
+
 }
 
 #include "faust_Transform_gpu.hpp"

src/faust_linear_operator/GPU2/faust_Transform_gpu.hpp

@@ -324,7 +324,7 @@ namespace Faust
 	template<typename FPP>
 		Vect<FPP,GPU2> Transform<FPP,GPU2>::multiply(const Vect<FPP,GPU2>& x, const char opThis/*='N'*/)
 		{
-			MatDense<FPP, GPU2> out = this->multiply((const MatDense<FPP,GPU2>)x, opThis);
+			MatDense<FPP, GPU2> out = this->multiply(*dynamic_cast<const MatDense<FPP, GPU2>*>(&x), opThis);
 			Vect<FPP,GPU2> v_out;
 			v_out.gpu_mat = out.gpu_mat;
 			out.gpu_mat = nullptr; // avoid freeing v_out.gpu_mat when out of scope
@@ -513,4 +513,95 @@ namespace Faust
 				data[i]->transpose();
 		}
 
+
+	template<typename FPP>
+		void Faust::Transform<FPP,GPU2>::adjoint()
+		{
+			transpose();
+			for(auto it = data.begin(); it != data.end(); it++){
+				(*it)->conjugate();
+			}
+		}
+
+	// compute the largest eigenvalue of A, A must be positive semi-definite
+	template<typename FPP>
+	template<typename FPP2>
+		FPP Transform<FPP, GPU2>::power_iteration(const faust_unsigned_int nbr_iter_max, const FPP2 threshold, int & flag, const bool rand_init/*=true*/) const
+		{
+			auto A = *this;
+
+			const int nb_col = A.getNbCol();
+			int i = 0;
+			flag = 0;
+
+			if (nbr_iter_max <= 0)
+			{
+				handleError("linear_algebra "," power_iteration :  nbr_iter_max <= 0");
+			}
+			if (nb_col != A.getNbRow())
+			{
+				handleError("linear_algebra "," power_iteration : Faust::Transform<FPP,GPU2> must be a square matrix");
+			}
+			Faust::Vect<FPP,GPU2> xk(nb_col);
+			if(rand_init)
+			{
+				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,GPU2> xk_norm(nb_col);
+			FPP lambda_old = 1.0;
+			FPP lambda = 0.0;
+			FPP alpha = 1.0;
+			FPP beta = 0.0;
+			while((Faust::fabs(lambda_old-lambda)> Faust::fabs(threshold) || Faust::fabs(lambda) <= Faust::fabs(threshold)) && i<nbr_iter_max)
+			{
+				i++;
+				lambda_old = lambda;
+				xk_norm = xk;
+				xk_norm.normalize();
+				xk = A.multiply(xk_norm);
+//				if(xk.isZero()) // A is most likely zero, lambda is zero //TODO
+//				{
+//					std::cerr << "WARNING: power_iteration product Ax leads to zero vector, A is most likely zero, lambda should be zero too." << std::endl;
+//					return FPP(0);
+//				}
+				lambda = xk_norm.dot(xk);
+				//std::cout << "i = " << i << " ; lambda=" << lambda << std::endl;
+			}
+			flag = (i<nbr_iter_max)?i:-1;
+
+			return lambda;
+		}
+
+
+	template<typename FPP>
+		Real<FPP> Transform<FPP,GPU2>::spectralNorm(const int nbr_iter_max, float threshold, int &flag) const // TODO: factorize with CPU code
+		{
+			if (size() == 0)
+			{
+				return 1; // TODO: why?
+			}else
+			{
+				// if(this->is_zero) // The Faust is zero by at least one of its factors
+					//return 0; //TODO
+				// The Faust can still be zero (without any of its factor being)
+				// this case will be detected in power_iteration
+				Transform<FPP,GPU2> AtA((*this));
+				AtA.adjoint();
+				if (getNbCol() < getNbRow())
+				{
+					AtA.multiply(*this);
+				}else
+				{
+					AtA.multiplyLeft(*this);
+				}
+				FPP maxAbsValue = std::sqrt(AtA.power_iteration(nbr_iter_max, threshold, flag));
+				return absValue(maxAbsValue);
+
+			}
+		}
+
+
 }