summaryrefslogtreecommitdiff
path: root/packages/hmatrix/src/Numeric/GSL
diff options
context:
space:
mode:
Diffstat (limited to 'packages/hmatrix/src/Numeric/GSL')
-rw-r--r--packages/hmatrix/src/Numeric/GSL/Differentiation.hs87
-rw-r--r--packages/hmatrix/src/Numeric/GSL/Fitting.hs179
-rw-r--r--packages/hmatrix/src/Numeric/GSL/Fourier.hs47
-rw-r--r--packages/hmatrix/src/Numeric/GSL/Integration.hs250
-rw-r--r--packages/hmatrix/src/Numeric/GSL/Internal.hs76
-rw-r--r--packages/hmatrix/src/Numeric/GSL/Minimization.hs227
-rw-r--r--packages/hmatrix/src/Numeric/GSL/ODE.hs138
-rw-r--r--packages/hmatrix/src/Numeric/GSL/Polynomials.hs58
-rw-r--r--packages/hmatrix/src/Numeric/GSL/Root.hs199
-rw-r--r--packages/hmatrix/src/Numeric/GSL/Vector.hs328
-rw-r--r--packages/hmatrix/src/Numeric/GSL/gsl-aux.c1541
-rw-r--r--packages/hmatrix/src/Numeric/GSL/gsl-ode.c182
12 files changed, 3312 insertions, 0 deletions
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 @@
1{-# OPTIONS #-}
2-----------------------------------------------------------------------------
3{- |
4Module : Numeric.GSL.Differentiation
5Copyright : (c) Alberto Ruiz 2006
6License : GPL-style
7
8Maintainer : Alberto Ruiz (aruiz at um dot es)
9Stability : provisional
10Portability : uses ffi
11
12Numerical differentiation.
13
14<http://www.gnu.org/software/gsl/manual/html_node/Numerical-Differentiation.html#Numerical-Differentiation>
15
16From 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.\"
17-}
18-----------------------------------------------------------------------------
19module Numeric.GSL.Differentiation (
20 derivCentral,
21 derivForward,
22 derivBackward
23) where
24
25import Foreign.C.Types
26import Foreign.Marshal.Alloc(malloc, free)
27import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr)
28import Foreign.Storable(peek)
29import Data.Packed.Internal(check,(//))
30import System.IO.Unsafe(unsafePerformIO)
31
32derivGen ::
33 CInt -- ^ type: 0 central, 1 forward, 2 backward
34 -> Double -- ^ initial step size
35 -> (Double -> Double) -- ^ function
36 -> Double -- ^ point where the derivative is taken
37 -> (Double, Double) -- ^ result and error
38derivGen c h f x = unsafePerformIO $ do
39 r <- malloc
40 e <- malloc
41 fp <- mkfun (\y _ -> f y)
42 c_deriv c fp x h r e // check "deriv"
43 vr <- peek r
44 ve <- peek e
45 let result = (vr,ve)
46 free r
47 free e
48 freeHaskellFunPtr fp
49 return result
50
51foreign import ccall safe "gsl-aux.h deriv"
52 c_deriv :: CInt -> FunPtr (Double -> Ptr () -> Double) -> Double -> Double
53 -> Ptr Double -> Ptr Double -> IO CInt
54
55
56{- | Adaptive central difference algorithm, /gsl_deriv_central/. For example:
57
58>>> let deriv = derivCentral 0.01
59>>> deriv sin (pi/4)
60(0.7071067812000676,1.0600063101654055e-10)
61>>> cos (pi/4)
620.7071067811865476
63
64-}
65derivCentral :: Double -- ^ initial step size
66 -> (Double -> Double) -- ^ function
67 -> Double -- ^ point where the derivative is taken
68 -> (Double, Double) -- ^ result and absolute error
69derivCentral = derivGen 0
70
71-- | 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.
72derivForward :: Double -- ^ initial step size
73 -> (Double -> Double) -- ^ function
74 -> Double -- ^ point where the derivative is taken
75 -> (Double, Double) -- ^ result and absolute error
76derivForward = derivGen 1
77
78-- | Adaptive backward difference algorithm, /gsl_deriv_backward/.
79derivBackward ::Double -- ^ initial step size
80 -> (Double -> Double) -- ^ function
81 -> Double -- ^ point where the derivative is taken
82 -> (Double, Double) -- ^ result and absolute error
83derivBackward = derivGen 2
84
85{- | conversion of Haskell functions into function pointers that can be used in the C side
86-}
87foreign 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 @@
1{- |
2Module : Numeric.GSL.Fitting
3Copyright : (c) Alberto Ruiz 2010
4License : GPL
5
6Maintainer : Alberto Ruiz (aruiz at um dot es)
7Stability : provisional
8Portability : uses ffi
9
10Nonlinear Least-Squares Fitting
11
12<http://www.gnu.org/software/gsl/manual/html_node/Nonlinear-Least_002dSquares-Fitting.html>
13
14The example program in the GSL manual (see examples/fitting.hs):
15
16@
17dat = [
18 ([0.0],([6.0133918608118675],0.1)),
19 ([1.0],([5.5153769909966535],0.1)),
20 ([2.0],([5.261094606015287],0.1)),
21 ...
22 ([39.0],([1.0619821710802808],0.1))]
23
24expModel [a,lambda,b] [t] = [a * exp (-lambda * t) + b]
25
26expModelDer [a,lambda,b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]]
27
28(sol,path) = fitModelScaled 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0]
29@
30
31>>> path
32(6><5)
33 [ 1.0, 76.45780563978782, 1.6465931240727802, 1.8147715267618197e-2, 0.6465931240727797
34 , 2.0, 37.683816318260355, 2.858760367632973, 8.092094813253975e-2, 1.4479636296208662
35 , 3.0, 9.5807893736187, 4.948995119561291, 0.11942927999921617, 1.0945766509238248
36 , 4.0, 5.630494933603935, 5.021755718065913, 0.10287787128056883, 1.0338835440862608
37 , 5.0, 5.443976278682909, 5.045204331329302, 0.10405523433131504, 1.019416067207375
38 , 6.0, 5.4439736648994685, 5.045357818922331, 0.10404905846029407, 1.0192487112786812 ]
39>>> sol
40[(5.045357818922331,6.027976702418132e-2),
41(0.10404905846029407,3.157045047172834e-3),
42(1.0192487112786812,3.782067731353722e-2)]
43
44-}
45-----------------------------------------------------------------------------
46
47module Numeric.GSL.Fitting (
48 -- * Levenberg-Marquardt
49 nlFitting, FittingMethod(..),
50 -- * Utilities
51 fitModelScaled, fitModel
52) where
53
54import Data.Packed.Internal
55import Numeric.LinearAlgebra
56import Numeric.GSL.Internal
57
58import Foreign.Ptr(FunPtr, freeHaskellFunPtr)
59import Foreign.C.Types
60import System.IO.Unsafe(unsafePerformIO)
61
62-------------------------------------------------------------------------
63
64data 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.
65 | 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.
66 deriving (Enum,Eq,Show,Bounded)
67
68
69-- | Nonlinear multidimensional least-squares fitting.
70nlFitting :: FittingMethod
71 -> Double -- ^ absolute tolerance
72 -> Double -- ^ relative tolerance
73 -> Int -- ^ maximum number of iterations allowed
74 -> (Vector Double -> Vector Double) -- ^ function to be minimized
75 -> (Vector Double -> Matrix Double) -- ^ Jacobian
76 -> Vector Double -- ^ starting point
77 -> (Vector Double, Matrix Double) -- ^ solution vector and optimization path
78
79nlFitting method epsabs epsrel maxit fun jac xinit = nlFitGen (fi (fromEnum method)) fun jac xinit epsabs epsrel maxit
80
81nlFitGen m f jac xiv epsabs epsrel maxit = unsafePerformIO $ do
82 let p = dim xiv
83 n = dim (f xiv)
84 fp <- mkVecVecfun (aux_vTov (checkdim1 n p . f))
85 jp <- mkVecMatfun (aux_vTom (checkdim2 n p . jac))
86 rawpath <- createMatrix RowMajor maxit (2+p)
87 app2 (c_nlfit m fp jp epsabs epsrel (fi maxit) (fi n)) vec xiv mat rawpath "c_nlfit"
88 let it = round (rawpath @@> (maxit-1,0))
89 path = takeRows it rawpath
90 [sol] = toRows $ dropRows (it-1) path
91 freeHaskellFunPtr fp
92 freeHaskellFunPtr jp
93 return (subVector 2 p sol, path)
94
95foreign import ccall safe "nlfit"
96 c_nlfit:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> Double -> CInt -> CInt -> TVM
97
98-------------------------------------------------------
99
100checkdim1 n _p v
101 | dim v == n = v
102 | otherwise = error $ "Error: "++ show n
103 ++ " components expected in the result of the function supplied to nlFitting"
104
105checkdim2 n p m
106 | rows m == n && cols m == p = m
107 | otherwise = error $ "Error: "++ show n ++ "x" ++ show p
108 ++ " Jacobian expected in nlFitting"
109
110------------------------------------------------------------
111
112err (model,deriv) dat vsol = zip sol errs where
113 sol = toList vsol
114 c = max 1 (chi/sqrt (fromIntegral dof))
115 dof = length dat - (rows cov)
116 chi = norm2 (fromList $ cost (resMs model) dat sol)
117 js = fromLists $ jacobian (resDs deriv) dat sol
118 cov = inv $ trans js <> js
119 errs = toList $ scalar c * sqrt (takeDiag cov)
120
121
122
123-- | Higher level interface to 'nlFitting' 'LevenbergMarquardtScaled'. The optimization function and
124-- Jacobian are automatically built from a model f vs x = y and its derivatives, and a list of
125-- instances (x, (y,sigma)) to be fitted.
126
127fitModelScaled
128 :: Double -- ^ absolute tolerance
129 -> Double -- ^ relative tolerance
130 -> Int -- ^ maximum number of iterations allowed
131 -> ([Double] -> x -> [Double], [Double] -> x -> [[Double]]) -- ^ (model, derivatives)
132 -> [(x, ([Double], Double))] -- ^ instances
133 -> [Double] -- ^ starting point
134 -> ([(Double, Double)], Matrix Double) -- ^ (solution, error) and optimization path
135fitModelScaled epsabs epsrel maxit (model,deriv) dt xin = (err (model,deriv) dt sol, path) where
136 (sol,path) = nlFitting LevenbergMarquardtScaled epsabs epsrel maxit
137 (fromList . cost (resMs model) dt . toList)
138 (fromLists . jacobian (resDs deriv) dt . toList)
139 (fromList xin)
140
141
142
143-- | Higher level interface to 'nlFitting' 'LevenbergMarquardt'. The optimization function and
144-- Jacobian are automatically built from a model f vs x = y and its derivatives, and a list of
145-- instances (x,y) to be fitted.
146
147fitModel :: Double -- ^ absolute tolerance
148 -> Double -- ^ relative tolerance
149 -> Int -- ^ maximum number of iterations allowed
150 -> ([Double] -> x -> [Double], [Double] -> x -> [[Double]]) -- ^ (model, derivatives)
151 -> [(x, [Double])] -- ^ instances
152 -> [Double] -- ^ starting point
153 -> ([Double], Matrix Double) -- ^ solution and optimization path
154fitModel epsabs epsrel maxit (model,deriv) dt xin = (toList sol, path) where
155 (sol,path) = nlFitting LevenbergMarquardt epsabs epsrel maxit
156 (fromList . cost (resM model) dt . toList)
157 (fromLists . jacobian (resD deriv) dt . toList)
158 (fromList xin)
159
160cost model ds vs = concatMap (model vs) ds
161
162jacobian modelDer ds vs = concatMap (modelDer vs) ds
163
164-- | Model-to-residual for association pairs with sigma, to be used with 'fitModel'.
165resMs :: ([Double] -> x -> [Double]) -> [Double] -> (x, ([Double], Double)) -> [Double]
166resMs m v = \(x,(ys,s)) -> zipWith (g s) (m v x) ys where g s a b = (a-b)/s
167
168-- | Associated derivative for 'resMs'.
169resDs :: ([Double] -> x -> [[Double]]) -> [Double] -> (x, ([Double], Double)) -> [[Double]]
170resDs m v = \(x,(_,s)) -> map (map (/s)) (m v x)
171
172-- | Model-to-residual for association pairs, to be used with 'fitModel'. It is equivalent
173-- to 'resMs' with all sigmas = 1.
174resM :: ([Double] -> x -> [Double]) -> [Double] -> (x, [Double]) -> [Double]
175resM m v = \(x,ys) -> zipWith g (m v x) ys where g a b = a-b
176
177-- | Associated derivative for 'resM'.
178resD :: ([Double] -> x -> [[Double]]) -> [Double] -> (x, [Double]) -> [[Double]]
179resD 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 @@
1{-# LANGUAGE ForeignFunctionInterface #-}
2-----------------------------------------------------------------------------
3{- |
4Module : Numeric.GSL.Fourier
5Copyright : (c) Alberto Ruiz 2006
6License : GPL-style
7
8Maintainer : Alberto Ruiz (aruiz at um dot es)
9Stability : provisional
10Portability : uses ffi
11
12Fourier Transform.
13
14<http://www.gnu.org/software/gsl/manual/html_node/Fast-Fourier-Transforms.html#Fast-Fourier-Transforms>
15
16-}
17-----------------------------------------------------------------------------
18module Numeric.GSL.Fourier (
19 fft,
20 ifft
21) where
22
23import Data.Packed.Internal
24import Data.Complex
25import Foreign.C.Types
26import System.IO.Unsafe (unsafePerformIO)
27
28genfft code v = unsafePerformIO $ do
29 r <- createVector (dim v)
30 app2 (c_fft code) vec v vec r "fft"
31 return r
32
33foreign import ccall unsafe "gsl-aux.h fft" c_fft :: CInt -> TCVCV
34
35
36{- | Fast 1D Fourier transform of a 'Vector' @(@'Complex' 'Double'@)@ using /gsl_fft_complex_forward/. It uses the same scaling conventions as GNU Octave.
37
38>>> fft (fromList [1,2,3,4])
39fromList [10.0 :+ 0.0,(-2.0) :+ 2.0,(-2.0) :+ 0.0,(-2.0) :+ (-2.0)]
40
41-}
42fft :: Vector (Complex Double) -> Vector (Complex Double)
43fft = genfft 0
44
45-- | The inverse of 'fft', using /gsl_fft_complex_inverse/.
46ifft :: Vector (Complex Double) -> Vector (Complex Double)
47ifft = 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 @@
1{-# OPTIONS #-}
2-----------------------------------------------------------------------------
3{- |
4Module : Numeric.GSL.Integration
5Copyright : (c) Alberto Ruiz 2006
6License : GPL-style
7
8Maintainer : Alberto Ruiz (aruiz at um dot es)
9Stability : provisional
10Portability : uses ffi
11
12Numerical integration routines.
13
14<http://www.gnu.org/software/gsl/manual/html_node/Numerical-Integration.html#Numerical-Integration>
15-}
16-----------------------------------------------------------------------------
17
18module Numeric.GSL.Integration (
19 integrateQNG,
20 integrateQAGS,
21 integrateQAGI,
22 integrateQAGIU,
23 integrateQAGIL,
24 integrateCQUAD
25) where
26
27import Foreign.C.Types
28import Foreign.Marshal.Alloc(malloc, free)
29import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr)
30import Foreign.Storable(peek)
31import Data.Packed.Internal(check,(//))
32import System.IO.Unsafe(unsafePerformIO)
33
34eps = 1e-12
35
36{- | conversion of Haskell functions into function pointers that can be used in the C side
37-}
38foreign import ccall safe "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double))
39
40--------------------------------------------------------------------
41{- | Numerical integration using /gsl_integration_qags/ (adaptive integration with singularities). For example:
42
43>>> let quad = integrateQAGS 1E-9 1000
44>>> let f a x = x**(-0.5) * log (a*x)
45>>> quad (f 1) 0 1
46(-3.999999999999974,4.871658632055187e-13)
47
48-}
49
50integrateQAGS :: Double -- ^ precision (e.g. 1E-9)
51 -> Int -- ^ size of auxiliary workspace (e.g. 1000)
52 -> (Double -> Double) -- ^ function to be integrated on the interval (a,b)
53 -> Double -- ^ a
54 -> Double -- ^ b
55 -> (Double, Double) -- ^ result of the integration and error
56integrateQAGS prec n f a b = unsafePerformIO $ do
57 r <- malloc
58 e <- malloc
59 fp <- mkfun (\x _ -> f x)
60 c_integrate_qags fp a b eps prec (fromIntegral n) r e // check "integrate_qags"
61 vr <- peek r
62 ve <- peek e
63 let result = (vr,ve)
64 free r
65 free e
66 freeHaskellFunPtr fp
67 return result
68
69foreign import ccall safe "integrate_qags" c_integrate_qags
70 :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
71 -> Double -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt
72
73-----------------------------------------------------------------
74{- | Numerical integration using /gsl_integration_qng/ (useful for fast integration of smooth functions). For example:
75
76>>> let quad = integrateQNG 1E-6
77>>> quad (\x -> 4/(1+x*x)) 0 1
78(3.141592653589793,3.487868498008632e-14)
79
80-}
81
82integrateQNG :: Double -- ^ precision (e.g. 1E-9)
83 -> (Double -> Double) -- ^ function to be integrated on the interval (a,b)
84 -> Double -- ^ a
85 -> Double -- ^ b
86 -> (Double, Double) -- ^ result of the integration and error
87integrateQNG prec f a b = unsafePerformIO $ do
88 r <- malloc
89 e <- malloc
90 fp <- mkfun (\x _ -> f x)
91 c_integrate_qng fp a b eps prec r e // check "integrate_qng"
92 vr <- peek r
93 ve <- peek e
94 let result = (vr,ve)
95 free r
96 free e
97 freeHaskellFunPtr fp
98 return result
99
100foreign import ccall safe "integrate_qng" c_integrate_qng
101 :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
102 -> Double -> Double -> Ptr Double -> Ptr Double -> IO CInt
103
104--------------------------------------------------------------------
105{- | Numerical integration using /gsl_integration_qagi/ (integration over the infinite integral -Inf..Inf using QAGS).
106For example:
107
108>>> let quad = integrateQAGI 1E-9 1000
109>>> let f a x = exp(-a * x^2)
110>>> quad (f 0.5)
111(2.5066282746310002,6.229215880648858e-11)
112
113-}
114
115integrateQAGI :: Double -- ^ precision (e.g. 1E-9)
116 -> Int -- ^ size of auxiliary workspace (e.g. 1000)
117 -> (Double -> Double) -- ^ function to be integrated on the interval (-Inf,Inf)
118 -> (Double, Double) -- ^ result of the integration and error
119integrateQAGI prec n f = unsafePerformIO $ do
120 r <- malloc
121 e <- malloc
122 fp <- mkfun (\x _ -> f x)
123 c_integrate_qagi fp eps prec (fromIntegral n) r e // check "integrate_qagi"
124 vr <- peek r
125 ve <- peek e
126 let result = (vr,ve)
127 free r
128 free e
129 freeHaskellFunPtr fp
130 return result
131
132foreign import ccall safe "integrate_qagi" c_integrate_qagi
133 :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
134 -> CInt -> Ptr Double -> Ptr Double -> IO CInt
135
136--------------------------------------------------------------------
137{- | Numerical integration using /gsl_integration_qagiu/ (integration over the semi-infinite integral a..Inf).
138For example:
139
140>>> let quad = integrateQAGIU 1E-9 1000
141>>> let f a x = exp(-a * x^2)
142>>> quad (f 0.5) 0
143(1.2533141373155001,3.114607940324429e-11)
144
145-}
146
147integrateQAGIU :: Double -- ^ precision (e.g. 1E-9)
148 -> Int -- ^ size of auxiliary workspace (e.g. 1000)
149 -> (Double -> Double) -- ^ function to be integrated on the interval (a,Inf)
150 -> Double -- ^ a
151 -> (Double, Double) -- ^ result of the integration and error
152integrateQAGIU prec n f a = unsafePerformIO $ do
153 r <- malloc
154 e <- malloc
155 fp <- mkfun (\x _ -> f x)
156 c_integrate_qagiu fp a eps prec (fromIntegral n) r e // check "integrate_qagiu"
157 vr <- peek r
158 ve <- peek e
159 let result = (vr,ve)
160 free r
161 free e
162 freeHaskellFunPtr fp
163 return result
164
165foreign import ccall safe "integrate_qagiu" c_integrate_qagiu
166 :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
167 -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt
168
169--------------------------------------------------------------------
170{- | Numerical integration using /gsl_integration_qagil/ (integration over the semi-infinite integral -Inf..b).
171For example:
172
173>>> let quad = integrateQAGIL 1E-9 1000
174>>> let f a x = exp(-a * x^2)
175>>> quad (f 0.5) 0
176(1.2533141373155001,3.114607940324429e-11)
177
178-}
179
180integrateQAGIL :: Double -- ^ precision (e.g. 1E-9)
181 -> Int -- ^ size of auxiliary workspace (e.g. 1000)
182 -> (Double -> Double) -- ^ function to be integrated on the interval (a,Inf)
183 -> Double -- ^ b
184 -> (Double, Double) -- ^ result of the integration and error
185integrateQAGIL prec n f b = unsafePerformIO $ do
186 r <- malloc
187 e <- malloc
188 fp <- mkfun (\x _ -> f x)
189 c_integrate_qagil fp b eps prec (fromIntegral n) r e // check "integrate_qagil"
190 vr <- peek r
191 ve <- peek e
192 let result = (vr,ve)
193 free r
194 free e
195 freeHaskellFunPtr fp
196 return result
197
198foreign import ccall safe "gsl-aux.h integrate_qagil" c_integrate_qagil
199 :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
200 -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt
201
202
203--------------------------------------------------------------------
204{- | Numerical integration using /gsl_integration_cquad/ (quadrature
205for general integrands). From the GSL manual:
206
207@CQUAD is a new doubly-adaptive general-purpose quadrature routine
208which can handle most types of singularities, non-numerical function
209values such as Inf or NaN, as well as some divergent integrals. It
210generally requires more function evaluations than the integration
211routines in QUADPACK, yet fails less often for difficult integrands.@
212
213For example:
214
215>>> let quad = integrateCQUAD 1E-12 1000
216>>> let f a x = exp(-a * x^2)
217>>> quad (f 0.5) 2 5
218(5.7025405463957006e-2,9.678874441303705e-16,95)
219
220Unlike other quadrature methods, integrateCQUAD also returns the
221number of function evaluations required.
222
223-}
224
225integrateCQUAD :: Double -- ^ precision (e.g. 1E-9)
226 -> Int -- ^ size of auxiliary workspace (e.g. 1000)
227 -> (Double -> Double) -- ^ function to be integrated on the interval (a, b)
228 -> Double -- ^ a
229 -> Double -- ^ b
230 -> (Double, Double, Int) -- ^ result of the integration, error and number of function evaluations performed
231integrateCQUAD prec n f a b = unsafePerformIO $ do
232 r <- malloc
233 e <- malloc
234 neval <- malloc
235 fp <- mkfun (\x _ -> f x)
236 c_integrate_cquad fp a b eps prec (fromIntegral n) r e neval // check "integrate_cquad"
237 vr <- peek r
238 ve <- peek e
239 vneval <- peek neval
240 let result = (vr,ve,vneval)
241 free r
242 free e
243 free neval
244 freeHaskellFunPtr fp
245 return result
246
247foreign import ccall safe "integrate_cquad" c_integrate_cquad
248 :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
249 -> Double -> Double -> CInt -> Ptr Double -> Ptr Double -> Ptr Int -> IO CInt
250
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 @@
1-- Module : Numeric.GSL.Internal
2-- Copyright : (c) Alberto Ruiz 2009
3-- License : GPL
4--
5-- Maintainer : Alberto Ruiz (aruiz at um dot es)
6-- Stability : provisional
7-- Portability : uses ffi
8--
9-- Auxiliary functions.
10--
11-- #hide
12
13module Numeric.GSL.Internal where
14
15import Data.Packed.Internal
16
17import Foreign.Marshal.Array(copyArray)
18import Foreign.Ptr(Ptr, FunPtr)
19import Foreign.C.Types
20import System.IO.Unsafe(unsafePerformIO)
21
22iv :: (Vector Double -> Double) -> (CInt -> Ptr Double -> Double)
23iv f n p = f (createV (fromIntegral n) copy "iv") where
24 copy n' q = do
25 copyArray q p (fromIntegral n')
26 return 0
27
28-- | conversion of Haskell functions into function pointers that can be used in the C side
29foreign import ccall safe "wrapper"
30 mkVecfun :: (CInt -> Ptr Double -> Double)
31 -> IO( FunPtr (CInt -> Ptr Double -> Double))
32
33foreign import ccall safe "wrapper"
34 mkVecVecfun :: TVV -> IO (FunPtr TVV)
35
36foreign import ccall safe "wrapper"
37 mkDoubleVecVecfun :: (Double -> TVV) -> IO (FunPtr (Double -> TVV))
38
39foreign import ccall safe "wrapper"
40 mkDoublefun :: (Double -> Double) -> IO (FunPtr (Double -> Double))
41
42aux_vTov :: (Vector Double -> Vector Double) -> TVV
43aux_vTov f n p nr r = g where
44 v = f x
45 x = createV (fromIntegral n) copy "aux_vTov"
46 copy n' q = do
47 copyArray q p (fromIntegral n')
48 return 0
49 g = do unsafeWith v $ \p' -> copyArray r p' (fromIntegral nr)
50 return 0
51
52foreign import ccall safe "wrapper"
53 mkVecMatfun :: TVM -> IO (FunPtr TVM)
54
55foreign import ccall safe "wrapper"
56 mkDoubleVecMatfun :: (Double -> TVM) -> IO (FunPtr (Double -> TVM))
57
58aux_vTom :: (Vector Double -> Matrix Double) -> TVM
59aux_vTom f n p rr cr r = g where
60 v = flatten $ f x
61 x = createV (fromIntegral n) copy "aux_vTov"
62 copy n' q = do
63 copyArray q p (fromIntegral n')
64 return 0
65 g = do unsafeWith v $ \p' -> copyArray r p' (fromIntegral $ rr*cr)
66 return 0
67
68createV n fun msg = unsafePerformIO $ do
69 r <- createVector n
70 app1 fun vec r msg
71 return r
72
73createMIO r c fun msg = do
74 res <- createMatrix RowMajor r c
75 app1 fun mat res msg
76 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 @@
1{-# LANGUAGE ForeignFunctionInterface #-}
2-----------------------------------------------------------------------------
3{- |
4Module : Numeric.GSL.Minimization
5Copyright : (c) Alberto Ruiz 2006-9
6License : GPL-style
7
8Maintainer : Alberto Ruiz (aruiz at um dot es)
9Stability : provisional
10Portability : uses ffi
11
12Minimization of a multidimensional function using some of the algorithms described in:
13
14<http://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Minimization.html>
15
16The example in the GSL manual:
17
18@
19f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30
20
21main = do
22 let (s,p) = minimize NMSimplex2 1E-2 30 [1,1] f [5,7]
23 print s
24 print p
25@
26
27>>> main
28[0.9920430849306288,1.9969168063253182]
29 0.000 512.500 1.130 6.500 5.000
30 1.000 290.625 1.409 5.250 4.000
31 2.000 290.625 1.409 5.250 4.000
32 3.000 252.500 1.409 5.500 1.000
33 ...
3422.000 30.001 0.013 0.992 1.997
3523.000 30.001 0.008 0.992 1.997
36
37The path to the solution can be graphically shown by means of:
38
39@'Graphics.Plot.mplot' $ drop 3 ('toColumns' p)@
40
41Taken from the GSL manual:
42
43The 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.
44
45The 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).
46
47The 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.
48
49-}
50
51-----------------------------------------------------------------------------
52module Numeric.GSL.Minimization (
53 minimize, minimizeV, MinimizeMethod(..),
54 minimizeD, minimizeVD, MinimizeMethodD(..),
55 uniMinimize, UniMinimizeMethod(..),
56
57 minimizeNMSimplex,
58 minimizeConjugateGradient,
59 minimizeVectorBFGS2
60) where
61
62
63import Data.Packed.Internal
64import Data.Packed.Matrix
65import Numeric.GSL.Internal
66
67import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr)
68import Foreign.C.Types
69import System.IO.Unsafe(unsafePerformIO)
70
71------------------------------------------------------------------------
72
73{-# DEPRECATED minimizeNMSimplex "use minimize NMSimplex2 eps maxit sizes f xi" #-}
74minimizeNMSimplex f xi szs eps maxit = minimize NMSimplex eps maxit szs f xi
75
76{-# DEPRECATED minimizeConjugateGradient "use minimizeD ConjugateFR eps maxit step tol f g xi" #-}
77minimizeConjugateGradient step tol eps maxit f g xi = minimizeD ConjugateFR eps maxit step tol f g xi
78
79{-# DEPRECATED minimizeVectorBFGS2 "use minimizeD VectorBFGS2 eps maxit step tol f g xi" #-}
80minimizeVectorBFGS2 step tol eps maxit f g xi = minimizeD VectorBFGS2 eps maxit step tol f g xi
81
82-------------------------------------------------------------------------
83
84data UniMinimizeMethod = GoldenSection
85 | BrentMini
86 | QuadGolden
87 deriving (Enum, Eq, Show, Bounded)
88
89-- | Onedimensional minimization.
90
91uniMinimize :: UniMinimizeMethod -- ^ The method used.
92 -> Double -- ^ desired precision of the solution
93 -> Int -- ^ maximum number of iterations allowed
94 -> (Double -> Double) -- ^ function to minimize
95 -> Double -- ^ guess for the location of the minimum
96 -> Double -- ^ lower bound of search interval
97 -> Double -- ^ upper bound of search interval
98 -> (Double, Matrix Double) -- ^ solution and optimization path
99
100uniMinimize method epsrel maxit fun xmin xl xu = uniMinimizeGen (fi (fromEnum method)) fun xmin xl xu epsrel maxit
101
102uniMinimizeGen m f xmin xl xu epsrel maxit = unsafePerformIO $ do
103 fp <- mkDoublefun f
104 rawpath <- createMIO maxit 4
105 (c_uniMinize m fp epsrel (fi maxit) xmin xl xu)
106 "uniMinimize"
107 let it = round (rawpath @@> (maxit-1,0))
108 path = takeRows it rawpath
109 [sol] = toLists $ dropRows (it-1) path
110 freeHaskellFunPtr fp
111 return (sol !! 1, path)
112
113
114foreign import ccall safe "uniMinimize"
115 c_uniMinize:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> Double -> TM
116
117data MinimizeMethod = NMSimplex
118 | NMSimplex2
119 deriving (Enum,Eq,Show,Bounded)
120
121-- | Minimization without derivatives
122minimize :: MinimizeMethod
123 -> Double -- ^ desired precision of the solution (size test)
124 -> Int -- ^ maximum number of iterations allowed
125 -> [Double] -- ^ sizes of the initial search box
126 -> ([Double] -> Double) -- ^ function to minimize
127 -> [Double] -- ^ starting point
128 -> ([Double], Matrix Double) -- ^ solution vector and optimization path
129
130-- | Minimization without derivatives (vector version)
131minimizeV :: MinimizeMethod
132 -> Double -- ^ desired precision of the solution (size test)
133 -> Int -- ^ maximum number of iterations allowed
134 -> Vector Double -- ^ sizes of the initial search box
135 -> (Vector Double -> Double) -- ^ function to minimize
136 -> Vector Double -- ^ starting point
137 -> (Vector Double, Matrix Double) -- ^ solution vector and optimization path
138
139minimize method eps maxit sz f xi = v2l $ minimizeV method eps maxit (fromList sz) (f.toList) (fromList xi)
140 where v2l (v,m) = (toList v, m)
141
142ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2
143
144minimizeV method eps maxit szv f xiv = unsafePerformIO $ do
145 let n = dim xiv
146 fp <- mkVecfun (iv f)
147 rawpath <- ww2 vec xiv vec szv $ \xiv' szv' ->
148 createMIO maxit (n+3)
149 (c_minimize (fi (fromEnum method)) fp eps (fi maxit) // xiv' // szv')
150 "minimize"
151 let it = round (rawpath @@> (maxit-1,0))
152 path = takeRows it rawpath
153 sol = cdat $ dropColumns 3 $ dropRows (it-1) path
154 freeHaskellFunPtr fp
155 return (sol, path)
156
157
158foreign import ccall safe "gsl-aux.h minimize"
159 c_minimize:: CInt -> FunPtr (CInt -> Ptr Double -> Double) -> Double -> CInt -> TVVM
160
161----------------------------------------------------------------------------------
162
163
164data MinimizeMethodD = ConjugateFR
165 | ConjugatePR
166 | VectorBFGS
167 | VectorBFGS2
168 | SteepestDescent
169 deriving (Enum,Eq,Show,Bounded)
170
171-- | Minimization with derivatives.
172minimizeD :: MinimizeMethodD
173 -> Double -- ^ desired precision of the solution (gradient test)
174 -> Int -- ^ maximum number of iterations allowed
175 -> Double -- ^ size of the first trial step
176 -> Double -- ^ tol (precise meaning depends on method)
177 -> ([Double] -> Double) -- ^ function to minimize
178 -> ([Double] -> [Double]) -- ^ gradient
179 -> [Double] -- ^ starting point
180 -> ([Double], Matrix Double) -- ^ solution vector and optimization path
181
182-- | Minimization with derivatives (vector version)
183minimizeVD :: MinimizeMethodD
184 -> Double -- ^ desired precision of the solution (gradient test)
185 -> Int -- ^ maximum number of iterations allowed
186 -> Double -- ^ size of the first trial step
187 -> Double -- ^ tol (precise meaning depends on method)
188 -> (Vector Double -> Double) -- ^ function to minimize
189 -> (Vector Double -> Vector Double) -- ^ gradient
190 -> Vector Double -- ^ starting point
191 -> (Vector Double, Matrix Double) -- ^ solution vector and optimization path
192
193minimizeD method eps maxit istep tol f df xi = v2l $ minimizeVD
194 method eps maxit istep tol (f.toList) (fromList.df.toList) (fromList xi)
195 where v2l (v,m) = (toList v, m)
196
197
198minimizeVD method eps maxit istep tol f df xiv = unsafePerformIO $ do
199 let n = dim xiv
200 f' = f
201 df' = (checkdim1 n . df)
202 fp <- mkVecfun (iv f')
203 dfp <- mkVecVecfun (aux_vTov df')
204 rawpath <- vec xiv $ \xiv' ->
205 createMIO maxit (n+2)
206 (c_minimizeD (fi (fromEnum method)) fp dfp istep tol eps (fi maxit) // xiv')
207 "minimizeD"
208 let it = round (rawpath @@> (maxit-1,0))
209 path = takeRows it rawpath
210 sol = cdat $ dropColumns 2 $ dropRows (it-1) path
211 freeHaskellFunPtr fp
212 freeHaskellFunPtr dfp
213 return (sol,path)
214
215foreign import ccall safe "gsl-aux.h minimizeD"
216 c_minimizeD :: CInt
217 -> FunPtr (CInt -> Ptr Double -> Double)
218 -> FunPtr TVV
219 -> Double -> Double -> Double -> CInt
220 -> TVM
221
222---------------------------------------------------------------------
223
224checkdim1 n v
225 | dim v == n = v
226 | otherwise = error $ "Error: "++ show n
227 ++ " 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 @@
1{- |
2Module : Numeric.GSL.ODE
3Copyright : (c) Alberto Ruiz 2010
4License : GPL
5
6Maintainer : Alberto Ruiz (aruiz at um dot es)
7Stability : provisional
8Portability : uses ffi
9
10Solution of ordinary differential equation (ODE) initial value problems.
11
12<http://www.gnu.org/software/gsl/manual/html_node/Ordinary-Differential-Equations.html>
13
14A simple example:
15
16@
17import Numeric.GSL.ODE
18import Numeric.LinearAlgebra
19import Numeric.LinearAlgebra.Util(mplot)
20
21xdot t [x,v] = [v, -0.95*x - 0.1*v]
22
23ts = linspace 100 (0,20 :: Double)
24
25sol = odeSolve xdot [10,0] ts
26
27main = mplot (ts : toColumns sol)
28@
29
30-}
31-----------------------------------------------------------------------------
32
33module Numeric.GSL.ODE (
34 odeSolve, odeSolveV, ODEMethod(..), Jacobian
35) where
36
37import Data.Packed.Internal
38import Numeric.GSL.Internal
39
40import Foreign.Ptr(FunPtr, nullFunPtr, freeHaskellFunPtr)
41import Foreign.C.Types
42import System.IO.Unsafe(unsafePerformIO)
43
44-------------------------------------------------------------------------
45
46type Jacobian = Double -> Vector Double -> Matrix Double
47
48-- | Stepping functions
49data ODEMethod = RK2 -- ^ Embedded Runge-Kutta (2, 3) method.
50 | 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.
51 | RKf45 -- ^ Embedded Runge-Kutta-Fehlberg (4, 5) method. This method is a good general-purpose integrator.
52 | RKck -- ^ Embedded Runge-Kutta Cash-Karp (4, 5) method.
53 | RK8pd -- ^ Embedded Runge-Kutta Prince-Dormand (8,9) method.
54 | RK2imp Jacobian -- ^ Implicit 2nd order Runge-Kutta at Gaussian points.
55 | RK4imp Jacobian -- ^ Implicit 4th order Runge-Kutta at Gaussian points.
56 | BSimp Jacobian -- ^ Implicit Bulirsch-Stoer method of Bader and Deuflhard. The method is generally suitable for stiff problems.
57 | 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.
58 | 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.
59 | 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.
60
61
62-- | A version of 'odeSolveV' with reasonable default parameters and system of equations defined using lists.
63odeSolve
64 :: (Double -> [Double] -> [Double]) -- ^ xdot(t,x)
65 -> [Double] -- ^ initial conditions
66 -> Vector Double -- ^ desired solution times
67 -> Matrix Double -- ^ solution
68odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) (fromList xi) ts
69 where hi = (ts@>1 - ts@>0)/100
70 epsAbs = 1.49012e-08
71 epsRel = 1.49012e-08
72 l2v f = \t -> fromList . f t . toList
73
74-- | Evolution of the system with adaptive step-size control.
75odeSolveV
76 :: ODEMethod
77 -> Double -- ^ initial step size
78 -> Double -- ^ absolute tolerance for the state vector
79 -> Double -- ^ relative tolerance for the state vector
80 -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x)
81 -> Vector Double -- ^ initial conditions
82 -> Vector Double -- ^ desired solution times
83 -> Matrix Double -- ^ solution
84odeSolveV RK2 = odeSolveV' 0 Nothing
85odeSolveV RK4 = odeSolveV' 1 Nothing
86odeSolveV RKf45 = odeSolveV' 2 Nothing
87odeSolveV RKck = odeSolveV' 3 Nothing
88odeSolveV RK8pd = odeSolveV' 4 Nothing
89odeSolveV (RK2imp jac) = odeSolveV' 5 (Just jac)
90odeSolveV (RK4imp jac) = odeSolveV' 6 (Just jac)
91odeSolveV (BSimp jac) = odeSolveV' 7 (Just jac)
92odeSolveV (RK1imp jac) = odeSolveV' 8 (Just jac)
93odeSolveV MSAdams = odeSolveV' 9 Nothing
94odeSolveV (MSBDF jac) = odeSolveV' 10 (Just jac)
95
96
97odeSolveV'
98 :: CInt
99 -> Maybe (Double -> Vector Double -> Matrix Double) -- ^ optional jacobian
100 -> Double -- ^ initial step size
101 -> Double -- ^ absolute tolerance for the state vector
102 -> Double -- ^ relative tolerance for the state vector
103 -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x)
104 -> Vector Double -- ^ initial conditions
105 -> Vector Double -- ^ desired solution times
106 -> Matrix Double -- ^ solution
107odeSolveV' method mbjac h epsAbs epsRel f xiv ts = unsafePerformIO $ do
108 let n = dim xiv
109 fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t))
110 jp <- case mbjac of
111 Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t))
112 Nothing -> return nullFunPtr
113 sol <- vec xiv $ \xiv' ->
114 vec (checkTimes ts) $ \ts' ->
115 createMIO (dim ts) n
116 (ode_c (method) h epsAbs epsRel fp jp // xiv' // ts' )
117 "ode"
118 freeHaskellFunPtr fp
119 return sol
120
121foreign import ccall safe "ode"
122 ode_c :: CInt -> Double -> Double -> Double -> FunPtr (Double -> TVV) -> FunPtr (Double -> TVM) -> TVVM
123
124-------------------------------------------------------
125
126checkdim1 n v
127 | dim v == n = v
128 | otherwise = error $ "Error: "++ show n
129 ++ " components expected in the result of the function supplied to odeSolve"
130
131checkdim2 n m
132 | rows m == n && cols m == n = m
133 | otherwise = error $ "Error: "++ show n ++ "x" ++ show n
134 ++ " Jacobian expected in odeSolve"
135
136checkTimes ts | dim ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts
137 | otherwise = error "odeSolve requires increasing times"
138 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 @@
1{-# LANGUAGE CPP, ForeignFunctionInterface #-}
2-----------------------------------------------------------------------------
3{- |
4Module : Numeric.GSL.Polynomials
5Copyright : (c) Alberto Ruiz 2006
6License : GPL-style
7
8Maintainer : Alberto Ruiz (aruiz at um dot es)
9Stability : provisional
10Portability : uses ffi
11
12Polynomials.
13
14<http://www.gnu.org/software/gsl/manual/html_node/General-Polynomial-Equations.html#General-Polynomial-Equations>
15
16-}
17-----------------------------------------------------------------------------
18module Numeric.GSL.Polynomials (
19 polySolve
20) where
21
22import Data.Packed.Internal
23import Data.Complex
24import System.IO.Unsafe (unsafePerformIO)
25
26#if __GLASGOW_HASKELL__ >= 704
27import Foreign.C.Types (CInt(..))
28#endif
29
30{- | Solution of general polynomial equations, using /gsl_poly_complex_solve/.
31
32For example, the three solutions of x^3 + 8 = 0
33
34>>> polySolve [8,0,0,1]
35[(-2.0) :+ 0.0,1.0 :+ 1.7320508075688776,1.0 :+ (-1.7320508075688776)]
36
37
38The example in the GSL manual: To find the roots of x^5 -1 = 0:
39
40>>> polySolve [-1, 0, 0, 0, 0, 1]
41[(-0.8090169943749472) :+ 0.5877852522924731,
42(-0.8090169943749472) :+ (-0.5877852522924731),
430.30901699437494756 :+ 0.9510565162951535,
440.30901699437494756 :+ (-0.9510565162951535),
451.0000000000000002 :+ 0.0]
46
47-}
48polySolve :: [Double] -> [Complex Double]
49polySolve = toList . polySolve' . fromList
50
51polySolve' :: Vector Double -> Vector (Complex Double)
52polySolve' v | dim v > 1 = unsafePerformIO $ do
53 r <- createVector (dim v-1)
54 app2 c_polySolve vec v vec r "polySolve"
55 return r
56 | otherwise = error "polySolve on a polynomial of degree zero"
57
58foreign 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 @@
1{- |
2Module : Numeric.GSL.Root
3Copyright : (c) Alberto Ruiz 2009
4License : GPL
5
6Maintainer : Alberto Ruiz (aruiz at um dot es)
7Stability : provisional
8Portability : uses ffi
9
10Multidimensional root finding.
11
12<http://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Root_002dFinding.html>
13
14The example in the GSL manual:
15
16>>> let rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ]
17>>> let (sol,path) = root Hybrids 1E-7 30 (rosenbrock 1 10) [-10,-5]
18>>> sol
19[1.0,1.0]
20>>> disp 3 path
2111x5
22 1.000 -10.000 -5.000 11.000 -1050.000
23 2.000 -3.976 24.827 4.976 90.203
24 3.000 -3.976 24.827 4.976 90.203
25 4.000 -3.976 24.827 4.976 90.203
26 5.000 -1.274 -5.680 2.274 -73.018
27 6.000 -1.274 -5.680 2.274 -73.018
28 7.000 0.249 0.298 0.751 2.359
29 8.000 0.249 0.298 0.751 2.359
30 9.000 1.000 0.878 -0.000 -1.218
3110.000 1.000 0.989 -0.000 -0.108
3211.000 1.000 1.000 0.000 0.000
33
34-}
35-----------------------------------------------------------------------------
36
37module Numeric.GSL.Root (
38 uniRoot, UniRootMethod(..),
39 uniRootJ, UniRootMethodJ(..),
40 root, RootMethod(..),
41 rootJ, RootMethodJ(..),
42) where
43
44import Data.Packed.Internal
45import Data.Packed.Matrix
46import Numeric.GSL.Internal
47import Foreign.Ptr(FunPtr, freeHaskellFunPtr)
48import Foreign.C.Types
49import System.IO.Unsafe(unsafePerformIO)
50
51-------------------------------------------------------------------------
52
53data UniRootMethod = Bisection
54 | FalsePos
55 | Brent
56 deriving (Enum, Eq, Show, Bounded)
57
58uniRoot :: UniRootMethod
59 -> Double
60 -> Int
61 -> (Double -> Double)
62 -> Double
63 -> Double
64 -> (Double, Matrix Double)
65uniRoot method epsrel maxit fun xl xu = uniRootGen (fi (fromEnum method)) fun xl xu epsrel maxit
66
67uniRootGen m f xl xu epsrel maxit = unsafePerformIO $ do
68 fp <- mkDoublefun f
69 rawpath <- createMIO maxit 4
70 (c_root m fp epsrel (fi maxit) xl xu)
71 "root"
72 let it = round (rawpath @@> (maxit-1,0))
73 path = takeRows it rawpath
74 [sol] = toLists $ dropRows (it-1) path
75 freeHaskellFunPtr fp
76 return (sol !! 1, path)
77
78foreign import ccall safe "root"
79 c_root:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> TM
80
81-------------------------------------------------------------------------
82data UniRootMethodJ = UNewton
83 | Secant
84 | Steffenson
85 deriving (Enum, Eq, Show, Bounded)
86
87uniRootJ :: UniRootMethodJ
88 -> Double
89 -> Int
90 -> (Double -> Double)
91 -> (Double -> Double)
92 -> Double
93 -> (Double, Matrix Double)
94uniRootJ method epsrel maxit fun dfun x = uniRootJGen (fi (fromEnum method)) fun
95 dfun x epsrel maxit
96
97uniRootJGen m f df x epsrel maxit = unsafePerformIO $ do
98 fp <- mkDoublefun f
99 dfp <- mkDoublefun df
100 rawpath <- createMIO maxit 2
101 (c_rootj m fp dfp epsrel (fi maxit) x)
102 "rootj"
103 let it = round (rawpath @@> (maxit-1,0))
104 path = takeRows it rawpath
105 [sol] = toLists $ dropRows (it-1) path
106 freeHaskellFunPtr fp
107 return (sol !! 1, path)
108
109foreign import ccall safe "rootj"
110 c_rootj :: CInt -> FunPtr (Double -> Double) -> FunPtr (Double -> Double)
111 -> Double -> CInt -> Double -> TM
112
113-------------------------------------------------------------------------
114
115data RootMethod = Hybrids
116 | Hybrid
117 | DNewton
118 | Broyden
119 deriving (Enum,Eq,Show,Bounded)
120
121-- | Nonlinear multidimensional root finding using algorithms that do not require
122-- any derivative information to be supplied by the user.
123-- Any derivatives needed are approximated by finite differences.
124root :: RootMethod
125 -> Double -- ^ maximum residual
126 -> Int -- ^ maximum number of iterations allowed
127 -> ([Double] -> [Double]) -- ^ function to minimize
128 -> [Double] -- ^ starting point
129 -> ([Double], Matrix Double) -- ^ solution vector and optimization path
130
131root method epsabs maxit fun xinit = rootGen (fi (fromEnum method)) fun xinit epsabs maxit
132
133rootGen m f xi epsabs maxit = unsafePerformIO $ do
134 let xiv = fromList xi
135 n = dim xiv
136 fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
137 rawpath <- vec xiv $ \xiv' ->
138 createMIO maxit (2*n+1)
139 (c_multiroot m fp epsabs (fi maxit) // xiv')
140 "multiroot"
141 let it = round (rawpath @@> (maxit-1,0))
142 path = takeRows it rawpath
143 [sol] = toLists $ dropRows (it-1) path
144 freeHaskellFunPtr fp
145 return (take n $ drop 1 sol, path)
146
147
148foreign import ccall safe "multiroot"
149 c_multiroot:: CInt -> FunPtr TVV -> Double -> CInt -> TVM
150
151-------------------------------------------------------------------------
152
153data RootMethodJ = HybridsJ
154 | HybridJ
155 | Newton
156 | GNewton
157 deriving (Enum,Eq,Show,Bounded)
158
159-- | Nonlinear multidimensional root finding using both the function and its derivatives.
160rootJ :: RootMethodJ
161 -> Double -- ^ maximum residual
162 -> Int -- ^ maximum number of iterations allowed
163 -> ([Double] -> [Double]) -- ^ function to minimize
164 -> ([Double] -> [[Double]]) -- ^ Jacobian
165 -> [Double] -- ^ starting point
166 -> ([Double], Matrix Double) -- ^ solution vector and optimization path
167
168rootJ method epsabs maxit fun jac xinit = rootJGen (fi (fromEnum method)) fun jac xinit epsabs maxit
169
170rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do
171 let xiv = fromList xi
172 n = dim xiv
173 fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
174 jp <- mkVecMatfun (aux_vTom (checkdim2 n . fromLists . jac . toList))
175 rawpath <- vec xiv $ \xiv' ->
176 createMIO maxit (2*n+1)
177 (c_multirootj m fp jp epsabs (fi maxit) // xiv')
178 "multiroot"
179 let it = round (rawpath @@> (maxit-1,0))
180 path = takeRows it rawpath
181 [sol] = toLists $ dropRows (it-1) path
182 freeHaskellFunPtr fp
183 freeHaskellFunPtr jp
184 return (take n $ drop 1 sol, path)
185
186foreign import ccall safe "multirootj"
187 c_multirootj:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> CInt -> TVM
188
189-------------------------------------------------------
190
191checkdim1 n v
192 | dim v == n = v
193 | otherwise = error $ "Error: "++ show n
194 ++ " components expected in the result of the function supplied to root"
195
196checkdim2 n m
197 | rows m == n && cols m == n = m
198 | otherwise = error $ "Error: "++ show n ++ "x" ++ show n
199 ++ " 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 @@
1-----------------------------------------------------------------------------
2-- |
3-- Module : Numeric.GSL.Vector
4-- Copyright : (c) Alberto Ruiz 2007
5-- License : GPL-style
6--
7-- Maintainer : Alberto Ruiz <aruiz@um.es>
8-- Stability : provisional
9-- Portability : portable (uses FFI)
10--
11-- Low level interface to vector operations.
12--
13-----------------------------------------------------------------------------
14
15module Numeric.GSL.Vector (
16 sumF, sumR, sumQ, sumC,
17 prodF, prodR, prodQ, prodC,
18 dotF, dotR, dotQ, dotC,
19 FunCodeS(..), toScalarR, toScalarF, toScalarC, toScalarQ,
20 FunCodeV(..), vectorMapR, vectorMapC, vectorMapF, vectorMapQ,
21 FunCodeSV(..), vectorMapValR, vectorMapValC, vectorMapValF, vectorMapValQ,
22 FunCodeVV(..), vectorZipR, vectorZipC, vectorZipF, vectorZipQ,
23 RandDist(..), randomVector
24) where
25
26import Data.Packed.Internal.Common
27import Data.Packed.Internal.Signatures
28import Data.Packed.Internal.Vector
29
30import Data.Complex
31import Foreign.Marshal.Alloc(free)
32import Foreign.Marshal.Array(newArray)
33import Foreign.Ptr(Ptr)
34import Foreign.C.Types
35import System.IO.Unsafe(unsafePerformIO)
36import Control.Monad(when)
37
38fromei x = fromIntegral (fromEnum x) :: CInt
39
40data FunCodeV = Sin
41 | Cos
42 | Tan
43 | Abs
44 | ASin
45 | ACos
46 | ATan
47 | Sinh
48 | Cosh
49 | Tanh
50 | ASinh
51 | ACosh
52 | ATanh
53 | Exp
54 | Log
55 | Sign
56 | Sqrt
57 deriving Enum
58
59data FunCodeSV = Scale
60 | Recip
61 | AddConstant
62 | Negate
63 | PowSV
64 | PowVS
65 deriving Enum
66
67data FunCodeVV = Add
68 | Sub
69 | Mul
70 | Div
71 | Pow
72 | ATan2
73 deriving Enum
74
75data FunCodeS = Norm2
76 | AbsSum
77 | MaxIdx
78 | Max
79 | MinIdx
80 | Min
81 deriving Enum
82
83------------------------------------------------------------------
84
85-- | sum of elements
86sumF :: Vector Float -> Float
87sumF x = unsafePerformIO $ do
88 r <- createVector 1
89 app2 c_sumF vec x vec r "sumF"
90 return $ r @> 0
91
92-- | sum of elements
93sumR :: Vector Double -> Double
94sumR x = unsafePerformIO $ do
95 r <- createVector 1
96 app2 c_sumR vec x vec r "sumR"
97 return $ r @> 0
98
99-- | sum of elements
100sumQ :: Vector (Complex Float) -> Complex Float
101sumQ x = unsafePerformIO $ do
102 r <- createVector 1
103 app2 c_sumQ vec x vec r "sumQ"
104 return $ r @> 0
105
106-- | sum of elements
107sumC :: Vector (Complex Double) -> Complex Double
108sumC x = unsafePerformIO $ do
109 r <- createVector 1
110 app2 c_sumC vec x vec r "sumC"
111 return $ r @> 0
112
113foreign import ccall unsafe "gsl-aux.h sumF" c_sumF :: TFF
114foreign import ccall unsafe "gsl-aux.h sumR" c_sumR :: TVV
115foreign import ccall unsafe "gsl-aux.h sumQ" c_sumQ :: TQVQV
116foreign import ccall unsafe "gsl-aux.h sumC" c_sumC :: TCVCV
117
118-- | product of elements
119prodF :: Vector Float -> Float
120prodF x = unsafePerformIO $ do
121 r <- createVector 1
122 app2 c_prodF vec x vec r "prodF"
123 return $ r @> 0
124
125-- | product of elements
126prodR :: Vector Double -> Double
127prodR x = unsafePerformIO $ do
128 r <- createVector 1
129 app2 c_prodR vec x vec r "prodR"
130 return $ r @> 0
131
132-- | product of elements
133prodQ :: Vector (Complex Float) -> Complex Float
134prodQ x = unsafePerformIO $ do
135 r <- createVector 1
136 app2 c_prodQ vec x vec r "prodQ"
137 return $ r @> 0
138
139-- | product of elements
140prodC :: Vector (Complex Double) -> Complex Double
141prodC x = unsafePerformIO $ do
142 r <- createVector 1
143 app2 c_prodC vec x vec r "prodC"
144 return $ r @> 0
145
146foreign import ccall unsafe "gsl-aux.h prodF" c_prodF :: TFF
147foreign import ccall unsafe "gsl-aux.h prodR" c_prodR :: TVV
148foreign import ccall unsafe "gsl-aux.h prodQ" c_prodQ :: TQVQV
149foreign import ccall unsafe "gsl-aux.h prodC" c_prodC :: TCVCV
150
151-- | dot product
152dotF :: Vector Float -> Vector Float -> Float
153dotF x y = unsafePerformIO $ do
154 r <- createVector 1
155 app3 c_dotF vec x vec y vec r "dotF"
156 return $ r @> 0
157
158-- | dot product
159dotR :: Vector Double -> Vector Double -> Double
160dotR x y = unsafePerformIO $ do
161 r <- createVector 1
162 app3 c_dotR vec x vec y vec r "dotR"
163 return $ r @> 0
164
165-- | dot product
166dotQ :: Vector (Complex Float) -> Vector (Complex Float) -> Complex Float
167dotQ x y = unsafePerformIO $ do
168 r <- createVector 1
169 app3 c_dotQ vec x vec y vec r "dotQ"
170 return $ r @> 0
171
172-- | dot product
173dotC :: Vector (Complex Double) -> Vector (Complex Double) -> Complex Double
174dotC x y = unsafePerformIO $ do
175 r <- createVector 1
176 app3 c_dotC vec x vec y vec r "dotC"
177 return $ r @> 0
178
179foreign import ccall unsafe "gsl-aux.h dotF" c_dotF :: TFFF
180foreign import ccall unsafe "gsl-aux.h dotR" c_dotR :: TVVV
181foreign import ccall unsafe "gsl-aux.h dotQ" c_dotQ :: TQVQVQV
182foreign import ccall unsafe "gsl-aux.h dotC" c_dotC :: TCVCVCV
183
184------------------------------------------------------------------
185
186toScalarAux fun code v = unsafePerformIO $ do
187 r <- createVector 1
188 app2 (fun (fromei code)) vec v vec r "toScalarAux"
189 return (r `at` 0)
190
191vectorMapAux fun code v = unsafePerformIO $ do
192 r <- createVector (dim v)
193 app2 (fun (fromei code)) vec v vec r "vectorMapAux"
194 return r
195
196vectorMapValAux fun code val v = unsafePerformIO $ do
197 r <- createVector (dim v)
198 pval <- newArray [val]
199 app2 (fun (fromei code) pval) vec v vec r "vectorMapValAux"
200 free pval
201 return r
202
203vectorZipAux fun code u v = unsafePerformIO $ do
204 r <- createVector (dim u)
205 when (dim u > 0) $ app3 (fun (fromei code)) vec u vec v vec r "vectorZipAux"
206 return r
207
208---------------------------------------------------------------------
209
210-- | obtains different functions of a vector: norm1, norm2, max, min, posmax, posmin, etc.
211toScalarR :: FunCodeS -> Vector Double -> Double
212toScalarR oper = toScalarAux c_toScalarR (fromei oper)
213
214foreign import ccall unsafe "gsl-aux.h toScalarR" c_toScalarR :: CInt -> TVV
215
216-- | obtains different functions of a vector: norm1, norm2, max, min, posmax, posmin, etc.
217toScalarF :: FunCodeS -> Vector Float -> Float
218toScalarF oper = toScalarAux c_toScalarF (fromei oper)
219
220foreign import ccall unsafe "gsl-aux.h toScalarF" c_toScalarF :: CInt -> TFF
221
222-- | obtains different functions of a vector: only norm1, norm2
223toScalarC :: FunCodeS -> Vector (Complex Double) -> Double
224toScalarC oper = toScalarAux c_toScalarC (fromei oper)
225
226foreign import ccall unsafe "gsl-aux.h toScalarC" c_toScalarC :: CInt -> TCVV
227
228-- | obtains different functions of a vector: only norm1, norm2
229toScalarQ :: FunCodeS -> Vector (Complex Float) -> Float
230toScalarQ oper = toScalarAux c_toScalarQ (fromei oper)
231
232foreign import ccall unsafe "gsl-aux.h toScalarQ" c_toScalarQ :: CInt -> TQVF
233
234------------------------------------------------------------------
235
236-- | map of real vectors with given function
237vectorMapR :: FunCodeV -> Vector Double -> Vector Double
238vectorMapR = vectorMapAux c_vectorMapR
239
240foreign import ccall unsafe "gsl-aux.h mapR" c_vectorMapR :: CInt -> TVV
241
242-- | map of complex vectors with given function
243vectorMapC :: FunCodeV -> Vector (Complex Double) -> Vector (Complex Double)
244vectorMapC oper = vectorMapAux c_vectorMapC (fromei oper)
245
246foreign import ccall unsafe "gsl-aux.h mapC" c_vectorMapC :: CInt -> TCVCV
247
248-- | map of real vectors with given function
249vectorMapF :: FunCodeV -> Vector Float -> Vector Float
250vectorMapF = vectorMapAux c_vectorMapF
251
252foreign import ccall unsafe "gsl-aux.h mapF" c_vectorMapF :: CInt -> TFF
253
254-- | map of real vectors with given function
255vectorMapQ :: FunCodeV -> Vector (Complex Float) -> Vector (Complex Float)
256vectorMapQ = vectorMapAux c_vectorMapQ
257
258foreign import ccall unsafe "gsl-aux.h mapQ" c_vectorMapQ :: CInt -> TQVQV
259
260-------------------------------------------------------------------
261
262-- | map of real vectors with given function
263vectorMapValR :: FunCodeSV -> Double -> Vector Double -> Vector Double
264vectorMapValR oper = vectorMapValAux c_vectorMapValR (fromei oper)
265
266foreign import ccall unsafe "gsl-aux.h mapValR" c_vectorMapValR :: CInt -> Ptr Double -> TVV
267
268-- | map of complex vectors with given function
269vectorMapValC :: FunCodeSV -> Complex Double -> Vector (Complex Double) -> Vector (Complex Double)
270vectorMapValC = vectorMapValAux c_vectorMapValC
271
272foreign import ccall unsafe "gsl-aux.h mapValC" c_vectorMapValC :: CInt -> Ptr (Complex Double) -> TCVCV
273
274-- | map of real vectors with given function
275vectorMapValF :: FunCodeSV -> Float -> Vector Float -> Vector Float
276vectorMapValF oper = vectorMapValAux c_vectorMapValF (fromei oper)
277
278foreign import ccall unsafe "gsl-aux.h mapValF" c_vectorMapValF :: CInt -> Ptr Float -> TFF
279
280-- | map of complex vectors with given function
281vectorMapValQ :: FunCodeSV -> Complex Float -> Vector (Complex Float) -> Vector (Complex Float)
282vectorMapValQ oper = vectorMapValAux c_vectorMapValQ (fromei oper)
283
284foreign import ccall unsafe "gsl-aux.h mapValQ" c_vectorMapValQ :: CInt -> Ptr (Complex Float) -> TQVQV
285
286-------------------------------------------------------------------
287
288-- | elementwise operation on real vectors
289vectorZipR :: FunCodeVV -> Vector Double -> Vector Double -> Vector Double
290vectorZipR = vectorZipAux c_vectorZipR
291
292foreign import ccall unsafe "gsl-aux.h zipR" c_vectorZipR :: CInt -> TVVV
293
294-- | elementwise operation on complex vectors
295vectorZipC :: FunCodeVV -> Vector (Complex Double) -> Vector (Complex Double) -> Vector (Complex Double)
296vectorZipC = vectorZipAux c_vectorZipC
297
298foreign import ccall unsafe "gsl-aux.h zipC" c_vectorZipC :: CInt -> TCVCVCV
299
300-- | elementwise operation on real vectors
301vectorZipF :: FunCodeVV -> Vector Float -> Vector Float -> Vector Float
302vectorZipF = vectorZipAux c_vectorZipF
303
304foreign import ccall unsafe "gsl-aux.h zipF" c_vectorZipF :: CInt -> TFFF
305
306-- | elementwise operation on complex vectors
307vectorZipQ :: FunCodeVV -> Vector (Complex Float) -> Vector (Complex Float) -> Vector (Complex Float)
308vectorZipQ = vectorZipAux c_vectorZipQ
309
310foreign import ccall unsafe "gsl-aux.h zipQ" c_vectorZipQ :: CInt -> TQVQVQV
311
312-----------------------------------------------------------------------
313
314data RandDist = Uniform -- ^ uniform distribution in [0,1)
315 | Gaussian -- ^ normal distribution with mean zero and standard deviation one
316 deriving Enum
317
318-- | Obtains a vector of pseudorandom elements from the the mt19937 generator in GSL, with a given seed. Use randomIO to get a random seed.
319randomVector :: Int -- ^ seed
320 -> RandDist -- ^ distribution
321 -> Int -- ^ vector size
322 -> Vector Double
323randomVector seed dist n = unsafePerformIO $ do
324 r <- createVector n
325 app1 (c_random_vector (fi seed) ((fi.fromEnum) dist)) vec r "randomVector"
326 return r
327
328foreign 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 @@
1#include <gsl/gsl_complex.h>
2
3#define RVEC(A) int A##n, double*A##p
4#define RMAT(A) int A##r, int A##c, double* A##p
5#define KRVEC(A) int A##n, const double*A##p
6#define KRMAT(A) int A##r, int A##c, const double* A##p
7
8#define CVEC(A) int A##n, gsl_complex*A##p
9#define CMAT(A) int A##r, int A##c, gsl_complex* A##p
10#define KCVEC(A) int A##n, const gsl_complex*A##p
11#define KCMAT(A) int A##r, int A##c, const gsl_complex* A##p
12
13#define FVEC(A) int A##n, float*A##p
14#define FMAT(A) int A##r, int A##c, float* A##p
15#define KFVEC(A) int A##n, const float*A##p
16#define KFMAT(A) int A##r, int A##c, const float* A##p
17
18#define QVEC(A) int A##n, gsl_complex_float*A##p
19#define QMAT(A) int A##r, int A##c, gsl_complex_float* A##p
20#define KQVEC(A) int A##n, const gsl_complex_float*A##p
21#define KQMAT(A) int A##r, int A##c, const gsl_complex_float* A##p
22
23#include <gsl/gsl_blas.h>
24#include <gsl/gsl_math.h>
25#include <gsl/gsl_errno.h>
26#include <gsl/gsl_fft_complex.h>
27#include <gsl/gsl_integration.h>
28#include <gsl/gsl_deriv.h>
29#include <gsl/gsl_poly.h>
30#include <gsl/gsl_multimin.h>
31#include <gsl/gsl_multiroots.h>
32#include <gsl/gsl_min.h>
33#include <gsl/gsl_complex_math.h>
34#include <gsl/gsl_rng.h>
35#include <gsl/gsl_randist.h>
36#include <gsl/gsl_roots.h>
37#include <gsl/gsl_multifit_nlin.h>
38#include <string.h>
39#include <stdio.h>
40
41#define MACRO(B) do {B} while (0)
42#define ERROR(CODE) MACRO(return CODE;)
43#define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);})
44#define OK return 0;
45
46#define MIN(A,B) ((A)<(B)?(A):(B))
47#define MAX(A,B) ((A)>(B)?(A):(B))
48
49#ifdef DBG
50#define DEBUGMSG(M) printf("*** calling aux C function: %s\n",M);
51#else
52#define DEBUGMSG(M)
53#endif
54
55#define CHECK(RES,CODE) MACRO(if(RES) return CODE;)
56
57#ifdef DBG
58#define DEBUGMAT(MSG,X) printf(MSG" = \n"); gsl_matrix_fprintf(stdout,X,"%f"); printf("\n");
59#else
60#define DEBUGMAT(MSG,X)
61#endif
62
63#ifdef DBG
64#define DEBUGVEC(MSG,X) printf(MSG" = \n"); gsl_vector_fprintf(stdout,X,"%f"); printf("\n");
65#else
66#define DEBUGVEC(MSG,X)
67#endif
68
69#define DVVIEW(A) gsl_vector_view A = gsl_vector_view_array(A##p,A##n)
70#define DMVIEW(A) gsl_matrix_view A = gsl_matrix_view_array(A##p,A##r,A##c)
71#define CVVIEW(A) gsl_vector_complex_view A = gsl_vector_complex_view_array((double*)A##p,A##n)
72#define CMVIEW(A) gsl_matrix_complex_view A = gsl_matrix_complex_view_array((double*)A##p,A##r,A##c)
73#define KDVVIEW(A) gsl_vector_const_view A = gsl_vector_const_view_array(A##p,A##n)
74#define KDMVIEW(A) gsl_matrix_const_view A = gsl_matrix_const_view_array(A##p,A##r,A##c)
75#define KCVVIEW(A) gsl_vector_complex_const_view A = gsl_vector_complex_const_view_array((double*)A##p,A##n)
76#define KCMVIEW(A) gsl_matrix_complex_const_view A = gsl_matrix_complex_const_view_array((double*)A##p,A##r,A##c)
77
78#define FVVIEW(A) gsl_vector_float_view A = gsl_vector_float_view_array(A##p,A##n)
79#define FMVIEW(A) gsl_matrix_float_view A = gsl_matrix_float_view_array(A##p,A##r,A##c)
80#define QVVIEW(A) gsl_vector_complex_float_view A = gsl_vector_float_complex_view_array((float*)A##p,A##n)
81#define QMVIEW(A) gsl_matrix_complex_float_view A = gsl_matrix_float_complex_view_array((float*)A##p,A##r,A##c)
82#define KFVVIEW(A) gsl_vector_float_const_view A = gsl_vector_float_const_view_array(A##p,A##n)
83#define KFMVIEW(A) gsl_matrix_float_const_view A = gsl_matrix_float_const_view_array(A##p,A##r,A##c)
84#define KQVVIEW(A) gsl_vector_complex_float_const_view A = gsl_vector_complex_float_const_view_array((float*)A##p,A##n)
85#define KQMVIEW(A) gsl_matrix_complex_float_const_view A = gsl_matrix_complex_float_const_view_array((float*)A##p,A##r,A##c)
86
87#define V(a) (&a.vector)
88#define M(a) (&a.matrix)
89
90#define GCVEC(A) int A##n, gsl_complex*A##p
91#define KGCVEC(A) int A##n, const gsl_complex*A##p
92
93#define GQVEC(A) int A##n, gsl_complex_float*A##p
94#define KGQVEC(A) int A##n, const gsl_complex_float*A##p
95
96#define BAD_SIZE 2000
97#define BAD_CODE 2001
98#define MEM 2002
99#define BAD_FILE 2003
100
101
102void no_abort_on_error() {
103 gsl_set_error_handler_off();
104}
105
106
107int sumF(KFVEC(x),FVEC(r)) {
108 DEBUGMSG("sumF");
109 REQUIRES(rn==1,BAD_SIZE);
110 int i;
111 float res = 0;
112 for (i = 0; i < xn; i++) res += xp[i];
113 rp[0] = res;
114 OK
115}
116
117int sumR(KRVEC(x),RVEC(r)) {
118 DEBUGMSG("sumR");
119 REQUIRES(rn==1,BAD_SIZE);
120 int i;
121 double res = 0;
122 for (i = 0; i < xn; i++) res += xp[i];
123 rp[0] = res;
124 OK
125}
126
127int sumQ(KQVEC(x),QVEC(r)) {
128 DEBUGMSG("sumQ");
129 REQUIRES(rn==1,BAD_SIZE);
130 int i;
131 gsl_complex_float res;
132 res.dat[0] = 0;
133 res.dat[1] = 0;
134 for (i = 0; i < xn; i++) {
135 res.dat[0] += xp[i].dat[0];
136 res.dat[1] += xp[i].dat[1];
137 }
138 rp[0] = res;
139 OK
140}
141
142int sumC(KCVEC(x),CVEC(r)) {
143 DEBUGMSG("sumC");
144 REQUIRES(rn==1,BAD_SIZE);
145 int i;
146 gsl_complex res;
147 res.dat[0] = 0;
148 res.dat[1] = 0;
149 for (i = 0; i < xn; i++) {
150 res.dat[0] += xp[i].dat[0];
151 res.dat[1] += xp[i].dat[1];
152 }
153 rp[0] = res;
154 OK
155}
156
157int prodF(KFVEC(x),FVEC(r)) {
158 DEBUGMSG("prodF");
159 REQUIRES(rn==1,BAD_SIZE);
160 int i;
161 float res = 1;
162 for (i = 0; i < xn; i++) res *= xp[i];
163 rp[0] = res;
164 OK
165}
166
167int prodR(KRVEC(x),RVEC(r)) {
168 DEBUGMSG("prodR");
169 REQUIRES(rn==1,BAD_SIZE);
170 int i;
171 double res = 1;
172 for (i = 0; i < xn; i++) res *= xp[i];
173 rp[0] = res;
174 OK
175}
176
177int prodQ(KQVEC(x),QVEC(r)) {
178 DEBUGMSG("prodQ");
179 REQUIRES(rn==1,BAD_SIZE);
180 int i;
181 gsl_complex_float res;
182 float temp;
183 res.dat[0] = 1;
184 res.dat[1] = 0;
185 for (i = 0; i < xn; i++) {
186 temp = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1];
187 res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0];
188 res.dat[0] = temp;
189 }
190 rp[0] = res;
191 OK
192}
193
194int prodC(KCVEC(x),CVEC(r)) {
195 DEBUGMSG("prodC");
196 REQUIRES(rn==1,BAD_SIZE);
197 int i;
198 gsl_complex res;
199 double temp;
200 res.dat[0] = 1;
201 res.dat[1] = 0;
202 for (i = 0; i < xn; i++) {
203 temp = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1];
204 res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0];
205 res.dat[0] = temp;
206 }
207 rp[0] = res;
208 OK
209}
210
211int dotF(KFVEC(x), KFVEC(y), FVEC(r)) {
212 DEBUGMSG("dotF");
213 REQUIRES(xn==yn,BAD_SIZE);
214 REQUIRES(rn==1,BAD_SIZE);
215 DEBUGMSG("dotF");
216 KFVVIEW(x);
217 KFVVIEW(y);
218 gsl_blas_sdot(V(x),V(y),rp);
219 OK
220}
221
222int dotR(KRVEC(x), KRVEC(y), RVEC(r)) {
223 DEBUGMSG("dotR");
224 REQUIRES(xn==yn,BAD_SIZE);
225 REQUIRES(rn==1,BAD_SIZE);
226 DEBUGMSG("dotR");
227 KDVVIEW(x);
228 KDVVIEW(y);
229 gsl_blas_ddot(V(x),V(y),rp);
230 OK
231}
232
233int dotQ(KQVEC(x), KQVEC(y), QVEC(r)) {
234 DEBUGMSG("dotQ");
235 REQUIRES(xn==yn,BAD_SIZE);
236 REQUIRES(rn==1,BAD_SIZE);
237 DEBUGMSG("dotQ");
238 KQVVIEW(x);
239 KQVVIEW(y);
240 gsl_blas_cdotu(V(x),V(y),rp);
241 OK
242}
243
244int dotC(KCVEC(x), KCVEC(y), CVEC(r)) {
245 DEBUGMSG("dotC");
246 REQUIRES(xn==yn,BAD_SIZE);
247 REQUIRES(rn==1,BAD_SIZE);
248 DEBUGMSG("dotC");
249 KCVVIEW(x);
250 KCVVIEW(y);
251 gsl_blas_zdotu(V(x),V(y),rp);
252 OK
253}
254
255int toScalarR(int code, KRVEC(x), RVEC(r)) {
256 REQUIRES(rn==1,BAD_SIZE);
257 DEBUGMSG("toScalarR");
258 KDVVIEW(x);
259 double res;
260 switch(code) {
261 case 0: { res = gsl_blas_dnrm2(V(x)); break; }
262 case 1: { res = gsl_blas_dasum(V(x)); break; }
263 case 2: { res = gsl_vector_max_index(V(x)); break; }
264 case 3: { res = gsl_vector_max(V(x)); break; }
265 case 4: { res = gsl_vector_min_index(V(x)); break; }
266 case 5: { res = gsl_vector_min(V(x)); break; }
267 default: ERROR(BAD_CODE);
268 }
269 rp[0] = res;
270 OK
271}
272
273int toScalarF(int code, KFVEC(x), FVEC(r)) {
274 REQUIRES(rn==1,BAD_SIZE);
275 DEBUGMSG("toScalarF");
276 KFVVIEW(x);
277 float res;
278 switch(code) {
279 case 0: { res = gsl_blas_snrm2(V(x)); break; }
280 case 1: { res = gsl_blas_sasum(V(x)); break; }
281 case 2: { res = gsl_vector_float_max_index(V(x)); break; }
282 case 3: { res = gsl_vector_float_max(V(x)); break; }
283 case 4: { res = gsl_vector_float_min_index(V(x)); break; }
284 case 5: { res = gsl_vector_float_min(V(x)); break; }
285 default: ERROR(BAD_CODE);
286 }
287 rp[0] = res;
288 OK
289}
290
291
292int toScalarC(int code, KCVEC(x), RVEC(r)) {
293 REQUIRES(rn==1,BAD_SIZE);
294 DEBUGMSG("toScalarC");
295 KCVVIEW(x);
296 double res;
297 switch(code) {
298 case 0: { res = gsl_blas_dznrm2(V(x)); break; }
299 case 1: { res = gsl_blas_dzasum(V(x)); break; }
300 default: ERROR(BAD_CODE);
301 }
302 rp[0] = res;
303 OK
304}
305
306int toScalarQ(int code, KQVEC(x), FVEC(r)) {
307 REQUIRES(rn==1,BAD_SIZE);
308 DEBUGMSG("toScalarQ");
309 KQVVIEW(x);
310 float res;
311 switch(code) {
312 case 0: { res = gsl_blas_scnrm2(V(x)); break; }
313 case 1: { res = gsl_blas_scasum(V(x)); break; }
314 default: ERROR(BAD_CODE);
315 }
316 rp[0] = res;
317 OK
318}
319
320
321inline double sign(double x) {
322 if(x>0) {
323 return +1.0;
324 } else if (x<0) {
325 return -1.0;
326 } else {
327 return 0.0;
328 }
329}
330
331inline float float_sign(float x) {
332 if(x>0) {
333 return +1.0;
334 } else if (x<0) {
335 return -1.0;
336 } else {
337 return 0.0;
338 }
339}
340
341inline gsl_complex complex_abs(gsl_complex z) {
342 gsl_complex r;
343 r.dat[0] = gsl_complex_abs(z);
344 r.dat[1] = 0;
345 return r;
346}
347
348inline gsl_complex complex_signum(gsl_complex z) {
349 gsl_complex r;
350 double mag;
351 if (z.dat[0] == 0 && z.dat[1] == 0) {
352 r.dat[0] = 0;
353 r.dat[1] = 0;
354 } else {
355 mag = gsl_complex_abs(z);
356 r.dat[0] = z.dat[0]/mag;
357 r.dat[1] = z.dat[1]/mag;
358 }
359 return r;
360}
361
362#define OP(C,F) case C: { for(k=0;k<xn;k++) rp[k] = F(xp[k]); OK }
363#define OPV(C,E) case C: { for(k=0;k<xn;k++) rp[k] = E; OK }
364int mapR(int code, KRVEC(x), RVEC(r)) {
365 int k;
366 REQUIRES(xn == rn,BAD_SIZE);
367 DEBUGMSG("mapR");
368 switch (code) {
369 OP(0,sin)
370 OP(1,cos)
371 OP(2,tan)
372 OP(3,fabs)
373 OP(4,asin)
374 OP(5,acos)
375 OP(6,atan) /* atan2 mediante vectorZip */
376 OP(7,sinh)
377 OP(8,cosh)
378 OP(9,tanh)
379 OP(10,gsl_asinh)
380 OP(11,gsl_acosh)
381 OP(12,gsl_atanh)
382 OP(13,exp)
383 OP(14,log)
384 OP(15,sign)
385 OP(16,sqrt)
386 default: ERROR(BAD_CODE);
387 }
388}
389
390int mapF(int code, KFVEC(x), FVEC(r)) {
391 int k;
392 REQUIRES(xn == rn,BAD_SIZE);
393 DEBUGMSG("mapF");
394 switch (code) {
395 OP(0,sin)
396 OP(1,cos)
397 OP(2,tan)
398 OP(3,fabs)
399 OP(4,asin)
400 OP(5,acos)
401 OP(6,atan) /* atan2 mediante vectorZip */
402 OP(7,sinh)
403 OP(8,cosh)
404 OP(9,tanh)
405 OP(10,gsl_asinh)
406 OP(11,gsl_acosh)
407 OP(12,gsl_atanh)
408 OP(13,exp)
409 OP(14,log)
410 OP(15,sign)
411 OP(16,sqrt)
412 default: ERROR(BAD_CODE);
413 }
414}
415
416
417int mapCAux(int code, KGCVEC(x), GCVEC(r)) {
418 int k;
419 REQUIRES(xn == rn,BAD_SIZE);
420 DEBUGMSG("mapC");
421 switch (code) {
422 OP(0,gsl_complex_sin)
423 OP(1,gsl_complex_cos)
424 OP(2,gsl_complex_tan)
425 OP(3,complex_abs)
426 OP(4,gsl_complex_arcsin)
427 OP(5,gsl_complex_arccos)
428 OP(6,gsl_complex_arctan)
429 OP(7,gsl_complex_sinh)
430 OP(8,gsl_complex_cosh)
431 OP(9,gsl_complex_tanh)
432 OP(10,gsl_complex_arcsinh)
433 OP(11,gsl_complex_arccosh)
434 OP(12,gsl_complex_arctanh)
435 OP(13,gsl_complex_exp)
436 OP(14,gsl_complex_log)
437 OP(15,complex_signum)
438 OP(16,gsl_complex_sqrt)
439
440 // gsl_complex_arg
441 // gsl_complex_abs
442 default: ERROR(BAD_CODE);
443 }
444}
445
446int mapC(int code, KCVEC(x), CVEC(r)) {
447 return mapCAux(code, xn, (gsl_complex*)xp, rn, (gsl_complex*)rp);
448}
449
450
451gsl_complex_float complex_float_math_fun(gsl_complex (*cf)(gsl_complex), gsl_complex_float a)
452{
453 gsl_complex c;
454 gsl_complex r;
455
456 gsl_complex_float float_r;
457
458 c.dat[0] = a.dat[0];
459 c.dat[1] = a.dat[1];
460
461 r = (*cf)(c);
462
463 float_r.dat[0] = r.dat[0];
464 float_r.dat[1] = r.dat[1];
465
466 return float_r;
467}
468
469gsl_complex_float complex_float_math_op(gsl_complex (*cf)(gsl_complex,gsl_complex),
470 gsl_complex_float a,gsl_complex_float b)
471{
472 gsl_complex c1;
473 gsl_complex c2;
474 gsl_complex r;
475
476 gsl_complex_float float_r;
477
478 c1.dat[0] = a.dat[0];
479 c1.dat[1] = a.dat[1];
480
481 c2.dat[0] = b.dat[0];
482 c2.dat[1] = b.dat[1];
483
484 r = (*cf)(c1,c2);
485
486 float_r.dat[0] = r.dat[0];
487 float_r.dat[1] = r.dat[1];
488
489 return float_r;
490}
491
492#define OPC(C,F) case C: { for(k=0;k<xn;k++) rp[k] = complex_float_math_fun(&F,xp[k]); OK }
493#define OPCA(C,F,A,B) case C: { for(k=0;k<xn;k++) rp[k] = complex_float_math_op(&F,A,B); OK }
494int mapQAux(int code, KGQVEC(x), GQVEC(r)) {
495 int k;
496 REQUIRES(xn == rn,BAD_SIZE);
497 DEBUGMSG("mapQ");
498 switch (code) {
499 OPC(0,gsl_complex_sin)
500 OPC(1,gsl_complex_cos)
501 OPC(2,gsl_complex_tan)
502 OPC(3,complex_abs)
503 OPC(4,gsl_complex_arcsin)
504 OPC(5,gsl_complex_arccos)
505 OPC(6,gsl_complex_arctan)
506 OPC(7,gsl_complex_sinh)
507 OPC(8,gsl_complex_cosh)
508 OPC(9,gsl_complex_tanh)
509 OPC(10,gsl_complex_arcsinh)
510 OPC(11,gsl_complex_arccosh)
511 OPC(12,gsl_complex_arctanh)
512 OPC(13,gsl_complex_exp)
513 OPC(14,gsl_complex_log)
514 OPC(15,complex_signum)
515 OPC(16,gsl_complex_sqrt)
516
517 // gsl_complex_arg
518 // gsl_complex_abs
519 default: ERROR(BAD_CODE);
520 }
521}
522
523int mapQ(int code, KQVEC(x), QVEC(r)) {
524 return mapQAux(code, xn, (gsl_complex_float*)xp, rn, (gsl_complex_float*)rp);
525}
526
527
528int mapValR(int code, double* pval, KRVEC(x), RVEC(r)) {
529 int k;
530 double val = *pval;
531 REQUIRES(xn == rn,BAD_SIZE);
532 DEBUGMSG("mapValR");
533 switch (code) {
534 OPV(0,val*xp[k])
535 OPV(1,val/xp[k])
536 OPV(2,val+xp[k])
537 OPV(3,val-xp[k])
538 OPV(4,pow(val,xp[k]))
539 OPV(5,pow(xp[k],val))
540 default: ERROR(BAD_CODE);
541 }
542}
543
544int mapValF(int code, float* pval, KFVEC(x), FVEC(r)) {
545 int k;
546 float val = *pval;
547 REQUIRES(xn == rn,BAD_SIZE);
548 DEBUGMSG("mapValF");
549 switch (code) {
550 OPV(0,val*xp[k])
551 OPV(1,val/xp[k])
552 OPV(2,val+xp[k])
553 OPV(3,val-xp[k])
554 OPV(4,pow(val,xp[k]))
555 OPV(5,pow(xp[k],val))
556 default: ERROR(BAD_CODE);
557 }
558}
559
560int mapValCAux(int code, gsl_complex* pval, KGCVEC(x), GCVEC(r)) {
561 int k;
562 gsl_complex val = *pval;
563 REQUIRES(xn == rn,BAD_SIZE);
564 DEBUGMSG("mapValC");
565 switch (code) {
566 OPV(0,gsl_complex_mul(val,xp[k]))
567 OPV(1,gsl_complex_div(val,xp[k]))
568 OPV(2,gsl_complex_add(val,xp[k]))
569 OPV(3,gsl_complex_sub(val,xp[k]))
570 OPV(4,gsl_complex_pow(val,xp[k]))
571 OPV(5,gsl_complex_pow(xp[k],val))
572 default: ERROR(BAD_CODE);
573 }
574}
575
576int mapValC(int code, gsl_complex* val, KCVEC(x), CVEC(r)) {
577 return mapValCAux(code, val, xn, (gsl_complex*)xp, rn, (gsl_complex*)rp);
578}
579
580
581int mapValQAux(int code, gsl_complex_float* pval, KQVEC(x), GQVEC(r)) {
582 int k;
583 gsl_complex_float val = *pval;
584 REQUIRES(xn == rn,BAD_SIZE);
585 DEBUGMSG("mapValQ");
586 switch (code) {
587 OPCA(0,gsl_complex_mul,val,xp[k])
588 OPCA(1,gsl_complex_div,val,xp[k])
589 OPCA(2,gsl_complex_add,val,xp[k])
590 OPCA(3,gsl_complex_sub,val,xp[k])
591 OPCA(4,gsl_complex_pow,val,xp[k])
592 OPCA(5,gsl_complex_pow,xp[k],val)
593 default: ERROR(BAD_CODE);
594 }
595}
596
597int mapValQ(int code, gsl_complex_float* val, KQVEC(x), QVEC(r)) {
598 return mapValQAux(code, val, xn, (gsl_complex_float*)xp, rn, (gsl_complex_float*)rp);
599}
600
601
602#define OPZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = E(ap[k],bp[k]); OK }
603#define OPZV(C,msg,E) case C: {DEBUGMSG(msg) res = E(V(r),V(b)); CHECK(res,res); OK }
604int zipR(int code, KRVEC(a), KRVEC(b), RVEC(r)) {
605 REQUIRES(an == bn && an == rn, BAD_SIZE);
606 int k;
607 switch(code) {
608 OPZE(4,"zipR Pow",pow)
609 OPZE(5,"zipR ATan2",atan2)
610 }
611 KDVVIEW(a);
612 KDVVIEW(b);
613 DVVIEW(r);
614 gsl_vector_memcpy(V(r),V(a));
615 int res;
616 switch(code) {
617 OPZV(0,"zipR Add",gsl_vector_add)
618 OPZV(1,"zipR Sub",gsl_vector_sub)
619 OPZV(2,"zipR Mul",gsl_vector_mul)
620 OPZV(3,"zipR Div",gsl_vector_div)
621 default: ERROR(BAD_CODE);
622 }
623}
624
625
626int zipF(int code, KFVEC(a), KFVEC(b), FVEC(r)) {
627 REQUIRES(an == bn && an == rn, BAD_SIZE);
628 int k;
629 switch(code) {
630 OPZE(4,"zipF Pow",pow)
631 OPZE(5,"zipF ATan2",atan2)
632 }
633 KFVVIEW(a);
634 KFVVIEW(b);
635 FVVIEW(r);
636 gsl_vector_float_memcpy(V(r),V(a));
637 int res;
638 switch(code) {
639 OPZV(0,"zipF Add",gsl_vector_float_add)
640 OPZV(1,"zipF Sub",gsl_vector_float_sub)
641 OPZV(2,"zipF Mul",gsl_vector_float_mul)
642 OPZV(3,"zipF Div",gsl_vector_float_div)
643 default: ERROR(BAD_CODE);
644 }
645}
646
647
648int zipCAux(int code, KGCVEC(a), KGCVEC(b), GCVEC(r)) {
649 REQUIRES(an == bn && an == rn, BAD_SIZE);
650 int k;
651 switch(code) {
652 OPZE(0,"zipC Add",gsl_complex_add)
653 OPZE(1,"zipC Sub",gsl_complex_sub)
654 OPZE(2,"zipC Mul",gsl_complex_mul)
655 OPZE(3,"zipC Div",gsl_complex_div)
656 OPZE(4,"zipC Pow",gsl_complex_pow)
657 //OPZE(5,"zipR ATan2",atan2)
658 }
659 //KCVVIEW(a);
660 //KCVVIEW(b);
661 //CVVIEW(r);
662 //gsl_vector_memcpy(V(r),V(a));
663 //int res;
664 switch(code) {
665 default: ERROR(BAD_CODE);
666 }
667}
668
669
670int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r)) {
671 return zipCAux(code, an, (gsl_complex*)ap, bn, (gsl_complex*)bp, rn, (gsl_complex*)rp);
672}
673
674
675#define OPCZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = complex_float_math_op(&E,ap[k],bp[k]); OK }
676int zipQAux(int code, KGQVEC(a), KGQVEC(b), GQVEC(r)) {
677 REQUIRES(an == bn && an == rn, BAD_SIZE);
678 int k;
679 switch(code) {
680 OPCZE(0,"zipQ Add",gsl_complex_add)
681 OPCZE(1,"zipQ Sub",gsl_complex_sub)
682 OPCZE(2,"zipQ Mul",gsl_complex_mul)
683 OPCZE(3,"zipQ Div",gsl_complex_div)
684 OPCZE(4,"zipQ Pow",gsl_complex_pow)
685 //OPZE(5,"zipR ATan2",atan2)
686 }
687 //KCVVIEW(a);
688 //KCVVIEW(b);
689 //CVVIEW(r);
690 //gsl_vector_memcpy(V(r),V(a));
691 //int res;
692 switch(code) {
693 default: ERROR(BAD_CODE);
694 }
695}
696
697
698int zipQ(int code, KQVEC(a), KQVEC(b), QVEC(r)) {
699 return zipQAux(code, an, (gsl_complex_float*)ap, bn, (gsl_complex_float*)bp, rn, (gsl_complex_float*)rp);
700}
701
702
703
704int fft(int code, KCVEC(X), CVEC(R)) {
705 REQUIRES(Xn == Rn,BAD_SIZE);
706 DEBUGMSG("fft");
707 int s = Xn;
708 gsl_fft_complex_wavetable * wavetable = gsl_fft_complex_wavetable_alloc (s);
709 gsl_fft_complex_workspace * workspace = gsl_fft_complex_workspace_alloc (s);
710 gsl_vector_const_view X = gsl_vector_const_view_array((double*)Xp, 2*Xn);
711 gsl_vector_view R = gsl_vector_view_array((double*)Rp, 2*Rn);
712 gsl_blas_dcopy(&X.vector,&R.vector);
713 if(code==0) {
714 gsl_fft_complex_forward ((double*)Rp, 1, s, wavetable, workspace);
715 } else {
716 gsl_fft_complex_inverse ((double*)Rp, 1, s, wavetable, workspace);
717 }
718 gsl_fft_complex_wavetable_free (wavetable);
719 gsl_fft_complex_workspace_free (workspace);
720 OK
721}
722
723
724int deriv(int code, double f(double, void*), double x, double h, double * result, double * abserr)
725{
726 gsl_function F;
727 F.function = f;
728 F.params = 0;
729
730 if(code==0) return gsl_deriv_central (&F, x, h, result, abserr);
731
732 if(code==1) return gsl_deriv_forward (&F, x, h, result, abserr);
733
734 if(code==2) return gsl_deriv_backward (&F, x, h, result, abserr);
735
736 return 0;
737}
738
739
740int integrate_qng(double f(double, void*), double a, double b, double aprec, double prec,
741 double *result, double*error) {
742 DEBUGMSG("integrate_qng");
743 gsl_function F;
744 F.function = f;
745 F.params = NULL;
746 size_t neval;
747 int res = gsl_integration_qng (&F, a,b, aprec, prec, result, error, &neval);
748 CHECK(res,res);
749 OK
750}
751
752int integrate_qags(double f(double,void*), double a, double b, double aprec, double prec, int w,
753 double *result, double* error) {
754 DEBUGMSG("integrate_qags");
755 gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
756 gsl_function F;
757 F.function = f;
758 F.params = NULL;
759 int res = gsl_integration_qags (&F, a,b, aprec, prec, w,wk, result, error);
760 CHECK(res,res);
761 gsl_integration_workspace_free (wk);
762 OK
763}
764
765int integrate_qagi(double f(double,void*), double aprec, double prec, int w,
766 double *result, double* error) {
767 DEBUGMSG("integrate_qagi");
768 gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
769 gsl_function F;
770 F.function = f;
771 F.params = NULL;
772 int res = gsl_integration_qagi (&F, aprec, prec, w,wk, result, error);
773 CHECK(res,res);
774 gsl_integration_workspace_free (wk);
775 OK
776}
777
778
779int integrate_qagiu(double f(double,void*), double a, double aprec, double prec, int w,
780 double *result, double* error) {
781 DEBUGMSG("integrate_qagiu");
782 gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
783 gsl_function F;
784 F.function = f;
785 F.params = NULL;
786 int res = gsl_integration_qagiu (&F, a, aprec, prec, w,wk, result, error);
787 CHECK(res,res);
788 gsl_integration_workspace_free (wk);
789 OK
790}
791
792
793int integrate_qagil(double f(double,void*), double b, double aprec, double prec, int w,
794 double *result, double* error) {
795 DEBUGMSG("integrate_qagil");
796 gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
797 gsl_function F;
798 F.function = f;
799 F.params = NULL;
800 int res = gsl_integration_qagil (&F, b, aprec, prec, w,wk, result, error);
801 CHECK(res,res);
802 gsl_integration_workspace_free (wk);
803 OK
804}
805
806int integrate_cquad(double f(double,void*), double a, double b, double aprec, double prec,
807 int w, double *result, double* error, int *neval) {
808 DEBUGMSG("integrate_cquad");
809 gsl_integration_cquad_workspace * wk = gsl_integration_cquad_workspace_alloc (w);
810 gsl_function F;
811 F.function = f;
812 F.params = NULL;
813 size_t * sneval = NULL;
814 int res = gsl_integration_cquad (&F, a, b, aprec, prec, wk, result, error, sneval);
815 *neval = *sneval;
816 CHECK(res,res);
817 gsl_integration_cquad_workspace_free (wk);
818 OK
819}
820
821
822int polySolve(KRVEC(a), CVEC(z)) {
823 DEBUGMSG("polySolve");
824 REQUIRES(an>1,BAD_SIZE);
825 gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc (an);
826 int res = gsl_poly_complex_solve ((double*)ap, an, w, (double*)zp);
827 CHECK(res,res);
828 gsl_poly_complex_workspace_free (w);
829 OK;
830}
831
832int vector_fscanf(char*filename, RVEC(a)) {
833 DEBUGMSG("gsl_vector_fscanf");
834 DVVIEW(a);
835 FILE * f = fopen(filename,"r");
836 CHECK(!f,BAD_FILE);
837 int res = gsl_vector_fscanf(f,V(a));
838 CHECK(res,res);
839 fclose (f);
840 OK
841}
842
843int vector_fprintf(char*filename, char*fmt, RVEC(a)) {
844 DEBUGMSG("gsl_vector_fprintf");
845 DVVIEW(a);
846 FILE * f = fopen(filename,"w");
847 CHECK(!f,BAD_FILE);
848 int res = gsl_vector_fprintf(f,V(a),fmt);
849 CHECK(res,res);
850 fclose (f);
851 OK
852}
853
854int vector_fread(char*filename, RVEC(a)) {
855 DEBUGMSG("gsl_vector_fread");
856 DVVIEW(a);
857 FILE * f = fopen(filename,"r");
858 CHECK(!f,BAD_FILE);
859 int res = gsl_vector_fread(f,V(a));
860 CHECK(res,res);
861 fclose (f);
862 OK
863}
864
865int vector_fwrite(char*filename, RVEC(a)) {
866 DEBUGMSG("gsl_vector_fwrite");
867 DVVIEW(a);
868 FILE * f = fopen(filename,"w");
869 CHECK(!f,BAD_FILE);
870 int res = gsl_vector_fwrite(f,V(a));
871 CHECK(res,res);
872 fclose (f);
873 OK
874}
875
876int matrix_fprintf(char*filename, char*fmt, int ro, RMAT(m)) {
877 DEBUGMSG("matrix_fprintf");
878 FILE * f = fopen(filename,"w");
879 CHECK(!f,BAD_FILE);
880 int i,j,sr,sc;
881 if (ro==1) { sr = mc; sc = 1;} else { sr = 1; sc = mr;}
882 #define AT(M,r,c) (M##p[(r)*sr+(c)*sc])
883 for (i=0; i<mr; i++) {
884 for (j=0; j<mc-1; j++) {
885 fprintf(f,fmt,AT(m,i,j));
886 fprintf(f," ");
887 }
888 fprintf(f,fmt,AT(m,i,j));
889 fprintf(f,"\n");
890 }
891 fclose (f);
892 OK
893}
894
895//---------------------------------------------------------------
896
897typedef double Trawfun(int, double*);
898
899double only_f_aux_min(const gsl_vector*x, void *pars) {
900 Trawfun * f = (Trawfun*) pars;
901 double* p = (double*)calloc(x->size,sizeof(double));
902 int k;
903 for(k=0;k<x->size;k++) {
904 p[k] = gsl_vector_get(x,k);
905 }
906 double res = f(x->size,p);
907 free(p);
908 return res;
909}
910
911double only_f_aux_root(double x, void *pars);
912int uniMinimize(int method, double f(double),
913 double epsrel, int maxit, double min,
914 double xl, double xu, RMAT(sol)) {
915 REQUIRES(solr == maxit && solc == 4,BAD_SIZE);
916 DEBUGMSG("minimize_only_f");
917 gsl_function my_func;
918 my_func.function = only_f_aux_root;
919 my_func.params = f;
920 size_t iter = 0;
921 int status;
922 const gsl_min_fminimizer_type *T;
923 gsl_min_fminimizer *s;
924 // Starting point
925 switch(method) {
926 case 0 : {T = gsl_min_fminimizer_goldensection; break; }
927 case 1 : {T = gsl_min_fminimizer_brent; break; }
928 case 2 : {T = gsl_min_fminimizer_quad_golden; break; }
929 default: ERROR(BAD_CODE);
930 }
931 s = gsl_min_fminimizer_alloc (T);
932 gsl_min_fminimizer_set (s, &my_func, min, xl, xu);
933 do {
934 double current_min, current_lo, current_hi;
935 status = gsl_min_fminimizer_iterate (s);
936 current_min = gsl_min_fminimizer_x_minimum (s);
937 current_lo = gsl_min_fminimizer_x_lower (s);
938 current_hi = gsl_min_fminimizer_x_upper (s);
939 solp[iter*solc] = iter + 1;
940 solp[iter*solc+1] = current_min;
941 solp[iter*solc+2] = current_lo;
942 solp[iter*solc+3] = current_hi;
943 iter++;
944 if (status) /* check if solver is stuck */
945 break;
946
947 status =
948 gsl_min_test_interval (current_lo, current_hi, 0, epsrel);
949 }
950 while (status == GSL_CONTINUE && iter < maxit);
951 int i;
952 for (i=iter; i<solr; i++) {
953 solp[i*solc+0] = iter;
954 solp[i*solc+1]=0.;
955 solp[i*solc+2]=0.;
956 solp[i*solc+3]=0.;
957 }
958 gsl_min_fminimizer_free(s);
959 OK
960}
961
962
963
964// this version returns info about intermediate steps
965int minimize(int method, double f(int, double*), double tolsize, int maxit,
966 KRVEC(xi), KRVEC(sz), RMAT(sol)) {
967 REQUIRES(xin==szn && solr == maxit && solc == 3+xin,BAD_SIZE);
968 DEBUGMSG("minimizeList (nmsimplex)");
969 gsl_multimin_function my_func;
970 // extract function from pars
971 my_func.f = only_f_aux_min;
972 my_func.n = xin;
973 my_func.params = f;
974 size_t iter = 0;
975 int status;
976 double size;
977 const gsl_multimin_fminimizer_type *T;
978 gsl_multimin_fminimizer *s = NULL;
979 // Initial vertex size vector
980 KDVVIEW(sz);
981 // Starting point
982 KDVVIEW(xi);
983 // Minimizer nmsimplex, without derivatives
984 switch(method) {
985 case 0 : {T = gsl_multimin_fminimizer_nmsimplex; break; }
986#ifdef GSL110
987 case 1 : {T = gsl_multimin_fminimizer_nmsimplex; break; }
988#else
989 case 1 : {T = gsl_multimin_fminimizer_nmsimplex2; break; }
990#endif
991 default: ERROR(BAD_CODE);
992 }
993 s = gsl_multimin_fminimizer_alloc (T, my_func.n);
994 gsl_multimin_fminimizer_set (s, &my_func, V(xi), V(sz));
995 do {
996 status = gsl_multimin_fminimizer_iterate (s);
997 size = gsl_multimin_fminimizer_size (s);
998
999 solp[iter*solc+0] = iter+1;
1000 solp[iter*solc+1] = s->fval;
1001 solp[iter*solc+2] = size;
1002
1003 int k;
1004 for(k=0;k<xin;k++) {
1005 solp[iter*solc+k+3] = gsl_vector_get(s->x,k);
1006 }
1007 iter++;
1008 if (status) break;
1009 status = gsl_multimin_test_size (size, tolsize);
1010 } while (status == GSL_CONTINUE && iter < maxit);
1011 int i,j;
1012 for (i=iter; i<solr; i++) {
1013 solp[i*solc+0] = iter;
1014 for(j=1;j<solc;j++) {
1015 solp[i*solc+j]=0.;
1016 }
1017 }
1018 gsl_multimin_fminimizer_free(s);
1019 OK
1020}
1021
1022// working with the gradient
1023
1024typedef struct {double (*f)(int, double*); int (*df)(int, double*, int, double*);} Tfdf;
1025
1026double f_aux_min(const gsl_vector*x, void *pars) {
1027 Tfdf * fdf = ((Tfdf*) pars);
1028 double* p = (double*)calloc(x->size,sizeof(double));
1029 int k;
1030 for(k=0;k<x->size;k++) {
1031 p[k] = gsl_vector_get(x,k);
1032 }
1033 double res = fdf->f(x->size,p);
1034 free(p);
1035 return res;
1036}
1037
1038
1039void df_aux_min(const gsl_vector * x, void * pars, gsl_vector * g) {
1040 Tfdf * fdf = ((Tfdf*) pars);
1041 double* p = (double*)calloc(x->size,sizeof(double));
1042 double* q = (double*)calloc(g->size,sizeof(double));
1043 int k;
1044 for(k=0;k<x->size;k++) {
1045 p[k] = gsl_vector_get(x,k);
1046 }
1047
1048 fdf->df(x->size,p,g->size,q);
1049
1050 for(k=0;k<x->size;k++) {
1051 gsl_vector_set(g,k,q[k]);
1052 }
1053 free(p);
1054 free(q);
1055}
1056
1057void fdf_aux_min(const gsl_vector * x, void * pars, double * f, gsl_vector * g) {
1058 *f = f_aux_min(x,pars);
1059 df_aux_min(x,pars,g);
1060}
1061
1062
1063int minimizeD(int method, double f(int, double*), int df(int, double*, int, double*),
1064 double initstep, double minimpar, double tolgrad, int maxit,
1065 KRVEC(xi), RMAT(sol)) {
1066 REQUIRES(solr == maxit && solc == 2+xin,BAD_SIZE);
1067 DEBUGMSG("minimizeWithDeriv (conjugate_fr)");
1068 gsl_multimin_function_fdf my_func;
1069 // extract function from pars
1070 my_func.f = f_aux_min;
1071 my_func.df = df_aux_min;
1072 my_func.fdf = fdf_aux_min;
1073 my_func.n = xin;
1074 Tfdf stfdf;
1075 stfdf.f = f;
1076 stfdf.df = df;
1077 my_func.params = &stfdf;
1078 size_t iter = 0;
1079 int status;
1080 const gsl_multimin_fdfminimizer_type *T;
1081 gsl_multimin_fdfminimizer *s = NULL;
1082 // Starting point
1083 KDVVIEW(xi);
1084 // conjugate gradient fr
1085 switch(method) {
1086 case 0 : {T = gsl_multimin_fdfminimizer_conjugate_fr; break; }
1087 case 1 : {T = gsl_multimin_fdfminimizer_conjugate_pr; break; }
1088 case 2 : {T = gsl_multimin_fdfminimizer_vector_bfgs; break; }
1089 case 3 : {T = gsl_multimin_fdfminimizer_vector_bfgs2; break; }
1090 case 4 : {T = gsl_multimin_fdfminimizer_steepest_descent; break; }
1091 default: ERROR(BAD_CODE);
1092 }
1093 s = gsl_multimin_fdfminimizer_alloc (T, my_func.n);
1094 gsl_multimin_fdfminimizer_set (s, &my_func, V(xi), initstep, minimpar);
1095 do {
1096 status = gsl_multimin_fdfminimizer_iterate (s);
1097 solp[iter*solc+0] = iter+1;
1098 solp[iter*solc+1] = s->f;
1099 int k;
1100 for(k=0;k<xin;k++) {
1101 solp[iter*solc+k+2] = gsl_vector_get(s->x,k);
1102 }
1103 iter++;
1104 if (status) break;
1105 status = gsl_multimin_test_gradient (s->gradient, tolgrad);
1106 } while (status == GSL_CONTINUE && iter < maxit);
1107 int i,j;
1108 for (i=iter; i<solr; i++) {
1109 solp[i*solc+0] = iter;
1110 for(j=1;j<solc;j++) {
1111 solp[i*solc+j]=0.;
1112 }
1113 }
1114 gsl_multimin_fdfminimizer_free(s);
1115 OK
1116}
1117
1118//---------------------------------------------------------------
1119
1120double only_f_aux_root(double x, void *pars) {
1121 double (*f)(double) = (double (*)(double)) pars;
1122 return f(x);
1123}
1124
1125int root(int method, double f(double),
1126 double epsrel, int maxit,
1127 double xl, double xu, RMAT(sol)) {
1128 REQUIRES(solr == maxit && solc == 4,BAD_SIZE);
1129 DEBUGMSG("root_only_f");
1130 gsl_function my_func;
1131 // extract function from pars
1132 my_func.function = only_f_aux_root;
1133 my_func.params = f;
1134 size_t iter = 0;
1135 int status;
1136 const gsl_root_fsolver_type *T;
1137 gsl_root_fsolver *s;
1138 // Starting point
1139 switch(method) {
1140 case 0 : {T = gsl_root_fsolver_bisection; printf("7\n"); break; }
1141 case 1 : {T = gsl_root_fsolver_falsepos; break; }
1142 case 2 : {T = gsl_root_fsolver_brent; break; }
1143 default: ERROR(BAD_CODE);
1144 }
1145 s = gsl_root_fsolver_alloc (T);
1146 gsl_root_fsolver_set (s, &my_func, xl, xu);
1147 do {
1148 double best, current_lo, current_hi;
1149 status = gsl_root_fsolver_iterate (s);
1150 best = gsl_root_fsolver_root (s);
1151 current_lo = gsl_root_fsolver_x_lower (s);
1152 current_hi = gsl_root_fsolver_x_upper (s);
1153 solp[iter*solc] = iter + 1;
1154 solp[iter*solc+1] = best;
1155 solp[iter*solc+2] = current_lo;
1156 solp[iter*solc+3] = current_hi;
1157 iter++;
1158 if (status) /* check if solver is stuck */
1159 break;
1160
1161 status =
1162 gsl_root_test_interval (current_lo, current_hi, 0, epsrel);
1163 }
1164 while (status == GSL_CONTINUE && iter < maxit);
1165 int i;
1166 for (i=iter; i<solr; i++) {
1167 solp[i*solc+0] = iter;
1168 solp[i*solc+1]=0.;
1169 solp[i*solc+2]=0.;
1170 solp[i*solc+3]=0.;
1171 }
1172 gsl_root_fsolver_free(s);
1173 OK
1174}
1175
1176typedef struct {
1177 double (*f)(double);
1178 double (*jf)(double);
1179} uniTfjf;
1180
1181double f_aux_uni(double x, void *pars) {
1182 uniTfjf * fjf = ((uniTfjf*) pars);
1183 return (fjf->f)(x);
1184}
1185
1186double jf_aux_uni(double x, void * pars) {
1187 uniTfjf * fjf = ((uniTfjf*) pars);
1188 return (fjf->jf)(x);
1189}
1190
1191void fjf_aux_uni(double x, void * pars, double * f, double * g) {
1192 *f = f_aux_uni(x,pars);
1193 *g = jf_aux_uni(x,pars);
1194}
1195
1196int rootj(int method, double f(double),
1197 double df(double),
1198 double epsrel, int maxit,
1199 double x, RMAT(sol)) {
1200 REQUIRES(solr == maxit && solc == 2,BAD_SIZE);
1201 DEBUGMSG("root_fjf");
1202 gsl_function_fdf my_func;
1203 // extract function from pars
1204 my_func.f = f_aux_uni;
1205 my_func.df = jf_aux_uni;
1206 my_func.fdf = fjf_aux_uni;
1207 uniTfjf stfjf;
1208 stfjf.f = f;
1209 stfjf.jf = df;
1210 my_func.params = &stfjf;
1211 size_t iter = 0;
1212 int status;
1213 const gsl_root_fdfsolver_type *T;
1214 gsl_root_fdfsolver *s;
1215 // Starting point
1216 switch(method) {
1217 case 0 : {T = gsl_root_fdfsolver_newton;; break; }
1218 case 1 : {T = gsl_root_fdfsolver_secant; break; }
1219 case 2 : {T = gsl_root_fdfsolver_steffenson; break; }
1220 default: ERROR(BAD_CODE);
1221 }
1222 s = gsl_root_fdfsolver_alloc (T);
1223
1224 gsl_root_fdfsolver_set (s, &my_func, x);
1225
1226 do {
1227 double x0;
1228 status = gsl_root_fdfsolver_iterate (s);
1229 x0 = x;
1230 x = gsl_root_fdfsolver_root(s);
1231 solp[iter*solc+0] = iter+1;
1232 solp[iter*solc+1] = x;
1233
1234 iter++;
1235 if (status) /* check if solver is stuck */
1236 break;
1237
1238 status =
1239 gsl_root_test_delta (x, x0, 0, epsrel);
1240 }
1241 while (status == GSL_CONTINUE && iter < maxit);
1242
1243 int i;
1244 for (i=iter; i<solr; i++) {
1245 solp[i*solc+0] = iter;
1246 solp[i*solc+1]=0.;
1247 }
1248 gsl_root_fdfsolver_free(s);
1249 OK
1250}
1251
1252
1253//---------------------------------------------------------------
1254
1255typedef void TrawfunV(int, double*, int, double*);
1256
1257int only_f_aux_multiroot(const gsl_vector*x, void *pars, gsl_vector*y) {
1258 TrawfunV * f = (TrawfunV*) pars;
1259 double* p = (double*)calloc(x->size,sizeof(double));
1260 double* q = (double*)calloc(y->size,sizeof(double));
1261 int k;
1262 for(k=0;k<x->size;k++) {
1263 p[k] = gsl_vector_get(x,k);
1264 }
1265 f(x->size,p,y->size,q);
1266 for(k=0;k<y->size;k++) {
1267 gsl_vector_set(y,k,q[k]);
1268 }
1269 free(p);
1270 free(q);
1271 return 0; //hmmm
1272}
1273
1274int multiroot(int method, void f(int, double*, int, double*),
1275 double epsabs, int maxit,
1276 KRVEC(xi), RMAT(sol)) {
1277 REQUIRES(solr == maxit && solc == 1+2*xin,BAD_SIZE);
1278 DEBUGMSG("root_only_f");
1279 gsl_multiroot_function my_func;
1280 // extract function from pars
1281 my_func.f = only_f_aux_multiroot;
1282 my_func.n = xin;
1283 my_func.params = f;
1284 size_t iter = 0;
1285 int status;
1286 const gsl_multiroot_fsolver_type *T;
1287 gsl_multiroot_fsolver *s;
1288 // Starting point
1289 KDVVIEW(xi);
1290 switch(method) {
1291 case 0 : {T = gsl_multiroot_fsolver_hybrids;; break; }
1292 case 1 : {T = gsl_multiroot_fsolver_hybrid; break; }
1293 case 2 : {T = gsl_multiroot_fsolver_dnewton; break; }
1294 case 3 : {T = gsl_multiroot_fsolver_broyden; break; }
1295 default: ERROR(BAD_CODE);
1296 }
1297 s = gsl_multiroot_fsolver_alloc (T, my_func.n);
1298 gsl_multiroot_fsolver_set (s, &my_func, V(xi));
1299
1300 do {
1301 status = gsl_multiroot_fsolver_iterate (s);
1302
1303 solp[iter*solc+0] = iter+1;
1304
1305 int k;
1306 for(k=0;k<xin;k++) {
1307 solp[iter*solc+k+1] = gsl_vector_get(s->x,k);
1308 }
1309 for(k=xin;k<2*xin;k++) {
1310 solp[iter*solc+k+1] = gsl_vector_get(s->f,k-xin);
1311 }
1312
1313 iter++;
1314 if (status) /* check if solver is stuck */
1315 break;
1316
1317 status =
1318 gsl_multiroot_test_residual (s->f, epsabs);
1319 }
1320 while (status == GSL_CONTINUE && iter < maxit);
1321
1322 int i,j;
1323 for (i=iter; i<solr; i++) {
1324 solp[i*solc+0] = iter;
1325 for(j=1;j<solc;j++) {
1326 solp[i*solc+j]=0.;
1327 }
1328 }
1329 gsl_multiroot_fsolver_free(s);
1330 OK
1331}
1332
1333// working with the jacobian
1334
1335typedef struct {int (*f)(int, double*, int, double *);
1336 int (*jf)(int, double*, int, int, double*);} Tfjf;
1337
1338int f_aux(const gsl_vector*x, void *pars, gsl_vector*y) {
1339 Tfjf * fjf = ((Tfjf*) pars);
1340 double* p = (double*)calloc(x->size,sizeof(double));
1341 double* q = (double*)calloc(y->size,sizeof(double));
1342 int k;
1343 for(k=0;k<x->size;k++) {
1344 p[k] = gsl_vector_get(x,k);
1345 }
1346 (fjf->f)(x->size,p,y->size,q);
1347 for(k=0;k<y->size;k++) {
1348 gsl_vector_set(y,k,q[k]);
1349 }
1350 free(p);
1351 free(q);
1352 return 0;
1353}
1354
1355int jf_aux(const gsl_vector * x, void * pars, gsl_matrix * jac) {
1356 Tfjf * fjf = ((Tfjf*) pars);
1357 double* p = (double*)calloc(x->size,sizeof(double));
1358 double* q = (double*)calloc((jac->size1)*(jac->size2),sizeof(double));
1359 int i,j,k;
1360 for(k=0;k<x->size;k++) {
1361 p[k] = gsl_vector_get(x,k);
1362 }
1363
1364 (fjf->jf)(x->size,p,jac->size1,jac->size2,q);
1365
1366 k=0;
1367 for(i=0;i<jac->size1;i++) {
1368 for(j=0;j<jac->size2;j++){
1369 gsl_matrix_set(jac,i,j,q[k++]);
1370 }
1371 }
1372 free(p);
1373 free(q);
1374 return 0;
1375}
1376
1377int fjf_aux(const gsl_vector * x, void * pars, gsl_vector * f, gsl_matrix * g) {
1378 f_aux(x,pars,f);
1379 jf_aux(x,pars,g);
1380 return 0;
1381}
1382
1383int multirootj(int method, int f(int, double*, int, double*),
1384 int jac(int, double*, int, int, double*),
1385 double epsabs, int maxit,
1386 KRVEC(xi), RMAT(sol)) {
1387 REQUIRES(solr == maxit && solc == 1+2*xin,BAD_SIZE);
1388 DEBUGMSG("root_fjf");
1389 gsl_multiroot_function_fdf my_func;
1390 // extract function from pars
1391 my_func.f = f_aux;
1392 my_func.df = jf_aux;
1393 my_func.fdf = fjf_aux;
1394 my_func.n = xin;
1395 Tfjf stfjf;
1396 stfjf.f = f;
1397 stfjf.jf = jac;
1398 my_func.params = &stfjf;
1399 size_t iter = 0;
1400 int status;
1401 const gsl_multiroot_fdfsolver_type *T;
1402 gsl_multiroot_fdfsolver *s;
1403 // Starting point
1404 KDVVIEW(xi);
1405 switch(method) {
1406 case 0 : {T = gsl_multiroot_fdfsolver_hybridsj;; break; }
1407 case 1 : {T = gsl_multiroot_fdfsolver_hybridj; break; }
1408 case 2 : {T = gsl_multiroot_fdfsolver_newton; break; }
1409 case 3 : {T = gsl_multiroot_fdfsolver_gnewton; break; }
1410 default: ERROR(BAD_CODE);
1411 }
1412 s = gsl_multiroot_fdfsolver_alloc (T, my_func.n);
1413
1414 gsl_multiroot_fdfsolver_set (s, &my_func, V(xi));
1415
1416 do {
1417 status = gsl_multiroot_fdfsolver_iterate (s);
1418
1419 solp[iter*solc+0] = iter+1;
1420
1421 int k;
1422 for(k=0;k<xin;k++) {
1423 solp[iter*solc+k+1] = gsl_vector_get(s->x,k);
1424 }
1425 for(k=xin;k<2*xin;k++) {
1426 solp[iter*solc+k+1] = gsl_vector_get(s->f,k-xin);
1427 }
1428
1429 iter++;
1430 if (status) /* check if solver is stuck */
1431 break;
1432
1433 status =
1434 gsl_multiroot_test_residual (s->f, epsabs);
1435 }
1436 while (status == GSL_CONTINUE && iter < maxit);
1437
1438 int i,j;
1439 for (i=iter; i<solr; i++) {
1440 solp[i*solc+0] = iter;
1441 for(j=1;j<solc;j++) {
1442 solp[i*solc+j]=0.;
1443 }
1444 }
1445 gsl_multiroot_fdfsolver_free(s);
1446 OK
1447}
1448
1449//-------------- non linear least squares fitting -------------------
1450
1451int nlfit(int method, int f(int, double*, int, double*),
1452 int jac(int, double*, int, int, double*),
1453 double epsabs, double epsrel, int maxit, int p,
1454 KRVEC(xi), RMAT(sol)) {
1455 REQUIRES(solr == maxit && solc == 2+xin,BAD_SIZE);
1456 DEBUGMSG("nlfit");
1457 const gsl_multifit_fdfsolver_type *T;
1458 gsl_multifit_fdfsolver *s;
1459 gsl_multifit_function_fdf my_f;
1460 // extract function from pars
1461 my_f.f = f_aux;
1462 my_f.df = jf_aux;
1463 my_f.fdf = fjf_aux;
1464 my_f.n = p;
1465 my_f.p = xin; // !!!!
1466 Tfjf stfjf;
1467 stfjf.f = f;
1468 stfjf.jf = jac;
1469 my_f.params = &stfjf;
1470 size_t iter = 0;
1471 int status;
1472
1473 KDVVIEW(xi);
1474 //DMVIEW(cov);
1475
1476 switch(method) {
1477 case 0 : { T = gsl_multifit_fdfsolver_lmsder; break; }
1478 case 1 : { T = gsl_multifit_fdfsolver_lmder; break; }
1479 default: ERROR(BAD_CODE);
1480 }
1481
1482 s = gsl_multifit_fdfsolver_alloc (T, my_f.n, my_f.p);
1483 gsl_multifit_fdfsolver_set (s, &my_f, V(xi));
1484
1485 do { status = gsl_multifit_fdfsolver_iterate (s);
1486
1487 solp[iter*solc+0] = iter+1;
1488 solp[iter*solc+1] = gsl_blas_dnrm2 (s->f);
1489
1490 int k;
1491 for(k=0;k<xin;k++) {
1492 solp[iter*solc+k+2] = gsl_vector_get(s->x,k);
1493 }
1494
1495 iter++;
1496 if (status) /* check if solver is stuck */
1497 break;
1498
1499 status = gsl_multifit_test_delta (s->dx, s->x, epsabs, epsrel);
1500 }
1501 while (status == GSL_CONTINUE && iter < maxit);
1502
1503 int i,j;
1504 for (i=iter; i<solr; i++) {
1505 solp[i*solc+0] = iter;
1506 for(j=1;j<solc;j++) {
1507 solp[i*solc+j]=0.;
1508 }
1509 }
1510
1511 //gsl_multifit_covar (s->J, 0.0, M(cov));
1512
1513 gsl_multifit_fdfsolver_free (s);
1514 OK
1515}
1516
1517
1518//////////////////////////////////////////////////////
1519
1520
1521#define RAN(C,F) case C: { for(k=0;k<rn;k++) { rp[k]= F(gen); }; OK }
1522
1523int random_vector(int seed, int code, RVEC(r)) {
1524 DEBUGMSG("random_vector")
1525 static gsl_rng * gen = NULL;
1526 if (!gen) { gen = gsl_rng_alloc (gsl_rng_mt19937);}
1527 gsl_rng_set (gen, seed);
1528 int k;
1529 switch (code) {
1530 RAN(0,gsl_rng_uniform)
1531 RAN(1,gsl_ran_ugaussian)
1532 default: ERROR(BAD_CODE);
1533 }
1534}
1535#undef RAN
1536
1537//////////////////////////////////////////////////////
1538
1539#include "gsl-ode.c"
1540
1541//////////////////////////////////////////////////////
diff --git a/packages/hmatrix/src/Numeric/GSL/gsl-ode.c b/packages/hmatrix/src/Numeric/GSL/gsl-ode.c
new file mode 100644
index 0000000..3f2771b
--- /dev/null
+++ b/packages/hmatrix/src/Numeric/GSL/gsl-ode.c
@@ -0,0 +1,182 @@
1
2#ifdef GSLODE1
3
4////////////////////////////// ODE V1 //////////////////////////////////////////
5
6#include <gsl/gsl_odeiv.h>
7
8typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode;
9
10int odefunc (double t, const double y[], double f[], void *params) {
11 Tode * P = (Tode*) params;
12 (P->f)(t,P->n,y,P->n,f);
13 return GSL_SUCCESS;
14}
15
16int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) {
17 Tode * P = ((Tode*) params);
18 (P->j)(t,P->n,y,P->n,P->n,dfdy);
19 int j;
20 for (j=0; j< P->n; j++)
21 dfdt[j] = 0.0;
22 return GSL_SUCCESS;
23}
24
25
26int ode(int method, double h, double eps_abs, double eps_rel,
27 int f(double, int, const double*, int, double*),
28 int jac(double, int, const double*, int, int, double*),
29 KRVEC(xi), KRVEC(ts), RMAT(sol)) {
30
31 const gsl_odeiv_step_type * T;
32
33 switch(method) {
34 case 0 : {T = gsl_odeiv_step_rk2; break; }
35 case 1 : {T = gsl_odeiv_step_rk4; break; }
36 case 2 : {T = gsl_odeiv_step_rkf45; break; }
37 case 3 : {T = gsl_odeiv_step_rkck; break; }
38 case 4 : {T = gsl_odeiv_step_rk8pd; break; }
39 case 5 : {T = gsl_odeiv_step_rk2imp; break; }
40 case 6 : {T = gsl_odeiv_step_rk4imp; break; }
41 case 7 : {T = gsl_odeiv_step_bsimp; break; }
42 case 8 : { printf("Sorry: ODE rk1imp not available in this GSL version\n"); exit(0); }
43 case 9 : { printf("Sorry: ODE msadams not available in this GSL version\n"); exit(0); }
44 case 10: { printf("Sorry: ODE msbdf not available in this GSL version\n"); exit(0); }
45 default: ERROR(BAD_CODE);
46 }
47
48 gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, xin);
49 gsl_odeiv_control * c = gsl_odeiv_control_y_new (eps_abs, eps_rel);
50 gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (xin);
51
52 Tode P;
53 P.f = f;
54 P.j = jac;
55 P.n = xin;
56
57 gsl_odeiv_system sys = {odefunc, odejac, xin, &P};
58
59 double t = tsp[0];
60
61 double* y = (double*)calloc(xin,sizeof(double));
62 int i,j;
63 for(i=0; i< xin; i++) {
64 y[i] = xip[i];
65 solp[i] = xip[i];
66 }
67
68 for (i = 1; i < tsn ; i++)
69 {
70 double ti = tsp[i];
71 while (t < ti)
72 {
73 gsl_odeiv_evolve_apply (e, c, s,
74 &sys,
75 &t, ti, &h,
76 y);
77 // if (h < hmin) h = hmin;
78 }
79 for(j=0; j<xin; j++) {
80 solp[i*xin + j] = y[j];
81 }
82 }
83
84 free(y);
85 gsl_odeiv_evolve_free (e);
86 gsl_odeiv_control_free (c);
87 gsl_odeiv_step_free (s);
88 return 0;
89}
90
91#else
92
93///////////////////// ODE V2 ///////////////////////////////////////////////////
94
95#include <gsl/gsl_odeiv2.h>
96
97typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode;
98
99int odefunc (double t, const double y[], double f[], void *params) {
100 Tode * P = (Tode*) params;
101 (P->f)(t,P->n,y,P->n,f);
102 return GSL_SUCCESS;
103}
104
105int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) {
106 Tode * P = ((Tode*) params);
107 (P->j)(t,P->n,y,P->n,P->n,dfdy);
108 int j;
109 for (j=0; j< P->n; j++)
110 dfdt[j] = 0.0;
111 return GSL_SUCCESS;
112}
113
114
115int ode(int method, double h, double eps_abs, double eps_rel,
116 int f(double, int, const double*, int, double*),
117 int jac(double, int, const double*, int, int, double*),
118 KRVEC(xi), KRVEC(ts), RMAT(sol)) {
119
120 const gsl_odeiv2_step_type * T;
121
122 switch(method) {
123 case 0 : {T = gsl_odeiv2_step_rk2; break; }
124 case 1 : {T = gsl_odeiv2_step_rk4; break; }
125 case 2 : {T = gsl_odeiv2_step_rkf45; break; }
126 case 3 : {T = gsl_odeiv2_step_rkck; break; }
127 case 4 : {T = gsl_odeiv2_step_rk8pd; break; }
128 case 5 : {T = gsl_odeiv2_step_rk2imp; break; }
129 case 6 : {T = gsl_odeiv2_step_rk4imp; break; }
130 case 7 : {T = gsl_odeiv2_step_bsimp; break; }
131 case 8 : {T = gsl_odeiv2_step_rk1imp; break; }
132 case 9 : {T = gsl_odeiv2_step_msadams; break; }
133 case 10: {T = gsl_odeiv2_step_msbdf; break; }
134 default: ERROR(BAD_CODE);
135 }
136
137 Tode P;
138 P.f = f;
139 P.j = jac;
140 P.n = xin;
141
142 gsl_odeiv2_system sys = {odefunc, odejac, xin, &P};
143
144 gsl_odeiv2_driver * d =
145 gsl_odeiv2_driver_alloc_y_new (&sys, T, h, eps_abs, eps_rel);
146
147 double t = tsp[0];
148
149 double* y = (double*)calloc(xin,sizeof(double));
150 int i,j;
151 int status=0;
152 for(i=0; i< xin; i++) {
153 y[i] = xip[i];
154 solp[i] = xip[i];
155 }
156
157 for (i = 1; i < tsn ; i++)
158 {
159 double ti = tsp[i];
160
161 status = gsl_odeiv2_driver_apply (d, &t, ti, y);
162
163 if (status != GSL_SUCCESS) {
164 printf ("error in ode, return value=%d\n", status);
165 break;
166 }
167
168// printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
169
170 for(j=0; j<xin; j++) {
171 solp[i*xin + j] = y[j];
172 }
173 }
174
175 free(y);
176 gsl_odeiv2_driver_free (d);
177
178 return status;
179}
180
181#endif
182
generated by cgit v1.2.3 (git 2.39.1) at 2024-07-05 17:18:15 +0000