Add the new matrix type MatBSR with its uncompleted test (test_MatBSR).

Only the Frobenius norm is tested in testMatBSR so far.

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 faust_sparse_prox_sp palm4msa_2020 hierarchical2020 hierarchical2020Hadamard hierarchicalFactorizationHadamard hierarchicalFactorizationButterfly)
+	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 hierarchicalFactorizationButterfly test_MatBSR)
 
 	if(FAUST_TORCH)
 		list(APPEND tests faust_torch)

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

+#include "faust_constant.h"
+#include "faust_MatDense.h"
+#include "faust_MatSparse.h"
+#include "faust_MatBSR.h"
+
+
+typedef @TEST_FPP@ FPP;
+
+using namespace std;
+using namespace Faust;
+
+void test_fro_norm(const MatBSR<FPP, Cpu>& bmat)
+{
+	cout << "testing bsr fro. norm:" << std::endl;
+	auto dmat = bmat.to_dense();
+	cout << "bsr mat fro-norm:" << bmat.norm() << std::endl;
+	cout << "dense mat fro-norm:" << dmat.norm() << std::endl;
+	assert(std::abs(bmat.norm()-dmat.norm()) < 1e-6);
+	cout << "fro-norm test is OK." << std::endl;
+}
+
+int main(int argc, char** argv)
+{
+	int m, n, bm, bn, bnnz;
+	m = 10;
+	n = 10;
+	bm = 2;
+	bn = 2;
+	bnnz = 8;
+	auto bmat = MatBSR<FPP, Cpu>::randMat(m, n, bm, bn, bnnz);
+	bmat->Display();
+	test_fro_norm(*bmat);
+	delete bmat;
+	return EXIT_SUCCESS;
+}
+
+

src/faust_linear_operator/CPU/faust_MatBSR.h

