Update/Fix a few elements to validate GivensFGFTParallel C++ impl.

- Fix the deletion of pivot candidates in GivensFGFTParallel: when pivot (p,q) is selected, all pivot candidates (p',q') such that p'==p or q'==p or p'==q or q'==q must be deleted from candidates (as in the matlab script).
- Change the matlab ref. script and the C++ impl. to store all pivots (only the last one was kept).
- Change the Laplacian matrix again in test data mat file (return to the same used for GivensFGFT parent algorithm), store the good pivots (since we store all of them now).
- Update output for pivots in test GivensFGFTParallel.cpp.in.

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

@@ -52,7 +52,7 @@ int main()
 	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);
+		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());

src/algorithm/factorization/faust_GivensFGFTParallel.hpp

@@ -7,6 +7,7 @@ GivensFGFTParallel<FPP,DEVICE,FPP2>::GivensFGFTParallel(Faust::MatDense<FPP,DEVI
 {
 	this->facts.resize(round(J/(float)t));
 	this->always_theta2 = true;
+	this->coord_choices.resize(0);
 }
 
 template<typename FPP, Device DEVICE, typename FPP2>
@@ -26,7 +27,7 @@ void GivensFGFTParallel<FPP,DEVICE,FPP2>::max_L()
 
 #ifdef DEBUG_GIVENS
 	for(auto &p : fact_nz_inds)
-		cout << "GivensFGFTParallel::max_L() before sort :" << L_low(p.first, p.second) << endl;
+		cout << "GivensFGFTParallel::max_L() before sort (" << p.first+1 <<  "," << p.second+1 << ") :" << 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.
@@ -36,7 +37,7 @@ void GivensFGFTParallel<FPP,DEVICE,FPP2>::max_L()
 			});
 #ifdef DEBUG_GIVENS
 	for(auto &p : fact_nz_inds)
-		cout << "GivensFGFTParallel::max_L() after sort :" << L_low(p.first, p.second) << endl;
+		cout << "GivensFGFTParallel::max_L() after sort (" << p.first+1 <<  "," << p.second+1 << ") :" << L_low(p.first, p.second) << endl;
 #endif
 }
 
@@ -53,11 +54,15 @@ void GivensFGFTParallel<FPP,DEVICE,FPP2>::update_fact_nz_inds()
 	//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)
+		if((*i).first == this->p || (*i).second == this->q || (*i).first == this->q || (*i).second == this->p)
 			i = fact_nz_inds.erase(i);
 		else
 			i++;
 	}
+#ifdef DEBUG_GIVENS
+	for(auto &p : fact_nz_inds)
+		cout << "GivensFGFTParallel::update_fact_nz_inds() after purge (" << p.first+1 <<  "," << p.second+1 << ") :" << endl;
+#endif
 }
 
 template<typename FPP, Device DEVICE, typename FPP2>
@@ -105,9 +110,10 @@ void GivensFGFTParallel<FPP,DEVICE,FPP2>::choose_pivot()
 	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);
+	this->coord_choices.push_back(pair<int,int>(this->p,this->q));
+#ifdef DEBUG_GIVENS
 	cout << "choose_pivot() p: " << this->p+1 << " q:" << this->q+1 << " " << "L(p,q): " << this->L(this->p,this->q) << " nrots: " << fact_nrots << endl;
+#endif
 }
 
 
@@ -147,8 +153,8 @@ 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
+#ifdef DEBUG_GIVENS
 	cout << "GivensFGFTParallel::finish_fact() ite: " << this->ite << " fact norm: " << this->facts[this->ite].norm() << endl;
 	this->facts[this->ite].Display();
-//#endif
+#endif
 }

wrapper/matlab/tools/diagonalization_givens_parall.m

@@ -51,7 +51,8 @@ n=size(Lap,1);
 L=sparse(Lap);
 C = zeros(n);
 err=zeros(1,round(J/t));
-coord_choices = zeros(2,J);
+%coord_choices = zeros(2,J);
+coord_choices = []
 %N_edges = (nnz(Lap)-n)/2;
 
 
@@ -64,8 +65,7 @@ for j=1:round(J/t)
     %%%%%%%%% Un tri
     ind_nnz = find(L_low);
     %disp(['nnz in L :' num2str(numel(ind_nnz))])
-    
-    [~,ind_sorted] = sort(L_low(ind_nnz),'descend');
+	[~,ind_sorted] = sort(L_low(ind_nnz),'descend');
     [Rvec, Svec] = ind2sub([n,n],ind_nnz(ind_sorted));
     RSmat = [Rvec, Svec];
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -116,8 +116,9 @@ for j=1:round(J/t)
             RSmatnew = RSmat(~tokill,:);
             RSmat = RSmatnew;
             % if  ~any(abs(p-chosen)<0.5) && ~any(abs(q-chosen)<0.5)%numel(find(S(p,:)))==1 && numel(find(S(:,q)))==1
-            coord_choices(1,j) = p;
-            coord_choices(2,j) = q;
+            %coord_choices(1,j) = p;
+            %coord_choices(2,j) = q;
+			coord_choices = [ coord_choices, [p;q]];
             theta = atan2(L(q,q) - L(p,p),2*L(p,q))/2 + pi/4;
             S(p,p) = cos(theta); S(p,q) = -sin(theta);
             S(q,p) = sin(theta); S(q,q) = cos(theta);
@@ -130,7 +131,8 @@ for j=1:round(J/t)
     end
     L = S'*L*S;
     facts{j} = sparse(S);
-    D = spdiag(diag(L));
+	D = spdiag(diag(L)); 
+    % D = sparse(diag(diag(L)));
     %    errchg = err - norm(D-L,'fro')^2
     if mod(j,1)==0
         %err(j) = norm(D-L,'fro')^2/norm(L,'fro')^2;