module Numeric.GSL.Root (
root, RootMethod(..),
rootJ, RootMethodJ(..),
) where
import Data.Packed.Internal
import Data.Packed.Matrix
import Foreign
import Foreign.C.Types(CInt)
import Numeric.GSL.Internal
data RootMethod = Hybrids
| Hybrid
| DNewton
| Broyden
deriving (Enum,Eq,Show,Bounded)
root :: RootMethod
-> Double
-> Int
-> ([Double] -> [Double])
-> [Double]
-> ([Double], Matrix Double)
root method epsabs maxit fun xinit = rootGen (fi (fromEnum method)) fun xinit epsabs maxit
rootGen m f xi epsabs maxit = unsafePerformIO $ do
let xiv = fromList xi
n = dim xiv
fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
rawpath <- withVector xiv $ \xiv' ->
createMIO maxit (2*n+1)
(c_root m fp epsabs (fi maxit) // xiv')
"root"
let it = round (rawpath @@> (maxit1,0))
path = takeRows it rawpath
[sol] = toLists $ dropRows (it1) path
freeHaskellFunPtr fp
return (take n $ drop 1 sol, path)
foreign import ccall "root"
c_root:: CInt -> FunPtr TVV -> Double -> CInt -> TVM
data RootMethodJ = HybridsJ
| HybridJ
| Newton
| GNewton
deriving (Enum,Eq,Show,Bounded)
rootJ :: RootMethodJ
-> Double
-> Int
-> ([Double] -> [Double])
-> ([Double] -> [[Double]])
-> [Double]
-> ([Double], Matrix Double)
rootJ method epsabs maxit fun jac xinit = rootJGen (fi (fromEnum method)) fun jac xinit epsabs maxit
rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do
let xiv = fromList xi
n = dim xiv
fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
jp <- mkVecMatfun (aux_vTom (checkdim2 n . fromLists . jac . toList))
rawpath <- withVector xiv $ \xiv' ->
createMIO maxit (2*n+1)
(c_rootj m fp jp epsabs (fi maxit) // xiv')
"root"
let it = round (rawpath @@> (maxit1,0))
path = takeRows it rawpath
[sol] = toLists $ dropRows (it1) path
freeHaskellFunPtr fp
freeHaskellFunPtr jp
return (take n $ drop 1 sol, path)
foreign import ccall "rootj"
c_rootj:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> CInt -> TVM
checkdim1 n v
| dim v == n = v
| otherwise = error $ "Error: "++ show n
++ " components expected in the result of the function supplied to root"
checkdim2 n m
| rows m == n && cols m == n = m
| otherwise = error $ "Error: "++ show n ++ "x" ++ show n
++ " Jacobian expected in rootJ"