Add MatDense::nonzeros_indices() and a unit test.

misc/test/CMakeLists.txt

@@ -156,7 +156,7 @@ endif()
 
 
 foreach(TEST_FPP float double)
-	foreach(FILE faust_mult faust_mult_cplx test_Vect_min test_MatDense_min test_MatDense_get_row test_MatDense_lower_upper_tri)
+	foreach(FILE faust_mult faust_mult_cplx test_Vect_min test_MatDense_min test_MatDense_get_row test_MatDense_lower_upper_tri test_MatDense_nonzeros_indices)
 		set(TEST_BIN_FILE ${FILE}_${TEST_FPP})
 		set(TEST_FILE_CPP ${TEST_BIN_FILE}.cpp)
 		message(STATUS ${TEST_FILE_CPP})

misc/test/src/C++/unit/test_MatDense_nonzeros_indices.cpp.in

+#include "faust_MatDense.h"
+
+
+#include <complex>
+#include <random>
+
+using namespace Faust;
+typedef @TEST_FPP@ FPP;
+
+
+int main()
+{
+	int mat_nrows = 5, mat_ncols = 5;
+	int nz = 3;
+	int nnz = mat_nrows*mat_ncols-nz;
+	int * zero_inds = new int[nz*2];
+
+	//ele 0
+	zero_inds[0] = 1;
+	zero_inds[1] = 3;
+	//ele 1
+	zero_inds[1*2+0] = 0;
+	zero_inds[1*2+1] = 4;
+	//ele 2
+	zero_inds[2*2+0] = 3;
+	zero_inds[2*2+1] = 2;
+
+
+	MatDense<FPP,Cpu>* M = MatDense<FPP,Cpu>::randMat(mat_nrows, mat_ncols);
+	M->scalarMultiply(100);
+
+	assert(M->getNonZeros() == mat_nrows*mat_ncols);
+
+	for(int i=0;i<nz;i++)
+		M->getData()[zero_inds[i*2+1]*M->getNbRow()+zero_inds[i*2+0]]= FPP(0);
+
+	assert(M->getNonZeros() == nnz);
+
+	vector<pair<int,int>> ref_nz_inds;
+	vector<pair<int,int>> nz_inds = M->nonzeros_indices();
+
+	for(auto indp : nz_inds)
+		cout <<  "(" << indp.first << "," << indp.second << ") ";
+	cout << endl;
+
+
+	bool keep;
+	for(int i=0;i<mat_nrows;i++)
+		for(int j=0;j<mat_ncols;j++)
+		{
+			keep = true;
+			for(int k=0;k<nz && keep;k++)
+				if(i == zero_inds[k*2+0] && j == zero_inds[k*2+1])
+					keep = false;
+			if(keep) ref_nz_inds.push_back(make_pair(i,j));
+		}
+
+	assert(ref_nz_inds.size() == nnz);
+	assert(nz_inds.size() == nnz);
+
+	//assert we find all nonzeros indices and that zero indices are really excluded
+	for(auto ref_pair : ref_nz_inds)
+	{
+//		cout << ref_pair.first << " " << ref_pair.second << endl;
+		assert(find(begin(nz_inds), end(nz_inds), ref_pair) != end(nz_inds));
+	}
+
+	for(int i=0;i<nz;i++)
+	{
+		auto p = make_pair(zero_inds[i*2+0], zero_inds[i*2+1]);
+		assert(find(begin(nz_inds), end(nz_inds), p) == end(nz_inds));
+	}
+
+
+	M->Display();
+
+
+	delete [] zero_inds;
+	return 0;
+}

src/faust_linear_operator/CPU/faust_MatDense.h

@@ -461,6 +461,11 @@ void spgemm(const Faust::MatSparse<FPP,Cpu> & A,const Faust::MatDense<FPP,Cpu> &
 		 */
 		Faust::MatDense<FPP,Cpu> upper_tri(const bool diag=true) const;
 
+		/**
+		 * \brief Returns the nonzeros indices.
+		 */
+		vector<pair<int,int>> nonzeros_indices() const;
+
 		private:
 		Eigen::Matrix<FPP, Eigen::Dynamic, Eigen::Dynamic> mat;
 		bool isIdentity;

src/faust_linear_operator/CPU/faust_MatDense.hpp

@@ -191,6 +191,29 @@ Faust::MatDense<FPP,Cpu> Faust::MatDense<FPP,Cpu>::upper_tri(const bool diag) co
 	return tri;
 }
 
+
+template<typename FPP>
+vector<pair<int,int>> Faust::MatDense<FPP,Cpu>::nonzeros_indices() const
+{
+	vector<pair<int,int>> nz_inds;
+	if(isIdentity)
+		for(int i=0;i<min(this->dim1, this->dim2);i++)
+			nz_inds.push_back(make_pair(i,i));
+	else if(! isZeros)
+	{
+		int i,j;
+		for(int k=0;k<this->dim1*this->dim2;k++)
+			if(mat(k) != FPP(0))
+			{
+				j = k/this->dim1;
+				i = k-j*this->dim1;
+				nz_inds.push_back(make_pair(i,j));
+			}
+	}
+	return nz_inds;
+}
+
+
 	template<typename FPP>
 void Faust::MatDense<FPP,Cpu>::setOnes()
 {