Her, LU, LDL, Linear, Additive

11 files changed, 239 insertions, 133 deletions

diff --git a/packages/base/src/Internal/Algorithms.hs b/packages/base/src/Internal/Algorithms.hs
index 3d25491..d2f17f4 100644
--- a/packages/base/src/Internal/Algorithms.hs
+++ b/packages/base/src/Internal/Algorithms.hs
@@ -4,6 +4,12 @@
 {-# LANGUAGE UndecidableInstances #-}
 {-# LANGUAGE TypeFamilies #-}
 
+{-# LANGUAGE FlexibleContexts, FlexibleInstances #-}
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE UndecidableInstances #-}
+{-# LANGUAGE TypeFamilies #-}
+
 -----------------------------------------------------------------------------
 {- |
 Module      :  Internal.Algorithms
@@ -32,6 +38,7 @@ import Data.List(foldl1')
 import qualified Data.Array as A
 import Internal.ST
 import Internal.Vectorized(range)
+import Control.DeepSeq
 
 {- | Generic linear algebra functions for double precision real and complex matrices.
 
@@ -43,6 +50,10 @@ class (Numeric t,
        Normed Matrix t,
        Normed Vector t,
        Floating t,
+       Linear t Vector,
+       Linear t Matrix,
+       Additive (Vector t),
+       Additive (Matrix 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)
@@ -306,25 +317,38 @@ leftSV m  | vertical m = let (u,s,_) = svd m     in (u,s)
 
 --------------------------------------------------------------
 
+-- | LU decomposition of a matrix in a compact format.
+data LU t = LU (Matrix t) [Int] deriving Show
+
+instance (NFData t, Numeric t) => NFData (LU t)
+  where
+    rnf (LU m _) = rnf m
+
 -- | 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'
+luPacked :: Field t => Matrix t -> LU t
+luPacked x = {-# SCC "luPacked" #-} LU m p
+  where
+    (m,p) = luPacked' x
 
 -- | 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'
+luSolve :: Field t => LU t -> Matrix t -> Matrix t
+luSolve (LU m p) = {-# SCC "luSolve" #-} luSolve' (m,p)
 
 -- | 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'. 
+-- | 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
+    :: Field t
+    => Matrix t -- ^ Cholesky decomposition of the coefficient matrix
+    -> Matrix t -- ^ right hand sides
+    -> Matrix t -- ^ solution
 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.
@@ -338,20 +362,28 @@ linearSolveLS = {-# SCC "linearSolveLS" #-} linearSolveLS'
 
 --------------------------------------------------------------------------------
 
+-- | LDL decomposition of a complex Hermitian or real symmetric matrix in a compact format.
+data LDL t = LDL (Matrix t) [Int] deriving Show
+
+instance (NFData t, Numeric t) => NFData (LDL t)
+  where
+    rnf (LDL m _) = rnf m
+
 -- | Similar to 'ldlPacked', without checking that the input matrix is hermitian or symmetric. It works with the lower triangular part.
-ldlPackedSH :: Field t => Matrix t -> (Matrix t, [Int])
-ldlPackedSH = {-# SCC "ldlPacked" #-} ldlPacked'
+ldlPackedSH :: Field t => Matrix t -> LDL t
+ldlPackedSH x = {-# SCC "ldlPacked" #-} LDL m p
+  where
+   (m,p) = ldlPacked' x
 
 -- | Obtains the LDL decomposition of a matrix in a compact data structure suitable for 'ldlSolve'.
-ldlPacked :: Field t => Matrix t -> (Matrix t, [Int])
-ldlPacked m
-    | exactHermitian m = {-# SCC "ldlPacked" #-} ldlPackedSH m
-    | otherwise = error "ldlPacked requires complex Hermitian or real symmetrix matrix"
-
+ldlPacked :: Field t => Her t -> LDL t
+ldlPacked (Her m) = ldlPackedSH m
 
--- | Solution of a linear system (for several right hand sides) from the precomputed LDL factorization obtained by 'ldlPacked'.
-ldlSolve :: Field t => (Matrix t, [Int]) -> Matrix t -> Matrix t
-ldlSolve = {-# SCC "ldlSolve" #-} ldlSolve'
+-- | Solution of a linear system (for several right hand sides) from a precomputed LDL factorization obtained by 'ldlPacked'.
+--
+--   Note: this can be slower than the general solver based on the LU decomposition.
+ldlSolve :: Field t => LDL t -> Matrix t -> Matrix t
+ldlSolve (LDL m p) = {-# SCC "ldlSolve" #-} ldlSolve' (m,p)
 
 --------------------------------------------------------------
 
@@ -429,14 +461,12 @@ fromList [11.344814282762075,0.17091518882717918,-0.5157294715892575]
 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"
+eigSH :: Field t => Her t -> (Vector Double, Matrix t)
+eigSH (Her m) = eigSH' m
 
 -- | 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"
+eigenvaluesSH :: Field t => Her t -> Vector Double
+eigenvaluesSH (Her m) = eigenvaluesSH' m
 
 --------------------------------------------------------------
 
@@ -490,14 +520,18 @@ 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'
+cholSH = cholSH'
 
 -- | Cholesky factorization of a positive definite hermitian or symmetric matrix.
 --
 -- If @c = chol m@ then @c@ is upper triangular and @m == tr 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"
+chol :: Field t => Her t ->  Matrix t
+chol (Her m) = {-# SCC "chol" #-} cholSH' m
+
+-- | Similar to 'chol', but instead of an error (e.g., caused by a matrix not positive definite) it returns 'Nothing'.
+mbChol :: Field t => Her t -> Maybe (Matrix t)
+mbChol (Her m) = {-# SCC "mbChol" #-} mbCholSH' m
+
 
 
 -- | Joint computation of inverse and logarithm of determinant of a square matrix.
@@ -507,7 +541,7 @@ invlndet :: Field t
 invlndet m | square m = (im,(ladm,sdm))
            | otherwise = error $ "invlndet of nonsquare "++ shSize m ++ " matrix"
   where
-    lp@(lup,perm) = luPacked m
+    lp@(LU lup perm) = luPacked m
     s = signlp (rows m) perm
     dg = toList $ takeDiag $ lup
     ladm = sum $ map (log.abs) dg
@@ -519,8 +553,9 @@ invlndet m | square m = (im,(ladm,sdm))
 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
+    where
+      LU lup perm = luPacked m
+      s = signlp (rows m) perm
 
 -- | Explicit LU factorization of a general matrix.
 --
@@ -720,7 +755,7 @@ diagonalize m = if rank v == n
                     else Nothing
     where n = rows m
           (l,v) = if exactHermitian m
-                    then let (l',v') = eigSH m in (real l', v')
+                    then let (l',v') = eigSH (trustSym m) in (real l', v')
                     else eig m
 
 -- | Generic matrix functions for diagonalizable matrices. For instance:
@@ -835,8 +870,9 @@ fixPerm' s = res $ mutable f s0
 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)
+luFact (LU l_u perm)
+    | r <= c    = (l ,u ,p, s)
+    | otherwise = (l',u',p, s)
   where
     r = rows l_u
     c = cols l_u
@@ -929,7 +965,13 @@ relativeError norm a b = r
 ----------------------------------------------------------------------
 
 -- | Generalized symmetric positive definite eigensystem Av = lBv,
--- for A and B symmetric, B positive definite (conditions not checked).
+-- for A and B symmetric, B positive definite.
+geigSH :: Field t
+        => Her t -- ^ A
+        -> Her t -- ^ B
+        -> (Vector Double, Matrix t)
+geigSH (Her a) (Her b) = geigSH' a b
+
 geigSH' :: Field t
         => Matrix t -- ^ A
         -> Matrix t -- ^ B
@@ -943,3 +985,33 @@ geigSH' a b = (l,v')
     v' = iu <> v
     (<>) = mXm
 
+--------------------------------------------------------------------------------
+
+-- | A matrix that, by construction, it is known to be complex Hermitian or real symmetric.
+--
+--   It can be created using 'sym', 'xTx', or 'trustSym', and the matrix can be extracted using 'her'.
+data Her t = Her (Matrix t) deriving Show
+
+-- | Extract the general matrix from a 'Her' structure, forgetting its symmetric or Hermitian property.
+her :: Her t -> Matrix t
+her (Her x) = x
+
+-- | Compute the complex Hermitian or real symmetric part of a square matrix (@(x + tr x)/2@).
+sym :: Field t => Matrix t -> Her t
+sym x = Her (scale 0.5 (tr x `add` x))
+
+-- | Compute the contraction @tr x <> x@ of a general matrix.
+xTx :: Numeric t => Matrix t -> Her t
+xTx x = Her (tr x `mXm` x)
+
+instance Field t => Linear t Her where
+    scale  x (Her m) = Her (scale x m)
+
+instance Field t => Additive (Her t) where
+    add    (Her a) (Her b) = Her (a `add` b)
+
+-- | At your own risk, declare that a matrix is complex Hermitian or real symmetric
+--   for usage in 'chol', 'eigSH', etc. Only a triangular part of the matrix will be used.
+trustSym :: Matrix t -> Her t
+trustSym x = (Her x)
+
diff --git a/packages/base/src/Internal/CG.hs b/packages/base/src/Internal/CG.hs
index f0142cd..cc10ad8 100644
--- a/packages/base/src/Internal/CG.hs
+++ b/packages/base/src/Internal/CG.hs
@@ -32,11 +32,11 @@ v /// b = debugMat b 2 asRow v
 type V = Vector R
 
 data CGState = CGState
-    { cgp  :: V  -- ^ conjugate gradient
-    , cgr  :: V  -- ^ residual
-    , cgr2 :: R  -- ^ squared norm of residual
-    , cgx  :: V  -- ^ current solution
-    , cgdx :: R  -- ^ normalized size of correction
+    { cgp  :: Vector R  -- ^ conjugate gradient
+    , cgr  :: Vector R  -- ^ residual
+    , cgr2 :: R         -- ^ squared norm of residual
+    , cgx  :: Vector R  -- ^ current solution
+    , cgdx :: R         -- ^ normalized size of correction
     }
 
 cg :: Bool -> (V -> V) -> (V -> V) -> CGState -> CGState
@@ -89,23 +89,25 @@ takeUntil q xs = a++ take 1 b
   where
     (a,b) = break q xs
 
+-- | Solve a sparse linear system using the conjugate gradient method with default parameters.
 cgSolve
   :: Bool          -- ^ is symmetric
   -> GMatrix       -- ^ coefficient matrix
-  -> Vector Double -- ^ right-hand side
-  -> Vector Double -- ^ solution
+  -> Vector R      -- ^ right-hand side
+  -> Vector R      -- ^ solution
 cgSolve sym a b  = cgx $ last $ cgSolve' sym 1E-4 1E-3 n a b 0
   where
     n = max 10 (round $ sqrt (fromIntegral (dim b) :: Double))
 
+-- | Solve a sparse linear system using the conjugate gradient method with default parameters.
 cgSolve'
   :: Bool      -- ^ symmetric
   -> R         -- ^ relative tolerance for the residual (e.g. 1E-4)
   -> R         -- ^ relative tolerance for δx (e.g. 1E-3)
   -> Int       -- ^ maximum number of iterations
   -> GMatrix   -- ^ coefficient matrix
-  -> V         -- ^ initial solution
-  -> V         -- ^ right-hand side
+  -> Vector R  -- ^ initial solution
+  -> Vector R  -- ^ right-hand side
   -> [CGState] -- ^ solution
 cgSolve' sym er es n a b x = take n $ conjugrad sym a b x er es
 
diff --git a/packages/base/src/Internal/Modular.hs b/packages/base/src/Internal/Modular.hs
index 64ed2bb..a3421a8 100644
--- a/packages/base/src/Internal/Modular.hs
+++ b/packages/base/src/Internal/Modular.hs
@@ -33,7 +33,7 @@ import Internal.Element
 import Internal.Container
 import Internal.Vectorized (prodI,sumI,prodL,sumL)
 import Internal.LAPACK (multiplyI, multiplyL)
-import Internal.Algorithms(luFact)
+import Internal.Algorithms(luFact,LU(..))
 import Internal.Util(Normed(..),Indexable(..),
                      gaussElim, gaussElim_1, gaussElim_2,
                      luST, luSolve', luPacked', magnit, invershur)
@@ -169,7 +169,7 @@ instance forall m . KnownNat m => Container Vector (Mod m I)
     size' = dim
     scale' s x = vmod (scale (unMod s) (f2i x))
     addConstant c x = vmod (addConstant (unMod c) (f2i x))
-    add a b = vmod (add (f2i a) (f2i b))
+    add' a b = vmod (add' (f2i a) (f2i b))
     sub a b = vmod (sub (f2i a) (f2i b))
     mul a b = vmod (mul (f2i a) (f2i b))
     equal u v = equal (f2i u) (f2i v)
@@ -209,7 +209,7 @@ instance forall m . KnownNat m => Container Vector (Mod m Z)
     size' = dim
     scale' s x = vmod (scale (unMod s) (f2i x))
     addConstant c x = vmod (addConstant (unMod c) (f2i x))
-    add a b = vmod (add (f2i a) (f2i b))
+    add' a b = vmod (add' (f2i a) (f2i b))
     sub a b = vmod (sub (f2i a) (f2i b))
     mul a b = vmod (mul (f2i a) (f2i b))
     equal u v = equal (f2i u) (f2i v)
@@ -371,7 +371,9 @@ test = (ok, info)
 
     checkLU okf t = norm_Inf $ flatten (l <> u <> p - t)
       where
-        (l,u,p,_ :: Int) = luFact $ mutable (luST okf) t
+        (l,u,p,_ :: Int) = luFact (LU x' p')
+          where
+            (x',p') = mutable (luST okf) t
 
     checkSolve aa = norm_Inf $ flatten (aa <> x - bb)
        where
diff --git a/packages/base/src/Internal/Numeric.hs b/packages/base/src/Internal/Numeric.hs
index a8ae2bb..e8c7440 100644
--- a/packages/base/src/Internal/Numeric.hs
+++ b/packages/base/src/Internal/Numeric.hs
@@ -49,7 +49,7 @@ class Element e => Container c e
     scalar'      :: e -> c e
     scale'       :: e -> c e -> c e
     addConstant :: e -> c e -> c e
-    add         :: c e -> c e -> c e
+    add'        :: c e -> c e -> c e
     sub         :: c e -> c e -> c e
     -- | element by element multiplication
     mul         :: c e -> c e -> c e
@@ -100,7 +100,7 @@ instance Container Vector I
     size' = dim
     scale' = vectorMapValI Scale
     addConstant = vectorMapValI AddConstant
-    add = vectorZipI Add
+    add' = vectorZipI Add
     sub = vectorZipI Sub
     mul = vectorZipI Mul
     equal u v = dim u == dim v && maxElement' (vectorMapI Abs (sub u v)) == 0
@@ -139,7 +139,7 @@ instance Container Vector Z
     size' = dim
     scale' = vectorMapValL Scale
     addConstant = vectorMapValL AddConstant
-    add = vectorZipL Add
+    add' = vectorZipL Add
     sub = vectorZipL Sub
     mul = vectorZipL Mul
     equal u v = dim u == dim v && maxElement' (vectorMapL Abs (sub u v)) == 0
@@ -179,7 +179,7 @@ instance Container Vector Float
     size' = dim
     scale' = vectorMapValF Scale
     addConstant = vectorMapValF AddConstant
-    add = vectorZipF Add
+    add' = vectorZipF Add
     sub = vectorZipF Sub
     mul = vectorZipF Mul
     equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0
@@ -216,7 +216,7 @@ instance Container Vector Double
     size' = dim
     scale' = vectorMapValR Scale
     addConstant = vectorMapValR AddConstant
-    add = vectorZipR Add
+    add' = vectorZipR Add
     sub = vectorZipR Sub
     mul = vectorZipR Mul
     equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0
@@ -253,7 +253,7 @@ instance Container Vector (Complex Double)
     size' = dim
     scale' = vectorMapValC Scale
     addConstant = vectorMapValC AddConstant
-    add = vectorZipC Add
+    add' = vectorZipC Add
     sub = vectorZipC Sub
     mul = vectorZipC Mul
     equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
@@ -289,7 +289,7 @@ instance Container Vector (Complex Float)
     size' = dim
     scale' = vectorMapValQ Scale
     addConstant = vectorMapValQ AddConstant
-    add = vectorZipQ Add
+    add' = vectorZipQ Add
     sub = vectorZipQ Sub
     mul = vectorZipQ Mul
     equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
@@ -327,7 +327,7 @@ instance (Num a, Element a, Container Vector a) => Container Matrix a
     size' = size
     scale' x = liftMatrix (scale' x)
     addConstant x = liftMatrix (addConstant x)
-    add = liftMatrix2 add
+    add' = liftMatrix2 add'
     sub = liftMatrix2 sub
     mul = liftMatrix2 mul
     equal a b = cols a == cols b && flatten a `equal` flatten b
@@ -387,9 +387,6 @@ scalar = scalar'
 conj :: Container c e => c e -> c e
 conj = conj'
 
--- | multiplication by scalar
-scale :: Container c e => e -> c e -> c e
-scale = scale'
 
 arctan2 :: (Fractional e, Container c e) => c e -> c e -> c e
 arctan2 = arctan2'
@@ -581,6 +578,10 @@ class ( Container Vector t
       , Konst t (Int,Int) Matrix
       , CTrans t
       , Product t
+      , Additive (Vector t)
+      , Additive (Matrix t)
+      , Linear t Vector
+      , Linear t Matrix
       ) => Numeric t
 
 instance Numeric Double
@@ -912,11 +913,30 @@ instance (CTrans t, Container Vector t) => Transposable (Matrix t) (Matrix t)
     tr  = ctrans
     tr' = trans
 
-class Linear t v
+class Additive c
   where
-    scalarL :: t -> v
-    addL    :: v -> v -> v
-    scaleL  :: t -> v -> v
+    add    :: c -> c -> c
+
+class Linear t c
+  where
+    scale  :: t -> c t -> c t
+
+
+instance Container Vector t => Linear t Vector
+  where
+    scale = scale'
+
+instance Container Matrix t => Linear t Matrix
+  where
+    scale = scale'
+
+instance Container Vector t => Additive (Vector t)
+  where
+    add = add'
+
+instance Container Matrix t => Additive (Matrix t)
+  where
+    add = add'
 
 
 class Testable t
diff --git a/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs
index 4123e6c..36b7855 100644
--- a/packages/base/src/Internal/Util.hs
+++ b/packages/base/src/Internal/Util.hs
@@ -458,12 +458,12 @@ rowOuters a b = a' * b'
 --------------------------------------------------------------------------------
 
 -- | solution of overconstrained homogeneous linear system
-null1 :: Matrix Double -> Vector Double
+null1 :: Matrix R -> Vector R
 null1 = last . toColumns . snd . rightSV
 
 -- | solution of overconstrained homogeneous symmetric linear system
-null1sym :: Matrix Double -> Vector Double
-null1sym = last . toColumns . snd . eigSH'
+null1sym :: Her R -> Vector R
+null1sym = last . toColumns . snd . eigSH
 
 --------------------------------------------------------------------------------
 
@@ -712,7 +712,9 @@ luST ok (r,_) x = do
  ,  0,  0,  0,  0,  1 ]
 
 -}
-luPacked' x = mutable (luST (magnit 0)) x
+luPacked' x = LU m p
+  where
+    (m,p) = mutable (luST (magnit 0)) x
 
 --------------------------------------------------------------------------------
 
@@ -782,7 +784,7 @@ forwSust' lup rhs = foldl' f (rhs?[]) ls
                 (b - l<>x)
 
 
-luSolve'' (lup,p) b = backSust' lup (forwSust' lup pb)
+luSolve'' (LU lup p) b = backSust' lup (forwSust' lup pb)
   where
     pb = b ?? (Pos (fixPerm' p), All)
 
@@ -827,7 +829,7 @@ backSust lup rhs = fst $ mutable f rhs
  , 7, 10, 6 ]
 
 -}
-luSolve' (lup,p) b = backSust lup (forwSust lup pb)
+luSolve' (LU lup p) b = backSust lup (forwSust lup pb)
   where
     pb = b ?? (Pos (fixPerm' p), All)
 
diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs
index 9a924e0..7be2600 100644
--- a/packages/base/src/Numeric/LinearAlgebra.hs
+++ b/packages/base/src/Numeric/LinearAlgebra.hs
@@ -53,11 +53,11 @@ module Numeric.LinearAlgebra (
     --
 
     -- * Products
-    -- ** dot
+    -- ** Dot
     dot, (<.>),
-    -- ** matrix-vector
+    -- ** Matrix-vector
     (#>), (<#), (!#>),
-    -- ** matrix-matrix
+    -- ** Matrix-matrix
     (<>),
     -- | The matrix product is also implemented in the "Data.Monoid" instance, where
     -- single-element matrices (created from numeric literals or using 'scalar')
@@ -73,20 +73,25 @@ module Numeric.LinearAlgebra (
     -- 'mconcat' uses 'optimiseMult' to get the optimal association order.
 
 
-    -- ** other
+    -- ** Other
     outer, kronecker, cross,
-    scale,
+    scale, add,
     sumElements, prodElements,
 
     -- * Linear systems
+    -- ** General
     (<\>),
-    linearSolve,
     linearSolveLS,
     linearSolveSVD,
-    luSolve,
-    luSolve',
+    -- ** Determined
+    linearSolve,
+    luSolve, luPacked,
+    luSolve', luPacked',
+    -- ** Symmetric indefinite
+    ldlSolve, ldlPacked,
+    -- ** Positive definite
     cholSolve,
-    ldlSolve,
+    -- ** Sparse
     cgSolve,
     cgSolve',
 
@@ -113,21 +118,18 @@ module Numeric.LinearAlgebra (
     leftSV, rightSV,
 
     -- * Eigendecomposition
-    eig, eigSH, eigSH',
-    eigenvalues, eigenvaluesSH, eigenvaluesSH',
-    geigSH',
+    eig, eigSH,
+    eigenvalues, eigenvaluesSH,
+    geigSH,
 
     -- * QR
     qr, rq, qrRaw, qrgr,
 
     -- * Cholesky
-    chol, cholSH, mbCholSH,
+    chol, mbChol,
 
     -- * LU
-    lu, luPacked, luPacked', luFact,
-
-    -- * LDL
-    ldlPacked, ldlPackedSH,
+    lu, luFact,
 
     -- * Hessenberg
     hess,
@@ -150,14 +152,16 @@ module Numeric.LinearAlgebra (
     -- * Misc
     meanCov, rowOuters, pairwiseD2, unitary, peps, relativeError, magnit,
     haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv,
-    iC,
+    iC, sym, xTx, trustSym, her,
     -- * Auxiliary classes
-    Element, Container, Product, Numeric, LSDiv,
+    Element, Container, Product, Numeric, LSDiv, Her,
     Complexable, RealElement,
     RealOf, ComplexOf, SingleOf, DoubleOf,
     IndexOf,
-    Field,
+    Field, Linear(), Additive(),
     Transposable,
+    LU(..),
+    LDL(..),
     CGState(..),
     Testable(..)
 ) where
@@ -169,7 +173,7 @@ import Numeric.Vector()
 import Internal.Matrix
 import Internal.Container hiding ((<>))
 import Internal.Numeric hiding (mul)
-import Internal.Algorithms hiding (linearSolve,Normed,orth,luPacked',linearSolve',luSolve')
+import Internal.Algorithms hiding (linearSolve,Normed,orth,luPacked',linearSolve',luSolve',ldlPacked')
 import qualified Internal.Algorithms as A
 import Internal.Util
 import Internal.Random
@@ -246,4 +250,3 @@ nullspace m = nullspaceSVD (Left (1*eps)) m (rightSV m)
 -- | return an orthonormal basis of the range space of a matrix. See also 'orthSVD'.
 orth m = orthSVD (Left (1*eps)) m (leftSV m)
 
-
diff --git a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs
index bac1c0c..5ce529c 100644
--- a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs
@@ -13,11 +13,12 @@ compatibility with previous version, to be removed
 
 module Numeric.LinearAlgebra.HMatrix (
     module Numeric.LinearAlgebra,
-    (¦),(——),ℝ,ℂ,(<·>),app,mul
+    (¦),(——),ℝ,ℂ,(<·>),app,mul, cholSH, mbCholSH, eigSH', eigenvaluesSH', geigSH'
 ) where
 
 import Numeric.LinearAlgebra
 import Internal.Util
+import Internal.Algorithms(cholSH, mbCholSH, eigSH', eigenvaluesSH', geigSH')
 
 infixr 8 <·>
 (<·>) :: Numeric t => Vector t -> Vector t -> t
diff --git a/packages/base/src/Numeric/LinearAlgebra/Static.hs b/packages/base/src/Numeric/LinearAlgebra/Static.hs
index 0dab0e6..ded69fa 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Static.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Static.hs
@@ -63,9 +63,9 @@ import GHC.TypeLits
 import Numeric.LinearAlgebra hiding (
     (<>),(#>),(<.>),Konst(..),diag, disp,(===),(|||),
     row,col,vector,matrix,linspace,toRows,toColumns,
-    (<\>),fromList,takeDiag,svd,eig,eigSH,eigSH',
-    eigenvalues,eigenvaluesSH,eigenvaluesSH',build,
-    qr,size,dot,chol,range,R,C)
+    (<\>),fromList,takeDiag,svd,eig,eigSH,
+    eigenvalues,eigenvaluesSH,build,
+    qr,size,dot,chol,range,R,C,Her,her,sym)
 import qualified Numeric.LinearAlgebra as LA
 import Data.Proxy(Proxy)
 import Internal.Static
@@ -292,10 +292,10 @@ her m = Her $ (m + LA.tr m)/2
 
 instance KnownNat n => Eigen (Sym n) (R n) (L n n)
   where
-    eigenvalues (Sym (extract -> m)) =  mkR . LA.eigenvaluesSH' $ m
+    eigenvalues (Sym (extract -> m)) =  mkR . LA.eigenvaluesSH . LA.trustSym $ m
     eigensystem (Sym (extract -> m)) = (mkR l, mkL v)
       where
-        (l,v) = LA.eigSH' m
+        (l,v) = LA.eigSH . LA.trustSym $ m
 
 instance KnownNat n => Eigen (Sq n) (C n) (M n n)
   where
@@ -305,7 +305,7 @@ instance KnownNat n => Eigen (Sq n) (C n) (M n n)
         (l,v) = LA.eig m
 
 chol :: KnownNat n => Sym n -> Sq n
-chol (extract . unSym -> m) = mkL $ LA.cholSH m
+chol (extract . unSym -> m) = mkL $ LA.chol $ LA.trustSym m
 
 --------------------------------------------------------------------------------
 
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
index 2ff1580..30480d7 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
@@ -127,8 +127,8 @@ expmTest2 = expm nd2 :~15~: (2><2)
 mbCholTest = utest "mbCholTest" (ok1 && ok2) where
     m1 = (2><2) [2,5,5,8 :: Double]
     m2 = (2><2) [3,5,5,9 :: Complex Double]
-    ok1 = mbCholSH m1 == Nothing
-    ok2 = mbCholSH m2 == Just (chol m2)
+    ok1 = mbChol (trustSym m1) == Nothing
+    ok2 = mbChol (trustSym m2) == Just (chol $ trustSym m2)
 
 ---------------------------------------------------------------------
 
@@ -403,8 +403,8 @@ indexProp g f x = a1 == g a2 && a2 == a3 && b1 == g b2 && b2 == b3
 --------------------------------------------------------------------------------
 
 sliceTest = utest "slice test" $ and
-    [ testSlice chol  (gen 5 :: Matrix R)
-    , testSlice chol  (gen 5 :: Matrix C)
+    [ testSlice (chol . trustSym)  (gen 5 :: Matrix R)
+    , testSlice (chol . trustSym)  (gen 5 :: Matrix C)
     , testSlice qr    (rec :: Matrix R)
     , testSlice qr    (rec :: Matrix C)
     , testSlice hess  (agen 5 :: Matrix R)
@@ -420,12 +420,12 @@ sliceTest = utest "slice test" $ and
 
     , testSlice eig   (agen 5 :: Matrix R)
     , testSlice eig   (agen 5 :: Matrix C)
-    , testSlice eigSH (gen 5 :: Matrix R)
-    , testSlice eigSH (gen 5 :: Matrix C)
+    , testSlice (eigSH . trustSym) (gen 5 :: Matrix R)
+    , testSlice (eigSH . trustSym) (gen 5 :: Matrix C)
     , testSlice eigenvalues   (agen 5 :: Matrix R)
     , testSlice eigenvalues   (agen 5 :: Matrix C)
-    , testSlice eigenvaluesSH (gen 5 :: Matrix R)
-    , testSlice eigenvaluesSH (gen 5 :: Matrix C)
+    , testSlice (eigenvaluesSH . trustSym) (gen 5 :: Matrix R)
+    , testSlice (eigenvaluesSH . trustSym) (gen 5 :: Matrix C)
 
     , testSlice svd           (rec :: Matrix R)
     , testSlice thinSVD       (rec :: Matrix R)
@@ -489,10 +489,10 @@ sliceTest = utest "slice test" $ and
     , testSlice ((<>) (ogen 5:: Matrix (Z ./. 7))) (gen 5)
     , testSlice (flip (<>) (gen 5:: Matrix (Z ./. 7))) (ogen 5)
 
-    , testSlice (flip cholSolve (agen 5:: Matrix R)) (chol $ gen 5)
-    , testSlice (flip cholSolve (agen 5:: Matrix C)) (chol $ gen 5)
-    , testSlice (cholSolve (chol $ gen 5:: Matrix R)) (agen 5)
-    , testSlice (cholSolve (chol $ gen 5:: Matrix C)) (agen 5)
+    , testSlice (flip cholSolve (agen 5:: Matrix R)) (chol $ trustSym $ gen 5)
+    , testSlice (flip cholSolve (agen 5:: Matrix C)) (chol $ trustSym $ gen 5)
+    , testSlice (cholSolve (chol $ trustSym $ gen 5:: Matrix R)) (agen 5)
+    , testSlice (cholSolve (chol $ trustSym $ gen 5:: Matrix C)) (agen 5)
 
     , ok_qrgr        (rec :: Matrix R)
     , ok_qrgr        (rec :: Matrix C)
@@ -515,8 +515,8 @@ sliceTest = utest "slice test" $ and
 
     test_lus m = testSlice f lup
       where
-        f x = luSolve (x,p) m
-        (lup,p) = luPacked m
+        f x = luSolve (LU x p) m
+        (LU lup p) = luPacked m
 
     gen :: Numeric t => Int -> Matrix t
     gen n = diagRect 1 (konst 5 n) n n
@@ -588,11 +588,11 @@ runTests n = do
     test (linearSolveProp (luSolve.luPacked) . rSqWC)
     test (linearSolveProp (luSolve.luPacked) . cSqWC)
     putStrLn "------ ldlSolve"
-    test (linearSolveProp (ldlSolve.ldlPacked) . rSymWC)
-    test (linearSolveProp (ldlSolve.ldlPacked) . cSymWC)
+    test (linearSolvePropH (ldlSolve.ldlPacked) . rSymWC)
+    test (linearSolvePropH (ldlSolve.ldlPacked) . cSymWC)
     putStrLn "------ cholSolve"
-    test (linearSolveProp (cholSolve.chol) . rPosDef)
-    test (linearSolveProp (cholSolve.chol) . cPosDef)
+    test (linearSolveProp (cholSolve.chol.trustSym) . rPosDef)
+    test (linearSolveProp (cholSolve.chol.trustSym) . cPosDef)
     putStrLn "------ luSolveLS"
     test (linearSolveProp linearSolveLS . rSqWC)
     test (linearSolveProp linearSolveLS . cSqWC)
@@ -865,8 +865,8 @@ eigBench = do
     let m = reshape 1000 (randomVector 777 Uniform (1000*1000))
         s = m + tr m
     m `seq` s `seq` putStrLn ""
-    time "eigenvalues  symmetric 1000x1000" (eigenvaluesSH' m)
-    time "eigenvectors symmetric 1000x1000" (snd $ eigSH' m)
+    time "eigenvalues  symmetric 1000x1000" (eigenvaluesSH (trustSym m))
+    time "eigenvectors symmetric 1000x1000" (snd $ eigSH (trustSym m))
     time "eigenvalues  general   1000x1000" (eigenvalues m)
     time "eigenvectors general   1000x1000" (snd $ eig m)
 
@@ -893,12 +893,14 @@ solveBenchN n = do
     time ("svd solve " ++ show n) (linearSolveSVD a b)
     time (" ls solve " ++ show n) (linearSolveLS a b)
     time ("    solve " ++ show n) (linearSolve a b)
-    time ("cholSolve " ++ show n) (cholSolve (chol a) b)
+--    time (" LU solve " ++ show n) (luSolve (luPacked a) b)
+    time ("LDL solve " ++ show n) (ldlSolve (ldlPacked (trustSym a)) b)
+    time ("cholSolve " ++ show n) (cholSolve (chol $ trustSym a) b)
 
 solveBench = do
     solveBenchN 500
     solveBenchN 1000
-    -- solveBenchN 1500
+    solveBenchN 1500
 
 --------------------------------
 
@@ -906,7 +908,7 @@ cholBenchN n = do
     let x = uniformSample 777 (2*n) (replicate n (-1,1))
         a = tr x <> x
     a `seq` putStr ""
-    time ("chol " ++ show n) (chol a)
+    time ("chol " ++ show n) (chol $ trustSym a)
 
 cholBench = do
     putStrLn ""
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
index 7c54535..4704989 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
@@ -14,7 +14,7 @@ Arbitrary instances for vectors, matrices.
 module Numeric.LinearAlgebra.Tests.Instances(
     Sq(..),     rSq,cSq,
     Rot(..),    rRot,cRot,
-    Her(..),    rHer,cHer,
+                rHer,cHer,
     WC(..),     rWC,cWC,
     SqWC(..),   rSqWC, cSqWC, rSymWC, cSymWC,
     PosDef(..), rPosDef, cPosDef,
@@ -81,12 +81,12 @@ instance (Field a, Arbitrary a) => Arbitrary (Rot a) where
 
 
 -- a complex hermitian or real symmetric matrix
-newtype (Her a) = Her (Matrix a) deriving Show
+--newtype (Her a) = Her (Matrix a) deriving Show
 instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Her a) where
     arbitrary = do
         Sq m <- arbitrary
         let m' = m/2
-        return $ Her (m' + tr m')
+        return $ sym m'
 
 
 class (Field a, Arbitrary a, Element (RealOf a), Random (RealOf a)) => ArbitraryField a
@@ -125,9 +125,9 @@ newtype (PosDef a) = PosDef (Matrix a) deriving Show
 instance (Numeric a, ArbitraryField a, Num (Vector a)) 
     => Arbitrary (PosDef a) where
     arbitrary = do
-        Her m <- arbitrary
+        m <- arbitrary
         let (_,v) = eigSH m
-            n = rows m
+            n = rows (her m)
         l <- replicateM n (choose (0,100))
         let s = diag (fromList l)
             p = v <> real s <> tr v
@@ -161,8 +161,8 @@ fM m = m :: FM
 zM m = m :: ZM
 
 
-rHer (Her m) = m :: RM
-cHer (Her m) = m :: CM
+rHer m = her m :: RM
+cHer m = her m :: CM
 
 rRot (Rot m) = m :: RM
 cRot (Rot m) = m :: CM
@@ -176,8 +176,8 @@ cWC (WC m) = m :: CM
 rSqWC (SqWC m) = m :: RM
 cSqWC (SqWC m) = m :: CM
 
-rSymWC (SqWC m) = m + tr m :: RM
-cSymWC (SqWC m) = m + tr m :: CM
+rSymWC (SqWC m) = sym m :: Her R
+cSymWC (SqWC m) = sym m :: Her C
 
 rPosDef (PosDef m) = m :: RM
 cPosDef (PosDef m) = m :: CM
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
index 207a303..2ac3588 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -39,7 +39,7 @@ module Numeric.LinearAlgebra.Tests.Properties (
     expmDiagProp,
     multProp1, multProp2,
     subProp,
-    linearSolveProp, linearSolveProp2
+    linearSolveProp, linearSolvePropH, linearSolveProp2
 ) where
 
 import Numeric.LinearAlgebra.HMatrix hiding (Testable,unitary)
@@ -209,11 +209,11 @@ eigProp m = complex m <> v |~| v <> diag s
 eigSHProp m = m <> v |~| v <> real (diag s)
               && unitary v
               && m |~| v <> real (diag s) <> tr v
-    where (s, v) = eigSH m
+    where (s, v) = eigSH' m
 
 eigProp2 m = fst (eig m) |~| eigenvalues m
 
-eigSHProp2 m = fst (eigSH m) |~| eigenvaluesSH m
+eigSHProp2 m = fst (eigSH' m) |~| eigenvaluesSH' m
 
 ------------------------------------------------------------------
 
@@ -246,9 +246,9 @@ schurProp2 m = m |~| u <> s <> tr u && unitary u && upperHessenberg s -- fixme
     where (u,s) = schur m
 
 cholProp m = m |~| tr c <> c && upperTriang c
-    where c = chol m
+    where c = chol (trustSym m)
 
-exactProp m = chol m == chol (m+0)
+exactProp m = chol (trustSym m) == chol (trustSym (m+0))
 
 expmDiagProp m = expm (logm m) :~ 7 ~: complex m
     where logm = matFunc log
@@ -263,6 +263,8 @@ multProp2 p (a,b) = (tr (a <> b)) :~p~: (tr b <> tr a)
 
 linearSolveProp f m = f m m |~| ident (rows m)
 
+linearSolvePropH f m = f m (her m) |~| ident (rows (her m))
+
 linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b)
     where q = min (rows a) (cols a)
           b = a <> x