From 551cf7498c33bc0948bb4cb8444ae6f8af7278ea Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Thu, 8 May 2014 13:43:07 +0200 Subject: separation ok --- packages/hmatrix/hmatrix.cabal | 37 ++---- packages/hmatrix/src/Numeric/Container.hs | 14 +-- packages/hmatrix/src/Numeric/ContainerBoot.hs | 12 +- packages/hmatrix/src/Numeric/GSL.hs | 8 +- .../hmatrix/src/Numeric/GSL/Differentiation.hs | 12 +- packages/hmatrix/src/Numeric/GSL/Fitting.hs | 9 +- packages/hmatrix/src/Numeric/GSL/Fourier.hs | 15 +-- packages/hmatrix/src/Numeric/GSL/Integration.hs | 12 +- packages/hmatrix/src/Numeric/GSL/Internal.hs | 62 +++++++++- packages/hmatrix/src/Numeric/GSL/Minimization.hs | 25 ++-- packages/hmatrix/src/Numeric/GSL/ODE.hs | 10 +- packages/hmatrix/src/Numeric/GSL/Polynomials.hs | 17 ++- packages/hmatrix/src/Numeric/GSL/Root.hs | 14 +-- packages/hmatrix/src/Numeric/GSL/Vector.hs | 133 +++++++++++++++++++-- packages/hmatrix/src/Numeric/IO.hs | 3 +- .../src/Numeric/LinearAlgebra/Algorithms.hs | 4 +- 16 files changed, 256 insertions(+), 131 deletions(-) (limited to 'packages/hmatrix') diff --git a/packages/hmatrix/hmatrix.cabal b/packages/hmatrix/hmatrix.cabal index afd678f..8496663 100644 --- a/packages/hmatrix/hmatrix.cabal +++ b/packages/hmatrix/hmatrix.cabal @@ -6,11 +6,11 @@ Author: Alberto Ruiz Maintainer: Alberto Ruiz Stability: provisional Homepage: https://github.com/albertoruiz/hmatrix -Synopsis: Linear algebra and numerical computation +Synopsis: Numerical computation Description: Purely functional interface to basic linear algebra and other numerical computations, internally implemented using GSL, BLAS and LAPACK. - + Category: Math tested-with: GHC ==7.8 @@ -45,8 +45,7 @@ extra-source-files: examples/deriv.hs examples/bool.hs examples/multiply.hs -extra-source-files: src/Numeric/LinearAlgebra/LAPACK/lapack-aux.h, - src/Numeric/GSL/gsl-ode.c +extra-source-files: src/Numeric/GSL/gsl-ode.c flag dd description: svd = zgesdd @@ -74,23 +73,15 @@ flag debugnan library - Build-Depends: base, hmatrix-base, - array, - storable-complex, - process, random, - vector >= 0.8, - binary, - deepseq + Build-Depends: base, hmatrix-base, array, vector, + process, random + Extensions: ForeignFunctionInterface, CPP hs-source-dirs: src - Exposed-modules: Data.Packed, - Data.Packed.Vector, - Data.Packed.Matrix, - Data.Packed.Foreign, - Numeric.GSL.Differentiation, + Exposed-modules: Numeric.GSL.Differentiation, Numeric.GSL.Integration, Numeric.GSL.Fourier, Numeric.GSL.Polynomials, @@ -101,24 +92,15 @@ library Numeric.GSL, Numeric.Container, Numeric.LinearAlgebra, - Numeric.LinearAlgebra.LAPACK, Numeric.LinearAlgebra.Algorithms, Numeric.LinearAlgebra.Util, - Data.Packed.ST, - Data.Packed.Development Graphics.Plot, Numeric.HMatrix, Numeric.HMatrix.Data, Numeric.HMatrix.Devel - other-modules: Data.Packed.Internal, - Data.Packed.Internal.Common, - Data.Packed.Internal.Signatures, - Data.Packed.Internal.Vector, - Data.Packed.Internal.Matrix, - Numeric.Random, + other-modules: Numeric.Random, Numeric.GSL.Internal, Numeric.GSL.Vector, - Numeric.Conversion, Numeric.ContainerBoot, Numeric.IO, Numeric.Chain, @@ -127,8 +109,7 @@ library Numeric.LinearAlgebra.Util.Convolution - C-sources: src/Numeric/LinearAlgebra/LAPACK/lapack-aux.c, - src/Numeric/GSL/gsl-aux.c + C-sources: src/Numeric/GSL/gsl-aux.c cc-options: -O4 -msse2 -Wall diff --git a/packages/hmatrix/src/Numeric/Container.hs b/packages/hmatrix/src/Numeric/Container.hs index ff649b7..3a8dd94 100644 --- a/packages/hmatrix/src/Numeric/Container.hs +++ b/packages/hmatrix/src/Numeric/Container.hs @@ -30,7 +30,7 @@ module Numeric.Container ( -- * Basic functions module Data.Packed, konst, build, - constant, linspace, + linspace, diag, ident, ctrans, -- * Generic operations @@ -64,8 +64,7 @@ module Numeric.Container ( fscanfVector, fprintfVector, freadVector, fwriteVector, ) where -import Data.Packed -import Data.Packed.Internal(constantD) +import Data.Packed hiding (stepD, stepF, condD, condF, conjugateC, conjugateQ) import Numeric.ContainerBoot import Numeric.Chain import Numeric.IO @@ -75,15 +74,6 @@ import Numeric.Random ------------------------------------------------------------------ -{- | creates a vector with a given number of equal components: - -@> constant 2 7 -7 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@ --} -constant :: Element a => a -> Int -> Vector a --- constant x n = runSTVector (newVector x n) -constant = constantD-- about 2x faster - {- | Creates a real vector containing a range of values: >>> linspace 5 (-3,7::Double) diff --git a/packages/hmatrix/src/Numeric/ContainerBoot.hs b/packages/hmatrix/src/Numeric/ContainerBoot.hs index ea4262c..ef21763 100644 --- a/packages/hmatrix/src/Numeric/ContainerBoot.hs +++ b/packages/hmatrix/src/Numeric/ContainerBoot.hs @@ -42,7 +42,7 @@ module Numeric.ContainerBoot ( import Data.Packed import Data.Packed.ST as ST import Numeric.Conversion -import Data.Packed.Internal +import Data.Packed.Development import Numeric.GSL.Vector import Data.Complex import Control.Applicative((<*>)) @@ -201,7 +201,7 @@ instance Container Vector Float where equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0 arctan2 = vectorZipF ATan2 scalar x = fromList [x] - konst' = constantD + konst' = constant build' = buildV conj = id cmap = mapVector @@ -229,7 +229,7 @@ instance Container Vector Double where equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0 arctan2 = vectorZipR ATan2 scalar x = fromList [x] - konst' = constantD + konst' = constant build' = buildV conj = id cmap = mapVector @@ -257,7 +257,7 @@ instance Container Vector (Complex Double) where equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0 arctan2 = vectorZipC ATan2 scalar x = fromList [x] - konst' = constantD + konst' = constant build' = buildV conj = conjugateC cmap = mapVector @@ -285,7 +285,7 @@ instance Container Vector (Complex Float) where equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0 arctan2 = vectorZipQ ATan2 scalar x = fromList [x] - konst' = constantD + konst' = constant build' = buildV conj = conjugateQ cmap = mapVector @@ -569,7 +569,7 @@ diag v = diagRect 0 v n n where n = dim v -- | creates the identity matrix of given dimension ident :: (Num a, Element a) => Int -> Matrix a -ident n = diag (constantD 1 n) +ident n = diag (constant 1 n) -------------------------------------------------------- diff --git a/packages/hmatrix/src/Numeric/GSL.hs b/packages/hmatrix/src/Numeric/GSL.hs index 5f39a3e..61b8d7b 100644 --- a/packages/hmatrix/src/Numeric/GSL.hs +++ b/packages/hmatrix/src/Numeric/GSL.hs @@ -1,12 +1,11 @@ {- | Module : Numeric.GSL -Copyright : (c) Alberto Ruiz 2006-7 -License : GPL-style +Copyright : (c) Alberto Ruiz 2006-14 +License : GPL -Maintainer : Alberto Ruiz (aruiz at um dot es) +Maintainer : Alberto Ruiz Stability : provisional -Portability : uses -fffi and -fglasgow-exts This module reexports all available GSL functions. @@ -41,3 +40,4 @@ import Data.Complex -- | This action removes the GSL default error handler (which aborts the program), so that -- GSL errors can be handled by Haskell (using Control.Exception) and ghci doesn't abort. foreign import ccall unsafe "GSL/gsl-aux.h no_abort_on_error" setErrorHandlerOff :: IO () + diff --git a/packages/hmatrix/src/Numeric/GSL/Differentiation.hs b/packages/hmatrix/src/Numeric/GSL/Differentiation.hs index 93c5007..0fb58ef 100644 --- a/packages/hmatrix/src/Numeric/GSL/Differentiation.hs +++ b/packages/hmatrix/src/Numeric/GSL/Differentiation.hs @@ -1,13 +1,10 @@ -{-# OPTIONS #-} ------------------------------------------------------------------------------ {- | Module : Numeric.GSL.Differentiation Copyright : (c) Alberto Ruiz 2006 -License : GPL-style +License : GPL -Maintainer : Alberto Ruiz (aruiz at um dot es) +Maintainer : Alberto Ruiz Stability : provisional -Portability : uses ffi Numerical differentiation. @@ -15,7 +12,8 @@ Numerical differentiation. From the GSL manual: \"The functions described in this chapter compute numerical derivatives by finite differencing. An adaptive algorithm is used to find the best choice of finite difference and to estimate the error in the derivative.\" -} ------------------------------------------------------------------------------ + + module Numeric.GSL.Differentiation ( derivCentral, derivForward, @@ -26,7 +24,7 @@ import Foreign.C.Types import Foreign.Marshal.Alloc(malloc, free) import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr) import Foreign.Storable(peek) -import Data.Packed.Internal(check,(//)) +import Numeric.GSL.Internal import System.IO.Unsafe(unsafePerformIO) derivGen :: diff --git a/packages/hmatrix/src/Numeric/GSL/Fitting.hs b/packages/hmatrix/src/Numeric/GSL/Fitting.hs index c4f3a91..0a92373 100644 --- a/packages/hmatrix/src/Numeric/GSL/Fitting.hs +++ b/packages/hmatrix/src/Numeric/GSL/Fitting.hs @@ -3,9 +3,8 @@ Module : Numeric.GSL.Fitting Copyright : (c) Alberto Ruiz 2010 License : GPL -Maintainer : Alberto Ruiz (aruiz at um dot es) +Maintainer : Alberto Ruiz Stability : provisional -Portability : uses ffi Nonlinear Least-Squares Fitting @@ -42,7 +41,7 @@ expModelDer [a,lambda,b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1] (1.0192487112786812,3.782067731353722e-2)] -} ------------------------------------------------------------------------------ + module Numeric.GSL.Fitting ( -- * Levenberg-Marquardt @@ -51,7 +50,6 @@ module Numeric.GSL.Fitting ( fitModelScaled, fitModel ) where -import Data.Packed.Internal import Numeric.LinearAlgebra import Numeric.GSL.Internal @@ -61,6 +59,9 @@ import System.IO.Unsafe(unsafePerformIO) ------------------------------------------------------------------------- +type TVV = TV (TV Res) +type TVM = TV (TM Res) + data FittingMethod = LevenbergMarquardtScaled -- ^ Interface to gsl_multifit_fdfsolver_lmsder. This is a robust and efficient version of the Levenberg-Marquardt algorithm as implemented in the scaled lmder routine in minpack. Minpack was written by Jorge J. More, Burton S. Garbow and Kenneth E. Hillstrom. | LevenbergMarquardt -- ^ This is an unscaled version of the lmder algorithm. The elements of the diagonal scaling matrix D are set to 1. This algorithm may be useful in circumstances where the scaled version of lmder converges too slowly, or the function is already scaled appropriately. deriving (Enum,Eq,Show,Bounded) diff --git a/packages/hmatrix/src/Numeric/GSL/Fourier.hs b/packages/hmatrix/src/Numeric/GSL/Fourier.hs index 86aedd6..734325b 100644 --- a/packages/hmatrix/src/Numeric/GSL/Fourier.hs +++ b/packages/hmatrix/src/Numeric/GSL/Fourier.hs @@ -1,26 +1,23 @@ -{-# LANGUAGE ForeignFunctionInterface #-} ------------------------------------------------------------------------------ {- | Module : Numeric.GSL.Fourier Copyright : (c) Alberto Ruiz 2006 -License : GPL-style - -Maintainer : Alberto Ruiz (aruiz at um dot es) +License : GPL +Maintainer : Alberto Ruiz Stability : provisional -Portability : uses ffi Fourier Transform. -} ------------------------------------------------------------------------------ + module Numeric.GSL.Fourier ( fft, ifft ) where -import Data.Packed.Internal +import Data.Packed +import Numeric.GSL.Internal import Data.Complex import Foreign.C.Types import System.IO.Unsafe (unsafePerformIO) @@ -30,7 +27,7 @@ genfft code v = unsafePerformIO $ do app2 (c_fft code) vec v vec r "fft" return r -foreign import ccall unsafe "gsl-aux.h fft" c_fft :: CInt -> TCVCV +foreign import ccall unsafe "gsl-aux.h fft" c_fft :: CInt -> TCV (TCV Res) {- | Fast 1D Fourier transform of a 'Vector' @(@'Complex' 'Double'@)@ using /gsl_fft_complex_forward/. It uses the same scaling conventions as GNU Octave. diff --git a/packages/hmatrix/src/Numeric/GSL/Integration.hs b/packages/hmatrix/src/Numeric/GSL/Integration.hs index 5f0a415..9c1d43a 100644 --- a/packages/hmatrix/src/Numeric/GSL/Integration.hs +++ b/packages/hmatrix/src/Numeric/GSL/Integration.hs @@ -1,19 +1,15 @@ -{-# OPTIONS #-} ------------------------------------------------------------------------------ {- | Module : Numeric.GSL.Integration Copyright : (c) Alberto Ruiz 2006 -License : GPL-style - -Maintainer : Alberto Ruiz (aruiz at um dot es) +License : GPL +Maintainer : Alberto Ruiz Stability : provisional -Portability : uses ffi Numerical integration routines. -} ------------------------------------------------------------------------------ + module Numeric.GSL.Integration ( integrateQNG, @@ -28,7 +24,7 @@ import Foreign.C.Types import Foreign.Marshal.Alloc(malloc, free) import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr) import Foreign.Storable(peek) -import Data.Packed.Internal(check,(//)) +import Numeric.GSL.Internal import System.IO.Unsafe(unsafePerformIO) eps = 1e-12 diff --git a/packages/hmatrix/src/Numeric/GSL/Internal.hs b/packages/hmatrix/src/Numeric/GSL/Internal.hs index 69a9750..a1c4e0c 100644 --- a/packages/hmatrix/src/Numeric/GSL/Internal.hs +++ b/packages/hmatrix/src/Numeric/GSL/Internal.hs @@ -1,23 +1,43 @@ +-- | -- Module : Numeric.GSL.Internal -- Copyright : (c) Alberto Ruiz 2009 -- License : GPL --- --- Maintainer : Alberto Ruiz (aruiz at um dot es) +-- Maintainer : Alberto Ruiz -- Stability : provisional --- Portability : uses ffi +-- -- -- Auxiliary functions. -- --- #hide -module Numeric.GSL.Internal where -import Data.Packed.Internal +module Numeric.GSL.Internal( + iv, + mkVecfun, + mkVecVecfun, + mkDoubleVecVecfun, + mkDoublefun, + aux_vTov, + mkVecMatfun, + mkDoubleVecMatfun, + aux_vTom, + createV, + createMIO, + module Data.Packed.Development, + check, + Res,TV,TM,TCV,TCM +) where + +import Data.Packed +import Data.Packed.Development hiding (check) +import Data.Complex import Foreign.Marshal.Array(copyArray) import Foreign.Ptr(Ptr, FunPtr) import Foreign.C.Types +import Foreign.C.String(peekCString) import System.IO.Unsafe(unsafePerformIO) +import Data.Vector.Storable(unsafeWith) +import Control.Monad(when) iv :: (Vector Double -> Double) -> (CInt -> Ptr Double -> Double) iv f n p = f (createV (fromIntegral n) copy "iv") where @@ -74,3 +94,33 @@ createMIO r c fun msg = do res <- createMatrix RowMajor r c app1 fun mat res msg return res + +-------------------------------------------------------------------------------- + +-- | check the error code +check :: String -> IO CInt -> IO () +check msg f = do + err <- f + when (err/=0) $ do + ps <- gsl_strerror err + s <- peekCString ps + error (msg++": "++s) + return () + +-- | description of GSL error codes +foreign import ccall unsafe "gsl_strerror" gsl_strerror :: CInt -> IO (Ptr CChar) + +type PF = Ptr Float +type PD = Ptr Double +type PQ = Ptr (Complex Float) +type PC = Ptr (Complex Double) + +type Res = IO CInt +type TV x = CInt -> PD -> x +type TM x = CInt -> CInt -> PD -> x +type TCV x = CInt -> PC -> x +type TCM x = CInt -> CInt -> PC -> x + +type TVV = TV (TV Res) +type TVM = TV (TM Res) + diff --git a/packages/hmatrix/src/Numeric/GSL/Minimization.hs b/packages/hmatrix/src/Numeric/GSL/Minimization.hs index 1879dab..056d463 100644 --- a/packages/hmatrix/src/Numeric/GSL/Minimization.hs +++ b/packages/hmatrix/src/Numeric/GSL/Minimization.hs @@ -1,13 +1,9 @@ -{-# LANGUAGE ForeignFunctionInterface #-} ------------------------------------------------------------------------------ {- | Module : Numeric.GSL.Minimization Copyright : (c) Alberto Ruiz 2006-9 -License : GPL-style - -Maintainer : Alberto Ruiz (aruiz at um dot es) +License : GPL +Maintainer : Alberto Ruiz Stability : provisional -Portability : uses ffi Minimization of a multidimensional function using some of the algorithms described in: @@ -48,7 +44,7 @@ The nmsimplex2 version is a new O(N) implementation of the earlier O(N^2) nmsimp -} ------------------------------------------------------------------------------ + module Numeric.GSL.Minimization ( minimize, minimizeV, MinimizeMethod(..), minimizeD, minimizeVD, MinimizeMethodD(..), @@ -60,8 +56,7 @@ module Numeric.GSL.Minimization ( ) where -import Data.Packed.Internal -import Data.Packed.Matrix +import Data.Packed import Numeric.GSL.Internal import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr) @@ -112,7 +107,7 @@ uniMinimizeGen m f xmin xl xu epsrel maxit = unsafePerformIO $ do foreign import ccall safe "uniMinimize" - c_uniMinize:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> Double -> TM + c_uniMinize:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> Double -> TM Res data MinimizeMethod = NMSimplex | NMSimplex2 @@ -150,13 +145,13 @@ minimizeV method eps maxit szv f xiv = unsafePerformIO $ do "minimize" let it = round (rawpath @@> (maxit-1,0)) path = takeRows it rawpath - sol = cdat $ dropColumns 3 $ dropRows (it-1) path + sol = flatten $ dropColumns 3 $ dropRows (it-1) path freeHaskellFunPtr fp return (sol, path) foreign import ccall safe "gsl-aux.h minimize" - c_minimize:: CInt -> FunPtr (CInt -> Ptr Double -> Double) -> Double -> CInt -> TVVM + c_minimize:: CInt -> FunPtr (CInt -> Ptr Double -> Double) -> Double -> CInt -> TV(TV(TM Res)) ---------------------------------------------------------------------------------- @@ -207,7 +202,7 @@ minimizeVD method eps maxit istep tol f df xiv = unsafePerformIO $ do "minimizeD" let it = round (rawpath @@> (maxit-1,0)) path = takeRows it rawpath - sol = cdat $ dropColumns 2 $ dropRows (it-1) path + sol = flatten $ dropColumns 2 $ dropRows (it-1) path freeHaskellFunPtr fp freeHaskellFunPtr dfp return (sol,path) @@ -215,9 +210,9 @@ minimizeVD method eps maxit istep tol f df xiv = unsafePerformIO $ do foreign import ccall safe "gsl-aux.h minimizeD" c_minimizeD :: CInt -> FunPtr (CInt -> Ptr Double -> Double) - -> FunPtr TVV + -> FunPtr (TV (TV Res)) -> Double -> Double -> Double -> CInt - -> TVM + -> TV (TM Res) --------------------------------------------------------------------- diff --git a/packages/hmatrix/src/Numeric/GSL/ODE.hs b/packages/hmatrix/src/Numeric/GSL/ODE.hs index 9a29085..7549a65 100644 --- a/packages/hmatrix/src/Numeric/GSL/ODE.hs +++ b/packages/hmatrix/src/Numeric/GSL/ODE.hs @@ -2,10 +2,8 @@ Module : Numeric.GSL.ODE Copyright : (c) Alberto Ruiz 2010 License : GPL - -Maintainer : Alberto Ruiz (aruiz at um dot es) +Maintainer : Alberto Ruiz Stability : provisional -Portability : uses ffi Solution of ordinary differential equation (ODE) initial value problems. @@ -34,7 +32,7 @@ module Numeric.GSL.ODE ( odeSolve, odeSolveV, ODEMethod(..), Jacobian ) where -import Data.Packed.Internal +import Data.Packed import Numeric.GSL.Internal import Foreign.Ptr(FunPtr, nullFunPtr, freeHaskellFunPtr) @@ -43,6 +41,10 @@ import System.IO.Unsafe(unsafePerformIO) ------------------------------------------------------------------------- +type TVV = TV (TV Res) +type TVM = TV (TM Res) +type TVVM = TV (TV (TM Res)) + type Jacobian = Double -> Vector Double -> Matrix Double -- | Stepping functions diff --git a/packages/hmatrix/src/Numeric/GSL/Polynomials.hs b/packages/hmatrix/src/Numeric/GSL/Polynomials.hs index 290c615..b1be85d 100644 --- a/packages/hmatrix/src/Numeric/GSL/Polynomials.hs +++ b/packages/hmatrix/src/Numeric/GSL/Polynomials.hs @@ -1,25 +1,23 @@ -{-# LANGUAGE CPP, ForeignFunctionInterface #-} ------------------------------------------------------------------------------ {- | Module : Numeric.GSL.Polynomials Copyright : (c) Alberto Ruiz 2006 -License : GPL-style - -Maintainer : Alberto Ruiz (aruiz at um dot es) +License : GPL +Maintainer : Alberto Ruiz Stability : provisional -Portability : uses ffi Polynomials. -} ------------------------------------------------------------------------------ + + module Numeric.GSL.Polynomials ( polySolve ) where -import Data.Packed.Internal +import Data.Packed +import Numeric.GSL.Internal import Data.Complex import System.IO.Unsafe (unsafePerformIO) @@ -55,4 +53,5 @@ polySolve' v | dim v > 1 = unsafePerformIO $ do return r | otherwise = error "polySolve on a polynomial of degree zero" -foreign import ccall unsafe "gsl-aux.h polySolve" c_polySolve:: TVCV +foreign import ccall unsafe "gsl-aux.h polySolve" c_polySolve:: TV (TCV Res) + diff --git a/packages/hmatrix/src/Numeric/GSL/Root.hs b/packages/hmatrix/src/Numeric/GSL/Root.hs index 9d561c4..b9f3b94 100644 --- a/packages/hmatrix/src/Numeric/GSL/Root.hs +++ b/packages/hmatrix/src/Numeric/GSL/Root.hs @@ -2,10 +2,8 @@ Module : Numeric.GSL.Root Copyright : (c) Alberto Ruiz 2009 License : GPL - -Maintainer : Alberto Ruiz (aruiz at um dot es) +Maintainer : Alberto Ruiz Stability : provisional -Portability : uses ffi Multidimensional root finding. @@ -41,14 +39,16 @@ module Numeric.GSL.Root ( rootJ, RootMethodJ(..), ) where -import Data.Packed.Internal -import Data.Packed.Matrix +import Data.Packed 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 @@ -76,7 +76,7 @@ uniRootGen m f xl xu epsrel maxit = unsafePerformIO $ do return (sol !! 1, path) foreign import ccall safe "root" - c_root:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> TM + c_root:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> TM Res ------------------------------------------------------------------------- data UniRootMethodJ = UNewton @@ -108,7 +108,7 @@ uniRootJGen m f df x epsrel maxit = unsafePerformIO $ do foreign import ccall safe "rootj" c_rootj :: CInt -> FunPtr (Double -> Double) -> FunPtr (Double -> Double) - -> Double -> CInt -> Double -> TM + -> Double -> CInt -> Double -> TM Res ------------------------------------------------------------------------- diff --git a/packages/hmatrix/src/Numeric/GSL/Vector.hs b/packages/hmatrix/src/Numeric/GSL/Vector.hs index 6204b8e..29c8bb7 100644 --- a/packages/hmatrix/src/Numeric/GSL/Vector.hs +++ b/packages/hmatrix/src/Numeric/GSL/Vector.hs @@ -2,11 +2,9 @@ -- | -- Module : Numeric.GSL.Vector -- Copyright : (c) Alberto Ruiz 2007 --- License : GPL-style --- --- Maintainer : Alberto Ruiz +-- License : GPL +-- Maintainer : Alberto Ruiz -- Stability : provisional --- Portability : portable (uses FFI) -- -- Low level interface to vector operations. -- @@ -20,18 +18,20 @@ module Numeric.GSL.Vector ( FunCodeV(..), vectorMapR, vectorMapC, vectorMapF, vectorMapQ, FunCodeSV(..), vectorMapValR, vectorMapValC, vectorMapValF, vectorMapValQ, FunCodeVV(..), vectorZipR, vectorZipC, vectorZipF, vectorZipQ, - RandDist(..), randomVector + RandDist(..), randomVector, + saveMatrix, + fwriteVector, freadVector, fprintfVector, fscanfVector ) where -import Data.Packed.Internal.Common -import Data.Packed.Internal.Signatures -import Data.Packed.Internal.Vector +import Data.Packed +import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM) import Data.Complex import Foreign.Marshal.Alloc(free) import Foreign.Marshal.Array(newArray) import Foreign.Ptr(Ptr) import Foreign.C.Types +import Foreign.C.String(newCString) import System.IO.Unsafe(unsafePerformIO) import Control.Monad(when) @@ -186,7 +186,7 @@ foreign import ccall unsafe "gsl-aux.h dotC" c_dotC :: TCVCVCV toScalarAux fun code v = unsafePerformIO $ do r <- createVector 1 app2 (fun (fromei code)) vec v vec r "toScalarAux" - return (r `at` 0) + return (r @> 0) vectorMapAux fun code v = unsafePerformIO $ do r <- createVector (dim v) @@ -326,3 +326,118 @@ randomVector seed dist n = unsafePerformIO $ do return r foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> TV + +-------------------------------------------------------------------------------- + +-- | Saves a matrix as 2D ASCII table. +saveMatrix :: FilePath + -> String -- ^ format (%f, %g, %e) + -> Matrix Double + -> IO () +saveMatrix filename fmt m = do + charname <- newCString filename + charfmt <- newCString fmt + let o = if orderOf m == RowMajor then 1 else 0 + app1 (matrix_fprintf charname charfmt o) mat m "matrix_fprintf" + free charname + free charfmt + +foreign import ccall unsafe "matrix_fprintf" matrix_fprintf :: Ptr CChar -> Ptr CChar -> CInt -> TM + +-------------------------------------------------------------------------------- + +-- | Loads a vector from an ASCII file (the number of elements must be known in advance). +fscanfVector :: FilePath -> Int -> IO (Vector Double) +fscanfVector filename n = do + charname <- newCString filename + res <- createVector n + app1 (gsl_vector_fscanf charname) vec res "gsl_vector_fscanf" + free charname + return res + +foreign import ccall unsafe "vector_fscanf" gsl_vector_fscanf:: Ptr CChar -> TV + +-- | Saves the elements of a vector, with a given format (%f, %e, %g), to an ASCII file. +fprintfVector :: FilePath -> String -> Vector Double -> IO () +fprintfVector filename fmt v = do + charname <- newCString filename + charfmt <- newCString fmt + app1 (gsl_vector_fprintf charname charfmt) vec v "gsl_vector_fprintf" + free charname + free charfmt + +foreign import ccall unsafe "vector_fprintf" gsl_vector_fprintf :: Ptr CChar -> Ptr CChar -> TV + +-- | Loads a vector from a binary file (the number of elements must be known in advance). +freadVector :: FilePath -> Int -> IO (Vector Double) +freadVector filename n = do + charname <- newCString filename + res <- createVector n + app1 (gsl_vector_fread charname) vec res "gsl_vector_fread" + free charname + return res + +foreign import ccall unsafe "vector_fread" gsl_vector_fread:: Ptr CChar -> TV + +-- | Saves the elements of a vector to a binary file. +fwriteVector :: FilePath -> Vector Double -> IO () +fwriteVector filename v = do + charname <- newCString filename + app1 (gsl_vector_fwrite charname) vec v "gsl_vector_fwrite" + free charname + +foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV + +type PF = Ptr Float -- +type PD = Ptr Double -- +type PQ = Ptr (Complex Float) -- +type PC = Ptr (Complex Double) -- +type TF = CInt -> PF -> IO CInt -- +type TFF = CInt -> PF -> TF -- +type TFV = CInt -> PF -> TV -- +type TVF = CInt -> PD -> TF -- +type TFFF = CInt -> PF -> TFF -- +type TV = CInt -> PD -> IO CInt -- +type TVV = CInt -> PD -> TV -- +type TVVV = CInt -> PD -> TVV -- +type TFM = CInt -> CInt -> PF -> IO CInt -- +type TFMFM = CInt -> CInt -> PF -> TFM -- +type TFMFMFM = CInt -> CInt -> PF -> TFMFM -- +type TM = CInt -> CInt -> PD -> IO CInt -- +type TMM = CInt -> CInt -> PD -> TM -- +type TVMM = CInt -> PD -> TMM -- +type TMVMM = CInt -> CInt -> PD -> TVMM -- +type TMMM = CInt -> CInt -> PD -> TMM -- +type TVM = CInt -> PD -> TM -- +type TVVM = CInt -> PD -> TVM -- +type TMV = CInt -> CInt -> PD -> TV -- +type TMMV = CInt -> CInt -> PD -> TMV -- +type TMVM = CInt -> CInt -> PD -> TVM -- +type TMMVM = CInt -> CInt -> PD -> TMVM -- +type TCM = CInt -> CInt -> PC -> IO CInt -- +type TCVCM = CInt -> PC -> TCM -- +type TCMCVCM = CInt -> CInt -> PC -> TCVCM -- +type TMCMCVCM = CInt -> CInt -> PD -> TCMCVCM -- +type TCMCMCVCM = CInt -> CInt -> PC -> TCMCVCM -- +type TCMCM = CInt -> CInt -> PC -> TCM -- +type TVCM = CInt -> PD -> TCM -- +type TCMVCM = CInt -> CInt -> PC -> TVCM -- +type TCMCMVCM = CInt -> CInt -> PC -> TCMVCM -- +type TCMCMCM = CInt -> CInt -> PC -> TCMCM -- +type TCV = CInt -> PC -> IO CInt -- +type TCVCV = CInt -> PC -> TCV -- +type TCVCVCV = CInt -> PC -> TCVCV -- +type TCVV = CInt -> PC -> TV -- +type TQV = CInt -> PQ -> IO CInt -- +type TQVQV = CInt -> PQ -> TQV -- +type TQVQVQV = CInt -> PQ -> TQVQV -- +type TQVF = CInt -> PQ -> TF -- +type TQM = CInt -> CInt -> PQ -> IO CInt -- +type TQMQM = CInt -> CInt -> PQ -> TQM -- +type TQMQMQM = CInt -> CInt -> PQ -> TQMQM -- +type TCMCV = CInt -> CInt -> PC -> TCV -- +type TVCV = CInt -> PD -> TCV -- +type TCVM = CInt -> PC -> TM -- +type TMCVM = CInt -> CInt -> PD -> TCVM -- +type TMMCVM = CInt -> CInt -> PD -> TMCVM -- + diff --git a/packages/hmatrix/src/Numeric/IO.hs b/packages/hmatrix/src/Numeric/IO.hs index 836f352..58fa2b4 100644 --- a/packages/hmatrix/src/Numeric/IO.hs +++ b/packages/hmatrix/src/Numeric/IO.hs @@ -20,11 +20,12 @@ module Numeric.IO ( ) where import Data.Packed -import Data.Packed.Internal +import Data.Packed.Development import System.Process(readProcess) import Text.Printf(printf) import Data.List(intersperse) import Data.Complex +import Numeric.GSL.Vector {- | Creates a string from a matrix given a separator and a function to show each entry. Using this function the user can easily define any desired display function: diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/Algorithms.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/Algorithms.hs index 8c4b610..0a6eaa1 100644 --- a/packages/hmatrix/src/Numeric/LinearAlgebra/Algorithms.hs +++ b/packages/hmatrix/src/Numeric/LinearAlgebra/Algorithms.hs @@ -79,8 +79,8 @@ module Numeric.LinearAlgebra.Algorithms ( ) where -import Data.Packed.Internal hiding ((//)) -import Data.Packed.Matrix +import Data.Packed.Development hiding ((//)) +import Data.Packed import Numeric.LinearAlgebra.LAPACK as LAPACK import Data.List(foldl1') import Data.Array -- cgit v1.2.3