17 files changed, 3021 insertions, 0 deletions

diff --git a/packages/gsl/src/Graphics/Plot.hs b/packages/gsl/src/Graphics/Plot.hs
new file mode 100644
index 0000000..0ea41ac
--- /dev/null
+++ b/packages/gsl/src/Graphics/Plot.hs
@@ -0,0 +1,184 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Graphics.Plot
+-- Copyright   :  (c) Alberto Ruiz 2005-8
+-- License     :  GPL-style
+--
+-- Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+-- Stability   :  provisional
+-- Portability :  uses gnuplot and ImageMagick
+--
+-- This module is deprecated. It can be replaced by improved drawing tools
+-- available in the plot\\plot-gtk packages by Vivian McPhail or Gnuplot by Henning Thielemann.
+-----------------------------------------------------------------------------
+{-# OPTIONS_HADDOCK hide #-}
+
+module Graphics.Plot(
+
+    mplot,
+
+    plot, parametricPlot, 
+
+    splot, mesh, meshdom,
+
+    matrixToPGM, imshow,
+
+    gnuplotX, gnuplotpdf, gnuplotWin
+
+) where
+
+import Numeric.Container
+import Data.List(intersperse)
+import System.Process (system)
+
+-- | From vectors x and y, it generates a pair of matrices to be used as x and y arguments for matrix functions.
+meshdom :: Vector Double -> Vector Double -> (Matrix Double , Matrix Double)
+meshdom r1 r2 = (outer r1 (constant 1 (dim r2)), outer (constant 1 (dim r1)) r2)
+
+
+{- | Draws a 3D surface representation of a real matrix.
+
+> > mesh $ build (10,10) (\\i j -> i + (j-5)^2)
+
+In certain versions you can interactively rotate the graphic using the mouse.
+
+-}
+mesh :: Matrix Double -> IO ()
+mesh m = gnuplotX (command++dat) where
+    command = "splot "++datafollows++" matrix with lines\n"
+    dat = prep $ toLists m
+
+{- | Draws the surface represented by the function f in the desired ranges and number of points, internally using 'mesh'.
+
+> > let f x y = cos (x + y) 
+> > splot f (0,pi) (0,2*pi) 50    
+
+-}
+splot :: (Matrix Double->Matrix Double->Matrix Double) -> (Double,Double) -> (Double,Double) -> Int -> IO () 
+splot f rx ry n = mesh z where
+    (x,y) = meshdom (linspace n rx) (linspace n ry)
+    z = f x y
+
+{- | plots several vectors against the first one 
+
+> > let t = linspace 100 (-3,3) in mplot [t, sin t, exp (-t^2)]
+
+-}
+mplot :: [Vector Double] -> IO ()
+mplot m = gnuplotX (commands++dats) where
+    commands = if length m == 1 then command1 else commandmore
+    command1 = "plot "++datafollows++" with lines\n" ++ dat
+    commandmore = "plot " ++ plots ++ "\n"
+    plots = concat $ intersperse ", " (map cmd [2 .. length m])
+    cmd k = datafollows++" using 1:"++show k++" with lines"
+    dat = prep $ toLists $ fromColumns m
+    dats = concat (replicate (length m-1) dat)
+
+
+{- | Draws a list of functions over a desired range and with a desired number of points 
+
+> > plot [sin, cos, sin.(3*)] (0,2*pi) 1000
+
+-}
+plot :: [Vector Double->Vector Double] -> (Double,Double) -> Int -> IO ()
+plot fs rx n = mplot (x: mapf fs x)
+    where x = linspace n rx
+          mapf gs y = map ($ y) gs
+
+{- | Draws a parametric curve. For instance, to draw a spiral we can do something like:
+
+> > parametricPlot (\t->(t * sin t, t * cos t)) (0,10*pi) 1000
+
+-}
+parametricPlot :: (Vector Double->(Vector Double,Vector Double)) -> (Double, Double) -> Int -> IO ()
+parametricPlot f rt n = mplot [fx, fy]
+    where t = linspace n rt
+          (fx,fy) = f t
+
+
+-- | writes a matrix to pgm image file
+matrixToPGM :: Matrix Double -> String
+matrixToPGM m = header ++ unlines (map unwords ll) where
+    c = cols m
+    r = rows m
+    header = "P2 "++show c++" "++show r++" "++show (round maxgray :: Int)++"\n"
+    maxgray = 255.0
+    maxval = maxElement m
+    minval = minElement m
+    scale' = if maxval == minval
+        then 0.0
+        else maxgray / (maxval - minval)
+    f x = show ( round ( scale' *(x - minval) ) :: Int )
+    ll = map (map f) (toLists m)
+
+-- | imshow shows a representation of a matrix as a gray level image using ImageMagick's display.
+imshow :: Matrix Double -> IO ()
+imshow m = do
+    _ <- system $ "echo \""++ matrixToPGM m ++"\"| display -antialias -resize 300 - &"
+    return ()
+
+----------------------------------------------------
+
+gnuplotX :: String -> IO ()
+gnuplotX command = do { _ <- system cmdstr; return()} where
+    cmdstr = "echo \""++command++"\" | gnuplot -persist"
+
+datafollows = "\\\"-\\\""
+
+prep = (++"e\n\n") . unlines . map (unwords . map show)
+
+
+gnuplotpdf :: String -> String -> [([[Double]], String)] -> IO ()
+gnuplotpdf title command ds = gnuplot (prelude ++ command ++" "++ draw) >> postproc where
+    prelude = "set terminal epslatex color; set output '"++title++".tex';"
+    (dats,defs) = unzip ds
+    draw = concat (intersperse ", " (map ("\"-\" "++) defs)) ++ "\n" ++
+           concatMap pr dats
+    postproc = do
+        _ <- system $ "epstopdf "++title++".eps"
+        mklatex
+        _ <- system $ "pdflatex "++title++"aux.tex > /dev/null"
+        _ <- system $ "pdfcrop "++title++"aux.pdf > /dev/null"
+        _ <- system $ "mv "++title++"aux-crop.pdf "++title++".pdf"
+        _ <- system $ "rm "++title++"aux.* "++title++".eps "++title++".tex"
+        return ()
+
+    mklatex = writeFile (title++"aux.tex") $
+       "\\documentclass{article}\n"++
+       "\\usepackage{graphics}\n"++
+       "\\usepackage{nopageno}\n"++
+       "\\usepackage{txfonts}\n"++
+       "\\renewcommand{\\familydefault}{phv}\n"++
+       "\\usepackage[usenames]{color}\n"++
+
+       "\\begin{document}\n"++
+
+       "\\begin{center}\n"++
+       "  \\input{./"++title++".tex}\n"++
+       "\\end{center}\n"++
+
+       "\\end{document}"
+
+    pr = (++"e\n") . unlines . map (unwords . map show)
+
+    gnuplot cmd = do
+        writeFile "gnuplotcommand" cmd
+        _ <- system "gnuplot gnuplotcommand"
+        _ <- system "rm gnuplotcommand"
+        return ()
+
+gnuplotWin :: String -> String -> [([[Double]], String)] -> IO ()
+gnuplotWin title command ds = gnuplot (prelude ++ command ++" "++ draw) where
+    (dats,defs) = unzip ds
+    draw = concat (intersperse ", " (map ("\"-\" "++) defs)) ++ "\n" ++
+           concatMap pr dats
+
+    pr = (++"e\n") . unlines . map (unwords . map show)
+
+    prelude = "set title \""++title++"\";"
+
+    gnuplot cmd = do
+        writeFile "gnuplotcommand" cmd
+        _ <- system "gnuplot -persist gnuplotcommand"
+        _ <- system "rm gnuplotcommand"
+        return ()
diff --git a/packages/gsl/src/Numeric/GSL.hs b/packages/gsl/src/Numeric/GSL.hs
new file mode 100644
index 0000000..61b8d7b
--- /dev/null
+++ b/packages/gsl/src/Numeric/GSL.hs
@@ -0,0 +1,43 @@
+{- |
+
+Module      :  Numeric.GSL
+Copyright   :  (c) Alberto Ruiz 2006-14
+License     :  GPL
+
+Maintainer  :  Alberto Ruiz
+Stability   :  provisional
+
+This module reexports all available GSL functions.
+
+The GSL special functions are in the separate package hmatrix-special.
+
+-}
+
+module Numeric.GSL (
+  module Numeric.GSL.Integration
+, module Numeric.GSL.Differentiation
+, module Numeric.GSL.Fourier
+, module Numeric.GSL.Polynomials
+, module Numeric.GSL.Minimization
+, module Numeric.GSL.Root
+, module Numeric.GSL.ODE
+, module Numeric.GSL.Fitting
+, module Data.Complex
+, setErrorHandlerOff
+) where
+
+import Numeric.GSL.Integration
+import Numeric.GSL.Differentiation
+import Numeric.GSL.Fourier
+import Numeric.GSL.Polynomials
+import Numeric.GSL.Minimization
+import Numeric.GSL.Root
+import Numeric.GSL.ODE
+import Numeric.GSL.Fitting
+import Data.Complex
+
+
+-- | This action removes the GSL default error handler (which aborts the program), so that
+-- GSL errors can be handled by Haskell (using Control.Exception) and ghci doesn't abort.
+foreign import ccall unsafe "GSL/gsl-aux.h no_abort_on_error" setErrorHandlerOff :: IO ()
+
diff --git a/packages/gsl/src/Numeric/GSL/Differentiation.hs b/packages/gsl/src/Numeric/GSL/Differentiation.hs
new file mode 100644
index 0000000..0fb58ef
--- /dev/null
+++ b/packages/gsl/src/Numeric/GSL/Differentiation.hs
@@ -0,0 +1,85 @@
+{- |
+Module      :  Numeric.GSL.Differentiation
+Copyright   :  (c) Alberto Ruiz 2006
+License     :  GPL
+
+Maintainer  :  Alberto Ruiz
+Stability   :  provisional
+
+Numerical differentiation.
+
+<http://www.gnu.org/software/gsl/manual/html_node/Numerical-Differentiation.html#Numerical-Differentiation>
+
+From the GSL manual: \"The functions described in this chapter compute numerical derivatives by finite differencing. An adaptive algorithm is used to find the best choice of finite difference and to estimate the error in the derivative.\"
+-}
+
+
+module Numeric.GSL.Differentiation (
+    derivCentral,
+    derivForward,
+    derivBackward
+) where
+
+import Foreign.C.Types
+import Foreign.Marshal.Alloc(malloc, free)
+import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr)
+import Foreign.Storable(peek)
+import Numeric.GSL.Internal
+import System.IO.Unsafe(unsafePerformIO)
+
+derivGen ::
+    CInt                   -- ^ type: 0 central, 1 forward, 2 backward
+    -> Double             -- ^ initial step size
+    -> (Double -> Double) -- ^ function
+    -> Double             -- ^ point where the derivative is taken
+    -> (Double, Double)   -- ^ result and error
+derivGen c h f x = unsafePerformIO $ do
+    r <- malloc
+    e <- malloc
+    fp <- mkfun (\y _ -> f y) 
+    c_deriv c fp x h r e // check "deriv"
+    vr <- peek r
+    ve <- peek e
+    let result = (vr,ve)
+    free r
+    free e
+    freeHaskellFunPtr fp
+    return result
+
+foreign import ccall safe "gsl-aux.h deriv" 
+ c_deriv :: CInt -> FunPtr (Double -> Ptr () -> Double) -> Double -> Double 
+                    -> Ptr Double -> Ptr Double -> IO CInt
+
+
+{- | Adaptive central difference algorithm, /gsl_deriv_central/. For example:
+
+>>> let deriv = derivCentral 0.01
+>>> deriv sin (pi/4)
+(0.7071067812000676,1.0600063101654055e-10)
+>>> cos (pi/4)
+0.7071067811865476
+
+-}
+derivCentral :: Double                  -- ^ initial step size
+                -> (Double -> Double)   -- ^ function 
+                -> Double               -- ^ point where the derivative is taken
+                -> (Double, Double)     -- ^ result and absolute error
+derivCentral = derivGen 0
+
+-- | Adaptive forward difference algorithm, /gsl_deriv_forward/. The function is evaluated only at points greater than x, and never at x itself. The derivative is returned in result and an estimate of its absolute error is returned in abserr. This function should be used if f(x) has a discontinuity at x, or is undefined for values less than x. A backward derivative can be obtained using a negative step.
+derivForward :: Double                  -- ^ initial step size
+                -> (Double -> Double)   -- ^ function 
+                -> Double               -- ^ point where the derivative is taken
+                -> (Double, Double)     -- ^ result and absolute error
+derivForward = derivGen 1
+
+-- | Adaptive backward difference algorithm, /gsl_deriv_backward/. 
+derivBackward ::Double                  -- ^ initial step size
+                -> (Double -> Double)   -- ^ function 
+                -> Double               -- ^ point where the derivative is taken
+                -> (Double, Double)     -- ^ result and absolute error
+derivBackward = derivGen 2
+
+{- | conversion of Haskell functions into function pointers that can be used in the C side
+-}
+foreign import ccall safe "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double)) 
diff --git a/packages/gsl/src/Numeric/GSL/Fitting.hs b/packages/gsl/src/Numeric/GSL/Fitting.hs
new file mode 100644
index 0000000..93fb281
--- /dev/null
+++ b/packages/gsl/src/Numeric/GSL/Fitting.hs
@@ -0,0 +1,180 @@
+{- |
+Module      :  Numeric.GSL.Fitting
+Copyright   :  (c) Alberto Ruiz 2010
+License     :  GPL
+
+Maintainer  :  Alberto Ruiz
+Stability   :  provisional
+
+Nonlinear Least-Squares Fitting
+
+<http://www.gnu.org/software/gsl/manual/html_node/Nonlinear-Least_002dSquares-Fitting.html>
+
+The example program in the GSL manual (see examples/fitting.hs):
+
+@
+dat = [
+ ([0.0],([6.0133918608118675],0.1)),
+ ([1.0],([5.5153769909966535],0.1)),
+ ([2.0],([5.261094606015287],0.1)),
+ ...
+ ([39.0],([1.0619821710802808],0.1))]
+
+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]
+@
+
+>>> path
+(6><5)
+ [ 1.0,  76.45780563978782, 1.6465931240727802, 1.8147715267618197e-2, 0.6465931240727797
+ , 2.0, 37.683816318260355,  2.858760367632973,  8.092094813253975e-2, 1.4479636296208662
+ , 3.0,    9.5807893736187,  4.948995119561291,   0.11942927999921617, 1.0945766509238248
+ , 4.0,  5.630494933603935,  5.021755718065913,   0.10287787128056883, 1.0338835440862608
+ , 5.0,  5.443976278682909,  5.045204331329302,   0.10405523433131504,  1.019416067207375
+ , 6.0, 5.4439736648994685,  5.045357818922331,   0.10404905846029407, 1.0192487112786812 ]
+>>> sol
+[(5.045357818922331,6.027976702418132e-2),
+(0.10404905846029407,3.157045047172834e-3),
+(1.0192487112786812,3.782067731353722e-2)]
+
+-}
+
+
+module Numeric.GSL.Fitting (
+    -- * Levenberg-Marquardt
+    nlFitting, FittingMethod(..),
+    -- * Utilities
+    fitModelScaled, fitModel
+) where
+
+import Numeric.LinearAlgebra
+import Numeric.GSL.Internal
+
+import Foreign.Ptr(FunPtr, freeHaskellFunPtr)
+import Foreign.C.Types
+import System.IO.Unsafe(unsafePerformIO)
+
+-------------------------------------------------------------------------
+
+type TVV = TV (TV Res)
+type TVM = TV (TM Res)
+
+data FittingMethod = LevenbergMarquardtScaled -- ^ Interface to gsl_multifit_fdfsolver_lmsder. This is a robust and efficient version of the Levenberg-Marquardt algorithm as implemented in the scaled lmder routine in minpack. Minpack was written by Jorge J. More, Burton S. Garbow and Kenneth E. Hillstrom.
+                   | LevenbergMarquardt -- ^ This is an unscaled version of the lmder algorithm. The elements of the diagonal scaling matrix D are set to 1. This algorithm may be useful in circumstances where the scaled version of lmder converges too slowly, or the function is already scaled appropriately.
+        deriving (Enum,Eq,Show,Bounded)
+
+
+-- | Nonlinear multidimensional least-squares fitting.
+nlFitting :: FittingMethod
+      -> Double                     -- ^ absolute tolerance
+      -> Double                     -- ^ relative tolerance
+      -> Int                        -- ^ maximum number of iterations allowed
+      -> (Vector Double -> Vector Double)     -- ^ function to be minimized
+      -> (Vector Double -> Matrix Double)   -- ^ Jacobian
+      -> Vector Double                   -- ^ starting point
+      -> (Vector Double, Matrix Double)  -- ^ solution vector and optimization path
+
+nlFitting method epsabs epsrel maxit fun jac xinit = nlFitGen (fi (fromEnum method)) fun jac xinit epsabs epsrel maxit
+
+nlFitGen m f jac xiv epsabs epsrel maxit = unsafePerformIO $ do
+    let p   = dim xiv
+        n   = dim (f xiv)
+    fp <- mkVecVecfun (aux_vTov (checkdim1 n p . f))
+    jp <- mkVecMatfun (aux_vTom (checkdim2 n p . jac))
+    rawpath <- createMatrix RowMajor maxit (2+p)
+    app2 (c_nlfit m fp jp epsabs epsrel (fi maxit) (fi n)) vec xiv mat rawpath "c_nlfit"
+    let it = round (rawpath @@> (maxit-1,0))
+        path = takeRows it rawpath
+        [sol] = toRows $ dropRows (it-1) path
+    freeHaskellFunPtr fp
+    freeHaskellFunPtr jp
+    return (subVector 2 p sol, path)
+
+foreign import ccall safe "nlfit"
+    c_nlfit:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> Double -> CInt -> CInt -> TVM
+
+-------------------------------------------------------
+
+checkdim1 n _p v
+    | dim v == n = v
+    | otherwise = error $ "Error: "++ show n
+                        ++ " components expected in the result of the function supplied to nlFitting"
+
+checkdim2 n p m
+    | rows m == n && cols m == p = m
+    | otherwise = error $ "Error: "++ show n ++ "x" ++ show p
+                        ++ " Jacobian expected in nlFitting"
+
+------------------------------------------------------------
+
+err (model,deriv) dat vsol = zip sol errs where
+    sol = toList vsol
+    c = max 1 (chi/sqrt (fromIntegral dof))
+    dof = length dat - (rows cov)
+    chi = norm2 (fromList $ cost (resMs model) dat sol)
+    js = fromLists $ jacobian (resDs deriv) dat sol
+    cov = inv $ trans js <.> js
+    errs = toList $ scalar c * sqrt (takeDiag cov)
+
+
+
+-- | Higher level interface to 'nlFitting' 'LevenbergMarquardtScaled'. The optimization function and
+-- Jacobian are automatically built from a model f vs x = y and its derivatives, and a list of
+-- instances (x, (y,sigma)) to be fitted.
+
+fitModelScaled
+         :: Double -- ^ absolute tolerance
+         -> Double -- ^ relative tolerance
+         -> Int    -- ^ maximum number of iterations allowed
+         -> ([Double] -> x -> [Double], [Double] -> x -> [[Double]]) -- ^ (model, derivatives)
+         -> [(x, ([Double], Double))] -- ^ instances
+         -> [Double] -- ^ starting point
+         -> ([(Double, Double)], Matrix Double) -- ^ (solution, error) and optimization path
+fitModelScaled epsabs epsrel maxit (model,deriv) dt xin = (err (model,deriv) dt sol, path) where
+    (sol,path) = nlFitting LevenbergMarquardtScaled epsabs epsrel maxit
+                (fromList . cost (resMs model) dt . toList)
+                (fromLists . jacobian (resDs deriv) dt . toList)
+                (fromList xin)
+
+
+
+-- | Higher level interface to 'nlFitting' 'LevenbergMarquardt'. The optimization function and
+-- Jacobian are automatically built from a model f vs x = y and its derivatives, and a list of
+-- instances (x,y) to be fitted.
+
+fitModel :: Double -- ^ absolute tolerance
+         -> Double -- ^ relative tolerance
+         -> Int    -- ^ maximum number of iterations allowed
+         -> ([Double] -> x -> [Double], [Double] -> x -> [[Double]]) -- ^ (model, derivatives)
+         -> [(x, [Double])]  -- ^ instances
+         -> [Double] -- ^ starting point
+         -> ([Double], Matrix Double) -- ^ solution and optimization path
+fitModel epsabs epsrel maxit (model,deriv) dt xin = (toList sol, path) where
+    (sol,path) = nlFitting LevenbergMarquardt epsabs epsrel maxit
+                (fromList . cost (resM model) dt . toList)
+                (fromLists . jacobian (resD deriv) dt . toList)
+                (fromList xin)
+
+cost model ds vs = concatMap (model vs) ds
+
+jacobian modelDer ds vs = concatMap (modelDer vs) ds
+
+-- | Model-to-residual for association pairs with sigma, to be used with 'fitModel'.
+resMs :: ([Double] -> x -> [Double]) -> [Double] -> (x, ([Double], Double)) -> [Double]
+resMs m v = \(x,(ys,s)) -> zipWith (g s) (m v x) ys where g s a b = (a-b)/s
+
+-- | Associated derivative for 'resMs'.
+resDs :: ([Double] -> x -> [[Double]]) -> [Double] -> (x, ([Double], Double)) -> [[Double]]
+resDs m v = \(x,(_,s)) -> map (map (/s)) (m v x)
+
+-- | Model-to-residual for association pairs, to be used with 'fitModel'. It is equivalent
+-- to 'resMs' with all sigmas = 1.
+resM :: ([Double] -> x -> [Double]) -> [Double] -> (x, [Double]) -> [Double]
+resM m v = \(x,ys) -> zipWith g (m v x) ys where g a b = a-b
+
+-- | Associated derivative for 'resM'.
+resD :: ([Double] -> x -> [[Double]]) -> [Double] -> (x, [Double]) -> [[Double]]
+resD m v = \(x,_) -> m v x
diff --git a/packages/gsl/src/Numeric/GSL/Fourier.hs b/packages/gsl/src/Numeric/GSL/Fourier.hs
new file mode 100644
index 0000000..734325b
--- /dev/null
+++ b/packages/gsl/src/Numeric/GSL/Fourier.hs
@@ -0,0 +1,44 @@
+{- |
+Module      : Numeric.GSL.Fourier
+Copyright   :  (c) Alberto Ruiz 2006
+License     :  GPL
+Maintainer  :  Alberto Ruiz
+Stability   :  provisional
+
+Fourier Transform.
+
+<http://www.gnu.org/software/gsl/manual/html_node/Fast-Fourier-Transforms.html#Fast-Fourier-Transforms>
+
+-}
+
+module Numeric.GSL.Fourier (
+    fft,
+    ifft
+) where
+
+import Data.Packed
+import Numeric.GSL.Internal
+import Data.Complex
+import Foreign.C.Types
+import System.IO.Unsafe (unsafePerformIO)
+
+genfft code v = unsafePerformIO $ do
+    r <- createVector (dim v)
+    app2 (c_fft code) vec v vec r "fft"
+    return r
+
+foreign import ccall unsafe "gsl-aux.h fft" c_fft ::  CInt -> TCV (TCV Res)
+
+
+{- | Fast 1D Fourier transform of a 'Vector' @(@'Complex' 'Double'@)@ using /gsl_fft_complex_forward/. It uses the same scaling conventions as GNU Octave.
+
+>>> fft (fromList [1,2,3,4])
+fromList [10.0 :+ 0.0,(-2.0) :+ 2.0,(-2.0) :+ 0.0,(-2.0) :+ (-2.0)]
+
+-}
+fft :: Vector (Complex Double) -> Vector (Complex Double)
+fft = genfft 0
+
+-- | The inverse of 'fft', using /gsl_fft_complex_inverse/.
+ifft :: Vector (Complex Double) -> Vector (Complex Double)
+ifft = genfft 1
diff --git a/packages/gsl/src/Numeric/GSL/IO.hs b/packages/gsl/src/Numeric/GSL/IO.hs
new file mode 100644
index 0000000..0d6031a
--- /dev/null
+++ b/packages/gsl/src/Numeric/GSL/IO.hs
@@ -0,0 +1,42 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.GSL.IO
+-- Copyright   :  (c) Alberto Ruiz 2007-14
+-- License     :  GPL
+-- Maintainer  :  Alberto Ruiz
+-- Stability   :  provisional
+--
+-----------------------------------------------------------------------------
+
+module Numeric.GSL.IO (
+    saveMatrix,
+    fwriteVector, freadVector, fprintfVector, fscanfVector,
+    fileDimensions, loadMatrix, fromFile
+) where
+
+import Data.Packed
+import Numeric.GSL.Vector
+import System.Process(readProcess)
+
+
+{- |  obtains the number of rows and columns in an ASCII data file
+      (provisionally using unix's wc).
+-}
+fileDimensions :: FilePath -> IO (Int,Int)
+fileDimensions fname = do
+    wcres <- readProcess "wc" ["-w",fname] ""
+    contents <- readFile fname
+    let tot = read . head . words $ wcres
+        c   = length . head . dropWhile null . map words . lines $ contents
+    if tot > 0
+        then return (tot `div` c, c)
+        else return (0,0)
+
+-- | Loads a matrix from an ASCII file formatted as a 2D table.
+loadMatrix :: FilePath -> IO (Matrix Double)
+loadMatrix file = fromFile file =<< fileDimensions file
+
+-- | Loads a matrix from an ASCII file (the number of rows and columns must be known in advance).
+fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double)
+fromFile filename (r,c) = reshape c `fmap` fscanfVector filename (r*c)
+
diff --git a/packages/gsl/src/Numeric/GSL/Integration.hs b/packages/gsl/src/Numeric/GSL/Integration.hs
new file mode 100644
index 0000000..9c1d43a
--- /dev/null
+++ b/packages/gsl/src/Numeric/GSL/Integration.hs
@@ -0,0 +1,246 @@
+{- |
+Module      :  Numeric.GSL.Integration
+Copyright   :  (c) Alberto Ruiz 2006
+License     :  GPL
+Maintainer  :  Alberto Ruiz
+Stability   :  provisional
+
+Numerical integration routines.
+
+<http://www.gnu.org/software/gsl/manual/html_node/Numerical-Integration.html#Numerical-Integration>
+-}
+
+
+module Numeric.GSL.Integration (
+    integrateQNG,
+    integrateQAGS,
+    integrateQAGI,
+    integrateQAGIU,
+    integrateQAGIL,
+    integrateCQUAD
+) where
+
+import Foreign.C.Types
+import Foreign.Marshal.Alloc(malloc, free)
+import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr)
+import Foreign.Storable(peek)
+import Numeric.GSL.Internal
+import System.IO.Unsafe(unsafePerformIO)
+
+eps = 1e-12
+
+{- | conversion of Haskell functions into function pointers that can be used in the C side
+-}
+foreign import ccall safe "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double))
+
+--------------------------------------------------------------------
+{- | Numerical integration using /gsl_integration_qags/ (adaptive integration with singularities). For example:
+
+>>> let quad = integrateQAGS 1E-9 1000
+>>> let f a x = x**(-0.5) * log (a*x)
+>>> quad (f 1) 0 1
+(-3.999999999999974,4.871658632055187e-13)
+
+-}
+
+integrateQAGS :: Double               -- ^ precision (e.g. 1E-9)
+                 -> Int               -- ^ size of auxiliary workspace (e.g. 1000)
+                 -> (Double -> Double) -- ^ function to be integrated on the interval (a,b)
+                 -> Double           -- ^ a
+                 -> Double           -- ^ b
+                 -> (Double, Double) -- ^ result of the integration and error
+integrateQAGS prec n f a b = unsafePerformIO $ do
+    r <- malloc
+    e <- malloc
+    fp <- mkfun (\x _ -> f x)
+    c_integrate_qags fp a b eps prec (fromIntegral n) r e // check "integrate_qags"
+    vr <- peek r
+    ve <- peek e
+    let result = (vr,ve)
+    free r
+    free e
+    freeHaskellFunPtr fp
+    return result
+
+foreign import ccall safe "integrate_qags" c_integrate_qags
+    :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
+    -> Double -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt
+
+-----------------------------------------------------------------
+{- | Numerical integration using /gsl_integration_qng/ (useful for fast integration of smooth functions). For example:
+
+>>> let quad = integrateQNG 1E-6
+>>> quad (\x -> 4/(1+x*x)) 0 1
+(3.141592653589793,3.487868498008632e-14)
+
+-}
+
+integrateQNG :: Double               -- ^ precision (e.g. 1E-9)
+                 -> (Double -> Double) -- ^ function to be integrated on the interval (a,b)
+                 -> Double           -- ^ a
+                 -> Double           -- ^ b
+                 -> (Double, Double) -- ^ result of the integration and error
+integrateQNG prec f a b = unsafePerformIO $ do
+    r <- malloc
+    e <- malloc
+    fp <- mkfun (\x _ -> f x)
+    c_integrate_qng fp a b eps prec r e  // check "integrate_qng"
+    vr <- peek r
+    ve <- peek e
+    let result = (vr,ve)
+    free r
+    free e
+    freeHaskellFunPtr fp
+    return result
+
+foreign import ccall safe "integrate_qng" c_integrate_qng
+    :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
+    -> Double -> Double -> Ptr Double -> Ptr Double -> IO CInt
+
+--------------------------------------------------------------------
+{- | Numerical integration using /gsl_integration_qagi/ (integration over the infinite integral -Inf..Inf using QAGS).
+For example:
+
+>>> let quad = integrateQAGI 1E-9 1000
+>>> let f a x = exp(-a * x^2)
+>>> quad (f 0.5)
+(2.5066282746310002,6.229215880648858e-11)
+
+-}
+
+integrateQAGI :: Double               -- ^ precision (e.g. 1E-9)
+                 -> Int               -- ^ size of auxiliary workspace (e.g. 1000)
+                 -> (Double -> Double) -- ^ function to be integrated on the interval (-Inf,Inf)
+                 -> (Double, Double) -- ^ result of the integration and error
+integrateQAGI prec n f = unsafePerformIO $ do
+    r <- malloc
+    e <- malloc
+    fp <- mkfun (\x _ -> f x)
+    c_integrate_qagi fp eps prec (fromIntegral n) r e // check "integrate_qagi"
+    vr <- peek r
+    ve <- peek e
+    let result = (vr,ve)
+    free r
+    free e
+    freeHaskellFunPtr fp
+    return result
+
+foreign import ccall safe "integrate_qagi" c_integrate_qagi
+    :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
+    -> CInt -> Ptr Double -> Ptr Double -> IO CInt
+
+--------------------------------------------------------------------
+{- | Numerical integration using /gsl_integration_qagiu/ (integration over the semi-infinite integral a..Inf).
+For example:
+
+>>> let quad = integrateQAGIU 1E-9 1000
+>>> let f a x = exp(-a * x^2)
+>>> quad (f 0.5) 0
+(1.2533141373155001,3.114607940324429e-11)
+
+-}
+
+integrateQAGIU :: Double               -- ^ precision (e.g. 1E-9)
+                 -> Int               -- ^ size of auxiliary workspace (e.g. 1000)
+                 -> (Double -> Double) -- ^ function to be integrated on the interval (a,Inf)
+                 -> Double             -- ^ a
+                 -> (Double, Double) -- ^ result of the integration and error
+integrateQAGIU prec n f a = unsafePerformIO $ do
+    r <- malloc
+    e <- malloc
+    fp <- mkfun (\x _ -> f x)
+    c_integrate_qagiu fp a eps prec (fromIntegral n) r e // check "integrate_qagiu"
+    vr <- peek r
+    ve <- peek e
+    let result = (vr,ve)
+    free r
+    free e
+    freeHaskellFunPtr fp
+    return result
+
+foreign import ccall safe "integrate_qagiu" c_integrate_qagiu
+    :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
+    -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt
+
+--------------------------------------------------------------------
+{- | Numerical integration using /gsl_integration_qagil/ (integration over the semi-infinite integral -Inf..b).
+For example:
+
+>>> let quad = integrateQAGIL 1E-9 1000
+>>> let f a x = exp(-a * x^2)
+>>> quad (f 0.5) 0
+(1.2533141373155001,3.114607940324429e-11)
+
+-}
+
+integrateQAGIL :: Double               -- ^ precision (e.g. 1E-9)
+                 -> Int               -- ^ size of auxiliary workspace (e.g. 1000)
+                 -> (Double -> Double) -- ^ function to be integrated on the interval (a,Inf)
+                 -> Double             -- ^ b
+                 -> (Double, Double) -- ^ result of the integration and error
+integrateQAGIL prec n f b = unsafePerformIO $ do
+    r <- malloc
+    e <- malloc
+    fp <- mkfun (\x _ -> f x)
+    c_integrate_qagil fp b eps prec (fromIntegral n) r e // check "integrate_qagil"
+    vr <- peek r
+    ve <- peek e
+    let result = (vr,ve)
+    free r
+    free e
+    freeHaskellFunPtr fp
+    return result
+
+foreign import ccall safe "gsl-aux.h integrate_qagil" c_integrate_qagil
+    :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
+    -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt
+
+
+--------------------------------------------------------------------
+{- | Numerical integration using /gsl_integration_cquad/ (quadrature
+for general integrands).  From the GSL manual:
+
+@CQUAD is a new doubly-adaptive general-purpose quadrature routine
+which can handle most types of singularities, non-numerical function
+values such as Inf or NaN, as well as some divergent integrals. It
+generally requires more function evaluations than the integration
+routines in QUADPACK, yet fails less often for difficult integrands.@
+
+For example:
+
+>>> let quad = integrateCQUAD 1E-12 1000
+>>> let f a x = exp(-a * x^2)
+>>> quad (f 0.5) 2 5
+(5.7025405463957006e-2,9.678874441303705e-16,95)
+
+Unlike other quadrature methods, integrateCQUAD also returns the
+number of function evaluations required.
+
+-}
+
+integrateCQUAD :: Double               -- ^ precision (e.g. 1E-9)
+                 -> Int               -- ^ size of auxiliary workspace (e.g. 1000)
+                 -> (Double -> Double) -- ^ function to be integrated on the interval (a, b)
+                 -> Double             -- ^ a
+                 -> Double             -- ^ b
+                 -> (Double, Double, Int) -- ^ result of the integration, error and number of function evaluations performed
+integrateCQUAD prec n f a b = unsafePerformIO $ do
+    r <- malloc
+    e <- malloc
+    neval <- malloc
+    fp <- mkfun (\x _ -> f x)
+    c_integrate_cquad fp a b eps prec (fromIntegral n) r e neval // check "integrate_cquad"
+    vr <- peek r
+    ve <- peek e
+    vneval <- peek neval
+    let result = (vr,ve,vneval)
+    free r
+    free e
+    free neval
+    freeHaskellFunPtr fp
+    return result
+
+foreign import ccall safe "integrate_cquad" c_integrate_cquad
+    :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
+    -> Double -> Double -> CInt -> Ptr Double -> Ptr Double -> Ptr Int -> IO CInt
+
diff --git a/packages/gsl/src/Numeric/GSL/Internal.hs b/packages/gsl/src/Numeric/GSL/Internal.hs
new file mode 100644
index 0000000..a1c4e0c
--- /dev/null
+++ b/packages/gsl/src/Numeric/GSL/Internal.hs
@@ -0,0 +1,126 @@
+-- |
+-- Module      :  Numeric.GSL.Internal
+-- Copyright   :  (c) Alberto Ruiz 2009
+-- License     :  GPL
+-- Maintainer  :  Alberto Ruiz
+-- Stability   :  provisional
+--
+--
+-- Auxiliary functions.
+--
+
+
+module Numeric.GSL.Internal(
+    iv,
+    mkVecfun,
+    mkVecVecfun,
+    mkDoubleVecVecfun,
+    mkDoublefun,
+    aux_vTov,
+    mkVecMatfun,
+    mkDoubleVecMatfun,
+    aux_vTom,
+    createV,
+    createMIO,
+    module Data.Packed.Development,
+    check,
+    Res,TV,TM,TCV,TCM
+) where
+
+import Data.Packed
+import Data.Packed.Development hiding (check)
+import Data.Complex
+
+import Foreign.Marshal.Array(copyArray)
+import Foreign.Ptr(Ptr, FunPtr)
+import Foreign.C.Types
+import Foreign.C.String(peekCString)
+import System.IO.Unsafe(unsafePerformIO)
+import Data.Vector.Storable(unsafeWith)
+import Control.Monad(when)
+
+iv :: (Vector Double -> Double) -> (CInt -> Ptr Double -> Double)
+iv f n p = f (createV (fromIntegral n) copy "iv") where
+    copy n' q = do
+        copyArray q p (fromIntegral n')
+        return 0
+
+-- | conversion of Haskell functions into function pointers that can be used in the C side
+foreign import ccall safe "wrapper"
+    mkVecfun :: (CInt -> Ptr Double -> Double)
+             -> IO( FunPtr (CInt -> Ptr Double -> Double))
+
+foreign import ccall safe "wrapper"
+    mkVecVecfun :: TVV -> IO (FunPtr TVV)
+
+foreign import ccall safe "wrapper"
+    mkDoubleVecVecfun :: (Double -> TVV) -> IO (FunPtr (Double -> TVV))
+
+foreign import ccall safe "wrapper"
+    mkDoublefun :: (Double -> Double) -> IO (FunPtr (Double -> Double))
+
+aux_vTov :: (Vector Double -> Vector Double) -> TVV
+aux_vTov f n p nr r = g where
+    v = f x
+    x = createV (fromIntegral n) copy "aux_vTov"
+    copy n' q = do
+        copyArray q p (fromIntegral n')
+        return 0
+    g = do unsafeWith v $ \p' -> copyArray r p' (fromIntegral nr)
+           return 0
+
+foreign import ccall safe "wrapper"
+    mkVecMatfun :: TVM -> IO (FunPtr TVM)
+
+foreign import ccall safe "wrapper"
+    mkDoubleVecMatfun :: (Double -> TVM) -> IO (FunPtr (Double -> TVM))
+
+aux_vTom :: (Vector Double -> Matrix Double) -> TVM
+aux_vTom f n p rr cr r = g where
+    v = flatten $ f x
+    x = createV (fromIntegral n) copy "aux_vTov"
+    copy n' q = do
+        copyArray q p (fromIntegral n')
+        return 0
+    g = do unsafeWith v $ \p' -> copyArray r p' (fromIntegral $ rr*cr)
+           return 0
+
+createV n fun msg = unsafePerformIO $ do
+    r <- createVector n
+    app1 fun vec r msg
+    return r
+
+createMIO r c fun msg = do
+    res <- createMatrix RowMajor r c
+    app1 fun mat res msg
+    return res
+
+--------------------------------------------------------------------------------
+
+-- | check the error code
+check :: String -> IO CInt -> IO ()
+check msg f = do
+    err <- f
+    when (err/=0) $ do
+        ps <- gsl_strerror err
+        s <- peekCString ps
+        error (msg++": "++s)
+    return ()
+
+-- | description of GSL error codes
+foreign import ccall unsafe "gsl_strerror" gsl_strerror :: CInt -> IO (Ptr CChar)
+
+type PF = Ptr Float
+type PD = Ptr Double
+type PQ = Ptr (Complex Float)
+type PC = Ptr (Complex Double)
+
+type Res = IO CInt
+type TV x  = CInt -> PD -> x
+type TM x  = CInt -> CInt -> PD -> x
+type TCV x = CInt -> PC -> x
+type TCM x = CInt -> CInt -> PC -> x
+
+type TVV = TV (TV Res)
+type TVM = TV (TM Res)
+
diff --git a/packages/gsl/src/Numeric/GSL/LinearAlgebra.hs b/packages/gsl/src/Numeric/GSL/LinearAlgebra.hs
new file mode 100644
index 0000000..17e2258
--- /dev/null
+++ b/packages/gsl/src/Numeric/GSL/LinearAlgebra.hs
@@ -0,0 +1,135 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.GSL.LinearAlgebra
+-- Copyright   :  (c) Alberto Ruiz 2007-14
+-- License     :  GPL
+-- Maintainer  :  Alberto Ruiz
+-- Stability   :  provisional
+--
+-----------------------------------------------------------------------------
+
+module Numeric.GSL.LinearAlgebra (
+    RandDist(..), randomVector,
+    saveMatrix,
+    fwriteVector, freadVector, fprintfVector, fscanfVector,
+    fileDimensions, loadMatrix, fromFile
+) where
+
+import Data.Packed
+import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM)
+
+import Foreign.Marshal.Alloc(free)
+import Foreign.Ptr(Ptr)
+import Foreign.C.Types
+import Foreign.C.String(newCString)
+import System.IO.Unsafe(unsafePerformIO)
+import System.Process(readProcess)
+
+fromei x = fromIntegral (fromEnum x) :: CInt
+
+-----------------------------------------------------------------------
+
+data RandDist = Uniform  -- ^ uniform distribution in [0,1)
+              | Gaussian -- ^ normal distribution with mean zero and standard deviation one
+              deriving Enum
+
+-- | Obtains a vector of pseudorandom elements from the the mt19937 generator in GSL, with a given seed. Use randomIO to get a random seed.
+randomVector :: Int      -- ^ seed
+             -> RandDist -- ^ distribution
+             -> Int      -- ^ vector size
+             -> Vector Double
+randomVector seed dist n = unsafePerformIO $ do
+    r <- createVector n
+    app1 (c_random_vector (fi seed) ((fi.fromEnum) dist)) vec r "randomVector"
+    return r
+
+foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> TV
+
+--------------------------------------------------------------------------------
+
+-- | Saves a matrix as 2D ASCII table.
+saveMatrix :: FilePath
+           -> String     -- ^ format (%f, %g, %e)
+           -> Matrix Double
+           -> IO ()
+saveMatrix filename fmt m = do
+    charname <- newCString filename
+    charfmt <- newCString fmt
+    let o = if orderOf m == RowMajor then 1 else 0
+    app1 (matrix_fprintf charname charfmt o) mat m "matrix_fprintf"
+    free charname
+    free charfmt
+
+foreign import ccall unsafe "matrix_fprintf" matrix_fprintf :: Ptr CChar -> Ptr CChar -> CInt -> TM
+
+--------------------------------------------------------------------------------
+
+-- | Loads a vector from an ASCII file (the number of elements must be known in advance).
+fscanfVector :: FilePath -> Int -> IO (Vector Double)
+fscanfVector filename n = do
+    charname <- newCString filename
+    res <- createVector n
+    app1 (gsl_vector_fscanf charname) vec res "gsl_vector_fscanf"
+    free charname
+    return res
+
+foreign import ccall unsafe "vector_fscanf" gsl_vector_fscanf:: Ptr CChar -> TV
+
+-- | Saves the elements of a vector, with a given format (%f, %e, %g), to an ASCII file.
+fprintfVector :: FilePath -> String -> Vector Double -> IO ()
+fprintfVector filename fmt v = do
+    charname <- newCString filename
+    charfmt <- newCString fmt
+    app1 (gsl_vector_fprintf charname charfmt) vec v "gsl_vector_fprintf"
+    free charname
+    free charfmt
+
+foreign import ccall unsafe "vector_fprintf" gsl_vector_fprintf :: Ptr CChar -> Ptr CChar -> TV
+
+-- | Loads a vector from a binary file (the number of elements must be known in advance).
+freadVector :: FilePath -> Int -> IO (Vector Double)
+freadVector filename n = do
+    charname <- newCString filename
+    res <- createVector n
+    app1 (gsl_vector_fread charname) vec res "gsl_vector_fread"
+    free charname
+    return res
+
+foreign import ccall unsafe "vector_fread" gsl_vector_fread:: Ptr CChar -> TV
+
+-- | Saves the elements of a vector to a binary file.
+fwriteVector :: FilePath -> Vector Double -> IO ()
+fwriteVector filename v = do
+    charname <- newCString filename
+    app1 (gsl_vector_fwrite charname) vec v "gsl_vector_fwrite"
+    free charname
+
+foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV
+
+type PD = Ptr Double                            --
+type TV = CInt -> PD -> IO CInt                 --
+type TM = CInt -> CInt -> PD -> IO CInt         --
+
+--------------------------------------------------------------------------------
+
+{- |  obtains the number of rows and columns in an ASCII data file
+      (provisionally using unix's wc).
+-}
+fileDimensions :: FilePath -> IO (Int,Int)
+fileDimensions fname = do
+    wcres <- readProcess "wc" ["-w",fname] ""
+    contents <- readFile fname
+    let tot = read . head . words $ wcres
+        c   = length . head . dropWhile null . map words . lines $ contents
+    if tot > 0
+        then return (tot `div` c, c)
+        else return (0,0)
+
+-- | Loads a matrix from an ASCII file formatted as a 2D table.
+loadMatrix :: FilePath -> IO (Matrix Double)
+loadMatrix file = fromFile file =<< fileDimensions file
+
+-- | Loads a matrix from an ASCII file (the number of rows and columns must be known in advance).
+fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double)
+fromFile filename (r,c) = reshape c `fmap` fscanfVector filename (r*c)
+
diff --git a/packages/gsl/src/Numeric/GSL/Minimization.hs b/packages/gsl/src/Numeric/GSL/Minimization.hs
new file mode 100644
index 0000000..056d463
--- /dev/null
+++ b/packages/gsl/src/Numeric/GSL/Minimization.hs
@@ -0,0 +1,222 @@
+{- |
+Module      :  Numeric.GSL.Minimization
+Copyright   :  (c) Alberto Ruiz 2006-9
+License     :  GPL
+Maintainer  :  Alberto Ruiz
+Stability   :  provisional
+
+Minimization of a multidimensional function using some of the algorithms described in:
+
+<http://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Minimization.html>
+
+The example in the GSL manual:
+
+@
+f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30
+
+main = do
+    let (s,p) = minimize NMSimplex2 1E-2 30 [1,1] f [5,7]
+    print s
+    print p
+@
+
+>>> main
+[0.9920430849306288,1.9969168063253182]
+ 0.000  512.500  1.130  6.500  5.000
+ 1.000  290.625  1.409  5.250  4.000
+ 2.000  290.625  1.409  5.250  4.000
+ 3.000  252.500  1.409  5.500  1.000
+ ...
+22.000   30.001  0.013  0.992  1.997
+23.000   30.001  0.008  0.992  1.997
+
+The path to the solution can be graphically shown by means of:
+
+@'Graphics.Plot.mplot' $ drop 3 ('toColumns' p)@
+
+Taken from the GSL manual:
+
+The vector Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm is a quasi-Newton method which builds up an approximation to the second derivatives of the function f using the difference between successive gradient vectors. By combining the first and second derivatives the algorithm is able to take Newton-type steps towards the function minimum, assuming quadratic behavior in that region.
+
+The bfgs2 version of this minimizer is the most efficient version available, and is a faithful implementation of the line minimization scheme described in Fletcher's Practical Methods of Optimization, Algorithms 2.6.2 and 2.6.4. It supercedes the original bfgs routine and requires substantially fewer function and gradient evaluations. The user-supplied tolerance tol corresponds to the parameter \sigma used by Fletcher. A value of 0.1 is recommended for typical use (larger values correspond to less accurate line searches).
+
+The nmsimplex2 version is a new O(N) implementation of the earlier O(N^2) nmsimplex minimiser. It calculates the size of simplex as the rms distance of each vertex from the center rather than the mean distance, which has the advantage of allowing a linear update.
+
+-}
+
+
+module Numeric.GSL.Minimization (
+    minimize, minimizeV, MinimizeMethod(..),
+    minimizeD, minimizeVD, MinimizeMethodD(..),
+    uniMinimize, UniMinimizeMethod(..),
+
+    minimizeNMSimplex,
+    minimizeConjugateGradient,
+    minimizeVectorBFGS2
+) where
+
+
+import Data.Packed
+import Numeric.GSL.Internal
+
+import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr)
+import Foreign.C.Types
+import System.IO.Unsafe(unsafePerformIO)
+
+------------------------------------------------------------------------
+
+{-# DEPRECATED minimizeNMSimplex "use minimize NMSimplex2 eps maxit sizes f xi" #-}
+minimizeNMSimplex f xi szs eps maxit = minimize NMSimplex eps maxit szs f xi
+
+{-# DEPRECATED minimizeConjugateGradient "use minimizeD ConjugateFR eps maxit step tol f g xi" #-}
+minimizeConjugateGradient step tol eps maxit f g xi = minimizeD ConjugateFR eps maxit step tol f g xi
+
+{-# DEPRECATED minimizeVectorBFGS2 "use minimizeD VectorBFGS2 eps maxit step tol f g xi" #-}
+minimizeVectorBFGS2 step tol eps maxit f g xi = minimizeD VectorBFGS2 eps maxit step tol f g xi
+
+-------------------------------------------------------------------------
+
+data UniMinimizeMethod = GoldenSection
+                       | BrentMini
+                       | QuadGolden
+                       deriving (Enum, Eq, Show, Bounded)
+
+-- | Onedimensional minimization.
+
+uniMinimize :: UniMinimizeMethod -- ^ The method used.
+        -> Double               -- ^ desired precision of the solution
+        -> Int                  -- ^ maximum number of iterations allowed
+        -> (Double -> Double)   -- ^ function to minimize
+        -> Double               -- ^ guess for the location of the minimum
+        -> Double               -- ^ lower bound of search interval
+        -> Double               -- ^ upper bound of search interval
+        -> (Double, Matrix Double) -- ^ solution and optimization path
+
+uniMinimize method epsrel maxit fun xmin xl xu = uniMinimizeGen (fi (fromEnum method)) fun xmin xl xu epsrel maxit
+
+uniMinimizeGen m f xmin xl xu epsrel maxit = unsafePerformIO $ do
+    fp <- mkDoublefun f
+    rawpath <- createMIO maxit 4
+                         (c_uniMinize m fp epsrel (fi maxit) xmin xl xu)
+                         "uniMinimize"
+    let it = round (rawpath @@> (maxit-1,0))
+        path = takeRows it rawpath
+        [sol] = toLists $ dropRows (it-1) path
+    freeHaskellFunPtr fp
+    return (sol !! 1, path)
+
+
+foreign import ccall safe "uniMinimize"
+    c_uniMinize:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> Double -> TM Res
+
+data MinimizeMethod = NMSimplex
+                    | NMSimplex2
+                    deriving (Enum,Eq,Show,Bounded)
+
+-- | Minimization without derivatives
+minimize :: MinimizeMethod
+         -> Double              -- ^ desired precision of the solution (size test)
+         -> Int                 -- ^ maximum number of iterations allowed
+         -> [Double]            -- ^ sizes of the initial search box
+         -> ([Double] -> Double) -- ^ function to minimize
+         -> [Double]            -- ^ starting point
+         -> ([Double], Matrix Double) -- ^ solution vector and optimization path
+
+-- | Minimization without derivatives (vector version)
+minimizeV :: MinimizeMethod
+         -> Double              -- ^ desired precision of the solution (size test)
+         -> Int                 -- ^ maximum number of iterations allowed
+         -> Vector Double       -- ^ sizes of the initial search box
+         -> (Vector Double -> Double) -- ^ function to minimize
+         -> Vector Double            -- ^ starting point
+         -> (Vector Double, Matrix Double) -- ^ solution vector and optimization path
+
+minimize method eps maxit sz f xi = v2l $ minimizeV method eps maxit (fromList sz) (f.toList) (fromList xi)
+    where v2l (v,m) = (toList v, m)
+
+ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2
+
+minimizeV method eps maxit szv f xiv = unsafePerformIO $ do
+    let n   = dim xiv
+    fp <- mkVecfun (iv f)
+    rawpath <- ww2 vec xiv vec szv $ \xiv' szv' ->
+                   createMIO maxit (n+3)
+                         (c_minimize (fi (fromEnum method)) fp eps (fi maxit) // xiv' // szv')
+                         "minimize"
+    let it = round (rawpath @@> (maxit-1,0))
+        path = takeRows it rawpath
+        sol = flatten $ dropColumns 3 $ dropRows (it-1) path
+    freeHaskellFunPtr fp
+    return (sol, path)
+
+
+foreign import ccall safe "gsl-aux.h minimize"
+    c_minimize:: CInt -> FunPtr (CInt -> Ptr Double -> Double) -> Double -> CInt -> TV(TV(TM Res))
+
+----------------------------------------------------------------------------------
+
+
+data MinimizeMethodD = ConjugateFR
+                     | ConjugatePR
+                     | VectorBFGS
+                     | VectorBFGS2
+                     | SteepestDescent
+                     deriving (Enum,Eq,Show,Bounded)
+
+-- | Minimization with derivatives.
+minimizeD :: MinimizeMethodD
+    -> Double                 -- ^ desired precision of the solution (gradient test)
+    -> Int                    -- ^ maximum number of iterations allowed
+    -> Double                 -- ^ size of the first trial step
+    -> Double                 -- ^ tol (precise meaning depends on method)
+    -> ([Double] -> Double)   -- ^ function to minimize
+    -> ([Double] -> [Double]) -- ^ gradient
+    -> [Double]               -- ^ starting point
+    -> ([Double], Matrix Double) -- ^ solution vector and optimization path
+
+-- | Minimization with derivatives (vector version)
+minimizeVD :: MinimizeMethodD
+    -> Double                 -- ^ desired precision of the solution (gradient test)
+    -> Int                    -- ^ maximum number of iterations allowed
+    -> Double                 -- ^ size of the first trial step
+    -> Double                 -- ^ tol (precise meaning depends on method)
+    -> (Vector Double -> Double)   -- ^ function to minimize
+    -> (Vector Double -> Vector Double) -- ^ gradient
+    -> Vector Double               -- ^ starting point
+    -> (Vector Double, Matrix Double) -- ^ solution vector and optimization path
+
+minimizeD method eps maxit istep tol f df xi = v2l $ minimizeVD
+          method eps maxit istep tol (f.toList) (fromList.df.toList) (fromList xi)
+    where v2l (v,m) = (toList v, m)
+
+
+minimizeVD method eps maxit istep tol f df xiv = unsafePerformIO $ do
+    let n = dim xiv
+        f' = f
+        df' = (checkdim1 n . df)
+    fp <- mkVecfun (iv f')
+    dfp <- mkVecVecfun (aux_vTov df')
+    rawpath <- vec xiv $ \xiv' ->
+                    createMIO maxit (n+2)
+                         (c_minimizeD (fi (fromEnum method)) fp dfp istep tol eps (fi maxit) // xiv')
+                         "minimizeD"
+    let it = round (rawpath @@> (maxit-1,0))
+        path = takeRows it rawpath
+        sol = flatten $ dropColumns 2 $ dropRows (it-1) path
+    freeHaskellFunPtr fp
+    freeHaskellFunPtr dfp
+    return (sol,path)
+
+foreign import ccall safe "gsl-aux.h minimizeD"
+    c_minimizeD :: CInt
+                -> FunPtr (CInt -> Ptr Double -> Double)
+                -> FunPtr (TV (TV Res))
+                -> Double -> Double -> Double -> CInt
+                -> TV (TM Res)
+
+---------------------------------------------------------------------
+
+checkdim1 n v
+    | dim v == n = v
+    | otherwise = error $ "Error: "++ show n
+                        ++ " components expected in the result of the gradient supplied to minimizeD"
diff --git a/packages/gsl/src/Numeric/GSL/ODE.hs b/packages/gsl/src/Numeric/GSL/ODE.hs
new file mode 100644
index 0000000..7549a65
--- /dev/null
+++ b/packages/gsl/src/Numeric/GSL/ODE.hs
@@ -0,0 +1,140 @@
+{- |
+Module      :  Numeric.GSL.ODE
+Copyright   :  (c) Alberto Ruiz 2010
+License     :  GPL
+Maintainer  :  Alberto Ruiz
+Stability   :  provisional
+
+Solution of ordinary differential equation (ODE) initial value problems.
+
+<http://www.gnu.org/software/gsl/manual/html_node/Ordinary-Differential-Equations.html>
+
+A simple example:
+
+@
+import Numeric.GSL.ODE
+import Numeric.LinearAlgebra
+import Numeric.LinearAlgebra.Util(mplot)
+
+xdot t [x,v] = [v, -0.95*x - 0.1*v]
+
+ts = linspace 100 (0,20 :: Double)
+
+sol = odeSolve xdot [10,0] ts
+
+main = mplot (ts : toColumns sol)
+@
+
+-}
+-----------------------------------------------------------------------------
+
+module Numeric.GSL.ODE (
+    odeSolve, odeSolveV, ODEMethod(..), Jacobian
+) where
+
+import Data.Packed
+import Numeric.GSL.Internal
+
+import Foreign.Ptr(FunPtr, nullFunPtr, freeHaskellFunPtr)
+import Foreign.C.Types
+import System.IO.Unsafe(unsafePerformIO)
+
+-------------------------------------------------------------------------
+
+type TVV  = TV (TV Res)
+type TVM  = TV (TM Res)
+type TVVM = TV (TV (TM Res))
+
+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 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 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
+    :: (Double -> [Double] -> [Double])        -- ^ xdot(t,x)
+    -> [Double]        -- ^ initial conditions
+    -> Vector Double   -- ^ desired solution times
+    -> Matrix Double   -- ^ solution
+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
+          l2v f = \t -> fromList  . f t . toList
+
+-- | Evolution of the system with adaptive step-size control.
+odeSolveV
+    :: ODEMethod
+    -> 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 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 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
+        Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t))
+        Nothing  -> return nullFunPtr
+    sol <- vec xiv $ \xiv' ->
+            vec (checkTimes ts) $ \ts' ->
+             createMIO (dim ts) n
+              (ode_c (method) h epsAbs epsRel fp jp // xiv' // ts' )
+              "ode"
+    freeHaskellFunPtr fp
+    return sol
+
+foreign import ccall safe "ode"
+    ode_c :: CInt -> Double -> Double -> Double -> FunPtr (Double -> TVV) -> FunPtr (Double -> TVM) -> TVVM
+
+-------------------------------------------------------
+
+checkdim1 n v
+    | dim v == n = v
+    | otherwise = error $ "Error: "++ show n
+                        ++ " components expected in the result of the function supplied to odeSolve"
+
+checkdim2 n m
+    | rows m == n && cols m == n = m
+    | otherwise = error $ "Error: "++ show n ++ "x" ++ show n
+                        ++ " Jacobian expected in odeSolve"
+
+checkTimes ts | dim ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts
+              | otherwise = error "odeSolve requires increasing times"
+    where ts' = toList ts
diff --git a/packages/gsl/src/Numeric/GSL/Polynomials.hs b/packages/gsl/src/Numeric/GSL/Polynomials.hs
new file mode 100644
index 0000000..b1be85d
--- /dev/null
+++ b/packages/gsl/src/Numeric/GSL/Polynomials.hs
@@ -0,0 +1,57 @@
+{- |
+Module      :  Numeric.GSL.Polynomials
+Copyright   :  (c) Alberto Ruiz 2006
+License     :  GPL
+Maintainer  :  Alberto Ruiz
+Stability   :  provisional
+
+Polynomials.
+
+<http://www.gnu.org/software/gsl/manual/html_node/General-Polynomial-Equations.html#General-Polynomial-Equations>
+
+-}
+
+
+module Numeric.GSL.Polynomials (
+    polySolve
+) where
+
+import Data.Packed
+import Numeric.GSL.Internal
+import Data.Complex
+import System.IO.Unsafe (unsafePerformIO)
+
+#if __GLASGOW_HASKELL__ >= 704
+import Foreign.C.Types (CInt(..))
+#endif
+
+{- | Solution of general polynomial equations, using /gsl_poly_complex_solve/. 
+
+For example, the three solutions of x^3 + 8 = 0
+
+>>> polySolve [8,0,0,1]
+[(-2.0) :+ 0.0,1.0 :+ 1.7320508075688776,1.0 :+ (-1.7320508075688776)]
+
+
+The example in the GSL manual: To find the roots of x^5 -1 = 0:
+
+>>> polySolve [-1, 0, 0, 0, 0, 1]
+[(-0.8090169943749472) :+ 0.5877852522924731,
+(-0.8090169943749472) :+ (-0.5877852522924731),
+0.30901699437494756 :+ 0.9510565162951535,
+0.30901699437494756 :+ (-0.9510565162951535),
+1.0000000000000002 :+ 0.0]
+
+-}  
+polySolve :: [Double] -> [Complex Double]
+polySolve = toList . polySolve' . fromList
+
+polySolve' :: Vector Double -> Vector (Complex Double)
+polySolve' v | dim v > 1 = unsafePerformIO $ do
+    r <- createVector (dim v-1)
+    app2 c_polySolve vec v vec r "polySolve"
+    return r
+             | otherwise = error "polySolve on a polynomial of degree zero"
+
+foreign import ccall unsafe "gsl-aux.h polySolve" c_polySolve:: TV (TCV Res)
+
diff --git a/packages/gsl/src/Numeric/GSL/Random.hs b/packages/gsl/src/Numeric/GSL/Random.hs
new file mode 100644
index 0000000..2872b17
--- /dev/null
+++ b/packages/gsl/src/Numeric/GSL/Random.hs
@@ -0,0 +1,84 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.GSL.Random
+-- Copyright   :  (c) Alberto Ruiz 2009-14
+-- License     :  GPL
+--
+-- Maintainer  :  Alberto Ruiz
+-- Stability   :  provisional
+--
+-- Random vectors and matrices.
+--
+-----------------------------------------------------------------------------
+
+module Numeric.GSL.Random (
+    Seed,
+    RandDist(..),
+    randomVector,
+    gaussianSample,
+    uniformSample,
+    rand, randn
+) where
+
+import Numeric.GSL.Vector
+import Numeric.Container
+import Numeric.LinearAlgebra(Seed,RandDist(..),cholSH)
+import System.Random(randomIO)
+
+
+
+
+-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
+-- Gaussian distribution.
+gaussianSample :: Seed
+               -> Int -- ^ number of rows
+               -> Vector Double -- ^ mean vector
+               -> Matrix Double -- ^ covariance matrix
+               -> Matrix Double -- ^ result
+gaussianSample seed n med cov = m where
+    c = dim med
+    meds = konst' 1 n `outer` med
+    rs = reshape c $ randomVector seed Gaussian (c * n)
+    m = rs `mXm` cholSH cov `add` meds
+
+-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
+-- uniform distribution.
+uniformSample :: Seed
+               -> Int -- ^ number of rows
+               -> [(Double,Double)] -- ^ ranges for each column
+               -> Matrix Double -- ^ result
+uniformSample seed n rgs = m where
+    (as,bs) = unzip rgs
+    a = fromList as
+    cs = zipWith subtract as bs
+    d = dim a
+    dat = toRows $ reshape n $ randomVector seed Uniform (n*d)
+    am = konst' 1 n `outer` a
+    m = fromColumns (zipWith scale cs dat) `add` am
+
+-- | pseudorandom matrix with uniform elements between 0 and 1
+randm :: RandDist
+     -> Int -- ^ rows
+     -> Int -- ^ columns
+     -> IO (Matrix Double)
+randm d r c = do
+    seed <- randomIO
+    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
+
+>>> x <- randn 3 5
+>>> disp 3 x
+3x5
+0.386  -1.141   0.491  -0.510   1.512
+0.069  -0.919   1.022  -0.181   0.745
+0.313  -0.670  -0.097  -1.575  -0.583
+
+-}
+randn :: Int -> Int -> IO (Matrix Double)
+randn = randm Gaussian
+
diff --git a/packages/gsl/src/Numeric/GSL/Root.hs b/packages/gsl/src/Numeric/GSL/Root.hs
new file mode 100644
index 0000000..b9f3b94
--- /dev/null
+++ b/packages/gsl/src/Numeric/GSL/Root.hs
@@ -0,0 +1,199 @@
+{- |
+Module      :  Numeric.GSL.Root
+Copyright   :  (c) Alberto Ruiz 2009
+License     :  GPL
+Maintainer  :  Alberto Ruiz
+Stability   :  provisional
+
+Multidimensional root finding.
+
+<http://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Root_002dFinding.html>
+
+The example in the GSL manual:
+
+>>> let rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ]
+>>> let (sol,path) = root Hybrids 1E-7 30 (rosenbrock 1 10) [-10,-5]
+>>> sol
+[1.0,1.0]
+>>> disp 3 path
+11x5
+ 1.000  -10.000  -5.000  11.000  -1050.000
+ 2.000   -3.976  24.827   4.976     90.203
+ 3.000   -3.976  24.827   4.976     90.203
+ 4.000   -3.976  24.827   4.976     90.203
+ 5.000   -1.274  -5.680   2.274    -73.018
+ 6.000   -1.274  -5.680   2.274    -73.018
+ 7.000    0.249   0.298   0.751      2.359
+ 8.000    0.249   0.298   0.751      2.359
+ 9.000    1.000   0.878  -0.000     -1.218
+10.000    1.000   0.989  -0.000     -0.108
+11.000    1.000   1.000   0.000      0.000
+
+-}
+-----------------------------------------------------------------------------
+
+module Numeric.GSL.Root (
+    uniRoot, UniRootMethod(..),
+    uniRootJ, UniRootMethodJ(..),
+    root, RootMethod(..),
+    rootJ, RootMethodJ(..),
+) where
+
+import Data.Packed
+import Numeric.GSL.Internal
+import Foreign.Ptr(FunPtr, freeHaskellFunPtr)
+import Foreign.C.Types
+import System.IO.Unsafe(unsafePerformIO)
+
+-------------------------------------------------------------------------
+type TVV = TV (TV Res)
+type TVM = TV (TM Res)
+
+
+data UniRootMethod = Bisection
+                   | FalsePos
+                   | Brent
+                   deriving (Enum, Eq, Show, Bounded)
+
+uniRoot :: UniRootMethod
+        -> Double
+        -> Int
+        -> (Double -> Double)
+        -> Double
+        -> Double
+        -> (Double, Matrix Double)
+uniRoot method epsrel maxit fun xl xu = uniRootGen (fi (fromEnum method)) fun xl xu epsrel maxit
+
+uniRootGen m f xl xu epsrel maxit = unsafePerformIO $ do
+    fp <- mkDoublefun f
+    rawpath <- createMIO maxit 4
+                         (c_root m fp epsrel (fi maxit) xl xu)
+                         "root"
+    let it = round (rawpath @@> (maxit-1,0))
+        path = takeRows it rawpath
+        [sol] = toLists $ dropRows (it-1) path
+    freeHaskellFunPtr fp
+    return (sol !! 1, path)
+
+foreign import ccall safe "root"
+    c_root:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> TM Res
+
+-------------------------------------------------------------------------
+data UniRootMethodJ = UNewton
+                    | Secant
+                    | Steffenson
+                    deriving (Enum, Eq, Show, Bounded)
+
+uniRootJ :: UniRootMethodJ
+        -> Double
+        -> Int
+        -> (Double -> Double)
+        -> (Double -> Double)
+        -> Double
+        -> (Double, Matrix Double)
+uniRootJ method epsrel maxit fun dfun x = uniRootJGen (fi (fromEnum method)) fun
+    dfun x epsrel maxit
+
+uniRootJGen m f df x epsrel maxit = unsafePerformIO $ do
+    fp <- mkDoublefun f
+    dfp <- mkDoublefun df
+    rawpath <- createMIO maxit 2
+                         (c_rootj m fp dfp epsrel (fi maxit) x)
+                         "rootj"
+    let it = round (rawpath @@> (maxit-1,0))
+        path = takeRows it rawpath
+        [sol] = toLists $ dropRows (it-1) path
+    freeHaskellFunPtr fp
+    return (sol !! 1, path)
+
+foreign import ccall safe "rootj"
+    c_rootj :: CInt -> FunPtr (Double -> Double) -> FunPtr (Double -> Double)
+            -> Double -> CInt -> Double -> TM Res
+
+-------------------------------------------------------------------------
+
+data RootMethod = Hybrids
+                | Hybrid
+                | DNewton
+                | Broyden
+                deriving (Enum,Eq,Show,Bounded)
+
+-- | Nonlinear multidimensional root finding using algorithms that do not require 
+-- any derivative information to be supplied by the user.
+-- Any derivatives needed are approximated by finite differences.
+root :: RootMethod
+     -> Double                     -- ^ maximum residual
+     -> Int                        -- ^ maximum number of iterations allowed
+     -> ([Double] -> [Double])     -- ^ function to minimize
+     -> [Double]                   -- ^ starting point
+     -> ([Double], Matrix Double)  -- ^ solution vector and optimization path
+
+root method epsabs maxit fun xinit = rootGen (fi (fromEnum method)) fun xinit epsabs maxit
+
+rootGen m f xi epsabs maxit = unsafePerformIO $ do
+    let xiv = fromList xi
+        n   = dim xiv
+    fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
+    rawpath <- vec xiv $ \xiv' ->
+                   createMIO maxit (2*n+1)
+                         (c_multiroot m fp epsabs (fi maxit) // xiv')
+                         "multiroot"
+    let it = round (rawpath @@> (maxit-1,0))
+        path = takeRows it rawpath
+        [sol] = toLists $ dropRows (it-1) path
+    freeHaskellFunPtr fp
+    return (take n $ drop 1 sol, path)
+
+
+foreign import ccall safe "multiroot"
+    c_multiroot:: CInt -> FunPtr TVV -> Double -> CInt -> TVM
+
+-------------------------------------------------------------------------
+
+data RootMethodJ = HybridsJ
+                 | HybridJ
+                 | Newton
+                 | GNewton
+                deriving (Enum,Eq,Show,Bounded)
+
+-- | Nonlinear multidimensional root finding using both the function and its derivatives.
+rootJ :: RootMethodJ
+      -> Double                     -- ^ maximum residual
+      -> Int                        -- ^ maximum number of iterations allowed
+      -> ([Double] -> [Double])     -- ^ function to minimize
+      -> ([Double] -> [[Double]])   -- ^ Jacobian
+      -> [Double]                   -- ^ starting point
+      -> ([Double], Matrix Double)  -- ^ solution vector and optimization path
+
+rootJ method epsabs maxit fun jac xinit = rootJGen (fi (fromEnum method)) fun jac xinit epsabs maxit
+
+rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do
+    let xiv = fromList xi
+        n   = dim xiv
+    fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
+    jp <- mkVecMatfun (aux_vTom (checkdim2 n . fromLists . jac . toList))
+    rawpath <- vec xiv $ \xiv' ->
+                   createMIO maxit (2*n+1)
+                         (c_multirootj m fp jp epsabs (fi maxit) // xiv')
+                         "multiroot"
+    let it = round (rawpath @@> (maxit-1,0))
+        path = takeRows it rawpath
+        [sol] = toLists $ dropRows (it-1) path
+    freeHaskellFunPtr fp
+    freeHaskellFunPtr jp
+    return (take n $ drop 1 sol, path)
+
+foreign import ccall safe "multirootj"
+    c_multirootj:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> CInt -> TVM
+
+-------------------------------------------------------
+
+checkdim1 n v
+    | dim v == n = v
+    | otherwise = error $ "Error: "++ show n
+                        ++ " components expected in the result of the function supplied to root"
+
+checkdim2 n m
+    | rows m == n && cols m == n = m
+    | otherwise = error $ "Error: "++ show n ++ "x" ++ show n
+                        ++ " Jacobian expected in rootJ"
diff --git a/packages/gsl/src/Numeric/GSL/Vector.hs b/packages/gsl/src/Numeric/GSL/Vector.hs
new file mode 100644
index 0000000..af79f32
--- /dev/null
+++ b/packages/gsl/src/Numeric/GSL/Vector.hs
@@ -0,0 +1,107 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.GSL.Vector
+-- Copyright   :  (c) Alberto Ruiz 2007-14
+-- License     :  GPL
+-- Maintainer  :  Alberto Ruiz
+-- Stability   :  provisional
+--
+-----------------------------------------------------------------------------
+
+module Numeric.GSL.Vector (
+    randomVector,
+    saveMatrix,
+    fwriteVector, freadVector, fprintfVector, fscanfVector
+) where
+
+import Data.Packed
+import Numeric.LinearAlgebra(RandDist(..))
+import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM)
+
+import Foreign.Marshal.Alloc(free)
+import Foreign.Ptr(Ptr)
+import Foreign.C.Types
+import Foreign.C.String(newCString)
+import System.IO.Unsafe(unsafePerformIO)
+
+fromei x = fromIntegral (fromEnum x) :: CInt
+
+-----------------------------------------------------------------------
+
+-- | Obtains a vector of pseudorandom elements from the the mt19937 generator in GSL, with a given seed. Use randomIO to get a random seed.
+randomVector :: Int      -- ^ seed
+             -> RandDist -- ^ distribution
+             -> Int      -- ^ vector size
+             -> Vector Double
+randomVector seed dist n = unsafePerformIO $ do
+    r <- createVector n
+    app1 (c_random_vector_GSL (fi seed) ((fi.fromEnum) dist)) vec r "randomVectorGSL"
+    return r
+
+foreign import ccall unsafe "random_vector_GSL" c_random_vector_GSL :: CInt -> CInt -> TV
+
+--------------------------------------------------------------------------------
+
+-- | Saves a matrix as 2D ASCII table.
+saveMatrix :: FilePath
+           -> String     -- ^ format (%f, %g, %e)
+           -> Matrix Double
+           -> IO ()
+saveMatrix filename fmt m = do
+    charname <- newCString filename
+    charfmt <- newCString fmt
+    let o = if orderOf m == RowMajor then 1 else 0
+    app1 (matrix_fprintf charname charfmt o) mat m "matrix_fprintf"
+    free charname
+    free charfmt
+
+foreign import ccall unsafe "matrix_fprintf" matrix_fprintf :: Ptr CChar -> Ptr CChar -> CInt -> TM
+
+--------------------------------------------------------------------------------
+
+-- | Loads a vector from an ASCII file (the number of elements must be known in advance).
+fscanfVector :: FilePath -> Int -> IO (Vector Double)
+fscanfVector filename n = do
+    charname <- newCString filename
+    res <- createVector n
+    app1 (gsl_vector_fscanf charname) vec res "gsl_vector_fscanf"
+    free charname
+    return res
+
+foreign import ccall unsafe "vector_fscanf" gsl_vector_fscanf:: Ptr CChar -> TV
+
+-- | Saves the elements of a vector, with a given format (%f, %e, %g), to an ASCII file.
+fprintfVector :: FilePath -> String -> Vector Double -> IO ()
+fprintfVector filename fmt v = do
+    charname <- newCString filename
+    charfmt <- newCString fmt
+    app1 (gsl_vector_fprintf charname charfmt) vec v "gsl_vector_fprintf"
+    free charname
+    free charfmt
+
+foreign import ccall unsafe "vector_fprintf" gsl_vector_fprintf :: Ptr CChar -> Ptr CChar -> TV
+
+-- | Loads a vector from a binary file (the number of elements must be known in advance).
+freadVector :: FilePath -> Int -> IO (Vector Double)
+freadVector filename n = do
+    charname <- newCString filename
+    res <- createVector n
+    app1 (gsl_vector_fread charname) vec res "gsl_vector_fread"
+    free charname
+    return res
+
+foreign import ccall unsafe "vector_fread" gsl_vector_fread:: Ptr CChar -> TV
+
+-- | Saves the elements of a vector to a binary file.
+fwriteVector :: FilePath -> Vector Double -> IO ()
+fwriteVector filename v = do
+    charname <- newCString filename
+    app1 (gsl_vector_fwrite charname) vec v "gsl_vector_fwrite"
+    free charname
+
+foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV
+
+type PD = Ptr Double                            --
+type TV = CInt -> PD -> IO CInt                 --
+type TM = CInt -> CInt -> PD -> IO CInt         --
+
diff --git a/packages/gsl/src/Numeric/GSL/gsl-aux.c b/packages/gsl/src/Numeric/GSL/gsl-aux.c
new file mode 100644
index 0000000..ffc5c20
--- /dev/null
+++ b/packages/gsl/src/Numeric/GSL/gsl-aux.c
@@ -0,0 +1,945 @@
+#include <gsl/gsl_complex.h>
+
+#define RVEC(A) int A##n, double*A##p
+#define RMAT(A) int A##r, int A##c, double* A##p
+#define KRVEC(A) int A##n, const double*A##p
+#define KRMAT(A) int A##r, int A##c, const double* A##p
+
+#define CVEC(A) int A##n, gsl_complex*A##p
+#define CMAT(A) int A##r, int A##c, gsl_complex* A##p
+#define KCVEC(A) int A##n, const gsl_complex*A##p
+#define KCMAT(A) int A##r, int A##c, const gsl_complex* A##p
+
+#define FVEC(A) int A##n, float*A##p
+#define FMAT(A) int A##r, int A##c, float* A##p
+#define KFVEC(A) int A##n, const float*A##p
+#define KFMAT(A) int A##r, int A##c, const float* A##p
+
+#define QVEC(A) int A##n, gsl_complex_float*A##p
+#define QMAT(A) int A##r, int A##c, gsl_complex_float* A##p
+#define KQVEC(A) int A##n, const gsl_complex_float*A##p
+#define KQMAT(A) int A##r, int A##c, const gsl_complex_float* A##p
+
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_fft_complex.h>
+#include <gsl/gsl_integration.h>
+#include <gsl/gsl_deriv.h>
+#include <gsl/gsl_poly.h>
+#include <gsl/gsl_multimin.h>
+#include <gsl/gsl_multiroots.h>
+#include <gsl/gsl_min.h>
+#include <gsl/gsl_complex_math.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+#include <gsl/gsl_roots.h>
+#include <gsl/gsl_multifit_nlin.h>
+#include <string.h>
+#include <stdio.h>
+
+#define MACRO(B) do {B} while (0)
+#define ERROR(CODE) MACRO(return CODE;)
+#define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);})
+#define OK return 0;
+
+#define MIN(A,B) ((A)<(B)?(A):(B))
+#define MAX(A,B) ((A)>(B)?(A):(B))
+
+#ifdef DBG
+#define DEBUGMSG(M) printf("*** calling aux C function: %s\n",M);
+#else
+#define DEBUGMSG(M)
+#endif
+
+#define CHECK(RES,CODE) MACRO(if(RES) return CODE;)
+
+#ifdef DBG
+#define DEBUGMAT(MSG,X) printf(MSG" = \n"); gsl_matrix_fprintf(stdout,X,"%f"); printf("\n");
+#else
+#define DEBUGMAT(MSG,X)
+#endif
+
+#ifdef DBG
+#define DEBUGVEC(MSG,X) printf(MSG" = \n"); gsl_vector_fprintf(stdout,X,"%f"); printf("\n");
+#else
+#define DEBUGVEC(MSG,X)
+#endif
+
+#define DVVIEW(A) gsl_vector_view A = gsl_vector_view_array(A##p,A##n)
+#define DMVIEW(A) gsl_matrix_view A = gsl_matrix_view_array(A##p,A##r,A##c)
+#define CVVIEW(A) gsl_vector_complex_view A = gsl_vector_complex_view_array((double*)A##p,A##n)
+#define CMVIEW(A) gsl_matrix_complex_view A = gsl_matrix_complex_view_array((double*)A##p,A##r,A##c)
+#define KDVVIEW(A) gsl_vector_const_view A = gsl_vector_const_view_array(A##p,A##n)
+#define KDMVIEW(A) gsl_matrix_const_view A = gsl_matrix_const_view_array(A##p,A##r,A##c)
+#define KCVVIEW(A) gsl_vector_complex_const_view A = gsl_vector_complex_const_view_array((double*)A##p,A##n)
+#define KCMVIEW(A) gsl_matrix_complex_const_view A = gsl_matrix_complex_const_view_array((double*)A##p,A##r,A##c)
+
+#define FVVIEW(A) gsl_vector_float_view A = gsl_vector_float_view_array(A##p,A##n)
+#define FMVIEW(A) gsl_matrix_float_view A = gsl_matrix_float_view_array(A##p,A##r,A##c)
+#define QVVIEW(A) gsl_vector_complex_float_view A = gsl_vector_float_complex_view_array((float*)A##p,A##n)
+#define QMVIEW(A) gsl_matrix_complex_float_view A = gsl_matrix_float_complex_view_array((float*)A##p,A##r,A##c)
+#define KFVVIEW(A) gsl_vector_float_const_view A = gsl_vector_float_const_view_array(A##p,A##n)
+#define KFMVIEW(A) gsl_matrix_float_const_view A = gsl_matrix_float_const_view_array(A##p,A##r,A##c)
+#define KQVVIEW(A) gsl_vector_complex_float_const_view A = gsl_vector_complex_float_const_view_array((float*)A##p,A##n)
+#define KQMVIEW(A) gsl_matrix_complex_float_const_view A = gsl_matrix_complex_float_const_view_array((float*)A##p,A##r,A##c)
+
+#define V(a) (&a.vector)
+#define M(a) (&a.matrix)
+
+#define GCVEC(A) int A##n, gsl_complex*A##p
+#define KGCVEC(A) int A##n, const gsl_complex*A##p
+
+#define GQVEC(A) int A##n, gsl_complex_float*A##p
+#define KGQVEC(A) int A##n, const gsl_complex_float*A##p
+
+#define BAD_SIZE 2000
+#define BAD_CODE 2001
+#define MEM      2002
+#define BAD_FILE 2003
+
+
+void no_abort_on_error() {
+    gsl_set_error_handler_off();
+}
+
+
+
+int fft(int code, KCVEC(X), CVEC(R)) {
+    REQUIRES(Xn == Rn,BAD_SIZE);
+    DEBUGMSG("fft");
+    int s = Xn;
+    gsl_fft_complex_wavetable * wavetable = gsl_fft_complex_wavetable_alloc (s);
+    gsl_fft_complex_workspace * workspace = gsl_fft_complex_workspace_alloc (s);
+    gsl_vector_const_view X = gsl_vector_const_view_array((double*)Xp, 2*Xn);
+    gsl_vector_view R = gsl_vector_view_array((double*)Rp, 2*Rn);
+    gsl_blas_dcopy(&X.vector,&R.vector);
+    if(code==0) {
+        gsl_fft_complex_forward ((double*)Rp, 1, s, wavetable, workspace);
+    } else {
+        gsl_fft_complex_inverse ((double*)Rp, 1, s, wavetable, workspace);
+    }
+    gsl_fft_complex_wavetable_free (wavetable);
+    gsl_fft_complex_workspace_free (workspace);
+    OK
+}
+
+
+int deriv(int code, double f(double, void*), double x, double h, double * result, double * abserr)
+{
+    gsl_function F;
+    F.function = f;
+    F.params = 0;
+
+    if(code==0) return gsl_deriv_central (&F, x, h, result, abserr);
+
+    if(code==1) return gsl_deriv_forward (&F, x, h, result, abserr);
+
+    if(code==2) return gsl_deriv_backward (&F, x, h, result, abserr);
+
+    return 0;
+}
+
+
+int integrate_qng(double f(double, void*), double a, double b, double aprec, double prec,
+                   double *result, double*error) {
+    DEBUGMSG("integrate_qng");
+    gsl_function F;
+    F.function = f;
+    F.params = NULL;
+    size_t neval;
+    int res = gsl_integration_qng (&F, a,b, aprec, prec, result, error, &neval);
+    CHECK(res,res);
+    OK
+}
+
+int integrate_qags(double f(double,void*), double a, double b, double aprec, double prec, int w,
+               double *result, double* error) {
+    DEBUGMSG("integrate_qags");
+    gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
+    gsl_function F;
+    F.function = f;
+    F.params = NULL;
+    int res = gsl_integration_qags (&F, a,b, aprec, prec, w,wk, result, error);
+    CHECK(res,res);
+    gsl_integration_workspace_free (wk); 
+    OK
+}
+
+int integrate_qagi(double f(double,void*), double aprec, double prec, int w,
+               double *result, double* error) {
+    DEBUGMSG("integrate_qagi");
+    gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
+    gsl_function F;
+    F.function = f;
+    F.params = NULL;
+    int res = gsl_integration_qagi (&F, aprec, prec, w,wk, result, error);
+    CHECK(res,res);
+    gsl_integration_workspace_free (wk); 
+    OK
+}
+
+
+int integrate_qagiu(double f(double,void*), double a, double aprec, double prec, int w,
+               double *result, double* error) {
+    DEBUGMSG("integrate_qagiu");
+    gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
+    gsl_function F;
+    F.function = f;
+    F.params = NULL;
+    int res = gsl_integration_qagiu (&F, a, aprec, prec, w,wk, result, error);
+    CHECK(res,res);
+    gsl_integration_workspace_free (wk); 
+    OK
+}
+
+
+int integrate_qagil(double f(double,void*), double b, double aprec, double prec, int w,
+               double *result, double* error) {
+    DEBUGMSG("integrate_qagil");
+    gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
+    gsl_function F;
+    F.function = f;
+    F.params = NULL;
+    int res = gsl_integration_qagil (&F, b, aprec, prec, w,wk, result, error);
+    CHECK(res,res);
+    gsl_integration_workspace_free (wk); 
+    OK
+}
+
+int integrate_cquad(double f(double,void*), double a, double b, double aprec, double prec,
+                    int w, double *result, double* error, int *neval) {
+    DEBUGMSG("integrate_cquad");
+    gsl_integration_cquad_workspace * wk = gsl_integration_cquad_workspace_alloc (w);
+    gsl_function F;
+    F.function = f;
+    F.params = NULL;
+    size_t * sneval = NULL;
+    int res = gsl_integration_cquad (&F, a, b, aprec, prec, wk, result, error, sneval);
+    *neval = *sneval;
+    CHECK(res,res);
+    gsl_integration_cquad_workspace_free (wk); 
+    OK
+}
+
+
+int polySolve(KRVEC(a), CVEC(z)) {
+    DEBUGMSG("polySolve");
+    REQUIRES(an>1,BAD_SIZE);
+    gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc (an);
+    int res = gsl_poly_complex_solve ((double*)ap, an, w, (double*)zp);
+    CHECK(res,res);
+    gsl_poly_complex_workspace_free (w);
+    OK;
+}
+
+int vector_fscanf(char*filename, RVEC(a)) {
+    DEBUGMSG("gsl_vector_fscanf");
+    DVVIEW(a);
+    FILE * f = fopen(filename,"r");
+    CHECK(!f,BAD_FILE);
+    int res = gsl_vector_fscanf(f,V(a));
+    CHECK(res,res);
+    fclose (f);
+    OK
+}
+
+int vector_fprintf(char*filename, char*fmt, RVEC(a)) {
+    DEBUGMSG("gsl_vector_fprintf");
+    DVVIEW(a);
+    FILE * f = fopen(filename,"w");
+    CHECK(!f,BAD_FILE);
+    int res = gsl_vector_fprintf(f,V(a),fmt);
+    CHECK(res,res);
+    fclose (f);
+    OK
+}
+
+int vector_fread(char*filename, RVEC(a)) {
+    DEBUGMSG("gsl_vector_fread");
+    DVVIEW(a);
+    FILE * f = fopen(filename,"r");
+    CHECK(!f,BAD_FILE);
+    int res = gsl_vector_fread(f,V(a));
+    CHECK(res,res);
+    fclose (f);
+    OK
+}
+
+int vector_fwrite(char*filename, RVEC(a)) {
+    DEBUGMSG("gsl_vector_fwrite");
+    DVVIEW(a);
+    FILE * f = fopen(filename,"w");
+    CHECK(!f,BAD_FILE);
+    int res = gsl_vector_fwrite(f,V(a));
+    CHECK(res,res);
+    fclose (f);
+    OK
+}
+
+int matrix_fprintf(char*filename, char*fmt, int ro, RMAT(m)) {
+    DEBUGMSG("matrix_fprintf");
+    FILE * f = fopen(filename,"w");
+    CHECK(!f,BAD_FILE);
+    int i,j,sr,sc;
+    if (ro==1) { sr = mc; sc = 1;} else { sr = 1; sc = mr;}
+    #define AT(M,r,c) (M##p[(r)*sr+(c)*sc])
+    for (i=0; i<mr; i++) {
+        for (j=0; j<mc-1; j++) {
+            fprintf(f,fmt,AT(m,i,j));
+            fprintf(f," ");
+        }
+        fprintf(f,fmt,AT(m,i,j));
+        fprintf(f,"\n");
+    }
+    fclose (f);
+    OK
+}
+
+//---------------------------------------------------------------
+
+typedef double Trawfun(int, double*);
+
+double only_f_aux_min(const gsl_vector*x, void *pars) {
+    Trawfun * f = (Trawfun*) pars;  
+    double* p = (double*)calloc(x->size,sizeof(double));
+    int k;
+    for(k=0;k<x->size;k++) {
+        p[k] = gsl_vector_get(x,k);
+    }  
+    double res = f(x->size,p);
+    free(p);
+    return res;
+}
+
+double only_f_aux_root(double x, void *pars);
+int uniMinimize(int method, double f(double),
+                double epsrel, int maxit, double min,
+                double xl, double xu, RMAT(sol)) {
+   REQUIRES(solr == maxit && solc == 4,BAD_SIZE);
+   DEBUGMSG("minimize_only_f");
+   gsl_function my_func;
+   my_func.function = only_f_aux_root;
+   my_func.params = f;
+   size_t iter = 0;
+   int status;
+   const gsl_min_fminimizer_type *T;
+   gsl_min_fminimizer *s;
+   // Starting point
+   switch(method) {
+     case 0 : {T = gsl_min_fminimizer_goldensection; break; }
+     case 1 : {T = gsl_min_fminimizer_brent; break; }
+     case 2 : {T = gsl_min_fminimizer_quad_golden; break; }
+     default: ERROR(BAD_CODE);
+   }
+   s = gsl_min_fminimizer_alloc (T);
+   gsl_min_fminimizer_set (s, &my_func, min, xl, xu);
+   do {
+       double current_min, current_lo, current_hi;
+       status = gsl_min_fminimizer_iterate (s);
+       current_min = gsl_min_fminimizer_x_minimum (s);
+       current_lo = gsl_min_fminimizer_x_lower (s);
+       current_hi = gsl_min_fminimizer_x_upper (s);
+       solp[iter*solc] = iter + 1;
+       solp[iter*solc+1] = current_min;
+       solp[iter*solc+2] = current_lo;
+       solp[iter*solc+3] = current_hi;
+       iter++;
+       if (status)   /* check if solver is stuck */
+          break;
+       
+       status =
+         gsl_min_test_interval (current_lo, current_hi, 0, epsrel);
+   }
+   while (status == GSL_CONTINUE && iter < maxit);
+   int i;
+   for (i=iter; i<solr; i++) {
+       solp[i*solc+0] = iter;
+       solp[i*solc+1]=0.;
+       solp[i*solc+2]=0.;
+       solp[i*solc+3]=0.;
+   }
+   gsl_min_fminimizer_free(s);
+   OK
+}
+
+   
+
+// this version returns info about intermediate steps
+int minimize(int method, double f(int, double*), double tolsize, int maxit, 
+                 KRVEC(xi), KRVEC(sz), RMAT(sol)) {
+    REQUIRES(xin==szn && solr == maxit && solc == 3+xin,BAD_SIZE);
+    DEBUGMSG("minimizeList (nmsimplex)");
+    gsl_multimin_function my_func;
+    // extract function from pars
+    my_func.f = only_f_aux_min;
+    my_func.n = xin; 
+    my_func.params = f;
+    size_t iter = 0;
+    int status;
+    double size;
+    const gsl_multimin_fminimizer_type *T;
+    gsl_multimin_fminimizer *s = NULL;
+    // Initial vertex size vector 
+    KDVVIEW(sz);
+    // Starting point
+    KDVVIEW(xi);
+    // Minimizer nmsimplex, without derivatives
+    switch(method) {
+        case 0 : {T = gsl_multimin_fminimizer_nmsimplex; break; }
+#ifdef GSL110
+        case 1 : {T = gsl_multimin_fminimizer_nmsimplex; break; }
+#else
+        case 1 : {T = gsl_multimin_fminimizer_nmsimplex2; break; }
+#endif
+        default: ERROR(BAD_CODE);
+    }
+    s = gsl_multimin_fminimizer_alloc (T, my_func.n);
+    gsl_multimin_fminimizer_set (s, &my_func, V(xi), V(sz));
+    do {
+        status = gsl_multimin_fminimizer_iterate (s);
+        size = gsl_multimin_fminimizer_size (s);
+
+        solp[iter*solc+0] = iter+1;
+        solp[iter*solc+1] = s->fval;
+        solp[iter*solc+2] = size;
+
+        int k;
+        for(k=0;k<xin;k++) {
+            solp[iter*solc+k+3] = gsl_vector_get(s->x,k);
+        }
+        iter++;
+        if (status) break;
+        status = gsl_multimin_test_size (size, tolsize);
+    } while (status == GSL_CONTINUE && iter < maxit);
+    int i,j;
+    for (i=iter; i<solr; i++) {
+        solp[i*solc+0] = iter;
+        for(j=1;j<solc;j++) {
+            solp[i*solc+j]=0.;
+        }
+    }
+    gsl_multimin_fminimizer_free(s);
+    OK
+}
+
+// working with the gradient
+
+typedef struct {double (*f)(int, double*); int (*df)(int, double*, int, double*);} Tfdf;
+
+double f_aux_min(const gsl_vector*x, void *pars) {
+    Tfdf * fdf = ((Tfdf*) pars);
+    double* p = (double*)calloc(x->size,sizeof(double));
+    int k;
+    for(k=0;k<x->size;k++) {
+        p[k] = gsl_vector_get(x,k);
+    }
+    double res = fdf->f(x->size,p);
+    free(p);
+    return res;
+}
+
+
+void df_aux_min(const gsl_vector * x, void * pars, gsl_vector * g) {
+    Tfdf * fdf = ((Tfdf*) pars);
+    double* p = (double*)calloc(x->size,sizeof(double));
+    double* q = (double*)calloc(g->size,sizeof(double));
+    int k;
+    for(k=0;k<x->size;k++) {
+        p[k] = gsl_vector_get(x,k);
+    }
+
+    fdf->df(x->size,p,g->size,q);
+
+    for(k=0;k<x->size;k++) {
+        gsl_vector_set(g,k,q[k]);
+    }
+    free(p);
+    free(q);
+}
+
+void fdf_aux_min(const gsl_vector * x, void * pars, double * f, gsl_vector * g) {
+    *f = f_aux_min(x,pars);
+    df_aux_min(x,pars,g);
+}
+
+
+int minimizeD(int method, double f(int, double*), int df(int, double*, int, double*),
+              double initstep, double minimpar, double tolgrad, int maxit, 
+              KRVEC(xi), RMAT(sol)) {
+    REQUIRES(solr == maxit && solc == 2+xin,BAD_SIZE);
+    DEBUGMSG("minimizeWithDeriv (conjugate_fr)");
+    gsl_multimin_function_fdf my_func;
+    // extract function from pars
+    my_func.f = f_aux_min;
+    my_func.df = df_aux_min;
+    my_func.fdf = fdf_aux_min;
+    my_func.n = xin; 
+    Tfdf stfdf;
+    stfdf.f = f;
+    stfdf.df = df;
+    my_func.params = &stfdf;
+    size_t iter = 0;
+    int status;
+    const gsl_multimin_fdfminimizer_type *T;
+    gsl_multimin_fdfminimizer *s = NULL;
+    // Starting point
+    KDVVIEW(xi);
+    // conjugate gradient fr
+    switch(method) {
+        case 0 : {T = gsl_multimin_fdfminimizer_conjugate_fr; break; }
+        case 1 : {T = gsl_multimin_fdfminimizer_conjugate_pr; break; }
+        case 2 : {T = gsl_multimin_fdfminimizer_vector_bfgs; break; }
+        case 3 : {T = gsl_multimin_fdfminimizer_vector_bfgs2; break; }
+        case 4 : {T = gsl_multimin_fdfminimizer_steepest_descent; break; }
+        default: ERROR(BAD_CODE);
+    }
+    s = gsl_multimin_fdfminimizer_alloc (T, my_func.n);
+    gsl_multimin_fdfminimizer_set (s, &my_func, V(xi), initstep, minimpar);
+    do {
+        status = gsl_multimin_fdfminimizer_iterate (s);
+        solp[iter*solc+0] = iter+1;
+        solp[iter*solc+1] = s->f;
+        int k;
+        for(k=0;k<xin;k++) {
+            solp[iter*solc+k+2] = gsl_vector_get(s->x,k);
+        }
+        iter++;
+        if (status) break;
+        status = gsl_multimin_test_gradient (s->gradient, tolgrad);
+    } while (status == GSL_CONTINUE && iter < maxit);
+    int i,j;
+    for (i=iter; i<solr; i++) {
+        solp[i*solc+0] = iter;
+        for(j=1;j<solc;j++) {
+            solp[i*solc+j]=0.;
+        }
+    }
+    gsl_multimin_fdfminimizer_free(s);
+    OK
+}
+
+//---------------------------------------------------------------
+
+double only_f_aux_root(double x, void *pars) {
+    double (*f)(double) = (double (*)(double)) pars;
+    return f(x);
+}
+
+int root(int method, double f(double),
+         double epsrel, int maxit,
+         double xl, double xu, RMAT(sol)) {
+    REQUIRES(solr == maxit && solc == 4,BAD_SIZE);
+    DEBUGMSG("root_only_f");
+    gsl_function my_func;
+    // extract function from pars
+    my_func.function = only_f_aux_root;
+    my_func.params = f;
+    size_t iter = 0;
+    int status;
+    const gsl_root_fsolver_type *T;
+    gsl_root_fsolver *s;
+    // Starting point
+    switch(method) {
+        case 0 : {T = gsl_root_fsolver_bisection; printf("7\n"); break; }
+        case 1 : {T = gsl_root_fsolver_falsepos; break; }
+        case 2 : {T = gsl_root_fsolver_brent; break; }
+        default: ERROR(BAD_CODE);
+    }
+    s = gsl_root_fsolver_alloc (T);
+    gsl_root_fsolver_set (s, &my_func, xl, xu);
+    do {
+           double best, current_lo, current_hi;
+           status = gsl_root_fsolver_iterate (s);
+           best = gsl_root_fsolver_root (s);
+           current_lo = gsl_root_fsolver_x_lower (s);
+           current_hi = gsl_root_fsolver_x_upper (s);
+           solp[iter*solc] = iter + 1;
+           solp[iter*solc+1] = best;
+           solp[iter*solc+2] = current_lo;
+           solp[iter*solc+3] = current_hi;
+           iter++;
+           if (status)   /* check if solver is stuck */
+             break;
+
+           status =
+               gsl_root_test_interval (current_lo, current_hi, 0, epsrel);
+        }
+        while (status == GSL_CONTINUE && iter < maxit);
+    int i;
+    for (i=iter; i<solr; i++) {
+        solp[i*solc+0] = iter;
+        solp[i*solc+1]=0.;
+        solp[i*solc+2]=0.;
+        solp[i*solc+3]=0.;
+    }
+    gsl_root_fsolver_free(s);
+    OK
+}
+
+typedef struct {
+    double (*f)(double);
+    double (*jf)(double);
+} uniTfjf;
+
+double f_aux_uni(double x, void *pars) {
+    uniTfjf * fjf = ((uniTfjf*) pars);
+    return (fjf->f)(x);
+}
+
+double jf_aux_uni(double x, void * pars) {
+    uniTfjf * fjf = ((uniTfjf*) pars);
+    return (fjf->jf)(x);
+}
+
+void fjf_aux_uni(double x, void * pars, double * f, double * g) {
+    *f = f_aux_uni(x,pars);
+    *g = jf_aux_uni(x,pars);
+}
+
+int rootj(int method, double f(double),
+          double df(double),
+         double epsrel, int maxit,
+         double x, RMAT(sol)) {
+    REQUIRES(solr == maxit && solc == 2,BAD_SIZE);
+    DEBUGMSG("root_fjf");
+    gsl_function_fdf my_func;
+    // extract function from pars
+    my_func.f = f_aux_uni;
+    my_func.df = jf_aux_uni;
+    my_func.fdf = fjf_aux_uni;
+    uniTfjf stfjf;
+    stfjf.f = f;
+    stfjf.jf = df;
+    my_func.params = &stfjf;
+    size_t iter = 0;
+    int status;
+    const gsl_root_fdfsolver_type *T;
+    gsl_root_fdfsolver *s;
+    // Starting point
+    switch(method) {
+        case 0 : {T = gsl_root_fdfsolver_newton;; break; }
+        case 1 : {T = gsl_root_fdfsolver_secant; break; }
+        case 2 : {T = gsl_root_fdfsolver_steffenson; break; }
+        default: ERROR(BAD_CODE);
+    }
+    s = gsl_root_fdfsolver_alloc (T);
+
+    gsl_root_fdfsolver_set (s, &my_func, x);
+
+    do {
+           double x0;
+           status = gsl_root_fdfsolver_iterate (s);
+           x0 = x;
+           x = gsl_root_fdfsolver_root(s);
+           solp[iter*solc+0] = iter+1;
+           solp[iter*solc+1] = x;
+
+           iter++;
+           if (status)   /* check if solver is stuck */
+             break;
+
+           status =
+               gsl_root_test_delta (x, x0, 0, epsrel);
+        }
+        while (status == GSL_CONTINUE && iter < maxit);
+
+    int i;
+    for (i=iter; i<solr; i++) {
+        solp[i*solc+0] = iter;
+        solp[i*solc+1]=0.;
+    }
+    gsl_root_fdfsolver_free(s);
+    OK
+}
+
+
+//---------------------------------------------------------------
+
+typedef void TrawfunV(int, double*, int, double*);
+
+int only_f_aux_multiroot(const gsl_vector*x, void *pars, gsl_vector*y) {
+    TrawfunV * f = (TrawfunV*) pars;
+    double* p = (double*)calloc(x->size,sizeof(double));
+    double* q = (double*)calloc(y->size,sizeof(double));
+    int k;
+    for(k=0;k<x->size;k++) {
+        p[k] = gsl_vector_get(x,k);
+    }
+    f(x->size,p,y->size,q);
+    for(k=0;k<y->size;k++) {
+        gsl_vector_set(y,k,q[k]);
+    }
+    free(p);
+    free(q);
+    return 0; //hmmm
+}
+
+int multiroot(int method, void f(int, double*, int, double*),
+         double epsabs, int maxit,
+         KRVEC(xi), RMAT(sol)) {
+    REQUIRES(solr == maxit && solc == 1+2*xin,BAD_SIZE);
+    DEBUGMSG("root_only_f");
+    gsl_multiroot_function my_func;
+    // extract function from pars
+    my_func.f = only_f_aux_multiroot;
+    my_func.n = xin;
+    my_func.params = f;
+    size_t iter = 0;
+    int status;
+    const gsl_multiroot_fsolver_type *T;
+    gsl_multiroot_fsolver *s;
+    // Starting point
+    KDVVIEW(xi);
+    switch(method) {
+        case 0 : {T = gsl_multiroot_fsolver_hybrids;; break; }
+        case 1 : {T = gsl_multiroot_fsolver_hybrid; break; }
+        case 2 : {T = gsl_multiroot_fsolver_dnewton; break; }
+        case 3 : {T = gsl_multiroot_fsolver_broyden; break; }
+        default: ERROR(BAD_CODE);
+    }
+    s = gsl_multiroot_fsolver_alloc (T, my_func.n);
+    gsl_multiroot_fsolver_set (s, &my_func, V(xi));
+
+    do {
+           status = gsl_multiroot_fsolver_iterate (s);
+
+           solp[iter*solc+0] = iter+1;
+
+           int k;
+           for(k=0;k<xin;k++) {
+               solp[iter*solc+k+1] = gsl_vector_get(s->x,k);
+           }
+           for(k=xin;k<2*xin;k++) {
+               solp[iter*solc+k+1] = gsl_vector_get(s->f,k-xin);
+           }
+
+           iter++;
+           if (status)   /* check if solver is stuck */
+             break;
+
+           status =
+             gsl_multiroot_test_residual (s->f, epsabs);
+        }
+        while (status == GSL_CONTINUE && iter < maxit);
+
+    int i,j;
+    for (i=iter; i<solr; i++) {
+        solp[i*solc+0] = iter;
+        for(j=1;j<solc;j++) {
+            solp[i*solc+j]=0.;
+        }
+    }
+    gsl_multiroot_fsolver_free(s);
+    OK
+}
+
+// working with the jacobian
+
+typedef struct {int (*f)(int, double*, int, double *);
+                int (*jf)(int, double*, int, int, double*);} Tfjf;
+
+int f_aux(const gsl_vector*x, void *pars, gsl_vector*y) {
+    Tfjf * fjf = ((Tfjf*) pars);
+    double* p = (double*)calloc(x->size,sizeof(double));
+    double* q = (double*)calloc(y->size,sizeof(double));
+    int k;
+    for(k=0;k<x->size;k++) {
+        p[k] = gsl_vector_get(x,k);
+    }
+    (fjf->f)(x->size,p,y->size,q);
+    for(k=0;k<y->size;k++) {
+        gsl_vector_set(y,k,q[k]);
+    }
+    free(p);
+    free(q);
+    return 0;
+}
+
+int jf_aux(const gsl_vector * x, void * pars, gsl_matrix * jac) {
+    Tfjf * fjf = ((Tfjf*) pars);
+    double* p = (double*)calloc(x->size,sizeof(double));
+    double* q = (double*)calloc((jac->size1)*(jac->size2),sizeof(double));
+    int i,j,k;
+    for(k=0;k<x->size;k++) {
+        p[k] = gsl_vector_get(x,k);
+    }
+
+    (fjf->jf)(x->size,p,jac->size1,jac->size2,q);
+
+    k=0;
+    for(i=0;i<jac->size1;i++) {
+        for(j=0;j<jac->size2;j++){
+            gsl_matrix_set(jac,i,j,q[k++]);
+        }
+    }
+    free(p);
+    free(q);
+    return 0;
+}
+
+int fjf_aux(const gsl_vector * x, void * pars, gsl_vector * f, gsl_matrix * g) {
+    f_aux(x,pars,f);
+    jf_aux(x,pars,g);
+    return 0;
+}
+
+int multirootj(int method, int f(int, double*, int, double*),
+                      int jac(int, double*, int, int, double*),
+         double epsabs, int maxit,
+         KRVEC(xi), RMAT(sol)) {
+    REQUIRES(solr == maxit && solc == 1+2*xin,BAD_SIZE);
+    DEBUGMSG("root_fjf");
+    gsl_multiroot_function_fdf my_func;
+    // extract function from pars
+    my_func.f = f_aux;
+    my_func.df = jf_aux;
+    my_func.fdf = fjf_aux;
+    my_func.n = xin;
+    Tfjf stfjf;
+    stfjf.f = f;
+    stfjf.jf = jac;
+    my_func.params = &stfjf;
+    size_t iter = 0;
+    int status;
+    const gsl_multiroot_fdfsolver_type *T;
+    gsl_multiroot_fdfsolver *s;
+    // Starting point
+    KDVVIEW(xi);
+    switch(method) {
+        case 0 : {T = gsl_multiroot_fdfsolver_hybridsj;; break; }
+        case 1 : {T = gsl_multiroot_fdfsolver_hybridj; break; }
+        case 2 : {T = gsl_multiroot_fdfsolver_newton; break; }
+        case 3 : {T = gsl_multiroot_fdfsolver_gnewton; break; }
+        default: ERROR(BAD_CODE);
+    }
+    s = gsl_multiroot_fdfsolver_alloc (T, my_func.n);
+
+    gsl_multiroot_fdfsolver_set (s, &my_func, V(xi));
+
+    do {
+           status = gsl_multiroot_fdfsolver_iterate (s);
+
+           solp[iter*solc+0] = iter+1;
+
+           int k;
+           for(k=0;k<xin;k++) {
+               solp[iter*solc+k+1] = gsl_vector_get(s->x,k);
+           }
+           for(k=xin;k<2*xin;k++) {
+               solp[iter*solc+k+1] = gsl_vector_get(s->f,k-xin);
+           }
+
+           iter++;
+           if (status)   /* check if solver is stuck */
+             break;
+
+           status =
+             gsl_multiroot_test_residual (s->f, epsabs);
+        }
+        while (status == GSL_CONTINUE && iter < maxit);
+
+    int i,j;
+    for (i=iter; i<solr; i++) {
+        solp[i*solc+0] = iter;
+        for(j=1;j<solc;j++) {
+            solp[i*solc+j]=0.;
+        }
+    }
+    gsl_multiroot_fdfsolver_free(s);
+    OK
+}
+
+//-------------- non linear least squares fitting -------------------
+
+int nlfit(int method, int f(int, double*, int, double*),
+                      int jac(int, double*, int, int, double*),
+         double epsabs, double epsrel, int maxit, int p,
+         KRVEC(xi), RMAT(sol)) {
+    REQUIRES(solr == maxit && solc == 2+xin,BAD_SIZE);
+    DEBUGMSG("nlfit");
+    const gsl_multifit_fdfsolver_type *T;
+    gsl_multifit_fdfsolver *s;
+    gsl_multifit_function_fdf my_f;
+    // extract function from pars
+    my_f.f = f_aux;
+    my_f.df = jf_aux;
+    my_f.fdf = fjf_aux;
+    my_f.n = p;
+    my_f.p = xin;  // !!!!
+    Tfjf stfjf;
+    stfjf.f = f;
+    stfjf.jf = jac;
+    my_f.params = &stfjf;
+    size_t iter = 0;
+    int status;
+
+    KDVVIEW(xi);
+    //DMVIEW(cov);
+
+    switch(method) {
+        case 0 : { T = gsl_multifit_fdfsolver_lmsder; break; }
+        case 1 : { T = gsl_multifit_fdfsolver_lmder; break; }
+        default: ERROR(BAD_CODE);
+    }
+
+    s = gsl_multifit_fdfsolver_alloc (T, my_f.n, my_f.p);
+    gsl_multifit_fdfsolver_set (s, &my_f, V(xi));
+
+    do {   status = gsl_multifit_fdfsolver_iterate (s);
+
+           solp[iter*solc+0] = iter+1;
+           solp[iter*solc+1] = gsl_blas_dnrm2 (s->f);
+
+           int k;
+           for(k=0;k<xin;k++) {
+               solp[iter*solc+k+2] = gsl_vector_get(s->x,k);
+           }
+
+           iter++;
+           if (status)   /* check if solver is stuck */
+             break;
+
+           status = gsl_multifit_test_delta (s->dx, s->x, epsabs, epsrel);
+        }
+        while (status == GSL_CONTINUE && iter < maxit);
+
+    int i,j;
+    for (i=iter; i<solr; i++) {
+        solp[i*solc+0] = iter;
+        for(j=1;j<solc;j++) {
+            solp[i*solc+j]=0.;
+        }
+    }
+
+    //gsl_multifit_covar (s->J, 0.0, M(cov));
+
+    gsl_multifit_fdfsolver_free (s);
+    OK
+}
+
+
+//////////////////////////////////////////////////////
+
+
+#define RAN(C,F) case C: { for(k=0;k<rn;k++) { rp[k]= F(gen); }; OK }
+
+int random_vector_GSL(int seed, int code, RVEC(r)) {
+    DEBUGMSG("random_vector_GSL")
+    static gsl_rng * gen = NULL;
+    if (!gen) { gen = gsl_rng_alloc (gsl_rng_mt19937);}
+    gsl_rng_set (gen, seed);
+    int k;
+    switch (code) {
+        RAN(0,gsl_rng_uniform)
+        RAN(1,gsl_ran_ugaussian)
+        default: ERROR(BAD_CODE);
+    }
+}
+#undef RAN
+
+//////////////////////////////////////////////////////
+
+#include "gsl-ode.c"
+
+//////////////////////////////////////////////////////
diff --git a/packages/gsl/src/Numeric/GSL/gsl-ode.c b/packages/gsl/src/Numeric/GSL/gsl-ode.c
new file mode 100644
index 0000000..3f2771b
--- /dev/null
+++ b/packages/gsl/src/Numeric/GSL/gsl-ode.c
@@ -0,0 +1,182 @@
+
+#ifdef GSLODE1
+
+////////////////////////////// ODE V1 //////////////////////////////////////////
+
+#include <gsl/gsl_odeiv.h>
+
+typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode;
+
+int odefunc (double t, const double y[], double f[], void *params) { 
+    Tode * P = (Tode*) params;
+    (P->f)(t,P->n,y,P->n,f);
+    return GSL_SUCCESS;
+}
+
+int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) {
+     Tode * P = ((Tode*) params);
+     (P->j)(t,P->n,y,P->n,P->n,dfdy);
+     int j;
+     for (j=0; j< P->n; j++)
+        dfdt[j] = 0.0;
+     return GSL_SUCCESS;
+}
+
+
+int ode(int method, double h, double eps_abs, double eps_rel,
+        int f(double, int, const double*, int, double*),
+        int jac(double, int, const double*, int, int, double*),
+        KRVEC(xi), KRVEC(ts), RMAT(sol)) {
+
+    const gsl_odeiv_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 : { printf("Sorry: ODE rk1imp not available in this GSL version\n"); exit(0); }
+        case 9 : { printf("Sorry: ODE msadams not available in this GSL version\n"); exit(0); }
+        case 10: { printf("Sorry: ODE msbdf not available in this GSL version\n"); exit(0); }
+        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};
+
+    double t = tsp[0];
+
+    double* y = (double*)calloc(xin,sizeof(double));
+    int i,j;
+    for(i=0; i< xin; i++) {
+        y[i] = xip[i];
+        solp[i] = xip[i];
+    }
+
+       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;
+             }
+           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;
+}
+
+#else
+
+///////////////////// ODE V2 ///////////////////////////////////////////////////
+                   
+#include <gsl/gsl_odeiv2.h>
+
+typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode;
+
+int odefunc (double t, const double y[], double f[], void *params) { 
+    Tode * P = (Tode*) params;
+    (P->f)(t,P->n,y,P->n,f);
+    return GSL_SUCCESS;
+}
+
+int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) {
+     Tode * P = ((Tode*) params);
+     (P->j)(t,P->n,y,P->n,P->n,dfdy);
+     int j;
+     for (j=0; j< P->n; j++)
+        dfdt[j] = 0.0;
+     return GSL_SUCCESS;
+}
+
+
+int ode(int method, double h, double eps_abs, double eps_rel,
+        int f(double, int, const double*, int, double*),
+        int jac(double, int, const double*, int, int, double*),
+        KRVEC(xi), KRVEC(ts), RMAT(sol)) {
+
+    const gsl_odeiv2_step_type * T;
+
+    switch(method) {
+        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);
+    }
+
+    Tode P;
+    P.f = f;
+    P.j = jac;
+    P.n = xin;
+
+    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=0;
+    for(i=0; i< xin; i++) {
+        y[i] = xip[i];
+        solp[i] = xip[i];
+    }
+
+       for (i = 1; i < tsn ; i++)
+         {
+           double ti = tsp[i];
+           
+           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_odeiv2_driver_free (d);
+    
+    return status;
+}
+
+#endif
+