8 files changed, 145 insertions, 48 deletions

diff --git a/CHANGES.md b/CHANGES.md
index de466e5..10be500 100644
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -3,6 +3,9 @@
 
 - integration over infinite intervals
 
+- msadams and msbdf methods for ode
+
+
 0.13.0.0
 --------
 
diff --git a/examples/Real.hs b/examples/Real.hs
index 02eba9f..9083b87 100644
--- a/examples/Real.hs
+++ b/examples/Real.hs
@@ -11,10 +11,10 @@ module Real(
     row,
     col,
     (#),(&), (//), blocks,
-    rand
+    rand, randn
 ) where
 
-import Numeric.LinearAlgebra hiding ((<>), (<|>), (<->), (<\>), (.*), (*/))
+import Numeric.LinearAlgebra hiding ((<>), (<\>))
 import System.Random(randomIO)
 
 infixl 7 <>
@@ -44,12 +44,21 @@ infixl 7 \>
 m \> v = flatten (m <\> reshape 1 v)
 
 -- | Pseudorandom matrix with uniform elements between 0 and 1.
-rand :: Int -- ^ rows
+randm :: RandDist
+     -> Int -- ^ rows
      -> Int -- ^ columns
      -> IO (Matrix Double)
-rand r c = do
+randm d r c = do
     seed <- randomIO
-    return (reshape c $ randomVector seed Uniform (r*c))
+    return (reshape c $ randomVector seed d (r*c))
+
+-- | Pseudorandom matrix with uniform elements between 0 and 1.
+rand :: Int -> Int -> IO (Matrix Double)
+rand = randm Uniform
+
+-- | Pseudorandom matrix with normal elements
+randn :: Int -> Int -> IO (Matrix Double)
+randn = randm Gaussian
 
 -- | Real identity matrix.
 eye :: Int -> Matrix Double
diff --git a/examples/ode.hs b/examples/ode.hs
index 082c46c..dc6e0ec 100644
--- a/examples/ode.hs
+++ b/examples/ode.hs
@@ -1,6 +1,9 @@
+{-# 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)]
@@ -32,3 +35,16 @@ main = do
     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/lib/Numeric/GSL/ODE.hs b/lib/Numeric/GSL/ODE.hs
index 2251acd..c243f4b 100644
--- a/lib/Numeric/GSL/ODE.hs
+++ b/lib/Numeric/GSL/ODE.hs
@@ -29,7 +29,7 @@ main = mplot (ts : toColumns sol)@
 -----------------------------------------------------------------------------
 
 module Numeric.GSL.ODE (
-    odeSolve, odeSolveV, ODEMethod(..)
+    odeSolve, odeSolveV, ODEMethod(..), Jacobian
 ) where
 
 import Data.Packed.Internal
@@ -41,18 +41,21 @@ import System.IO.Unsafe(unsafePerformIO)
 
 -------------------------------------------------------------------------
 
+type Jacobian = Double -> Vector Double -> Matrix Double
+
 -- | Stepping functions
 data ODEMethod = RK2 -- ^ Embedded Runge-Kutta (2, 3) method.
-               | RK4 -- ^ 4th order (classical) Runge-Kutta. The error estimate is obtained by halving the step-size. For more efficient estimate of the error, use 'RKf45'.
+               | RK4 -- ^ 4th order (classical) Runge-Kutta. The error estimate is obtained by halving the step-size. For more efficient estimate of the error, use the embedded methods.
                | RKf45 -- ^ Embedded Runge-Kutta-Fehlberg (4, 5) method. This method is a good general-purpose integrator.
                | RKck -- ^ Embedded Runge-Kutta Cash-Karp (4, 5) method.
                | RK8pd -- ^ Embedded Runge-Kutta Prince-Dormand (8,9) method.
-               | RK2imp -- ^ Implicit 2nd order Runge-Kutta at Gaussian points.
-               | RK4imp -- ^ Implicit 4th order Runge-Kutta at Gaussian points.
-               | BSimp -- ^ Implicit Bulirsch-Stoer method of Bader and Deuflhard. This algorithm requires the Jacobian.
-               | Gear1 -- ^ M=1 implicit Gear method.
-               | Gear2 -- ^ M=2 implicit Gear method.
-               deriving (Enum,Eq,Show,Bounded)
+               | RK2imp Jacobian -- ^ Implicit 2nd order Runge-Kutta at Gaussian points.
+               | RK4imp Jacobian -- ^ Implicit 4th order Runge-Kutta at Gaussian points.
+               | BSimp Jacobian -- ^ Implicit Bulirsch-Stoer method of Bader and Deuflhard. The method is generally suitable for stiff problems.
+               | RK1imp Jacobian -- ^ Implicit Gaussian first order Runge-Kutta. Also known as implicit Euler or backward Euler method. Error estimation is carried out by the step doubling method.
+               | MSAdams -- ^ A variable-coefficient linear multistep Adams method in Nordsieck form. This stepper uses explicit Adams-Bashforth (predictor) and implicit Adams-Moulton (corrector) methods in P(EC)^m functional iteration mode. Method order varies dynamically between 1 and 12. 
+               | MSBDF Jacobian -- ^ A variable-coefficient linear multistep backward differentiation formula (BDF) method in Nordsieck form. This stepper uses the explicit BDF formula as predictor and implicit BDF formula as corrector. A modified Newton iteration method is used to solve the system of non-linear equations. Method order varies dynamically between 1 and 5. The method is generally suitable for stiff problems.
+
 
 -- | A version of 'odeSolveV' with reasonable default parameters and system of equations defined using lists.
 odeSolve
@@ -60,7 +63,7 @@ odeSolve
     -> [Double]        -- ^ initial conditions
     -> Vector Double   -- ^ desired solution times
     -> Matrix Double   -- ^ solution
-odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) Nothing (fromList xi) ts
+odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) (fromList xi) ts
     where hi = (ts@>1 - ts@>0)/100
           epsAbs = 1.49012e-08
           epsRel = 1.49012e-08
