3 files changed, 89 insertions, 34 deletions

diff --git a/packages/base/src/Numeric/HMatrix.hs b/packages/base/src/Numeric/HMatrix.hs
index 9d34658..ec96bfc 100644
--- a/packages/base/src/Numeric/HMatrix.hs
+++ b/packages/base/src/Numeric/HMatrix.hs
@@ -17,10 +17,10 @@ module Numeric.HMatrix (
     -- |
     -- The standard numeric classes are defined elementwise:
     --
-    -- >>> fromList [1,2,3] * fromList [3,0,-2 :: Double]
+    -- >>>  vect [1,2,3] * vect [3,0,-2]
     -- fromList [3.0,0.0,-6.0]
     --
-    -- >>> (3><3) [1..9] * ident 3 :: Matrix Double
+    -- >>> mat 3 [1..9] * ident 3
     -- (3><3)
     --  [ 1.0, 0.0, 0.0
     --  , 0.0, 5.0, 0.0
@@ -36,6 +36,12 @@ module Numeric.HMatrix (
     --  , 5.0, 7.0, 5.0
     --  , 5.0, 5.0, 7.0 ]
     --
+    -- >>> mat 3 [1..9] + mat 1 [10,20,30]
+    -- (3><3)
+    --  [ 11.0, 12.0, 13.0
+    --  , 24.0, 25.0, 26.0
+    --  , 37.0, 38.0, 39.0 ]
+    --
 
     -- * Products
     -- ** dot
@@ -48,11 +54,12 @@ module Numeric.HMatrix (
     -- single-element matrices (created from numeric literals or using 'scalar')
     -- are used for scaling.
     --
-    -- >>> let m = (2><3)[1..] :: Matrix Double
-    -- >>> m <> 2 <> diagl[0.5,1,0]
+    -- >>> import Data.Monoid as M
+    -- >>>  let m = mat 3 [1..6]
+    -- >>> m M.<> 2 M.<> diagl[0.5,1,0]
     -- (2><3)
-    -- [ 1.0,  4.0, 0.0
-    -- , 4.0, 10.0, 0.0 ]
+    --  [ 1.0,  4.0, 0.0
+    --  , 4.0, 10.0, 0.0 ]
     --
     -- 'mconcat' uses 'optimiseMult' to get the optimal association order.
 
@@ -76,10 +83,18 @@ module Numeric.HMatrix (
     inv, pinv, pinvTol,
 
     -- * Determinant and rank
-    rcond, rank, ranksv,
+    rcond, rank,
     det, invlndet,
 
-    -- * Singular value decomposition
+    -- * Norms
+    Normed(..),
+    norm_Frob, norm_nuclear,
+
+    -- * Nullspace and range
+    orth,
+    nullspace, null1, null1sym,
+
+    -- * SVD
     svd,
     fullSVD,
     thinSVD,
@@ -112,18 +127,6 @@ module Numeric.HMatrix (
     sqrtm,
     matFunc,
 
-    -- * Nullspace
-    nullspacePrec,
-    nullVector,
-    nullspaceSVD,
-    null1, null1sym,
-
-    orth,
-
-    -- * Norms
-    Normed(..),
-    norm_Frob, norm_nuclear,
-
     -- * Correlation and convolution
     corr, conv, corrMin, corr2, conv2,
 
@@ -132,7 +135,8 @@ module Numeric.HMatrix (
     Seed, RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample,
 
     -- * Misc
-    meanCov, peps, relativeError, haussholder, optimiseMult, udot,
+    meanCov, peps, relativeError, haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv,
+    ℝ,ℂ,iC,
     -- * Auxiliary classes
     Element, Container, Product, Numeric, LSDiv,
     Complexable, RealElement,
@@ -142,8 +146,7 @@ module Numeric.HMatrix (
 --    Normed,
     Transposable,
     CGState(..),
-    Testable(..),
-    ℕ,ℤ,ℝ,ℂ, i_C
+    Testable(..)
 ) where
 
 import Numeric.LinearAlgebra.Data
@@ -151,14 +154,38 @@ import Numeric.LinearAlgebra.Data
 import Numeric.Matrix()
 import Numeric.Vector()
 import Data.Packed.Numeric hiding ((<>))
-import Numeric.LinearAlgebra.Algorithms hiding (linearSolve,Normed)
+import Numeric.LinearAlgebra.Algorithms hiding (linearSolve,Normed,orth)
 import qualified Numeric.LinearAlgebra.Algorithms as A
 import Numeric.LinearAlgebra.Util
 import Numeric.LinearAlgebra.Random
 import Numeric.Sparse((!#>))
 import Numeric.LinearAlgebra.Util.CG
 
--- | matrix product
+{- | dense matrix product
+
+>>> let a = (3><5) [1..]
+>>> a
+(3><5)
+ [  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 ]
+
+>>> let b = (5><2) [1,3, 0,2, -1,5, 7,7, 6,0]
+>>> b
+(5><2)
+ [  1.0, 3.0
+ ,  0.0, 2.0
+ , -1.0, 5.0
+ ,  7.0, 7.0
+ ,  6.0, 0.0 ]
+
+>>> a <> b
+(3><2)
+ [  56.0,  50.0
+ , 121.0, 135.0
+ , 186.0, 220.0 ]
+
+-}
 (<>) :: Numeric t => Matrix t -> Matrix t -> Matrix t
 (<>) = mXm
 infixr 8 <>
@@ -166,3 +193,9 @@ infixr 8 <>
 -- | 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'. 
 linearSolve m b = A.mbLinearSolve m b
 
+-- | return an orthonormal basis of the null space of a matrix. See also 'nullspaceSVD'.
+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/Util.hs b/packages/base/src/Numeric/LinearAlgebra/Util.hs
index 324fb44..4824af4 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Util.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Util.hs
@@ -32,7 +32,7 @@ module Numeric.LinearAlgebra.Util(
     rand, randn,
     cross,
     norm,
-    ℕ,ℤ,ℝ,ℂ,𝑖,i_C, --ℍ
+    ℕ,ℤ,ℝ,ℂ,iC,
     Normed(..), norm_Frob, norm_nuclear,
     unitary,
     mt,
@@ -72,13 +72,10 @@ type ℝ = Double
 type ℕ = Int
 type ℤ = Int
 type ℂ = Complex Double
---type ℝn = Vector ℝ
---type ℂn = Vector ℂ
---newtype ℍ m = H m
 
-i_C, 𝑖 :: ℂ
-𝑖 = 0:+1
-i_C = 𝑖
+-- | imaginary unit
+iC :: ℂ
+iC = 0:+1
 
 {- | create a real vector
 
diff --git a/packages/base/src/Numeric/Sparse.hs b/packages/base/src/Numeric/Sparse.hs
index 1b8a7b3..f495e3a 100644
--- a/packages/base/src/Numeric/Sparse.hs
+++ b/packages/base/src/Numeric/Sparse.hs
@@ -62,7 +62,26 @@ mkCSR sm' = CSR{..}
     csrNRows = dim csrRows - 1
     csrNCols = fromIntegral (V.maximum csrCols)
 
-
+{- | General matrix with specialized internal representations for
+     dense, sparse, diagonal, banded, and constant elements.
+
+>>> let m = mkSparse [((0,999),1.0),((1,1999),2.0)]
+>>> m
+SparseR {gmCSR = CSR {csrVals = fromList [1.0,2.0],
+                      csrCols = fromList [1000,2000],
+                      csrRows = fromList [1,2,3],
+                      csrNRows = 2,
+                      csrNCols = 2000},
+                      nRows = 2,
+                      nCols = 2000}
+
+>>> let m = mkDense (mat 2 [1..4])
+>>> m
+Dense {gmDense = (2><2)
+ [ 1.0, 2.0
+ , 3.0, 4.0 ], nRows = 2, nCols = 2}
+
+-}
 data GMatrix
     = SparseR
         { gmCSR   :: CSR
@@ -146,7 +165,13 @@ gmXv Dense{..} v
                                  nRows nCols (dim v)
 
 
--- | general matrix - vector product
+{- | general matrix - vector product
+
+>>> let m = mkSparse [((0,999),1.0),((1,1999),2.0)]
+>>> m !#> vect[1..2000]
+fromList [1000.0,4000.0]
+
+-}
 infixr 8 !#>
 (!#>) :: GMatrix -> Vector Double -> Vector Double
 (!#>) = gmXv