Add a first implementation of HierarchicalFactFFT (along with ParamsFFT and a test which compiles).

- New baseline test (just a copy of basis HierarchicalFact except that we use child classes HierarchicalFactFFT and ParamsFFT). A real functional test will come later but we already know that Palm4MSAFFT is working so it shouldn't be difficult.
- Change Faust::HierarchicalFact::compute_errors() access modifier from private to protected (in order to override it later in FFT child), likewise all private fields move to protected access.
- New constructor in Palm4MSAFFT to initialize from a ParamsFFT object (like Palm4MSA can initialize from a Params object).
- Other minor changes (function decl. order, typo.).

misc/test/CMakeLists.txt

@@ -192,7 +192,7 @@ if(MATIO_LIB_FILE AND MATIO_INC_DIR AND BUILD_READ_MAT_FILE) # AND HDF5_LIB_FILE
 	# faust_multiplication : time comparison between Faust-vector product and Dense matrix-vector product
  
 	foreach(TEST_FPP float double complex<float> complex<double>)	
-		foreach(testin hierarchicalFactorization test_palm4MSA test_palm4MSAFFT faust_multiplication faust_matdense_conjugate)
+		foreach(testin hierarchicalFactorization hierarchicalFactorizationFFT test_palm4MSA test_palm4MSAFFT faust_multiplication faust_matdense_conjugate)
 		string(REGEX REPLACE "<|>" "" TEST_FILE_CPP ${TEST_FPP})
 		set(TEST_CPP_NAME ${testin}_${TEST_FILE_CPP})
 		set(TEST_FILE_CPP ${TEST_CPP_NAME}.cpp)

misc/test/src/C++/hierarchicalFactorizationFFT.cpp.in

