ranksv and tests

2 files changed, 18 insertions, 10 deletions

diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index f749cc1..954b214 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -63,7 +63,7 @@ module Numeric.LinearAlgebra.Algorithms (
     haussholder,
     unpackQR, unpackHess,
     pinvTol,
-    rankSVD,
+    rankSVD, ranksv,
     nullspaceSVD
 ) where
 
@@ -262,15 +262,18 @@ rankSVD :: Element t
         -> Matrix t -- ^ input matrix m
         -> (Matrix t, Vector Double, Matrix t) -- ^ 'svd' of m
         -> Int      -- ^ rank of m
-rankSVD teps m (_,s,_) = k where
-    sl@(g:_) = toList s
-    r = rows m
-    c = cols m
-    tol = fromIntegral (max r c) * g * teps
-    s' = fromList . filter (>tol) $ sl
-    k = if g > teps
-        then dim s'
-             else 0
+rankSVD teps m (_,s,_) = ranksv teps (max (rows m) (cols m)) (toList s)
+
+-- | Numeric rank of a matrix from its singular values.
+ranksv ::  Double   -- ^ numeric zero (e.g. 1*'eps')
+        -> Int      -- ^ maximum dimension of the matrix
+        -> [Double] -- ^ singular values
+        -> Int      -- ^ rank of m
+ranksv teps maxdim s = k where
+    g = maximum s
+    tol = fromIntegral maxdim * g * teps
+    s' = filter (>tol) s
+    k = if g > teps then length s' else 0
 
 -- | The machine precision of a Double: @eps = 2.22044604925031e-16@ (the value used by GNU-Octave).
 eps :: Double
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs
index 4495396..bf29ee9 100644
--- a/lib/Numeric/LinearAlgebra/Tests.hs
+++ b/lib/Numeric/LinearAlgebra/Tests.hs
@@ -243,6 +243,11 @@ runTests n = do
         , rootFindingTest
         , utest "randomGaussian" randomTestGaussian
         , utest "randomUniform" randomTestUniform
+        , utest "buildVector/Matrix" $
+                        comp (10 |> [0::Double ..]) == buildVector 10 fromIntegral
+                     && ident 5 == buildMatrix 5 5 (\(r,c) -> if r==c then 1::Double else 0)
+        , utest "rank" $  rank ((2><3)[1,0,0,1,6*eps,0]) == 1
+                       && rank ((2><3)[1,0,0,1,7*eps,0]) == 2
         ]
     return ()