Replace use of MatSparse by Vect for D in Faust::GivensFGFT* algorithms.

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

@@ -73,7 +73,7 @@ int main()
 
 	cout << "fourier_diag fro norm: " << fourier_diag.norm() << endl;
 
-	const Faust::MatSparse<FPP,Cpu>& D = algo.get_D();
+	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;
 
@@ -91,7 +91,7 @@ int main()
 	assert(fourier_diag2.norm()==0);
 	assert(ordered_fourier_diag2.norm()==0);
 
-	Faust::MatDense<FPP,Cpu> ordered_D(algo.get_D(true));
+	Faust::MatDense<FPP,Cpu> ordered_D(algo.get_Dspm(true));
 
 	cout << "orderded D fro norm: " << ordered_D.norm() << endl;
 	cout << "ordered_D:" << endl;

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

@@ -63,7 +63,7 @@ int main()
 
 	cout << "fourier_diag fro norm: " << fourier_diag.norm() << endl;
 
-	const Faust::MatSparse<FPP,Cpu> D = algo.get_D();
+	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;
 
@@ -81,7 +81,7 @@ int main()
 	assert(fourier_diag2.norm()==0);
 	assert(ordered_fourier_diag2.norm()==0);
 
-	Faust::MatDense<FPP,Cpu> ordered_D(algo.get_D(true));
+	Faust::MatDense<FPP,Cpu> ordered_D(algo.get_Dspm(true));
 
 	cout << "orderded D fro norm: " << ordered_D.norm() << endl;
 	cout << "ordered_D:" << endl;

src/algorithm/factorization/faust_GivensFGFT.h

