blob: ad1b72c27017f6a861106c16d9f6ffe383bc0aa5 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
|
{- |
Module : Numeric.GSL.Root
Copyright : (c) Alberto Ruiz 2009
License : GPL
Maintainer : Alberto Ruiz (aruiz at um dot es)
Stability : provisional
Portability : uses ffi
Multidimensional root finding.
<http://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Root_002dFinding.html>
The example in the GSL manual:
@import Numeric.GSL
import Numeric.LinearAlgebra(format)
import Text.Printf(printf)
rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ]
disp = putStrLn . format \" \" (printf \"%.3f\")
main = do
let (sol,path) = root Hybrids 1E-7 30 (rosenbrock 1 10) [-10,-5]
print sol
disp path
\> main
[1.0,1.0]
0.000 -10.000 -5.000 11.000 -1050.000
1.000 -3.976 24.827 4.976 90.203
2.000 -3.976 24.827 4.976 90.203
3.000 -3.976 24.827 4.976 90.203
4.000 -1.274 -5.680 2.274 -73.018
5.000 -1.274 -5.680 2.274 -73.018
6.000 0.249 0.298 0.751 2.359
7.000 0.249 0.298 0.751 2.359
8.000 1.000 0.878 -0.000 -1.218
9.000 1.000 0.989 -0.000 -0.108
10.000 1.000 1.000 0.000 0.000
@
-}
-----------------------------------------------------------------------------
module Numeric.GSL.Root (
root, RootMethod(..)
) where
import Data.Packed.Internal
import Data.Packed.Matrix
import Foreign
import Foreign.C.Types(CInt)
-------------------------------------------------------------------------
data RootMethod = Hybrids
| Hybrid
| DNewton
| Broyden
deriving (Enum,Eq,Show)
-- | Nonlinear multidimensional root finding using algorithms that do not require
-- any derivative information to be supplied by the user.
-- Any derivatives needed are approximated by finite differences.
root :: RootMethod
-> Double -- ^ maximum residual
-> Int -- ^ maximum number of iterations allowed
-> ([Double] -> [Double]) -- ^ function to minimize
-> [Double] -- ^ starting point
-> ([Double], Matrix Double) -- ^ solution vector and optimization path
root method epsabs maxit fun xinit = rootGen (fi (fromEnum method)) fun xinit epsabs maxit
rootGen m f xi epsabs maxit = unsafePerformIO $ do
let xiv = fromList xi
n = dim xiv
fp <- mkVecVecfun (aux_vTov (fromList.f.toList))
rawpath <- withVector xiv $ \xiv' ->
createMIO maxit (2*n+1)
(c_root m fp epsabs (fi maxit) // xiv')
"root"
let it = round (rawpath @@> (maxit-1,0))
path = takeRows it rawpath
[sol] = toLists $ dropRows (it-1) path
freeHaskellFunPtr fp
return (take n $ drop 1 sol, path)
foreign import ccall "root"
c_root:: CInt -> FunPtr (CInt -> Ptr Double -> Ptr Double -> IO ()) -> Double -> CInt -> TVM
---------------------------------------------------------------------
foreign import ccall "wrapper"
mkVecVecfun :: (CInt -> Ptr Double -> Ptr Double -> IO ())
-> IO (FunPtr (CInt -> Ptr Double -> Ptr Double->IO()))
aux_vTov :: (Vector Double -> Vector Double) -> (CInt -> Ptr Double -> Ptr Double -> IO())
aux_vTov f n p r = g where
V {fptr = pr} = f x
x = createV (fromIntegral n) copy "aux_vTov"
copy n' q = do
copyArray q p (fromIntegral n')
return 0
g = withForeignPtr pr $ \p' -> copyArray r p' (fromIntegral n)
createV n fun msg = unsafePerformIO $ do
r <- createVector n
app1 fun vec r msg
return r
createMIO r c fun msg = do
res <- createMatrix RowMajor r c
app1 fun mat res msg
return res
|