+#include "faust_MatSparse.h"
+#include "faust_HierarchicalFactFFT.h"
+#include "faust_Timer.h"
+#include "faust_Transform.h"
+#include "faust_init_from_matio_params.h"
+#include "faust_init_from_matio_core.h"
+#include "faust_ParamsFFT.h"
+#include <string>
+#include <sstream>
+#include "faust_BlasHandle.h"
+#include "faust_SpBlasHandle.h"
+
+
+#include <iostream>
+#include <iomanip>
+
+typedef @TEST_FPP@ FPP;
+typedef @TEST_FPP2@ FPP2;
+
+int main(int argc, char* argv[])
+{
+	if (typeid(FPP) == typeid(double))
+  	{
+		cout<<"floating point precision == double"<<endl;
+  	}
+
+  	if (typeid(FPP) == typeid(float))
+  	{
+		cout<<"floating point precision == float"<<endl;
+  	}
+
+	//TODO: this test is for now totally irelevant, we must get a Laplacian and a Fourier graph matrix to test properly HierarchicalFactFFT	
+	//default value
+	string configFilename = "@FAUST_CONFIG_MAT_DIR@/config_hierarchical_fact.mat";
+	string MatrixFilename = "@FAUST_DATA_MAT_DIR@/matrix_hierarchical_fact.mat";
+
+	if (argc >= 3)
+	{	
+		MatrixFilename = argv[1];
+		configFilename = argv[2];
+		
+	}
+
+	FPP expectedLambda = 0;
+	if (argc >= 4)
+		expectedLambda = atof(argv[3]);
+
+
+	FPP2 epsilon = 0.0001;
+	if (argc >= 5)
+		epsilon = atof(argv[4]);
+
+	char transposedMatrix='N';
+	if (argc >= 6)
+		transposedMatrix=argv[5][0];
+
+
+	if ((transposedMatrix != 'N') && (transposedMatrix != 'T'))
+	{
+		cerr << "transposedMatrix value is "<<transposedMatrix<< endl;		
+		cerr << "transposedMatrix parameter must be equal to ''N'' or ''T'' " << endl;
+		exit(EXIT_FAILURE);
+	}
+
+
+	size_t ind = configFilename.find_last_of(".");
+
+	if(ind<=0 || ind>= configFilename.size())
+	{
+		cerr << "Le nom du fichier est incorrect" << endl;
+		exit(EXIT_FAILURE);
+	}
+	string configFileExtension(configFilename, ind);
+	if(configFileExtension.compare(".mat") != 0)
+	{
+		cerr << "Le nom du fichier doit se terminer par \".mat\"" << endl;
+		exit(EXIT_FAILURE);
+	}
+	string configFileBodyTmp(configFilename, 0, ind);
+
+
+	string configFileBodyDir, configFileBodyFile;
+
+	ind = configFileBodyTmp.find_last_of("/");
+	if(ind<=0 || ind>= configFileBodyTmp.size())
+	{
+		configFileBodyDir = string("");
+		configFileBodyFile = configFileBodyTmp;
+	}
+	else
+	{
+		configFileBodyDir = string(configFileBodyTmp, 0, ind+1);
+		configFileBodyFile = string(configFileBodyTmp, ind+1);
+	}
+
+
+
+	// useless for CPU but use for compatibility with GPU
+	Faust::BlasHandle<Cpu> blasHandle;
+	Faust::SpBlasHandle<Cpu> spblasHandle;
+
+	// parameter setting
+	Faust::ParamsFFT<FPP,Cpu,FPP2> params;
+	init_params_from_matiofile<FPP,Cpu,FPP2>(params,configFilename.c_str(),"params");
+
+	params.Display();
+
+	// matrix to be factorized
+	Faust::MatDense<FPP,Cpu> matrix;
+	init_faust_mat_from_matio(matrix,MatrixFilename.c_str(),"matrix");
+	
+	// transposed the matrix if needed
+	if (transposedMatrix == 'T')
+		matrix.transpose();
+	
+	//algorithm
+	Faust::HierarchicalFactFFT<FPP,Cpu,FPP2> hierFact(matrix,matrix,params,blasHandle,spblasHandle);
+	Faust::Timer t1;
+	t1.start();
+
+	hierFact.compute_facts();
+
+	t1.stop();
+	#ifdef __COMPILE_TIMERS__
+		hierFact.print_timers();
+		//hierFact.print_prox_timers();
+	#endif
+	cout <<"total hierarchical fact = "<<t1.get_time()<<endl;
+
+	vector<Faust::MatSparse<FPP,Cpu> > facts;
+	hierFact.get_facts(facts);
+	FPP lambda = hierFact.get_lambda();
+	if (argc >= 3)
+	{
+		if (Faust::fabs(lambda - expectedLambda) > epsilon)
+		{
+			std::cerr<<"invalid lambda, must be equal to "<<std::setprecision(20)<<std::setprecision(20)<<expectedLambda<<" in the precision of "<<epsilon<<std::endl;
+			std::cerr<<"current value is "<<std::setprecision(20)<<lambda<<std::endl;
+			exit(EXIT_FAILURE);
+		}
+
+	}
+	(facts[0]) *= hierFact.get_lambda();
+	// transform the sparse matrix into generic one
+	std::vector<Faust::MatGeneric<FPP,Cpu> *> list_fact_generic;
+	list_fact_generic.resize(facts.size());
+	for (int i=0;i<list_fact_generic.size();i++)
+		list_fact_generic[i]=facts[i].Clone();
+
+	Faust::Transform<FPP,Cpu> hierFactCore(list_fact_generic);
+
+	for (int i=0;i<list_fact_generic.size();i++)
+		delete list_fact_generic[i];
+
+	char nomFichier[100];
+	string outputFile="@FAUST_BIN_TEST_OUTPUT_DIR@/hier_fact_factorisation.dat";
+	
+	//WARNING no implemented	
+	// hierFactCore.print_file(outputFile.c_str());
+
+	//write the given factorisation into a mat file
+	stringstream outputFilename;
+	outputFilename<<"@FAUST_BIN_TEST_OUTPUT_DIR@/"<<configFileBodyFile<<"_factorisation.mat";
+	std::cout<<"**************** WRITING FACTORISATION INTO ****************"<<std::endl;
+	std::cout<<"output filename : "<<outputFilename.str();
+	//WARNING no implemented	
+	// hierFactCore.print_file(outputFile.c_str());	
+	// modif NB : v1102 not implemented
+	//write_faust_core_into_matfile(hierFactCore,outputFilename.str().c_str(),"fact");
+
+
+	//relativeError
+	Faust::MatDense<FPP,Cpu> faustProduct;
+	faustProduct=hierFactCore.get_product();
+	faustProduct-=matrix;
+	FPP2 relativeError = Faust::fabs(faustProduct.norm()/matrix.norm());
+
+	std::cout<<std::endl;
+	std::cout<<"**************** RELATIVE ERROR BETWEEN FAUST AND DATA MATRIX **************** "<<std::endl;
+	std::cout<<"		"<<relativeError<<std::endl<<std::endl;
+
+
+	//time comparison between matrix vector product and faust-vector product
+	int niterTimeComp = 10;
+	if (niterTimeComp > 0)
+	{
+
+		Faust::Timer tdense;
+		Faust::Timer tfaust;
+		Faust::Vect<FPP,Cpu> x(matrix.getNbCol());
+		Faust::Vect<FPP,Cpu> ydense(matrix.getNbRow());
+		Faust::Vect<FPP,Cpu> yfaust(hierFactCore.getNbRow());
+		for (int i=0;i<niterTimeComp;i++)
+		{
+			//random initilisation of vector x
+	 		for (int j=0;j<x.size();j++)
+			{
+				x[j]=std::rand()*2.0/RAND_MAX-1.0;
+			}
+
+			tdense.start();
+			ydense = matrix * x;
+			tdense.stop();
+
+	 		tfaust.start();
+			yfaust = hierFactCore * x;
+			tfaust.stop();
+
+	 	}
+		std::cout<<std::endl;
+
+		std::cout<<"**************** TIME COMPARISON MATRIX VECTOR PRODUCT **************** "<<std::endl;
+		std::cout<<"	TIME  SPEED-UP : "<<tdense.get_time()/tfaust.get_time()<<std::endl;
+		std::cout<<"	MEAN TIME dense : "<<tdense.get_time()/((float) niterTimeComp)<<std::endl;
+		std::cout<<"	MEAN TIME faust : "<<tfaust.get_time()/((float) niterTimeComp)<<std::endl;
+		cout<<"lambda="<<std::setprecision(20)<<hierFact.get_lambda()<<endl;
+	}
+
+	blasHandle.Destroy();
+	spblasHandle.Destroy();
+
+
+	return 0;
+
+}

