Update mex C++ svdtj to sort singular values in descendant order.

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

@@ -29,15 +29,87 @@
 %>
 %====================================================================
 function [U,S,V] = svdtj(M, maxiter, varargin)
-	[W1,D1] = matfaust.fact.eigtj(M*M', maxiter, varargin{:});
-	[W2,D2] = matfaust.fact.eigtj(M'*M, maxiter, varargin{:});
-	S = diag(W1'*M*W2);
-	[~,I] = sort(abs(S), 'descend');
-	S = sparse(diag(S(I)));
-	sign_S = sign(S);
-	S = S*sign_S;
-	Id = eye(size(S));
-	U = W1(:,1:size(Id,1))*matfaust.Faust({Id(:,I),sign_S});
-	V = W2(:,1:size(Id,1))*matfaust.Faust(Id(:,I));
+%	[W1,D1] = matfaust.fact.eigtj(M*M', maxiter, varargin{:});
+%	[W2,D2] = matfaust.fact.eigtj(M'*M, maxiter, varargin{:});
+%	S = diag(W1'*M*W2);
+%	[~,I] = sort(abs(S), 'descend');
+%	S = sparse(diag(S(I)));
+%	sign_S = sign(S);
+%	S = S*sign_S;
+%	Id = eye(size(S));
+%	U = W1(:,1:size(Id,1))*matfaust.Faust({Id(:,I),sign_S});
+%	V = W2(:,1:size(Id,1))*matfaust.Faust(Id(:,I));
+%	TODO: factorize argument parsing code with fgft_givens
+	import matfaust.Faust
+	nGivens_per_fac = 1; % default value
+	verbosity = 0; % default value
+%	if(~ ismatrix(M) || ~ isreal(M))
+%		error('M must be a real matrix.')
+%	end
+	if(size(M,1) ~= size(M,2))
+		error('M must be square')
+	end
+	if(~ isnumeric(maxiter) || maxiter-floor(maxiter) > 0 || maxiter <= 0)
+		error('maxiter must be a positive integer.')
+	end
+	bad_arg_err = 'bad number of arguments.';
+	tol = 0;
+	relerr = true;
+	verbosity = 0;
+	argc = length(varargin);
+	order = 1; % ascending order
+	if(argc > 0)
+		for i=1:argc
+			switch(varargin{i})
+				case 'tol'
+					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
+					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};
+					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}));
+					end
+				case '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);
+					end
+				case 'order'
+					if(argc == i || (~ strcmp(varargin{i+1}, 'ascend') && ~ strcmp(varargin{i+1}, 'descend') && ~ strcmp(varargin{i+1}, 'undef')))
+						error('order must be followed by a char array among ''ascend'', ''descend'' or ''undef''.')
+					else
+						order = varargin{i+1};
+						if(order(1) == 'a')
+							order = 1
+						elseif(order(1) == 'd')
+							order = -1
+						else
+							order = 0
+						end
+					end
+				otherwise
+					if(isstr(varargin{i}) && (~ strcmp(varargin{i}, 'ascend') && ~ strcmp(varargin{i}, 'descend') && ~ strcmp(varargin{i}, 'undef')) )
+						error([ varargin{i} ' unrecognized argument'])
+					end
+			end
+		end
+	end
+	[core_obj1, S, core_obj2] = mexsvdtjReal(M, maxiter, nGivens_per_fac, verbosity, tol, relerr, order);
+	S = sparse(diag(real(S)));
+	U = Faust(core_obj1, isreal(M));
+	V = Faust(core_obj2, isreal(M));
 end
 % experimental block end

wrapper/matlab/src/mexsvdtj.cpp.in

@@ -120,7 +120,7 @@ void svdtj(const mxArray* matlab_matrix, int J, int t, double tol, unsigned int
 			else
 			{
 				algoW1 = new GivensFGFTParallel<SCALAR,Cpu,FPP2>(dMM_, J, t, verbosity, tol, relErr);
-				algoW2 = new GivensFGFT<SCALAR,Cpu,FPP2>(dM_M, J, verbosity, tol, relErr);
+				algoW2 = new GivensFGFTParallel<SCALAR,Cpu,FPP2>(dM_M, J, t, verbosity, tol, relErr);
 			}
 		}else
 		{
@@ -136,12 +136,12 @@ void svdtj(const mxArray* matlab_matrix, int J, int t, double tol, unsigned int
 			else
 			{
 				algoW1 = new GivensFGFTParallel<SCALAR,Cpu,FPP2>(dMM_, J, t, verbosity, tol, relErr);
-				algoW2 = new GivensFGFT<SCALAR,Cpu,FPP2>(dM_M, J, verbosity, tol, relErr);
+				algoW2 = new GivensFGFTParallel<SCALAR,Cpu,FPP2>(dM_M, J, t, verbosity, tol, relErr);
 			}
 
 		}
 
-		//TODO: parallelize
+		//TODO: parallelize with OpenMP
 		algoW1->compute_facts();
 		algoW2->compute_facts();
 
@@ -153,21 +153,41 @@ void svdtj(const mxArray* matlab_matrix, int J, int t, double tol, unsigned int
 
 		plhs[0] = convertPtr2Mat<Faust::TransformHelper<SCALAR, Cpu>>(thW1);
 
-		Faust::Transform<SCALAR,Cpu> transW2 = std::move(algoW1->get_transform(order));
+		Faust::Transform<SCALAR,Cpu> transW2 = std::move(algoW2->get_transform(order));
 		TransformHelper<SCALAR,Cpu> *thW2 = new TransformHelper<SCALAR,Cpu>(transW2, 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[2] = convertPtr2Mat<Faust::TransformHelper<SCALAR, Cpu>>(thW2);
 
 		// compute S = W1'*M*W2 = W1'*(W2^T*M)^T
-		thW2->transpose();
+		dM.transpose();
 		Faust::MatDense<SCALAR,Cpu> MW2 = thW2->multiply(dM, /* transpose */ true);
-		thW1->transpose();
+		MW2.transpose();
 		Faust::MatDense<SCALAR,Cpu> W1_MW2 = thW1->multiply(MW2, /* transpose */ true);
+
 		// create diagonal vector
-		for(int i=0;i<S.size();i++)
+		for(int i=0;i<S.size();i++){
 			S.getData()[i] = W1_MW2(i,i);
+			std::cout << "% " << S.getData()[i] << std::endl;
+		}
 
-		plhs[1] = FaustVec2mxArray(S);
+		//order D descendently according to the abs value
+		vector<int> ord_indices;
+		Faust::Vect<SCALAR,Cpu> ordered_S = Faust::Vect<SCALAR,Cpu>(S.size());
+		ord_indices.resize(0);
+		order = 1;
+		for(int i=0;i<S.size();i++)
+			ord_indices.push_back(i);
+		sort(ord_indices.begin(), ord_indices.end(), [S, &order](int i, int j) {
+				return Faust::fabs(S.getData()[i]) > Faust::fabs(S.getData()[j])?1:0;
+				});
+		for(int i=0;i<ord_indices.size();i++)
+		{
+			ordered_S.getData()[i] = S.getData()[ord_indices[i]];
+			std::cout << "* " << ordered_S.getData()[i] << std::endl;
+		}
+		// TODO: and change the sign when the value is negative
+		// it gives a signed permutation matrix to append to W1
+		plhs[1] = FaustVec2mxArray(ordered_S);
 
 		delete algoW1;
 		delete algoW2;