diff --git a/Numeric/LinearAlgebra/SVD/SVDLIBC.hs b/Numeric/LinearAlgebra/SVD/SVDLIBC.hs
--- a/Numeric/LinearAlgebra/SVD/SVDLIBC.hs
+++ b/Numeric/LinearAlgebra/SVD/SVDLIBC.hs
@@ -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
diff --git a/Test.hs b/Test.hs
deleted file mode 100644
--- a/Test.hs
+++ /dev/null
@@ -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
diff --git a/cbits/glue.c b/cbits/glue.c
--- a/cbits/glue.c
+++ b/cbits/glue.c
@@ -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; }
diff --git a/hmatrix-svdlibc.cabal b/hmatrix-svdlibc.cabal
--- a/hmatrix-svdlibc.cabal
+++ b/hmatrix-svdlibc.cabal
@@ -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
diff --git a/include/glue.h b/include/glue.h
--- a/include/glue.h
+++ b/include/glue.h
@@ -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);
diff --git a/test/Numeric/LinearAlgebra/SVD/Benchmarks.hs b/test/Numeric/LinearAlgebra/SVD/Benchmarks.hs
new file mode 100644
--- /dev/null
+++ b/test/Numeric/LinearAlgebra/SVD/Benchmarks.hs
@@ -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
+  ]
diff --git a/test/Numeric/LinearAlgebra/SVD/Spec.hs b/test/Numeric/LinearAlgebra/SVD/Spec.hs
new file mode 100644
--- /dev/null
+++ b/test/Numeric/LinearAlgebra/SVD/Spec.hs
@@ -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