From 3a058b3707eecaac8ee3d964baf3e1ea1faabf51 Mon Sep 17 00:00:00 2001
From: Alberto Ruiz <aruiz@um.es>
Date: Wed, 25 Jul 2007 09:59:44 +0000
Subject: nullspace

---
 lib/LinearAlgebra/Algorithms.hs | 19 +++++++++++++++++++
 1 file changed, 19 insertions(+)

(limited to 'lib/LinearAlgebra/Algorithms.hs')

diff --git a/lib/LinearAlgebra/Algorithms.hs b/lib/LinearAlgebra/Algorithms.hs
index 126549a..4ac6dde 100644
--- a/lib/LinearAlgebra/Algorithms.hs
+++ b/lib/LinearAlgebra/Algorithms.hs
@@ -19,6 +19,8 @@ module LinearAlgebra.Algorithms (
     pinv,
     pinvTol,
     pinvTolg,
+    nullspacePrec,
+    nullVector,
     Normed(..), NormType(..),
     det,
     eps, i
@@ -240,3 +242,20 @@ instance Normed (Matrix Double) where
 
 instance Normed (Matrix (Complex Double)) where
     pnorm = pnormCM
+
+-----------------------------------------------------------------------
+
+-- | The nullspace of a real matrix from its SVD decomposition.
+nullspacePrec :: Double          -- ^ relative tolerance in 'eps' units
+              -> Matrix Double   -- ^ input matrix
+              -> [Vector Double] -- ^ list of unitary vectors spanning the nullspace
+nullspacePrec t m = ns where
+    (_,s,v) = svdR' m
+    sl@(g:_) = toList s
+    tol = (fromIntegral (max (rows m) (cols m)) * g * t * eps)
+    rank = length (filter (> g*tol) sl)
+    ns = drop rank (toColumns v)
+
+-- | The nullspace of a real matrix, assumed to be one-dimensional, with default tolerance (shortcut for @last . nullspacePrec 1@).
+nullVector :: Matrix Double -> Vector Double
+nullVector = last . nullspacePrec 1
-- 
cgit v1.2.3