14 files changed, 271 insertions, 45 deletions

diff --git a/packages/base/CHANGELOG b/packages/base/CHANGELOG
index 35ccdbc..c137285 100644
--- a/packages/base/CHANGELOG
+++ b/packages/base/CHANGELOG
@@ -1,3 +1,10 @@
+0.16.1.0
+--------
+
+    * Added (|||) and (===) for "besides" and "above"
+
+    * rowOuters
+
 0.16.0.0
 --------
 
@@ -11,9 +18,11 @@
         Numeric.LinearAlgebra.Data
         Numeric.LinearAlgebra.Devel
 
-      The documentation is now hidden for Data.Packed.*, Numeric.Container, and
+      For normal usage we only need to import Numeric.LinearAlgebra.HMatrix.
+
+      (The documentation is now hidden for Data.Packed.*, Numeric.Container, and
       the other Numeric.LinearAlgebra.* modules, but they continue to be exposed
-      for backwards compatibility.
+      for backwards compatibility.)
 
     * Added support for empty arrays, extending automatic conformability
       (very useful for construction of block matrices).
diff --git a/packages/base/THANKS.md b/packages/base/THANKS.md
index b1417a6..fdbbe14 100644
--- a/packages/base/THANKS.md
+++ b/packages/base/THANKS.md
@@ -103,7 +103,7 @@ module reorganization, monadic mapVectorM, and many other improvements.
   deprecation warnings, used more explicit imports, and updated to ghc-7.4.
 
 - Tom Nielsen discovered a problem in Config.hs, exposed by link problems
-  in Ubuntu 11.10 beta.
+  in Ubuntu 11.10 beta, and fixed the link options on freebsd.
 
 - Daniel Fischer reported some Haddock markup errors.
 
@@ -120,7 +120,8 @@ module reorganization, monadic mapVectorM, and many other improvements.
 
 - Takano Akio fixed off-by-one errors in gsl-aux.c producing segfaults.
 
-- Alex Lang implemented uniRoot and uniRootJ for one-dimensional root-finding.
+- Alex Lang implemented uniRoot and uniRootJ for one-dimensional root-finding, and
+  fixed asRow and asColumn for empty vectors.
 
 - Mike Ledger contributed alternative FFI helpers for matrix interoperation with C
 
@@ -139,7 +140,8 @@ module reorganization, monadic mapVectorM, and many other improvements.
 
 - Clemens Lang updated the MacPort installation instructions.
 
-- Henning Thielemann reported the pinv inefficient implementation.
+- Henning Thielemann reported the pinv inefficient implementation and the need of
+  pkgconfig-depends.
 
 - bdoering reported the problem of zero absolute tolerance in the integration functions.
 
@@ -155,3 +157,25 @@ module reorganization, monadic mapVectorM, and many other improvements.
 
 - Samium Gromoff reported a build failure caused by a size_t - int mismatch.
 
+- Denis Laxalde separated the gsl tests from the base ones.
+
+- Dominic Steinitz (idontgetoutmuch) reported a bug in the static diagonal creation functions.
+
+- Dylan Thurston reported an error in the glpk documentation and ambiguity in
+  the description of linearSolve.
+
+- Adrian Victor Crisciu developed an installation method for platforms which
+  don't provide shared lapack libraries.
+
+- Ian Ross reported the max/minIndex bug.
+
+- Niklas Hambüchen improved the documentation.
+
+- "erdeszt" optimized "conv" using a direct vector reverse.
+
+- John Shahbazian added support for openBLAS.
+
+- "yongqli" reported the bug in randomVector (rand() is not thread-safe and drand48_r() is not portable).
+
+- Kiwamu Ishikura improved randomVector for OSX
+
diff --git a/packages/base/hmatrix.cabal b/packages/base/hmatrix.cabal
index ba719a1..8ba3e06 100644
--- a/packages/base/hmatrix.cabal
+++ b/packages/base/hmatrix.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix
-Version:            0.16.0.3
+Version:            0.16.1.3
 License:            BSD3
 License-file:       LICENSE
 Author:             Alberto Ruiz
@@ -14,6 +14,12 @@ Description:        Linear algebra based on BLAS and LAPACK.
                     ["Numeric.LinearAlgebra.HMatrix"] Starting point and recommended import module for most applications.
                     .
                     ["Numeric.LinearAlgebra.Static"] Experimental alternative interface.
+                    .
+                    ["Numeric.LinearAlgebra.Devel"] Tools for extending the library.
+                    .
+                    (Other modules are exposed with hidden documentation for backwards compatibility.)
+                    .
+                    Code examples: <http://dis.um.es/~alberto/hmatrix/hmatrix.html>
 
 Category:           Math
 tested-with:        GHC==7.8
@@ -26,6 +32,10 @@ extra-source-files: THANKS.md CHANGELOG
 
 extra-source-files: src/C/lapack-aux.h
 
