Add a test for GivensFGFTParallelComplex (upcoming code).

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

+
+#include "faust_GivensFGFTParallelComplex.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<complex<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)
+	GivensFGFTParallelComplex<complex<FPP>,Cpu, FPP> 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 eq=%d\n", j, p_choices[j], q_choices[j], coord_choices[j].first+1, coord_choices[j].second+1, p_choices[j] == coord_choices[j].first+1 && q_choices[j] == coord_choices[j].second+1);
+
+
+	Faust::MatDense<complex<FPP>,Cpu> fourier_diag(Lap.getNbRow(), Lap.getNbCol());
+	fourier_diag.setEyes();
+	const vector<Faust::MatSparse<complex<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<complex<FPP>,Cpu> D = algo.get_Dspm();
+	for(int i=0; i < Lap.getNbRow(); i++)
+	{
+		cout << "-- " << D(i,i) << endl;
+	}
+	Faust::MatDense<complex<FPP>,Cpu> full_D = Faust::MatDense<complex<FPP>,Cpu>(D);
+	cout << "D fro norm:" << D.norm() << endl;
+
+	Faust::MatDense<complex<FPP>,Cpu> fourier_diag2 = algo.compute_fourier();
+	Faust::MatDense<complex<FPP>,Cpu> ordered_fourier_diag2 = algo.compute_fourier(true);
+
+	Faust::MatDense<complex<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()==FPP(0));
+	assert(ordered_fourier_diag2.norm()==FPP(0));
+
+	Faust::MatDense<complex<FPP>,Cpu> ordered_D(algo.get_Dspm(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<complex<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;
+}
+