+#ifndef __BSR_MAT__
+#define __BSR_MAT__
+#include <Eigen/SparseCore>
+#include "matio.h"
+#include <Eigen/Dense>
+#include <list>
+#include <utility>
+#include "faust_MatGeneric.h"
+
+//template<typename FPP, FDevice DEVICE> class MatGeneric;
+
+template<typename T>
+using DenseMat = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>;
+
+template<typename T>
+using Vec = Eigen::Matrix<T, Eigen::Dynamic, 1>;
+
+template<typename T>
+using DenseMatMap = Eigen::Map<DenseMat<T>>;
+
+template<typename T>
+using SparseMat = Eigen::SparseMatrix<T,Eigen::RowMajor>;
+
+#ifndef FAUST_CONSTANT_H
+template<typename T>
+using Real = typename Eigen::NumTraits<T>::Real;
+#endif
+
+template<typename T,int BlockStorageOrder=0> class BSRMat;
+
+namespace Faust
+{
+	template<typename FPP>
+		class MatBSR<FPP,Cpu> : public MatGeneric<FPP,Cpu>
+		{
+			BSRMat<FPP> bmat; // low-level BSRMat
+			MatBSR() : MatGeneric<FPP, Cpu>() {}
+			MatBSR(BSRMat<FPP>& mat);
+			public:
+			MatGeneric<FPP,Cpu>* Clone(const bool isOptimize=false) const;
+			void multiply(Vect<FPP,Cpu> & vec, char opThis='N') const;
+			Vect<FPP,Cpu> multiply(const Vect<FPP,Cpu> &v) const; // from LinearOperator
+			void multiply(MatDense<FPP,Cpu> & M, char opThis) const;
+			void multiply(MatSparse<FPP, Cpu>& M, char opThis) const;
+			void multiplyRight(MatSparse<FPP, Cpu> const& M) ;
+			void faust_gemm(const MatDense<FPP,Cpu> & B, MatDense<FPP,Cpu> & C,const FPP & alpha, const FPP & beta, char typeA, char typeB)const; // from LinearOperator
+			void transpose();
+			void conjugate(const bool eval=true);
+			void adjoint();
+			faust_unsigned_int getNonZeros()const;
+			size_t getNBytes() const;
+			MatType getType() const;
+			void operator*=(const FPP alpha);
+			std::string to_string(MatType type, const bool transpose=false, const bool displaying_small_mat_elts=false) const;
+			std::string to_string(const bool transpose=false, const bool displaying_small_mat_elts=false) const;
+			matvar_t* toMatIOVar(bool transpose, bool conjugate) const;
+			Real<FPP> normL1(const bool transpose) const;
+			Real<FPP> norm() const;
+			Real<FPP> normL1(faust_unsigned_int& col_id, const bool transpose) const;
+			Vect<FPP,Cpu> get_col(faust_unsigned_int id) const;
+			MatGeneric<FPP,Cpu>* get_cols(faust_unsigned_int col_id_start, faust_unsigned_int num_cols) const;
+			MatGeneric<FPP,Cpu>* get_rows(faust_unsigned_int row_id_start, faust_unsigned_int num_rows) const;
+			MatGeneric<FPP,Cpu>* get_cols(const faust_unsigned_int* col_ids, faust_unsigned_int num_cols) const;
+			MatGeneric<FPP,Cpu>* get_rows(const faust_unsigned_int* row_ids, faust_unsigned_int num_rows) const;
+			std::list<std::pair<int,int>> nonzeros_indices() const;
+			void setZeros();
+			bool containsNaN();
+			const FPP& operator()(faust_unsigned_int i, faust_unsigned_int j)const ;
+			MatDense<FPP, Cpu> to_dense() const;
+			virtual ~MatBSR();
+			static MatBSR<FPP, Cpu>* randMat(int m, int n, int bm, int bn, int bnnz);
+		};
+}
+
+
+//TODO: put all this class and member defs in Faust namespace after unit tests integration
+template<typename T,int BlockStorageOrder>
+class BSRMat
+{
+
+	friend class Faust::MatBSR<T, Cpu>;
+	// \brief the data of nonzero blocks (size bnnz*bm*bn).
+	T* data;
+	// \brief column indices of the nonzero blocks (size bnnz).
+	int* bcolinds;
+	// \brief the buffer of cumulative number of blocks per block row (size b_per_coldim+1).
+	int* browptr;
+	// \brief the number of nonzero blocks.
+	int bnnz;
+	// \brief the number of matrix rows.
+	int m;
+	// \brief the number of matrix columns.
+	int n;
+	// \brief the number of block rows.
+	int bm;
+	// \brief the number of block columns.
+	int bn;
+	// \brief the number of blocks along the row dimension of the matrix (m/bm).
+	int b_per_rowdim;
+	// \brief the number of blocks along the column dimension of the matrix (n/bn).
+	int b_per_coldim;
+
+	// private default constructor
+	BSRMat(): data(nullptr), bcolinds(nullptr), browptr(nullptr), bnnz(0), m(0), n(0), bm(0), bn(0), b_per_rowdim(0), b_per_coldim(0) {}
+	public:
+	/** Copy constructor */
+	BSRMat(const BSRMat<T, BlockStorageOrder>& src_bmat);
+	/** Move constructor */
+	BSRMat(BSRMat<T, BlockStorageOrder>&& src_bmat);
+
+	/** Copy operator */
+	BSRMat<T, BlockStorageOrder>& operator=(const BSRMat<T, BlockStorageOrder>& src_bmat);
+	/** Move operator */
+	BSRMat<T, BlockStorageOrder>& operator=(BSRMat<T, BlockStorageOrder>&& src_bmat);
+
+
+	/**
+	 * Prints internal buffers of this: browptr, bcolinds, data.
+	 */
+	void print_bufs();
+	/**
+	 * Returns true if this contains NaN elements.
+	 */
+	bool contains_nan() const;
+	/**
+	 * Returns the (i,j) entry of this.
+	 */
+	T operator()(unsigned int i, unsigned int j) const;
+	/**
+	 * Returns the indices on the nonzeros of this.
+	 */
+	std::list<std::pair<int,int>> nonzeros_indices() const;
+	/**
+	 * Returns the number of bytes consumed by this.
+	 */
+	size_t nbytes() const;
+	/**
+	 * Returns the nnz of this.
+	 */
+	size_t nnz() const;
+	/**
+	 * Returns the density of this.
+	 */
+	Real<T> density() const;
+	/**
+	 * Returns one column of this as a vector.
+	 *
+	 * \param col_id: the column index.
+	 */
+	Vec<T> get_col(unsigned int col_id) const;
+	/**
+	 * Returns a sequence of num_rows of this into a SparseMat.
+	 *
+	 * \param start_row_id: the first row to return.
+	 * \param num_rows: the number of rows returned.
+	 */
+	SparseMat<T> get_rows(unsigned int start_row_id, unsigned int num_rows) const;
+	/**
+	 * Returns a sequence of num_cols columns of this into a SparseMat.
+	 *
+	 * \param start_col_id: the first column to return.
+	 * \param num_cols: the number of columns returned.
+	 */
+	SparseMat<T> get_cols(unsigned int start_col_id, unsigned int num_cols) const;
+	/**
+	 * Returns a sequence of rows of this into a SparseMat.
+	 *
+	 * \param row_ids: the indices of the rows in the desired order.
+	 * \param num_rows: the size of row_ids.
+	 */
+	SparseMat<T> get_rows(const unsigned long int* row_ids, unsigned int num_rows) const;
+	/**
+	 * Returns a sequence of columns of this into a SparseMat.
+	 *
+	 * \param col_ids: the indices of the columns in the desired order.
+	 * \param num_cols: the size of col_ids.
+	 */
+	SparseMat<T> get_cols(const unsigned long int* col_ids, unsigned int num_cols) const;
+	/**
+	 * Converts this to a DenseMat.
+	 */
+	DenseMat<T> to_dense() const;
+	/**
+	 * Converts this to a SparseMat.
+	 */
+	SparseMat<T> to_sparse() const;
+	/**
+	 * Multiplies this by a DenseMat and returns the result as a DenseMat.
+	 */
+	DenseMat<T> mul(const DenseMat<T>&) const;
+	/**
+	 * Multiplies this by a SparseMat and returns the result as a DenseMat.
+	 */
+	DenseMat<T> mul(const SparseMat<T>&) const;
+	/**
+	 * Multiplies this by a vector (passed directly as a buffer) and returns the result as a Vec.
+	 */
+	Vec<T> mul(const T* vec_data) const;
+	/**
+	 * Multiplies this by a BSRMat and returns the result as DenseMat.
+	 */
+	DenseMat<T> mul(const BSRMat<T, BlockStorageOrder>&) const;
+	/**
+	 * Multiplies this by a scalar (in place).
+	 */
+	void mul(const T& scal) const;
+	/**
+	 * Iterates on the matrix blocks, row by row.
+	 * For each block the argument function is called with the row and column indices of the upper left corner of the block and the block_offset in data and bcolinds.
+	 */
+	void iter_block(std::function<void(int /* mat_row_id */, int /* mat_col_id */, int /* block_offset */)> op) const;
+	/**
+	 * Frobenius norm.
+	 */
+	Real<T> norm() const;
+
+	/**
+	 * Inf-norm.
+	 */
+	Real<T> normInf() const;
+
+	/**
+	 * 1-norm
+	 */
+	Real<T> normL1() const;
+	/**
+	 * Transpose
+	 */
+	BSRMat<T, BlockStorageOrder> transpose(const bool inplace=false);
+	BSRMat<T, BlockStorageOrder> conjugate(const bool inplace=false);
+	/**
+	 * \param m: matrix number of rows.
+	 * \param n: matrix number of columns.
+	 * \param bm: number of rows of blocks (must divide m evenly).
+	 * \param bn: number of columns of blocks (must divide n evenly).
+	 * \param bnnz: number of blocks that contain nonzeros among the m*n/bm/bn blocks.
+	 * \param min_bdensity: the inner minimal density of each block.
+	 * \param max_bdensity: the inner maximal density of each block.
+	 */
+	static BSRMat<T, BlockStorageOrder> rand(int m, int n, int bm, int bn, int bnnz);
+
+	void free_bufs();
+	~BSRMat();
+};
+#include "faust_MatBSR.hpp"
+#endif

