10 files changed, 95 insertions, 72 deletions

diff --git a/packages/base/CHANGELOG b/packages/base/CHANGELOG
index 0336a28..581d2ac 100644
--- a/packages/base/CHANGELOG
+++ b/packages/base/CHANGELOG
@@ -1,24 +1,32 @@
 0.17.0.0
 --------
 
-    * improved matrix extraction (??) and rectangular matrix slices without data copy
+    * eigSH, chol, and other functions that work with Hermitian or symmetric matrices
+      now take a special "Herm" argument that can be created by means of "sym"
+      or "mTm". The unchecked versions of those functions have been removed and we
+      use "trustSym" to create the Herm type when the matrix is known to be Hermitian/symmetric.
+
+    * Improved matrix extraction (??) and rectangular matrix slices without data copy
+
+    * Basic support of Int32 and Int64 elements
 
-    * basic support of Int32 and Int64 elements
-    
     * remap, more general cond, sortIndex
 
-    * experimental support of type safe modular arithmetic, including linear
-      systems and lu factorization
-    
-    * elementary row operations and inplace matrix slice products in the ST monad
+    * Experimental support of type safe modular arithmetic, including linear
+      system solver and LU factorization
+
+    * Elementary row operations and inplace matrix slice products in the ST monad
 
     * Improved development tools.
 
-    * old compatibility modules removed, simpler organization of internal modules
+    * Old compatibility modules removed, simpler organization of internal modules
 
     * unitary, pairwiseD2, tr'
 
-    * ldlPacked, ldlSolve
+    * ldlPacked, ldlSolve for indefinite symmetric systems (apparently not faster
+      than the general solver based on the LU)
+
+    * LU, LDL, and QR types for these compact decompositions.
 
 0.16.1.0
 --------
diff --git a/packages/base/hmatrix.cabal b/packages/base/hmatrix.cabal
index 93ca3d7..e248886 100644
--- a/packages/base/hmatrix.cabal
+++ b/packages/base/hmatrix.cabal
@@ -9,9 +9,9 @@ Homepage:           https://github.com/albertoruiz/hmatrix
 Synopsis:           Numeric Linear Algebra
 Description:        Linear systems, matrix decompositions, and other numerical computations based on BLAS and LAPACK.
                     .
-                    The standard interface is provided by the module "Numeric.LinearAlgebra".
+                    Standard interface: "Numeric.LinearAlgebra".
                     .
-                    A safer interface with statically checked dimensions is provided by "Numeric.LinearAlgebra.Static".
+                    Safer interface with statically checked dimensions: "Numeric.LinearAlgebra.Static".
                     .
                     Code examples: <http://dis.um.es/~alberto/hmatrix/hmatrix.html>
 
