24 files changed, 1107 insertions, 0 deletions

diff --git a/packages/hmatrix/examples/bool.hs b/packages/hmatrix/examples/bool.hs
new file mode 100644
index 0000000..679b8bf
--- /dev/null
+++ b/packages/hmatrix/examples/bool.hs
@@ -0,0 +1,54 @@
+-- vectorized boolean operations defined in terms of step or cond
+
+import Numeric.LinearAlgebra
+
+infix  4  .==., ./=., .<., .<=., .>=., .>.
+infixr 3  .&&.
+infixr 2  .||.
+
+a .<.  b = step (b-a)
+a .<=. b = cond a b 1 1 0
+a .==. b = cond a b 0 1 0
+a ./=. b = cond a b 1 0 1
+a .>=. b = cond a b 0 1 1
+a .>.  b = step (a-b)
+
+a .&&. b  = step (a*b)
+a .||. b  = step (a+b)
+no a      = 1-a
+xor a b   = a ./=. b
+equiv a b = a .==. b
+imp a b   = no a .||. b
+
+taut x = minElement x == 1
+
+minEvery a b = cond a b a a b
+maxEvery a b = cond a b b b a
+
+-- examples
+
+clip a b x = cond y b y y b where y = cond x a a x x
+
+disp = putStr . dispf 3
+
+eye n = ident n :: Matrix Double
+row = asRow . fromList    :: [Double] -> Matrix Double
+col = asColumn . fromList :: [Double] -> Matrix Double
+
+m = (3><4) [1..] :: Matrix Double
+
+p = row [0,0,1,1]
+q = row [0,1,0,1]
+
+main = do
+    print $ find (>6) m
+    disp $ assoc (6,8) 7 $ zip (find (/=0) (eye 5)) [10..]
+    disp $ accum (eye 5) (+) [((0,2),3), ((3,1),7), ((1,1),1)]
+    disp $ m .>=. 10  .||.  m .<. 4
+    (disp . fromColumns . map flatten) [p, q, p.&&.q, p .||.q, p `xor` q, p `equiv` q, p `imp` q]
+    print $ taut $ (p `imp` q ) `equiv` (no q `imp` no p)
+    print $ taut $ (xor p q) `equiv` (p .&&. no q .||. no p .&&. q)
+    disp $ clip 3 8 m
+    disp $ col [1..7] .<=. row [1..5]
+    disp $ cond (col [1..3]) (row [1..4]) m 50 (3*m)
+
diff --git a/packages/hmatrix/examples/data.txt b/packages/hmatrix/examples/data.txt
new file mode 100644
index 0000000..2b9bfea
--- /dev/null
+++ b/packages/hmatrix/examples/data.txt
@@ -0,0 +1,10 @@
+ 0.9    1.1
+ 2.1    3.9
+ 3.1    9.2
+ 4.0   51.8
+ 4.9   25.3
+ 6.1   35.7
+ 7.0   49.4
+ 7.9    3.6
+ 9.1   81.5
+10.2   99.5
\ No newline at end of file
diff --git a/packages/hmatrix/examples/deriv.hs b/packages/hmatrix/examples/deriv.hs
new file mode 100644
index 0000000..c9456d1
--- /dev/null
+++ b/packages/hmatrix/examples/deriv.hs
@@ -0,0 +1,8 @@
+-- Numerical differentiation
+
+import Numeric.GSL
+
+d :: (Double -> Double) -> (Double -> Double)
+d f x = fst $ derivCentral 0.01 f x
+
+main = print $ d (\x-> x * d (\y-> x+y) 1) 1
diff --git a/packages/hmatrix/examples/devel/ej1/functions.c b/packages/hmatrix/examples/devel/ej1/functions.c
new file mode 100644
index 0000000..02a4cdd
--- /dev/null
+++ b/packages/hmatrix/examples/devel/ej1/functions.c
@@ -0,0 +1,35 @@
+/* assuming row order */
+
+typedef struct { double r, i; } doublecomplex;
+
+#define DVEC(A) int A##n, double*A##p
+#define CVEC(A) int A##n, doublecomplex*A##p
+#define DMAT(A) int A##r, int A##c, double*A##p
+#define CMAT(A) int A##r, int A##c, doublecomplex*A##p
+
+#define AT(M,row,col) (M##p[(row)*M##c + (col)])
+
+/*-----------------------------------------------------*/
+
+int c_scale_vector(double s, DVEC(x), DVEC(y)) {
+    int k;
+    for (k=0; k<=yn; k++) {
+        yp[k] = s*xp[k];
+    }
+    return 0;
+}
+
+/*-----------------------------------------------------*/
+
+int c_diag(DMAT(m),DVEC(y),DMAT(z)) {
+    int i,j;
+    for (j=0; j<yn; j++) {
+        yp[j] = AT(m,j,j);
+    }
+    for (i=0; i<mr; i++) {
+        for (j=0; j<mc; j++) {
+            AT(z,i,j) = i==j?yp[i]:0;
+        }
+    }
+    return 0;
+}
diff --git a/packages/hmatrix/examples/devel/ej1/wrappers.hs b/packages/hmatrix/examples/devel/ej1/wrappers.hs
new file mode 100644
index 0000000..a88f74b
--- /dev/null
+++ b/packages/hmatrix/examples/devel/ej1/wrappers.hs
@@ -0,0 +1,44 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+
+-- $ ghc -O2 --make wrappers.hs functions.c
+
+import Numeric.LinearAlgebra
+import Data.Packed.Development
+import Foreign(Ptr,unsafePerformIO)
+import Foreign.C.Types(CInt)
+
+-----------------------------------------------------
+
+main = do
+    print $ myScale 3.0 (fromList [1..10])
+    print $ myDiag $ (3><5) [1..]
+
+-----------------------------------------------------
+
+foreign import ccall unsafe "c_scale_vector"
+    cScaleVector :: Double                -- scale
+                 -> CInt -> Ptr Double    -- argument
+                 -> CInt -> Ptr Double    -- result
+                 -> IO CInt               -- exit code
+
+myScale s x = unsafePerformIO $ do
+    y <- createVector (dim x)
+    app2 (cScaleVector s) vec x vec y "cScaleVector"
+    return y
+
+-----------------------------------------------------
+-- forcing row order
+
+foreign import ccall unsafe "c_diag"
+    cDiag :: CInt -> CInt -> Ptr Double  -- argument
+          -> CInt -> Ptr Double          -- result1
+          -> CInt -> CInt -> Ptr Double  -- result2
+          -> IO CInt                     -- exit code
+
+myDiag m = unsafePerformIO $ do
+    y <- createVector (min r c)
+    z <- createMatrix RowMajor r c
+    app3 cDiag mat (cmat m) vec y mat z "cDiag"
+    return (y,z)
+  where r = rows m
+        c = cols m
diff --git a/packages/hmatrix/examples/devel/ej2/functions.c b/packages/hmatrix/examples/devel/ej2/functions.c
new file mode 100644
index 0000000..4dcd377
--- /dev/null
+++ b/packages/hmatrix/examples/devel/ej2/functions.c
@@ -0,0 +1,24 @@
+/*  general element order   */
+
+typedef struct { double r, i; } doublecomplex;
+
+#define DVEC(A) int A##n, double*A##p
+#define CVEC(A) int A##n, doublecomplex*A##p
+#define DMAT(A) int A##r, int A##c, double*A##p
+#define CMAT(A) int A##r, int A##c, doublecomplex*A##p
+
+#define AT(M,r,c) (M##p[(r)*sr+(c)*sc])
+
+int c_diag(int ro, DMAT(m),DVEC(y),DMAT(z)) {
+    int i,j,sr,sc;
+    if (ro==1) { sr = mc; sc = 1;} else { sr = 1; sc = mr;}
+    for (j=0; j<yn; j++) {
+        yp[j] = AT(m,j,j);
+    }
+    for (i=0; i<mr; i++) {
+        for (j=0; j<mc; j++) {
+            AT(z,i,j) = i==j?yp[i]:0;
+        }
+    }
+    return 0;
+}
diff --git a/packages/hmatrix/examples/devel/ej2/wrappers.hs b/packages/hmatrix/examples/devel/ej2/wrappers.hs
new file mode 100644
index 0000000..1c02a24
--- /dev/null
+++ b/packages/hmatrix/examples/devel/ej2/wrappers.hs
@@ -0,0 +1,32 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+
+-- $ ghc -O2 --make wrappers.hs functions.c
+
+import Numeric.LinearAlgebra
+import Data.Packed.Development
+import Foreign(Ptr,unsafePerformIO)
+import Foreign.C.Types(CInt)
+
+-----------------------------------------------------
+
+main = do
+    print $ myDiag $ (3><5) [1..]
+
+-----------------------------------------------------
+-- arbitrary data order
+
+foreign import ccall unsafe "c_diag"
+    cDiag :: CInt                        -- matrix order
+          -> CInt -> CInt -> Ptr Double  -- argument
+          -> CInt -> Ptr Double          -- result1
+          -> CInt -> CInt -> Ptr Double  -- result2
+          -> IO CInt                     -- exit code
+
+myDiag m = unsafePerformIO $ do
+    y <- createVector (min r c)
+    z <- createMatrix (orderOf m) r c
+    app3 (cDiag o) mat m vec y mat z "cDiag"
+    return (y,z)
+  where r = rows m
+        c = cols m
+        o = if orderOf m == RowMajor then 1 else 0
diff --git a/packages/hmatrix/examples/error.hs b/packages/hmatrix/examples/error.hs
new file mode 100644
index 0000000..5efae7c
--- /dev/null
+++ b/packages/hmatrix/examples/error.hs
@@ -0,0 +1,21 @@
+import Numeric.GSL
+import Numeric.GSL.Special
+import Numeric.LinearAlgebra
+import Prelude hiding (catch)
+import Control.Exception
+
+test x = catch
+       (print x)
+       (\e -> putStrLn $ "captured ["++ show (e :: SomeException) ++"]")
+
+main = do
+    setErrorHandlerOff
+
+    test $ log_e (-1)
+    test $ 5 + (fst.exp_e) 1000
+    test $ bessel_zero_Jnu_e (-0.3) 2
+
+    test $ (linearSolve 0 4 :: Matrix Double)
+    test $ (linearSolve 5 (sqrt (-1)) :: Matrix Double)
+
+    putStrLn "Bye"
\ No newline at end of file
diff --git a/packages/hmatrix/examples/fitting.hs b/packages/hmatrix/examples/fitting.hs
new file mode 100644
index 0000000..a8f6b1c
--- /dev/null
+++ b/packages/hmatrix/examples/fitting.hs
@@ -0,0 +1,24 @@
+-- nonlinear least-squares fitting
+
+import Numeric.GSL.Fitting
+import Numeric.LinearAlgebra
+
+xs = map return [0 .. 39]
+sigma = 0.1
+ys = map return $ toList $ fromList (map (head . expModel [5,0.1,1]) xs)
+            + scalar sigma * (randomVector 0 Gaussian 40)
+
+dat :: [([Double],([Double],Double))]
+
+dat = zip xs (zip ys (repeat sigma))
+
+expModel [a,lambda,b] [t] = [a * exp (-lambda * t) + b]
+
+expModelDer [a,lambda,b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]]
+
+(sol,path) = fitModelScaled 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0]
+
+main = do
+    print dat
+    print path
+    print sol
diff --git a/packages/hmatrix/examples/inplace.hs b/packages/hmatrix/examples/inplace.hs
new file mode 100644
index 0000000..dcfff56
--- /dev/null
+++ b/packages/hmatrix/examples/inplace.hs
@@ -0,0 +1,152 @@
+-- some tests of the interface for pure
+-- computations with inplace updates
+
+import Numeric.LinearAlgebra
+import Data.Packed.ST
+import Data.Packed.Convert
+
+import Data.Array.Unboxed
+import Data.Array.ST
+import Control.Monad.ST
+import Control.Monad
+
+main = sequence_[
+    print test1,
+    print test2,
+    print test3,
+    print test4,
+    test5,
+    test6,
+    print test7,
+    test8,
+    test0]
+
+-- helper functions
+vector l = fromList l :: Vector Double
+norm v = pnorm PNorm2 v
+
+-- hmatrix vector and matrix
+v = vector [1..10]
+m = (5><10) [1..50::Double]
+
+----------------------------------------------------------------------
+
+-- vector creation by in-place updates on a copy of the argument
+test1 = fun v
+
+fun :: Element t => Vector t -> Vector t
+fun x = runSTVector $ do
+    a <- thawVector x
+    mapM_ (flip (modifyVector a) (+57)) [0 .. dim x `div` 2 - 1]
+    return a
+
+-- another example: creation of an antidiagonal matrix from a list
+test2 = antiDiag 5 8 [1..] :: Matrix Double
+
+antiDiag :: (Element b) => Int -> Int -> [b] -> Matrix b
+antiDiag r c l = runSTMatrix $ do
+    m <- newMatrix 0 r c
+    let d = min r c - 1
+    sequence_ $ zipWith (\i v -> writeMatrix m i (c-1-i) v) [0..d] l
+    return m
+
+-- using vector or matrix functions on mutable objects requires freezing:
+test3 = g1 v
+
+g1 x = runST $ do
+    a <- thawVector x
+    writeVector a (dim x -1) 0
+    b <- freezeVector a
+    return (norm b)
+
+--  another possibility:
+test4 = g2 v
+
+g2 x = runST $ do
+    a <- thawVector x
+    writeVector a (dim x -1) 0
+    t <- liftSTVector norm a
+    return t
+
+--------------------------------------------------------------
+
+-- haskell arrays
+hv = listArray (0,9) [1..10::Double]
+hm = listArray ((0,0),(4,9)) [1..50::Double]
+
+
+
+-- conversion from standard Haskell arrays
+test5 = do
+    print $ norm (vectorFromArray hv)
+    print $ norm v
+    print $ rcond (matrixFromArray hm)
+    print $ rcond m
+
+
+-- conversion to mutable ST arrays
+test6 = do
+    let y = clearColumn m 1
+    print y
+    print (matrixFromArray y)
+
+clearColumn x c = runSTUArray $ do
+    a <- mArrayFromMatrix x
+    forM_ [0..rows x-1] $ \i->
+        writeArray a (i,c) (0::Double)
+    return a
+
+-- hmatrix functions applied to mutable ST arrays
+test7 = unitary (listArray (1,4) [3,5,7,2] :: UArray Int Double)
+
+unitary v = runSTUArray $ do
+    a <- thaw v
+    n <- norm `fmap` vectorFromMArray a
+    b <- mapArray (/n) a
+    return b
+
+-------------------------------------------------
+
+-- (just to check that they are not affected)
+test0 = do
+    print v
+    print m
+    --print hv
+    --print hm
+
+-------------------------------------------------
+
+histogram n ds = runSTVector $ do
+    h <- newVector (0::Double) n -- number of bins
+    let inc k = modifyVector h k (+1)
+    mapM_ inc ds
+    return h
+
+-- check that newVector is really called with a fresh new array
+histoCheck ds = runSTVector $ do
+    h <- newVector (0::Double) 15 -- > constant for this test
+    let inc k = modifyVector h k (+1)
+    mapM_ inc ds
+    return h
+
+hc = fromList [1 .. 15::Double]
+
+-- check that thawVector creates a new array
+histoCheck2 ds = runSTVector $ do
+    h <- thawVector hc
+    let inc k = modifyVector h k (+1)
+    mapM_ inc ds
+    return h
+
+test8 = do
+    let ds = [0..14]
+    print $ histogram 15 ds
+    print $ histogram 15 ds
+    print $ histogram 15 ds
+    print $ histoCheck ds
+    print $ histoCheck ds
+    print $ histoCheck ds
+    print $ histoCheck2 ds
+    print $ histoCheck2 ds
+    print $ histoCheck2 ds
+    putStrLn "----------------------"
diff --git a/packages/hmatrix/examples/integrate.hs b/packages/hmatrix/examples/integrate.hs
new file mode 100644
index 0000000..3a03da6
--- /dev/null
+++ b/packages/hmatrix/examples/integrate.hs
@@ -0,0 +1,24 @@
+-- Numerical integration
+import Numeric.GSL
+
+quad f a b = fst $ integrateQAGS 1E-9 100 f a b  
+
+-- A multiple integral can be easily defined using partial application
+quad2 f y1 y2 g1 g2 = quad h y1 y2
+  where
+    h y = quad (flip f y) (g1 y) (g2 y)
+
+volSphere r = 8 * quad2 (\x y -> sqrt (r*r-x*x-y*y)) 
+                        0 r (const 0) (\x->sqrt (r*r-x*x))
+
+-- wikipedia example
+exw = quad2 f 7 10 (const 11) (const 14)
+  where
+    f x y = x**2 + 4*y
+
+main = do
+    print $ quad (\x -> 4/(x^2+1)) 0 1
+    print pi
+    print $ volSphere 2.5
+    print $ 4/3*pi*2.5**3
+    print $ exw
diff --git a/packages/hmatrix/examples/kalman.hs b/packages/hmatrix/examples/kalman.hs
new file mode 100644
index 0000000..7fbe3d2
--- /dev/null
+++ b/packages/hmatrix/examples/kalman.hs
@@ -0,0 +1,51 @@
+import Numeric.LinearAlgebra
+import Graphics.Plot
+
+vector l = fromList l :: Vector Double
+matrix ls = fromLists ls :: Matrix Double
+diagl = diag . vector
+
+f = matrix [[1,0,0,0],
+            [1,1,0,0],
+            [0,0,1,0],
+            [0,0,0,1]]
+
+h = matrix [[0,-1,1,0],
+            [0,-1,0,1]]
+
+q = diagl [1,1,0,0]
+
+r = diagl [2,2]
+
+s0 = State (vector [0, 0, 10, -10]) (diagl [10,0, 100, 100])
+
+data System = System {kF, kH, kQ, kR :: Matrix Double}
+data State = State {sX :: Vector Double , sP :: Matrix Double}
+type Measurement = Vector Double
+
+kalman :: System -> State -> Measurement -> State
+kalman (System f h q r) (State x p) z = State x' p' where
+    px = f <> x                            -- prediction
+    pq = f <> p <> trans f + q             -- its covariance
+    y  = z - h <> px                       -- residue
+    cy = h <> pq <> trans h + r            -- its covariance
+    k  = pq <> trans h <> inv cy           -- kalman gain
+    x' = px + k <> y                       -- new state
+    p' = (ident (dim x) - k <> h) <> pq    -- its covariance
+
+sys = System f h q r
+
+zs = [vector [15-k,-20-k] | k <- [0..]]
+
+xs = s0 : zipWith (kalman sys) xs zs
+
+des = map (sqrt.takeDiag.sP) xs
+
+evolution n (xs,des) =
+    vector [1.. fromIntegral n]:(toColumns $ fromRows $ take n (zipWith (-) (map sX xs) des)) ++
+    (toColumns $ fromRows $ take n (zipWith (+) (map sX xs) des))
+
+main = do
+    print $ fromRows $ take 10 (map sX xs)
+    mapM_ (print . sP) $ take 10 xs
+    mplot (evolution 20 (xs,des))
diff --git a/packages/hmatrix/examples/lie.hs b/packages/hmatrix/examples/lie.hs
new file mode 100644
index 0000000..db21ea8
--- /dev/null
+++ b/packages/hmatrix/examples/lie.hs
@@ -0,0 +1,65 @@
+-- The magic of Lie Algebra
+
+import Numeric.LinearAlgebra
+
+disp = putStrLn . dispf 5
+
+rot1 :: Double -> Matrix Double
+rot1 a = (3><3)
+ [ 1, 0, 0
+ , 0, c, s
+ , 0,-s, c ]
+    where c = cos a
+          s = sin a
+
+g1,g2,g3 :: Matrix Double
+
+g1 = (3><3) [0, 0,0
+            ,0, 0,1
+            ,0,-1,0]
+
+rot2 :: Double -> Matrix Double
+rot2 a = (3><3)
+ [ c, 0, s
+ , 0, 1, 0
+ ,-s, 0, c ]
+    where c = cos a
+          s = sin a
+
+g2 = (3><3) [ 0,0,1
+            , 0,0,0
+            ,-1,0,0]
+
+rot3 :: Double -> Matrix Double
+rot3 a = (3><3)
+ [ c, s, 0
+ ,-s, c, 0
+ , 0, 0, 1 ]
+    where c = cos a
+          s = sin a
+
+g3 = (3><3) [ 0,1,0
+            ,-1,0,0
+            , 0,0,0]
+
+deg=pi/180
+
+-- commutator
+infix 8 &
+a & b = a <> b - b <> a
+
+infixl 6 |+|
+a |+| b = a + b + a&b /2  + (a-b)&(a & b) /12
+
+main = do
+   let a =  45*deg
+       b =  50*deg
+       c = -30*deg
+       exact = rot3 a <> rot1 b <> rot2 c
+       lie = scalar a * g3 |+| scalar b * g1 |+| scalar c * g2
+   putStrLn "position in the tangent space:"
+   disp lie
+   putStrLn "exponential map back to the group (2 terms):"
+   disp (expm lie)
+   putStrLn "exact position:"
+   disp exact
diff --git a/packages/hmatrix/examples/minimize.hs b/packages/hmatrix/examples/minimize.hs
new file mode 100644
index 0000000..19b2cb3
--- /dev/null
+++ b/packages/hmatrix/examples/minimize.hs
@@ -0,0 +1,50 @@
+-- the multidimensional minimization example in the GSL manual
+import Numeric.GSL
+import Numeric.LinearAlgebra
+import Graphics.Plot
+import Text.Printf(printf)
+
+-- the function to be minimized
+f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30
+
+-- exact gradient
+df [x,y] = [20*(x-1), 40*(y-2)]
+
+-- a minimization algorithm which does not require the gradient
+minimizeS f xi = minimize NMSimplex2 1E-2 100 (replicate (length xi) 1) f xi
+
+-- Numerical estimation of the gradient
+gradient f v = [partialDerivative k f v | k <- [0 .. length v -1]]
+
+partialDerivative n f v = fst (derivCentral 0.01 g (v!!n)) where
+    g x = f (concat [a,x:b])
+    (a,_:b) = splitAt n v
+
+disp = putStrLn . format "  " (printf "%.3f")
+
+allMethods :: (Enum a, Bounded a) => [a]
+allMethods = [minBound .. maxBound]
+
+test method = do
+    print method
+    let (s,p) = minimize method 1E-2 30 [1,1] f [5,7]
+    print s
+    disp p
+
+testD method = do
+    print method
+    let (s,p) = minimizeD method 1E-3 30 1E-2 1E-4 f df [5,7]
+    print s
+    disp p
+
+testD' method = do
+    putStrLn $ show method ++ " with estimated gradient"
+    let (s,p) = minimizeD method 1E-3 30 1E-2 1E-4 f (gradient f) [5,7]
+    print s
+    disp p
+
+main = do
+    mapM_ test [NMSimplex, NMSimplex2]
+    mapM_ testD allMethods
+    testD' ConjugateFR
+    mplot $ drop 3 . toColumns . snd $ minimizeS f [5,7]
diff --git a/packages/hmatrix/examples/monadic.hs b/packages/hmatrix/examples/monadic.hs
new file mode 100644
index 0000000..7c6e0f4
--- /dev/null
+++ b/packages/hmatrix/examples/monadic.hs
@@ -0,0 +1,118 @@
+-- monadic computations
+-- (contributed by Vivian McPhail)
+
+import Numeric.LinearAlgebra
+import Control.Monad.State.Strict
+import Control.Monad.Maybe
+import Foreign.Storable(Storable)
+import System.Random(randomIO)
+
+-------------------------------------------
+
+-- an instance of MonadIO, a monad transformer
+type VectorMonadT = StateT Int IO
+
+test1 :: Vector Int -> IO (Vector Int)
+test1 = mapVectorM $ \x -> do
+    putStr $ (show x) ++ " "
+    return (x + 1)
+
+-- we can have an arbitrary monad AND do IO
+addInitialM :: Vector Int -> VectorMonadT ()
+addInitialM = mapVectorM_ $ \x -> do
+    i <- get
+    liftIO $ putStr $ (show $ x + i) ++ " "
+    put $ x + i
+
+-- sum the values of the even indiced elements
+sumEvens :: Vector Int -> Int
+sumEvens = foldVectorWithIndex (\x a b -> if x `mod` 2 == 0 then a + b else b) 0
+
+-- sum and print running total of evens
+sumEvensAndPrint :: Vector Int -> VectorMonadT ()
+sumEvensAndPrint = mapVectorWithIndexM_ $ \ i x -> do
+    when (i `mod` 2 == 0) $ do
+        v <- get
+        put $ v + x
+        v' <- get
+        liftIO $ putStr $ (show v') ++ " "
+
+
+indexPlusSum :: Vector Int -> VectorMonadT ()
+indexPlusSum v' = do
+    let f i x = do
+            s <- get
+            let inc = x+s
+            liftIO $ putStr $ show (i,inc) ++ " "
+            put inc
+            return inc
+    v <- mapVectorWithIndexM f v'
+    liftIO $ do
+        putStrLn ""
+        putStrLn $ show v
+
+-------------------------------------------
+
+-- short circuit
+monoStep :: Double -> MaybeT (State Double) ()
+monoStep d = do
+    dp <- get
+    when (d < dp) (fail "negative difference")
+    put d
+{-# INLINE monoStep #-}
+
+isMonotoneIncreasing :: Vector Double -> Bool
+isMonotoneIncreasing v =
+    let res = evalState (runMaybeT $ (mapVectorM_ monoStep v)) (v @> 0)
+     in case res of
+        Nothing -> False
+        Just _  -> True
+
+
+-------------------------------------------
+
+-- | apply a test to successive elements of a vector, evaluates to true iff test passes for all pairs
+successive_ :: Storable a => (a -> a -> Bool) -> Vector a -> Bool
+successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ step (subVector 1 (dim v - 1) v))) (v @> 0)
+   where step e = do
+                  ep <- lift $ get
+                  if t e ep
+                     then lift $ put e
+                     else (fail "successive_ test failed")
+
+-- | operate on successive elements of a vector and return the resulting vector, whose length 1 less than that of the input
+successive :: (Storable a, Storable b) => (a -> a -> b) -> Vector a -> Vector b
+successive f v = evalState (mapVectorM step (subVector 1 (dim v - 1) v)) (v @> 0)
+   where step e = do
+                  ep <- get
+                  put e
+                  return $ f ep e
+
+-------------------------------------------
+
+v :: Vector Int
+v = 10 |> [0..]
+
+w = fromList ([1..10]++[10,9..1]) :: Vector Double
+
+
+main = do
+    v' <- test1 v
+    putStrLn ""
+    putStrLn $ show v'
+    evalStateT (addInitialM v) 0
+    putStrLn ""
+    putStrLn $ show (sumEvens v)
+    evalStateT (sumEvensAndPrint v) 0
+    putStrLn ""
+    evalStateT (indexPlusSum v) 0
+    putStrLn "-----------------------"
+    mapVectorM_ print v
+    print =<< (mapVectorM (const randomIO) v :: IO (Vector Double))
+    print =<< (mapVectorM (\a -> fmap (+a) randomIO) (5|>[0,100..1000]) :: IO (Vector Double))
+    putStrLn "-----------------------"
+    print $ isMonotoneIncreasing w
+    print $ isMonotoneIncreasing (subVector 0 7 w)
+    print $ successive_ (>) v
+    print $ successive_ (>) w
+    print $ successive (+) v
diff --git a/packages/hmatrix/examples/multiply.hs b/packages/hmatrix/examples/multiply.hs
new file mode 100644
index 0000000..572961c
--- /dev/null
+++ b/packages/hmatrix/examples/multiply.hs
@@ -0,0 +1,104 @@
+{-# LANGUAGE UnicodeSyntax
+           , MultiParamTypeClasses
+           , FunctionalDependencies
+           , FlexibleInstances
+           , FlexibleContexts
+--           , OverlappingInstances
+           , UndecidableInstances #-}
+
+import Numeric.LinearAlgebra
+
+class Scaling a b c | a b -> c where
+ -- ^ 0x22C5    8901    DOT OPERATOR, scaling
+ infixl 7 ⋅
+ (⋅) :: a -> b -> c
+
+instance (Num t) => Scaling t t t where
+    (⋅) = (*)
+
+instance Container Vector t => Scaling t (Vector t) (Vector t) where
+    (⋅) = scale
+
+instance Container Vector t => Scaling (Vector t) t (Vector t) where
+    (⋅) = flip scale
+
+instance Container Vector t => Scaling t (Matrix t) (Matrix t) where
+    (⋅) = scale
+
+instance Container Vector t => Scaling (Matrix t) t (Matrix t) where
+    (⋅) = flip scale
+
+
+class Mul a b c | a b -> c, a c -> b, b c -> a where
+ -- ^ 0x00D7    215     MULTIPLICATION SIGN     ×, contraction
+ infixl 7 ×
+ (×) :: a -> b -> c
+
+
+-------
+
+
+
+instance Product t => Mul (Vector t) (Vector t) t where
+    (×) = udot
+
+instance Product t => Mul (Matrix t) (Vector t) (Vector t) where
+    (×) = mXv
+
+instance Product t => Mul (Vector t) (Matrix t) (Vector t) where
+    (×) = vXm
+
+instance Product t => Mul (Matrix t) (Matrix t) (Matrix t) where
+    (×) = mXm
+
+
+--instance Scaling a b c => Contraction a b c where
+--    (×) = (⋅)
+
+--------------------------------------------------------------------------------
+
+class Outer a
+  where
+    infixl 7 ⊗
+    -- | unicode 0x2297 8855 CIRCLED TIMES      ⊗
+    --
+    -- vector outer product and matrix Kronecker product
+    (⊗) :: Product t => a t -> a t -> Matrix t
+
+instance Outer Vector where
+    (⊗) = outer
+
+instance Outer Matrix where
+    (⊗) = kronecker
+
+--------------------------------------------------------------------------------
+
+
+v = 3 |> [1..] :: Vector Double
+
+m = (3 >< 3) [1..] :: Matrix Double
+
+s = 3 :: Double
+
+a = s ⋅ v × m × m × v ⋅ s
+
+--b = (v ⊗ m) ⊗ (v ⊗ m)
+
+--c = v ⊗ m ⊗ v ⊗ m
+
+d = s ⋅ (3 |> [10,20..] :: Vector Double)
+
+u = fromList [3,0,5]
+w = konst 1 (2,3) :: Matrix Double
+
+main = do
+    print $ (scale s v <> m) `udot` v
+    print $ scale s v `udot` (m <> v)
+    print $ s * ((v <> m) `udot` v)
+    print $ s ⋅ v × m × v
+    print a
+--    print (b == c)
+    print d
+    print $ asColumn u ⊗ w
+    print $ w ⊗ asColumn u
+
diff --git a/packages/hmatrix/examples/ode.hs b/packages/hmatrix/examples/ode.hs
new file mode 100644
index 0000000..dc6e0ec
--- /dev/null
+++ b/packages/hmatrix/examples/ode.hs
@@ -0,0 +1,50 @@
+{-# LANGUAGE ViewPatterns #-}
+import Numeric.GSL.ODE
+import Numeric.LinearAlgebra
+import Graphics.Plot
+import Debug.Trace(trace)
+debug x = trace (show x) x
+
+vanderpol mu = do
+    let xdot mu t [x,v] = [v, -x + mu * v * (1-x^2)]
+        ts = linspace 1000 (0,50)
+        sol = toColumns $ odeSolve (xdot mu) [1,0] ts
+    mplot (ts : sol)
+    mplot sol
+
+
+harmonic w d = do
+    let xdot w d t [x,v] = [v, a*x + b*v] where a = -w^2; b = -2*d*w
+        ts = linspace 100 (0,20)
+        sol = odeSolve (xdot w d) [1,0] ts
+    mplot (ts : toColumns sol)
+
+
+kepler v a = mplot (take 2 $ toColumns sol) where
+    xdot t [x,y,vx,vy] = [vx,vy,x*k,y*k]
+        where g=1
+              k=(-g)*(x*x+y*y)**(-1.5)
+    ts = linspace 100 (0,30)
+    sol = odeSolve xdot [4, 0, v * cos (a*degree), v * sin (a*degree)] ts
+    degree = pi/180
+
+
+main = do
+    vanderpol 2
+    harmonic 1 0
+    harmonic 1 0.1
+    kepler 0.3 60
+    kepler 0.4 70
+    vanderpol' 2
+
+-- example of odeSolveV with jacobian
+vanderpol' mu = do
+    let xdot mu t (toList->[x,v]) = fromList [v, -x + mu * v * (1-x^2)]
+        jac t (toList->[x,v]) = (2><2) [ 0          ,          1
+                                       , -1-2*x*v*mu, mu*(1-x**2) ]
+        ts = linspace 1000 (0,50)
+        hi = (ts@>1 - ts@>0)/100
+        sol = toColumns $ odeSolveV (MSBDF jac) hi 1E-8 1E-8 (xdot mu) (fromList [1,0]) ts
+    mplot sol
+
+
diff --git a/packages/hmatrix/examples/parallel.hs b/packages/hmatrix/examples/parallel.hs
new file mode 100644
index 0000000..e875407
--- /dev/null
+++ b/packages/hmatrix/examples/parallel.hs
@@ -0,0 +1,28 @@
+-- $ ghc --make -O -rtsopts -threaded parallel.hs
+-- $ ./parallel 3000 +RTS -N4 -s -A200M
+
+import System.Environment(getArgs)
+import Numeric.LinearAlgebra
+import Control.Parallel.Strategies
+import System.Time
+
+inParallel = parMap rwhnf id
+
+-- matrix product decomposed into p parallel subtasks
+parMul p x y = fromBlocks [ inParallel ( map (x <>) ys ) ]
+    where [ys] = toBlocksEvery (rows y) (cols y `div` p) y
+
+main = do
+    n <- (read . head) `fmap` getArgs
+    let m = ident n :: Matrix Double
+    time $ print $ maxElement $ takeDiag $ m <> m
+    time $ print $ maxElement $ takeDiag $ parMul 2 m m
+    time $ print $ maxElement $ takeDiag $ parMul 4 m m
+    time $ print $ maxElement $ takeDiag $ parMul 8 m m
+
+time act = do
+    t0 <- getClockTime
+    act
+    t1 <- getClockTime
+    print $ tdSec $ normalizeTimeDiff $ diffClockTimes t1 t0
+
diff --git a/packages/hmatrix/examples/pca1.hs b/packages/hmatrix/examples/pca1.hs
new file mode 100644
index 0000000..a11eba9
--- /dev/null
+++ b/packages/hmatrix/examples/pca1.hs
@@ -0,0 +1,46 @@
+-- Principal component analysis
+
+import Numeric.LinearAlgebra
+import System.Directory(doesFileExist)
+import System.Process(system)
+import Control.Monad(when)
+
+type Vec = Vector Double
+type Mat = Matrix Double
+
+
+-- Vector with the mean value of the columns of a matrix
+mean a = constant (recip . fromIntegral . rows $ a) (rows a) <> a
+
+-- covariance matrix of a list of observations stored as rows
+cov x = (trans xc <> xc) / fromIntegral (rows x - 1)
+    where xc = x - asRow (mean x)
+
+
+-- creates the compression and decompression functions from the desired number of components
+pca :: Int -> Mat -> (Vec -> Vec , Vec -> Vec)
+pca n dataSet = (encode,decode)
+  where
+    encode x = vp <> (x - m)
+    decode x = x <> vp + m
+    m = mean dataSet
+    c = cov dataSet
+    (_,v) = eigSH' c
+    vp = takeRows n (trans v)
+
+norm = pnorm PNorm2
+
+main = do
+    ok <- doesFileExist ("mnist.txt")
+    when (not ok)  $ do
+        putStrLn "\nTrying to download test datafile..."
+        system("wget -nv http://dis.um.es/~alberto/material/sp/mnist.txt.gz")
+        system("gunzip mnist.txt.gz")
+        return ()
+    m <- loadMatrix "mnist.txt" -- fromFile "mnist.txt" (5000,785)
+    let xs = takeColumns (cols m -1) m -- the last column is the digit type (class label)
+    let x = toRows xs !! 4  -- an arbitrary test Vec
+    let (pe,pd) = pca 10 xs
+    let y = pe x
+    print y  -- compressed version
+    print $ norm (x - pd y) / norm x --reconstruction quality
diff --git a/packages/hmatrix/examples/pca2.hs b/packages/hmatrix/examples/pca2.hs
new file mode 100644
index 0000000..e7ea95f
--- /dev/null
+++ b/packages/hmatrix/examples/pca2.hs
@@ -0,0 +1,65 @@
+-- Improved PCA, including illustrative graphics
+
+import Numeric.LinearAlgebra
+import Graphics.Plot
+import System.Directory(doesFileExist)
+import System.Process(system)
+import Control.Monad(when)
+
+type Vec = Vector Double
+type Mat = Matrix Double
+
+-- Vector with the mean value of the columns of a matrix
+mean a = constant (recip . fromIntegral . rows $ a) (rows a) <> a
+
+-- covariance matrix of a list of observations stored as rows
+cov x = (trans xc <> xc) / fromIntegral (rows x - 1)
+    where xc = x - asRow (mean x)
+
+
+type Stat = (Vec, [Double], Mat)
+-- 1st and 2nd order statistics of a dataset (mean, eigenvalues and eigenvectors of cov)
+stat :: Mat -> Stat
+stat x = (m, toList s, trans v) where
+    m = mean x
+    (s,v) = eigSH' (cov x)
+
+-- creates the compression and decompression functions from the desired reconstruction
+-- quality and the statistics of a data set
+pca :: Double -> Stat -> (Vec -> Vec , Vec -> Vec)
+pca prec (m,s,v) = (encode,decode)    
+  where    
+    encode x = vp <> (x - m)
+    decode x = x <> vp + m
+    vp = takeRows n v
+    n = 1 + (length $ fst $ span (< (prec'*sum s)) $ cumSum s)
+    cumSum = tail . scanl (+) 0.0
+    prec' = if prec <=0.0 || prec >= 1.0
+                then error "the precision in pca must be 0<prec<1"
+                else prec
+
+shdigit :: Vec -> IO ()
+shdigit v = imshow (reshape 28 (-v))
+
+-- shows the effect of a given reconstruction quality on a test vector
+test :: Stat -> Double -> Vec -> IO ()
+test st prec x = do
+    let (pe,pd) = pca prec st
+    let y = pe x
+    print $ dim y
+    shdigit (pd y)
+
+main = do
+    ok <- doesFileExist ("mnist.txt")
+    when (not ok)  $ do
+        putStrLn "\nTrying to download test datafile..."
+        system("wget -nv http://dis.um.es/~alberto/material/sp/mnist.txt.gz")
+        system("gunzip mnist.txt.gz")
+        return ()
+    m <- loadMatrix "mnist.txt"
+    let xs = takeColumns (cols m -1) m
+    let x = toRows xs !! 4  -- an arbitrary test vector
+    shdigit x
+    let st = stat xs
+    test st 0.90 x
+    test st 0.50 x
diff --git a/packages/hmatrix/examples/pinv.hs b/packages/hmatrix/examples/pinv.hs
new file mode 100644
index 0000000..7de50b8
--- /dev/null
+++ b/packages/hmatrix/examples/pinv.hs
@@ -0,0 +1,20 @@
+import Numeric.LinearAlgebra
+import Graphics.Plot
+import Text.Printf(printf)
+
+expand :: Int -> Vector Double -> Matrix Double
+expand n x = fromColumns $ map (x^) [0 .. n]
+
+polynomialModel :: Vector Double -> Vector Double -> Int
+                -> (Vector Double -> Vector Double)
+polynomialModel x y n = f where
+    f z = expand n z <> ws
+    ws  = expand n x <\> y
+
+main = do
+    [x,y] <- (toColumns . readMatrix) `fmap` readFile "data.txt"
+    let pol = polynomialModel x y
+    let view = [x, y, pol 1 x, pol 2 x, pol 3 x]
+    putStrLn $ "  x      y      p 1    p 2    p 3"
+    putStrLn $ format "  " (printf "%.2f") $ fromColumns view
+    mplot view
diff --git a/packages/hmatrix/examples/plot.hs b/packages/hmatrix/examples/plot.hs
new file mode 100644
index 0000000..f950aa5
--- /dev/null
+++ b/packages/hmatrix/examples/plot.hs
@@ -0,0 +1,20 @@
+import Numeric.LinearAlgebra
+import Graphics.Plot
+import Numeric.GSL.Special(erf_Z, erf)
+
+sombrero n = f x y where 
+    (x,y) = meshdom range range
+    range = linspace n (-2,2)
+    f x y = exp (-r2) * cos (2*r2) where 
+        r2 = x*x+y*y
+
+f x =  sin x + 0.5 * sin (5*x)
+
+gaussianPDF = erf_Z
+cumdist x = 0.5 * (1+ erf (x/sqrt 2))
+
+main = do
+    let x = linspace 1000 (-4,4)
+    mplot [f x]
+    mplot [x, mapVector cumdist x, mapVector gaussianPDF x]
+    mesh (sombrero 40)
\ No newline at end of file
diff --git a/packages/hmatrix/examples/root.hs b/packages/hmatrix/examples/root.hs
new file mode 100644
index 0000000..8546ff5
--- /dev/null
+++ b/packages/hmatrix/examples/root.hs
@@ -0,0 +1,31 @@
+-- root finding examples
+import Numeric.GSL
+import Numeric.LinearAlgebra
+import Text.Printf(printf)
+
+rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ]
+
+test method = do
+    print method
+    let (s,p) = root method 1E-7 30 (rosenbrock 1 10) [-10,-5]
+    print s -- solution
+    disp p -- evolution of the algorithm
+
+jacobian a b [x,y] = [ [-a    , 0]
+                     , [-2*b*x, b] ]
+
+testJ method = do
+    print method
+    let (s,p) = rootJ method 1E-7 30 (rosenbrock 1 10) (jacobian 1 10) [-10,-5]
+    print s
+    disp p
+
+disp = putStrLn . format "  " (printf "%.3f")
+
+main = do
+    test Hybrids
+    test Hybrid
+    test DNewton
+    test Broyden
+
+    mapM_ testJ [HybridsJ .. GNewton]
diff --git a/packages/hmatrix/examples/vector.hs b/packages/hmatrix/examples/vector.hs
new file mode 100644
index 0000000..f531cbd
--- /dev/null
+++ b/packages/hmatrix/examples/vector.hs
@@ -0,0 +1,31 @@
+-- conversion to/from Data.Vector.Storable
+-- from Roman Leshchinskiy "vector" package
+--
+-- In the future Data.Packed.Vector will be replaced by Data.Vector.Storable
+
+-------------------------------------------
+
+import Numeric.LinearAlgebra as H
+import Data.Packed.Development(unsafeFromForeignPtr, unsafeToForeignPtr)
+import Foreign.Storable
+import qualified Data.Vector.Storable as V
+
+fromVector :: Storable t => V.Vector t -> H.Vector t
+fromVector v = unsafeFromForeignPtr p i n where
+    (p,i,n) = V.unsafeToForeignPtr v
+
+toVector :: Storable t => H.Vector t -> V.Vector t
+toVector v = V.unsafeFromForeignPtr p i n where
+    (p,i,n) = unsafeToForeignPtr v
+
+-------------------------------------------
+
+v = V.slice 5 10 (V.fromList [1 .. 10::Double] V.++ V.replicate 10 7)
+
+w = subVector 2 3 (linspace 5 (0,1)) :: Vector Double
+
+main = do
+    print v
+    print $ fromVector v
+    print w
+    print $ toVector w