Add to matfaust eigtj/fgft_givens an order argument to define the eigenvalues sort.

src/algorithm/factorization/faust_GivensFGFT.h

@@ -69,6 +69,7 @@ namespace Faust {
 				Faust::Vect<FPP,DEVICE> ordered_D;
 				/** \brief true if D has already been ordered (order_D() was called). */
 				bool is_D_ordered;
+				int D_order_dir;
 
 				/** \brief The level of verbosity (0 for nothing, 1 for iteration numbers,...) */
 				unsigned int verbosity;
@@ -196,10 +197,18 @@ namespace Faust {
 				void update_err();
 
 				/**
-				 * Order D into ordered_D and keeps ordered indices in ord_indices.
+				 * Order (sort) D ascendantly into ordered_D and keeps ordered indices in ord_indices.
 				 */
 				void order_D();
 
+				/**
+				 * Order D into ascendantly or descendantly into order_D and keeps ordered indices in ord_indices.
+				 *
+				 * \param order -1 to order descendantly, 1 to order ascendantly.
+				 */
+				void order_D(int order);
+
+
 			public:
 
 				/**
@@ -225,6 +234,12 @@ namespace Faust {
 				 */
 				const Faust::Vect<FPP,DEVICE>& get_D(const bool ord=false);
 
+				/**
+				 * \param ord 0 (no sort), 1 (ascendant sort), -1 (descendant sort).
+				 */
+				const Faust::Vect<FPP,DEVICE>& get_D(const int ord=0);
+
+
 				/** \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).
@@ -239,6 +254,12 @@ namespace Faust {
 				 */
 				void get_D(FPP* diag_data, const bool ord=false);
 
+				/**
+				 *
+				 * \param ord: 0 for no sort, -1 for descending order, 1 for descending order.
+				 */
+				void get_D(FPP* diag_data, const int ord=0);
+
 
 				/**
 				 * Computes and returns the Fourier matrix approximate from the Givens factors computed up to this time.
@@ -283,6 +304,10 @@ namespace Faust {
 				 * \param ord true to get the Transform's facts ordering the last one columns according to ascendant eigenvalues, false to let facts as they went out from the algorithm (without reordering).
 				 */
 				Faust::Transform<FPP,DEVICE> get_transform(bool ord);
+				/**
+				 * \param ord -1 for descending order, 1 for ascending order.
+				 */
+				Faust::Transform<FPP,DEVICE> get_transform(int ord);
 
 
 		};

src/algorithm/factorization/faust_GivensFGFT.hpp

@@ -402,17 +402,24 @@ void GivensFGFT<FPP,DEVICE,FPP2>::update_err()
 
 template<typename FPP, Device DEVICE, typename FPP2>
 void GivensFGFT<FPP,DEVICE,FPP2>::order_D()
+{
+	order_D(1);
+}
+
+template<typename FPP, Device DEVICE, typename FPP2>
+void GivensFGFT<FPP,DEVICE,FPP2>::order_D(const int order /* -1 for descending order, 1 for ascending order */)
 {
 	ordered_D = Faust::Vect<FPP,DEVICE>(D.size());
 	ord_indices.resize(0);
 	for(int i=0;i<D.size();i++)
 		ord_indices.push_back(i);
-	sort(ord_indices.begin(), ord_indices.end(), [this](int i, int j) {
-			return D.getData()[i] < D.getData()[j];
+	sort(ord_indices.begin(), ord_indices.end(), [this, &order](int i, int j) {
+			return order>0?D.getData()[i] < D.getData()[j]:D.getData()[i] > D.getData()[j];
 			});
 	for(int i=0;i<ord_indices.size();i++)
 		ordered_D.getData()[i] = D.getData()[ord_indices[i]];
 	is_D_ordered = true;
+	D_order_dir = order;
 }
 
 template<typename FPP, Device DEVICE, typename FPP2>
@@ -528,7 +535,18 @@ const Faust::Vect<FPP,DEVICE>& GivensFGFT<FPP,DEVICE,FPP2>::get_D(const bool ord
 		if(!is_D_ordered)
 			order_D();
 		return ordered_D;
+	}
+	return D;
+}
 
+template<typename FPP, Device DEVICE, typename FPP2>
+const Faust::Vect<FPP,DEVICE>& GivensFGFT<FPP,DEVICE,FPP2>::get_D(const int ord /* default to false */)
+{
+	if(ord != 0)
+	{
+		if(!is_D_ordered || ord != D_order_dir)
+			order_D(ord);
+		return ordered_D;
 	}
 	return D;
 }
@@ -550,13 +568,18 @@ const Faust::MatSparse<FPP,DEVICE> GivensFGFT<FPP,DEVICE,FPP2>::get_Dspm(const b
 
 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(diag_data, ord?1:0);
+}
+
+template<typename FPP, Device DEVICE, typename FPP2>
+void GivensFGFT<FPP,DEVICE,FPP2>::get_D(FPP* diag_data, const int ord /* default to false */)
 {
 	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());
 }
 
-
 template<typename FPP, Device DEVICE, typename FPP2>
 const Faust::MatDense<FPP,DEVICE> GivensFGFT<FPP,DEVICE,FPP2>::compute_fourier(const bool ord /* default to false */)
 {
@@ -615,13 +638,20 @@ const vector<Faust::MatSparse<FPP,DEVICE>>& GivensFGFT<FPP,DEVICE,FPP2>::get_fac
 
 template<typename FPP, Device DEVICE, typename FPP2>
 Faust::Transform<FPP,DEVICE> GivensFGFT<FPP,DEVICE,FPP2>::get_transform(bool ord)
+{
+	return get_transform(ord?1:0);
+}
+
+
+template<typename FPP, Device DEVICE, typename FPP2>
+Faust::Transform<FPP,DEVICE> GivensFGFT<FPP,DEVICE,FPP2>::get_transform(int ord)
 {
 	//TODO: an optimization is possible by changing type of facts to vector<MatGeneric*> it would avoid copying facts into Transform and rather use them directly. It will need a destructor that deletes them eventually if they weren't transfered to a Transform object before.
 	if(ord)
 	{
 		// add a permutation factor to reorder according to ordered D
-		if(!is_D_ordered)
-			order_D();
+		if(ord != 0 && (!is_D_ordered || ord != D_order_dir))
+			order_D(ord);
 		MatSparse<FPP,DEVICE> & last_fact = *(facts.end()-1);
 		MatSparse<FPP,DEVICE> P(last_fact.getNbCol(), last_fact.getNbCol()); //last_fact is a sqr. mat.
 		P.set_orthogonal(true);

wrapper/matlab/+matfaust/+fact/eigtj.m

@@ -22,11 +22,12 @@
 %> @param 'relerr', true (optional) For a stopping criterion based on the relative squared error (this is the default error).
 %> @param 'relerr', false (optional) For a stopping criterion based on the absolute squared error.
 %> @param 'verbosity', integer (optional) the level of verbosity, the greater the value the more info. is displayed.
+%> @param 'order', integer (optional) -1 for a descending order of eigenvalues, 1 for an ascending order.
 %>
 %> @retval [V,D]
 %> - V the Faust object representing the approximate eigenvector transform. The column V(:, i) is the eigenvector corresponding to the eigenvalue D(i,i).
 %> The last factor of V is a permutation matrix. The goal of this factor is to apply to the columns of V the same order as eigenvalues set in D.
-%> - D the approximate sparse diagonal matrix of the eigenvalues (in ascendant order along the rows/columns).
+%> - D the approximate sparse diagonal matrix of the eigenvalues (by default in ascendant order along the rows/columns).
 %>
 %> @b Example
 %> @code

wrapper/matlab/+matfaust/+fact/fgft_givens.m

@@ -11,10 +11,11 @@
 %> @param 'tol', number see fact.eigtj
 %> @param 'relerr', bool see fact.eigtj
 %> @param 'verbosity', integer see fact.eigtj
+%> @param 'order', integer see fact.eigtj
 %>
 %> @retval [FGFT,D]:
 %> - with FGFT being the Faust object representing the Fourier transform and,
-%> -  D as a sparse diagonal matrix of the eigenvalues in ascendant order along the rows/columns.
+%> -  D as a sparse diagonal matrix of the eigenvalues by default in ascendant order along the rows/columns.
 %>
 %>
 %> @b References:
@@ -45,6 +46,7 @@ function [FGFT,D] = fgft_givens(Lap, maxiter, varargin)
 	relerr = true;
 	verbosity = 0;
 	argc = length(varargin);
+	order = 1; % ascending order
 	if(argc > 0)
 		for i=1:argc
 			switch(varargin{i})
@@ -52,27 +54,33 @@ function [FGFT,D] = fgft_givens(Lap, maxiter, varargin)
 					if(argc == i || ~ isscalar(varargin{i+1}))
 						error('tol keyword arg. is not followed by a number')
 					else
-						tol = real(varargin{i+1}) % real in case of cplx num
+						tol = real(varargin{i+1}); % real in case of cplx num
 					end
 				case 'relerr'
 					if(argc == i || ~ islogical(varargin{i+1}))
 						error('relerr keyword argument is not followed by a logical')
 					else
-						relerr = varargin{i+1}
+						relerr = varargin{i+1};
 					end
 				case 'verbosity'
 					if(argc == i || ~ isscalar(varargin{i+1}))
 						error('verbose keyword argument is not followed by a number')
 					else
-						verbosity = floor(real(varargin{i+1}))
+						verbosity = floor(real(varargin{i+1}));
 					end
 				case 'nGivens_per_fac'
-					if(argc == i || ~ isscalar(varargin{i+1}) || ~ isnumeric(nGivens_per_fac))
+					if(argc == i || ~ isscalar(varargin{i+1}) || ~ isnumeric(varargin{i+1}))
 						error('nGivens_per_fac must be followed by a positive integer.')
 					else
 						nGivens_per_fac = floor(abs(real(varargin{i+1})));
 						nGivens_per_fac = min(nGivens_per_fac, maxiter);
-						nGivens_per_fac = max(1, nGivens_per_fac)
+						nGivens_per_fac = max(1, nGivens_per_fac);
+					end
+				case 'order'
+					if(argc == i || ~ isscalar(varargin{i+1}) || ~ isnumeric(varargin{i+1}))
+						error('order must be followed by an integer.')
+					else
+						order = varargin{i+1};
 					end
 				otherwise
 					if(isstr(varargin{i}))
@@ -81,7 +89,7 @@ function [FGFT,D] = fgft_givens(Lap, maxiter, varargin)
 			end
 		end
 	end
-	[core_obj, D] = mexfgftgivensReal(Lap, maxiter, nGivens_per_fac, verbosity, tol, relerr);
+	[core_obj, D] = mexfgftgivensReal(Lap, maxiter, nGivens_per_fac, verbosity, tol, relerr, order);
 	D = sparse(diag(D));
 	FGFT = Faust(core_obj, true);
 end

wrapper/matlab/src/mexfgftgivens.cpp.in

@@ -59,8 +59,9 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 	unsigned int verbosity = 0; //default verbosity (no info. displayed)
 	double tol = 0;
 	bool relErr = true;
+	int order;
 
-	if(nrhs < 2 || nrhs > 6)
+	if(nrhs < 2 || nrhs > 7)
 		mexErrMsgTxt("Bad Number of inputs arguments");
 
 	J = (int) mxGetScalar(prhs[1]);
@@ -72,6 +73,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 		tol = (double) mxGetScalar(prhs[4]);
 	if(nrhs >= 6)
 		relErr = (bool) mxGetScalar(prhs[5]);
+	if(nrhs >= 7)
+		order = (int) mxGetScalar(prhs[6]); //eigenvalues order
 
 	const mxArray* matlab_matrix = prhs[0]; // Laplacian
 
@@ -112,11 +115,11 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 
 		Faust::Vect<SCALAR,Cpu> D(Lap->getNbRow());
 		// true below is for ascendant order of eigenvalues
-		algo->get_D(const_cast<SCALAR*>(D.getData()), true); // don't respect constness for optimization (saving a vector copy)
+		algo->get_D(const_cast<SCALAR*>(D.getData()), order); // don't respect constness for optimization (saving a vector copy)
 
 		plhs[1] = FaustVec2mxArray(D);
 
-		Faust::Transform<SCALAR,Cpu> trans = std::move(algo->get_transform(true)); // true for same order of "eigenvectors" as eigenvalues
+		Faust::Transform<SCALAR,Cpu> trans = std::move(algo->get_transform(order)); // true for same order of "eigenvectors" as eigenvalues
 		TransformHelper<SCALAR,Cpu> *th = new TransformHelper<SCALAR,Cpu>(trans, true); // true is for moving and not copying the Transform object into TransformHelper (optimization possible cause we know the original object won't be used later)
 
 		plhs[0] = convertPtr2Mat<Faust::TransformHelper<SCALAR, Cpu>>(th);