Faust-sparse matrix multiplication available also for matfaust.

0496252e added it for pyfaust (same C++ backend).
Related to issue #24.

src/faust_linear_operator/CPU/faust_TransformHelper.h

@@ -92,9 +92,9 @@ TransformHelper(Transform<FPP,Cpu> &t, const bool moving=false);
 
 			Vect<FPP,Cpu> multiply(const Vect<FPP,Cpu> x) const;
 			Vect<FPP,Cpu> multiply(const Vect<FPP,Cpu> x, const bool transpose);
-			MatDense<FPP,Cpu> multiply(const MatDense<FPP,Cpu> A) const;
-			MatDense<FPP, Cpu> multiply(const MatDense<FPP,Cpu> A, const bool transpose);
-			MatSparse<FPP, Cpu> multiply(const MatSparse<FPP,Cpu> A, const bool transpose=false) const;
+//			MatDense<FPP,Cpu> multiply(const MatDense<FPP,Cpu> A) const;
+			MatDense<FPP, Cpu> multiply(const MatDense<FPP,Cpu> A, const bool transpose=false);
+			MatDense<FPP, Cpu> multiply(const MatSparse<FPP,Cpu> A, const bool transpose=false);
 
 			TransformHelper<FPP, Cpu>* multiply(TransformHelper<FPP, Cpu>*);
 			TransformHelper<FPP, Cpu>* multiply(FPP& scalar);

src/faust_linear_operator/CPU/faust_TransformHelper.hpp

@@ -160,18 +160,11 @@ namespace Faust {
 		}
 
 	template<typename FPP>
-		MatDense<FPP,Cpu> TransformHelper<FPP,Cpu>::multiply(const MatDense<FPP,Cpu> A) const
+		MatDense<FPP,Cpu> TransformHelper<FPP,Cpu>::multiply(const MatSparse<FPP,Cpu> A, const bool transpose /* deft to false */)
 		{
+			is_transposed ^= transpose;
 			MatDense<FPP,Cpu> M = this->transform->multiply(A, isTransposed2char());
-			if(is_conjugate) M.conjugate();
-			return M;
-		}
-
-	template<typename FPP>
-		MatSparse<FPP,Cpu> TransformHelper<FPP,Cpu>::multiply(const MatSparse<FPP,Cpu> A, const bool transpose /* deft to false */) const
-		{
-			MatSparse<FPP,Cpu> M = this->transform->multiply(A, isTransposed2char());
-			if(is_conjugate) M.conjugate();
+			is_transposed ^= transpose;
 			return M;
 		}
 

wrapper/matlab/+matfaust/@Faust/Faust.m

@@ -376,7 +376,7 @@ classdef Faust
 		%> @b Usage
 		%>
 		%> &nbsp;&nbsp;&nbsp; @b G = mtimes(F, A)<br/>
-		%> &nbsp;&nbsp;&nbsp; @b G = F*A with A and G being full matrices.<br/>
+		%> &nbsp;&nbsp;&nbsp; @b G = F*A with A and G being full matrices. A could also be a sparse matrix but note that often the Faust-sparse matrix multiplication is slower than performing F*full(A) multiplication. In some cases though, it stays quicker: moreover when the Faust is composed of a small number of factors.<br/>
 		%> &nbsp;&nbsp;&nbsp; @b G = F*s with s a scalar and G a Faust.<br/>
 		%> &nbsp;&nbsp;&nbsp; @b G = s*F with s a scalar and G a Faust.<br/>
 		%> &nbsp;&nbsp;&nbsp; @b G = s*F' with s a scalar multiplying the conjugate-transpose of F.<br/>
@@ -409,21 +409,6 @@ classdef Faust
 		%> @endcode
 		%>
 		%> @b Errors
-		%> - The multiplying operand A is a sparse matrix:
-		%>
-		%> @code
-		%>>> issparse(S)
-		%>
-		%>ans =
-		%>
-		%>  logical
-		%>
-		%>     1
-		%>
-		%>>> F*S
-		%> Error using matfaust
-		%> Faust multiplication to a sparse matrix isn't supported.
-		%> @endcode
 		%>
 		%> - F is real but A is a complex scalar.
 		%>
