check hermitian on chol and eigSH

3 files changed, 28 insertions, 19 deletions

diff --git a/lib/LinearAlgebra/Algorithms.hs b/lib/LinearAlgebra/Algorithms.hs
index 481dbd6..682b17f 100644
--- a/lib/LinearAlgebra/Algorithms.hs
+++ b/lib/LinearAlgebra/Algorithms.hs
@@ -28,7 +28,7 @@ module LinearAlgebra.Algorithms (
     svd,
     full, economy,
 -- ** Eigensystems
-    eig, LinearAlgebra.Algorithms.eigS, LinearAlgebra.Algorithms.eigH,
+    eig, eigSH,
 -- ** Other 
     LinearAlgebra.Algorithms.qr, chol,
 -- * Nullspace
@@ -38,7 +38,7 @@ module LinearAlgebra.Algorithms (
     eps, i,
     ctrans,
     Normed(..), NormType(..),
-    GenMat(linearSolveSVD,lu,eigSH)
+    GenMat(linearSolveSVD,lu,eigSH')
 ) where
 
 
@@ -50,15 +50,15 @@ import LinearAlgebra.LAPACK as LAPACK
 import Complex
 import LinearAlgebra.Linear
 
--- | matrix computations available for both real and complex matrices
+-- | Auxiliary typeclass used to define generic computations for both real and complex matrices.
 class (Linear Matrix t) => GenMat t where
     svd         :: Matrix t -> (Matrix t, Vector Double, Matrix t)
     lu          :: Matrix t -> (Matrix t, Matrix t, [Int], t)
     linearSolve :: Matrix t -> Matrix t -> Matrix t
     linearSolveSVD :: Matrix t -> Matrix t -> Matrix t
     eig         :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
-    eigSH       :: Matrix t -> (Vector Double, Matrix t)
-    chol        :: Matrix t -> Matrix t
+    eigSH'      :: Matrix t -> (Vector Double, Matrix t)
+    cholSH      :: Matrix t -> Matrix t
     -- | conjugate transpose
     ctrans :: Matrix t -> Matrix t
 
@@ -69,8 +69,8 @@ instance GenMat Double where
     linearSolveSVD = linearSolveSVDR Nothing
     ctrans = trans
     eig = eigR
-    eigSH = LAPACK.eigS
-    chol = cholS
+    eigSH' = eigS
+    cholSH = cholS
 
 instance GenMat (Complex Double) where
     svd = svdC
@@ -79,16 +79,18 @@ instance GenMat (Complex Double) where
     linearSolveSVD = linearSolveSVDC Nothing
     ctrans = conjTrans
     eig = eigC
-    eigSH =  LAPACK.eigH
-    chol = cholH
-
--- | eigensystem of a symmetric matrix
-eigS :: Matrix Double -> (Vector Double, Matrix Double)
-eigS = LAPACK.eigS
-
--- | eigensystem of a hermitian matrix
-eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
-eigH = LAPACK.eigH
+    eigSH' = eigH
+    cholSH = cholH
+
+-- | eigensystem of complex hermitian or real symmetric matrix
+eigSH :: GenMat t => Matrix t -> (Vector Double, Matrix t)
+eigSH m | m `equal` ctrans m = eigSH' m
+        | otherwise = error "eigSH requires complex hermitian or real symmetric matrix"
+
+-- | Cholesky factorization of a positive definite hermitian or symmetric matrix
+chol :: GenMat t => Matrix t ->  Matrix t
+chol m | m `equal` ctrans m = cholSH m
+       | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix"
 
 qr :: Matrix Double -> (Matrix Double, Matrix Double)
 qr = GSL.Matrix.qr
diff --git a/lib/LinearAlgebra/LAPACK.hs b/lib/LinearAlgebra/LAPACK.hs
index 2861061..c6908e3 100644
--- a/lib/LinearAlgebra/LAPACK.hs
+++ b/lib/LinearAlgebra/LAPACK.hs
@@ -15,7 +15,7 @@
 
 module LinearAlgebra.LAPACK (
     svdR, svdRdd, svdC,
-    eigC, eigR, eigS, eigH,
+    eigC, eigR, eigS, eigH, eigS', eigH',
     linearSolveR, linearSolveC,
     linearSolveLSR, linearSolveLSC,
     linearSolveSVDR, linearSolveSVDC,
@@ -166,7 +166,7 @@ foreign import ccall "LAPACK/lapack-aux.h eig_l_H" zheev :: TCMVCM
 -- if @(l,v)=eigH m@ then @m \<\> s v = v \<\> diag l@.
 --
 -- The eigenvectors are the columns of v.
--- The eigenvalues are sorted in descending order.
+-- The eigenvalues are sorted in descending order (use eigH' for ascending order).
 eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
 eigH m = (s', fliprl v)
     where (s,v) = eigH' m
diff --git a/lib/LinearAlgebra/Linear.hs b/lib/LinearAlgebra/Linear.hs
index 3e6c55d..7cd1a52 100644
--- a/lib/LinearAlgebra/Linear.hs
+++ b/lib/LinearAlgebra/Linear.hs
@@ -37,6 +37,8 @@ class (Container c e) => Linear c e where
     divide      :: c e -> c e -> c e
     -- | scale the element by element reciprocal of the object: @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@
     scaleRecip  :: e -> c e -> c e
+    equal       :: c e -> c e -> Bool
+--  numequal    :: Double -> c e -> c e -> Bool
 
 instance Linear Vector Double where
     scale = vectorMapValR Scale
@@ -46,6 +48,7 @@ instance Linear Vector Double where
     sub = vectorZipR Sub
     mul = vectorZipR Mul
     divide = vectorZipR Div
+    equal u v = dim u == dim v && vectorMax (vectorMapR Abs (sub u v)) == 0.0
 
 instance Linear Vector (Complex Double) where
     scale = vectorMapValC Scale
@@ -55,6 +58,7 @@ instance Linear Vector (Complex Double) where
     sub = vectorZipC Sub
     mul = vectorZipC Mul
     divide = vectorZipC Div
+    equal u v = dim u == dim v && vectorMax (liftVector magnitude (sub u v)) == 0.0
 
 instance Linear Matrix Double where
     scale x = liftMatrix (scale x)
@@ -64,6 +68,8 @@ instance Linear Matrix Double where
     sub = liftMatrix2 sub
     mul = liftMatrix2 mul
     divide = liftMatrix2 divide
+    equal a b = cols a == cols b && cdat a `equal` cdat b
+
 
 instance Linear Matrix (Complex Double) where
     scale x = liftMatrix (scale x)
@@ -73,6 +79,7 @@ instance Linear Matrix (Complex Double) where
     sub = liftMatrix2 sub
     mul = liftMatrix2 mul
     divide = liftMatrix2 divide
+    equal a b = cols a == cols b && cdat a `equal` cdat b
 
 --------------------------------------------------