src/faust_linear_operator/CPU/faust_MatDense.h

@@ -131,11 +131,15 @@ namespace Faust
 	template<typename FPP>
 		class MatDiag;
 
+	template<typename FPP,FDevice DEVICE>
+		class MatBSR;
+
 	template<typename FPP>
 		class MatDense<FPP,Cpu> : public MatGeneric<FPP,Cpu>
 		{
 
 			friend class MatSparse<FPP,Cpu>;
+			friend class MatBSR<FPP, Cpu>;
 			friend TransformHelper<FPP,Cpu>; // TODO: limit to needed member functions only
 			friend Transform<FPP,Cpu>; //TODO: limit to needed member functions only (multiply)
 			friend void  MatDiag<FPP>::multiply(MatDense<FPP,Cpu> & M, char opThis) const;
@@ -156,8 +160,8 @@ namespace Faust
 																										 \tparam nbRow : number of row of the matrix
 																										 \tparam nbCol : number of column of the matrix
 																										 */
-			MatDense(const FPP  *data_,const faust_unsigned_int nbRow, const faust_unsigned_int nbCol );
-			MatDense(const faust_unsigned_int nbRow, const faust_unsigned_int nbCol, const FPP  *data_);
+			explicit MatDense(const FPP  *data_,const faust_unsigned_int nbRow, const faust_unsigned_int nbCol );
+			explicit MatDense(const faust_unsigned_int nbRow, const faust_unsigned_int nbCol, const FPP  *data_);
 			MatDense() : MatGeneric<FPP,Cpu>(), mat(0,0), isZeros(false) {}
 			/*!
 			 *  \brief Copy Constructor of MatDense
@@ -170,8 +174,8 @@ namespace Faust
 				MatDense(const MatSparse<FPP1,Cpu> & A)  {this->operator=(A);this->is_ortho = A.is_ortho;}
 			MatDense(const MatSparse<FPP,Cpu> & A)  {this->operator=(A);this->is_ortho = A.is_ortho;}
 
-			MatDense(const faust_unsigned_int nbRow, const faust_unsigned_int nbCol) : MatGeneric<FPP,Cpu>(nbRow,nbCol), mat(nbRow,nbCol), isZeros(false){}
-			MatDense(const faust_unsigned_int nbRow) : MatGeneric<FPP,Cpu>(nbRow,nbRow), mat(nbRow,nbRow), isZeros(false){}
+			explicit MatDense(const faust_unsigned_int nbRow, const faust_unsigned_int nbCol) : MatGeneric<FPP,Cpu>(nbRow,nbCol), mat(nbRow,nbCol), isZeros(false){}
+			explicit MatDense(const faust_unsigned_int nbRow) : MatGeneric<FPP,Cpu>(nbRow,nbRow), mat(nbRow,nbRow), isZeros(false){}
 			/// Destructor of MatDense
 			~MatDense(){resize(0,0);/*std::cout<<"destructor dense mat"<<std::endl;*/}
 

