Add MatSparse::vstack and MatSparse::hstack along with their unit tests.

misc/test/CMakeLists.txt

@@ -150,7 +150,7 @@ endif()
 
 if(NOT NOCPPTESTS)
 	foreach(TEST_FPP float double)
-		foreach(FILE faust_mult2 faust_mult faust_mult_cplx test_Vect_min test_MatDense_get_row test_MatDense_lower_upper_tri test_MatDense_nonzeros_indices test_Transform_move test_TransformHelper_and_Transform_copy_ctor test_TransformHelper_and_Transform_fac_iterato test_TransformHelper_variadic_template_ctor test_MatDense_min test_TH_pack_factors faust_transform_omp_mul faust_pruneout faust_transform_optimize_storage faust_transform_optimize faust_prox_blockdiag)
+		foreach(FILE faust_mult2 faust_mult faust_mult_cplx test_Vect_min test_MatDense_get_row test_MatDense_lower_upper_tri test_MatDense_nonzeros_indices test_Transform_move test_TransformHelper_and_Transform_copy_ctor test_TransformHelper_and_Transform_fac_iterato test_TransformHelper_variadic_template_ctor test_MatDense_min test_TH_pack_factors faust_transform_omp_mul faust_pruneout faust_transform_optimize_storage faust_transform_optimize faust_prox_blockdiag test_MatSparse_vstack test_MatSparse_hstack)
 			if(USE_GPU_MOD AND ${FILE} MATCHES test_TH_pack_factors|test_TransformHelper_variadic_template_ctor|faust_transform_optimize|test_TransformHelper_and_Transform_fac_iterato|faust_mult2|faust_transform_omp_mul|test_Transform_move|test_TransformHelper_and_Transform_copy_ctor AND ${TEST_FPP} MATCHES float)# gpu_mod handles float but not FaustGPU
 				message(STATUS ${FILE}_${TEST_FPP} " skipped")
 				continue()

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

+#include "faust_MatSparse.h"
+
+#include <complex>
+#include <random>
+#include <iostream>
+
+using namespace Faust;
+using namespace std;
+typedef @TEST_FPP@ FPP;
+
+
+int main()
+{
+
+	MatSparse<FPP,Cpu>* M1 = MatSparse<FPP,Cpu>::randMat(4, 3, .5);
+	MatSparse<FPP,Cpu>* M2 = MatSparse<FPP,Cpu>::randMat(4, 13, .5);
+	MatSparse<FPP,Cpu> M;
+	M.hstack(*M1, *M2);
+	cout << "M1:" << endl;
+	M1->Display();
+	cout << "M2:" << endl;
+	M2->Display();
+	cout << "M:" << endl;
+	M.Display();
+	assert(M.getNbRow() == M1->getNbRow() && M.getNbRow() == M2->getNbRow());
+	assert(M.getNbCol() == M1->getNbCol() + M2->getNbCol());
+	assert(M1->getNonZeros()+M2->getNonZeros() == M.getNonZeros());
+	FPP* buf = new FPP[M.getNonZeros()];
+	memcpy(buf, M1->getValuePtr(), sizeof(FPP)*M1->getNonZeros());
+	memcpy(buf+M1->getNonZeros(), M2->getValuePtr(), sizeof(FPP)*M2->getNonZeros());
+	MatDense<FPP,Cpu> Mds(buf, M.getNonZeros(), 1);
+	assert(abs(Mds.norm()-M.norm()) < 1e-1);
+	for(int i=0;i < M.getNonZeros(); i++)
+//		cout << M.getValuePtr()[i] << " " << buf[i] << endl;
+		assert(abs(M.getValuePtr()[i] - buf[i] < 1e-6));
+	delete [] buf;
+	return 0;
+}

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

+#include "faust_MatSparse.h"
+
+#include <complex>
+#include <random>
+#include <iostream>
+
+using namespace Faust;
+using namespace std;
+typedef @TEST_FPP@ FPP;
+
+
+int main()
+{
+
+	MatSparse<FPP,Cpu>* M1 = MatSparse<FPP,Cpu>::randMat(3, 4, .5);
+	MatSparse<FPP,Cpu>* M2 = MatSparse<FPP,Cpu>::randMat(13, 4, .5);
+	MatSparse<FPP,Cpu> M;
+	M.vstack(*M1, *M2);
+	cout << "M1:" << endl;
+	M1->Display();
+	cout << "M2:" << endl;
+	M2->Display();
+	cout << "M:" << endl;
+	M.Display();
+	assert(M.getNbCol() == M1->getNbCol() && M.getNbCol() == M2->getNbCol());
+	assert(M.getNbRow() == M1->getNbRow() + M2->getNbRow());
+	assert(M1->getNonZeros()+M2->getNonZeros() == M.getNonZeros());
+	FPP* buf = new FPP[M.getNonZeros()];
+	memcpy(buf, M1->getValuePtr(), sizeof(FPP)*M1->getNonZeros());
+	memcpy(buf+M1->getNonZeros(), M2->getValuePtr(), sizeof(FPP)*M2->getNonZeros());
+	MatDense<FPP,Cpu> Mds(buf, M.getNonZeros(), 1);
+	assert(abs(Mds.norm()-M.norm()) < 1e-1);
+	for(int i=0;i < M.getNonZeros(); i++)
+//		cout << M.getValuePtr()[i] << " " << buf[i] << endl;
+		assert(abs(M.getValuePtr()[i] - buf[i] < 1e-6));
+	delete [] buf;
+	return 0;
+}

