diff options
Diffstat (limited to 'lib/Numeric/GSL/Root.hs')
-rw-r--r-- | lib/Numeric/GSL/Root.hs | 49 |
1 files changed, 7 insertions, 42 deletions
diff --git a/lib/Numeric/GSL/Root.hs b/lib/Numeric/GSL/Root.hs index 6ce2c4c..840e8ee 100644 --- a/lib/Numeric/GSL/Root.hs +++ b/lib/Numeric/GSL/Root.hs | |||
@@ -53,6 +53,7 @@ import Data.Packed.Internal | |||
53 | import Data.Packed.Matrix | 53 | import Data.Packed.Matrix |
54 | import Foreign | 54 | import Foreign |
55 | import Foreign.C.Types(CInt) | 55 | import Foreign.C.Types(CInt) |
56 | import Numeric.GSL.Internal | ||
56 | 57 | ||
57 | ------------------------------------------------------------------------- | 58 | ------------------------------------------------------------------------- |
58 | 59 | ||
@@ -60,7 +61,7 @@ data RootMethod = Hybrids | |||
60 | | Hybrid | 61 | | Hybrid |
61 | | DNewton | 62 | | DNewton |
62 | | Broyden | 63 | | Broyden |
63 | deriving (Enum,Eq,Show) | 64 | deriving (Enum,Eq,Show,Bounded) |
64 | 65 | ||
65 | -- | Nonlinear multidimensional root finding using algorithms that do not require | 66 | -- | Nonlinear multidimensional root finding using algorithms that do not require |
66 | -- any derivative information to be supplied by the user. | 67 | -- any derivative information to be supplied by the user. |
@@ -98,7 +99,7 @@ data RootMethodJ = HybridsJ | |||
98 | | HybridJ | 99 | | HybridJ |
99 | | Newton | 100 | | Newton |
100 | | GNewton | 101 | | GNewton |
101 | deriving (Enum,Eq,Show) | 102 | deriving (Enum,Eq,Show,Bounded) |
102 | 103 | ||
103 | -- | Nonlinear multidimensional root finding using both the function and its derivatives. | 104 | -- | Nonlinear multidimensional root finding using both the function and its derivatives. |
104 | rootJ :: RootMethodJ | 105 | rootJ :: RootMethodJ |
@@ -124,57 +125,21 @@ rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do | |||
124 | path = takeRows it rawpath | 125 | path = takeRows it rawpath |
125 | [sol] = toLists $ dropRows (it-1) path | 126 | [sol] = toLists $ dropRows (it-1) path |
126 | freeHaskellFunPtr fp | 127 | freeHaskellFunPtr fp |
128 | freeHaskellFunPtr jp | ||
127 | return (take n $ drop 1 sol, path) | 129 | return (take n $ drop 1 sol, path) |
128 | 130 | ||
129 | 131 | ||
130 | foreign import ccall "rootj" | 132 | foreign import ccall "rootj" |
131 | c_rootj:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> CInt -> TVM | 133 | c_rootj:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> CInt -> TVM |
132 | 134 | ||
133 | 135 | ------------------------------------------------------- | |
134 | --------------------------------------------------------------------- | ||
135 | |||
136 | foreign import ccall "wrapper" | ||
137 | mkVecVecfun :: TVV -> IO (FunPtr TVV) | ||
138 | |||
139 | aux_vTov :: (Vector Double -> Vector Double) -> TVV | ||
140 | aux_vTov f n p nr r = g where | ||
141 | V {fptr = pr} = f x | ||
142 | x = createV (fromIntegral n) copy "aux_vTov" | ||
143 | copy n' q = do | ||
144 | copyArray q p (fromIntegral n') | ||
145 | return 0 | ||
146 | g = do withForeignPtr pr $ \p' -> copyArray r p' (fromIntegral nr) | ||
147 | return 0 | ||
148 | |||
149 | foreign import ccall "wrapper" | ||
150 | mkVecMatfun :: TVM -> IO (FunPtr TVM) | ||
151 | |||
152 | aux_vTom :: (Vector Double -> Matrix Double) -> TVM | ||
153 | aux_vTom f n p rr cr r = g where | ||
154 | V {fptr = pr} = flatten $ f x | ||
155 | x = createV (fromIntegral n) copy "aux_vTov" | ||
156 | copy n' q = do | ||
157 | copyArray q p (fromIntegral n') | ||
158 | return 0 | ||
159 | g = do withForeignPtr pr $ \p' -> copyArray r p' (fromIntegral $ rr*cr) | ||
160 | return 0 | ||
161 | |||
162 | createV n fun msg = unsafePerformIO $ do | ||
163 | r <- createVector n | ||
164 | app1 fun vec r msg | ||
165 | return r | ||
166 | |||
167 | createMIO r c fun msg = do | ||
168 | res <- createMatrix RowMajor r c | ||
169 | app1 fun mat res msg | ||
170 | return res | ||
171 | 136 | ||
172 | checkdim1 n v | 137 | checkdim1 n v |
173 | | dim v == n = v | 138 | | dim v == n = v |
174 | | otherwise = error $ "Error: "++ show n | 139 | | otherwise = error $ "Error: "++ show n |
175 | ++ " results expected in the function supplied to root" | 140 | ++ " components expected in the result of the function supplied to root" |
176 | 141 | ||
177 | checkdim2 n m | 142 | checkdim2 n m |
178 | | rows m == n && cols m == n = m | 143 | | rows m == n && cols m == n = m |
179 | | otherwise = error $ "Error: "++ show n ++ "x" ++ show n | 144 | | otherwise = error $ "Error: "++ show n ++ "x" ++ show n |
180 | ++ " Jacobian expected in root" | 145 | ++ " Jacobian expected in rootJ" |