@@ -38,7 +38,7 @@ namespace Faust {
 				/** \brief Fourier matrix factorization matrices (Givens matrix). */
 				vector<Faust::MatSparse<FPP,DEVICE>> facts;
 				/** \brief Diagonalization approximate of Laplacian. */
-				Faust::MatSparse<FPP,DEVICE> D;
+				Faust::Vect<FPP,DEVICE> D;
 				/** \brief Queue of errors (cf. calc_err()). */
 				vector<FPP2> errs;
 				/** \brief Pivot choices (p, q) coordinates. */
@@ -65,7 +65,7 @@ namespace Faust {
 				/** \brief Ordered indices of D to get increasing eigenvalues along the diagonal. */
 				vector<int> ord_indices;
 				/** \brief Cache for the ordered D. */
-				Faust::MatSparse<FPP,DEVICE> ordered_D;
+				Faust::Vect<FPP,DEVICE> ordered_D;
 				/** \brief true if D has already been ordered (order_D() was called). */
 				bool is_D_ordered;
 
@@ -197,9 +197,16 @@ namespace Faust {
 				FPP2 get_err(int j) const;
 
 				/**
-				 * Returns the matrix D in its current status (which is updated at each iteration).
+				 * Returns the diag. vector D in its current status (which is updated at each iteration).
 				 */
-				const MatSparse<FPP,DEVICE>& get_D(const bool ord=false);
+				const Vect<FPP,DEVICE>& get_D(const bool ord=false);
+
+				/** \brief Returns the diagonal vector as a sparse matrix.
+				 *
+				 * \note This function makes copies, is not intented for repeated use (use get_D() for optimized calls).
+				 *
+				 **/
+				const Faust::MatSparse<FPP,DEVICE> get_Dspm(const bool ord=false);
 
 				/**
 				 * Returns the diagonal by copying it in the buffer diag_data (should be allocated by the callee).

src/algorithm/factorization/faust_GivensFGFT.hpp

@@ -212,7 +212,7 @@ void GivensFGFT<FPP,DEVICE,FPP2>::update_D()
 {
 	// D = spdiag(diag(L))
 	for(int i=0;i<L.getNbRow();i++)
-		D.setCoeff(i,i, L(i,i));
+		D.getData()[i] = L(i,i);
 #ifdef DEBUG_GIVENS
 	D.Display();
 	cout << "D fro. norm: " << D.norm() << endl;
@@ -233,7 +233,7 @@ void GivensFGFT<FPP,DEVICE,FPP2>::update_err()
 	//
 	if(!((ite+1)%100) || verbosity > 1)
 	{
-		MatDense<FPP,Cpu> tmp = D;
+		MatDense<FPP,Cpu> tmp = this->get_Dspm(false);
 		tmp -= L;
 		FPP2 err = tmp.norm(), err_d;
 		err *= err;
@@ -249,18 +249,15 @@ void GivensFGFT<FPP,DEVICE,FPP2>::update_err()
 template<typename FPP, Device DEVICE, typename FPP2>
 void GivensFGFT<FPP,DEVICE,FPP2>::order_D()
 {
-	vector<FPP> D_ord_diag;
-	vector<int> nat_ord_indices;
-	for(int i=0;i<D.getNbRow();i++)
+	ordered_D = Vect<FPP,DEVICE>(D.size());
+	ord_indices.resize(0);
+	for(int i=0;i<D.size();i++)
 		ord_indices.push_back(i);
-	nat_ord_indices = ord_indices;
 	sort(ord_indices.begin(), ord_indices.end(), [this](int i, int j) {
-			return D.getValuePtr()[i] < D.getValuePtr()[j];
+			return D.getData()[i] < D.getData()[j];
 			});
-	for(int i=0;i<ord_indices.size();i++){
-		D_ord_diag.push_back(this->D.getValuePtr()[ord_indices[i]]);
-	}
-	ordered_D = MatSparse<FPP,Cpu>(nat_ord_indices, nat_ord_indices, D_ord_diag, nat_ord_indices.size(), nat_ord_indices.size());
+	for(int i=0;i<ord_indices.size();i++)
+		ordered_D.getData()[i] = D.getData()[ord_indices[i]];
 	is_D_ordered = true;
 }
 
@@ -286,7 +283,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, unsigned int verbosity /* deft val == 0 */) : Lap(Lap), facts(J), D(Lap.getNbRow(), Lap.getNbCol()), C(Lap.getNbRow(), Lap.getNbCol()), errs(0), coord_choices(J), L(Lap), q_candidates(new int[Lap.getNbCol()]), is_D_ordered(false), always_theta2(false), verbosity(verbosity)
+GivensFGFT<FPP,DEVICE,FPP2>::GivensFGFT(Faust::MatDense<FPP,DEVICE>& Lap, int J, unsigned int verbosity /* deft val == 0 */) : Lap(Lap), facts(J), D(Lap.getNbRow()), C(Lap.getNbRow(), Lap.getNbCol()), errs(0), coord_choices(J), L(Lap), q_candidates(new int[Lap.getNbCol()]), is_D_ordered(false), always_theta2(false), verbosity(verbosity)
 {
 	/* Matlab ref. code:
 	 *     facts = cell(1,J);
@@ -312,7 +309,7 @@ GivensFGFT<FPP,DEVICE,FPP2>::GivensFGFT(Faust::MatDense<FPP,DEVICE>& Lap, int J,
 	}
 
 	// init. D
-	D.setZeros();
+	memset(D.getData(), 0, sizeof(FPP)*Lap.getNbRow());
 }
 
 template<typename FPP, Device DEVICE, typename FPP2>
@@ -331,7 +328,7 @@ const vector<FPP2>& GivensFGFT<FPP,DEVICE,FPP2>::get_errs() const
 }
 
 template<typename FPP, Device DEVICE, typename FPP2>
-const Faust::MatSparse<FPP,DEVICE>& GivensFGFT<FPP,DEVICE,FPP2>::get_D(const bool ord /* default to false */)
+const Faust::Vect<FPP,DEVICE>& GivensFGFT<FPP,DEVICE,FPP2>::get_D(const bool ord /* default to false */)
 {
 	if(ord)
 	{
@@ -343,12 +340,27 @@ const Faust::MatSparse<FPP,DEVICE>& GivensFGFT<FPP,DEVICE,FPP2>::get_D(const boo
 	return D;
 }
 
+template<typename FPP, Device DEVICE, typename FPP2>
+const Faust::MatSparse<FPP,DEVICE> GivensFGFT<FPP,DEVICE,FPP2>::get_Dspm(const bool ord /* default to false */)
+{
+	const Faust::Vect<FPP,DEVICE>& D_ = this->get_D(ord);
+	vector<int> nat_ord_indices;
+	vector<FPP> diag_D;
+	for(int i=0;i<D_.size();i++)
+	{
+		nat_ord_indices.push_back(i);
+		diag_D.push_back(D_.getData()[i]);
+	}
+	MatSparse<FPP,DEVICE> spD(nat_ord_indices, nat_ord_indices, diag_D, nat_ord_indices.size(), nat_ord_indices.size());
+	return spD;
+}
+
 template<typename FPP, Device DEVICE, typename FPP2>
 void GivensFGFT<FPP,DEVICE,FPP2>::get_D(FPP* diag_data, const bool ord /* default to false */)
 {
-	get_D(ord);
-	FPP* src_data_ptr = ord?ordered_D.getValuePtr():D.getValuePtr();
-	memcpy(diag_data, src_data_ptr, sizeof(FPP)*D.getNbRow());
+	const Faust::Vect<FPP,DEVICE>& D_ = get_D(ord);
+	const FPP* src_data_ptr = D_.getData();
+	memcpy(diag_data, src_data_ptr, sizeof(FPP)*D_.size());
 }