From 1925c123d7d8184a1d2ddc0a413e0fd2776e1083 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Thu, 8 May 2014 08:48:12 +0200 Subject: empty hmatrix-base --- .../hmatrix/src/Numeric/GSL/Differentiation.hs | 87 ++ packages/hmatrix/src/Numeric/GSL/Fitting.hs | 179 +++ packages/hmatrix/src/Numeric/GSL/Fourier.hs | 47 + packages/hmatrix/src/Numeric/GSL/Integration.hs | 250 ++++ packages/hmatrix/src/Numeric/GSL/Internal.hs | 76 + packages/hmatrix/src/Numeric/GSL/Minimization.hs | 227 +++ packages/hmatrix/src/Numeric/GSL/ODE.hs | 138 ++ packages/hmatrix/src/Numeric/GSL/Polynomials.hs | 58 + packages/hmatrix/src/Numeric/GSL/Root.hs | 199 +++ packages/hmatrix/src/Numeric/GSL/Vector.hs | 328 +++++ packages/hmatrix/src/Numeric/GSL/gsl-aux.c | 1541 ++++++++++++++++++++ packages/hmatrix/src/Numeric/GSL/gsl-ode.c | 182 +++ 12 files changed, 3312 insertions(+) create mode 100644 packages/hmatrix/src/Numeric/GSL/Differentiation.hs create mode 100644 packages/hmatrix/src/Numeric/GSL/Fitting.hs create mode 100644 packages/hmatrix/src/Numeric/GSL/Fourier.hs create mode 100644 packages/hmatrix/src/Numeric/GSL/Integration.hs create mode 100644 packages/hmatrix/src/Numeric/GSL/Internal.hs create mode 100644 packages/hmatrix/src/Numeric/GSL/Minimization.hs create mode 100644 packages/hmatrix/src/Numeric/GSL/ODE.hs create mode 100644 packages/hmatrix/src/Numeric/GSL/Polynomials.hs create mode 100644 packages/hmatrix/src/Numeric/GSL/Root.hs create mode 100644 packages/hmatrix/src/Numeric/GSL/Vector.hs create mode 100644 packages/hmatrix/src/Numeric/GSL/gsl-aux.c create mode 100644 packages/hmatrix/src/Numeric/GSL/gsl-ode.c (limited to 'packages/hmatrix/src/Numeric/GSL') diff --git a/packages/hmatrix/src/Numeric/GSL/Differentiation.hs b/packages/hmatrix/src/Numeric/GSL/Differentiation.hs new file mode 100644 index 0000000..93c5007 --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/Differentiation.hs @@ -0,0 +1,87 @@ +{-# OPTIONS #-} +----------------------------------------------------------------------------- +{- | +Module : Numeric.GSL.Differentiation +Copyright : (c) Alberto Ruiz 2006 +License : GPL-style + +Maintainer : Alberto Ruiz (aruiz at um dot es) +Stability : provisional +Portability : uses ffi + +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 Data.Packed.Internal(check,(//)) +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/hmatrix/src/Numeric/GSL/Fitting.hs b/packages/hmatrix/src/Numeric/GSL/Fitting.hs new file mode 100644 index 0000000..c4f3a91 --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/Fitting.hs @@ -0,0 +1,179 @@ +{- | +Module : Numeric.GSL.Fitting +Copyright : (c) Alberto Ruiz 2010 +License : GPL + +Maintainer : Alberto Ruiz (aruiz at um dot es) +Stability : provisional +Portability : uses ffi + +Nonlinear Least-Squares Fitting + + + +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 Data.Packed.Internal +import Numeric.LinearAlgebra +import Numeric.GSL.Internal + +import Foreign.Ptr(FunPtr, freeHaskellFunPtr) +import Foreign.C.Types +import System.IO.Unsafe(unsafePerformIO) + +------------------------------------------------------------------------- + +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/hmatrix/src/Numeric/GSL/Fourier.hs b/packages/hmatrix/src/Numeric/GSL/Fourier.hs new file mode 100644 index 0000000..86aedd6 --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/Fourier.hs @@ -0,0 +1,47 @@ +{-# LANGUAGE ForeignFunctionInterface #-} +----------------------------------------------------------------------------- +{- | +Module : Numeric.GSL.Fourier +Copyright : (c) Alberto Ruiz 2006 +License : GPL-style + +Maintainer : Alberto Ruiz (aruiz at um dot es) +Stability : provisional +Portability : uses ffi + +Fourier Transform. + + + +-} +----------------------------------------------------------------------------- +module Numeric.GSL.Fourier ( + fft, + ifft +) where + +import Data.Packed.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 -> TCVCV + + +{- | 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/hmatrix/src/Numeric/GSL/Integration.hs b/packages/hmatrix/src/Numeric/GSL/Integration.hs new file mode 100644 index 0000000..5f0a415 --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/Integration.hs @@ -0,0 +1,250 @@ +{-# OPTIONS #-} +----------------------------------------------------------------------------- +{- | +Module : Numeric.GSL.Integration +Copyright : (c) Alberto Ruiz 2006 +License : GPL-style + +Maintainer : Alberto Ruiz (aruiz at um dot es) +Stability : provisional +Portability : uses ffi + +Numerical integration routines. + + +-} +----------------------------------------------------------------------------- + +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 Data.Packed.Internal(check,(//)) +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/hmatrix/src/Numeric/GSL/Internal.hs b/packages/hmatrix/src/Numeric/GSL/Internal.hs new file mode 100644 index 0000000..69a9750 --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/Internal.hs @@ -0,0 +1,76 @@ +-- Module : Numeric.GSL.Internal +-- Copyright : (c) Alberto Ruiz 2009 +-- License : GPL +-- +-- Maintainer : Alberto Ruiz (aruiz at um dot es) +-- Stability : provisional +-- Portability : uses ffi +-- +-- Auxiliary functions. +-- +-- #hide + +module Numeric.GSL.Internal where + +import Data.Packed.Internal + +import Foreign.Marshal.Array(copyArray) +import Foreign.Ptr(Ptr, FunPtr) +import Foreign.C.Types +import System.IO.Unsafe(unsafePerformIO) + +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 diff --git a/packages/hmatrix/src/Numeric/GSL/Minimization.hs b/packages/hmatrix/src/Numeric/GSL/Minimization.hs new file mode 100644 index 0000000..1879dab --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/Minimization.hs @@ -0,0 +1,227 @@ +{-# LANGUAGE ForeignFunctionInterface #-} +----------------------------------------------------------------------------- +{- | +Module : Numeric.GSL.Minimization +Copyright : (c) Alberto Ruiz 2006-9 +License : GPL-style + +Maintainer : Alberto Ruiz (aruiz at um dot es) +Stability : provisional +Portability : uses ffi + +Minimization of a multidimensional function using some of the algorithms described in: + + + +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.Internal +import Data.Packed.Matrix +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 + +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 = cdat $ 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 -> TVVM + +---------------------------------------------------------------------------------- + + +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 = cdat $ 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 TVV + -> Double -> Double -> Double -> CInt + -> TVM + +--------------------------------------------------------------------- + +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/hmatrix/src/Numeric/GSL/ODE.hs b/packages/hmatrix/src/Numeric/GSL/ODE.hs new file mode 100644 index 0000000..9a29085 --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/ODE.hs @@ -0,0 +1,138 @@ +{- | +Module : Numeric.GSL.ODE +Copyright : (c) Alberto Ruiz 2010 +License : GPL + +Maintainer : Alberto Ruiz (aruiz at um dot es) +Stability : provisional +Portability : uses ffi + +Solution of ordinary differential equation (ODE) initial value problems. + + + +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.Internal +import Numeric.GSL.Internal + +import Foreign.Ptr(FunPtr, nullFunPtr, freeHaskellFunPtr) +import Foreign.C.Types +import System.IO.Unsafe(unsafePerformIO) + +------------------------------------------------------------------------- + +type Jacobian = Double -> Vector Double -> Matrix Double + +-- | Stepping functions +data ODEMethod = RK2 -- ^ Embedded Runge-Kutta (2, 3) method. + | RK4 -- ^ 4th order (classical) Runge-Kutta. The error estimate is obtained by halving the step-size. For more efficient estimate of the error, use 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/hmatrix/src/Numeric/GSL/Polynomials.hs b/packages/hmatrix/src/Numeric/GSL/Polynomials.hs new file mode 100644 index 0000000..290c615 --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/Polynomials.hs @@ -0,0 +1,58 @@ +{-# LANGUAGE CPP, ForeignFunctionInterface #-} +----------------------------------------------------------------------------- +{- | +Module : Numeric.GSL.Polynomials +Copyright : (c) Alberto Ruiz 2006 +License : GPL-style + +Maintainer : Alberto Ruiz (aruiz at um dot es) +Stability : provisional +Portability : uses ffi + +Polynomials. + + + +-} +----------------------------------------------------------------------------- +module Numeric.GSL.Polynomials ( + polySolve +) where + +import Data.Packed.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:: TVCV diff --git a/packages/hmatrix/src/Numeric/GSL/Root.hs b/packages/hmatrix/src/Numeric/GSL/Root.hs new file mode 100644 index 0000000..9d561c4 --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/Root.hs @@ -0,0 +1,199 @@ +{- | +Module : Numeric.GSL.Root +Copyright : (c) Alberto Ruiz 2009 +License : GPL + +Maintainer : Alberto Ruiz (aruiz at um dot es) +Stability : provisional +Portability : uses ffi + +Multidimensional root finding. + + + +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.Internal +import Data.Packed.Matrix +import Numeric.GSL.Internal +import Foreign.Ptr(FunPtr, freeHaskellFunPtr) +import Foreign.C.Types +import System.IO.Unsafe(unsafePerformIO) + +------------------------------------------------------------------------- + +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 + +------------------------------------------------------------------------- +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 + +------------------------------------------------------------------------- + +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/hmatrix/src/Numeric/GSL/Vector.hs b/packages/hmatrix/src/Numeric/GSL/Vector.hs new file mode 100644 index 0000000..6204b8e --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/Vector.hs @@ -0,0 +1,328 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.GSL.Vector +-- Copyright : (c) Alberto Ruiz 2007 +-- License : GPL-style +-- +-- Maintainer : Alberto Ruiz +-- Stability : provisional +-- Portability : portable (uses FFI) +-- +-- Low level interface to vector operations. +-- +----------------------------------------------------------------------------- + +module Numeric.GSL.Vector ( + sumF, sumR, sumQ, sumC, + prodF, prodR, prodQ, prodC, + dotF, dotR, dotQ, dotC, + FunCodeS(..), toScalarR, toScalarF, toScalarC, toScalarQ, + FunCodeV(..), vectorMapR, vectorMapC, vectorMapF, vectorMapQ, + FunCodeSV(..), vectorMapValR, vectorMapValC, vectorMapValF, vectorMapValQ, + FunCodeVV(..), vectorZipR, vectorZipC, vectorZipF, vectorZipQ, + RandDist(..), randomVector +) where + +import Data.Packed.Internal.Common +import Data.Packed.Internal.Signatures +import Data.Packed.Internal.Vector + +import Data.Complex +import Foreign.Marshal.Alloc(free) +import Foreign.Marshal.Array(newArray) +import Foreign.Ptr(Ptr) +import Foreign.C.Types +import System.IO.Unsafe(unsafePerformIO) +import Control.Monad(when) + +fromei x = fromIntegral (fromEnum x) :: CInt + +data FunCodeV = Sin + | Cos + | Tan + | Abs + | ASin + | ACos + | ATan + | Sinh + | Cosh + | Tanh + | ASinh + | ACosh + | ATanh + | Exp + | Log + | Sign + | Sqrt + deriving Enum + +data FunCodeSV = Scale + | Recip + | AddConstant + | Negate + | PowSV + | PowVS + deriving Enum + +data FunCodeVV = Add + | Sub + | Mul + | Div + | Pow + | ATan2 + deriving Enum + +data FunCodeS = Norm2 + | AbsSum + | MaxIdx + | Max + | MinIdx + | Min + deriving Enum + +------------------------------------------------------------------ + +-- | sum of elements +sumF :: Vector Float -> Float +sumF x = unsafePerformIO $ do + r <- createVector 1 + app2 c_sumF vec x vec r "sumF" + return $ r @> 0 + +-- | sum of elements +sumR :: Vector Double -> Double +sumR x = unsafePerformIO $ do + r <- createVector 1 + app2 c_sumR vec x vec r "sumR" + return $ r @> 0 + +-- | sum of elements +sumQ :: Vector (Complex Float) -> Complex Float +sumQ x = unsafePerformIO $ do + r <- createVector 1 + app2 c_sumQ vec x vec r "sumQ" + return $ r @> 0 + +-- | sum of elements +sumC :: Vector (Complex Double) -> Complex Double +sumC x = unsafePerformIO $ do + r <- createVector 1 + app2 c_sumC vec x vec r "sumC" + return $ r @> 0 + +foreign import ccall unsafe "gsl-aux.h sumF" c_sumF :: TFF +foreign import ccall unsafe "gsl-aux.h sumR" c_sumR :: TVV +foreign import ccall unsafe "gsl-aux.h sumQ" c_sumQ :: TQVQV +foreign import ccall unsafe "gsl-aux.h sumC" c_sumC :: TCVCV + +-- | product of elements +prodF :: Vector Float -> Float +prodF x = unsafePerformIO $ do + r <- createVector 1 + app2 c_prodF vec x vec r "prodF" + return $ r @> 0 + +-- | product of elements +prodR :: Vector Double -> Double +prodR x = unsafePerformIO $ do + r <- createVector 1 + app2 c_prodR vec x vec r "prodR" + return $ r @> 0 + +-- | product of elements +prodQ :: Vector (Complex Float) -> Complex Float +prodQ x = unsafePerformIO $ do + r <- createVector 1 + app2 c_prodQ vec x vec r "prodQ" + return $ r @> 0 + +-- | product of elements +prodC :: Vector (Complex Double) -> Complex Double +prodC x = unsafePerformIO $ do + r <- createVector 1 + app2 c_prodC vec x vec r "prodC" + return $ r @> 0 + +foreign import ccall unsafe "gsl-aux.h prodF" c_prodF :: TFF +foreign import ccall unsafe "gsl-aux.h prodR" c_prodR :: TVV +foreign import ccall unsafe "gsl-aux.h prodQ" c_prodQ :: TQVQV +foreign import ccall unsafe "gsl-aux.h prodC" c_prodC :: TCVCV + +-- | dot product +dotF :: Vector Float -> Vector Float -> Float +dotF x y = unsafePerformIO $ do + r <- createVector 1 + app3 c_dotF vec x vec y vec r "dotF" + return $ r @> 0 + +-- | dot product +dotR :: Vector Double -> Vector Double -> Double +dotR x y = unsafePerformIO $ do + r <- createVector 1 + app3 c_dotR vec x vec y vec r "dotR" + return $ r @> 0 + +-- | dot product +dotQ :: Vector (Complex Float) -> Vector (Complex Float) -> Complex Float +dotQ x y = unsafePerformIO $ do + r <- createVector 1 + app3 c_dotQ vec x vec y vec r "dotQ" + return $ r @> 0 + +-- | dot product +dotC :: Vector (Complex Double) -> Vector (Complex Double) -> Complex Double +dotC x y = unsafePerformIO $ do + r <- createVector 1 + app3 c_dotC vec x vec y vec r "dotC" + return $ r @> 0 + +foreign import ccall unsafe "gsl-aux.h dotF" c_dotF :: TFFF +foreign import ccall unsafe "gsl-aux.h dotR" c_dotR :: TVVV +foreign import ccall unsafe "gsl-aux.h dotQ" c_dotQ :: TQVQVQV +foreign import ccall unsafe "gsl-aux.h dotC" c_dotC :: TCVCVCV + +------------------------------------------------------------------ + +toScalarAux fun code v = unsafePerformIO $ do + r <- createVector 1 + app2 (fun (fromei code)) vec v vec r "toScalarAux" + return (r `at` 0) + +vectorMapAux fun code v = unsafePerformIO $ do + r <- createVector (dim v) + app2 (fun (fromei code)) vec v vec r "vectorMapAux" + return r + +vectorMapValAux fun code val v = unsafePerformIO $ do + r <- createVector (dim v) + pval <- newArray [val] + app2 (fun (fromei code) pval) vec v vec r "vectorMapValAux" + free pval + return r + +vectorZipAux fun code u v = unsafePerformIO $ do + r <- createVector (dim u) + when (dim u > 0) $ app3 (fun (fromei code)) vec u vec v vec r "vectorZipAux" + return r + +--------------------------------------------------------------------- + +-- | obtains different functions of a vector: norm1, norm2, max, min, posmax, posmin, etc. +toScalarR :: FunCodeS -> Vector Double -> Double +toScalarR oper = toScalarAux c_toScalarR (fromei oper) + +foreign import ccall unsafe "gsl-aux.h toScalarR" c_toScalarR :: CInt -> TVV + +-- | obtains different functions of a vector: norm1, norm2, max, min, posmax, posmin, etc. +toScalarF :: FunCodeS -> Vector Float -> Float +toScalarF oper = toScalarAux c_toScalarF (fromei oper) + +foreign import ccall unsafe "gsl-aux.h toScalarF" c_toScalarF :: CInt -> TFF + +-- | obtains different functions of a vector: only norm1, norm2 +toScalarC :: FunCodeS -> Vector (Complex Double) -> Double +toScalarC oper = toScalarAux c_toScalarC (fromei oper) + +foreign import ccall unsafe "gsl-aux.h toScalarC" c_toScalarC :: CInt -> TCVV + +-- | obtains different functions of a vector: only norm1, norm2 +toScalarQ :: FunCodeS -> Vector (Complex Float) -> Float +toScalarQ oper = toScalarAux c_toScalarQ (fromei oper) + +foreign import ccall unsafe "gsl-aux.h toScalarQ" c_toScalarQ :: CInt -> TQVF + +------------------------------------------------------------------ + +-- | map of real vectors with given function +vectorMapR :: FunCodeV -> Vector Double -> Vector Double +vectorMapR = vectorMapAux c_vectorMapR + +foreign import ccall unsafe "gsl-aux.h mapR" c_vectorMapR :: CInt -> TVV + +-- | map of complex vectors with given function +vectorMapC :: FunCodeV -> Vector (Complex Double) -> Vector (Complex Double) +vectorMapC oper = vectorMapAux c_vectorMapC (fromei oper) + +foreign import ccall unsafe "gsl-aux.h mapC" c_vectorMapC :: CInt -> TCVCV + +-- | map of real vectors with given function +vectorMapF :: FunCodeV -> Vector Float -> Vector Float +vectorMapF = vectorMapAux c_vectorMapF + +foreign import ccall unsafe "gsl-aux.h mapF" c_vectorMapF :: CInt -> TFF + +-- | map of real vectors with given function +vectorMapQ :: FunCodeV -> Vector (Complex Float) -> Vector (Complex Float) +vectorMapQ = vectorMapAux c_vectorMapQ + +foreign import ccall unsafe "gsl-aux.h mapQ" c_vectorMapQ :: CInt -> TQVQV + +------------------------------------------------------------------- + +-- | map of real vectors with given function +vectorMapValR :: FunCodeSV -> Double -> Vector Double -> Vector Double +vectorMapValR oper = vectorMapValAux c_vectorMapValR (fromei oper) + +foreign import ccall unsafe "gsl-aux.h mapValR" c_vectorMapValR :: CInt -> Ptr Double -> TVV + +-- | map of complex vectors with given function +vectorMapValC :: FunCodeSV -> Complex Double -> Vector (Complex Double) -> Vector (Complex Double) +vectorMapValC = vectorMapValAux c_vectorMapValC + +foreign import ccall unsafe "gsl-aux.h mapValC" c_vectorMapValC :: CInt -> Ptr (Complex Double) -> TCVCV + +-- | map of real vectors with given function +vectorMapValF :: FunCodeSV -> Float -> Vector Float -> Vector Float +vectorMapValF oper = vectorMapValAux c_vectorMapValF (fromei oper) + +foreign import ccall unsafe "gsl-aux.h mapValF" c_vectorMapValF :: CInt -> Ptr Float -> TFF + +-- | map of complex vectors with given function +vectorMapValQ :: FunCodeSV -> Complex Float -> Vector (Complex Float) -> Vector (Complex Float) +vectorMapValQ oper = vectorMapValAux c_vectorMapValQ (fromei oper) + +foreign import ccall unsafe "gsl-aux.h mapValQ" c_vectorMapValQ :: CInt -> Ptr (Complex Float) -> TQVQV + +------------------------------------------------------------------- + +-- | elementwise operation on real vectors +vectorZipR :: FunCodeVV -> Vector Double -> Vector Double -> Vector Double +vectorZipR = vectorZipAux c_vectorZipR + +foreign import ccall unsafe "gsl-aux.h zipR" c_vectorZipR :: CInt -> TVVV + +-- | elementwise operation on complex vectors +vectorZipC :: FunCodeVV -> Vector (Complex Double) -> Vector (Complex Double) -> Vector (Complex Double) +vectorZipC = vectorZipAux c_vectorZipC + +foreign import ccall unsafe "gsl-aux.h zipC" c_vectorZipC :: CInt -> TCVCVCV + +-- | elementwise operation on real vectors +vectorZipF :: FunCodeVV -> Vector Float -> Vector Float -> Vector Float +vectorZipF = vectorZipAux c_vectorZipF + +foreign import ccall unsafe "gsl-aux.h zipF" c_vectorZipF :: CInt -> TFFF + +-- | elementwise operation on complex vectors +vectorZipQ :: FunCodeVV -> Vector (Complex Float) -> Vector (Complex Float) -> Vector (Complex Float) +vectorZipQ = vectorZipAux c_vectorZipQ + +foreign import ccall unsafe "gsl-aux.h zipQ" c_vectorZipQ :: CInt -> TQVQVQV + +----------------------------------------------------------------------- + +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 diff --git a/packages/hmatrix/src/Numeric/GSL/gsl-aux.c b/packages/hmatrix/src/Numeric/GSL/gsl-aux.c new file mode 100644 index 0000000..410d157 --- /dev/null +++ b/packages/hmatrix/src/Numeric/GSL/gsl-aux.c @@ -0,0 +1,1541 @@ +#include + +#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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#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 sumF(KFVEC(x),FVEC(r)) { + DEBUGMSG("sumF"); + REQUIRES(rn==1,BAD_SIZE); + int i; + float res = 0; + for (i = 0; i < xn; i++) res += xp[i]; + rp[0] = res; + OK +} + +int sumR(KRVEC(x),RVEC(r)) { + DEBUGMSG("sumR"); + REQUIRES(rn==1,BAD_SIZE); + int i; + double res = 0; + for (i = 0; i < xn; i++) res += xp[i]; + rp[0] = res; + OK +} + +int sumQ(KQVEC(x),QVEC(r)) { + DEBUGMSG("sumQ"); + REQUIRES(rn==1,BAD_SIZE); + int i; + gsl_complex_float res; + res.dat[0] = 0; + res.dat[1] = 0; + for (i = 0; i < xn; i++) { + res.dat[0] += xp[i].dat[0]; + res.dat[1] += xp[i].dat[1]; + } + rp[0] = res; + OK +} + +int sumC(KCVEC(x),CVEC(r)) { + DEBUGMSG("sumC"); + REQUIRES(rn==1,BAD_SIZE); + int i; + gsl_complex res; + res.dat[0] = 0; + res.dat[1] = 0; + for (i = 0; i < xn; i++) { + res.dat[0] += xp[i].dat[0]; + res.dat[1] += xp[i].dat[1]; + } + rp[0] = res; + OK +} + +int prodF(KFVEC(x),FVEC(r)) { + DEBUGMSG("prodF"); + REQUIRES(rn==1,BAD_SIZE); + int i; + float res = 1; + for (i = 0; i < xn; i++) res *= xp[i]; + rp[0] = res; + OK +} + +int prodR(KRVEC(x),RVEC(r)) { + DEBUGMSG("prodR"); + REQUIRES(rn==1,BAD_SIZE); + int i; + double res = 1; + for (i = 0; i < xn; i++) res *= xp[i]; + rp[0] = res; + OK +} + +int prodQ(KQVEC(x),QVEC(r)) { + DEBUGMSG("prodQ"); + REQUIRES(rn==1,BAD_SIZE); + int i; + gsl_complex_float res; + float temp; + res.dat[0] = 1; + res.dat[1] = 0; + for (i = 0; i < xn; i++) { + temp = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1]; + res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0]; + res.dat[0] = temp; + } + rp[0] = res; + OK +} + +int prodC(KCVEC(x),CVEC(r)) { + DEBUGMSG("prodC"); + REQUIRES(rn==1,BAD_SIZE); + int i; + gsl_complex res; + double temp; + res.dat[0] = 1; + res.dat[1] = 0; + for (i = 0; i < xn; i++) { + temp = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1]; + res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0]; + res.dat[0] = temp; + } + rp[0] = res; + OK +} + +int dotF(KFVEC(x), KFVEC(y), FVEC(r)) { + DEBUGMSG("dotF"); + REQUIRES(xn==yn,BAD_SIZE); + REQUIRES(rn==1,BAD_SIZE); + DEBUGMSG("dotF"); + KFVVIEW(x); + KFVVIEW(y); + gsl_blas_sdot(V(x),V(y),rp); + OK +} + +int dotR(KRVEC(x), KRVEC(y), RVEC(r)) { + DEBUGMSG("dotR"); + REQUIRES(xn==yn,BAD_SIZE); + REQUIRES(rn==1,BAD_SIZE); + DEBUGMSG("dotR"); + KDVVIEW(x); + KDVVIEW(y); + gsl_blas_ddot(V(x),V(y),rp); + OK +} + +int dotQ(KQVEC(x), KQVEC(y), QVEC(r)) { + DEBUGMSG("dotQ"); + REQUIRES(xn==yn,BAD_SIZE); + REQUIRES(rn==1,BAD_SIZE); + DEBUGMSG("dotQ"); + KQVVIEW(x); + KQVVIEW(y); + gsl_blas_cdotu(V(x),V(y),rp); + OK +} + +int dotC(KCVEC(x), KCVEC(y), CVEC(r)) { + DEBUGMSG("dotC"); + REQUIRES(xn==yn,BAD_SIZE); + REQUIRES(rn==1,BAD_SIZE); + DEBUGMSG("dotC"); + KCVVIEW(x); + KCVVIEW(y); + gsl_blas_zdotu(V(x),V(y),rp); + OK +} + +int toScalarR(int code, KRVEC(x), RVEC(r)) { + REQUIRES(rn==1,BAD_SIZE); + DEBUGMSG("toScalarR"); + KDVVIEW(x); + double res; + switch(code) { + case 0: { res = gsl_blas_dnrm2(V(x)); break; } + case 1: { res = gsl_blas_dasum(V(x)); break; } + case 2: { res = gsl_vector_max_index(V(x)); break; } + case 3: { res = gsl_vector_max(V(x)); break; } + case 4: { res = gsl_vector_min_index(V(x)); break; } + case 5: { res = gsl_vector_min(V(x)); break; } + default: ERROR(BAD_CODE); + } + rp[0] = res; + OK +} + +int toScalarF(int code, KFVEC(x), FVEC(r)) { + REQUIRES(rn==1,BAD_SIZE); + DEBUGMSG("toScalarF"); + KFVVIEW(x); + float res; + switch(code) { + case 0: { res = gsl_blas_snrm2(V(x)); break; } + case 1: { res = gsl_blas_sasum(V(x)); break; } + case 2: { res = gsl_vector_float_max_index(V(x)); break; } + case 3: { res = gsl_vector_float_max(V(x)); break; } + case 4: { res = gsl_vector_float_min_index(V(x)); break; } + case 5: { res = gsl_vector_float_min(V(x)); break; } + default: ERROR(BAD_CODE); + } + rp[0] = res; + OK +} + + +int toScalarC(int code, KCVEC(x), RVEC(r)) { + REQUIRES(rn==1,BAD_SIZE); + DEBUGMSG("toScalarC"); + KCVVIEW(x); + double res; + switch(code) { + case 0: { res = gsl_blas_dznrm2(V(x)); break; } + case 1: { res = gsl_blas_dzasum(V(x)); break; } + default: ERROR(BAD_CODE); + } + rp[0] = res; + OK +} + +int toScalarQ(int code, KQVEC(x), FVEC(r)) { + REQUIRES(rn==1,BAD_SIZE); + DEBUGMSG("toScalarQ"); + KQVVIEW(x); + float res; + switch(code) { + case 0: { res = gsl_blas_scnrm2(V(x)); break; } + case 1: { res = gsl_blas_scasum(V(x)); break; } + default: ERROR(BAD_CODE); + } + rp[0] = res; + OK +} + + +inline double sign(double x) { + if(x>0) { + return +1.0; + } else if (x<0) { + return -1.0; + } else { + return 0.0; + } +} + +inline float float_sign(float x) { + if(x>0) { + return +1.0; + } else if (x<0) { + return -1.0; + } else { + return 0.0; + } +} + +inline gsl_complex complex_abs(gsl_complex z) { + gsl_complex r; + r.dat[0] = gsl_complex_abs(z); + r.dat[1] = 0; + return r; +} + +inline gsl_complex complex_signum(gsl_complex z) { + gsl_complex r; + double mag; + if (z.dat[0] == 0 && z.dat[1] == 0) { + r.dat[0] = 0; + r.dat[1] = 0; + } else { + mag = gsl_complex_abs(z); + r.dat[0] = z.dat[0]/mag; + r.dat[1] = z.dat[1]/mag; + } + return r; +} + +#define OP(C,F) case C: { for(k=0;k1,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; isize,sizeof(double)); + int k; + for(k=0;ksize;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; ifval; + solp[iter*solc+2] = size; + + int k; + for(k=0;kx,k); + } + iter++; + if (status) break; + status = gsl_multimin_test_size (size, tolsize); + } while (status == GSL_CONTINUE && iter < maxit); + int i,j; + for (i=iter; isize,sizeof(double)); + int k; + for(k=0;ksize;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;ksize;k++) { + p[k] = gsl_vector_get(x,k); + } + + fdf->df(x->size,p,g->size,q); + + for(k=0;ksize;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;kx,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; if)(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; isize,sizeof(double)); + double* q = (double*)calloc(y->size,sizeof(double)); + int k; + for(k=0;ksize;k++) { + p[k] = gsl_vector_get(x,k); + } + f(x->size,p,y->size,q); + for(k=0;ksize;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;kx,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; isize,sizeof(double)); + double* q = (double*)calloc(y->size,sizeof(double)); + int k; + for(k=0;ksize;k++) { + p[k] = gsl_vector_get(x,k); + } + (fjf->f)(x->size,p,y->size,q); + for(k=0;ksize;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;ksize;k++) { + p[k] = gsl_vector_get(x,k); + } + + (fjf->jf)(x->size,p,jac->size1,jac->size2,q); + + k=0; + for(i=0;isize1;i++) { + for(j=0;jsize2;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;kx,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; if); + + int k; + for(k=0;kx,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; iJ, 0.0, M(cov)); + + gsl_multifit_fdfsolver_free (s); + OK +} + + +////////////////////////////////////////////////////// + + +#define RAN(C,F) case C: { for(k=0;k + +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 + +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