Add a new version of prox_sp that authorizes dynamically to output a MatSparse or MatDense image.

The test shows on an exemple that it is not slower:
./run_test/bin/faust_sparse_prox_sp_double
sparse sp prox time:0.340821
dense sp prox time:0.369167

This commit is issue #179 related.

misc/test/CMakeLists.txt

@@ -199,7 +199,7 @@ if(MATIO_LIB_FILE AND MATIO_INC_DIR AND BUILD_READ_MAT_FILE AND NOT NOCPPTESTS)
 	# faust_multiplication : time comparison between Faust-vector product and Dense matrix-vector product
 
 
-	list(APPEND tests hierarchicalFactorization hierarchicalFactorizationFFT test_palm4MSA test_palm4MSAFFT faust_multiplication  faust_matdense_conjugate GivensFGFT GivensFGFTSparse GivensFGFTParallel GivensFGFTParallelSparse test_MatDiag faust_matsparse_mul faust_matsparse_index_op GivensFGFTComplex GivensFGFTComplexSparse GivensFGFTParallelComplex faust_toeplitz faust_circ faust_hankel palm4msa_2020 hierarchical2020 hierarchical2020Hadamard hierarchicalFactorizationHadamard)
+	list(APPEND tests hierarchicalFactorization hierarchicalFactorizationFFT test_palm4MSA test_palm4MSAFFT faust_multiplication  faust_matdense_conjugate GivensFGFT GivensFGFTSparse GivensFGFTParallel GivensFGFTParallelSparse test_MatDiag faust_matsparse_mul faust_matsparse_index_op GivensFGFTComplex GivensFGFTComplexSparse GivensFGFTParallelComplex faust_toeplitz faust_circ faust_hankel faust_sparse_prox_sp palm4msa_2020 hierarchical2020 hierarchical2020Hadamard hierarchicalFactorizationHadamard)
 
 	if(FAUST_TORCH)
 		list(APPEND tests faust_torch)

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

+#include "faust_MatDense.h"
+#include "faust_MatSparse.h"
+#include "faust_prox.h"
+#include <iostream>
+
+/** \brief unitary test for MatDense conjugate
+*/
+
+typedef @TEST_FPP@ FPP;
+
+
+using namespace Faust;
+using namespace std;
+
+void print_mat_data(MatDense<FPP,Cpu>& mat)
+{
+	int nrows = mat.getNbRow(), ncols = mat.getNbCol();
+	for (int i=0;i < nrows; i++)
+	{
+		for(int j=0;j < ncols; j++)
+			cout << mat.getData()[j*nrows+i] << " ";
+		cout << endl;
+	}
+}
+
+int main(int argc, char* argv[])
+{
+	faust_unsigned_int dim1 = 204;
+	faust_unsigned_int dim2 = 204;
+
+	MatDense<FPP,Cpu>* M;
+	MatDense<FPP,Cpu> M_copy;
+	M = MatDense<FPP,Cpu>::randMat(dim1,dim2);
+	*M *= FPP(100);
+	MatDense<FPP,Cpu> dM = *M;
+	MatSparse<FPP,Cpu> spM;
+	MatGeneric<FPP,Cpu>* G;
+
+	auto num_its = 300;
+
+	auto k = 1632;
+
+
+	std::chrono::time_point<std::chrono::steady_clock> sparse_prox_start, sparse_prox_end, dense_prox_start, dense_prox_end;
+	std::chrono::duration<double> sparse_prox_dur, dense_prox_dur;
+
+	for(int i=0; i < num_its; i++)
+	{
+		sparse_prox_start = std::chrono::steady_clock::now();
+		G = prox_sp(dM, spM, k, /* normalized */ false);
+		sparse_prox_end = std::chrono::steady_clock::now();
+//		G->Display();
+//		cout << G->norm() << endl;
+		MatDense<FPP,Cpu> spM2d(*dynamic_cast<MatSparse<FPP,Cpu>*>(G));
+//		spM2d.Display();
+//		print_mat_data(spM2d);
+		sparse_prox_dur += sparse_prox_end-sparse_prox_start;
+		dense_prox_start = std::chrono::steady_clock::now();
+		dM = *M;
+		prox_sp(dM, k, /* normalized */ false);
+		spM = dM;
+		dense_prox_end = std::chrono::steady_clock::now();
+//		cout <<dM.norm() << endl;
+//		print_mat_data(dM);
+//		dM.Display();
+		dense_prox_dur += dense_prox_end-dense_prox_start;
+		dM -= spM2d;
+		auto err = dM.norm();
+//		cout << "err:" << err << endl;
+		assert(err < 1e-6);
+		dM = *M;
+	}
+
+	cout << "sparse sp prox time:" << sparse_prox_dur.count() << endl;
+	cout << "dense sp prox time:" << dense_prox_dur.count() << endl;
+
+	return 0;
+}

src/faust_linear_operator/CPU/faust_prox.h

@@ -100,7 +100,15 @@ namespace Faust
 	template<typename FPP>
 		void prox_skperm(MatDense<FPP, Cpu> & M,const unsigned int k,  const bool normalized=true, const bool pos=false);
 
