diff options
Diffstat (limited to 'lib')
-rw-r--r-- | lib/Data/Packed/Internal/Common.hs | 5 | ||||
-rw-r--r-- | lib/Data/Packed/Internal/Matrix.hs | 41 | ||||
-rw-r--r-- | lib/Data/Packed/Internal/Vector.hs | 34 | ||||
-rw-r--r-- | lib/Numeric/GSL/Differentiation.hs | 2 | ||||
-rw-r--r-- | lib/Numeric/GSL/Fourier.hs | 3 | ||||
-rw-r--r-- | lib/Numeric/GSL/Integration.hs | 4 | ||||
-rw-r--r-- | lib/Numeric/GSL/Matrix.hs | 52 | ||||
-rw-r--r-- | lib/Numeric/GSL/Minimization.hs | 34 | ||||
-rw-r--r-- | lib/Numeric/GSL/Polynomials.hs | 3 | ||||
-rw-r--r-- | lib/Numeric/GSL/Special/Internal.hs | 4 | ||||
-rw-r--r-- | lib/Numeric/GSL/Vector.hs | 12 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/LAPACK.hs | 34 |
12 files changed, 146 insertions, 82 deletions
diff --git a/lib/Data/Packed/Internal/Common.hs b/lib/Data/Packed/Internal/Common.hs index 2b3ec28..c3a733c 100644 --- a/lib/Data/Packed/Internal/Common.hs +++ b/lib/Data/Packed/Internal/Common.hs | |||
@@ -60,6 +60,11 @@ common f = commonval . map f where | |||
60 | infixl 0 // | 60 | infixl 0 // |
61 | (//) = flip ($) | 61 | (//) = flip ($) |
62 | 62 | ||
63 | -- hmm.. | ||
64 | ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2 | ||
65 | ww3 w1 o1 w2 o2 w3 o3 f = w1 o1 $ \a1 -> ww2 w2 o2 w3 o3 (f a1) | ||
66 | ww4 w1 o1 w2 o2 w3 o3 w4 o4 f = w1 o1 $ \a1 -> ww3 w2 o2 w3 o3 w4 o4 (f a1) | ||
67 | |||
63 | -- GSL error codes are <= 1024 | 68 | -- GSL error codes are <= 1024 |
64 | -- | error codes for the auxiliary functions required by the wrappers | 69 | -- | error codes for the auxiliary functions required by the wrappers |
65 | errorCode :: Int -> String | 70 | errorCode :: Int -> String |
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs index 5617996..90a96b5 100644 --- a/lib/Data/Packed/Internal/Matrix.hs +++ b/lib/Data/Packed/Internal/Matrix.hs | |||
@@ -79,9 +79,20 @@ cmat MF {rows = r, cols = c, fdat = d } = MC {rows = r, cols = c, cdat = transda | |||
79 | fmat m@MF{} = m | 79 | fmat m@MF{} = m |
80 | fmat MC {rows = r, cols = c, cdat = d } = MF {rows = r, cols = c, fdat = transdata c d r} | 80 | fmat MC {rows = r, cols = c, cdat = d } = MF {rows = r, cols = c, fdat = transdata c d r} |
81 | 81 | ||
82 | matc m f = f (rows m) (cols m) (ptr (cdat m)) | 82 | --matc m f = f (rows m) (cols m) (ptr (cdat m)) |
83 | matf m f = f (rows m) (cols m) (ptr (fdat m)) | 83 | --matf m f = f (rows m) (cols m) (ptr (fdat m)) |
84 | 84 | ||
85 | withMatrix MC {rows = r, cols = c, cdat = d } f = | ||
86 | withForeignPtr (fptr d) $ \p -> do | ||
87 | let m f = do | ||
88 | f r c p | ||
89 | f m | ||
90 | |||
91 | withMatrix MF {rows = r, cols = c, fdat = d } f = | ||
92 | withForeignPtr (fptr d) $ \p -> do | ||
93 | let m f = do | ||
94 | f r c p | ||
95 | f m | ||
85 | 96 | ||
86 | {- | Creates a vector by concatenation of rows | 97 | {- | Creates a vector by concatenation of rows |
87 | 98 | ||
@@ -236,7 +247,9 @@ transdataAux fun c1 d c2 = | |||
236 | then d | 247 | then d |
237 | else unsafePerformIO $ do | 248 | else unsafePerformIO $ do |
238 | v <- createVector (dim d) | 249 | v <- createVector (dim d) |
239 | fun r1 c1 (ptr d) r2 c2 (ptr v) // check "transdataAux" [d,v] | 250 | withForeignPtr (fptr d) $ \pd -> |
251 | withForeignPtr (fptr v) $ \pv -> | ||
252 | fun r1 c1 pd r2 c2 pv // check "transdataAux" | ||
240 | -- putStrLn $ "---> transdataAux" ++ show (toList d) ++ show (toList v) | 253 | -- putStrLn $ "---> transdataAux" ++ show (toList d) ++ show (toList v) |
241 | return v | 254 | return v |
242 | where r1 = dim d `div` c1 | 255 | where r1 = dim d `div` c1 |
@@ -250,8 +263,8 @@ foreign import ccall safe "auxi.h transC" | |||
250 | 263 | ||
251 | ------------------------------------------------------------------ | 264 | ------------------------------------------------------------------ |
252 | 265 | ||
253 | gmatC MF {rows = r, cols = c, fdat = d} f = f 1 c r (ptr d) | 266 | gmatC MF {rows = r, cols = c, fdat = d} p f = f 1 c r p |
254 | gmatC MC {rows = r, cols = c, cdat = d} f = f 0 r c (ptr d) | 267 | gmatC MC {rows = r, cols = c, cdat = d} p f = f 0 r c p |
255 | 268 | ||
256 | dtt MC { cdat = d } = d | 269 | dtt MC { cdat = d } = d |
257 | dtt MF { fdat = d } = d | 270 | dtt MF { fdat = d } = d |
@@ -260,7 +273,9 @@ multiplyAux fun a b = unsafePerformIO $ do | |||
260 | when (cols a /= rows b) $ error $ "inconsistent dimensions in contraction "++ | 273 | when (cols a /= rows b) $ error $ "inconsistent dimensions in contraction "++ |
261 | show (rows a,cols a) ++ " x " ++ show (rows b, cols b) | 274 | show (rows a,cols a) ++ " x " ++ show (rows b, cols b) |
262 | r <- createMatrix RowMajor (rows a) (cols b) | 275 | r <- createMatrix RowMajor (rows a) (cols b) |
263 | fun // gmatC a // gmatC b // matc r // check "multiplyAux" [dtt a, dtt b, cdat r] | 276 | withForeignPtr (fptr (dtt a)) $ \pa -> withForeignPtr (fptr (dtt b)) $ \pb -> |
277 | withMatrix r $ \r -> | ||
278 | fun // gmatC a pa // gmatC b pb // r // check "multiplyAux" | ||
264 | return r | 279 | return r |
265 | 280 | ||
266 | multiplyR = multiplyAux cmultiplyR | 281 | multiplyR = multiplyAux cmultiplyR |
@@ -293,7 +308,8 @@ subMatrixR :: (Int,Int) -> (Int,Int) -> Matrix Double -> Matrix Double | |||
293 | subMatrixR (r0,c0) (rt,ct) x' = unsafePerformIO $ do | 308 | subMatrixR (r0,c0) (rt,ct) x' = unsafePerformIO $ do |
294 | r <- createMatrix RowMajor rt ct | 309 | r <- createMatrix RowMajor rt ct |
295 | let x = cmat x' | 310 | let x = cmat x' |
296 | c_submatrixR r0 (r0+rt-1) c0 (c0+ct-1) // matc x // matc r // check "subMatrixR" [cdat x] | 311 | ww2 withMatrix x withMatrix r $ \x r -> |
312 | c_submatrixR r0 (r0+rt-1) c0 (c0+ct-1) // x // r // check "subMatrixR" | ||
297 | return r | 313 | return r |
298 | foreign import ccall "auxi.h submatrixR" c_submatrixR :: Int -> Int -> Int -> Int -> TMM | 314 | foreign import ccall "auxi.h submatrixR" c_submatrixR :: Int -> Int -> Int -> Int -> TMM |
299 | 315 | ||
@@ -317,8 +333,9 @@ subMatrix = subMatrixD | |||
317 | 333 | ||
318 | diagAux fun msg (v@V {dim = n}) = unsafePerformIO $ do | 334 | diagAux fun msg (v@V {dim = n}) = unsafePerformIO $ do |
319 | m <- createMatrix RowMajor n n | 335 | m <- createMatrix RowMajor n n |
320 | fun // vec v // matc m // check msg [v] | 336 | ww2 withVector v withMatrix m $ \v m -> |
321 | return m -- {tdat = dat m} | 337 | fun // v // m // check msg |
338 | return m | ||
322 | 339 | ||
323 | -- | diagonal matrix from a real vector | 340 | -- | diagonal matrix from a real vector |
324 | diagR :: Vector Double -> Matrix Double | 341 | diagR :: Vector Double -> Matrix Double |
@@ -339,7 +356,8 @@ diag = diagD | |||
339 | constantAux fun x n = unsafePerformIO $ do | 356 | constantAux fun x n = unsafePerformIO $ do |
340 | v <- createVector n | 357 | v <- createVector n |
341 | px <- newArray [x] | 358 | px <- newArray [x] |
342 | fun px // vec v // check "constantAux" [] | 359 | withVector v $ \v -> |
360 | fun px // v // check "constantAux" | ||
343 | free px | 361 | free px |
344 | return v | 362 | return v |
345 | 363 | ||
@@ -385,7 +403,8 @@ fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double) | |||
385 | fromFile filename (r,c) = do | 403 | fromFile filename (r,c) = do |
386 | charname <- newCString filename | 404 | charname <- newCString filename |
387 | res <- createMatrix RowMajor r c | 405 | res <- createMatrix RowMajor r c |
388 | c_gslReadMatrix charname // matc res // check "gslReadMatrix" [] | 406 | withMatrix res $ \res -> |
407 | c_gslReadMatrix charname // res // check "gslReadMatrix" | ||
389 | --free charname -- TO DO: free the auxiliary CString | 408 | --free charname -- TO DO: free the auxiliary CString |
390 | return res | 409 | return res |
391 | foreign import ccall "auxi.h matrix_fscanf" c_gslReadMatrix:: Ptr CChar -> TM | 410 | foreign import ccall "auxi.h matrix_fscanf" c_gslReadMatrix:: Ptr CChar -> TM |
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs index dc86484..386ebb5 100644 --- a/lib/Data/Packed/Internal/Vector.hs +++ b/lib/Data/Packed/Internal/Vector.hs | |||
@@ -31,11 +31,11 @@ data Vector t = V { dim :: Int -- ^ number of elements | |||
31 | , fptr :: ForeignPtr t -- ^ foreign pointer to the memory block | 31 | , fptr :: ForeignPtr t -- ^ foreign pointer to the memory block |
32 | } | 32 | } |
33 | 33 | ||
34 | ptr (V _ fptr) = unsafeForeignPtrToPtr fptr | 34 | --ptr (V _ fptr) = unsafeForeignPtrToPtr fptr |
35 | 35 | ||
36 | -- | check the error code and touch foreign ptr of vector arguments (if any) | 36 | -- | check the error code |
37 | check :: String -> [Vector a] -> IO Int -> IO () | 37 | check :: String -> IO Int -> IO () |
38 | check msg ls f = do | 38 | check msg f = do |
39 | err <- f | 39 | err <- f |
40 | when (err/=0) $ if err > 1024 | 40 | when (err/=0) $ if err > 1024 |
41 | then (error (msg++": "++errorCode err)) -- our errors | 41 | then (error (msg++": "++errorCode err)) -- our errors |
@@ -43,7 +43,6 @@ check msg ls f = do | |||
43 | ps <- gsl_strerror err | 43 | ps <- gsl_strerror err |
44 | s <- peekCString ps | 44 | s <- peekCString ps |
45 | error (msg++": "++s) | 45 | error (msg++": "++s) |
46 | mapM_ (touchForeignPtr . fptr) ls | ||
47 | return () | 46 | return () |
48 | 47 | ||
49 | -- | description of GSL error codes | 48 | -- | description of GSL error codes |
@@ -55,9 +54,14 @@ type Vc t s = Int -> Ptr t -> s | |||
55 | -- infixr 5 :> | 54 | -- infixr 5 :> |
56 | -- type t :> s = Vc t s | 55 | -- type t :> s = Vc t s |
57 | 56 | ||
58 | -- | adaptation of our vectors to be admitted by foreign functions: @f \/\/ vec v@ | 57 | --- | adaptation of our vectors to be admitted by foreign functions: @f \/\/ vec v@ |
59 | vec :: Vector t -> (Vc t s) -> s | 58 | --vec :: Vector t -> (Vc t s) -> s |
60 | vec v f = f (dim v) (ptr v) | 59 | --vec v f = f (dim v) (ptr v) |
60 | |||
61 | withVector (V n fp) f = withForeignPtr fp $ \p -> do | ||
62 | let v f = do | ||
63 | f n p | ||
64 | f v | ||
61 | 65 | ||
62 | -- | allocates memory for a new vector | 66 | -- | allocates memory for a new vector |
63 | createVector :: Storable a => Int -> IO (Vector a) | 67 | createVector :: Storable a => Int -> IO (Vector a) |
@@ -76,7 +80,8 @@ fromList :: Storable a => [a] -> Vector a | |||
76 | fromList l = unsafePerformIO $ do | 80 | fromList l = unsafePerformIO $ do |
77 | v <- createVector (length l) | 81 | v <- createVector (length l) |
78 | let f _ p = pokeArray p l >> return 0 | 82 | let f _ p = pokeArray p l >> return 0 |
79 | f // vec v // check "fromList" [] | 83 | withVector v $ \v -> |
84 | f // v // check "fromList" | ||
80 | return v | 85 | return v |
81 | 86 | ||
82 | safeRead v = unsafePerformIO . withForeignPtr (fptr v) | 87 | safeRead v = unsafePerformIO . withForeignPtr (fptr v) |
@@ -118,8 +123,9 @@ subVector k l (v@V {dim=n}) | |||
118 | | k<0 || k >= n || k+l > n || l < 0 = error "subVector out of range" | 123 | | k<0 || k >= n || k+l > n || l < 0 = error "subVector out of range" |
119 | | otherwise = unsafePerformIO $ do | 124 | | otherwise = unsafePerformIO $ do |
120 | r <- createVector l | 125 | r <- createVector l |
121 | let f = copyArray (ptr r) (advancePtr (ptr v) k) l >> return 0 | 126 | let f _ s _ d = copyArray d (advancePtr s k) l >> return 0 |
122 | f // check "subVector" [v,r] | 127 | ww2 withVector v withVector r $ \v r -> |
128 | f // v // r // check "subVector" | ||
123 | return r | 129 | return r |
124 | 130 | ||
125 | {- | Reads a vector position: | 131 | {- | Reads a vector position: |
@@ -144,12 +150,12 @@ join [] = error "joining zero vectors" | |||
144 | join as = unsafePerformIO $ do | 150 | join as = unsafePerformIO $ do |
145 | let tot = sum (map dim as) | 151 | let tot = sum (map dim as) |
146 | r@V {fptr = p} <- createVector tot | 152 | r@V {fptr = p} <- createVector tot |
147 | withForeignPtr p $ \_ -> | 153 | withForeignPtr p $ \ptr -> |
148 | joiner as tot (ptr r) | 154 | joiner as tot ptr |
149 | return r | 155 | return r |
150 | where joiner [] _ _ = return () | 156 | where joiner [] _ _ = return () |
151 | joiner (r@V {dim = n, fptr = b} : cs) _ p = do | 157 | joiner (r@V {dim = n, fptr = b} : cs) _ p = do |
152 | withForeignPtr b $ \_ -> copyArray p (ptr r) n | 158 | withForeignPtr b $ \pb -> copyArray p pb n |
153 | joiner cs 0 (advancePtr p n) | 159 | joiner cs 0 (advancePtr p n) |
154 | 160 | ||
155 | 161 | ||
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 | ||
27 | genfft code v = unsafePerformIO $ do | 27 | genfft 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 | ||
32 | foreign import ccall "gsl-aux.h fft" c_fft :: Int -> TCVCV | 33 | foreign 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 |
57 | foreign import ccall "gsl-aux.h eigensystemR" c_eigS :: TMVM | 58 | foreign 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 |
90 | foreign import ccall "gsl-aux.h eigensystemC" c_eigH :: TCMVCM | 92 | foreign 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 | |||
149 | qr' x = unsafePerformIO $ do | 152 | qr' 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 | |||
161 | qrPacked' x = unsafePerformIO $ do | 165 | qrPacked' 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 | ||
173 | unpackQR' (qr,tau) = unsafePerformIO $ do | 178 | unpackQR' (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 |
180 | foreign import ccall "gsl-aux.h QRunpack" c_qrUnpack :: TMVMM | 186 | foreign import ccall "gsl-aux.h QRunpack" c_qrUnpack :: TMVMM |
@@ -196,20 +202,22 @@ cholR :: Matrix Double -> Matrix Double | |||
196 | cholR = cholR' . cmat | 202 | cholR = cholR' . cmat |
197 | 203 | ||
198 | cholR' x = unsafePerformIO $ do | 204 | cholR' 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 | ||
203 | foreign import ccall "gsl-aux.h cholR" c_cholR :: TMM | 210 | foreign import ccall "gsl-aux.h cholR" c_cholR :: TMM |
204 | 211 | ||
205 | cholC :: Matrix (Complex Double) -> Matrix (Complex Double) | 212 | cholC :: Matrix (Complex Double) -> Matrix (Complex Double) |
206 | cholC = cholC' . cmat | 213 | cholC = cholC' . cmat |
207 | 214 | ||
208 | cholC' x = unsafePerformIO $ do | 215 | cholC' 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 | ||
213 | foreign import ccall "gsl-aux.h cholC" c_cholC :: TCMCM | 221 | foreign 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) | |||
223 | luSolveR' a b | 231 | luSolveR' 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) | |||
240 | luSolveC' a b | 249 | luSolveC' 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 | ||
257 | luRaux' x = unsafePerformIO $ do | 267 | luRaux' 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 | ||
270 | luCaux' x = unsafePerformIO $ do | 281 | luCaux' 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 | |||
26 | import Foreign | 26 | import Foreign |
27 | import Complex | 27 | import 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 | --------------------------------------------------------------------- |
171 | iv :: (Vector Double -> Double) -> (Int -> Ptr Double -> Double) | 172 | iv :: (Vector Double -> Double) -> (Int -> Ptr Double -> Double) |
172 | iv f n p = f (createV n copy "iv" []) where | 173 | iv 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" | |||
187 | aux_vTov :: (Vector Double -> Vector Double) -> (Int -> Ptr Double -> Ptr Double -> IO()) | 188 | aux_vTov :: (Vector Double -> Vector Double) -> (Int -> Ptr Double -> Ptr Double -> IO()) |
188 | aux_vTov f n p r = g where | 189 | aux_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 | ||
198 | createV n fun msg ptrs = unsafePerformIO $ do | 199 | |
200 | createV 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 | {- | ||
203 | createM r c fun msg ptrs = unsafePerformIO $ do | 207 | createM 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 | ||
208 | createMIO r c fun msg ptrs = do | 213 | createMIO 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 | |||
47 | polySolve' :: Vector Double -> Vector (Complex Double) | 47 | polySolve' :: Vector Double -> Vector (Complex Double) |
48 | polySolve' v | dim v > 1 = unsafePerformIO $ do | 48 | polySolve' 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 | |||
45 | createSFR :: Storable a => String -> (Ptr a -> IO Int) -> (a, a) | 45 | createSFR :: Storable a => String -> (Ptr a -> IO Int) -> (a, a) |
46 | createSFR s f = unsafePerformIO $ do | 46 | createSFR 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 | ||
74 | toScalarAux fun code v = unsafePerformIO $ do | 74 | toScalarAux 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 | ||
79 | vectorMapAux fun code v = unsafePerformIO $ do | 80 | vectorMapAux 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 | ||
84 | vectorMapValAux fun code val v = unsafePerformIO $ do | 86 | vectorMapValAux 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 | ||
91 | vectorZipAux fun code u v = unsafePerformIO $ do | 94 | vectorZipAux 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 | |||
177 | linearSolveSQAux f st a b | 182 | linearSolveSQAux 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 | ||
202 | linearSolveAux f st a b = unsafePerformIO $ do | 208 | linearSolveAux 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 | ||
252 | cholAux f st a = unsafePerformIO $ do | 259 | cholAux 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 | |||
270 | qrAux f st a = unsafePerformIO $ do | 278 | qrAux 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 | |||
292 | hessAux f st a = unsafePerformIO $ do | 300 | hessAux 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 | |||
313 | schurAux f st a = unsafePerformIO $ do | 322 | schurAux 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 | ||