diff --git a/packages/base/src/Internal/Algorithms.hs b/packages/base/src/Internal/Algorithms.hs
index d2f17f4..c4f1a60 100644
--- a/packages/base/src/Internal/Algorithms.hs
+++ b/packages/base/src/Internal/Algorithms.hs
@@ -4,12 +4,6 @@
 {-# LANGUAGE UndecidableInstances #-}
 {-# LANGUAGE TypeFamilies #-}
 
-{-# LANGUAGE FlexibleContexts, FlexibleInstances #-}
-{-# LANGUAGE CPP #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE UndecidableInstances #-}
-{-# LANGUAGE TypeFamilies #-}
-
 -----------------------------------------------------------------------------
 {- |
 Module      :  Internal.Algorithms
@@ -376,8 +370,8 @@ ldlPackedSH x = {-# SCC "ldlPacked" #-} LDL m p
    (m,p) = ldlPacked' x
 
 -- | Obtains the LDL decomposition of a matrix in a compact data structure suitable for 'ldlSolve'.
-ldlPacked :: Field t => Her t -> LDL t
-ldlPacked (Her m) = ldlPackedSH m
+ldlPacked :: Field t => Herm t -> LDL t
+ldlPacked (Herm m) = ldlPackedSH m
 
 -- | Solution of a linear system (for several right hand sides) from a precomputed LDL factorization obtained by 'ldlPacked'.
 --
@@ -461,26 +455,39 @@ fromList [11.344814282762075,0.17091518882717918,-0.5157294715892575]
 3.000  5.000  6.000
 
 -}
-eigSH :: Field t => Her t -> (Vector Double, Matrix t)
-eigSH (Her m) = eigSH' m
+eigSH :: Field t => Herm t -> (Vector Double, Matrix t)
+eigSH (Herm m) = eigSH' m
 
 -- | Eigenvalues (in descending order) of a complex hermitian or real symmetric matrix.
-eigenvaluesSH :: Field t => Her t -> Vector Double
-eigenvaluesSH (Her m) = eigenvaluesSH' m
+eigenvaluesSH :: Field t => Herm t -> Vector Double
+eigenvaluesSH (Herm m) = eigenvaluesSH' m
 
 --------------------------------------------------------------
 
+-- | QR decomposition of a matrix in compact form. (The orthogonal matrix is not explicitly formed.)
+data QR t = QR (Matrix t) (Vector t)
+
+instance (NFData t, Numeric t) => NFData (QR t)
+  where
+    rnf (QR m _) = rnf m
+
+
 -- | 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
+-- | Compute the QR decomposition of a matrix in compact form.
+qrRaw :: Field t => Matrix t -> QR t
+qrRaw m = QR x v
+  where
+    (x,v) = qr' m
 
-{- | generate a matrix with k orthogonal columns from the output of qrRaw
--}
-qrgr n (a,t)
+-- | generate a matrix with k orthogonal columns from the compact QR decomposition obtained by 'qrRaw'.
+--
+qrgr :: Field t => Int -> QR t -> Matrix t
+qrgr n (QR 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)
 
@@ -525,12 +532,12 @@ 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 => Her t ->  Matrix t
-chol (Her m) = {-# SCC "chol" #-} cholSH' m
+chol :: Field t => Herm t ->  Matrix t
+chol (Herm 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
+mbChol :: Field t => Herm t -> Maybe (Matrix t)
+mbChol (Herm m) = {-# SCC "mbChol" #-} mbCholSH' m
 
 
 
@@ -870,6 +877,8 @@ 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
 
+-- | Compute the explicit LU decomposition from the compact one obtained by 'luPacked'.
+luFact :: Numeric t => LU t -> (Matrix t, Matrix t, Matrix t, t)
 luFact (LU l_u perm)
     | r <= c    = (l ,u ,p, s)
     | otherwise = (l',u',p, s)
@@ -967,10 +976,10 @@ relativeError norm a b = r
 -- | Generalized symmetric positive definite eigensystem Av = lBv,
 -- for A and B symmetric, B positive definite.
 geigSH :: Field t
-        => Her t -- ^ A
-        -> Her t -- ^ B
+        => Herm t -- ^ A
+        -> Herm t -- ^ B
         -> (Vector Double, Matrix t)
-geigSH (Her a) (Her b) = geigSH' a b
+geigSH (Herm a) (Herm b) = geigSH' a b
 
 geigSH' :: Field t
         => Matrix t -- ^ A
@@ -989,29 +998,33 @@ geigSH' a b = (l,v')
 
 -- | 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
+--   It can be created using 'sym', 'mTm', or 'trustSym', and the matrix can be extracted using 'unSym'.
+newtype Herm t = Herm (Matrix t) deriving Show
+
+instance (NFData t, Numeric t) => NFData (Herm t)
+  where
+    rnf (Herm m) = rnf m
 
--- | Extract the general matrix from a 'Her' structure, forgetting its symmetric or Hermitian property.
-her :: Her t -> Matrix t
-her (Her x) = x
+-- | Extract the general matrix from a 'Herm' structure, forgetting its symmetric or Hermitian property.
+unSym :: Herm t -> Matrix t
+unSym (Herm 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))
+sym :: Field t => Matrix t -> Herm t
+sym x = Herm (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)
+mTm :: Numeric t => Matrix t -> Herm t
+mTm x = Herm (tr x `mXm` x)
 
-instance Field t => Linear t Her where
-    scale  x (Her m) = Her (scale x m)
+instance Field t => Linear t Herm where
+    scale  x (Herm m) = Herm (scale x m)
 
-instance Field t => Additive (Her t) where
-    add    (Her a) (Her b) = Her (a `add` b)
+instance Field t => Additive (Herm t) where
+    add (Herm a) (Herm b) = Herm (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)
+trustSym :: Matrix t -> Herm t
+trustSym x = (Herm x)
 
diff --git a/packages/base/src/Internal/Modular.hs b/packages/base/src/Internal/Modular.hs
index a3421a8..239c742 100644
--- a/packages/base/src/Internal/Modular.hs
+++ b/packages/base/src/Internal/Modular.hs
@@ -371,7 +371,7 @@ test = (ok, info)
 
     checkLU okf t = norm_Inf $ flatten (l <> u <> p - t)
       where
-        (l,u,p,_ :: Int) = luFact (LU x' p')
+        (l,u,p,_) = luFact (LU x' p')
           where
             (x',p') = mutable (luST okf) t
 
diff --git a/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs
index 36b7855..cf42961 100644
--- a/packages/base/src/Internal/Util.hs
+++ b/packages/base/src/Internal/Util.hs
@@ -462,7 +462,7 @@ null1 :: Matrix R -> Vector R
 null1 = last . toColumns . snd . rightSV
 
 -- | solution of overconstrained homogeneous symmetric linear system
-null1sym :: Her R -> Vector R
+null1sym :: Herm R -> Vector R
 null1sym = last . toColumns . snd . eigSH
 
 --------------------------------------------------------------------------------
diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs
index 7be2600..6a9c33a 100644
--- a/packages/base/src/Numeric/LinearAlgebra.hs
+++ b/packages/base/src/Numeric/LinearAlgebra.hs
@@ -35,8 +35,9 @@ module Numeric.LinearAlgebra (
 
     -- * Autoconformable dimensions
     -- |
-    -- In arithmetic operations single-element vectors and matrices
-    -- (created from numeric literals or using 'scalar') automatically
+    -- In most operations, single-element vectors and matrices
+    -- (created from numeric literals or using 'scalar'), and matrices
+    -- with just one row or column, automatically
     -- expand to match the dimensions of the other operand:
     --
     -- >>> 5 + 2*ident 3 :: Matrix Double
@@ -45,11 +46,12 @@ module Numeric.LinearAlgebra (
     --  , 5.0, 7.0, 5.0
     --  , 5.0, 5.0, 7.0 ]
     --
-    -- >>> matrix 3 [1..9] + matrix 1 [10,20,30]
-    -- (3><3)
-    --  [ 11.0, 12.0, 13.0
-    --  , 24.0, 25.0, 26.0
-    --  , 37.0, 38.0, 39.0 ]
+    -- >>> (4><3) [1..] + row [10,20,30]
+    -- (4><3)
+    --  [ 11.0, 22.0, 33.0
+    --  , 14.0, 25.0, 36.0
+    --  , 17.0, 28.0, 39.0
+    --  , 20.0, 31.0, 42.0 ]
     --
 
     -- * Products
@@ -152,9 +154,9 @@ module Numeric.LinearAlgebra (
     -- * Misc
     meanCov, rowOuters, pairwiseD2, unitary, peps, relativeError, magnit,
     haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv,
-    iC, sym, xTx, trustSym, her,
+    iC, sym, mTm, trustSym, unSym,
     -- * Auxiliary classes
-    Element, Container, Product, Numeric, LSDiv, Her,
+    Element, Container, Product, Numeric, LSDiv, Herm,
     Complexable, RealElement,
     RealOf, ComplexOf, SingleOf, DoubleOf,
     IndexOf,
@@ -162,6 +164,7 @@ module Numeric.LinearAlgebra (
     Transposable,
     LU(..),
     LDL(..),
+    QR(..),
     CGState(..),
     Testable(..)
 ) where
diff --git a/packages/base/src/Numeric/LinearAlgebra/Static.hs b/packages/base/src/Numeric/LinearAlgebra/Static.hs
index ded69fa..843c727 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Static.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Static.hs
@@ -65,7 +65,7 @@ import Numeric.LinearAlgebra hiding (
     row,col,vector,matrix,linspace,toRows,toColumns,
     (<\>),fromList,takeDiag,svd,eig,eigSH,
     eigenvalues,eigenvaluesSH,build,
-    qr,size,dot,chol,range,R,C,Her,her,sym)
+    qr,size,dot,chol,range,R,C,sym,mTm,unSym)
 import qualified Numeric.LinearAlgebra as LA
 import Data.Proxy(Proxy)
 import Internal.Static
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
index 30480d7..ecda6c1 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
@@ -500,10 +500,10 @@ sliceTest = utest "slice test" $ and
     , testSlice (test_qrgr 4 tau2) qrr2
     ]
   where
-    (qrr1,tau1) = qrRaw (rec :: Matrix R)
-    (qrr2,tau2) = qrRaw (rec :: Matrix C)
+    QR qrr1 tau1 = qrRaw (rec :: Matrix R)
+    QR qrr2 tau2 = qrRaw (rec :: Matrix C)
 
-    test_qrgr n t x = qrgr n (x,t)
+    test_qrgr n t x = qrgr n (QR x t)
 
     ok_qrgr x = simeq 1E-15 q q'
       where
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
index 4704989..3d5441d 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
@@ -81,8 +81,7 @@ instance (Field a, Arbitrary a) => Arbitrary (Rot a) where
 
 
 -- a complex hermitian or real symmetric matrix
---newtype (Her a) = Her (Matrix a) deriving Show
-instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Her a) where
+instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Herm a) where
     arbitrary = do
         Sq m <- arbitrary
         let m' = m/2
@@ -127,7 +126,7 @@ instance (Numeric a, ArbitraryField a, Num (Vector a))
     arbitrary = do
         m <- arbitrary
         let (_,v) = eigSH m
-            n = rows (her m)
+            n = rows (unSym m)
         l <- replicateM n (choose (0,100))
         let s = diag (fromList l)
             p = v <> real s <> tr v
@@ -161,8 +160,8 @@ fM m = m :: FM
 zM m = m :: ZM
 
 
-rHer m = her m :: RM
-cHer m = her m :: CM
+rHer m = unSym m :: RM
+cHer m = unSym m :: CM
 
 rRot (Rot m) = m :: RM
 cRot (Rot m) = m :: CM
@@ -176,8 +175,8 @@ cWC (WC m) = m :: CM
 rSqWC (SqWC m) = m :: RM
 cSqWC (SqWC m) = m :: CM
 
-rSymWC (SqWC m) = sym m :: Her R
-cSymWC (SqWC m) = sym m :: Her C
+rSymWC (SqWC m) = sym m :: Herm R
+cSymWC (SqWC m) = sym m :: Herm 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 2ac3588..720b7bd 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -263,7 +263,7 @@ 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))
+linearSolvePropH f m = f m (unSym m) |~| ident (rows (unSym m))
 
 linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b)
     where q = min (rows a) (cols a)