Refactor bit_reversal_factor functions and add a permutation parameter to...

Refactor bit_reversal_factor functions and add a permutation parameter to butterfly_support function.

src/algorithm/factorization/faust_butterfly.h

@@ -34,11 +34,14 @@ namespace Faust
 		void butterfly_support(int nfactors, std::vector<Faust::MatSparse<FPP, Cpu>*> out);
 
 	template<typename FPP>
-		std::vector<Faust::MatSparse<FPP, Cpu>*> support_DFT(int nfactors);
+		std::vector<Faust::MatSparse<FPP, Cpu>*> butterfly_support(int nfactors, Faust::MatSparse<FPP, Cpu>* P=nullptr);
 
 	template<typename FPP>
 		void bit_reversal_factor(int nfactors, std::vector<Faust::MatSparse<FPP, Cpu>*>& out);
 
+	template<typename FPP>
+		void bit_reversal_factor(int nfactors, Faust::MatSparse<FPP, Cpu>*& out);
+
 	template<typename FPP>
 		Faust::TransformHelper<FPP, Cpu>* butterfly_hierarchical(const Faust::MatDense<FPP, Cpu>& A, const std::vector<Faust::MatSparse<FPP, Cpu>*> &supports, const ButterflyFactDir& dir=RIGHT);
 

src/algorithm/factorization/faust_butterfly.hpp

@@ -326,9 +326,9 @@ namespace Faust
 		}
 
 	template<typename FPP>
-		std::vector<Faust::MatSparse<FPP, Cpu>*> support_DFT(int nfactors)
+		std::vector<Faust::MatSparse<FPP, Cpu>*> butterfly_support(int nfactors, Faust::MatSparse<FPP, Cpu>* P/*=nullptr*/)
 		{
-			// WARNING: free the MatSparse of out after calling this function
+			// WARNING: free the MatSparse after calling this function
 			std::vector<Faust::MatSparse<FPP, Cpu>*> out;
 			auto size = 1 << nfactors;
 			Faust::MatDense<FPP, Cpu> kernel, id;
@@ -343,28 +343,31 @@ namespace Faust
 				sp_support = new Faust::MatSparse<FPP, Cpu>(support);
 				out.push_back(sp_support);
 			}
-			if(! std::is_same<FPP, Real<FPP>>::value)
+			if(P == nullptr)
 			{
-				bit_reversal_factor(nfactors, out);
+				if(! std::is_same<FPP, Real<FPP>>::value)
+					bit_reversal_factor(nfactors, out);
 			}
+			else
+				out.push_back(P);
 			return out;
 		}
 
 	// TODO: specializations should be in .cpp file
 	template<>
-		void bit_reversal_factor<double>(int nfactors, std::vector<Faust::MatSparse<double, Cpu>*> &out)
+		void bit_reversal_factor<double>(int nfactors, Faust::MatSparse<double, Cpu>*& out)
 		{
 			//nothing to do for double matrices
 		}
 
 	template<>
-		void bit_reversal_factor<float>(int nfactors, std::vector<Faust::MatSparse<float, Cpu>*> &out)
+		void bit_reversal_factor<float>(int nfactors, Faust::MatSparse<float, Cpu>*& out)
 		{
 			//nothing to do for float matrices
 		}
 
 	template<>
-		void bit_reversal_factor<std::complex<double>>(int nfactors, std::vector<Faust::MatSparse<std::complex<double>, Cpu>*>& out)
+		void bit_reversal_factor<std::complex<double>>(int nfactors, Faust::MatSparse<std::complex<double>, Cpu>*& out)
 		{
 			// bit reversal permutation factor
 			unsigned int dim_size = 1u << nfactors, L, r, L_over_2, L_times_2;
@@ -377,7 +380,26 @@ namespace Faust
 			std::vector<std::complex<Real<double>>> ones(dim_size);
 			for(typename std::vector<std::complex<double>>::iterator it=ones.begin(); it != ones.end(); it++)
 				*it = std::complex<double>(1.0);
-			MatSparse<std::complex<double>,Cpu> *P = new MatSparse<std::complex<double>,Cpu>(index, new_index, ones, dim_size, dim_size);
+			out = new MatSparse<std::complex<double>,Cpu>(index, new_index, ones, dim_size, dim_size);
+		}
+
+	template<>
+		void bit_reversal_factor<double>(int nfactors, std::vector<Faust::MatSparse<double, Cpu>*> &out)
+		{
+			//nothing to do for double matrices
+		}
+
+	template<>
+		void bit_reversal_factor<float>(int nfactors, std::vector<Faust::MatSparse<float, Cpu>*> &out)
+		{
+			//nothing to do for float matrices
+		}
+
+	template<>
+		void bit_reversal_factor<std::complex<double>>(int nfactors, std::vector<Faust::MatSparse<std::complex<double>, Cpu>*>& out)
+		{
+			MatSparse<std::complex<double>,Cpu> *P;
+			bit_reversal_factor(nfactors, P);
 			out.push_back(P);
 		}
 
@@ -648,7 +670,7 @@ namespace Faust
 				throw std::runtime_error("The matrix to factorize must be square.");
 			if(log2_size - int(log2_size) > std::numeric_limits<Real<FPP>>::epsilon())
 				throw std::runtime_error("The matrix to factorize must be of a size equal to a power of two.");
-			auto support = support_DFT<FPP>((int) log2_size);
+			auto support = butterfly_support<FPP>((int) log2_size);
 			//	std::cout << "support norms" << std::endl;
 			//	for(auto s: support)
 			//		std::cout << s->norm() << std::endl;