diff --git a/misc/test/src/C++/EigTJParallelSparse.cpp.in b/misc/test/src/C++/EigTJParallelSparse.cpp.in index 07006b8c336beb198ba13baeba23e4422ac602fa..e8476392187819f14e99a0fb53097508ba2cd07d 100644 --- a/misc/test/src/C++/EigTJParallelSparse.cpp.in +++ b/misc/test/src/C++/EigTJParallelSparse.cpp.in @@ -1,114 +1,117 @@ - -#include "faust_EigTJParallel.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 EigTJ. - */ - - -int main() -{ - string configFilename = "@FAUST_DATA_MAT_DIR@/test_GivensDiagParallel_Lap_U_J_choices.mat"; - - Faust::MatDense<FPP,Cpu> Lap; - Faust::MatSparse<FPP,Cpu> sLap; - init_faust_mat_from_matio(Lap,configFilename.c_str(),"Lap"); - - sLap = Lap; - int J = 5472; //TODO: should be retrieved from .mat file (not be a hard-coded value) - EigTJParallel<FPP,Cpu, FPP> algo(sLap,J,/*t=*/ Lap.getNbRow()/2, - /*verbosity */ 0, /* stoppingError */ 0.0, /* errIsRel */ true, /* enable_large_Faust */ true); - - 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<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_Dspm(); - 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_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<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; -} - + +#include "faust_EigTJParallel.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 EigTJ. + */ + + +int main() +{ + string configFilename = "@FAUST_DATA_MAT_DIR@/test_GivensDiagParallel_Lap_U_J_choices.mat"; + + Faust::MatDense<FPP,Cpu> Lap; + Faust::MatSparse<FPP,Cpu> sLap; + init_faust_mat_from_matio(Lap,configFilename.c_str(),"Lap"); + + sLap = Lap; + + + 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; + } + auto mv_info = Mat_VarReadInfo(matfp, "choices"); + matvar_t* cell_choices = Mat_VarRead(matfp, "choices"); +// Mat_VarPrint(cell_choices, 0); + int J = cell_choices->dims[1]; + EigTJParallel<FPP,Cpu, FPP> algo(sLap,J,/*t=*/ Lap.getNbRow()/2, + /*verbosity */ 0, /* stoppingError */ 0.0, /* errIsRel */ true, /* enable_large_Faust */ true); + + algo.compute_facts(); + + vector<pair<int,int>> coord_choices = algo.get_coord_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]; + } + Mat_VarFree(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<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_Dspm(); + 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_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<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; +} +