11 files changed, 216 insertions, 37 deletions

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)