1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
|
-- |
-- Module : Numeric.GSL.Internal
-- Copyright : (c) Alberto Ruiz 2009
-- License : GPL
-- Maintainer : Alberto Ruiz
-- Stability : provisional
--
--
-- Auxiliary functions.
--
module Numeric.GSL.Internal(
iv,
mkVecfun,
mkVecVecfun,
mkDoubleVecVecfun,
mkDoublefun,
aux_vTov,
mkVecMatfun,
mkDoubleVecMatfun,
aux_vTom,
createV,
createMIO,
module Numeric.LinearAlgebra.Devel,
check,(#),vec, ww2,
Res,TV,TM,TCV,TCM
) where
import Numeric.LinearAlgebra.HMatrix
import Numeric.LinearAlgebra.Devel hiding (check)
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 as V (unsafeWith,length)
import Control.Monad(when)
iv :: (Vector Double -> Double) -> (CInt -> Ptr Double -> Double)
iv f n p = f (createV (fromIntegral n) copy "iv") where
copy n' q = do
copyArray q p (fromIntegral n')
return 0
-- | conversion of Haskell functions into function pointers that can be used in the C side
foreign import ccall safe "wrapper"
mkVecfun :: (CInt -> Ptr Double -> Double)
-> IO( FunPtr (CInt -> Ptr Double -> Double))
foreign import ccall safe "wrapper"
mkVecVecfun :: TVV -> IO (FunPtr TVV)
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
x = createV (fromIntegral n) copy "aux_vTov"
copy n' q = do
copyArray q p (fromIntegral n')
return 0
g = do unsafeWith v $ \p' -> copyArray r p' (fromIntegral nr)
return 0
foreign import ccall safe "wrapper"
mkVecMatfun :: TVM -> IO (FunPtr TVM)
foreign import ccall safe "wrapper"
mkDoubleVecMatfun :: (Double -> TVM) -> IO (FunPtr (Double -> TVM))
aux_vTom :: (Vector Double -> Matrix Double) -> TVM
aux_vTom f n p rr cr r = g where
v = flatten $ f x
x = createV (fromIntegral n) copy "aux_vTov"
copy n' q = do
copyArray q p (fromIntegral n')
return 0
g = do unsafeWith v $ \p' -> copyArray r p' (fromIntegral $ rr*cr)
return 0
createV n fun msg = unsafePerformIO $ do
r <- createVector n
fun # r #| msg
return r
createMIO r c fun msg = do
res <- createMatrix RowMajor r c
fun # 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)
ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2
vec x f = unsafeWith x $ \p -> do
let v g = do
g (fi $ V.length x) p
f v
{-# INLINE vec #-}
infixl 1 #
a # b = applyRaw a b
{-# INLINE (#) #-}
|
