summaryrefslogtreecommitdiff
path: root/lib/Numeric/GSL/Root.hs
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Numeric/GSL/Root.hs')
-rw-r--r--lib/Numeric/GSL/Root.hs199
1 files changed, 0 insertions, 199 deletions
diff --git a/lib/Numeric/GSL/Root.hs b/lib/Numeric/GSL/Root.hs
deleted file mode 100644
index 9d561c4..0000000
--- a/lib/Numeric/GSL/Root.hs
+++ /dev/null
@@ -1,199 +0,0 @@
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"
generated by cgit v1.2.3 (git 2.39.1) at 2024-09-15 09:20:01 +0000