move algorithms

1 files changed, 920 insertions, 0 deletions

diff --git a/packages/base/src/Internal/Algorithms.hs b/packages/base/src/Internal/Algorithms.hs
new file mode 100644
index 0000000..328af22
--- /dev/null
+++ b/packages/base/src/Internal/Algorithms.hs
@@ -0,0 +1,920 @@
+{-# LANGUAGE FlexibleContexts, FlexibleInstances #-}
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE UndecidableInstances #-}
+{-# LANGUAGE TypeFamilies #-}
+
+-----------------------------------------------------------------------------
+{- |
+Module      :  Internal.Algorithms
+Copyright   :  (c) Alberto Ruiz 2006-14
+License     :  BSD3
+Maintainer  :  Alberto Ruiz
+Stability   :  provisional
+
+High level generic interface to common matrix computations.
+
+Specific functions for particular base types can also be explicitly
+imported from "Numeric.LinearAlgebra.LAPACK".
+
+-}
+-----------------------------------------------------------------------------
+
+module Internal.Algorithms where
+
+import Internal.Vector
+import Internal.Matrix
+import Internal.Element
+import Internal.Conversion
+import Internal.LAPACK as LAPACK
+import Data.List(foldl1')
+import Data.Array
+import Internal.Numeric
+import Data.Vector.Storable(fromList)
+
+-- :i mul
+
+{- | Generic linear algebra functions for double precision real and complex matrices.
+
+(Single precision data can be converted using 'single' and 'double').
+
+-}
+class (Product t,
+       Convert t,
+       Container Vector t,
+       Container Matrix t,
+       Normed Matrix t,
+       Normed Vector t,
+       Floating t,
+       RealOf t ~ Double) => Field t where
+    svd'         :: Matrix t -> (Matrix t, Vector Double, Matrix t)
+    thinSVD'     :: Matrix t -> (Matrix t, Vector Double, Matrix t)
+    sv'          :: Matrix t -> Vector Double
+    luPacked'    :: Matrix t -> (Matrix t, [Int])
+    luSolve'     :: (Matrix t, [Int]) -> Matrix t -> Matrix t
+    mbLinearSolve' :: Matrix t -> Matrix t -> Maybe (Matrix t)
+    linearSolve' :: Matrix t -> Matrix t -> Matrix t
+    cholSolve'   :: Matrix t -> Matrix t -> Matrix t
+    linearSolveSVD' :: Matrix t -> Matrix t -> Matrix t
+    linearSolveLS'  :: Matrix t -> Matrix t -> Matrix t
+    eig'         :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
+    eigSH''      :: Matrix t -> (Vector Double, Matrix t)
+    eigOnly      :: Matrix t -> Vector (Complex Double)
+    eigOnlySH    :: Matrix t -> Vector Double
+    cholSH'      :: Matrix t -> Matrix t
+    mbCholSH'    :: Matrix t -> Maybe (Matrix t)
+    qr'          :: Matrix t -> (Matrix t, Vector t)
+    qrgr'        :: Int -> (Matrix t, Vector t) -> Matrix t
+    hess'        :: Matrix t -> (Matrix t, Matrix t)
+    schur'       :: Matrix t -> (Matrix t, Matrix t)
+
+
+instance Field Double where
+    svd' = svdRd
+    thinSVD' = thinSVDRd
+    sv' = svR
+    luPacked' = luR
+    luSolve' (l_u,perm) = lusR l_u perm
+    linearSolve' = linearSolveR                 -- (luSolve . luPacked) ??
+    mbLinearSolve' = mbLinearSolveR
+    cholSolve' = cholSolveR
+    linearSolveLS' = linearSolveLSR
+    linearSolveSVD' = linearSolveSVDR Nothing
+    eig' = eigR
+    eigSH'' = eigS
+    eigOnly = eigOnlyR
+    eigOnlySH = eigOnlyS
+    cholSH' = cholS
+    mbCholSH' = mbCholS
+    qr' = qrR
+    qrgr' = qrgrR
+    hess' = unpackHess hessR
+    schur' = schurR
+
+instance Field (Complex Double) where
+#ifdef NOZGESDD
+    svd' = svdC
+    thinSVD' = thinSVDC
+#else
+    svd' = svdCd
+    thinSVD' = thinSVDCd
+#endif
+    sv' = svC
+    luPacked' = luC
+    luSolve' (l_u,perm) = lusC l_u perm
+    linearSolve' = linearSolveC
+    mbLinearSolve' = mbLinearSolveC
+    cholSolve' = cholSolveC
+    linearSolveLS' = linearSolveLSC
+    linearSolveSVD' = linearSolveSVDC Nothing
+    eig' = eigC
+    eigOnly = eigOnlyC
+    eigSH'' = eigH
+    eigOnlySH = eigOnlyH
+    cholSH' = cholH
+    mbCholSH' = mbCholH
+    qr' = qrC
+    qrgr' = qrgrC
+    hess' = unpackHess hessC
+    schur' = schurC
+
+--------------------------------------------------------------
+
+square m = rows m == cols m
+
+vertical m = rows m >= cols m
+
+exactHermitian m = m `equal` ctrans m
+
+--------------------------------------------------------------
+
+{- | Full singular value decomposition.
+
+@
+a = (5><3)
+ [  1.0,  2.0,  3.0
+ ,  4.0,  5.0,  6.0
+ ,  7.0,  8.0,  9.0
+ , 10.0, 11.0, 12.0
+ , 13.0, 14.0, 15.0 ] :: Matrix Double
+@
+
+>>> let (u,s,v) = svd a
+
+>>> disp 3 u
+5x5
+-0.101   0.768   0.614   0.028  -0.149
+-0.249   0.488  -0.503   0.172   0.646
+-0.396   0.208  -0.405  -0.660  -0.449
+-0.543  -0.072  -0.140   0.693  -0.447
+-0.690  -0.352   0.433  -0.233   0.398
+
+>>> s
+fromList [35.18264833189422,1.4769076999800903,1.089145439970417e-15]
+
+>>> disp 3 v
+3x3
+-0.519  -0.751   0.408
+-0.576  -0.046  -0.816
+-0.632   0.659   0.408
+
+>>> let d = diagRect 0 s 5 3
+>>> disp 3 d
+5x3
+35.183  0.000  0.000
+ 0.000  1.477  0.000
+ 0.000  0.000  0.000
+ 0.000  0.000  0.000
+
+>>> disp 3 $ u <> d <> tr v
+5x3
+ 1.000   2.000   3.000
+ 4.000   5.000   6.000
+ 7.000   8.000   9.000
+10.000  11.000  12.000
+13.000  14.000  15.000
+
+-}
+svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
+svd = {-# SCC "svd" #-} g . svd'
+  where
+    g (u,s,v) = (u,s,tr v)
+
+{- | A version of 'svd' which returns only the @min (rows m) (cols m)@ singular vectors of @m@.
+
+If @(u,s,v) = thinSVD m@ then @m == u \<> diag s \<> tr v@.
+
+@
+a = (5><3)
+ [  1.0,  2.0,  3.0
+ ,  4.0,  5.0,  6.0
+ ,  7.0,  8.0,  9.0
+ , 10.0, 11.0, 12.0
+ , 13.0, 14.0, 15.0 ] :: Matrix Double
+@
+
+>>> let (u,s,v) = thinSVD a
+
+>>> disp 3 u
+5x3
+-0.101   0.768   0.614
+-0.249   0.488  -0.503
+-0.396   0.208  -0.405
+-0.543  -0.072  -0.140
+-0.690  -0.352   0.433
+
+>>> s
+fromList [35.18264833189422,1.4769076999800903,1.089145439970417e-15]
+
+>>> disp 3 v
+3x3
+-0.519  -0.751   0.408
+-0.576  -0.046  -0.816
+-0.632   0.659   0.408
+
+>>> disp 3 $ u <> diag s <> tr v
+5x3
+ 1.000   2.000   3.000
+ 4.000   5.000   6.000
+ 7.000   8.000   9.000
+10.000  11.000  12.000
+13.000  14.000  15.000
+
+-}
+thinSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
+thinSVD = {-# SCC "thinSVD" #-} g . thinSVD'
+  where
+    g (u,s,v) = (u,s,tr v)
+
+
+-- | Singular values only.
+singularValues :: Field t => Matrix t -> Vector Double
+singularValues = {-# SCC "singularValues" #-} sv'
+
+-- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values.
+--
+-- If @(u,d,v) = fullSVD m@ then @m == u \<> d \<> tr v@.
+fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t)
+fullSVD m = (u,d,v) where
+    (u,s,v) = svd m
+    d = diagRect 0 s r c
+    r = rows m
+    c = cols m
+
+{- | Similar to 'thinSVD', returning only the nonzero singular values and the corresponding singular vectors.
+
+@
+a = (5><3)
+ [  1.0,  2.0,  3.0
+ ,  4.0,  5.0,  6.0
+ ,  7.0,  8.0,  9.0
+ , 10.0, 11.0, 12.0
+ , 13.0, 14.0, 15.0 ] :: Matrix Double
+@
+
+>>> let (u,s,v) = compactSVD a
+
+>>> disp 3 u
+5x2
+-0.101   0.768
+-0.249   0.488
+-0.396   0.208
+-0.543  -0.072
+-0.690  -0.352
+
+>>> s
+fromList [35.18264833189422,1.4769076999800903]
+
+>>> disp 3 u
+5x2
+-0.101   0.768
+-0.249   0.488
+-0.396   0.208
+-0.543  -0.072
+-0.690  -0.352
+
+>>> disp 3 $ u <> diag s <> tr v
+5x3
+ 1.000   2.000   3.000
+ 4.000   5.000   6.000
+ 7.000   8.000   9.000
+10.000  11.000  12.000
+13.000  14.000  15.000
+
+-}
+compactSVD :: Field t  => Matrix t -> (Matrix t, Vector Double, Matrix t)
+compactSVD m = (u', subVector 0 d s, v') where
+    (u,s,v) = thinSVD m
+    d = rankSVD (1*eps) m s `max` 1
+    u' = takeColumns d u
+    v' = takeColumns d v
+
+
+-- | Singular values and all right singular vectors (as columns).
+rightSV :: Field t => Matrix t -> (Vector Double, Matrix t)
+rightSV m | vertical m = let (_,s,v) = thinSVD m in (s,v)
+          | otherwise  = let (_,s,v) = svd m     in (s,v)
+
+-- | Singular values and all left singular vectors (as columns).
+leftSV :: Field t => Matrix t -> (Matrix t, Vector Double)
+leftSV m  | vertical m = let (u,s,_) = svd m     in (u,s)
+          | otherwise  = let (u,s,_) = thinSVD m in (u,s)
+
+
+--------------------------------------------------------------
+
+-- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'.
+luPacked :: Field t => Matrix t -> (Matrix t, [Int])
+luPacked = {-# SCC "luPacked" #-} luPacked'
+
+-- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization obtained by 'luPacked'.
+luSolve :: Field t => (Matrix t, [Int]) -> Matrix t -> Matrix t
+luSolve = {-# SCC "luSolve" #-} luSolve'
+
+-- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'.
+-- It is similar to 'luSolve' . 'luPacked', but @linearSolve@ raises an error if called on a singular system.
+linearSolve :: Field t => Matrix t -> Matrix t -> Matrix t
+linearSolve = {-# SCC "linearSolve" #-} linearSolve'
+
+-- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, returning Nothing for a singular system. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'. 
+mbLinearSolve :: Field t => Matrix t -> Matrix t -> Maybe (Matrix t)
+mbLinearSolve = {-# SCC "linearSolve" #-} mbLinearSolve'
+
+-- | Solve a symmetric or Hermitian positive definite linear system using a precomputed Cholesky decomposition obtained by 'chol'.
+cholSolve :: Field t => Matrix t -> Matrix t -> Matrix t
+cholSolve = {-# SCC "cholSolve" #-} cholSolve'
+
+-- | Minimum norm solution of a general linear least squares problem Ax=B using the SVD. Admits rank-deficient systems but it is slower than 'linearSolveLS'. The effective rank of A is determined by treating as zero those singular valures which are less than 'eps' times the largest singular value.
+linearSolveSVD :: Field t => Matrix t -> Matrix t -> Matrix t
+linearSolveSVD = {-# SCC "linearSolveSVD" #-} linearSolveSVD'
+
+
+-- | Least squared error solution of an overconstrained linear system, or the minimum norm solution of an underconstrained system. For rank-deficient systems use 'linearSolveSVD'.
+linearSolveLS :: Field t => Matrix t -> Matrix t -> Matrix t
+linearSolveLS = {-# SCC "linearSolveLS" #-} linearSolveLS'
+
+--------------------------------------------------------------
+
+{- | Eigenvalues (not ordered) and eigenvectors (as columns) of a general square matrix.
+
+If @(s,v) = eig m@ then @m \<> v == v \<> diag s@
+
+@
+a = (3><3)
+ [ 3, 0, -2
+ , 4, 5, -1
+ , 3, 1,  0 ] :: Matrix Double
+@
+
+>>> let (l, v) = eig a
+
+>>> putStr . dispcf 3 . asRow $ l
+1x3
+1.925+1.523i  1.925-1.523i  4.151
+
+>>> putStr . dispcf 3 $ v
+3x3
+-0.455+0.365i  -0.455-0.365i   0.181
+        0.603          0.603  -0.978
+ 0.033+0.543i   0.033-0.543i  -0.104
+
+>>> putStr . dispcf 3 $ complex a <> v
+3x3
+-1.432+0.010i  -1.432-0.010i   0.753
+ 1.160+0.918i   1.160-0.918i  -4.059
+-0.763+1.096i  -0.763-1.096i  -0.433
+
+>>> putStr . dispcf 3 $ v <> diag l
+3x3
+-1.432+0.010i  -1.432-0.010i   0.753
+ 1.160+0.918i   1.160-0.918i  -4.059
+-0.763+1.096i  -0.763-1.096i  -0.433
+
+-}
+eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
+eig = {-# SCC "eig" #-} eig'
+
+-- | Eigenvalues (not ordered) of a general square matrix.
+eigenvalues :: Field t => Matrix t -> Vector (Complex Double)
+eigenvalues = {-# SCC "eigenvalues" #-} eigOnly
+
+-- | Similar to 'eigSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
+eigSH' :: Field t => Matrix t -> (Vector Double, Matrix t)
+eigSH' = {-# SCC "eigSH'" #-} eigSH''
+
+-- | Similar to 'eigenvaluesSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
+eigenvaluesSH' :: Field t => Matrix t -> Vector Double
+eigenvaluesSH' = {-# SCC "eigenvaluesSH'" #-} eigOnlySH
+
+{- | Eigenvalues and eigenvectors (as columns) of a complex hermitian or real symmetric matrix, in descending order.
+
+If @(s,v) = eigSH m@ then @m == v \<> diag s \<> tr v@
+
+@
+a = (3><3)
+ [ 1.0, 2.0, 3.0
+ , 2.0, 4.0, 5.0
+ , 3.0, 5.0, 6.0 ]
+@
+
+>>> let (l, v) = eigSH a
+
+>>> l
+fromList [11.344814282762075,0.17091518882717918,-0.5157294715892575]
+
+>>> disp 3 $ v <> diag l <> tr v
+3x3
+1.000  2.000  3.000
+2.000  4.000  5.000
+3.000  5.000  6.000
+
+-}
+eigSH :: Field t => Matrix t -> (Vector Double, Matrix t)
+eigSH m | exactHermitian m = eigSH' m
+        | otherwise = error "eigSH requires complex hermitian or real symmetric matrix"
+
+-- | Eigenvalues (in descending order) of a complex hermitian or real symmetric matrix.
+eigenvaluesSH :: Field t => Matrix t -> Vector Double
+eigenvaluesSH m | exactHermitian m = eigenvaluesSH' m
+                | otherwise = error "eigenvaluesSH requires complex hermitian or real symmetric matrix"
+
+--------------------------------------------------------------
+
+-- | QR factorization.
+--
+-- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular.
+qr :: Field t => Matrix t -> (Matrix t, Matrix t)
+qr = {-# SCC "qr" #-} unpackQR . qr'
+
+qrRaw m = qr' m
+
+{- | generate a matrix with k orthogonal columns from the output of qrRaw
+-}
+qrgr n (a,t)
+    | dim t > min (cols a) (rows a) || n < 0 || n > dim t = error "qrgr expects k <= min(rows,cols)"
+    | otherwise = qrgr' n (a,t)
+
+-- | RQ factorization.
+--
+-- If @(r,q) = rq m@ then @m == r \<> q@, where q is unitary and r is upper triangular.
+rq :: Field t => Matrix t -> (Matrix t, Matrix t)
+rq m =  {-# SCC "rq" #-} (r,q) where
+    (q',r') = qr $ trans $ rev1 m
+    r = rev2 (trans r')
+    q = rev2 (trans q')
+    rev1 = flipud . fliprl
+    rev2 = fliprl . flipud
+
+-- | Hessenberg factorization.
+--
+-- If @(p,h) = hess m@ then @m == p \<> h \<> ctrans p@, where p is unitary
+-- and h is in upper Hessenberg form (it has zero entries below the first subdiagonal).
+hess        :: Field t => Matrix t -> (Matrix t, Matrix t)
+hess = hess'
+
+-- | Schur factorization.
+--
+-- If @(u,s) = schur m@ then @m == u \<> s \<> ctrans u@, where u is unitary
+-- and s is a Shur matrix. A complex Schur matrix is upper triangular. A real Schur matrix is
+-- upper triangular in 2x2 blocks.
+--
+-- \"Anything that the Jordan decomposition can do, the Schur decomposition
+-- can do better!\" (Van Loan)
+schur       :: Field t => Matrix t -> (Matrix t, Matrix t)
+schur = schur'
+
+
+-- | Similar to 'cholSH', but instead of an error (e.g., caused by a matrix not positive definite) it returns 'Nothing'.
+mbCholSH :: Field t => Matrix t -> Maybe (Matrix t)
+mbCholSH = {-# SCC "mbCholSH" #-} mbCholSH'
+
+-- | Similar to 'chol', without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
+cholSH      :: Field t => Matrix t -> Matrix t
+cholSH = {-# SCC "cholSH" #-} cholSH'
+
+-- | Cholesky factorization of a positive definite hermitian or symmetric matrix.
+--
+-- If @c = chol m@ then @c@ is upper triangular and @m == ctrans c \<> c@.
+chol :: Field t => Matrix t ->  Matrix t
+chol m | exactHermitian m = cholSH m
+       | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix"
+
+
+-- | Joint computation of inverse and logarithm of determinant of a square matrix.
+invlndet :: Field t
+         => Matrix t
+         -> (Matrix t, (t, t)) -- ^ (inverse, (log abs det, sign or phase of det))
+invlndet m | square m = (im,(ladm,sdm))
+           | otherwise = error $ "invlndet of nonsquare "++ shSize m ++ " matrix"
+  where
+    lp@(lup,perm) = luPacked m
+    s = signlp (rows m) perm
+    dg = toList $ takeDiag $ lup
+    ladm = sum $ map (log.abs) dg
+    sdm = s* product (map signum dg)
+    im = luSolve lp (ident (rows m))
+
+
+-- | Determinant of a square matrix. To avoid possible overflow or underflow use 'invlndet'.
+det :: Field t => Matrix t -> t
+det m | square m = {-# SCC "det" #-} s * (product $ toList $ takeDiag $ lup)
+      | otherwise = error $ "det of nonsquare "++ shSize m ++ " matrix"
+    where (lup,perm) = luPacked m
+          s = signlp (rows m) perm
+
+-- | Explicit LU factorization of a general matrix.
+--
+-- If @(l,u,p,s) = lu m@ then @m == p \<> l \<> u@, where l is lower triangular,
+-- u is upper triangular, p is a permutation matrix and s is the signature of the permutation.
+lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t)
+lu = luFact . luPacked
+
+-- | Inverse of a square matrix. See also 'invlndet'.
+inv :: Field t => Matrix t -> Matrix t
+inv m | square m = m `linearSolve` ident (rows m)
+      | otherwise = error $ "inv of nonsquare "++ shSize m ++ " matrix"
+
+
+-- | Pseudoinverse of a general matrix with default tolerance ('pinvTol' 1, similar to GNU-Octave).
+pinv :: Field t => Matrix t -> Matrix t
+pinv = pinvTol 1
+
+{- | @pinvTol r@ computes the pseudoinverse of a matrix with tolerance @tol=r*g*eps*(max rows cols)@, where g is the greatest singular value.
+
+@
+m = (3><3) [ 1, 0,    0
+           , 0, 1,    0
+           , 0, 0, 1e-10] :: Matrix Double
+@
+
+>>> pinv m
+1. 0.           0.
+0. 1.           0.
+0. 0. 10000000000.
+
+>>> pinvTol 1E8 m
+1. 0. 0.
+0. 1. 0.
+0. 0. 1.
+
+-}
+
+pinvTol :: Field t => Double -> Matrix t -> Matrix t
+pinvTol t m = v' `mXm` diag s' `mXm` ctrans u' where
+    (u,s,v) = thinSVD m
+    sl@(g:_) = toList s
+    s' = real . fromList . map rec $ sl
+    rec x = if x <= g*tol then x else 1/x
+    tol = (fromIntegral (max r c) * g * t * eps)
+    r = rows m
+    c = cols m
+    d = dim s
+    u' = takeColumns d u
+    v' = takeColumns d v
+
+
+-- | Numeric rank of a matrix from the SVD decomposition.
+rankSVD :: Element t
+        => Double   -- ^ numeric zero (e.g. 1*'eps')
+        -> Matrix t -- ^ input matrix m
+        -> Vector Double -- ^ 'sv' of m
+        -> Int      -- ^ rank of m
+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
+eps =  2.22044604925031e-16
+
+
+-- | 1 + 0.5*peps == 1,  1 + 0.6*peps /= 1
+peps :: RealFloat x => x
+peps = x where x = 2.0 ** fromIntegral (1 - floatDigits x)
+
+
+-- | The imaginary unit: @i = 0.0 :+ 1.0@
+i :: Complex Double
+i = 0:+1
+
+-----------------------------------------------------------------------
+
+-- | The nullspace of a matrix from its precomputed SVD decomposition.
+nullspaceSVD :: Field t
+             => Either Double Int -- ^ Left \"numeric\" zero (eg. 1*'eps'),
+                                  --   or Right \"theoretical\" matrix rank.
+             -> Matrix t          -- ^ input matrix m
+             -> (Vector Double, Matrix t) -- ^ 'rightSV' of m
+             -> Matrix t          -- ^ nullspace
+nullspaceSVD hint a (s,v) = vs where
+    tol = case hint of
+        Left t -> t
+        _      -> eps
+    k = case hint of
+        Right t -> t
+        _       -> rankSVD tol a s
+    vs = dropColumns k v
+
+
+-- | The nullspace of a matrix. See also 'nullspaceSVD'.
+nullspacePrec :: Field t
+              => Double     -- ^ relative tolerance in 'eps' units (e.g., use 3 to get 3*'eps')
+              -> Matrix t   -- ^ input matrix
+              -> [Vector t] -- ^ list of unitary vectors spanning the nullspace
+nullspacePrec t m = toColumns $ nullspaceSVD (Left (t*eps)) m (rightSV m)
+
+-- | The nullspace of a matrix, assumed to be one-dimensional, with machine precision.
+nullVector :: Field t => Matrix t -> Vector t
+nullVector = last . nullspacePrec 1
+
+-- | The range space a matrix from its precomputed SVD decomposition.
+orthSVD :: Field t
+             => Either Double Int -- ^ Left \"numeric\" zero (eg. 1*'eps'),
+                                  --   or Right \"theoretical\" matrix rank.
+             -> Matrix t          -- ^ input matrix m
+             -> (Matrix t, Vector Double) -- ^ 'leftSV' of m
+             -> Matrix t          -- ^ orth
+orthSVD hint a (v,s) = vs where
+    tol = case hint of
+        Left t -> t
+        _      -> eps
+    k = case hint of
+        Right t -> t
+        _       -> rankSVD tol a s
+    vs = takeColumns k v
+
+
+orth :: Field t => Matrix t -> [Vector t]
+-- ^ Return an orthonormal basis of the range space of a matrix
+orth m = take r $ toColumns u
+  where
+    (u,s,_) = compactSVD m
+    r = ranksv eps (max (rows m) (cols m)) (toList s)
+
+------------------------------------------------------------------------
+
+-- many thanks, quickcheck!
+
+haussholder :: (Field a) => a -> Vector a -> Matrix a
+haussholder tau v = ident (dim v) `sub` (tau `scale` (w `mXm` ctrans w))
+    where w = asColumn v
+
+
+zh k v = fromList $ replicate (k-1) 0 ++ (1:drop k xs)
+              where xs = toList v
+
+zt 0 v = v
+zt k v = vjoin [subVector 0 (dim v - k) v, konst' 0 k]
+
+
+unpackQR :: (Field t) => (Matrix t, Vector t) -> (Matrix t, Matrix t)
+unpackQR (pq, tau) =  {-# SCC "unpackQR" #-} (q,r)
+    where cs = toColumns pq
+          m = rows pq
+          n = cols pq
+          mn = min m n
+          r = fromColumns $ zipWith zt ([m-1, m-2 .. 1] ++ repeat 0) cs
+          vs = zipWith zh [1..mn] cs
+          hs = zipWith haussholder (toList tau) vs
+          q = foldl1' mXm hs
+
+unpackHess :: (Field t) => (Matrix t -> (Matrix t,Vector t)) -> Matrix t -> (Matrix t, Matrix t)
+unpackHess hf m
+    | rows m == 1 = ((1><1)[1],m)
+    | otherwise = (uH . hf) m
+
+uH (pq, tau) = (p,h)
+    where cs = toColumns pq
+          m = rows pq
+          n = cols pq
+          mn = min m n
+          h = fromColumns $ zipWith zt ([m-2, m-3 .. 1] ++ repeat 0) cs
+          vs = zipWith zh [2..mn] cs
+          hs = zipWith haussholder (toList tau) vs
+          p = foldl1' mXm hs
+
+--------------------------------------------------------------------------
+
+-- | Reciprocal of the 2-norm condition number of a matrix, computed from the singular values.
+rcond :: Field t => Matrix t -> Double
+rcond m = last s / head s
+    where s = toList (singularValues m)
+
+-- | Number of linearly independent rows or columns. See also 'ranksv'
+rank :: Field t => Matrix t -> Int
+rank m = rankSVD eps m (singularValues m)
+
+{-
+expm' m = case diagonalize (complex m) of
+    Just (l,v) -> v `mXm` diag (exp l) `mXm` inv v
+    Nothing -> error "Sorry, expm not yet implemented for non-diagonalizable matrices"
+  where exp = vectorMapC Exp
+-}
+
+diagonalize m = if rank v == n
+                    then Just (l,v)
+                    else Nothing
+    where n = rows m
+          (l,v) = if exactHermitian m
+                    then let (l',v') = eigSH m in (real l', v')
+                    else eig m
+
+-- | Generic matrix functions for diagonalizable matrices. For instance:
+--
+-- @logm = matFunc log@
+--
+matFunc :: (Complex Double -> Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
+matFunc f m = case diagonalize m of
+    Just (l,v) -> v `mXm` diag (mapVector f l) `mXm` inv v
+    Nothing -> error "Sorry, matFunc requires a diagonalizable matrix"
+
+--------------------------------------------------------------
+
+golubeps :: Integer -> Integer -> Double
+golubeps p q = a * fromIntegral b / fromIntegral c where
+    a = 2^^(3-p-q)
+    b = fact p * fact q
+    c = fact (p+q) * fact (p+q+1)
+    fact n = product [1..n]
+
+epslist :: [(Int,Double)]
+epslist = [ (fromIntegral k, golubeps k k) | k <- [1..]]
+
+geps delta = head [ k | (k,g) <- epslist, g<delta]
+
+
+{- | Matrix exponential. It uses a direct translation of Algorithm 11.3.1 in Golub & Van Loan,
+     based on a scaled Pade approximation.
+-}
+expm :: Field t => Matrix t -> Matrix t
+expm = expGolub
+
+expGolub :: Field t => Matrix t -> Matrix t
+expGolub m = iterate msq f !! j
+    where j = max 0 $ floor $ logBase 2 $ pnorm Infinity m
+          a = m */ fromIntegral ((2::Int)^j)
+          q = geps eps -- 7 steps
+          eye = ident (rows m)
+          work (k,c,x,n,d) = (k',c',x',n',d')
+              where k' = k+1
+                    c' = c * fromIntegral (q-k+1) / fromIntegral ((2*q-k+1)*k)
+                    x' = a <> x
+                    n' = n |+| (c' .* x')
+                    d' = d |+| (((-1)^k * c') .* x')
+          (_,_,_,nf,df) = iterate work (1,1,eye,eye,eye) !! q
+          f = linearSolve df nf
+          msq x = x <> x
+
+          (<>) = multiply
+          v */ x = scale (recip x) v
+          (.*) = scale
+          (|+|) = add
+
+--------------------------------------------------------------
+
+{- | Matrix square root. Currently it uses a simple iterative algorithm described in Wikipedia.
+It only works with invertible matrices that have a real solution.
+
+@m = (2><2) [4,9
+           ,0,4] :: Matrix Double@
+
+>>> sqrtm m
+(2><2)
+ [ 2.0, 2.25
+ , 0.0,  2.0 ]
+
+For diagonalizable matrices you can try 'matFunc' @sqrt@:
+
+>>> matFunc sqrt ((2><2) [1,0,0,-1])
+(2><2)
+ [ 1.0 :+ 0.0, 0.0 :+ 0.0
+ , 0.0 :+ 0.0, 0.0 :+ 1.0 ]
+
+-}
+sqrtm ::  Field t => Matrix t -> Matrix t
+sqrtm = sqrtmInv
+
+sqrtmInv x = fst $ fixedPoint $ iterate f (x, ident (rows x))
+    where fixedPoint (a:b:rest) | pnorm PNorm1 (fst a |-| fst b) < peps   = a
+                                | otherwise = fixedPoint (b:rest)
+          fixedPoint _ = error "fixedpoint with impossible inputs"
+          f (y,z) = (0.5 .* (y |+| inv z),
+                     0.5 .* (inv y |+| z))
+          (.*) = scale
+          (|+|) = add
+          (|-|) = sub
+
+------------------------------------------------------------------
+
+signlp r vals = foldl f 1 (zip [0..r-1] vals)
+    where f s (a,b) | a /= b    = -s
+                    | otherwise =  s
+
+swap (arr,s) (a,b) | a /= b    = (arr // [(a, arr!b),(b,arr!a)],-s)
+                   | otherwise = (arr,s)
+
+fixPerm r vals = (fromColumns $ elems res, sign)
+    where v = [0..r-1]
+          s = toColumns (ident r)
+          (res,sign) = foldl swap (listArray (0,r-1) s, 1) (zip v vals)
+
+triang r c h v = (r><c) [el s t | s<-[0..r-1], t<-[0..c-1]]
+    where el p q = if q-p>=h then v else 1 - v
+
+luFact (l_u,perm) | r <= c    = (l ,u ,p, s)
+                  | otherwise = (l',u',p, s)
+  where
+    r = rows l_u
+    c = cols l_u
+    tu = triang r c 0 1
+    tl = triang r c 0 0
+    l = takeColumns r (l_u |*| tl) |+| diagRect 0 (konst' 1 r) r r
+    u = l_u |*| tu
+    (p,s) = fixPerm r perm
+    l' = (l_u |*| tl) |+| diagRect 0 (konst' 1 c) r c
+    u' = takeRows c (l_u |*| tu)
+    (|+|) = add
+    (|*|) = mul
+
+---------------------------------------------------------------------------
+
+data NormType = Infinity | PNorm1 | PNorm2 | Frobenius
+
+class (RealFloat (RealOf t)) => Normed c t where
+    pnorm :: NormType -> c t -> RealOf t
+
+instance Normed Vector Double where
+    pnorm PNorm1    = norm1
+    pnorm PNorm2    = norm2
+    pnorm Infinity  = normInf
+    pnorm Frobenius = norm2
+
+instance Normed Vector (Complex Double) where
+    pnorm PNorm1    = norm1
+    pnorm PNorm2    = norm2
+    pnorm Infinity  = normInf
+    pnorm Frobenius = pnorm PNorm2
+
+instance Normed Vector Float where
+    pnorm PNorm1    = norm1
+    pnorm PNorm2    = norm2
+    pnorm Infinity  = normInf
+    pnorm Frobenius = pnorm PNorm2
+
+instance Normed Vector (Complex Float) where
+    pnorm PNorm1    = norm1
+    pnorm PNorm2    = norm2
+    pnorm Infinity  = normInf
+    pnorm Frobenius = pnorm PNorm2
+
+
+instance Normed Matrix Double where
+    pnorm PNorm1    = maximum . map (pnorm PNorm1) . toColumns
+    pnorm PNorm2    = (@>0) . singularValues
+    pnorm Infinity  = pnorm PNorm1 . trans
+    pnorm Frobenius = pnorm PNorm2 . flatten
+
+instance Normed Matrix (Complex Double) where
+    pnorm PNorm1    = maximum . map (pnorm PNorm1) . toColumns
+    pnorm PNorm2    = (@>0) . singularValues
+    pnorm Infinity  = pnorm PNorm1 . trans
+    pnorm Frobenius = pnorm PNorm2 . flatten
+
+instance Normed Matrix Float where
+    pnorm PNorm1    = maximum . map (pnorm PNorm1) . toColumns
+    pnorm PNorm2    = realToFrac . (@>0) . singularValues . double
+    pnorm Infinity  = pnorm PNorm1 . trans
+    pnorm Frobenius = pnorm PNorm2 . flatten
+
+instance Normed Matrix (Complex Float) where
+    pnorm PNorm1    = maximum . map (pnorm PNorm1) . toColumns
+    pnorm PNorm2    = realToFrac . (@>0) . singularValues . double
+    pnorm Infinity  = pnorm PNorm1 . trans
+    pnorm Frobenius = pnorm PNorm2 . flatten
+
+-- | Approximate number of common digits in the maximum element.
+relativeError' :: (Normed c t, Container c t) => c t -> c t -> Int
+relativeError' x y = dig (norm (x `sub` y) / norm x)
+    where norm = pnorm Infinity
+          dig r = round $ -logBase 10 (realToFrac r :: Double)
+
+
+relativeError :: Num a => (a -> Double) -> a -> a -> Double
+relativeError norm a b = r
+  where
+    na = norm a
+    nb = norm b
+    nab = norm (a-b)
+    mx = max na nb
+    mn = min na nb
+    r = if mn < peps
+        then mx
+        else nab/mx
+
+
+----------------------------------------------------------------------
+
+-- | Generalized symmetric positive definite eigensystem Av = lBv,
+-- for A and B symmetric, B positive definite (conditions not checked).
+geigSH' :: Field t
+        => Matrix t -- ^ A
+        -> Matrix t -- ^ B
+        -> (Vector Double, Matrix t)
+geigSH' a b = (l,v')
+  where
+    u = cholSH b
+    iu = inv u
+    c = ctrans iu <> a <> iu
+    (l,v) = eigSH' c
+    v' = iu <> v
+    (<>) = mXm
+