src/faust_linear_operator/CPU/faust_MatSparse.h

@@ -91,6 +91,7 @@ namespace Faust
 
 
 	template<typename FPP> class MatDiag;
+	template<typename FPP, FDevice DEVICE> class MatBSR;
 
 	//! MatDense class template of dense matrix
 	template<typename FPP,FDevice DEVICE> class MatDense;
@@ -108,6 +109,7 @@ namespace Faust
 		class MatSparse<FPP,Cpu> : public MatGeneric<FPP,Cpu>
 		{
 
+			friend MatBSR<FPP, Cpu>;
 			friend GivensFGFT<FPP,Cpu, double>;
 			friend GivensFGFTParallel<FPP,Cpu, double>;
 			friend GivensFGFT<FPP,Cpu, float>;

src/faust_linear_operator/CPU/faust_Vect.h

@@ -91,6 +91,9 @@ namespace Faust
 	template<typename FPP>
 		class MatDiag;
 
+	template<typename FPP,FDevice DEVICE>
+		class MatBSR;
+
    // friend function of faust_linear_algebra.h
    template<typename FPP>
    void gemv(const MatDense<FPP,Cpu> & A,const Vect<FPP,Cpu> & x,Vect<FPP,Cpu> & y,const FPP & alpha, const FPP & beta, char typeA);
@@ -112,6 +115,7 @@ namespace Faust
 		friend class MatSparse<FPP,Cpu>; // TODO: limit friendship to the necessary functions
 //		friend void  MatDiag<FPP>::multiply(Vect<FPP,Cpu> & vec, char opThis) const;
 
+		friend class MatBSR<FPP, Cpu>;
 		friend TransformHelperPoly<FPP>; // TODO: limit to needed member functions only
 
         public :

src/faust_linear_operator/faust_MatGeneric.h

@@ -66,6 +66,9 @@ namespace Faust
     template<typename FPP,FDevice DEVICE>
     class LinearOperator;
 
+    template<typename FPP,FDevice DEVICE>
+    class MatBSR;
+
     template<typename FPP,FDevice DEVICE>
     class MatGeneric : public Faust::LinearOperator<FPP,DEVICE>
 	{
@@ -146,7 +149,7 @@ namespace Faust
 			//! \return value between 0 and 1
 			float density() const{return ((float) this->getNonZeros())/((float)this->getNbCol()*this->getNbRow());}
 
-			//! \brief get the dynamic type of the matrix (SPARSE or DENSE)
+			//! \brief get the dynamic type of the matrix (SPARSE, DENSE, etc.)
 			virtual MatType getType() const=0;
 
 			//! \brief multiply a matrix by the given scalar alpha
@@ -217,7 +220,7 @@ namespace Faust
 
 			//! \brief
 			//! \warning : declare a virtual destructor is mandatory for an abstract class
-			//! in order to allow descendant class destructor to clean up in case of pointer to the abstract class
+			//! in order to allow child class destructor to clean up in case of pointer to the abstract class
 			virtual ~ MatGeneric()=0;
 
 

src/faust_linear_operator/faust_MatGeneric.hpp

@@ -42,6 +42,7 @@
 
 // useful for optimize with multiplication
 #include "faust_Timer.h"
+#include "faust_MatBSR.h"
 
 
 template<typename FPP,FDevice DEVICE>
@@ -148,6 +149,8 @@ std::string Faust::MatGeneric<FPP,DEVICE>::to_string(const bool transpose /* set
 		type = Dense;
 	else if (dynamic_cast<const MatDiag<FPP>*>(this))
 		type = Diag;
+	else if (dynamic_cast<const MatBSR<FPP, DEVICE>*>(this))
+		type = BSR;
 	else
 		throw std::runtime_error("Unhandled matrix type in MatGeneric::to_string()"); // shouldn't happen
 	return this->to_string(getNbRow(), getNbCol(), transpose, this->density(), this->getNonZeros(), this->is_identity, type);
@@ -181,6 +184,8 @@ std::string Faust::MatGeneric<FPP,DEVICE>::to_string(int32_t nrows, int32_t ncol
 		str << "SPARSE,";
 	else if(type == Diag)
 		str << "DIAG,";
+	else if(type == BSR)
+		str << "BSR,";
 	else if(type == None)
 		str << "UNKNOWN MATRIX TYPE,";
 	else

src/faust_linear_operator/faust_constant.h.in

@@ -93,6 +93,7 @@ enum MatType
    Dense, 
    Sparse, 
    Diag,
+   BSR,
    None
 };