New prototype of svdtj in pyfaust wrapper for supporting (sparse) csr_matrix.

wrapper/python/pyfaust/fact.py

@@ -115,7 +115,10 @@ def svdtj(M, maxiter, tol=0, relerr=True,  nGivens_per_fac=None, verbosity=0):
      See also:
         eigtj
     """
-    Ucore, S, Vcore =  _FaustCorePy.FaustFact.svdtj(M, maxiter, nGivens_per_fac, verbosity, tol, relerr)
+    if(isinstance(M, np.ndarray)):
+        Ucore, S, Vcore =  _FaustCorePy.FaustFact.svdtj(M, maxiter, nGivens_per_fac, verbosity, tol, relerr)
+    elif(isinstance(M, csr_matrix)):
+        Ucore, S, Vcore =  _FaustCorePy.FaustFact.svdtj_sparse(M, maxiter, nGivens_per_fac, verbosity, tol, relerr)
     U = Faust(core_obj=Ucore)
     V = Faust(core_obj=Vcore)
     return U, S, V

wrapper/python/src/FaustCoreCy.pxd

@@ -226,3 +226,9 @@ cdef extern from "FaustFactGivensFGFT.h":
                          num_cols, unsigned int J, unsigned int t, unsigned int
                          verbosity, const double stoppingError, const bool errIsRel)
 
+    cdef void svdtj_sparse[FPP, FPP2](FaustCoreCpp[FPP]** U, FaustCoreCpp[FPP] **V, FPP* S,
+                                      const FPP* data, int* row_ptr, int* id_col, int
+                                      nnz, int nrows, int ncols, unsigned int J,
+                                      unsigned int t, unsigned int verbosity,
+                                      const double stoppingError, const bool errIsRel)
+

wrapper/python/src/FaustFactGivensFGFT.h

@@ -27,6 +27,9 @@ FaustCoreCpp<FPP>* fact_givens_fgft_generic_cplx(GivensFGFTComplex<FPP, Cpu, FPP
 template<typename FPP, typename FPP2 = float>
 void svdtj( FaustCoreCpp<FPP>** U, FaustCoreCpp<FPP> **V, FPP* S /** end of output parameters */, const FPP* M, unsigned int num_rows, unsigned int num_cols, unsigned int J, unsigned int t, unsigned int verbosity = 0, const double stoppingError = 0.0, const bool errIsRel = true /* end of input parameters*/);
 
+template<typename FPP, typename FPP2 = float>
+void svdtj_sparse(FaustCoreCpp<FPP>** U, FaustCoreCpp<FPP> **V, FPP* S, /*start of input parameters*/ const FPP* data, int* row_ptr, int* id_col, int nnz, int nrows, int ncols, unsigned int J, unsigned int t, unsigned int verbosity = 0, const double stoppingError = 0.0, const bool errIsRel = true);
+
 #include "FaustFactGivensFGFT.hpp"
 
 #endif

wrapper/python/src/FaustFactGivensFGFT.hpp

@@ -116,18 +116,30 @@ FaustCoreCpp<FPP>* fact_givens_fgft_generic_cplx(GivensFGFTComplex<FPP, Cpu, FPP
     return new FaustCoreCpp<FPP>(th);
 }
 
-template<typename FPP, typename FPP2 = float>
+template<typename FPP, typename FPP2>
 void svdtj(FaustCoreCpp<FPP>** U, FaustCoreCpp<FPP> **V, FPP* S, /*start of input parameters*/ const FPP* M_data, unsigned int num_rows, unsigned int num_cols, unsigned int J, unsigned int t, unsigned int verbosity, const double stoppingError, const bool errIsRel)
 {
     Faust::MatDense<FPP,Cpu> M(M_data, (faust_unsigned_int) num_rows, (faust_unsigned_int) num_cols);
-//MatDense<FPP, Cpu> & dM, int J, int t, double tol, unsigned int verbosity, bool relErr, int order, TransformHelper<FPP,Cpu> ** U, TransformHelper<FPP,Cpu> **V, Faust::Vect<FPP,Cpu> ** S_
     TransformHelper<FPP,Cpu> *U_,  *V_;
     Faust::Vect<FPP,Cpu> * S_;
     svdtj(M, J, t, stoppingError, verbosity, errIsRel, 1 /* order (useless) */, &U_, &V_, &S_);
     *U = new FaustCoreCpp<FPP>(U_);
     *V = new FaustCoreCpp<FPP>(V_);
     //TODO: avoid this copy by directly edit outside buffer S
-//    delete S;
+    memcpy(S, S_->getData(), sizeof(FPP)* S_->size());
+    delete S_;
+}
+
+template<typename FPP, typename FPP2>
+void svdtj_sparse(FaustCoreCpp<FPP>** U, FaustCoreCpp<FPP> **V, FPP* S, /*start of input parameters*/ const FPP* data, int* row_ptr, int* id_col, int nnz, int nrows, int ncols, unsigned int J, unsigned int t, unsigned int verbosity, const double stoppingError, const bool errIsRel)
+{
+    Faust::MatSparse<FPP, Cpu> M(nnz, nrows, ncols, data, id_col, row_ptr);
+    TransformHelper<FPP,Cpu> *U_,  *V_;
+    Faust::Vect<FPP,Cpu> * S_;
+    svdtj(M, J, t, stoppingError, verbosity, errIsRel, 1 /* order (useless) */, &U_, &V_, &S_);
+    *U = new FaustCoreCpp<FPP>(U_);
+    *V = new FaustCoreCpp<FPP>(V_);
+    //TODO factorize with svdtj
     memcpy(S, S_->getData(), sizeof(FPP)* S_->size());
     delete S_;
 }

wrapper/python/src/_FaustCorePy.pyx

@@ -1592,3 +1592,44 @@ cdef class FaustFact:
         #return core, S_spdiag
         return coreU, S, coreV
 
+    @staticmethod
+    def svdtj_sparse(M, J, t, verbosity=0, stoppingError = 0.0,
+                     errIsRel=True):
+        from scipy.sparse import spdiags
+        cdef double [:] data1d #only for csr mat factor
+        cdef int [:] indices # only for csr mat
+        cdef int [:] indptr # only for csr mat
+        cdef unsigned int M_num_rows=M.shape[0]
+        cdef unsigned int M_num_cols=M.shape[1]
+        cdef double[:] S_view
+
+
+        data1d = M.data.astype(float,'F')
+        indices = M.indices.astype(np.int32, 'F')
+        indptr = M.indptr.astype(np.int32, 'F')
+
+        S = np.empty(M.shape[0], dtype=M.dtype)
+        S_view = S
+
+        coreU = FaustCore(core=True)
+        coreV = FaustCore(core=True)
+        FaustCoreCy.svdtj_sparse[double,double](
+            &(coreU.core_faust_dbl),
+            &(coreV.core_faust_dbl),
+            &S_view[0],
+            &data1d[0],
+            &indices[0],
+            &indptr[0],
+            M.nnz,
+            M_num_rows,
+            M_num_cols, J, t,
+            verbosity,
+            stoppingError,
+            errIsRel)
+
+        coreU._isReal = coreV._isReal = True
+        #D_spdiag = spdiags(D, [0], M.shape[0], M.shape[0])
+        #return core, D_spdiag
+        return coreU, S, coreV
+
+