Add GivensFGFTParallel and a test (not validated yet).

- Update GivensFGFT (parent class) in the same goal.
- Update the matlab test (to test the same Laplacian).

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 hierarchicalFactorizationFFT test_palm4MSA test_palm4MSAFFT faust_multiplication faust_matdense_conjugate GivensFGFT)
+		foreach(testin hierarchicalFactorization hierarchicalFactorizationFFT test_palm4MSA test_palm4MSAFFT faust_multiplication faust_matdense_conjugate GivensFGFT GivensFGFTParallel)
 		if(${TEST_FPP} MATCHES complex AND ${testin} MATCHES GivensFGFT)
 			continue()
 			# GivensFGFT doesn't handle complex matrices

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

+
+#include "faust_GivensFGFTParallel.h"
+#include "faust_MatDense.h"
+#include "faust_init_from_matio_params.h"
+#include "faust_init_from_matio_core.h"
+
+#include <complex>
+
+using namespace Faust;
+typedef @TEST_FPP@ FPP;
+typedef @TEST_FPP2@ FPP2;
+
+/***
+ * This test results have to be compared with what misc/test/src/Matlab/test_GivensDiagPar.m outputs.
+ * That's a way to validate the C++ impl. of Parallel Givens Jacobi FGFT.
+ */
+
+
+int main()
+{
+	string configFilename = "@FAUST_DATA_MAT_DIR@/test_GivensDiagParallel_Lap_U_J_choices.mat";
+
+	Faust::MatDense<FPP,Cpu> Lap;
+	init_faust_mat_from_matio(Lap,configFilename.c_str(),"Lap");
+
+	int J = 5472; //TODO: should be retrieved from .mat file (not be a hard-coded value)
+	GivensFGFTParallel<FPP,Cpu> algo(Lap,J,/*t=*/ Lap.getNbRow()/2);
+
+	algo.compute_facts();
+
+	vector<pair<int,int>> coord_choices = algo.get_coord_choices();
+
+	mat_t *matfp;
+
+	matfp = Mat_Open(configFilename.c_str(),MAT_ACC_RDONLY);
+	if ( NULL == matfp ) {
+		fprintf(stderr,"Error opening MAT file %s", configFilename.c_str());
+		return EXIT_FAILURE;
+	}
+	matvar_t* cell_choices = Mat_VarRead(matfp, "choices");
+	int* p_choices = new int[J];
+	int* q_choices = new int[J];
+	cout << "=======================================" << endl;
+    if ( NULL != cell_choices) {
+        for (int i = 0; i < J*2; i+=2 )
+		{
+				p_choices[i/2] = (int)((double*) (cell_choices->data))[i];
+				q_choices[i/2] = (int)((double*) (cell_choices->data))[i+1];
+		}
+        free(cell_choices);
+	}
+	Mat_Close(matfp);
+
+	for(int j=0;j<J;j++)
+		printf("ite=%d (1-base index) ref. p=%d q=%d, algo. p=%d q=%d\n", j, p_choices[j], q_choices[j], coord_choices[j].first+1, coord_choices[j].second+1);
+
+
+	Faust::MatDense<FPP,Cpu> fourier_diag(Lap.getNbRow(), Lap.getNbCol());
+	fourier_diag.setEyes();
+	const vector<Faust::MatSparse<FPP,Cpu>>& givens_facts = algo.get_facts();
+	for(int i=givens_facts.size()-1;i>=0;i--)
+		givens_facts[i].multiply(fourier_diag, 'N');
+
+	cout << "fourier_diag fro norm: " << fourier_diag.norm() << endl;
+
+	const Faust::MatSparse<FPP,Cpu> D = algo.get_D();
+	Faust::MatDense<FPP,Cpu> full_D = Faust::MatDense<FPP,Cpu>(D);
+	cout << "D fro norm:" << D.norm() << endl;
+
+	Faust::MatDense<FPP,Cpu> fourier_diag2 = algo.compute_fourier();
+	Faust::MatDense<FPP,Cpu> ordered_fourier_diag2 = algo.compute_fourier(true);
+
+	Faust::MatDense<FPP,Cpu> * ordered_fourier_diag = fourier_diag.get_cols(algo.get_ord_indices());
+
+	fourier_diag2 -= fourier_diag;
+	ordered_fourier_diag2 -= *ordered_fourier_diag;
+	cout << "norm(fourier_diag2-fourier_diag): " << fourier_diag2.norm() << endl;
+	cout << "norm(ordered_fourier_diag2-ordered_fourier_diag): " << ordered_fourier_diag2.norm() << endl;
+
+	//verify approx. fourier is properly computed (when ordered or not)
+	assert(fourier_diag2.norm()==0);
+	assert(ordered_fourier_diag2.norm()==0);
+
+	Faust::MatDense<FPP,Cpu> ordered_D(algo.get_D(true));
+
+	cout << "orderded D fro norm: " << ordered_D.norm() << endl;
+	cout << "ordered_D:" << endl;
+	ordered_D.Display();
+
+	cout << "ordered D eles" << endl;
+	for(int i=0;i<ordered_D.getNbRow();i++)
+	{
+		cout << ordered_D.getData()[i*ordered_D.getNbRow()+i] << " " << full_D(i,i) << endl;
+	}
+
+	cout << "ordered fourier_diag fro norm: " << ordered_fourier_diag->norm() << endl;
+	cout << "Lap fro norm: " << Lap.norm() << endl;
+
+	Faust::MatDense<FPP,Cpu> tmp = *ordered_fourier_diag;
+	tmp.multiplyRight(ordered_D);
+	ordered_fourier_diag->transpose();
+	tmp.multiplyRight(*ordered_fourier_diag);
+	ordered_fourier_diag->transpose();
+	tmp -= Lap;
+	cout << tmp.norm()/Lap.norm() << endl;
+
+	delete []p_choices;
+	delete []q_choices;
+	return 0;
+}
+