+flag openblas
+    description:    Link with OpenBLAS (https://github.com/xianyi/OpenBLAS) optimized libraries.
+    default:        False
+
 library
 
     Build-Depends:      base >= 4 && < 5,
@@ -94,7 +104,11 @@ library
 
     cpp-options:        -DBINARY
 
-    extra-libraries:    blas lapack
+    if flag(openblas)
+        extra-lib-dirs:     /usr/lib/openblas/lib    
+        extra-libraries:    openblas
+    else 
+        extra-libraries:    blas lapack
 
     if os(OSX)
         extra-lib-dirs: /opt/local/lib/
@@ -108,7 +122,7 @@ library
     if os(freebsd)
        extra-lib-dirs: /usr/local/lib
        include-dirs: /usr/local/include
-       extra-libraries: blas lapack
+       extra-libraries: blas lapack gfortran
 
     if os(windows)
         extra-libraries: blas lapack
diff --git a/packages/base/src/C/vector-aux.c b/packages/base/src/C/vector-aux.c
index 2f47c8f..dda47cb 100644
--- a/packages/base/src/C/vector-aux.c
+++ b/packages/base/src/C/vector-aux.c
@@ -172,7 +172,7 @@ double vector_min(KDVEC(x)) {
 double vector_max_index(KDVEC(x)) {
     int k, r = 0;
     for (k = 1; k<xn; k++) {
-        if(xp[k]>xp[0]) {
+        if(xp[k]>xp[r]) {
             r = k;
         }
     }
@@ -182,7 +182,7 @@ double vector_max_index(KDVEC(x)) {
 double vector_min_index(KDVEC(x)) {
     int k, r = 0;
     for (k = 1; k<xn; k++) {
-        if(xp[k]<xp[0]) {
+        if(xp[k]<xp[r]) {
             r = k;
         }
     }
@@ -237,7 +237,7 @@ float vector_min_f(KFVEC(x)) {
 float vector_max_index_f(KFVEC(x)) {
     int k, r = 0;
     for (k = 1; k<xn; k++) {
-        if(xp[k]>xp[0]) {
+        if(xp[k]>xp[r]) {
             r = k;
         }
     }
@@ -247,7 +247,7 @@ float vector_max_index_f(KFVEC(x)) {
 float vector_min_index_f(KFVEC(x)) {
     int k, r = 0;
     for (k = 1; k<xn; k++) {
-        if(xp[k]<xp[0]) {
+        if(xp[k]<xp[r]) {
             r = k;
         }
     }
@@ -700,17 +700,86 @@ int saveMatrix(char * file, char * format, KDMAT(a)){
 
 ////////////////////////////////////////////////////////////////////////////////
 
+#ifdef __APPLE__
+
+#pragma message "randomVector is not thread-safe in OSX"
+
+inline double urandom() {
+    const long max_random = 2147483647; // 2**31 - 1
+    return (double)random() / (double)max_random;
+}
+
+double gaussrand(int *phase, double *pV1, double *pV2, double *pS)
+{
+        double V1=*pV1, V2=*pV2, S=*pS;
+        double X;
+
+        if(*phase == 0) {
+                do {
+            double U1 = urandom();
+                        double U2 = urandom();
+
+                        V1 = 2 * U1 - 1;
+                        V2 = 2 * U2 - 1;
+                        S = V1 * V1 + V2 * V2;
+                        } while(S >= 1 || S == 0);
+
+                X = V1 * sqrt(-2 * log(S) / S);
+        } else
+                X = V2 * sqrt(-2 * log(S) / S);
+
+        *phase = 1 - *phase;
+    *pV1=V1; *pV2=V2; *pS=S;
+
+        return X;
+
+}
+
+int random_vector(unsigned int seed, int code, DVEC(r)) {
+    int phase = 0;
+    double V1,V2,S;
+
+    srandom(seed);
+    
+    int k;
+    switch (code) {
+      case 0: { // uniform
+        for (k=0; k<rn; k++) {
+            rp[k] = urandom();
+        }
+        OK
+      }
+      case 1: { // gaussian
+        for (k=0; k<rn; k++) {
+            rp[k] = gaussrand(&phase,&V1,&V2,&S);
+        }
+        OK
+      }
+
+      default: ERROR(BAD_CODE);
+    }
+}
+
+#else
+
+inline double urandom(struct random_data * buffer) {
+    int32_t res;
+    random_r(buffer,&res);
+    return (double)res/RAND_MAX;
+}
+
+
 // http://c-faq.com/lib/gaussian.html
-double gaussrand()
+double gaussrand(struct random_data *buffer,
+                 int *phase, double *pV1, double *pV2, double *pS)
 {
-        static double V1, V2, S;
-        static int phase = 0;
+        double V1=*pV1, V2=*pV2, S=*pS;
         double X;
 
-        if(phase == 0) {
+        if(*phase == 0) {
                 do {
-                        double U1 = (double)rand() / RAND_MAX;
-                        double U2 = (double)rand() / RAND_MAX;
+            double U1 = urandom(buffer);
+                        double U2 = urandom(buffer);
 
                         V1 = 2 * U1 - 1;
                         V2 = 2 * U2 - 1;
@@ -721,24 +790,37 @@ double gaussrand()
         } else
                 X = V2 * sqrt(-2 * log(S) / S);
 
-        phase = 1 - phase;
+        *phase = 1 - *phase;
+    *pV1=V1; *pV2=V2; *pS=S;
 
         return X;
+
 }
 
-int random_vector(int seed, int code, DVEC(r)) {
-    srand(seed);
+int random_vector(unsigned int seed, int code, DVEC(r)) {
+    struct random_data buffer;
+    char   random_state[128];
+    memset(&buffer, 0, sizeof(struct random_data));
+    memset(random_state, 0, sizeof(random_state));
+    
+    initstate_r(seed,random_state,sizeof(random_state),&buffer);
+    // setstate_r(random_state,&buffer);
+    // srandom_r(seed,&buffer);
+
+    int phase = 0;
+    double V1,V2,S;
+
     int k;
     switch (code) {
       case 0: { // uniform
         for (k=0; k<rn; k++) {
-            rp[k] = (double)rand()/RAND_MAX;
+            rp[k] = urandom(&buffer);
         }
         OK
       }
       case 1: { // gaussian
         for (k=0; k<rn; k++) {
-            rp[k] = gaussrand();
+            rp[k] = gaussrand(&buffer,&phase,&V1,&V2,&S);
         }
         OK
       }
@@ -747,6 +829,8 @@ int random_vector(int seed, int code, DVEC(r)) {
     }
 }
 
+#endif
+
 ////////////////////////////////////////////////////////////////////////////////
 
 int smXv(KDVEC(vals),KIVEC(cols),KIVEC(rows),KDVEC(x),DVEC(r)) {
diff --git a/packages/base/src/Data/Packed/Internal/Matrix.hs b/packages/base/src/Data/Packed/Internal/Matrix.hs
index 91a9466..150b978 100644
--- a/packages/base/src/Data/Packed/Internal/Matrix.hs
+++ b/packages/base/src/Data/Packed/Internal/Matrix.hs
@@ -253,7 +253,8 @@ compat m1 m2 = rows m1 == rows m2 && cols m1 == cols m2
     operations for selected element types.
     It provides unoptimised defaults for any 'Storable' type,
     so you can create instances simply as:
-    @instance Element Foo@.
+
+    >instance Element Foo
 -}
 class (Storable a) => Element a where
     subMatrixD :: (Int,Int) -- ^ (r0,c0) starting position 
diff --git a/packages/base/src/Data/Packed/Matrix.hs b/packages/base/src/Data/Packed/Matrix.hs
index 6445ce4..70b9232 100644
--- a/packages/base/src/Data/Packed/Matrix.hs
+++ b/packages/base/src/Data/Packed/Matrix.hs
@@ -284,7 +284,7 @@ fromLists = fromRows . map fromList
 --   [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
 --
 asRow :: Storable a => Vector a -> Matrix a
-asRow v = reshape (dim v) v
+asRow = trans . asColumn
 
 -- | creates a 1-column matrix from a vector
 --
@@ -297,7 +297,7 @@ asRow v = reshape (dim v) v
 --  , 5.0 ]
 --
 asColumn :: Storable a => Vector a -> Matrix a
-asColumn = trans . asRow
+asColumn v = reshape 1 v
 
 
 
@@ -476,7 +476,7 @@ mapMatrixWithIndexM g m = liftM (reshape c) . mapVectorWithIndexM (mk c g) . fla
 
 {- |
 
->>> mapMatrixWithIndex (\\(i,j) v -> 100*v + 10*fromIntegral i + fromIntegral j) (ident 3:: Matrix Double)
+>>> mapMatrixWithIndex (\(i,j) v -> 100*v + 10*fromIntegral i + fromIntegral j) (ident 3:: Matrix Double)
 (3><3)
  [ 100.0,   1.0,   2.0
  ,  10.0, 111.0,  12.0
diff --git a/packages/base/src/Data/Packed/Numeric.hs b/packages/base/src/Data/Packed/Numeric.hs
index 7aa53f1..6027f43 100644
--- a/packages/base/src/Data/Packed/Numeric.hs
+++ b/packages/base/src/Data/Packed/Numeric.hs
@@ -160,7 +160,29 @@ instance Mul Vector Matrix Vector where
 
 --------------------------------------------------------------------------------
 
--- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD)
+{- | Least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD)
+
+@
+a = (3><2)
+ [ 1.0,  2.0
+ , 2.0,  4.0
+ , 2.0, -1.0 ]
+@
+
+@
+v = vector [13.0,27.0,1.0]
+@
+
+>>> let x = a <\> v
+>>> x
+fromList [3.0799999999999996,5.159999999999999]
+
+>>> a #> x
+fromList [13.399999999999999,26.799999999999997,1.0]
+
+It also admits multiple right-hand sides stored as columns in a matrix.
+
+-}
 infixl 7 <\>
 (<\>) :: (LSDiv c, Field t) => Matrix t -> c t -> c t
 (<\>) = linSolve
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)