@@ -502,9 +487,10 @@ classdef Faust
 			% TODO: take trans into account when mul F to a scal or a Faust
 			% it's not a serious issue because mtimes_trans() shouln't be called by final user
 			% (func's doc is filtered out in doxydoc)
-			if(issparse(A))
-				error('Faust multiplication to a sparse matrix isn''t supported.')
-			elseif(isa(A,'matfaust.Faust'))
+			%if(issparse(A))
+			%	error('Faust multiplication to a sparse matrix isn''t supported.')
+			%elseif(isa(A,'matfaust.Faust'))
+			if(isa(A,'matfaust.Faust'))
 				if (F.isReal)
 					if(isreal(A))
 						C = matfaust.Faust(F, mexFaustReal('mul_faust', F.matrix.objectHandle, A.matrix.objectHandle));

wrapper/matlab/src/mexFaust.cpp.in

@@ -592,57 +592,65 @@ void save(Faust::TransformHelper<SCALAR,Cpu>* core_ptr, int nargs, const mxArray
                 // input matrix or vector from MATLAB
                 const mxArray * inMatlabMatrix = prhs[2];
 
+				if(mxIsSparse(inMatlabMatrix))
+				{
+					Faust::MatSparse<SCALAR,Cpu> spA;
+					mxArray2FaustspMat<SCALAR>(inMatlabMatrix, spA);
+					Faust::MatDense<SCALAR,Cpu> B;
+					B = (*core_ptr).multiply(spA, transpose_flag);
+					plhs[0] = FaustMat2mxArray(B);
+				}
+				else
+				{
+					const size_t nbRowA = mxGetM(inMatlabMatrix);
+					const size_t nbColA = mxGetN(inMatlabMatrix);
+					faust_unsigned_int nbRowOp_,nbColOp_;
+					const size_t nbRowOp = core_ptr->getNbRow();
+					const size_t nbColOp = core_ptr->getNbCol();
+					const size_t nbRowB = nbRowOp;
+					const size_t nbColB = nbColA;
 
 
-                const size_t nbRowA = mxGetM(inMatlabMatrix);
-                const size_t nbColA = mxGetN(inMatlabMatrix);
-                faust_unsigned_int nbRowOp_,nbColOp_;
-                const size_t nbRowOp = core_ptr->getNbRow();
-                const size_t nbColOp = core_ptr->getNbCol();
-                const size_t nbRowB = nbRowOp;
-                const size_t nbColB = nbColA;
-
-
-                /** Check parameters **/
-
-                //check dimension match
-                if (mxGetNumberOfDimensions(inMatlabMatrix) != 2
-                        || nbRowA != nbColOp && (nbRowA != nbRowOp && transpose_flag))
-                    mexErrMsgTxt("Multiply : Wrong number of dimensions for the input vector or matrix (third argument).");
+					/** Check parameters **/
 
+					//check dimension match
+					if (mxGetNumberOfDimensions(inMatlabMatrix) != 2
+							|| nbRowA != nbColOp && (nbRowA != nbRowOp && transpose_flag))
+						mexErrMsgTxt("Multiply : Wrong number of dimensions for the input vector or matrix (third argument).");
 
-                SCALAR* ptr_data = NULL;
-                if ((core_ptr->isReal() ) && ( mxIsComplex(inMatlabMatrix) ) )
-                    mexErrMsgTxt("impossibility to multiply a real Faust with complex matrix");
 
-                mxArray2Ptr(inMatlabMatrix, ptr_data);
+					SCALAR* ptr_data = NULL;
+					if ((core_ptr->isReal() ) && ( mxIsComplex(inMatlabMatrix) ) )
+						mexErrMsgTxt("impossibility to multiply a real Faust with complex matrix");
 
-                // Si inMatlabMatrix est un vecteur
+					mxArray2Ptr(inMatlabMatrix, ptr_data);
 
-                if(nbColA == 1)
-                {
-                    // applying the Faust to a vector
-                    Faust::Vect<SCALAR,Cpu> A(nbRowA, ptr_data);
-                    Faust::Vect<SCALAR,Cpu> B(nbRowB);
-                    B = (*core_ptr).multiply(A, transpose_flag);
+					// Si inMatlabMatrix est un vecteur
 
+					if(nbColA == 1)
+					{
+						// applying the Faust to a vector
+						Faust::Vect<SCALAR,Cpu> A(nbRowA, ptr_data);
+						Faust::Vect<SCALAR,Cpu> B(nbRowB);
+						B = (*core_ptr).multiply(A, transpose_flag);
 
-                    plhs[0]=FaustVec2mxArray(B);
-                }
-                else
-                { // applying the Faust to a matrix
-                    Faust::MatDense<SCALAR,Cpu> A(ptr_data, nbRowA, nbColA);
-                    Faust::MatDense<SCALAR,Cpu> B(nbRowB, nbColA);
-                    B = (*core_ptr).multiply(A, transpose_flag);
 
-                    plhs[0]=FaustMat2mxArray(B);
+						plhs[0]=FaustVec2mxArray(B);
+					}
+					else
+					{ // applying the Faust to a matrix
+						Faust::MatDense<SCALAR,Cpu> A(ptr_data, nbRowA, nbColA);
+						Faust::MatDense<SCALAR,Cpu> B(nbRowB, nbColA);
+						B = (*core_ptr).multiply(A, transpose_flag);
 
-                }
+						plhs[0]=FaustMat2mxArray(B);
 
+					}
 
 
-                if(ptr_data) {delete [] ptr_data ; ptr_data = NULL;}
 
+					if(ptr_data) {delete [] ptr_data ; ptr_data = NULL;}
+				}
                 return;
             }
 

wrapper/python/pyfaust.py

@@ -551,10 +551,16 @@ class Faust:
 
         Args:
             F: the Faust object.
-            A: a Faust object or a 2D full matrix (numpy.ndarray, numpy.matrix).
+            A: a Faust object a scipy.sparse.csr_matrix or a 2D full matrix (numpy.ndarray, numpy.matrix).
             <br/> In the latter case, A must be Fortran contiguous (i.e. Column-major order;
                 `order' argument passed to np.ndararray() must be equal to str
                 'F').
+            <br/> Note that performing a Faust-csr_matrix product is often
+            slower than converting first the csr_matrix to a dense
+            representation (toarray()) and then proceed to the
+            Faust-dense_matrix multiplication. In some cases though, it stays
+            quicker: moreover when the Faust is composed of a small number of
+            factors.
 
         Returns: The result of the multiplication
             - as a numpy.ndarray if A is a ndarray,<br/>
@@ -861,9 +867,8 @@ class Faust:
 
         <b/>See also Faust.todense
         """
-        identity = np.eye(F.shape[1], F.shape[1])
-        F_dense = F*identity
-        return F_dense
+        return F.m_faust.get_product()
+
 
     def todense(F):
         """

wrapper/python/src/FaustCoreCpp.h

@@ -63,6 +63,7 @@ class FaustCoreCpp
     void push_back(FPP* data, int* row_ptr, int* id_col, int nnz, int nrows, int ncols, bool optimizedCopy=false);
     unsigned int getNbRow() const;
     unsigned int getNbCol() const;
+    void get_product(FPP* y_data, int y_nrows, int y_ncols);
     void multiply(FPP* y_data, int y_nrows, int y_ncols, FPP* x_data, int* x_row_ptr, int* x_id_col, int x_nnz, int x_nrows, int x_ncols);
     void multiply(FPP* value_y,int nbrow_y,int nbcol_y,FPP* value_x,int nbrow_x,int nbcol_x/*,bool isTranspose*/)const;
     FaustCoreCpp<FPP>* mul_faust(FaustCoreCpp<FPP>* right);

wrapper/python/src/FaustCoreCpp.hpp

@@ -108,6 +108,12 @@ FaustCoreCpp<FPP>* FaustCoreCpp<FPP>::mul_scal(FPP scal)
     FaustCoreCpp<FPP>* core = new FaustCoreCpp<FPP>(th);
     return core;
 }
+template<typename FPP>
+void FaustCoreCpp<FPP>::get_product(FPP* y_data, int y_nrows, int y_ncols)
+{
+    Faust::MatDense<FPP, Cpu> Y = this->transform->get_product();
+    memcpy(y_data, Y.getData(), sizeof(FPP)*y_ncols*y_nrows);
+}
 
 template<typename FPP>
 void FaustCoreCpp<FPP>::multiply(FPP* y_data, int y_nrows, int y_ncols, FPP* x_data, int* x_row_ptr, int* x_id_col, int x_nnz, int x_nrows, int x_ncols)

wrapper/python/src/FaustCoreCy.pxd

@@ -48,6 +48,7 @@ cdef extern from "FaustCoreCpp.h" :
                        bool optimizedCopy)
         void push_back(FPP* data, int* row_ptr, int* id_col, int nnz, int nrows,
                        int ncols, bool optimizedCopy)
+        void get_product(FPP* data, int nbrow, int nbcol);
         void multiply(FPP* value_y,int nbrow_y,int nbcol_y,FPP* value_x,
                       int nbrow_x, int nbcol_x);#,bool isTranspose*/);
         # Faust-by-csr product -> dense mat

wrapper/python/src/_FaustCorePy.pyx

@@ -245,7 +245,24 @@ cdef class FaustCore:
             raise Exception("complex Faust-csr_matrix mul not yet supported")
         return y_data_arr
 
+    def get_product(self):
+        cdef double [:,:] y_data
+        cdef complex [:,:] y_data_cplx
 
+        if(self._isReal):
+            y_arr = np.empty((self.core_faust_dbl.getNbRow(), self.core_faust_dbl.getNbCol()), dtype='d',
+                             order='F')
+            y_data = y_arr
+            self.core_faust_dbl.get_product(&y_data[0,0], y_arr.shape[0],
+                                            y_arr.shape[1])
+        else:
+            y_arr = np.empty((self.core_faust_cplx.getNbRow(), self.core_faust_cplx.getNbCol()),
+                             dtype=np.complex,
+                             order='F')
+            y_data_cplx = y_arr
+            self.core_faust_cplx.get_product(&y_data_cplx[0,0], y_arr.shape[0],
+                                            y_arr.shape[1])
+        return y_arr
 
     cdef multiply_faust(self, F):
         if(isinstance(F, FaustCore)):