misc/test/src/Matlab/test_GivensDiagPar.m

@@ -4,8 +4,8 @@
 
 p = mfilename('fullpath');
 [filepath, ~, ~ ] = fileparts(p);
-load([filepath, '/../../../data/mat/test_GivensDiag_Lap_U_J.mat']) % Lap, U, J
-ref_choices = choices
+load([filepath, '/../../../data/mat/test_GivensDiagParallel_Lap_U_J_choices.mat']) % Lap, U, J
+ref_choices = choices;
 
 
 [facts_givens,D,err,L,choices] = diagonalization_givens_parall(Lap,J,size(Lap,1)/2);
@@ -33,4 +33,6 @@ disp("Error: norm(Uhat*Dhat*Uhat'- Lap, 'fro')/norm(Lap, 'fro')")
 err0 = norm(Uhat_givens*full(diag(sorted_spectrum))*Uhat_givens' - Lap, 'fro')/norm(Lap, 'fro')
 %err2 = norm(Uhat_givens'*Lap*Uhat_givens - diag(diag(Uhat_givens'*Lap*Uhat_givens)),'fro')/norm(Lap,'fro') % equals err0
 %
-
+disp("Verifying that reference choices for pivots are respected by this exec.")
+disp("ref_choices == choices:")
+all(all(ref_choices == choices))

src/algorithm/factorization/faust_GivensFGFT.h

