summaryrefslogtreecommitdiff
path: root/lib/Numeric
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Numeric')
-rw-r--r--lib/Numeric/GSL/Differentiation.hs2
-rw-r--r--lib/Numeric/GSL/Fourier.hs3
-rw-r--r--lib/Numeric/GSL/Integration.hs4
-rw-r--r--lib/Numeric/GSL/Matrix.hs52
-rw-r--r--lib/Numeric/GSL/Minimization.hs34
-rw-r--r--lib/Numeric/GSL/Polynomials.hs3
-rw-r--r--lib/Numeric/GSL/Special/Internal.hs4
-rw-r--r--lib/Numeric/GSL/Vector.hs12
-rw-r--r--lib/Numeric/LinearAlgebra/LAPACK.hs34
9 files changed, 91 insertions, 57 deletions
diff --git a/lib/Numeric/GSL/Differentiation.hs b/lib/Numeric/GSL/Differentiation.hs
index e7fea92..09236bd 100644
--- a/lib/Numeric/GSL/Differentiation.hs
+++ b/lib/Numeric/GSL/Differentiation.hs
@@ -35,7 +35,7 @@ derivGen c h f x = unsafePerformIO $ do
35 r <- malloc 35 r <- malloc
36 e <- malloc 36 e <- malloc
37 fp <- mkfun (\x _ -> f x) 37 fp <- mkfun (\x _ -> f x)
38 c_deriv c fp x h r e // check "deriv" [] 38 c_deriv c fp x h r e // check "deriv"
39 vr <- peek r 39 vr <- peek r
40 ve <- peek e 40 ve <- peek e
41 let result = (vr,ve) 41 let result = (vr,ve)
diff --git a/lib/Numeric/GSL/Fourier.hs b/lib/Numeric/GSL/Fourier.hs
index e975fbf..4b08625 100644
--- a/lib/Numeric/GSL/Fourier.hs
+++ b/lib/Numeric/GSL/Fourier.hs
@@ -26,7 +26,8 @@ import Foreign
26 26
27genfft code v = unsafePerformIO $ do 27genfft code v = unsafePerformIO $ do
28 r <- createVector (dim v) 28 r <- createVector (dim v)
29 c_fft code // vec v // vec r // check "fft" [v] 29 ww2 withVector v withVector r $ \ v r ->
30 c_fft code // v // r // check "fft"
30 return r 31 return r
31 32
32foreign import ccall "gsl-aux.h fft" c_fft :: Int -> TCVCV 33foreign import ccall "gsl-aux.h fft" c_fft :: Int -> TCVCV
diff --git a/lib/Numeric/GSL/Integration.hs b/lib/Numeric/GSL/Integration.hs
index d756417..747b34c 100644
--- a/lib/Numeric/GSL/Integration.hs
+++ b/lib/Numeric/GSL/Integration.hs
@@ -44,7 +44,7 @@ integrateQAGS prec n f a b = unsafePerformIO $ do
44 r <- malloc 44 r <- malloc
45 e <- malloc 45 e <- malloc
46 fp <- mkfun (\x _ -> f x) 46 fp <- mkfun (\x _ -> f x)
47 c_integrate_qags fp a b prec n r e // check "integrate_qags" [] 47 c_integrate_qags fp a b prec n r e // check "integrate_qags"
48 vr <- peek r 48 vr <- peek r
49 ve <- peek e 49 ve <- peek e
50 let result = (vr,ve) 50 let result = (vr,ve)
@@ -75,7 +75,7 @@ integrateQNG prec f a b = unsafePerformIO $ do
75 r <- malloc 75 r <- malloc
76 e <- malloc 76 e <- malloc
77 fp <- mkfun (\x _ -> f x) 77 fp <- mkfun (\x _ -> f x)
78 c_integrate_qng fp a b prec r e // check "integrate_qng" [] 78 c_integrate_qng fp a b prec r e // check "integrate_qng"
79 vr <- peek r 79 vr <- peek r
80 ve <- peek e 80 ve <- peek e
81 let result = (vr,ve) 81 let result = (vr,ve)
diff --git a/lib/Numeric/GSL/Matrix.hs b/lib/Numeric/GSL/Matrix.hs
index e803c53..09a0be4 100644
--- a/lib/Numeric/GSL/Matrix.hs
+++ b/lib/Numeric/GSL/Matrix.hs
@@ -51,7 +51,8 @@ eigSg' m
51 | otherwise = unsafePerformIO $ do 51 | otherwise = unsafePerformIO $ do
52 l <- createVector r 52 l <- createVector r
53 v <- createMatrix RowMajor r r 53 v <- createMatrix RowMajor r r
54 c_eigS // matc m // vec l // matc v // check "eigSg" [cdat m] 54 ww3 withMatrix m withVector l withMatrix v $ \m l v ->
55 c_eigS // m // l // v // check "eigSg"
55 return (l,v) 56 return (l,v)
56 where r = rows m 57 where r = rows m
57foreign import ccall "gsl-aux.h eigensystemR" c_eigS :: TMVM 58foreign import ccall "gsl-aux.h eigensystemR" c_eigS :: TMVM
@@ -84,7 +85,8 @@ eigHg' m
84 | otherwise = unsafePerformIO $ do 85 | otherwise = unsafePerformIO $ do
85 l <- createVector r 86 l <- createVector r
86 v <- createMatrix RowMajor r r 87 v <- createMatrix RowMajor r r
87 c_eigH // matc m // vec l // matc v // check "eigHg" [cdat m] 88 ww3 withMatrix m withVector l withMatrix v $ \m l v ->
89 c_eigH // m // l // v // check "eigHg"
88 return (l,v) 90 return (l,v)
89 where r = rows m 91 where r = rows m
90foreign import ccall "gsl-aux.h eigensystemC" c_eigH :: TCMVCM 92foreign import ccall "gsl-aux.h eigensystemC" c_eigH :: TCMVCM
@@ -120,7 +122,8 @@ svd' x = unsafePerformIO $ do
120 u <- createMatrix RowMajor r c 122 u <- createMatrix RowMajor r c
121 s <- createVector c 123 s <- createVector c
122 v <- createMatrix RowMajor c c 124 v <- createMatrix RowMajor c c
123 c_svd // matc x // matc u // vec s // matc v // check "svdg" [cdat x] 125 ww4 withMatrix x withMatrix u withVector s withMatrix v $ \x u s v ->
126 c_svd // x // u // s // v // check "svdg"
124 return (u,s,v) 127 return (u,s,v)
125 where r = rows x 128 where r = rows x
126 c = cols x 129 c = cols x
@@ -149,7 +152,8 @@ qr = qr' . cmat
149qr' x = unsafePerformIO $ do 152qr' x = unsafePerformIO $ do
150 q <- createMatrix RowMajor r r 153 q <- createMatrix RowMajor r r
151 rot <- createMatrix RowMajor r c 154 rot <- createMatrix RowMajor r c
152 c_qr // matc x // matc q // matc rot // check "qr" [cdat x] 155 ww3 withMatrix x withMatrix q withMatrix rot $ \x q rot ->
156 c_qr // x // q // rot // check "qr"
153 return (q,rot) 157 return (q,rot)
154 where r = rows x 158 where r = rows x
155 c = cols x 159 c = cols x
@@ -161,7 +165,8 @@ qrPacked = qrPacked' . cmat
161qrPacked' x = unsafePerformIO $ do 165qrPacked' x = unsafePerformIO $ do
162 qr <- createMatrix RowMajor r c 166 qr <- createMatrix RowMajor r c
163 tau <- createVector (min r c) 167 tau <- createVector (min r c)
164 c_qrPacked // matc x // matc qr // vec tau // check "qrUnpacked" [cdat x] 168 ww3 withMatrix x withMatrix qr withVector tau $ \x qr tau ->
169 c_qrPacked // x // qr // tau // check "qrUnpacked"
165 return (qr,tau) 170 return (qr,tau)
166 where r = rows x 171 where r = rows x
167 c = cols x 172 c = cols x
@@ -172,9 +177,10 @@ unpackQR (qr,tau) = unpackQR' (cmat qr, tau)
172 177
173unpackQR' (qr,tau) = unsafePerformIO $ do 178unpackQR' (qr,tau) = unsafePerformIO $ do
174 q <- createMatrix RowMajor r r 179 q <- createMatrix RowMajor r r
175 rot <- createMatrix RowMajor r c 180 res <- createMatrix RowMajor r c
176 c_qrUnpack // matc qr // vec tau // matc q // matc rot // check "qrUnpack" [cdat qr,tau] 181 ww4 withMatrix qr withVector tau withMatrix q withMatrix res $ \qr tau q res ->
177 return (q,rot) 182 c_qrUnpack // qr // tau // q // res // check "qrUnpack"
183 return (q,res)
178 where r = rows qr 184 where r = rows qr
179 c = cols qr 185 c = cols qr
180foreign import ccall "gsl-aux.h QRunpack" c_qrUnpack :: TMVMM 186foreign import ccall "gsl-aux.h QRunpack" c_qrUnpack :: TMVMM
@@ -196,20 +202,22 @@ cholR :: Matrix Double -> Matrix Double
196cholR = cholR' . cmat 202cholR = cholR' . cmat
197 203
198cholR' x = unsafePerformIO $ do 204cholR' x = unsafePerformIO $ do
199 res <- createMatrix RowMajor r r 205 r <- createMatrix RowMajor n n
200 c_cholR // matc x // matc res // check "cholR" [cdat x] 206 ww2 withMatrix x withMatrix r $ \x r ->
201 return res 207 c_cholR // x // r // check "cholR"
202 where r = rows x 208 return r
209 where n = rows x
203foreign import ccall "gsl-aux.h cholR" c_cholR :: TMM 210foreign import ccall "gsl-aux.h cholR" c_cholR :: TMM
204 211
205cholC :: Matrix (Complex Double) -> Matrix (Complex Double) 212cholC :: Matrix (Complex Double) -> Matrix (Complex Double)
206cholC = cholC' . cmat 213cholC = cholC' . cmat
207 214
208cholC' x = unsafePerformIO $ do 215cholC' x = unsafePerformIO $ do
209 res <- createMatrix RowMajor r r 216 r <- createMatrix RowMajor n n
210 c_cholC // matc x // matc res // check "cholC" [cdat x] 217 ww2 withMatrix x withMatrix r $ \x r ->
211 return res 218 c_cholC // x // r // check "cholC"
212 where r = rows x 219 return r
220 where n = rows x
213foreign import ccall "gsl-aux.h cholC" c_cholC :: TCMCM 221foreign import ccall "gsl-aux.h cholC" c_cholC :: TCMCM
214 222
215 223
@@ -223,7 +231,8 @@ luSolveR a b = luSolveR' (cmat a) (cmat b)
223luSolveR' a b 231luSolveR' a b
224 | n1==n2 && n1==r = unsafePerformIO $ do 232 | n1==n2 && n1==r = unsafePerformIO $ do
225 s <- createMatrix RowMajor r c 233 s <- createMatrix RowMajor r c
226 c_luSolveR // matc a // matc b // matc s // check "luSolveR" [cdat a, cdat b] 234 ww3 withMatrix a withMatrix b withMatrix s $ \ a b s ->
235 c_luSolveR // a // b // s // check "luSolveR"
227 return s 236 return s
228 | otherwise = error "luSolveR of nonsquare matrix" 237 | otherwise = error "luSolveR of nonsquare matrix"
229 where n1 = rows a 238 where n1 = rows a
@@ -240,7 +249,8 @@ luSolveC a b = luSolveC' (cmat a) (cmat b)
240luSolveC' a b 249luSolveC' a b
241 | n1==n2 && n1==r = unsafePerformIO $ do 250 | n1==n2 && n1==r = unsafePerformIO $ do
242 s <- createMatrix RowMajor r c 251 s <- createMatrix RowMajor r c
243 c_luSolveC // matc a // matc b // matc s // check "luSolveC" [cdat a, cdat b] 252 ww3 withMatrix a withMatrix b withMatrix s $ \ a b s ->
253 c_luSolveC // a // b // s // check "luSolveC"
244 return s 254 return s
245 | otherwise = error "luSolveC of nonsquare matrix" 255 | otherwise = error "luSolveC of nonsquare matrix"
246 where n1 = rows a 256 where n1 = rows a
@@ -256,7 +266,8 @@ luRaux = luRaux' . cmat
256 266
257luRaux' x = unsafePerformIO $ do 267luRaux' x = unsafePerformIO $ do
258 res <- createVector (r*r+r+1) 268 res <- createVector (r*r+r+1)
259 c_luRaux // matc x // vec res // check "luRaux" [cdat x] 269 ww2 withMatrix x withVector res $ \x res ->
270 c_luRaux // x // res // check "luRaux"
260 return res 271 return res
261 where r = rows x 272 where r = rows x
262 c = cols x 273 c = cols x
@@ -269,7 +280,8 @@ luCaux = luCaux' . cmat
269 280
270luCaux' x = unsafePerformIO $ do 281luCaux' x = unsafePerformIO $ do
271 res <- createVector (r*r+r+1) 282 res <- createVector (r*r+r+1)
272 c_luCaux // matc x // vec res // check "luCaux" [cdat x] 283 ww2 withMatrix x withVector res $ \x res ->
284 c_luCaux // x // res // check "luCaux"
273 return res 285 return res
274 where r = rows x 286 where r = rows x
275 c = cols x 287 c = cols x
diff --git a/lib/Numeric/GSL/Minimization.hs b/lib/Numeric/GSL/Minimization.hs
index f523849..e44b2e5 100644
--- a/lib/Numeric/GSL/Minimization.hs
+++ b/lib/Numeric/GSL/Minimization.hs
@@ -26,7 +26,6 @@ import Data.Packed.Matrix
26import Foreign 26import Foreign
27import Complex 27import Complex
28 28
29
30------------------------------------------------------------------------- 29-------------------------------------------------------------------------
31 30
32{- | The method of Nelder and Mead, implemented by /gsl_multimin_fminimizer_nmsimplex/. The gradient of the function is not required. This is the example in the GSL manual: 31{- | The method of Nelder and Mead, implemented by /gsl_multimin_fminimizer_nmsimplex/. The gradient of the function is not required. This is the example in the GSL manual:
@@ -85,9 +84,10 @@ minimizeNMSimplex f xi sz tol maxit = unsafePerformIO $ do
85 szv = fromList sz 84 szv = fromList sz
86 n = dim xiv 85 n = dim xiv
87 fp <- mkVecfun (iv (f.toList)) 86 fp <- mkVecfun (iv (f.toList))
88 rawpath <- createMIO maxit (n+3) 87 rawpath <- ww2 withVector xiv withVector szv $ \xiv szv ->
89 (c_minimizeNMSimplex fp tol maxit // vec xiv // vec szv) 88 createMIO maxit (n+3)
90 "minimizeNMSimplex" [xiv,szv] 89 (c_minimizeNMSimplex fp tol maxit // xiv // szv)
90 "minimizeNMSimplex"
91 let it = round (rawpath @@> (maxit-1,0)) 91 let it = round (rawpath @@> (maxit-1,0))
92 path = takeRows it rawpath 92 path = takeRows it rawpath
93 [sol] = toLists $ dropRows (it-1) path 93 [sol] = toLists $ dropRows (it-1) path
@@ -150,9 +150,10 @@ minimizeConjugateGradient istep minimpar tol maxit f df xi = unsafePerformIO $ d
150 df' = (fromList . df . toList) 150 df' = (fromList . df . toList)
151 fp <- mkVecfun (iv f') 151 fp <- mkVecfun (iv f')
152 dfp <- mkVecVecfun (aux_vTov df') 152 dfp <- mkVecVecfun (aux_vTov df')
153 rawpath <- createMIO maxit (n+2) 153 rawpath <- withVector xiv $ \xiv ->
154 (c_minimizeConjugateGradient fp dfp istep minimpar tol maxit // vec xiv) 154 createMIO maxit (n+2)
155 "minimizeDerivV" [xiv] 155 (c_minimizeConjugateGradient fp dfp istep minimpar tol maxit // xiv)
156 "minimizeDerivV"
156 let it = round (rawpath @@> (maxit-1,0)) 157 let it = round (rawpath @@> (maxit-1,0))
157 path = takeRows it rawpath 158 path = takeRows it rawpath
158 sol = toList $ cdat $ dropColumns 2 $ dropRows (it-1) path 159 sol = toList $ cdat $ dropColumns 2 $ dropRows (it-1) path
@@ -169,7 +170,7 @@ foreign import ccall "gsl-aux.h minimizeWithDeriv"
169 170
170--------------------------------------------------------------------- 171---------------------------------------------------------------------
171iv :: (Vector Double -> Double) -> (Int -> Ptr Double -> Double) 172iv :: (Vector Double -> Double) -> (Int -> Ptr Double -> Double)
172iv f n p = f (createV n copy "iv" []) where 173iv f n p = f (createV n copy "iv") where
173 copy n q = do 174 copy n q = do
174 copyArray q p n 175 copyArray q p n
175 return 0 176 return 0
@@ -187,25 +188,30 @@ foreign import ccall "wrapper"
187aux_vTov :: (Vector Double -> Vector Double) -> (Int -> Ptr Double -> Ptr Double -> IO()) 188aux_vTov :: (Vector Double -> Vector Double) -> (Int -> Ptr Double -> Ptr Double -> IO())
188aux_vTov f n p r = g where 189aux_vTov f n p r = g where
189 v@V {fptr = pr} = f x 190 v@V {fptr = pr} = f x
190 x = createV n copy "aux_vTov" [] 191 x = createV n copy "aux_vTov"
191 copy n q = do 192 copy n q = do
192 copyArray q p n 193 copyArray q p n
193 return 0 194 return 0
194 g = withForeignPtr pr $ \_ -> copyArray r (ptr v) n 195 g = withForeignPtr pr $ \p -> copyArray r p n
195 196
196-------------------------------------------------------------------- 197--------------------------------------------------------------------
197 198
198createV n fun msg ptrs = unsafePerformIO $ do 199
200createV n fun msg = unsafePerformIO $ do
199 r <- createVector n 201 r <- createVector n
200 fun // vec r // check msg ptrs 202 withVector r $ \ r ->
203 fun // r // check msg
201 return r 204 return r
202 205
206{-
203createM r c fun msg ptrs = unsafePerformIO $ do 207createM r c fun msg ptrs = unsafePerformIO $ do
204 r <- createMatrix RowMajor r c 208 r <- createMatrix RowMajor r c
205 fun // matc r // check msg ptrs 209 fun // matc r // check msg ptrs
206 return r 210 return r
211-}
207 212
208createMIO r c fun msg ptrs = do 213createMIO r c fun msg = do
209 r <- createMatrix RowMajor r c 214 r <- createMatrix RowMajor r c
210 fun // matc r // check msg ptrs 215 withMatrix r $ \ r ->
216 fun // r // check msg
211 return r 217 return r
diff --git a/lib/Numeric/GSL/Polynomials.hs b/lib/Numeric/GSL/Polynomials.hs
index 42694f0..e663711 100644
--- a/lib/Numeric/GSL/Polynomials.hs
+++ b/lib/Numeric/GSL/Polynomials.hs
@@ -47,7 +47,8 @@ polySolve = toList . polySolve' . fromList
47polySolve' :: Vector Double -> Vector (Complex Double) 47polySolve' :: Vector Double -> Vector (Complex Double)
48polySolve' v | dim v > 1 = unsafePerformIO $ do 48polySolve' v | dim v > 1 = unsafePerformIO $ do
49 r <- createVector (dim v-1) 49 r <- createVector (dim v-1)
50 c_polySolve // vec v // vec r // check "polySolve" [v] 50 ww2 withVector v withVector r $ \ v r ->
51 c_polySolve // v // r // check "polySolve"
51 return r 52 return r
52 | otherwise = error "polySolve on a polynomial of degree zero" 53 | otherwise = error "polySolve on a polynomial of degree zero"
53 54
diff --git a/lib/Numeric/GSL/Special/Internal.hs b/lib/Numeric/GSL/Special/Internal.hs
index a08809b..ca36009 100644
--- a/lib/Numeric/GSL/Special/Internal.hs
+++ b/lib/Numeric/GSL/Special/Internal.hs
@@ -45,7 +45,7 @@ type Size_t = Int
45createSFR :: Storable a => String -> (Ptr a -> IO Int) -> (a, a) 45createSFR :: Storable a => String -> (Ptr a -> IO Int) -> (a, a)
46createSFR s f = unsafePerformIO $ do 46createSFR s f = unsafePerformIO $ do
47 p <- mallocArray 2 47 p <- mallocArray 2
48 f p // check s [] 48 f p // check s
49 [val,err] <- peekArray 2 p 49 [val,err] <- peekArray 2 p
50 free p 50 free p
51 return (val,err) 51 return (val,err)
@@ -60,7 +60,7 @@ createSFR_E10 s f = unsafePerformIO $ do
60 let sd = sizeOf (0::Double) 60 let sd = sizeOf (0::Double)
61 let si = sizeOf (0::Int) 61 let si = sizeOf (0::Int)
62 p <- mallocBytes (2*sd + si) 62 p <- mallocBytes (2*sd + si)
63 f p // check s [] 63 f p // check s
64 val <- peekByteOff p 0 64 val <- peekByteOff p 0
65 err <- peekByteOff p sd 65 err <- peekByteOff p sd
66 expo <- peekByteOff p (2*sd) 66 expo <- peekByteOff p (2*sd)
diff --git a/lib/Numeric/GSL/Vector.hs b/lib/Numeric/GSL/Vector.hs
index d94b377..65f3a2e 100644
--- a/lib/Numeric/GSL/Vector.hs
+++ b/lib/Numeric/GSL/Vector.hs
@@ -73,24 +73,28 @@ data FunCodeS = Norm2
73 73
74toScalarAux fun code v = unsafePerformIO $ do 74toScalarAux fun code v = unsafePerformIO $ do
75 r <- createVector 1 75 r <- createVector 1
76 fun (fromEnum code) // vec v // vec r // check "toScalarAux" [v] 76 ww2 withVector v withVector r $ \v r ->
77 fun (fromEnum code) // v // r // check "toScalarAux"
77 return (r `at` 0) 78 return (r `at` 0)
78 79
79vectorMapAux fun code v = unsafePerformIO $ do 80vectorMapAux fun code v = unsafePerformIO $ do
80 r <- createVector (dim v) 81 r <- createVector (dim v)
81 fun (fromEnum code) // vec v // vec r // check "vectorMapAux" [v] 82 ww2 withVector v withVector r $ \v r ->
83 fun (fromEnum code) // v // r // check "vectorMapAux"
82 return r 84 return r
83 85
84vectorMapValAux fun code val v = unsafePerformIO $ do 86vectorMapValAux fun code val v = unsafePerformIO $ do
85 r <- createVector (dim v) 87 r <- createVector (dim v)
86 pval <- newArray [val] 88 pval <- newArray [val]
87 fun (fromEnum code) pval // vec v // vec r // check "vectorMapValAux" [v] 89 ww2 withVector v withVector r $ \v r ->
90 fun (fromEnum code) pval // v // r // check "vectorMapValAux"
88 free pval 91 free pval
89 return r 92 return r
90 93
91vectorZipAux fun code u v = unsafePerformIO $ do 94vectorZipAux fun code u v = unsafePerformIO $ do
92 r <- createVector (dim u) 95 r <- createVector (dim u)
93 fun (fromEnum code) // vec u // vec v // vec r // check "vectorZipAux" [u,v] 96 ww3 withVector u withVector v withVector r $ \u v r ->
97 fun (fromEnum code) // u // v // r // check "vectorZipAux"
94 return r 98 return r
95 99
96--------------------------------------------------------------------- 100---------------------------------------------------------------------
diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs
index 315be17..19516e3 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK.hs
+++ b/lib/Numeric/LinearAlgebra/LAPACK.hs
@@ -61,7 +61,8 @@ svdAux f st x = unsafePerformIO $ do
61 u <- createMatrix ColumnMajor r r 61 u <- createMatrix ColumnMajor r r
62 s <- createVector (min r c) 62 s <- createVector (min r c)
63 v <- createMatrix ColumnMajor c c 63 v <- createMatrix ColumnMajor c c
64 f // matf x // matf u // vec s // matf v // check st [fdat x] 64 ww4 withMatrix x withMatrix u withVector s withMatrix v $ \x u s v ->
65 f // x // u // s // v // check st
65 return (u,s,trans v) 66 return (u,s,trans v)
66 where r = rows x 67 where r = rows x
67 c = cols x 68 c = cols x
@@ -73,7 +74,8 @@ eigAux f st m
73 l <- createVector r 74 l <- createVector r
74 v <- createMatrix ColumnMajor r r 75 v <- createMatrix ColumnMajor r r
75 dummy <- createMatrix ColumnMajor 1 1 76 dummy <- createMatrix ColumnMajor 1 1
76 f // matf m // matf dummy // vec l // matf v // check st [fdat m] 77 ww4 withMatrix m withMatrix dummy withVector l withMatrix v $ \m dummy l v ->
78 f // m // dummy // l // v // check st
77 return (l,v) 79 return (l,v)
78 where r = rows m 80 where r = rows m
79 81
@@ -115,7 +117,8 @@ eigRaux m
115 l <- createVector r 117 l <- createVector r
116 v <- createMatrix ColumnMajor r r 118 v <- createMatrix ColumnMajor r r
117 dummy <- createMatrix ColumnMajor 1 1 119 dummy <- createMatrix ColumnMajor 1 1
118 dgeev // matf m // matf dummy // vec l // matf v // check "eigR" [fdat m] 120 ww4 withMatrix m withMatrix dummy withVector l withMatrix v $ \m dummy l v ->
121 dgeev // m // dummy // l // v // check "eigR"
119 return (l,v) 122 return (l,v)
120 where r = rows m 123 where r = rows m
121 124
@@ -144,7 +147,8 @@ eigS' m
144 | otherwise = unsafePerformIO $ do 147 | otherwise = unsafePerformIO $ do
145 l <- createVector r 148 l <- createVector r
146 v <- createMatrix ColumnMajor r r 149 v <- createMatrix ColumnMajor r r
147 dsyev // matf m // vec l // matf v // check "eigS" [fdat m] 150 ww3 withMatrix m withVector l withMatrix v $ \m l v ->
151 dsyev // m // l // v // check "eigS"
148 return (l,v) 152 return (l,v)
149 where r = rows m 153 where r = rows m
150 154
@@ -166,7 +170,8 @@ eigH' m
166 | otherwise = unsafePerformIO $ do 170 | otherwise = unsafePerformIO $ do
167 l <- createVector r 171 l <- createVector r
168 v <- createMatrix ColumnMajor r r 172 v <- createMatrix ColumnMajor r r
169 zheev // matf m // vec l // matf v // check "eigH" [fdat m] 173 ww3 withMatrix m withVector l withMatrix v $ \m l v ->
174 zheev // m // l // v // check "eigH"
170 return (l,v) 175 return (l,v)
171 where r = rows m 176 where r = rows m
172 177
@@ -177,7 +182,8 @@ foreign import ccall "LAPACK/lapack-aux.h linearSolveC_l" zgesv :: TCMCMCM
177linearSolveSQAux f st a b 182linearSolveSQAux f st a b
178 | n1==n2 && n1==r = unsafePerformIO $ do 183 | n1==n2 && n1==r = unsafePerformIO $ do
179 s <- createMatrix ColumnMajor r c 184 s <- createMatrix ColumnMajor r c
180 f // matf a // matf b // matf s // check st [fdat a, fdat b] 185 ww3 withMatrix a withMatrix b withMatrix s $ \a b s ->
186 f // a // b // s // check st
181 return s 187 return s
182 | otherwise = error $ st ++ " of nonsquare matrix" 188 | otherwise = error $ st ++ " of nonsquare matrix"
183 where n1 = rows a 189 where n1 = rows a
@@ -201,7 +207,8 @@ foreign import ccall "LAPACK/lapack-aux.h linearSolveSVDC_l" zgelss :: Double ->
201 207
202linearSolveAux f st a b = unsafePerformIO $ do 208linearSolveAux f st a b = unsafePerformIO $ do
203 r <- createMatrix ColumnMajor (max m n) nrhs 209 r <- createMatrix ColumnMajor (max m n) nrhs
204 f // matf a // matf b // matf r // check st [fdat a, fdat b] 210 ww3 withMatrix a withMatrix b withMatrix r $ \a b r ->
211 f // a // b // r // check st
205 return r 212 return r
206 where m = rows a 213 where m = rows a
207 n = cols a 214 n = cols a
@@ -251,7 +258,8 @@ cholS = cholAux dpotrf "cholS" . fmat
251 258
252cholAux f st a = unsafePerformIO $ do 259cholAux f st a = unsafePerformIO $ do
253 r <- createMatrix ColumnMajor n n 260 r <- createMatrix ColumnMajor n n
254 f // matf a // matf r // check st [fdat a] 261 ww2 withMatrix a withMatrix r $ \a r ->
262 f // a // r // check st
255 return r 263 return r
256 where n = rows a 264 where n = rows a
257 265
@@ -270,8 +278,8 @@ qrC = qrAux zgeqr2 "qrC" . fmat
270qrAux f st a = unsafePerformIO $ do 278qrAux f st a = unsafePerformIO $ do
271 r <- createMatrix ColumnMajor m n 279 r <- createMatrix ColumnMajor m n
272 tau <- createVector mn 280 tau <- createVector mn
273 withForeignPtr (fptr $ fdat $ a) $ \p -> 281 ww3 withMatrix a withMatrix r withVector tau $ \ a r tau ->
274 f m n p // vec tau // matf r // check st [fdat a] 282 f // a // tau // r // check st
275 return (r,tau) 283 return (r,tau)
276 where m = rows a 284 where m = rows a
277 n = cols a 285 n = cols a
@@ -292,7 +300,8 @@ hessC = hessAux zgehrd "hessC" . fmat
292hessAux f st a = unsafePerformIO $ do 300hessAux f st a = unsafePerformIO $ do
293 r <- createMatrix ColumnMajor m n 301 r <- createMatrix ColumnMajor m n
294 tau <- createVector (mn-1) 302 tau <- createVector (mn-1)
295 f // matf a // vec tau // matf r // check st [fdat a] 303 ww3 withMatrix a withMatrix r withVector tau $ \ a r tau ->
304 f // a // tau // r // check st
296 return (r,tau) 305 return (r,tau)
297 where m = rows a 306 where m = rows a
298 n = cols a 307 n = cols a
@@ -313,7 +322,8 @@ schurC = schurAux zgees "schurC" . fmat
313schurAux f st a = unsafePerformIO $ do 322schurAux f st a = unsafePerformIO $ do
314 u <- createMatrix ColumnMajor n n 323 u <- createMatrix ColumnMajor n n
315 s <- createMatrix ColumnMajor n n 324 s <- createMatrix ColumnMajor n n
316 f // matf a // matf u // matf s // check st [fdat a] 325 ww3 withMatrix a withMatrix u withMatrix s $ \ a u s ->
326 f // a // u // s // check st
317 return (u,s) 327 return (u,s)
318 where n = rows a 328 where n = rows a
319 329