src/faust_linear_operator/CPU/faust_MatSparse.h

@@ -347,6 +347,16 @@ namespace Faust
 
 			std::list<std::pair<int,int>> nonzeros_indices() const;
 
+			/**
+			 * \brief Concatenates vertically top and bottom matrices, resizing this if necessary.
+			 */
+			void vstack(MatSparse<FPP, Cpu>& top, MatSparse<FPP, Cpu>& bottom);
+
+			/**
+			 * \brief Concatenates horizontally left and right matrices, resizing this if necessary.
+			 */
+			void hstack(MatSparse<FPP, Cpu>& left, MatSparse<FPP, Cpu>& right);
+
 			static MatSparse<FPP, Cpu>* randMat(faust_unsigned_int num_rows, faust_unsigned_int num_cols, double density);
 			//\param : per_row means the density applies for each line rather than globally for the matrix
 			static MatSparse<FPP, Cpu>* randMat(faust_unsigned_int num_rows, faust_unsigned_int num_cols, double density, bool per_row);

src/faust_linear_operator/CPU/faust_MatSparse.hpp

@@ -1213,5 +1213,82 @@ void Faust::MatSparse<FPP,Cpu>::copyBufs(U* out_rowptr, U* out_colind, FPP* out_
 	copyColInd(out_colind);
 }
 
+template<typename FPP>
+void Faust::MatSparse<FPP,Cpu>::vstack(MatSparse<FPP, Cpu>& top, MatSparse<FPP, Cpu>& bottom)
+{
+	auto tncols = top.getNbCol();
+	auto bncols = bottom.getNbCol();
+	auto tnrows = top.getNbRow();
+	auto bnrows = bottom.getNbRow();
+	auto tnnz = top.getNonZeros();
+	auto bnnz = bottom.getNonZeros();
+	auto nrows = tnrows+bnrows;
+	auto nnz = tnnz + bnnz;
+	if(tncols != bncols)
+		throw std::runtime_error("vstack error: dimensions must agree.");
+	if(this->getNbCol() != tncols || this->getNbRow() != nrows || this->getNonZeros() != nnz)
+		resize(nnz, nrows, tncols);
+	// copy column indices
+	memcpy(getColInd(), top.getColInd(), sizeof(int)*tnnz);
+	memcpy(getColInd()+tnnz, bottom.getColInd(), sizeof(int)*bnnz);
+	// copy values
+	memcpy(getValuePtr(), top.getValuePtr(), sizeof(FPP)*tnnz);
+	memcpy(getValuePtr()+tnnz, bottom.getValuePtr(), sizeof(FPP)*bnnz);
+	// build rowptr
+	memcpy(getRowPtr(), top.getRowPtr(), sizeof(FPP)*tnrows+1);
+	int *rowptr, *browptr, row_offset = *(top.getRowPtr()+tnrows);
+	int j = tnrows;
+	//TODO: openmp
+//#pragma omp parallel for
+	for(int i=1;i < bnrows; i++)
+	{
+		rowptr = getRowPtr()+tnrows+i;
+		browptr = bottom.getRowPtr()+i;
+		*rowptr = *browptr+row_offset;
+	}
+}
+
+template<typename FPP>
+void Faust::MatSparse<FPP,Cpu>::hstack(MatSparse<FPP, Cpu>& left, MatSparse<FPP, Cpu>& right)
+{
+	auto lncols = left.getNbCol();
+	auto rncols = right.getNbCol();
+	auto lnrows = left.getNbRow();
+	auto rnrows = right.getNbRow();
+	auto lnnz = left.getNonZeros();
+	auto rnnz = right.getNonZeros();
+	auto nrows = lnrows;
+	auto ncols = lncols + rncols;
+	auto nnz = lnnz + rnnz;
+	int *rowptr = getRowPtr(), *lrowptr = left.getRowPtr(), *rrowptr = right.getRowPtr();
+	if(lnrows != rnrows)
+		throw std::runtime_error("vstack error: dimensions must agree.");
+	if(this->getNbCol() != lncols || this->getNbRow() != nrows || this->getNonZeros() != nnz)
+		resize(nnz, nrows, ncols);
+	// copy column indices
+	memcpy(getColInd(), left.getColInd(), sizeof(int)*lnnz);
+	int *colind = getColInd()+lnnz;
+	//TODO: openmp
+//#pragma omp parallel for
+	for(int i=0;i<rnnz;i++)
+	{
+		*colind = *(right.getColInd()+i)+lncols;
+		colind++;
+	}
+	// copy values
+	memcpy(getValuePtr(), left.getValuePtr(), sizeof(FPP)*lnnz);
+	memcpy(getValuePtr()+lnnz, right.getValuePtr(), sizeof(FPP)*rnnz);
+	// build rowptr
+	//TODO: openmp
+//#pragma omp parallel for
+	for(int i=0;i <= nrows; i++)
+	{
+		*rowptr = *lrowptr + *rrowptr;
+		rowptr++;
+		lrowptr++;
+		rrowptr++;
+	}
+}
+
 
 #endif