diff --git a/Changelog.md b/Changelog.md new file mode 100644 --- /dev/null +++ b/Changelog.md @@ -0,0 +1,7 @@ +# hmatrix-svdlibc + +## 0.5.0 + +* Added `svd'`, `sparseSvd'`, `SVDParams`, and `defaultSVDParams` providing + greater control over decomposition parameters and convergence. + 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,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 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.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 diff --git a/svdlibc/las2.c b/svdlibc/las2.c --- a/svdlibc/las2.c +++ b/svdlibc/las2.c @@ -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);