hmatrix-svdlibc 0.4.1 → 0.5.0
raw patch · 4 files changed
+61/−16 lines, 4 filesdep ~hmatrix
Dependency ranges changed: hmatrix
Files
- Changelog.md +7/−0
- Numeric/LinearAlgebra/SVD/SVDLIBC.hs +47/−10
- hmatrix-svdlibc.cabal +6/−6
- svdlibc/las2.c +1/−0
+ Changelog.md view
@@ -0,0 +1,7 @@+# hmatrix-svdlibc++## 0.5.0++* Added `svd'`, `sparseSvd'`, `SVDParams`, and `defaultSVDParams` providing+ greater control over decomposition parameters and convergence.+
Numeric/LinearAlgebra/SVD/SVDLIBC.hs view
@@ -1,7 +1,14 @@ {-# LANGUAGE ForeignFunctionInterface #-} module Numeric.LinearAlgebra.SVD.SVDLIBC- (svd, sparsify, sparseSvd) where+ ( -- * Singular value decomposition+ svd, svd'+ -- * Sparse SVD+ , sparsify+ , sparseSvd, sparseSvd'+ -- * Parameters+ , SVDParams(..), defaultSVDParams+ ) where import Control.Applicative import qualified Numeric.LinearAlgebra.Data as P@@ -30,7 +37,8 @@ foreign import ccall unsafe "svdConvertStoD" _convertSToD :: Ptr SparseMatrix -> IO (Ptr DenseMatrix) newtype SVDRec = SVDRec (ForeignPtr SVDRec)-foreign import ccall unsafe "svdLAS2A" _svdLAS2 :: Ptr SparseMatrix -> CLong -> IO (Ptr SVDRec)+foreign import ccall unsafe "svdLAS2" _svdLAS2 :: Ptr SparseMatrix -> CLong+ -> CLong -> Ptr CDouble -> CDouble -> IO (Ptr SVDRec) foreign import ccall unsafe "get_svdrec_ut" getUt :: Ptr SVDRec -> IO (Ptr DenseMatrix) foreign import ccall unsafe "get_svdrec_s" getS :: Ptr SVDRec -> IO (Ptr Double)@@ -100,9 +108,27 @@ sMatrixToDMatrix (SMat fptr) = withForeignPtr fptr $ \ptr-> _convertSToD ptr >>= asDMat -runSvd :: Int -> SparseMatrix -> IO SVDRec-runSvd rank (SMat fptr) = withForeignPtr fptr $ \ptr->do- res <- _svdLAS2 ptr (fromIntegral rank)+data SVDParams = SVDParams { maxIters :: Maybe Int+ -- ^ Maximum iteration count+ , minEigenvalRange :: (Double, Double)+ -- ^ Eigenvalues in this range are considered uninteresting+ , kappa :: Double+ }++-- | No iteration limit, exclude eigenvalues in range \((-10^{-30}, +10^{-30})@,+-- kappa of @10^{-6}\).+defaultSVDParams :: SVDParams+defaultSVDParams = SVDParams { maxIters = Nothing+ , minEigenvalRange = (-1e-30, 1e-30)+ , kappa = 1e-6+ }++runSvd :: SVDParams -> Int -> SparseMatrix -> IO SVDRec+runSvd params rank (SMat fptr) = withForeignPtr fptr $ \ptr->do+ let iterLimit = maybe 0 fromIntegral (maxIters params)+ (a,b) = minEigenvalRange params+ res <- withArray [realToFrac a, realToFrac b] $ \eigenvalRangePtr ->+ _svdLAS2 ptr (fromIntegral rank) iterLimit eigenvalRangePtr (realToFrac $ kappa params) SVDRec <$> newForeignPtr finalizerFree res unpackSvdRec :: SVDRec -> IO (P.Matrix Double, P.Vector Double, P.Matrix Double)@@ -117,11 +143,16 @@ -- | @svd rank a@ is the sparse SVD of matrix @a@ with the given rank -- 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+svd = svd' defaultSVDParams++-- | Similar to 'svd' but allowing control over the 'SVDParams'.+svd' :: SVDParams -> Int -> P.Matrix Double+ -> (P.Matrix Double, P.Vector Double, P.Matrix Double)+svd' params rank m = unsafePerformIO $ do setVerbosity 0 >> matrixToDMatrix m >>= dMatrixToSMatrix- >>= runSvd rank+ >>= runSvd params rank >>= unpackSvdRec sparsify :: P.Matrix Double -> I.CSR@@ -143,12 +174,18 @@ I.csrCols = Vec.map (subtract 1) $ I.csrCols csr } --- | @svd rank a@ is the sparse SVD of matrix @a@ with the given rank+-- | @sparseSvd 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 $+sparseSvd = sparseSvd' defaultSVDParams++-- | Like 'sparseSvd' but providing greater control over the 'SVDParams' used+-- for the computation.+sparseSvd' :: SVDParams -> Int -> I.CSR+ -> (P.Matrix Double, P.Vector Double, P.Matrix Double)+sparseSvd' params rank m = unsafePerformIO $ setVerbosity 0 >> (wrapSMatrix $ shiftCSR m)- >>= runSvd rank+ >>= runSvd params rank >>= unpackSvdRec
hmatrix-svdlibc.cabal view
@@ -1,5 +1,5 @@ name: hmatrix-svdlibc-version: 0.4.1+version: 0.5.0 synopsis: SVDLIBC bindings for HMatrix description: Bindings for the sparse singular value decomposition@@ -13,8 +13,8 @@ category: Math build-type: Simple 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+extra-source-files: Changelog.md, cbits/*.c, include/*.h, svdlibc/*.c, svdlibc/*.h, svdlibc/README.md+tested-with: GHC ==8.0.2, GHC ==8.2.2, GHC ==8.4.3, GHC ==8.6.1 source-repository head type: git@@ -26,7 +26,7 @@ Include-dirs: include, svdlibc Includes: glue.h, svdlib.h build-depends: base >=4.6 && <5.0,- hmatrix >=0.17 && <0.19,+ hmatrix >=0.17 && <0.20, vector >= 0.11 default-language: Haskell2010 @@ -39,7 +39,7 @@ main-is: Numeric/LinearAlgebra/SVD/Spec.hs build-depends: hmatrix-svdlibc, base >=4.6 && <5.0,- hmatrix >=0.17 && <0.19,+ hmatrix >=0.17 && <0.20, vector >= 0.11, hspec >= 2.2, QuickCheck >= 2.8@@ -55,7 +55,7 @@ main-is: Numeric/LinearAlgebra/SVD/Benchmarks.hs build-depends: hmatrix-svdlibc, base >=4.6 && <5.0,- hmatrix >=0.17 && <0.19,+ hmatrix >=0.17 && <0.20, vector >= 0.11, criterion >= 1.1 ghc-options: -threaded -rtsopts -with-rtsopts=-N
svdlibc/las2.c view
@@ -331,6 +331,7 @@ SVDRec R = NULL; ierr = 0; // reset the global error flag + svdWriteSparseMatrix(A, "svd.in.m", SVD_F_DT); svdResetCounters(); m = svd_imin(A->rows, A->cols);