From edb20f3c9993d9c9c349f6a001317592949154f1 Mon Sep 17 00:00:00 2001 From: Alex Lang Date: Wed, 8 May 2013 10:36:47 +0900 Subject: Implemented uniRoot for one-dimensional root-finding without derivatives --- lib/Numeric/GSL/Internal.hs | 3 +++ lib/Numeric/GSL/Root.hs | 44 +++++++++++++++++++++++++----- lib/Numeric/GSL/gsl-aux.c | 65 ++++++++++++++++++++++++++++++++++++++++++--- 3 files changed, 101 insertions(+), 11 deletions(-) (limited to 'lib') diff --git a/lib/Numeric/GSL/Internal.hs b/lib/Numeric/GSL/Internal.hs index 84417ce..69a9750 100644 --- a/lib/Numeric/GSL/Internal.hs +++ b/lib/Numeric/GSL/Internal.hs @@ -36,6 +36,9 @@ foreign import ccall safe "wrapper" foreign import ccall safe "wrapper" mkDoubleVecVecfun :: (Double -> TVV) -> IO (FunPtr (Double -> TVV)) +foreign import ccall safe "wrapper" + mkDoublefun :: (Double -> Double) -> IO (FunPtr (Double -> Double)) + aux_vTov :: (Vector Double -> Vector Double) -> TVV aux_vTov f n p nr r = g where v = f x 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 ------------------------------------------------------- diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c index 61b6a54..27ffbdc 100644 --- a/lib/Numeric/GSL/gsl-aux.c +++ b/lib/Numeric/GSL/gsl-aux.c @@ -32,6 +32,7 @@ #include #include #include +#include #include #include #include @@ -1047,9 +1048,65 @@ int minimizeD(int method, double f(int, double*), int df(int, double*, int, doub //--------------------------------------------------------------- +double only_f_aux_root(double x, void *pars) { + double (*f)(double) = (double (*)(double)) pars; + return f(x); +} + +int root(int method, double f(double), + double epsrel, int maxit, + double xl, double xu, RMAT(sol)) { + REQUIRES(solr == maxit && solc == 4,BAD_SIZE); + DEBUGMSG("root_only_f"); + gsl_function my_func; + // extract function from pars + my_func.function = only_f_aux_root; + my_func.params = f; + size_t iter = 0; + int status; + const gsl_root_fsolver_type *T; + gsl_root_fsolver *s; + // Starting point + switch(method) { + case 0 : {T = gsl_root_fsolver_bisection; printf("7\n"); break; } + case 1 : {T = gsl_root_fsolver_falsepos; break; } + case 2 : {T = gsl_root_fsolver_brent; break; } + default: ERROR(BAD_CODE); + } + s = gsl_root_fsolver_alloc (T); + gsl_root_fsolver_set (s, &my_func, xl, xu); + do { + double best, current_lo, current_hi; + status = gsl_root_fsolver_iterate (s); + best = gsl_root_fsolver_root (s); + current_lo = gsl_root_fsolver_x_lower (s); + current_hi = gsl_root_fsolver_x_upper (s); + solp[iter*solc] = iter + 1; + solp[iter*solc+1] = best; + solp[iter*solc+2] = current_lo; + solp[iter*solc+3] = current_hi; + iter++; + if (status) /* check if solver is stuck */ + break; + + status = + gsl_root_test_interval (current_lo, current_hi, 0, epsrel); + } + while (status == GSL_CONTINUE && iter < maxit); + int i; + for (i=iter; isize,sizeof(double)); double* q = (double*)calloc(y->size,sizeof(double)); @@ -1066,14 +1123,14 @@ int only_f_aux_root(const gsl_vector*x, void *pars, gsl_vector*y) { return 0; //hmmm } -int root(int method, void f(int, double*, int, double*), +int multiroot(int method, void f(int, double*, int, double*), double epsabs, int maxit, KRVEC(xi), RMAT(sol)) { REQUIRES(solr == maxit && solc == 1+2*xin,BAD_SIZE); DEBUGMSG("root_only_f"); gsl_multiroot_function my_func; // extract function from pars - my_func.f = only_f_aux_root; + my_func.f = only_f_aux_multiroot; my_func.n = xin; my_func.params = f; size_t iter = 0; @@ -1175,7 +1232,7 @@ int fjf_aux(const gsl_vector * x, void * pars, gsl_vector * f, gsl_matrix * g) { return 0; } -int rootj(int method, int f(int, double*, int, double*), +int multirootj(int method, int f(int, double*, int, double*), int jac(int, double*, int, int, double*), double epsabs, int maxit, KRVEC(xi), RMAT(sol)) { -- cgit v1.2.3