From edb20f3c9993d9c9c349f6a001317592949154f1 Mon Sep 17 00:00:00 2001
From: Alex Lang <me@alang.ca>
Date: Wed, 8 May 2013 10:36:47 +0900
Subject: Implemented uniRoot for one-dimensional root-finding without
 derivatives

---
 lib/Numeric/GSL/Root.hs | 44 +++++++++++++++++++++++++++++++++++++-------
 1 file changed, 37 insertions(+), 7 deletions(-)

(limited to 'lib/Numeric/GSL/Root.hs')

diff --git a/lib/Numeric/GSL/Root.hs b/lib/Numeric/GSL/Root.hs
index cd2982a..25ec5f5 100644
--- a/lib/Numeric/GSL/Root.hs
+++ b/lib/Numeric/GSL/Root.hs
@@ -45,6 +45,7 @@ main = do
 -----------------------------------------------------------------------------
 
 module Numeric.GSL.Root (
+    uniRoot, UniRootMethod(..),
     root, RootMethod(..),
     rootJ, RootMethodJ(..),
 ) where
@@ -58,6 +59,36 @@ import System.IO.Unsafe(unsafePerformIO)
 
 -------------------------------------------------------------------------
 
+data UniRootMethod = Bisection
+                   | FalsePos
+                   | Brent
+                   deriving (Enum, Eq, Show, Bounded)
+
+uniRoot :: UniRootMethod
+        -> Double
+        -> Int
+        -> (Double -> Double)
+        -> Double
+        -> Double
+        -> (Double, Matrix Double)
+uniRoot method epsrel maxit fun xl xu = uniRootGen (fi (fromEnum method)) fun xl xu epsrel maxit
+
+uniRootGen m f xl xu epsrel maxit = unsafePerformIO $ do
+    fp <- mkDoublefun f
+    rawpath <- createMIO maxit 4
+                         (c_root m fp epsrel (fi maxit) xl xu)
+                         "root"
+    let it = round (rawpath @@> (maxit-1,0))
+        path = takeRows it rawpath
+        [sol] = toLists $ dropRows (it-1) path
+    freeHaskellFunPtr fp
+    return (sol !! 1, path)
+
+foreign import ccall safe "root"
+    c_root:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> TM
+
+-------------------------------------------------------------------------
+
 data RootMethod = Hybrids
                 | Hybrid
                 | DNewton
@@ -82,8 +113,8 @@ rootGen m f xi epsabs maxit = unsafePerformIO $ do
     fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
     rawpath <- vec xiv $ \xiv' ->
                    createMIO maxit (2*n+1)
-                         (c_root m fp epsabs (fi maxit) // xiv')
-                         "root"
+                         (c_multiroot m fp epsabs (fi maxit) // xiv')
+                         "multiroot"
     let it = round (rawpath @@> (maxit-1,0))
         path = takeRows it rawpath
         [sol] = toLists $ dropRows (it-1) path
@@ -91,8 +122,8 @@ rootGen m f xi epsabs maxit = unsafePerformIO $ do
     return (take n $ drop 1 sol, path)
 
 
-foreign import ccall safe "root"
-    c_root:: CInt -> FunPtr TVV -> Double -> CInt -> TVM
+foreign import ccall safe "multiroot"
+    c_multiroot:: CInt -> FunPtr TVV -> Double -> CInt -> TVM
 
 -------------------------------------------------------------------------
 
@@ -121,7 +152,7 @@ rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do
     rawpath <- vec xiv $ \xiv' ->
                    createMIO maxit (2*n+1)
                          (c_rootj m fp jp epsabs (fi maxit) // xiv')
-                         "root"
+                         "multiroot"
     let it = round (rawpath @@> (maxit-1,0))
         path = takeRows it rawpath
         [sol] = toLists $ dropRows (it-1) path
@@ -129,8 +160,7 @@ rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do
     freeHaskellFunPtr jp
     return (take n $ drop 1 sol, path)
 
-
-foreign import ccall safe "rootj"
+foreign import ccall safe "multirootj"
     c_rootj:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> CInt -> TVM
 
 -------------------------------------------------------
-- 
cgit v1.2.3