7 files changed, 89 insertions, 17 deletions

diff --git a/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs b/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs
index 25700bc..02ac6a0 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs
@@ -809,7 +809,7 @@ expGolub m = iterate msq f !! j
 --------------------------------------------------------------
 
 {- | Matrix square root. Currently it uses a simple iterative algorithm described in Wikipedia.
-It only works with invertible matrices that have a real solution. For diagonalizable matrices you can try @matFunc sqrt@.
+It only works with invertible matrices that have a real solution.
 
 @m = (2><2) [4,9
            ,0,4] :: Matrix Double@
@@ -819,6 +819,13 @@ It only works with invertible matrices that have a real solution. For diagonaliz
  [ 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
diff --git a/packages/base/src/Numeric/LinearAlgebra/Data.hs b/packages/base/src/Numeric/LinearAlgebra/Data.hs
index b1a31fc..6dea407 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Data.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Data.hs
@@ -40,7 +40,7 @@ module Numeric.LinearAlgebra.Data(
     takeRows, dropRows, takeColumns, dropColumns, subMatrix, (?), (¿), fliprl, flipud,
   
     -- * Block matrix
-    fromBlocks, (¦), (——), diagBlock, repmat, toBlocks, toBlocksEvery,
+    fromBlocks, (|||), (===), diagBlock, repmat, toBlocks, toBlocksEvery,
 
     -- * Mapping functions
     conj, cmap, step, cond,
@@ -66,7 +66,7 @@ module Numeric.LinearAlgebra.Data(
     arctan2,
     rows, cols,
     separable,
-
+    (¦),(——),
     module Data.Complex,
 
     Vector, Matrix, GMatrix, nRows, nCols
diff --git a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs
index d2cae6c..677f9ee 100644
--- a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs
@@ -134,7 +134,7 @@ module Numeric.LinearAlgebra.HMatrix (
     Seed, RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample,
 
     -- * Misc
-    meanCov, peps, relativeError, haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv,
+    meanCov, rowOuters, peps, relativeError, haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv,
     ℝ,ℂ,iC,
     -- * Auxiliary classes
     Element, Container, Product, Numeric, LSDiv,
@@ -194,7 +194,37 @@ mul :: Numeric t => Matrix t -> Matrix t -> Matrix t
 mul = mXm
 
 
--- | 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'.
+
+@
+a = (2><2)
+ [ 1.0, 2.0
+ , 3.0, 5.0 ]
+@
+
+@
+b = (2><3)
+ [  6.0, 1.0, 10.0
+ , 15.0, 3.0, 26.0 ]
+@
+
+>>> linearSolve a b
+Just (2><3)
+ [ -1.4802973661668753e-15,     0.9999999999999997, 1.999999999999997
+ ,       3.000000000000001, 1.6653345369377348e-16, 4.000000000000002 ]
+
+>>> let Just x = it
+>>> disp 5 x
+2x3
+-0.00000  1.00000  2.00000
+ 3.00000  0.00000  4.00000
+
+>>> a <> x
+(2><3)
+ [  6.0, 1.0, 10.0
+ , 15.0, 3.0, 26.0 ]
+
+-}
 linearSolve m b = A.mbLinearSolve m b
 
 -- | return an orthonormal basis of the null space of a matrix. See also 'nullspaceSVD'.
diff --git a/packages/base/src/Numeric/LinearAlgebra/Static.hs b/packages/base/src/Numeric/LinearAlgebra/Static.hs
index 5749c40..037396d 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Static.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Static.hs
@@ -512,7 +512,10 @@ crossC (extract -> x) (extract -> y) = mkC (LA.fromList [z1, z2, z3])
 --------------------------------------------------------------------------------
 
 diagRectR :: forall m n k . (KnownNat m, KnownNat n, KnownNat k) => ℝ -> R k -> L m n
-diagRectR x v = r
+diagRectR x v
+    | m' == 1 = mkL (LA.diagRect x ev m' n')
+    | m'*n' > 0 = r
+    | otherwise = matrix []
   where
     r = mkL (asRow (vjoin [scalar x, ev, zeros]))
     ev = extract v
@@ -521,7 +524,10 @@ diagRectR x v = r
 
 
 diagRectC :: forall m n k . (KnownNat m, KnownNat n, KnownNat k) => ℂ -> C k -> M m n
-diagRectC x v = r
+diagRectC x v
+    | m' == 1 = mkM (LA.diagRect x ev m' n')
+    | m'*n' > 0 = r
+    | otherwise = fromList []
   where
     r = mkM (asRow (vjoin [scalar x, ev, zeros]))
     ev = extract v
diff --git a/packages/base/src/Numeric/LinearAlgebra/Static/Internal.hs b/packages/base/src/Numeric/LinearAlgebra/Static/Internal.hs
index 339ef7d..ec02cf6 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Static/Internal.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Static/Internal.hs
@@ -150,7 +150,7 @@ gmat st xs'
     (xs,rest) = splitAt (m'*n') xs'
     v = LA.fromList xs
     x = reshape n' v
-    ok = rem (LA.size v) n' == 0 && LA.size x == (m',n') && null rest
+    ok = null rest && ((n' == 0 && dim v == 0) || n'> 0 && (rem (LA.size v) n' == 0) && LA.size x == (m',n'))
     m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
     n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
     abort info = error $ st ++" "++show m' ++ " " ++ show n'++" can't be created from elements " ++ info
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util.hs b/packages/base/src/Numeric/LinearAlgebra/Util.hs
index 6bb9d15..043aa21 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Util.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Util.hs
@@ -33,7 +33,7 @@ module Numeric.LinearAlgebra.Util(
     diagl,
     row,
     col,
-    (&), (¦), (——), (#),
+    (&), (¦), (|||), (——), (===), (#),
     (?), (¿),
     Indexable(..), size,
     Numeric,
@@ -157,25 +157,34 @@ a & b = vjoin [a,b]
 
 {- | horizontal concatenation of real matrices
 
- (unicode 0x00a6, broken bar)
-
->>> ident 3 ¦ konst 7 (3,4)
+>>> ident 3 ||| konst 7 (3,4)
 (3><7)
  [ 1.0, 0.0, 0.0, 7.0, 7.0, 7.0, 7.0
  , 0.0, 1.0, 0.0, 7.0, 7.0, 7.0, 7.0
  , 0.0, 0.0, 1.0, 7.0, 7.0, 7.0, 7.0 ]
 
 -}
+infixl 3 |||
+(|||) :: Matrix Double -> Matrix Double -> Matrix Double
+a ||| b = fromBlocks [[a,b]]
+
+-- | a synonym for ('|||') (unicode 0x00a6, broken bar)
 infixl 3 ¦
 (¦) :: Matrix Double -> Matrix Double -> Matrix Double
-a ¦ b = fromBlocks [[a,b]]
+(¦) = (|||)
+
 
 -- | vertical concatenation of real matrices
 --
--- (unicode 0x2014, em dash)
+(===) :: Matrix Double -> Matrix Double -> Matrix Double
+infixl 2 ===
+a === b = fromBlocks [[a],[b]]
+
+-- | a synonym for ('===') (unicode 0x2014, em dash)
 (——) :: Matrix Double -> Matrix Double -> Matrix Double
 infixl 2 ——
-a —— b = fromBlocks [[a],[b]]
+(——) = (===)
+
 
 (#) :: Matrix Double -> Matrix Double -> Matrix Double
 infixl 2 #
@@ -356,7 +365,26 @@ pairwiseD2 x y | ok = x2 `outer` oy + ox `outer` y2 - 2* x <> trans y
 
 --------------------------------------------------------------------------------
 
--- | outer products of rows
+{- | outer products of rows
+
+>>> a
+(3><2)
+ [   1.0,   2.0
+ ,  10.0,  20.0
+ , 100.0, 200.0 ]
+>>> b
+(3><3)
+ [ 1.0, 2.0, 3.0
+ , 4.0, 5.0, 6.0
+ , 7.0, 8.0, 9.0 ]
+
+>>> rowOuters a (b ||| 1)
+(3><8)
+ [   1.0,   2.0,   3.0,   1.0,    2.0,    4.0,    6.0,   2.0
+ ,  40.0,  50.0,  60.0,  10.0,   80.0,  100.0,  120.0,  20.0
+ , 700.0, 800.0, 900.0, 100.0, 1400.0, 1600.0, 1800.0, 200.0 ]
+
+-}
 rowOuters :: Matrix Double -> Matrix Double -> Matrix Double
 rowOuters a b = a' * b'
   where
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs b/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs
index c8c7536..c9e75de 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs
@@ -16,6 +16,7 @@ module Numeric.LinearAlgebra.Util.Convolution(
    corr2, conv2, separable
 ) where
 
+import qualified Data.Vector.Storable as SV
 import Data.Packed.Numeric
 
 
@@ -51,7 +52,7 @@ conv ker v
     | dim ker == 0 = konst 0 (dim v)
     | otherwise = corr ker' v'
   where
-    ker' = (flatten.fliprl.asRow) ker
+    ker' = SV.reverse ker
     v' = vjoin [z,v,z]
     z = konst 0 (dim ker -1)