src/algorithm/factorization/faust_HierarchicalFact.h

@@ -76,7 +76,7 @@ namespace Faust
 
        public:
 
-          HierarchicalFact(const Faust::MatDense<FPP,DEVICE>& M,const Faust::Params<FPP,DEVICE,FPP2>& params_, Faust::BlasHandle<DEVICE> cublasHandle, SpBlasHandle<DEVICE> cusparseHandle);
+          HierarchicalFact(const Faust::MatDense<FPP,DEVICE>& M, const Faust::Params<FPP,DEVICE,FPP2>& params_, Faust::BlasHandle<DEVICE> cublasHandle, SpBlasHandle<DEVICE> cusparseHandle);
           void get_facts(Faust::Transform<FPP,DEVICE> &)const;
           void get_facts(std::vector<Faust::MatSparse<FPP,DEVICE> >&)const;
           void get_facts(std::vector<Faust::MatDense<FPP,DEVICE> >& fact)const{fact = palm_global.get_facts();}
@@ -88,11 +88,12 @@ namespace Faust
         private:
           void init();
           void next_step();
+		protected:
           void compute_errors();
 
 
 
-        private:
+        protected:
           const std::vector< std::vector<const Faust::ConstraintGeneric*>> cons;
           bool m_isUpdateWayR2L;
           bool m_isFactSideLeft;

src/algorithm/factorization/faust_HierarchicalFact.hpp

