    {-# LANGUAGE  FlexibleContexts #-}

    {-# OPTIONS_GHC -fno-warn-missing-signatures #-}

    {- |
    Module      :  Numeric.GSL.Root
    Copyright   :  (c) Alberto Ruiz 2009
    License     :  GPL
    Maintainer  :  Alberto Ruiz
    Stability   :  provisional

    Multidimensional root finding.

    <http://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Root_002dFinding.html>

    The example in the GSL manual:

    >>> let rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ]
    >>> let (sol,path) = root Hybrids 1E-7 30 (rosenbrock 1 10) [-10,-5]
    >>> sol
    [1.0,1.0]
    >>> disp 3 path
    11x5
     1.000  -10.000  -5.000  11.000  -1050.000
     2.000   -3.976  24.827   4.976     90.203
     3.000   -3.976  24.827   4.976     90.203
     4.000   -3.976  24.827   4.976     90.203
     5.000   -1.274  -5.680   2.274    -73.018
     6.000   -1.274  -5.680   2.274    -73.018
     7.000    0.249   0.298   0.751      2.359
     8.000    0.249   0.298   0.751      2.359
     9.000    1.000   0.878  -0.000     -1.218
    10.000    1.000   0.989  -0.000     -0.108
    11.000    1.000   1.000   0.000      0.000

    -}
    -----------------------------------------------------------------------------

    module Numeric.GSL.Root (
        uniRoot, UniRootMethod(..),
        uniRootJ, UniRootMethodJ(..),
        root, RootMethod(..),
        rootJ, RootMethodJ(..),
    ) where

    import Numeric.LinearAlgebra.HMatrix
    import Numeric.GSL.Internal
    import Foreign.Ptr(FunPtr, freeHaskellFunPtr)
    import Foreign.C.Types
    import System.IO.Unsafe(unsafePerformIO)

    -------------------------------------------------------------------------
    type TVV = TV (TV Res)
    type TVM = TV (TM Res)

    data UniRootMethod = Bisection
                       | FalsePos
                       | Brent
                       deriving (Int -> UniRootMethod
    UniRootMethod -> Int
    UniRootMethod -> [UniRootMethod]
    UniRootMethod -> UniRootMethod
    UniRootMethod -> UniRootMethod -> [UniRootMethod]
    UniRootMethod -> UniRootMethod -> UniRootMethod -> [UniRootMethod]
    (UniRootMethod -> UniRootMethod)
    -> (UniRootMethod -> UniRootMethod)
    -> (Int -> UniRootMethod)
    -> (UniRootMethod -> Int)
    -> (UniRootMethod -> [UniRootMethod])
    -> (UniRootMethod -> UniRootMethod -> [UniRootMethod])
    -> (UniRootMethod -> UniRootMethod -> [UniRootMethod])
    -> (UniRootMethod
        -> UniRootMethod -> UniRootMethod -> [UniRootMethod])
    -> Enum UniRootMethod
    forall a.
    (a -> a)
    -> (a -> a)
    -> (Int -> a)
    -> (a -> Int)
    -> (a -> [a])
    -> (a -> a -> [a])
    -> (a -> a -> [a])
    -> (a -> a -> a -> [a])
    -> Enum a
    $csucc :: UniRootMethod -> UniRootMethod
    succ :: UniRootMethod -> UniRootMethod
    $cpred :: UniRootMethod -> UniRootMethod
    pred :: UniRootMethod -> UniRootMethod
    $ctoEnum :: Int -> UniRootMethod
    toEnum :: Int -> UniRootMethod
    $cfromEnum :: UniRootMethod -> Int
    fromEnum :: UniRootMethod -> Int
    $cenumFrom :: UniRootMethod -> [UniRootMethod]
    enumFrom :: UniRootMethod -> [UniRootMethod]
    $cenumFromThen :: UniRootMethod -> UniRootMethod -> [UniRootMethod]
    enumFromThen :: UniRootMethod -> UniRootMethod -> [UniRootMethod]
    $cenumFromTo :: UniRootMethod -> UniRootMethod -> [UniRootMethod]
    enumFromTo :: UniRootMethod -> UniRootMethod -> [UniRootMethod]
    $cenumFromThenTo :: UniRootMethod -> UniRootMethod -> UniRootMethod -> [UniRootMethod]
    enumFromThenTo :: UniRootMethod -> UniRootMethod -> UniRootMethod -> [UniRootMethod]
    Enum, UniRootMethod -> UniRootMethod -> Bool
    (UniRootMethod -> UniRootMethod -> Bool)
    -> (UniRootMethod -> UniRootMethod -> Bool) -> Eq UniRootMethod
    forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
    $c== :: UniRootMethod -> UniRootMethod -> Bool
    == :: UniRootMethod -> UniRootMethod -> Bool
    $c/= :: UniRootMethod -> UniRootMethod -> Bool
    /= :: UniRootMethod -> UniRootMethod -> Bool
    Eq, Int -> UniRootMethod -> ShowS
    [UniRootMethod] -> ShowS
    UniRootMethod -> String
    (Int -> UniRootMethod -> ShowS)
    -> (UniRootMethod -> String)
    -> ([UniRootMethod] -> ShowS)
    -> Show UniRootMethod
    forall a.
    (Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
    $cshowsPrec :: Int -> UniRootMethod -> ShowS
    showsPrec :: Int -> UniRootMethod -> ShowS
    $cshow :: UniRootMethod -> String
    show :: UniRootMethod -> String
    $cshowList :: [UniRootMethod] -> ShowS
    showList :: [UniRootMethod] -> ShowS
    Show, UniRootMethod
    UniRootMethod -> UniRootMethod -> Bounded UniRootMethod
    forall a. a -> a -> Bounded a
    $cminBound :: UniRootMethod
    minBound :: UniRootMethod
    $cmaxBound :: UniRootMethod
    maxBound :: UniRootMethod
    Bounded)

    uniRoot :: UniRootMethod
            -> Double
            -> Int
            -> (Double -> Double)
            -> Double
            -> Double
            -> (Double, Matrix Double)
    uniRoot :: UniRootMethod
    -> Double
    -> Int
    -> (Double -> Double)
    -> Double
    -> Double
    -> (Double, Matrix Double)
    uniRoot UniRootMethod
    method Double
    epsrel Int
    maxit Double -> Double
    fun Double
    xl Double
    xu = CInt
    -> (Double -> Double)
    -> Double
    -> Double
    -> Double
    -> Int
    -> (Double, Matrix Double)
    uniRootGen (Int -> CInt
    fi (UniRootMethod -> Int
    forall a. Enum a => a -> Int
    fromEnum UniRootMethod
    method)) Double -> Double
    fun Double
    xl Double
    xu Double
    epsrel Int
    maxit

    uniRootGen :: CInt
    -> (Double -> Double)
    -> Double
    -> Double
    -> Double
    -> Int
    -> (Double, Matrix Double)
    uniRootGen CInt
    m Double -> Double
    f Double
    xl Double
    xu Double
    epsrel Int
    maxit = IO (Double, Matrix Double) -> (Double, Matrix Double)
    forall a. IO a -> a
    unsafePerformIO (IO (Double, Matrix Double) -> (Double, Matrix Double))
    -> IO (Double, Matrix Double) -> (Double, Matrix Double)
    forall a b. (a -> b) -> a -> b
    $ do
        fp <- (Double -> Double) -> IO (FunPtr (Double -> Double))
    mkDoublefun Double -> Double
    f
        rawpath <- createMIO maxit 4
                             (c_root m fp epsrel (fi maxit) xl xu)
                             "root"
        let it = Double -> Int
    forall b. Integral b => Double -> b
    forall a b. (RealFrac a, Integral b) => a -> b
    round (Matrix Double
    rawpath Matrix Double -> IndexOf Matrix -> Double
    forall (c :: * -> *) e. Container c e => c e -> IndexOf c -> e
    `atIndex` (Int
    maxitInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1,Int
    0))
            path = Int -> Matrix Double -> Matrix Double
    forall t. Element t => Int -> Matrix t -> Matrix t
    takeRows Int
    it Matrix Double
    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 Res

    -------------------------------------------------------------------------
    data UniRootMethodJ = UNewton
                        | Secant
                        | Steffenson
                        deriving (Int -> UniRootMethodJ
    UniRootMethodJ -> Int
    UniRootMethodJ -> [UniRootMethodJ]
    UniRootMethodJ -> UniRootMethodJ
    UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ]
    UniRootMethodJ
    -> UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ]
    (UniRootMethodJ -> UniRootMethodJ)
    -> (UniRootMethodJ -> UniRootMethodJ)
    -> (Int -> UniRootMethodJ)
    -> (UniRootMethodJ -> Int)
    -> (UniRootMethodJ -> [UniRootMethodJ])
    -> (UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ])
    -> (UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ])
    -> (UniRootMethodJ
        -> UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ])
    -> Enum UniRootMethodJ
    forall a.
    (a -> a)
    -> (a -> a)
    -> (Int -> a)
    -> (a -> Int)
    -> (a -> [a])
    -> (a -> a -> [a])
    -> (a -> a -> [a])
    -> (a -> a -> a -> [a])
    -> Enum a
    $csucc :: UniRootMethodJ -> UniRootMethodJ
    succ :: UniRootMethodJ -> UniRootMethodJ
    $cpred :: UniRootMethodJ -> UniRootMethodJ
    pred :: UniRootMethodJ -> UniRootMethodJ
    $ctoEnum :: Int -> UniRootMethodJ
    toEnum :: Int -> UniRootMethodJ
    $cfromEnum :: UniRootMethodJ -> Int
    fromEnum :: UniRootMethodJ -> Int
    $cenumFrom :: UniRootMethodJ -> [UniRootMethodJ]
    enumFrom :: UniRootMethodJ -> [UniRootMethodJ]
    $cenumFromThen :: UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ]
    enumFromThen :: UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ]
    $cenumFromTo :: UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ]
    enumFromTo :: UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ]
    $cenumFromThenTo :: UniRootMethodJ
    -> UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ]
    enumFromThenTo :: UniRootMethodJ
    -> UniRootMethodJ -> UniRootMethodJ -> [UniRootMethodJ]
    Enum, UniRootMethodJ -> UniRootMethodJ -> Bool
    (UniRootMethodJ -> UniRootMethodJ -> Bool)
    -> (UniRootMethodJ -> UniRootMethodJ -> Bool) -> Eq UniRootMethodJ
    forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
    $c== :: UniRootMethodJ -> UniRootMethodJ -> Bool
    == :: UniRootMethodJ -> UniRootMethodJ -> Bool
    $c/= :: UniRootMethodJ -> UniRootMethodJ -> Bool
    /= :: UniRootMethodJ -> UniRootMethodJ -> Bool
    Eq, Int -> UniRootMethodJ -> ShowS
    [UniRootMethodJ] -> ShowS
    UniRootMethodJ -> String
    (Int -> UniRootMethodJ -> ShowS)
    -> (UniRootMethodJ -> String)
    -> ([UniRootMethodJ] -> ShowS)
    -> Show UniRootMethodJ
    forall a.
    (Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
    $cshowsPrec :: Int -> UniRootMethodJ -> ShowS
    showsPrec :: Int -> UniRootMethodJ -> ShowS
    $cshow :: UniRootMethodJ -> String
    show :: UniRootMethodJ -> String
    $cshowList :: [UniRootMethodJ] -> ShowS
    showList :: [UniRootMethodJ] -> ShowS
    Show, UniRootMethodJ
    UniRootMethodJ -> UniRootMethodJ -> Bounded UniRootMethodJ
    forall a. a -> a -> Bounded a
    $cminBound :: UniRootMethodJ
    minBound :: UniRootMethodJ
    $cmaxBound :: UniRootMethodJ
    maxBound :: UniRootMethodJ
    Bounded)

    uniRootJ :: UniRootMethodJ
            -> Double
            -> Int
            -> (Double -> Double)
            -> (Double -> Double)
            -> Double
            -> (Double, Matrix Double)
    uniRootJ :: UniRootMethodJ
    -> Double
    -> Int
    -> (Double -> Double)
    -> (Double -> Double)
    -> Double
    -> (Double, Matrix Double)
    uniRootJ UniRootMethodJ
    method Double
    epsrel Int
    maxit Double -> Double
    fun Double -> Double
    dfun Double
    x = CInt
    -> (Double -> Double)
    -> (Double -> Double)
    -> Double
    -> Double
    -> Int
    -> (Double, Matrix Double)
    uniRootJGen (Int -> CInt
    fi (UniRootMethodJ -> Int
    forall a. Enum a => a -> Int
    fromEnum UniRootMethodJ
    method)) Double -> Double
    fun
        Double -> Double
    dfun Double
    x Double
    epsrel Int
    maxit

    uniRootJGen :: CInt
    -> (Double -> Double)
    -> (Double -> Double)
    -> Double
    -> Double
    -> Int
    -> (Double, Matrix Double)
    uniRootJGen CInt
    m Double -> Double
    f Double -> Double
    df Double
    x Double
    epsrel Int
    maxit = IO (Double, Matrix Double) -> (Double, Matrix Double)
    forall a. IO a -> a
    unsafePerformIO (IO (Double, Matrix Double) -> (Double, Matrix Double))
    -> IO (Double, Matrix Double) -> (Double, Matrix Double)
    forall a b. (a -> b) -> a -> b
    $ do
        fp <- (Double -> Double) -> IO (FunPtr (Double -> Double))
    mkDoublefun Double -> Double
    f
        dfp <- mkDoublefun df
        rawpath <- createMIO maxit 2
                             (c_rootj m fp dfp epsrel (fi maxit) x)
                             "rootj"
        let it = Double -> Int
    forall b. Integral b => Double -> b
    forall a b. (RealFrac a, Integral b) => a -> b
    round (Matrix Double
    rawpath Matrix Double -> IndexOf Matrix -> Double
    forall (c :: * -> *) e. Container c e => c e -> IndexOf c -> e
    `atIndex` (Int
    maxitInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1,Int
    0))
            path = Int -> Matrix Double -> Matrix Double
    forall t. Element t => Int -> Matrix t -> Matrix t
    takeRows Int
    it Matrix Double
    rawpath
            [sol] = toLists $ dropRows (it-1) path
        freeHaskellFunPtr fp
        return (sol !! 1, path)

    foreign import ccall safe "rootj"
        c_rootj :: CInt -> FunPtr (Double -> Double) -> FunPtr (Double -> Double)
                -> Double -> CInt -> Double -> TM Res

    -------------------------------------------------------------------------

    data RootMethod = Hybrids
                    | Hybrid
                    | DNewton
                    | Broyden
                    deriving (Int -> RootMethod
    RootMethod -> Int
    RootMethod -> [RootMethod]
    RootMethod -> RootMethod
    RootMethod -> RootMethod -> [RootMethod]
    RootMethod -> RootMethod -> RootMethod -> [RootMethod]
    (RootMethod -> RootMethod)
    -> (RootMethod -> RootMethod)
    -> (Int -> RootMethod)
    -> (RootMethod -> Int)
    -> (RootMethod -> [RootMethod])
    -> (RootMethod -> RootMethod -> [RootMethod])
    -> (RootMethod -> RootMethod -> [RootMethod])
    -> (RootMethod -> RootMethod -> RootMethod -> [RootMethod])
    -> Enum RootMethod
    forall a.
    (a -> a)
    -> (a -> a)
    -> (Int -> a)
    -> (a -> Int)
    -> (a -> [a])
    -> (a -> a -> [a])
    -> (a -> a -> [a])
    -> (a -> a -> a -> [a])
    -> Enum a
    $csucc :: RootMethod -> RootMethod
    succ :: RootMethod -> RootMethod
    $cpred :: RootMethod -> RootMethod
    pred :: RootMethod -> RootMethod
    $ctoEnum :: Int -> RootMethod
    toEnum :: Int -> RootMethod
    $cfromEnum :: RootMethod -> Int
    fromEnum :: RootMethod -> Int
    $cenumFrom :: RootMethod -> [RootMethod]
    enumFrom :: RootMethod -> [RootMethod]
    $cenumFromThen :: RootMethod -> RootMethod -> [RootMethod]
    enumFromThen :: RootMethod -> RootMethod -> [RootMethod]
    $cenumFromTo :: RootMethod -> RootMethod -> [RootMethod]
    enumFromTo :: RootMethod -> RootMethod -> [RootMethod]
    $cenumFromThenTo :: RootMethod -> RootMethod -> RootMethod -> [RootMethod]
    enumFromThenTo :: RootMethod -> RootMethod -> RootMethod -> [RootMethod]
    Enum,RootMethod -> RootMethod -> Bool
    (RootMethod -> RootMethod -> Bool)
    -> (RootMethod -> RootMethod -> Bool) -> Eq RootMethod
    forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
    $c== :: RootMethod -> RootMethod -> Bool
    == :: RootMethod -> RootMethod -> Bool
    $c/= :: RootMethod -> RootMethod -> Bool
    /= :: RootMethod -> RootMethod -> Bool
    Eq,Int -> RootMethod -> ShowS
    [RootMethod] -> ShowS
    RootMethod -> String
    (Int -> RootMethod -> ShowS)
    -> (RootMethod -> String)
    -> ([RootMethod] -> ShowS)
    -> Show RootMethod
    forall a.
    (Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
    $cshowsPrec :: Int -> RootMethod -> ShowS
    showsPrec :: Int -> RootMethod -> ShowS
    $cshow :: RootMethod -> String
    show :: RootMethod -> String
    $cshowList :: [RootMethod] -> ShowS
    showList :: [RootMethod] -> ShowS
    Show,RootMethod
    RootMethod -> RootMethod -> Bounded RootMethod
    forall a. a -> a -> Bounded a
    $cminBound :: RootMethod
    minBound :: RootMethod
    $cmaxBound :: RootMethod
    maxBound :: RootMethod
    Bounded)

    -- | 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 :: RootMethod
    -> Double
    -> Int
    -> ([Double] -> [Double])
    -> [Double]
    -> ([Double], Matrix Double)
    root RootMethod
    method Double
    epsabs Int
    maxit [Double] -> [Double]
    fun [Double]
    xinit = CInt
    -> ([Double] -> [Double])
    -> [Double]
    -> Double
    -> Int
    -> ([Double], Matrix Double)
    rootGen (Int -> CInt
    fi (RootMethod -> Int
    forall a. Enum a => a -> Int
    fromEnum RootMethod
    method)) [Double] -> [Double]
    fun [Double]
    xinit Double
    epsabs Int
    maxit

    rootGen :: CInt
    -> ([Double] -> [Double])
    -> [Double]
    -> Double
    -> Int
    -> ([Double], Matrix Double)
    rootGen CInt
    m [Double] -> [Double]
    f [Double]
    xi Double
    epsabs Int
    maxit = IO ([Double], Matrix Double) -> ([Double], Matrix Double)
    forall a. IO a -> a
    unsafePerformIO (IO ([Double], Matrix Double) -> ([Double], Matrix Double))
    -> IO ([Double], Matrix Double) -> ([Double], Matrix Double)
    forall a b. (a -> b) -> a -> b
    $ do
        let xiv :: Vector Double
    xiv = [Double] -> Vector Double
    forall a. Storable a => [a] -> Vector a
    fromList [Double]
    xi
            n :: IndexOf Vector
    n   = Vector Double -> IndexOf Vector
    forall (c :: * -> *) t. Container c t => c t -> IndexOf c
    size Vector Double
    xiv
        fp <- TVV -> IO (FunPtr TVV)
    mkVecVecfun ((Vector Double -> Vector Double) -> TVV
    aux_vTov (IndexOf Vector -> Vector Double -> Vector Double
    forall {c :: * -> *} {t}.
    (Eq (IndexOf c), Container c t, Show (IndexOf c)) =>
    IndexOf c -> c t -> c t
    checkdim1 Int
    IndexOf Vector
    n (Vector Double -> Vector Double)
    -> (Vector Double -> Vector Double)
    -> Vector Double
    -> Vector Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . [Double] -> Vector Double
    forall a. Storable a => [a] -> Vector a
    fromList ([Double] -> Vector Double)
    -> (Vector Double -> [Double]) -> Vector Double -> Vector Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . [Double] -> [Double]
    f ([Double] -> [Double])
    -> (Vector Double -> [Double]) -> Vector Double -> [Double]
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector Double -> [Double]
    forall a. Storable a => Vector a -> [a]
    toList))
        rawpath <- vec xiv $ \(CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
    -> CInt -> CInt -> Ptr Double -> IO CInt
    xiv' ->
                       Int
    -> Int
    -> (CInt -> CInt -> Ptr Double -> IO CInt)
    -> String
    -> IO (Matrix Double)
    forall {a}.
    Storable a =>
    Int
    -> Int
    -> (CInt -> CInt -> Ptr a -> IO CInt)
    -> String
    -> IO (Matrix a)
    createMIO Int
    maxit (Int
    2Int -> Int -> Int
    forall a. Num a => a -> a -> a
    *Int
    nInt -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    1)
                             (CInt
    -> FunPtr TVV
    -> Double
    -> CInt
    -> CInt
    -> Ptr Double
    -> CInt
    -> CInt
    -> Ptr Double
    -> IO CInt
    c_multiroot CInt
    m FunPtr TVV
    fp Double
    epsabs (Int -> CInt
    fi Int
    maxit) (CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
    -> ((CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
        -> CInt -> CInt -> Ptr Double -> IO CInt)
    -> CInt
    -> CInt
    -> Ptr Double
    -> IO CInt
    forall x y. x -> (x -> y) -> y
    // (CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
    -> CInt -> CInt -> Ptr Double -> IO CInt
    xiv')
                             String
    "multiroot"
        let it = Double -> Int
    forall b. Integral b => Double -> b
    forall a b. (RealFrac a, Integral b) => a -> b
    round (Matrix Double
    rawpath Matrix Double -> IndexOf Matrix -> Double
    forall (c :: * -> *) e. Container c e => c e -> IndexOf c -> e
    `atIndex` (Int
    maxitInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1,Int
    0))
            path = Int -> Matrix Double -> Matrix Double
    forall t. Element t => Int -> Matrix t -> Matrix t
    takeRows Int
    it Matrix Double
    rawpath
            [sol] = toLists $ dropRows (it-1) path
        freeHaskellFunPtr fp
        return (take n $ drop 1 sol, path)

    foreign import ccall safe "multiroot"
        c_multiroot:: CInt -> FunPtr TVV -> Double -> CInt -> TVM

    -------------------------------------------------------------------------

    data RootMethodJ = HybridsJ
                     | HybridJ
                     | Newton
                     | GNewton
                    deriving (Int -> RootMethodJ
    RootMethodJ -> Int
    RootMethodJ -> [RootMethodJ]
    RootMethodJ -> RootMethodJ
    RootMethodJ -> RootMethodJ -> [RootMethodJ]
    RootMethodJ -> RootMethodJ -> RootMethodJ -> [RootMethodJ]
    (RootMethodJ -> RootMethodJ)
    -> (RootMethodJ -> RootMethodJ)
    -> (Int -> RootMethodJ)
    -> (RootMethodJ -> Int)
    -> (RootMethodJ -> [RootMethodJ])
    -> (RootMethodJ -> RootMethodJ -> [RootMethodJ])
    -> (RootMethodJ -> RootMethodJ -> [RootMethodJ])
    -> (RootMethodJ -> RootMethodJ -> RootMethodJ -> [RootMethodJ])
    -> Enum RootMethodJ
    forall a.
    (a -> a)
    -> (a -> a)
    -> (Int -> a)
    -> (a -> Int)
    -> (a -> [a])
    -> (a -> a -> [a])
    -> (a -> a -> [a])
    -> (a -> a -> a -> [a])
    -> Enum a
    $csucc :: RootMethodJ -> RootMethodJ
    succ :: RootMethodJ -> RootMethodJ
    $cpred :: RootMethodJ -> RootMethodJ
    pred :: RootMethodJ -> RootMethodJ
    $ctoEnum :: Int -> RootMethodJ
    toEnum :: Int -> RootMethodJ
    $cfromEnum :: RootMethodJ -> Int
    fromEnum :: RootMethodJ -> Int
    $cenumFrom :: RootMethodJ -> [RootMethodJ]
    enumFrom :: RootMethodJ -> [RootMethodJ]
    $cenumFromThen :: RootMethodJ -> RootMethodJ -> [RootMethodJ]
    enumFromThen :: RootMethodJ -> RootMethodJ -> [RootMethodJ]
    $cenumFromTo :: RootMethodJ -> RootMethodJ -> [RootMethodJ]
    enumFromTo :: RootMethodJ -> RootMethodJ -> [RootMethodJ]
    $cenumFromThenTo :: RootMethodJ -> RootMethodJ -> RootMethodJ -> [RootMethodJ]
    enumFromThenTo :: RootMethodJ -> RootMethodJ -> RootMethodJ -> [RootMethodJ]
    Enum,RootMethodJ -> RootMethodJ -> Bool
    (RootMethodJ -> RootMethodJ -> Bool)
    -> (RootMethodJ -> RootMethodJ -> Bool) -> Eq RootMethodJ
    forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
    $c== :: RootMethodJ -> RootMethodJ -> Bool
    == :: RootMethodJ -> RootMethodJ -> Bool
    $c/= :: RootMethodJ -> RootMethodJ -> Bool
    /= :: RootMethodJ -> RootMethodJ -> Bool
    Eq,Int -> RootMethodJ -> ShowS
    [RootMethodJ] -> ShowS
    RootMethodJ -> String
    (Int -> RootMethodJ -> ShowS)
    -> (RootMethodJ -> String)
    -> ([RootMethodJ] -> ShowS)
    -> Show RootMethodJ
    forall a.
    (Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
    $cshowsPrec :: Int -> RootMethodJ -> ShowS
    showsPrec :: Int -> RootMethodJ -> ShowS
    $cshow :: RootMethodJ -> String
    show :: RootMethodJ -> String
    $cshowList :: [RootMethodJ] -> ShowS
    showList :: [RootMethodJ] -> ShowS
    Show,RootMethodJ
    RootMethodJ -> RootMethodJ -> Bounded RootMethodJ
    forall a. a -> a -> Bounded a
    $cminBound :: RootMethodJ
    minBound :: RootMethodJ
    $cmaxBound :: RootMethodJ
    maxBound :: RootMethodJ
    Bounded)

    -- | Nonlinear multidimensional root finding using both the function and its derivatives.
    rootJ :: RootMethodJ
          -> Double                     -- ^ maximum residual
          -> Int                        -- ^ maximum number of iterations allowed
          -> ([Double] -> [Double])     -- ^ function to minimize
          -> ([Double] -> [[Double]])   -- ^ Jacobian
          -> [Double]                   -- ^ starting point
          -> ([Double], Matrix Double)  -- ^ solution vector and optimization path

    rootJ :: RootMethodJ
    -> Double
    -> Int
    -> ([Double] -> [Double])
    -> ([Double] -> [[Double]])
    -> [Double]
    -> ([Double], Matrix Double)
    rootJ RootMethodJ
    method Double
    epsabs Int
    maxit [Double] -> [Double]
    fun [Double] -> [[Double]]
    jac [Double]
    xinit = CInt
    -> ([Double] -> [Double])
    -> ([Double] -> [[Double]])
    -> [Double]
    -> Double
    -> Int
    -> ([Double], Matrix Double)
    rootJGen (Int -> CInt
    fi (RootMethodJ -> Int
    forall a. Enum a => a -> Int
    fromEnum RootMethodJ
    method)) [Double] -> [Double]
    fun [Double] -> [[Double]]
    jac [Double]
    xinit Double
    epsabs Int
    maxit

    rootJGen :: CInt
    -> ([Double] -> [Double])
    -> ([Double] -> [[Double]])
    -> [Double]
    -> Double
    -> Int
    -> ([Double], Matrix Double)
    rootJGen CInt
    m [Double] -> [Double]
    f [Double] -> [[Double]]
    jac [Double]
    xi Double
    epsabs Int
    maxit = IO ([Double], Matrix Double) -> ([Double], Matrix Double)
    forall a. IO a -> a
    unsafePerformIO (IO ([Double], Matrix Double) -> ([Double], Matrix Double))
    -> IO ([Double], Matrix Double) -> ([Double], Matrix Double)
    forall a b. (a -> b) -> a -> b
    $ do
        let xiv :: Vector Double
    xiv = [Double] -> Vector Double
    forall a. Storable a => [a] -> Vector a
    fromList [Double]
    xi
            n :: IndexOf Vector
    n   = Vector Double -> IndexOf Vector
    forall (c :: * -> *) t. Container c t => c t -> IndexOf c
    size Vector Double
    xiv
        fp <- TVV -> IO (FunPtr TVV)
    mkVecVecfun ((Vector Double -> Vector Double) -> TVV
    aux_vTov (IndexOf Vector -> Vector Double -> Vector Double
    forall {c :: * -> *} {t}.
    (Eq (IndexOf c), Container c t, Show (IndexOf c)) =>
    IndexOf c -> c t -> c t
    checkdim1 Int
    IndexOf Vector
    n (Vector Double -> Vector Double)
    -> (Vector Double -> Vector Double)
    -> Vector Double
    -> Vector Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . [Double] -> Vector Double
    forall a. Storable a => [a] -> Vector a
    fromList ([Double] -> Vector Double)
    -> (Vector Double -> [Double]) -> Vector Double -> Vector Double
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . [Double] -> [Double]
    f ([Double] -> [Double])
    -> (Vector Double -> [Double]) -> Vector Double -> [Double]
    forall b c a. (b -> c) -> (a -> b) -> a -> c
    . Vector Double -> [Double]
    forall a. Storable a => Vector a -> [a]
    toList))
        jp <- mkVecMatfun (aux_vTom (checkdim2 n . fromLists . jac . toList))
        rawpath <- vec xiv $ \(CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
    -> CInt -> CInt -> Ptr Double -> IO CInt
    xiv' ->
                       Int
    -> Int
    -> (CInt -> CInt -> Ptr Double -> IO CInt)
    -> String
    -> IO (Matrix Double)
    forall {a}.
    Storable a =>
    Int
    -> Int
    -> (CInt -> CInt -> Ptr a -> IO CInt)
    -> String
    -> IO (Matrix a)
    createMIO Int
    maxit (Int
    2Int -> Int -> Int
    forall a. Num a => a -> a -> a
    *Int
    nInt -> Int -> Int
    forall a. Num a => a -> a -> a
    +Int
    1)
                             (CInt
    -> FunPtr TVV
    -> FunPtr
         (CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
    -> Double
    -> CInt
    -> CInt
    -> Ptr Double
    -> CInt
    -> CInt
    -> Ptr Double
    -> IO CInt
    c_multirootj CInt
    m FunPtr TVV
    fp FunPtr
      (CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
    jp Double
    epsabs (Int -> CInt
    fi Int
    maxit) (CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
    -> ((CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
        -> CInt -> CInt -> Ptr Double -> IO CInt)
    -> CInt
    -> CInt
    -> Ptr Double
    -> IO CInt
    forall x y. x -> (x -> y) -> y
    // (CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt)
    -> CInt -> CInt -> Ptr Double -> IO CInt
    xiv')
                             String
    "multiroot"
        let it = Double -> Int
    forall b. Integral b => Double -> b
    forall a b. (RealFrac a, Integral b) => a -> b
    round (Matrix Double
    rawpath Matrix Double -> IndexOf Matrix -> Double
    forall (c :: * -> *) e. Container c e => c e -> IndexOf c -> e
    `atIndex` (Int
    maxitInt -> Int -> Int
    forall a. Num a => a -> a -> a
    -Int
    1,Int
    0))
            path = Int -> Matrix Double -> Matrix Double
    forall t. Element t => Int -> Matrix t -> Matrix t
    takeRows Int
    it Matrix Double
    rawpath
            [sol] = toLists $ dropRows (it-1) path
        freeHaskellFunPtr fp
        freeHaskellFunPtr jp
        return (take n $ drop 1 sol, path)

    foreign import ccall safe "multirootj"
        c_multirootj:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> CInt -> TVM

    -------------------------------------------------------

    checkdim1 :: IndexOf c -> c t -> c t
    checkdim1 IndexOf c
    n c t
    v
        | c t -> IndexOf c
    forall (c :: * -> *) t. Container c t => c t -> IndexOf c
    size c t
    v IndexOf c -> IndexOf c -> Bool
    forall a. Eq a => a -> a -> Bool
    == IndexOf c
    n = c t
    v
        | Bool
    otherwise = String -> c t
    forall a. HasCallStack => String -> a
    error (String -> c t) -> String -> c t
    forall a b. (a -> b) -> a -> b
    $ String
    "Error: "String -> ShowS
    forall a. [a] -> [a] -> [a]
    ++ IndexOf c -> String
    forall a. Show a => a -> String
    show IndexOf c
    n
                            String -> ShowS
    forall a. [a] -> [a] -> [a]
    ++ String
    " components expected in the result of the function supplied to root"

    checkdim2 :: Int -> Matrix t -> Matrix t
    checkdim2 Int
    n Matrix t
    m
        | Matrix t -> Int
    forall t. Matrix t -> Int
    rows Matrix t
    m Int -> Int -> Bool
    forall a. Eq a => a -> a -> Bool
    == Int
    n Bool -> Bool -> Bool
    && Matrix t -> Int
    forall t. Matrix t -> Int
    cols Matrix t
    m Int -> Int -> Bool
    forall a. Eq a => a -> a -> Bool
    == Int
    n = Matrix t
    m
        | Bool
    otherwise = String -> Matrix t
    forall a. HasCallStack => String -> a
    error (String -> Matrix t) -> String -> Matrix t
    forall a b. (a -> b) -> a -> b
    $ String
    "Error: "String -> ShowS
    forall a. [a] -> [a] -> [a]
    ++ Int -> String
    forall a. Show a => a -> String
    show Int
    n String -> ShowS
    forall a. [a] -> [a] -> [a]
    ++ String
    "x" String -> ShowS
    forall a. [a] -> [a] -> [a]
    ++ Int -> String
    forall a. Show a => a -> String
    show Int
    n
                            String -> ShowS
    forall a. [a] -> [a] -> [a]
    ++ String
    " Jacobian expected in rootJ"