-
+	template<typename FPP> faust_unsigned_int sparse_size(faust_unsigned_int nnz, faust_unsigned_int nrows);
+	template<typename FPP> faust_unsigned_int dense_size(faust_unsigned_int nrows, faust_unsigned_int ncols);
+	/**
+	 * Decides which output format to use when appliying the SP prox. op. to M. Either M or spM is the output, it depends on the byte size. The minimum memory fingerprint is targeted.
+	 *
+	 *  \param forcedType: used to choose explicitely the output format with values Sparse or Dense (MatSparse or MatDense).
+	 */
+	template<typename FPP>
+		MatGeneric<FPP,Cpu>* prox_sp(MatDense<FPP,Cpu> & M, MatSparse<FPP, Cpu> & spM, faust_unsigned_int k, const bool normalized=true, const bool pos=false, const MatType forcedType=None);
 }
 
 #include "faust_prox.hpp"

src/faust_linear_operator/CPU/faust_prox.hpp

@@ -94,6 +94,94 @@ void Faust::prox_sp(Faust::MatDense<FPP,Cpu> & M,faust_unsigned_int k, const boo
 	if(normalized) M.normalize();
 }
 
+template<typename FPP> faust_unsigned_int Faust::sparse_size(faust_unsigned_int nnz, faust_unsigned_int nrows)
+{
+	//TODO: move elsewhere
+	auto size = nnz*(sizeof(FPP)+sizeof(int)) + (nrows + 1) * sizeof(int);
+	return size;
+}
+
+template<typename FPP> faust_unsigned_int Faust::dense_size(faust_unsigned_int nrows, faust_unsigned_int ncols)
+{
+	//TODO: move elsewhere
+	auto size = nrows*ncols*sizeof(FPP);
+	return size;
+}
+
+template<typename FPP>
+Faust::MatGeneric<FPP,Cpu>* Faust::prox_sp(Faust::MatDense<FPP,Cpu> & M, Faust::MatSparse<FPP, Cpu> & spM, faust_unsigned_int k, const bool normalized /* true by deft */, const bool pos, const MatType forcedType/*=None*/)
+{
+	// M is the dense matrix on which to compute projection
+	// M is also the output matrix if k is high enough
+	// otherwise spM is used as output
+	// The output matrix is returned as a MatGeneric (as explained it could be MatSparse or MatDense)
+	const faust_unsigned_int dim1 = M.getNbRow();
+	const faust_unsigned_int dim2 = M.getNbCol();
+	const faust_unsigned_int nnz = M.getNonZeros();
+
+	auto out_is_dense = sparse_size<FPP>(k, dim1) > dense_size<FPP>(dim1, dim2) && forcedType == None || forcedType == Dense;
+
+	if(pos) Faust::pre_prox_pos(M);
+
+	const faust_unsigned_int numel = dim1*dim2;
+
+	if (k < numel /* && k < M.getNonZeros()*/)
+	{
+		const std::vector<FPP> vec(M.getData(), M.getData()+numel);
+		std::vector<int> index;
+		Faust::sort_idx(vec, index, k);
+		index.erase(index.begin()+k, index.end());
+
+		if(out_is_dense)
+		{
+			M.setZeros();
+			for (int i=0 ; i<index.size() ; i++)
+				M.getData()[index[i]] = vec[index[i]];
+
+			if(normalized) M.normalize();
+			return &M;
+		}
+		else
+		{
+			FPP *values = new FPP[index.size()];
+			unsigned int *row_ids = new unsigned int[index.size()];
+			unsigned int *col_ids = new unsigned int[index.size()];
+
+
+
+			std::vector<Eigen::Triplet<FPP> > tripletList;
+			for (int i=0 ; i<index.size() ; i++)
+			{
+				values[i] = vec[index[i]];
+				row_ids[i] = index[i]%dim1;
+				col_ids[i] = index[i]/dim1;
+			}
+
+			if(normalized)
+			{
+				Faust::Vect<FPP,Cpu> nvalues(index.size(), values);
+				nvalues.normalize();
+				values = nvalues.getData();
+			}
+			spM = MatSparse<FPP,Cpu>(row_ids, col_ids, values, dim1, dim2, index.size());
+			delete[] row_ids;
+			delete[] col_ids;
+			delete[] values;
+			return &spM;
+		}
+	}
+	else
+	{
+		if(sparse_size<FPP>(M.getNonZeros(), dim1) > dense_size<FPP>(dim1, dim2))
+			return &M;
+		else
+		{
+			spM = M;
+			return &spM;
+		}
+	}
+}
+
 template<typename FPP>
 void Faust::prox_spcol(Faust::MatDense<FPP,Cpu> & M,faust_unsigned_int k, const bool normalized /* deft to true */, const bool pos)
 {

src/faust_linear_operator/faust_constant.h.in

@@ -92,7 +92,8 @@ enum MatType
 {
    Dense, 
    Sparse, 
-   Diag
+   Diag,
+   None
 };