hmatrix-svdlibc 0.3.2 → 0.4.0
raw patch · 7 files changed
+227/−24 lines, 7 filesdep +QuickCheckdep +criteriondep +hmatrix-svdlibcdep ~hmatrix
Dependencies added: QuickCheck, criterion, hmatrix-svdlibc, hspec, vector
Dependency ranges changed: hmatrix
Files
- Numeric/LinearAlgebra/SVD/SVDLIBC.hs +60/−7
- Test.hs +0/−9
- cbits/glue.c +44/−0
- hmatrix-svdlibc.cabal +34/−8
- include/glue.h +10/−0
- test/Numeric/LinearAlgebra/SVD/Benchmarks.hs +19/−0
- test/Numeric/LinearAlgebra/SVD/Spec.hs +60/−0
Numeric/LinearAlgebra/SVD/SVDLIBC.hs view
@@ -1,15 +1,17 @@ {-# LANGUAGE ForeignFunctionInterface #-} module Numeric.LinearAlgebra.SVD.SVDLIBC- (svd) where+ (svd, sparsify, sparseSvd) where import Control.Applicative-import qualified Data.Packed as P-import qualified Data.Packed.Development as I+import qualified Numeric.LinearAlgebra.Data as P+import qualified Numeric.LinearAlgebra.Devel as I+import qualified Data.Vector.Storable as Vec import Foreign hiding (unsafePerformIO) import Foreign.C.Types import System.IO.Unsafe import Foreign.Marshal.Alloc+import System.IO newtype DenseMatrix = DMat (ForeignPtr DenseMatrix) deriving (Eq, Ord, Show)@@ -20,6 +22,7 @@ deriving (Eq, Ord, Show) foreign import ccall unsafe "svdNewSMat" _newSMat :: CInt -> CInt -> IO (Ptr SparseMatrix)+foreign import ccall unsafe "svd_new_smat_from_csr" _newSMatFromCSRT :: CInt -> CInt -> CInt -> Ptr CLong -> Ptr CLong -> Ptr Double -> IO (Ptr SparseMatrix) foreign import ccall unsafe "&svdFreeSMat" p_freeSMat :: FunPtr (Ptr SparseMatrix -> IO ()) foreign import ccall unsafe "svdTransposeS" _transposeSMat :: Ptr SparseMatrix -> IO (Ptr SparseMatrix) @@ -38,6 +41,12 @@ foreign import ccall unsafe "get_dmat_cols" getCols :: Ptr DenseMatrix -> IO CLong foreign import ccall unsafe "get_dmat_buffer" getBuffer :: Ptr DenseMatrix -> IO (Ptr Double) +foreign import ccall unsafe "get_smat_rows" getSRows :: Ptr SparseMatrix -> IO CLong+foreign import ccall unsafe "get_smat_cols" getSCols :: Ptr SparseMatrix -> IO CLong+foreign import ccall unsafe "get_smat_pointr" getSPointr :: Ptr SparseMatrix -> IO (Ptr CLong)+foreign import ccall unsafe "get_smat_rowind" getSRowind :: Ptr SparseMatrix -> IO (Ptr CLong) --long?+foreign import ccall unsafe "get_smat_buffer" getSBuffer :: Ptr SparseMatrix -> IO (Ptr Double)+ foreign import ccall unsafe "set_verbosity" setVerbosity :: CLong -> IO () -- Our approach to memory management for dmats isn't entirely future-proof as we currently@@ -49,9 +58,20 @@ asSMat :: Ptr SparseMatrix -> IO SparseMatrix asSMat ptr = SMat <$> newForeignPtr p_freeSMat ptr -createSMatrix :: Int -> Int -> IO SparseMatrix-createSMatrix rows cols = do- _newSMat (fromIntegral rows) (fromIntegral cols) >>= asSMat+wrapSMatrix :: I.CSR -> IO SparseMatrix+wrapSMatrix csr = do+ let tfst (x,_,_) = x+ Vec.unsafeWith (Vec.map fromIntegral $ I.csrRows csr) $ \rowvec ->+ Vec.unsafeWith (Vec.map fromIntegral $ I.csrCols csr) $ \colvec ->+ Vec.unsafeWith (I.csrVals csr) $ \valvec -> do+ smatptr <- _newSMatFromCSRT+ (fromIntegral $ I.csrNRows csr)+ (fromIntegral $ I.csrNCols csr)+ (fromIntegral $ P.size $ I.csrVals csr)+ rowvec+ colvec+ valvec+ asSMat smatptr transposeSMatrix :: SparseMatrix -> IO SparseMatrix transposeSMatrix (SMat fptr) = withForeignPtr fptr $ \ptr->@@ -98,4 +118,37 @@ -- This function handles the conversion to svdlibc's sparse representation. svd :: Int -> P.Matrix Double -> (P.Matrix Double, P.Vector Double, P.Matrix Double) svd rank m = unsafePerformIO $ do- setVerbosity 0 >> matrixToDMatrix m >>= dMatrixToSMatrix >>= runSvd rank >>= unpackSvdRec+ setVerbosity 0+ >> matrixToDMatrix m+ >>= dMatrixToSMatrix+ >>= runSvd rank+ >>= unpackSvdRec++sparsify :: P.Matrix Double -> I.CSR+sparsify mat = I.CSR {+ I.csrVals = P.flatten mat,+ -- The following are indexed from 1! This will be undone in sparseSvd(shiftCSR)+ -- but it needs to match HMatrix's standard until then+ I.csrCols = Vec.iterateN (rows*cols) (\x -> (x `mod` i32cols) + 1) 1,+ I.csrRows = Vec.iterateN (rows+1) (+i32cols) 1,+ I.csrNRows = rows,+ I.csrNCols = cols+ } where+ (rows, cols) = P.size mat+ (i32rows, i32cols) = (fromIntegral rows, fromIntegral cols)++shiftCSR :: I.CSR -> I.CSR+shiftCSR csr = csr {+ I.csrRows = Vec.map (subtract 1) $ I.csrRows csr,+ I.csrCols = Vec.map (subtract 1) $ I.csrCols csr+ }++-- | @svd rank a@ is the sparse SVD of matrix @a@ with the given rank+-- This function handles the conversion to svdlibc's sparse representation,+-- but does not require making the whole matrix dense first+sparseSvd :: Int -> I.CSR -> (P.Matrix Double, P.Vector Double, P.Matrix Double)+sparseSvd rank m = unsafePerformIO $+ setVerbosity 0+ >> (wrapSMatrix $ shiftCSR m)+ >>= runSvd rank+ >>= unpackSvdRec
− Test.hs
@@ -1,9 +0,0 @@-import Data.Packed as P-import Numeric.LinearAlgebra.SVD.SVDLIBC as SVD-import Numeric.LinearAlgebra--main = do- let m = ident 100- (u,s,v) = SVD.svd 50 m- print $ s- print $ u `mXm` diag s `mXm` trans v
cbits/glue.c view
@@ -2,6 +2,14 @@ #include <svdlib.h> #include <glue.h> +/* For reference, this is from svdlib.h+Row-major dense matrix. Rows are consecutive vectors.+struct dmat {+ long rows;+ long cols;+ double **value; // Accessed by [row][col]. Free value[0] and value to free.+}; */+ DMat get_svdrec_ut(SVDRec s) { return s->Ut; } double *get_svdrec_s(SVDRec s) { return s->S; } DMat get_svdrec_vt(SVDRec s) { return s->Vt; }@@ -15,5 +23,41 @@ free(d->value); free(d); }++/* For reference, this is in svdlib.h: Harwell-Boeing sparse matrix.+struct smat {+ long rows;+ long cols;+ long vals; // Total non-zero entries.+ long *pointr; // For each col (plus 1), index of first non-zero entry.+ long *rowind; // For each nz entry, the row index.+ double *value; // For each nz entry, the value.+}; */++/** Take a matrix in CSR format (which is what HMatrix gives) and turn it into+ CSC format, meanwhile making a copy so that GC is happy.*/+SMat svd_new_smat_from_csr(int rows, int cols, int vals, long *rowstart, long *colind, double *csr_value) {+ return svdTransposeS(& (struct smat){+ .rows = cols,+ .cols = rows,+ .vals = vals,+ .pointr = rowstart,+ .rowind = colind,+ .value = csr_value+ });+}++long get_smat_rows(SMat m) { return m->rows; };+long get_smat_cols(SMat m) { return m->cols; };+long *get_smat_pointr(SMat m) { return m->pointr; };+long *get_smat_rowind(SMat m) { return m->rowind; };+double *get_smat_value(SMat m) { return m->value; };++void free_smat(SMat m) {+ free(m->pointr);+ free(m->rowind);+ free(m->value);+ free(m);+}; void set_verbosity(long v) { SVDVerbosity = v; }
hmatrix-svdlibc.cabal view
@@ -1,5 +1,5 @@ name: hmatrix-svdlibc-version: 0.3.2+version: 0.4.0 synopsis: SVDLIBC bindings for HMatrix description: Bindings for the sparse singular value decomposition@@ -8,12 +8,13 @@ license: BSD3 license-file: LICENSE author: Ben Gamari-maintainer: bgamari.foss@gmail.com+maintainer: ben@smart-cactus.org copyright: (c) 2014 Ben Gamari category: Math build-type: Simple-cabal-version: >=1.8+cabal-version: >=1.10 extra-source-files: cbits/*.c, include/*.h, svdlibc/*.c, svdlibc/*.h, svdlibc/README.md+tested-with: GHC ==7.10.3, GHC ==8.0.1 source-repository head type: git@@ -25,12 +26,37 @@ Include-dirs: include, svdlibc Includes: glue.h, svdlib.h build-depends: base >=4.6 && <5.0,- hmatrix >=0.16 && <0.18+ hmatrix >=0.17 && <0.18,+ vector >= 0.11+ default-language: Haskell2010 -executable svdlibc-test- main-is: Test.hs+test-suite svdlibc-test+ type: exitcode-stdio-1.0+ hs-source-dirs: test c-sources: cbits/glue.c, svdlibc/svdlib.c, svdlibc/svdutil.c, svdlibc/las2.c include-dirs: include, svdlibc includes: glue.h, svdlib.h- build-depends: base >=4.6 && <5.0,- hmatrix >=0.16 && <0.18+ main-is: Numeric/LinearAlgebra/SVD/Spec.hs+ build-depends: hmatrix-svdlibc,+ base >=4.6 && <5.0,+ hmatrix >=0.17 && <0.18,+ vector >= 0.11,+ hspec >= 2.2,+ QuickCheck >= 2.8+ ghc-options: -threaded -rtsopts -with-rtsopts=-N+ default-language: Haskell2010++benchmark svdlibc-benchmarks+ type: exitcode-stdio-1.0+ hs-source-dirs: test+ c-sources: cbits/glue.c, svdlibc/svdlib.c, svdlibc/svdutil.c, svdlibc/las2.c+ include-dirs: include, svdlibc+ includes: glue.h, svdlib.h+ main-is: Numeric/LinearAlgebra/SVD/Benchmarks.hs+ build-depends: hmatrix-svdlibc,+ base >=4.6 && <5.0,+ hmatrix >=0.17 && <0.18,+ vector >= 0.11,+ criterion >= 1.1+ ghc-options: -threaded -rtsopts -with-rtsopts=-N+ default-language: Haskell2010
include/glue.h view
@@ -11,4 +11,14 @@ void free_dmat(DMat d); ++long get_smat_rows(SMat m);+long get_smat_cols(SMat m);+long *get_smat_pointr(SMat m);+long *get_smat_rowind(SMat m);+double *get_smat_value(SMat m);+SMat svd_new_smat_from_csr(int rows, int cols, int vals, long *pointr, long *rowind, double *value);++void free_smat(SMat m);+ void set_verbosity(long v);
+ test/Numeric/LinearAlgebra/SVD/Benchmarks.hs view
@@ -0,0 +1,19 @@+{-# LANGUAGE OverloadedStrings #-}+import Criterion.Main+import qualified Data.Vector.Storable as Vec+import Numeric.LinearAlgebra+import Numeric.LinearAlgebra.Data+import Numeric.LinearAlgebra.Devel+import Numeric.LinearAlgebra.SVD.SVDLIBC as SVD+import System.IO+++main :: IO ()+main = defaultMain+ [ bench "Factorizes a random dense matrix" $ nfIO $ do+ let (width, height) = (20, 10)+ SVD.svd (min width height) <$> randn width height+ , bench "Factorizes a random sparse matrix" $ nfIO $ do+ let (width, height) = (20, 10)+ SVD.sparseSvd (min width height) <$> SVD.sparsify <$> randn width height+ ]
+ test/Numeric/LinearAlgebra/SVD/Spec.hs view
@@ -0,0 +1,60 @@+{-# LANGUAGE OverloadedStrings #-}+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck hiding ((><))+import qualified Data.Vector.Storable as Vec+import Numeric.LinearAlgebra+import Numeric.LinearAlgebra.Data+import Numeric.LinearAlgebra.Devel+import Numeric.LinearAlgebra.SVD.SVDLIBC as SVD+import Debug.Trace+import System.IO++newtype RandHMat = RandHMat (Matrix Double) deriving Show++instance Arbitrary (RandHMat) where+ arbitrary = do+ x <- arbitrary+ y <- arbitrary+ vector <- Vec.replicateM (x*y) arbitrary+ return $ RandHMat $ reshape y vector++main :: IO ()+main = hspec $ do+ let summat = Vec.foldr (+) 0 . flatten+ validate u s v ref = summat (mconcat [tr u, diag s, v] - ref) / summat ref < 0.01++ describe "Dense matrices" $ do+ it "Factorizes a random dense matrix" $ do+ let (width, height) = (20, 10)+ mat <- randn width height+ short <- return $! min width height+ (u,s',v) <- return $! SVD.svd short mat+ s <- return $! vjoin [s', konst 0 (short - size s') ]+ validate u s v mat `shouldBe` True+++ prop "Factorizes a random dense matrix (quickcheck)" $ \(RandHMat mat) ->+ let (width, height) = size mat+ short = min width height+ (u,s',v) = SVD.svd short mat+ s = vjoin [s', konst 0 (short - size s') ]+ in (width > 1 && height > 1) ==> validate u s v mat++ describe "Sparse matrices" $ do+ it "Factorizes a random sparse matrix" $ do+ let (width, height) = (20, 10)+ m <- randn width height+ csr <- return $! SVD.sparsify m+ short <- return $! min width height+ (u,s',v) <- return $! SVD.sparseSvd short csr+ s <- return $! vjoin [s', konst 0 (short - size s') ]+ validate u s v m `shouldBe` True++ prop "Factorizes a random sparse matrix (quickcheck)" $ \(RandHMat mat) ->+ let (width, height) = size mat+ short = min width height+ csr = SVD.sparsify mat+ (u,s',v) = SVD.sparseSvd short csr+ s = vjoin [s', konst 0 (short - size s') ]+ in (width > 1 && height > 1) ==> validate u s v mat