@@ -13,6 +13,7 @@ namespace Faust {
 	template<typename FPP, Device DEVICE, typename FPP2 = float>
 		class GivensFGFT {
 
+			protected:
 			vector<Faust::MatSparse<FPP,DEVICE>> facts;
 			Faust::MatSparse<FPP,DEVICE> D;
 			Faust::MatDense<FPP,DEVICE> C;
@@ -35,7 +36,7 @@ namespace Faust {
 			bool is_D_ordered;
 
 			/**
-			 * Matrix pivot and column indices.
+			 * L pivot row and column indices.
 			 */
 			int p, q;
 
@@ -44,6 +45,9 @@ namespace Faust {
 			 */
 			unsigned int ite;
 
+			// in calc_theta() two values are calculated, set this bool to true to always choose theta2 (useful for GivensFGFTParallel)
+			bool always_theta2;
+
 			/**
 			 * Function pointer to any step of the algorithm.
 			 */
@@ -52,27 +56,27 @@ namespace Faust {
 
 			public:
 			GivensFGFT(Faust::MatDense<FPP,DEVICE>& Lap, int J);
-			~GivensFGFT() {delete[] q_candidates;};
+			virtual ~GivensFGFT() {delete[] q_candidates;};
 
 			/**
 			 * \brief Algo. main step.
 			 */
-			void next_step();
+			virtual void next_step();
 
 			/**
 			 * \brief Algo. main loop (facts.size() iterations).
 			 */
 			void compute_facts();
 
-			private:
+			protected:
 			/** \brief Algo. step 2.1.
 			*/
-			void choose_pivot();
+			virtual void choose_pivot();
 
 			/** \brief Algo. step 2.1.1
 			 *
 			 */
-			void max_L_into_C();
+			virtual void max_L();
 
 			/** \brief Algo. step 2.2.
 			*/
@@ -81,7 +85,7 @@ namespace Faust {
 			/**
 			 * \brief Algo. step 2.3.
 			 */
-			void update_fact();
+			virtual void update_fact();
 
 			/**
 			 * \brief Algo. step 2.4

src/algorithm/factorization/faust_GivensFGFT.hpp

@@ -7,7 +7,7 @@ void GivensFGFT<FPP,DEVICE,FPP2>::next_step()
 {
 
 	substep_fun substep[] = {
-		&GivensFGFT<FPP,DEVICE,FPP2>::max_L_into_C,
+		&GivensFGFT<FPP,DEVICE,FPP2>::max_L,
 		&GivensFGFT<FPP,DEVICE,FPP2>::choose_pivot,
 		&GivensFGFT<FPP,DEVICE,FPP2>::calc_theta,
 		&GivensFGFT<FPP,DEVICE,FPP2>::update_fact,
@@ -38,7 +38,7 @@ void GivensFGFT<FPP,DEVICE,FPP2>::choose_pivot()
 }
 
 template<typename FPP, Device DEVICE, typename FPP2>
-void GivensFGFT<FPP,DEVICE,FPP2>::max_L_into_C()
+void GivensFGFT<FPP,DEVICE,FPP2>::max_L()
 {
 	// Matlab ref. code:
 	//**************  at initialization
@@ -139,7 +139,7 @@ void GivensFGFT<FPP,DEVICE,FPP2>::calc_theta()
 	theta2 = theta1 + M_PI_4; // from cmath
 	err_theta1 = calc_err(theta1);
 	err_theta2 = calc_err(theta2);
-	if(err_theta1 < err_theta2)
+	if(err_theta1 < err_theta2 && !always_theta2)
 		theta = theta1;
 	else
 		theta = theta2;
@@ -277,7 +277,7 @@ void GivensFGFT<FPP,DEVICE,FPP2>::compute_facts()
 }
 
 template<typename FPP, Device DEVICE, typename FPP2>
-GivensFGFT<FPP,DEVICE,FPP2>::GivensFGFT(Faust::MatDense<FPP,DEVICE>& Lap, int J) : Lap(Lap), facts(J), D(Lap.getNbRow(), Lap.getNbCol()), C(Lap.getNbRow(), Lap.getNbCol()), errs(J), coord_choices(J), L(Lap), q_candidates(new int[Lap.getNbCol()]), is_D_ordered(false)
+GivensFGFT<FPP,DEVICE,FPP2>::GivensFGFT(Faust::MatDense<FPP,DEVICE>& Lap, int J) : Lap(Lap), facts(J), D(Lap.getNbRow(), Lap.getNbCol()), C(Lap.getNbRow(), Lap.getNbCol()), errs(J), coord_choices(J), L(Lap), q_candidates(new int[Lap.getNbCol()]), is_D_ordered(false), always_theta2(false)
 {
 	/** Matlab ref. code:
 	 *     facts = cell(1,J);

src/algorithm/factorization/faust_GivensFGFTParallel.h

+#include "faust_GivensFGFT.h"
+
+#include <list>
+
+#ifndef __GIVENS_FGFT_PARALLEL__
+#define __GIVENS_FGFT_PARALLEL__
+
+namespace Faust {
+
+
+	template<typename FPP, Device DEVICE, typename FPP2 = float>
+		class GivensFGFTParallel : public GivensFGFT<FPP,DEVICE,FPP2>
+	{
+
+		/** Maximum number of rotations per factor
+		* the effective number of rotations is the min(t, number_of_pivot_candidates)
+		*/
+		unsigned int t;
+		// Number of rotations already made for the current factor facts[ite]
+		unsigned int fact_nrots;
+		// current fact nonzero indices (descendantly sorted)
+		list<pair<int,int>> fact_nz_inds;
+
+		void max_L();
+		void update_fact_nz_inds();
+		void loop_update_fact();
+		void choose_pivot();
+		void update_fact();
+		void finish_fact();
+
+		/**
+		 * Function pointer to any step of the algorithm.
+		 */
+		typedef void (GivensFGFTParallel<FPP,DEVICE,FPP2>::*substep_fun)();
+
+
+		public:
+
+		void next_step();
+
+		GivensFGFTParallel(Faust::MatDense<FPP,DEVICE>& Lap, int J, int t);
+
+	};
+
+#include "faust_GivensFGFTParallel.hpp"
+
+
+}
+#endif

src/algorithm/factorization/faust_GivensFGFTParallel.hpp

+using namespace Faust;
+
+#include <cmath>
+
+template<typename FPP, Device DEVICE, typename FPP2>
+GivensFGFTParallel<FPP,DEVICE,FPP2>::GivensFGFTParallel(Faust::MatDense<FPP,DEVICE>& Lap, int J, int t) : GivensFGFT<FPP,DEVICE,FPP2>(Lap,J), t(t), fact_nrots(0)
+{
+	this->facts.resize(round(J/(float)t));
+	this->always_theta2 = true;
+}
+
+template<typename FPP, Device DEVICE, typename FPP2>
+void GivensFGFTParallel<FPP,DEVICE,FPP2>::max_L()
+{
+	// 	  matlab ref
+	// 	  L_low = abs(tril(L,-1));
+	//    ind_nnz = find(L_low);
+	//    [~,ind_sorted] = sort(L_low(ind_nnz),'descend');
+	//    [Rvec, Svec] = ind2sub([n,n],ind_nnz(ind_sorted));
+	//    RSmat = [Rvec, Svec];
+	MatDense<FPP,DEVICE> L_low = this->L;
+	L_low.abs();
+	L_low =	L_low.lower_tri(false);
+	cout << " GivensFGFTParallel::max_L() norm(L_low): " <<  L_low.norm() << endl;
+	fact_nz_inds = L_low.nonzeros_indices();
+
+#ifdef DEBUG_GIVENS
+	for(auto &p : fact_nz_inds)
+		cout << "GivensFGFTParallel::max_L() before sort :" << L_low(p.first, p.second) << endl;
+#endif
+
+	//ENOTE: can't use sort(it,it, lambda) as vector because list doesn't provide random access it.
+	fact_nz_inds.sort([&L_low](const pair<int,int> &a, const pair<int,int> &b)
+			{
+			return L_low(a.first, a.second) > L_low(b.first, b.second);
+			});
+#ifdef DEBUG_GIVENS
+	for(auto &p : fact_nz_inds)
+		cout << "GivensFGFTParallel::max_L() after sort :" << L_low(p.first, p.second) << endl;
+#endif
+}
+
+template<typename FPP, Device DEVICE, typename FPP2>
+void GivensFGFTParallel<FPP,DEVICE,FPP2>::update_fact_nz_inds()
+{
+	// matlab ref
+//	other = RSmat==p | RSmat==q;
+//	tokill = any(other,2);
+//	%RSmat( tokill, : ) = [];
+//	RSmatnew = RSmat(~tokill,:);
+//	RSmat = RSmatnew;
+
+	//remove all pairs containing p or q
+	for(auto i = fact_nz_inds.begin(); i != fact_nz_inds.end();)
+	{
+		if((*i).first == this->p || (*i).second == this->q)
+			i = fact_nz_inds.erase(i);
+		else
+			i++;
+	}
+}
+
+template<typename FPP, Device DEVICE, typename FPP2>
+void GivensFGFTParallel<FPP,DEVICE,FPP2>::next_step()
+{
+
+	substep_fun substep[] = {
+		&GivensFGFTParallel<FPP,DEVICE,FPP2>::max_L,
+		&GivensFGFTParallel<FPP,DEVICE,FPP2>::loop_update_fact, //responsible to call choose_pivot() and calc_theta(), update_fact()
+		&GivensFGFTParallel<FPP,DEVICE,FPP2>::update_L,
+		&GivensFGFTParallel<FPP,DEVICE,FPP2>::update_D,
+		&GivensFGFTParallel<FPP,DEVICE,FPP2>::update_err};
+
+	for(int i=0;i<sizeof(substep)/sizeof(substep_fun);i++)
+	{
+#ifdef DEBUG_GIVENS
+		cout << "GivensFGFTParallel ite=" << this->ite << " substep i=" << i << endl;
+#endif
+		(this->*substep[i])();
+	}
+}
+
+template<typename FPP, Device DEVICE, typename FPP2>
+void GivensFGFTParallel<FPP,DEVICE,FPP2>::loop_update_fact()
+{
+	fact_nrots = 0;
+	int k = 0;
+	while(k < t && k < fact_nz_inds.size())
+	{
+		choose_pivot();
+		update_fact_nz_inds();
+		this->calc_theta();
+		this->update_fact();
+		fact_nrots++;
+		k++;
+	}
+	finish_fact();
+}
+
+template<typename FPP, Device DEVICE, typename FPP2>
+void GivensFGFTParallel<FPP,DEVICE,FPP2>::choose_pivot()
+{
+	// fact_nz_inds is assumed to be sorted in max_L()
+	// the first one in the list maps the biggest abs val in L
+	pair<int,int> max_elt = *(fact_nz_inds.begin());
+	this->p = max_elt.first;
+	this->q = max_elt.second;
+	//keep only the last pivot selected as in Matlab version
+	this->coord_choices[this->ite] = pair<int,int>(this->p,this->q);
+	cout << "choose_pivot() p: " << this->p+1 << " q:" << this->q+1 << " " << "L(p,q): " << this->L(this->p,this->q) << " nrots: " << fact_nrots << endl;
+}
+
+
+template<typename FPP, Device DEVICE, typename FPP2>
+void GivensFGFTParallel<FPP,DEVICE,FPP2>::update_fact()
+{
+	if(fact_nrots == 0){
+		int n = this->Lap.getNbRow();
+		// forget previous rotation coeffs
+		// and keep identity part (n first coeffs)
+		this->fact_mod_row_ids.resize(n);
+		this->fact_mod_col_ids.resize(n);
+		this->fact_mod_values.resize(n);
+	}
+	// write new coeffs
+	// 1st one
+	this->fact_mod_row_ids.push_back(this->p);
+	this->fact_mod_col_ids.push_back(this->p);
+	this->fact_mod_values.push_back(cos(this->theta));
+	// 2nd
+	this->fact_mod_row_ids.push_back(this->p);
+	this->fact_mod_col_ids.push_back(this->q);
+	this->fact_mod_values.push_back(-sin(this->theta));
+	// 3rd
+	this->fact_mod_row_ids.push_back(this->q);
+	this->fact_mod_col_ids.push_back(this->p);
+	this->fact_mod_values.push_back(sin(this->theta));
+	// 4th
+	this->fact_mod_row_ids.push_back(this->q);
+	this->fact_mod_col_ids.push_back(this->q);
+	this->fact_mod_values.push_back(cos(this->theta));
+}
+
+
+template<typename FPP, Device DEVICE, typename FPP2>
+void GivensFGFTParallel<FPP,DEVICE,FPP2>::finish_fact()
+{
+	int n = this->Lap.getNbRow();
+	this->facts[this->ite] = MatSparse<FPP,DEVICE>(this->fact_mod_row_ids, this->fact_mod_col_ids, this->fact_mod_values, n, n);
+//#ifdef DEBUG_GIVENS
+	cout << "GivensFGFTParallel::finish_fact() ite: " << this->ite << " fact norm: " << this->facts[this->ite].norm() << endl;
+	this->facts[this->ite].Display();
+//#endif
+}