@@ -81,7 +81,7 @@ Faust::HierarchicalFact<FPP,DEVICE,FPP2>::HierarchicalFact(const Faust::MatDense
 {
 	   // check if the params and M are compatible
    if ((M.getNbRow() != params_.m_nbRow) |  (M.getNbCol() != params_.m_nbCol))
-		handleError(m_className,"constructor : params and matrix haven't compatible size"); 
+		handleError(m_className,"constructor : params and matrix haven't compatible size");
 }
 
 
@@ -251,9 +251,9 @@ void Faust::HierarchicalFact<FPP,DEVICE,FPP2>::get_facts(std::vector<Faust::MatS
    const std::vector<Faust::MatDense<FPP,DEVICE> >& full_facts = palm_global.get_facts();
    sparse_facts.resize(full_facts.size());
    for (int i=0 ; i<sparse_facts.size() ; i++)
-{
-      sparse_facts[i].init(full_facts[i],cusparse_handle);
-}
+   {
+	   sparse_facts[i].init(full_facts[i],cusparse_handle);
+   }
 }
 
 
@@ -269,7 +269,7 @@ void Faust::HierarchicalFact<FPP,DEVICE,FPP2>::compute_facts()
   init();
   for (int i=0 ; i<=nbFact-1 ; i++)
   {
-     cout << "Faust::HierarchicalFact<FPP,DEVICE,FPP2>::compute_facts : factorisation "<<i+1<<"/"<<nbFact <<endl;
+     cout << "Faust::HierarchicalFact<FPP,DEVICE,FPP2>::compute_facts : factorization "<<i+1<<"/"<<nbFact <<endl;
      next_step();
   }
 

src/algorithm/factorization/faust_HierarchicalFactFFT.h

+#ifndef __FAUST_HIERARCHICAL_FACT_FFT_H__
+#define __FAUST_HIERARCHICAL_FACT_FFT_H__
+#include "faust_Palm4MSAFFT.h"
+#include "faust_ParamsFFT.h"
+#include "faust_HierarchicalFact.h"
+
+using namespace Faust;
+
+namespace Faust
+{
+    template<typename FPP,Device DEVICE,typename FPP2 = double>
+    class HierarchicalFactFFT : public HierarchicalFact<FPP, DEVICE, FPP2>
+    {
+
+		Palm4MSAFFT<FPP,DEVICE,FPP2> palm_global;
+		static const char * m_className;
+
+		public:
+		//TODO: move def. code in .hpp
+		HierarchicalFactFFT(const MatDense<FPP,DEVICE>& U, const MatDense<FPP,DEVICE>& Lap,const ParamsFFT<FPP,DEVICE,FPP2>& params, BlasHandle<DEVICE> cublasHandle, SpBlasHandle<DEVICE> cusparseHandle): HierarchicalFact<FPP, DEVICE, FPP2>(U, params, cublasHandle, cusparseHandle), palm_global(Palm4MSAFFT<FPP,DEVICE,FPP2>(Lap, params, cublasHandle, true))
+																																			//TODO: verify if palm_global is really initialized after parent ctor call
+		{
+			if ((U.getNbRow() != params.m_nbRow) |  (U.getNbCol() != params.m_nbCol))
+				handleError(m_className,"constructor : params and Fourier matrix U haven't compatible size");
+			if((Lap.getNbRow() != params.m_nbRow) |  (Lap.getNbCol() != params.m_nbCol))
+				handleError(m_className,"constructor : params and Laplacian matrix Lap haven't compatible size");
+		}
+
+
+	};
+}
+
+#include "faust_HierarchicalFactFFT.hpp"
+
+#endif

src/algorithm/factorization/faust_HierarchicalFactFFT.hpp

+template<typename FPP,Device DEVICE,typename FPP2>
+const char * Faust::HierarchicalFactFFT<FPP,DEVICE,FPP2>::m_className="Faust::HierarchicalFactFFT";
+

src/algorithm/factorization/faust_Palm4MSA.hpp

@@ -743,7 +743,7 @@ t_local_init_fact_from_palm.start();
 
     if(!isConstraintSet)
     {
-        handleError(m_className,"init_fact_from_palm : constrainst must be set before calling init_fact_from_palm");
+        handleError(m_className,"init_fact_from_palm : constraints must be set before calling init_fact_from_palm");
     }
 
     if(isFactSideLeft)

src/algorithm/factorization/faust_Palm4MSAFFT.h

 #include "faust_ParamsPalmFFT.h"
 #include "faust_Palm4MSA.h"
+#include "faust_ParamsFFT.h"
 
 #ifndef __FAUST_PALM4MSA_FFT_H__
 #define __FAUST_PALM4MSA_FFT_H__
@@ -12,10 +13,11 @@ namespace Faust {
 		class Palm4MSAFFT : public Palm4MSA<FPP, DEVICE, FPP2>
 	{
 		MatDense<FPP, DEVICE> D; //TODO: later it will need to be Sparse (which needs to add a prototype overload for multiplication in faust_linear_algebra.h)
-		Faust::MatDense<FPP,DEVICE> D_grad_over_c; //TODO: move to sparse mat later
+		MatDense<FPP,DEVICE> D_grad_over_c; //TODO: move to sparse mat later
 		public:
 			//TODO: another ctor (like in Palm4MSA) for hierarchical algo. use
 			Palm4MSAFFT(const ParamsPalmFFT<FPP, DEVICE, FPP2>& params, const BlasHandle<DEVICE> blasHandle, const bool isGlobal=false);
+			Palm4MSAFFT(const MatDense<FPP,DEVICE>& Lap, const ParamsFFT<FPP,DEVICE,FPP2> & params, const BlasHandle<DEVICE> blasHandle, const bool isGlobal);
 			virtual void next_step();
 			const MatDense<FPP, DEVICE>& get_D();
 		private:

src/algorithm/factorization/faust_Palm4MSAFFT.hpp

@@ -6,6 +6,11 @@ Palm4MSAFFT<FPP,DEVICE,FPP2>::Palm4MSAFFT(const ParamsPalmFFT<FPP, DEVICE, FPP2>
 	//TODO: is there something to check additionally to what parent's ctor checks ?
 }
 
+template<typename FPP,Device DEVICE,typename FPP2>
+Palm4MSAFFT<FPP,DEVICE,FPP2>::Palm4MSAFFT(const MatDense<FPP,DEVICE>& Lap, const ParamsFFT<FPP,DEVICE,FPP2> & params, const BlasHandle<DEVICE> blasHandle, const bool isGlobal) : Palm4MSA<FPP,DEVICE,FPP2>(Lap, params, blasHandle, isGlobal), D(params.init_D)
+{
+
+}
 
 template <typename FPP, Device DEVICE, typename FPP2>
 void Palm4MSAFFT<FPP,DEVICE,FPP2>::compute_grad_over_c()

src/algorithm/factorization/faust_Params.h

@@ -136,8 +136,8 @@ namespace Faust
             const bool constant_step_size_ = defaultConstantStepSize,
             const FPP step_size_ = defaultStepSize);
 
-        void init_from_file(const char* filename);
         Params();
+        void init_from_file(const char* filename);
 
         void check_constraint_validity();
         void check_bool_validity();

src/algorithm/factorization/faust_ParamsFFT.h

+
+#ifndef __FAUST_PARAMSFFT_H__
+#define __FAUST_PARAMSFFT_H__
+
+#include "faust_Params.h"
+
+using namespace Faust;
+
+namespace Faust
+{
+	template<typename FPP, Device DEVICE, typename FPP2 = double>
+	class ParamsFFT : public Params<FPP, DEVICE, FPP2>
+	{
+
+		public:
+		MatDense<FPP, DEVICE> init_D; //TODO: convert to Sparse or Diag repres.
+		//TODO: does it really need to be public 
+		//TODO: move the ctor def into .hpp
+		ParamsFFT(
+				const faust_unsigned_int nbRow,
+				const faust_unsigned_int nbCol,
+				const unsigned int nbFact,
+				const std::vector<std::vector<const Faust::ConstraintGeneric*>> & cons,
+				const std::vector<Faust::MatDense<FPP,DEVICE> >& init_fact,
+				const MatDense<FPP, DEVICE>& init_D,
+				const Faust::StoppingCriterion<FPP2>& stop_crit_2facts = StoppingCriterion<FPP2>(Params<FPP,DEVICE,FPP2>::defaultNiter1),
+				const Faust::StoppingCriterion<FPP2>& stop_crit_global = StoppingCriterion<FPP2>(Params<FPP,DEVICE,FPP2>::defaultNiter2),
+				const bool isVerbose = Params<FPP,DEVICE,FPP2>::defaultVerbosity,
+				const bool isUpdateWayR2L = Params<FPP,DEVICE,FPP2>::defaultUpdateWayR2L,
+				const bool isFactSideLeft = Params<FPP,DEVICE,FPP2>::defaultFactSideLeft,
+				const FPP init_lambda = Params<FPP,DEVICE,FPP2>::defaultLambda,
+				const bool constant_step_size = Params<FPP,DEVICE,FPP2>::defaultConstantStepSize,
+				const FPP step_size = Params<FPP,DEVICE,FPP2>::defaultStepSize): Params<FPP, DEVICE, FPP2>(nbRow, nbCol, nbFact, cons, init_fact, stop_crit_2facts, stop_crit_global, isVerbose, isFactSideLeft, init_lambda, constant_step_size, step_size), init_D(init_D)
+		{
+
+		}
+
+		ParamsFFT() {}
+
+	};
+}
+
+#include "faust_ParamsFFT.hpp"
+
+#endif