@@ -73,11 +76,33 @@ odeSolveV
     -> Double -- ^ absolute tolerance for the state vector
     -> Double -- ^ relative tolerance for the state vector
     -> (Double -> Vector Double -> Vector Double)   -- ^ xdot(t,x)
+    -> Vector Double     -- ^ initial conditions
+    -> Vector Double     -- ^ desired solution times
+    -> Matrix Double     -- ^ solution
+odeSolveV RK2 = odeSolveV' 0 Nothing
+odeSolveV RK4 = odeSolveV' 1 Nothing
+odeSolveV RKf45 = odeSolveV' 2 Nothing
+odeSolveV RKck = odeSolveV' 3 Nothing
+odeSolveV RK8pd = odeSolveV' 4 Nothing
+odeSolveV (RK2imp jac) = odeSolveV' 5 (Just jac)
+odeSolveV (RK4imp jac) = odeSolveV' 6 (Just jac)
+odeSolveV (BSimp jac) = odeSolveV' 7 (Just jac)
+odeSolveV (RK1imp jac) = odeSolveV' 8 (Just jac)
+odeSolveV MSAdams = odeSolveV' 9 Nothing
+odeSolveV (MSBDF jac) = odeSolveV' 10 (Just jac)
+
+
+odeSolveV'
+    :: CInt
     -> Maybe (Double -> Vector Double -> Matrix Double)   -- ^ optional jacobian
+    -> Double -- ^ initial step size
+    -> Double -- ^ absolute tolerance for the state vector
+    -> Double -- ^ relative tolerance for the state vector
+    -> (Double -> Vector Double -> Vector Double)   -- ^ xdot(t,x)
     -> Vector Double     -- ^ initial conditions
     -> Vector Double     -- ^ desired solution times
     -> Matrix Double     -- ^ solution
-odeSolveV method h epsAbs epsRel f mbjac xiv ts = unsafePerformIO $ do
+odeSolveV' method mbjac h epsAbs epsRel f  xiv ts = unsafePerformIO $ do
     let n   = dim xiv
     fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t))
     jp <- case mbjac of
@@ -86,7 +111,7 @@ odeSolveV method h epsAbs epsRel f mbjac xiv ts = unsafePerformIO $ do
     sol <- vec xiv $ \xiv' ->
             vec (checkTimes ts) $ \ts' ->
              createMIO (dim ts) n
-              (ode_c (fi (fromEnum method)) h epsAbs epsRel fp jp // xiv' // ts' )
+              (ode_c (method) h epsAbs epsRel fp jp // xiv' // ts' )
               "ode"
     freeHaskellFunPtr fp
     return sol
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c
index 5c17836..7f1cf68 100644
--- a/lib/Numeric/GSL/gsl-aux.c
+++ b/lib/Numeric/GSL/gsl-aux.c
@@ -32,7 +32,7 @@
 #include <gsl/gsl_complex_math.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
-#include <gsl/gsl_odeiv.h>
+#include <gsl/gsl_odeiv2.h>
 #include <gsl/gsl_multifit_nlin.h>
 #include <string.h>
 #include <stdio.h>
@@ -1356,38 +1356,38 @@ int ode(int method, double h, double eps_abs, double eps_rel,
         int jac(double, int, const double*, int, int, double*),
         KRVEC(xi), KRVEC(ts), RMAT(sol)) {
 
-    const gsl_odeiv_step_type * T;
+    const gsl_odeiv2_step_type * T;
 
     switch(method) {
-        case 0 : {T = gsl_odeiv_step_rk2; break; }
-        case 1 : {T = gsl_odeiv_step_rk4; break; }
-        case 2 : {T = gsl_odeiv_step_rkf45; break; }
-        case 3 : {T = gsl_odeiv_step_rkck; break; }
-        case 4 : {T = gsl_odeiv_step_rk8pd; break; }
-        case 5 : {T = gsl_odeiv_step_rk2imp; break; }
-        case 6 : {T = gsl_odeiv_step_rk4imp; break; }
-        case 7 : {T = gsl_odeiv_step_bsimp; break; }
-        case 8 : {T = gsl_odeiv_step_gear1; break; }
-        case 9 : {T = gsl_odeiv_step_gear2; break; }
+        case 0 : {T = gsl_odeiv2_step_rk2; break; }
+        case 1 : {T = gsl_odeiv2_step_rk4; break; }
+        case 2 : {T = gsl_odeiv2_step_rkf45; break; }
+        case 3 : {T = gsl_odeiv2_step_rkck; break; }
+        case 4 : {T = gsl_odeiv2_step_rk8pd; break; }
+        case 5 : {T = gsl_odeiv2_step_rk2imp; break; }
+        case 6 : {T = gsl_odeiv2_step_rk4imp; break; }
+        case 7 : {T = gsl_odeiv2_step_bsimp; break; }
+        case 8 : {T = gsl_odeiv2_step_rk1imp; break; }
+        case 9 : {T = gsl_odeiv2_step_msadams; break; }
+        case 10: {T = gsl_odeiv2_step_msbdf; break; }
         default: ERROR(BAD_CODE);
     }
 
-
-    gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, xin);
-    gsl_odeiv_control * c = gsl_odeiv_control_y_new (eps_abs, eps_rel);
-    gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (xin);
-
     Tode P;
     P.f = f;
     P.j = jac;
     P.n = xin;
 
-    gsl_odeiv_system sys = {odefunc, odejac, xin, &P};
+    gsl_odeiv2_system sys = {odefunc, odejac, xin, &P};
+
+    gsl_odeiv2_driver * d =
+         gsl_odeiv2_driver_alloc_y_new (&sys, T, h, eps_abs, eps_rel);
 
     double t = tsp[0];
 
     double* y = (double*)calloc(xin,sizeof(double));
     int i,j;
+    int status;
     for(i=0; i< xin; i++) {
         y[i] = xip[i];
         solp[i] = xip[i];
@@ -1396,22 +1396,24 @@ int ode(int method, double h, double eps_abs, double eps_rel,
        for (i = 1; i < tsn ; i++)
          {
            double ti = tsp[i];
-           while (t < ti)
-             {
-               gsl_odeiv_evolve_apply (e, c, s,
-                                       &sys,
-                                       &t, ti, &h,
-                                       y);
-               // if (h < hmin) h = hmin;
-             }
+           
+           status = gsl_odeiv2_driver_apply (d, &t, ti, y);
+     
+           if (status != GSL_SUCCESS) {
+                  printf ("error in ode, return value=%d\n", status);
+                  break;
+                }
+
+//           printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
+           
            for(j=0; j<xin; j++) {
                solp[i*xin + j] = y[j];
            }
          }
 
     free(y);
-    gsl_odeiv_evolve_free (e);
-    gsl_odeiv_control_free (c);
-    gsl_odeiv_step_free (s);
-    return 0;
+    gsl_odeiv2_driver_free (d);
+    
+    return status;
 }
+
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index e2ecd4d..435cc5a 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -209,7 +209,7 @@ rightSV :: Field t => Matrix t -> (Vector Double, Matrix t)
 rightSV m | vertical m = let (_,s,v) = thinSVD m in (s,v)
           | otherwise  = let (_,s,v) = svd m     in (s,v)
 
--- | Singular values and all right singular vectors.
+-- | Singular values and all left singular vectors.
 leftSV :: Field t => Matrix t -> (Matrix t, Vector Double)
 leftSV m  | vertical m = let (u,s,_) = svd m     in (u,s)
           | otherwise  = let (u,s,_) = thinSVD m in (u,s)
diff --git a/packages/tests/hmatrix-tests.cabal b/packages/tests/hmatrix-tests.cabal
index 9f7bcdc..d1f449b 100644
--- a/packages/tests/hmatrix-tests.cabal
+++ b/packages/tests/hmatrix-tests.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix-tests
-Version:            0.1.0.0
+Version:            0.1.1
 License:            GPL
 License-file:       LICENSE
 Author:             Alberto Ruiz
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
index 8d402d0..7c41209 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
@@ -43,6 +43,8 @@ import Control.Arrow((***))
 import Debug.Trace
 import Control.Monad(when)
 
+import Data.Packed.ST
+
 import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector
                       ,sized,classify,Testable,Property
                       ,quickCheckWithResult,maxSize,stdArgs,shrink)
@@ -622,6 +624,7 @@ runBenchmarks :: IO ()
 runBenchmarks = do
     solveBench
     subBench
+    mkVecBench
     multBench
     cholBench
     svdBench
@@ -638,6 +641,14 @@ time msg act = do
     printf "%6.2f s CPU\n" $ (fromIntegral (t1 - t0) / (10^12 :: Double)) :: IO ()
     return ()
 
+timeR msg act = do
+    putStr (msg++" ")
+    t0 <- getCPUTime
+    putStr (show act)
+    t1 <- getCPUTime
+    printf "%6.2f s CPU\n" $ (fromIntegral (t1 - t0) / (10^12 :: Double)) :: IO ()
+    return ()
+
 --------------------------------
 
 manymult n = foldl1' (<>) (map rot2 angles) where
@@ -653,6 +664,37 @@ multb n = foldl1' (<>) (replicate (10^6) (ident n :: Matrix Double))
 
 --------------------------------
 
+manyvec0 xs = sum $ map (\x -> x + x**2 + x**3) xs
+manyvec1 xs = sumElements $ fromRows $ map (\x -> fromList [x,x**2,x**3]) xs
+manyvec5 xs = sumElements $ fromRows $ map (\x -> vec3 x (x**2) (x**3)) xs
+
+
+manyvec2 xs = sum $ map (\x -> sqrt(x^2 + (x**2)^2 +(x**3)^2)) xs
+manyvec3 xs = sum $ map (pnorm PNorm2 . (\x -> fromList [x,x**2,x**3])) xs
+
+manyvec4 xs = sum $ map (pnorm PNorm2 . (\x -> vec3 x (x**2) (x**3))) xs
+
+vec3 :: Double -> Double -> Double -> Vector Double
+vec3 a b c = runSTVector $ do
+    v <- newUndefinedVector 3
+    writeVector v 0 a
+    writeVector v 1 b
+    writeVector v 2 c
+    return v
+
+mkVecBench = do
+    let n = 1000000
+        xs = toList $ linspace n (0,1::Double)
+    putStr "\neval data... "; print (sum xs)
+    timeR "listproc        " $ manyvec0 xs
+    timeR "fromList matrix " $ manyvec1 xs
+    timeR "vec3 matrix     " $ manyvec5 xs
+    timeR "listproc norm   " $ manyvec2 xs
+    timeR "norm fromList   " $ manyvec3 xs
+    timeR "norm vec3       " $ manyvec4 xs
+
+--------------------------------
+
 subBench = do
     putStrLn ""
     let g = foldl1' (.) (replicate (10^5) (\v -> subVector 1 (dim v -1) v))