summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--INSTALL.md8
-rw-r--r--packages/base/THANKS.md15
-rw-r--r--packages/base/hmatrix.cabal1
-rw-r--r--packages/base/src/Internal/Algorithms.hs32
-rw-r--r--packages/base/src/Internal/Devel.hs5
-rw-r--r--packages/base/src/Internal/Foreign.hs102
-rw-r--r--packages/base/src/Internal/LAPACK.hs54
-rw-r--r--packages/base/src/Internal/Matrix.hs35
-rw-r--r--packages/base/src/Internal/Sparse.hs5
-rw-r--r--packages/base/src/Internal/Util.hs4
-rw-r--r--packages/base/src/Internal/Vector.hs24
-rw-r--r--packages/base/src/Internal/Vectorized.hs34
-rw-r--r--packages/base/src/Numeric/LinearAlgebra.hs3
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Devel.hs8
-rw-r--r--packages/gsl/src/Numeric/GSL/ODE.hs2
15 files changed, 131 insertions, 201 deletions
diff --git a/INSTALL.md b/INSTALL.md
index 014d614..e48064d 100644
--- a/INSTALL.md
+++ b/INSTALL.md
@@ -42,11 +42,11 @@ the dlls contributed by Gilberto Camara available in [gsl-lapack-windows.zip][wi
42 42
433.a) In a msys shell: 433.a) In a msys shell:
44 44
45 $ cabal install hmatrix-0.13.1.0 --extra-lib-dir=${GSL} --extra-include-dir=${GSL} 45 $ cabal install hmatrix-0.13.1.0 --extra-lib-dir=${GSL} --extra-include-dirs=${GSL}
46 46
473.b) In a normal windows cmd: 473.b) In a normal windows cmd:
48 48
49 > cabal install --extra-lib-dir=%GSL% --extra-include-dir=%GSL% 49 > cabal install --extra-lib-dir=%GSL% --extra-include-dirs=%GSL%
50 50
51It may be necessary to put the dlls in the search path. 51It may be necessary to put the dlls in the search path.
52 52
@@ -65,7 +65,7 @@ using this method.
65 65
663) In a normal Windows cmd: 663) In a normal Windows cmd:
67 67
68 > cabal install --flags=openblas --extra-lib-dirs=C:\...\OpenBLAS\lib --extra-include-dir=C:\...\OpenBLAS\include 68 > cabal install --flags=openblas --extra-lib-dirs=C:\...\OpenBLAS\lib --extra-include-dirs=C:\...\OpenBLAS\include
69 69
70### Stack-based Windows build 70### Stack-based Windows build
71 71
@@ -79,7 +79,7 @@ Similar should be build under other OSes, like Linux and OSX.
79 79
803) Example in a normal Windows cmd for building hmatrix base lib: 803) Example in a normal Windows cmd for building hmatrix base lib:
81 81
82 > stack install hmatrix --flag hmatrix:openblas --extra-lib-dirs=C:\...\OpenBLAS\lib --extra-include-dir=C:\...\OpenBLAS\include 82 > stack install hmatrix --flag hmatrix:openblas --extra-lib-dirs=C:\...\OpenBLAS\lib --extra-include-dirs=C:\...\OpenBLAS\include
83 83
84## Tests ############################################### 84## Tests ###############################################
85 85
diff --git a/packages/base/THANKS.md b/packages/base/THANKS.md
index c0ef1a3..3da4177 100644
--- a/packages/base/THANKS.md
+++ b/packages/base/THANKS.md
@@ -212,6 +212,17 @@ module reorganization, monadic mapVectorM, and many other improvements.
212- Justin Le added NFData instances for Static types, added mapping and outer product 212- Justin Le added NFData instances for Static types, added mapping and outer product
213 methods to Domain, and many other functions to the Static module. 213 methods to Domain, and many other functions to the Static module.
214 214
215- Sidharth Kapur added Normed and numeric instances for several Static types, and 215- Sidharth Kapur added Normed and numeric instances for several Static types,
216fixed the CPP issue in cabal files. 216fixed the CPP issue in cabal files, and made many other contributions.
217
218- Matt Renaud improved the documentation.
219
220- Joshua Moerman fixed cabal/stack flags for windows.
221
222- Francesco Mazzoli, Niklas Hambüchen, Patrick Chilton, and Andras Slemmer
223 discovered a serious and subtle bug in the wrapper helpers causing memory corruption.
224 Andras Slemmer fixed the bug. Thank you all.
225
226- Kevin Slagle implemented thinQR and thinRQ, much faster than the original qr,
227 and added compactSVDTol.
217 228
diff --git a/packages/base/hmatrix.cabal b/packages/base/hmatrix.cabal
index faf84e8..89e5423 100644
--- a/packages/base/hmatrix.cabal
+++ b/packages/base/hmatrix.cabal
@@ -55,7 +55,6 @@ library
55 Internal.Devel 55 Internal.Devel
56 Internal.Vectorized 56 Internal.Vectorized
57 Internal.Matrix 57 Internal.Matrix
58 Internal.Foreign
59 Internal.ST 58 Internal.ST
60 Internal.IO 59 Internal.IO
61 Internal.Element 60 Internal.Element
diff --git a/packages/base/src/Internal/Algorithms.hs b/packages/base/src/Internal/Algorithms.hs
index c4f1a60..7ead5b4 100644
--- a/packages/base/src/Internal/Algorithms.hs
+++ b/packages/base/src/Internal/Algorithms.hs
@@ -30,6 +30,7 @@ import Internal.LAPACK as LAPACK
30import Internal.Numeric 30import Internal.Numeric
31import Data.List(foldl1') 31import Data.List(foldl1')
32import qualified Data.Array as A 32import qualified Data.Array as A
33import qualified Data.Vector.Storable as Vector
33import Internal.ST 34import Internal.ST
34import Internal.Vectorized(range) 35import Internal.Vectorized(range)
35import Control.DeepSeq 36import Control.DeepSeq
@@ -291,9 +292,13 @@ fromList [35.18264833189422,1.4769076999800903]
291 292
292-} 293-}
293compactSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) 294compactSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
294compactSVD m = (u', subVector 0 d s, v') where 295compactSVD = compactSVDTol 1
296
297-- | @compactSVDTol r@ is similar to 'compactSVD', but uses tolerance @tol=r*g*eps*(max rows cols)@ to distinguish nonzero singular values, where @g@ is the greatest singular value.
298compactSVDTol :: Field t => Double -> Matrix t -> (Matrix t, Vector Double, Matrix t)
299compactSVDTol r m = (u', subVector 0 d s, v') where
295 (u,s,v) = thinSVD m 300 (u,s,v) = thinSVD m
296 d = rankSVD (1*eps) m s `max` 1 301 d = rankSVD (r*eps) m s `max` 1
297 u' = takeColumns d u 302 u' = takeColumns d u
298 v' = takeColumns d v 303 v' = takeColumns d v
299 304
@@ -475,9 +480,14 @@ instance (NFData t, Numeric t) => NFData (QR t)
475-- | QR factorization. 480-- | QR factorization.
476-- 481--
477-- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular. 482-- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular.
483-- Note: the current implementation is very slow for large matrices. 'thinQR' is much faster.
478qr :: Field t => Matrix t -> (Matrix t, Matrix t) 484qr :: Field t => Matrix t -> (Matrix t, Matrix t)
479qr = {-# SCC "qr" #-} unpackQR . qr' 485qr = {-# SCC "qr" #-} unpackQR . qr'
480 486
487-- | A version of 'qr' which returns only the @min (rows m) (cols m)@ columns of @q@ and rows of @r@.
488thinQR :: Field t => Matrix t -> (Matrix t, Matrix t)
489thinQR = {-# SCC "thinQR" #-} thinUnpackQR . qr'
490
481-- | Compute the QR decomposition of a matrix in compact form. 491-- | Compute the QR decomposition of a matrix in compact form.
482qrRaw :: Field t => Matrix t -> QR t 492qrRaw :: Field t => Matrix t -> QR t
483qrRaw m = QR x v 493qrRaw m = QR x v
@@ -494,9 +504,17 @@ qrgr n (QR a t)
494-- | RQ factorization. 504-- | RQ factorization.
495-- 505--
496-- If @(r,q) = rq m@ then @m == r \<> q@, where q is unitary and r is upper triangular. 506-- If @(r,q) = rq m@ then @m == r \<> q@, where q is unitary and r is upper triangular.
507-- Note: the current implementation is very slow for large matrices. 'thinRQ' is much faster.
497rq :: Field t => Matrix t -> (Matrix t, Matrix t) 508rq :: Field t => Matrix t -> (Matrix t, Matrix t)
498rq m = {-# SCC "rq" #-} (r,q) where 509rq = {-# SCC "rq" #-} rqFromQR qr
499 (q',r') = qr $ trans $ rev1 m 510
511-- | A version of 'rq' which returns only the @min (rows m) (cols m)@ columns of @r@ and rows of @q@.
512thinRQ :: Field t => Matrix t -> (Matrix t, Matrix t)
513thinRQ = {-# SCC "thinQR" #-} rqFromQR thinQR
514
515rqFromQR :: Field t => (Matrix t -> (Matrix t, Matrix t)) -> Matrix t -> (Matrix t, Matrix t)
516rqFromQR qr0 m = (r,q) where
517 (q',r') = qr0 $ trans $ rev1 m
500 r = rev2 (trans r') 518 r = rev2 (trans r')
501 q = rev2 (trans q') 519 q = rev2 (trans q')
502 rev1 = flipud . fliprl 520 rev1 = flipud . fliprl
@@ -724,6 +742,12 @@ unpackQR (pq, tau) = {-# SCC "unpackQR" #-} (q,r)
724 hs = zipWith haussholder (toList tau) vs 742 hs = zipWith haussholder (toList tau) vs
725 q = foldl1' mXm hs 743 q = foldl1' mXm hs
726 744
745thinUnpackQR :: (Field t) => (Matrix t, Vector t) -> (Matrix t, Matrix t)
746thinUnpackQR (pq, tau) = (q, r)
747 where mn = uncurry min $ size pq
748 q = qrgr mn $ QR pq tau
749 r = fromRows $ zipWith (\i v -> Vector.replicate i 0 Vector.++ Vector.drop i v) [0..mn-1] (toRows pq)
750
727unpackHess :: (Field t) => (Matrix t -> (Matrix t,Vector t)) -> Matrix t -> (Matrix t, Matrix t) 751unpackHess :: (Field t) => (Matrix t -> (Matrix t,Vector t)) -> Matrix t -> (Matrix t, Matrix t)
728unpackHess hf m 752unpackHess hf m
729 | rows m == 1 = ((1><1)[1],m) 753 | rows m == 1 = ((1><1)[1],m)
diff --git a/packages/base/src/Internal/Devel.hs b/packages/base/src/Internal/Devel.hs
index 92b5604..3887663 100644
--- a/packages/base/src/Internal/Devel.hs
+++ b/packages/base/src/Internal/Devel.hs
@@ -80,8 +80,8 @@ class TransArray c
80 where 80 where
81 type Trans c b 81 type Trans c b
82 type TransRaw c b 82 type TransRaw c b
83 apply :: (Trans c b) -> c -> b 83 apply :: c -> (b -> IO r) -> (Trans c b) -> IO r
84 applyRaw :: (TransRaw c b) -> c -> b 84 applyRaw :: c -> (b -> IO r) -> (TransRaw c b) -> IO r
85 infixl 1 `apply`, `applyRaw` 85 infixl 1 `apply`, `applyRaw`
86 86
87instance Storable t => TransArray (Vector t) 87instance Storable t => TransArray (Vector t)
@@ -92,4 +92,3 @@ instance Storable t => TransArray (Vector t)
92 {-# INLINE apply #-} 92 {-# INLINE apply #-}
93 applyRaw = avec 93 applyRaw = avec
94 {-# INLINE applyRaw #-} 94 {-# INLINE applyRaw #-}
95
diff --git a/packages/base/src/Internal/Foreign.hs b/packages/base/src/Internal/Foreign.hs
deleted file mode 100644
index ea071a4..0000000
--- a/packages/base/src/Internal/Foreign.hs
+++ /dev/null
@@ -1,102 +0,0 @@
1{-# LANGUAGE MagicHash, UnboxedTuples #-}
2-- | FFI and hmatrix helpers.
3--
4-- Sample usage, to upload a perspective matrix to a shader.
5--
6-- @ glUniformMatrix4fv 0 1 (fromIntegral gl_TRUE) \`appMatrix\` perspective 0.01 100 (pi\/2) (4\/3)
7-- @
8--
9
10module Internal.Foreign
11 ( app
12 , appVector, appVectorLen
13 , appMatrix, appMatrixLen, appMatrixRaw, appMatrixRawLen
14 , unsafeMatrixToVector, unsafeMatrixToForeignPtr
15 ) where
16
17import Foreign.C.Types(CInt)
18import Internal.Vector
19import Internal.Matrix
20import qualified Data.Vector.Storable as S
21import Foreign (Ptr, ForeignPtr, Storable)
22import GHC.Base (IO(..), realWorld#)
23
24{-# INLINE unsafeInlinePerformIO #-}
25-- | If we use unsafePerformIO, it may not get inlined, so in a function that returns IO (which are all safe uses of app* in this module), there would be
26-- unecessary calls to unsafePerformIO or its internals.
27unsafeInlinePerformIO :: IO a -> a
28unsafeInlinePerformIO (IO f) = case f realWorld# of
29 (# _, x #) -> x
30
31{-# INLINE app #-}
32-- | Only useful since it is left associated with a precedence of 1, unlike 'Prelude.$', which is right associative.
33-- e.g.
34--
35-- @
36-- someFunction
37-- \`appMatrixLen\` m
38-- \`appVectorLen\` v
39-- \`app\` other
40-- \`app\` arguments
41-- \`app\` go here
42-- @
43--
44-- One could also write:
45--
46-- @
47-- (someFunction
48-- \`appMatrixLen\` m
49-- \`appVectorLen\` v)
50-- other
51-- arguments
52-- (go here)
53-- @
54--
55app :: (a -> b) -> a -> b
56app f = f
57
58{-# INLINE appVector #-}
59appVector :: Storable a => (Ptr a -> b) -> Vector a -> b
60appVector f x = unsafeInlinePerformIO (S.unsafeWith x (return . f))
61
62{-# INLINE appVectorLen #-}
63appVectorLen :: Storable a => (CInt -> Ptr a -> b) -> Vector a -> b
64appVectorLen f x = unsafeInlinePerformIO (S.unsafeWith x (return . f (fromIntegral (S.length x))))
65
66{-# INLINE appMatrix #-}
67appMatrix :: Element a => (Ptr a -> b) -> Matrix a -> b
68appMatrix f x = unsafeInlinePerformIO (S.unsafeWith (flatten x) (return . f))
69
70{-# INLINE appMatrixLen #-}
71appMatrixLen :: Element a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b
72appMatrixLen f x = unsafeInlinePerformIO (S.unsafeWith (flatten x) (return . f r c))
73 where
74 r = fromIntegral (rows x)
75 c = fromIntegral (cols x)
76
77{-# INLINE appMatrixRaw #-}
78appMatrixRaw :: Storable a => (Ptr a -> b) -> Matrix a -> b
79appMatrixRaw f x = unsafeInlinePerformIO (S.unsafeWith (xdat x) (return . f))
80
81{-# INLINE appMatrixRawLen #-}
82appMatrixRawLen :: Element a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b
83appMatrixRawLen f x = unsafeInlinePerformIO (S.unsafeWith (xdat x) (return . f r c))
84 where
85 r = fromIntegral (rows x)
86 c = fromIntegral (cols x)
87
88infixl 1 `app`
89infixl 1 `appVector`
90infixl 1 `appMatrix`
91infixl 1 `appMatrixRaw`
92
93{-# INLINE unsafeMatrixToVector #-}
94-- | This will disregard the order of the matrix, and simply return it as-is.
95-- If the order of the matrix is RowMajor, this function is identical to 'flatten'.
96unsafeMatrixToVector :: Matrix a -> Vector a
97unsafeMatrixToVector = xdat
98
99{-# INLINE unsafeMatrixToForeignPtr #-}
100unsafeMatrixToForeignPtr :: Storable a => Matrix a -> (ForeignPtr a, Int)
101unsafeMatrixToForeignPtr m = S.unsafeToForeignPtr0 (xdat m)
102
diff --git a/packages/base/src/Internal/LAPACK.hs b/packages/base/src/Internal/LAPACK.hs
index c2c140b..231109a 100644
--- a/packages/base/src/Internal/LAPACK.hs
+++ b/packages/base/src/Internal/LAPACK.hs
@@ -18,7 +18,7 @@ module Internal.LAPACK where
18 18
19import Internal.Devel 19import Internal.Devel
20import Internal.Vector 20import Internal.Vector
21import Internal.Matrix hiding ((#)) 21import Internal.Matrix hiding ((#), (#!))
22import Internal.Conversion 22import Internal.Conversion
23import Internal.Element 23import Internal.Element
24import Foreign.Ptr(nullPtr) 24import Foreign.Ptr(nullPtr)
@@ -28,10 +28,13 @@ import System.IO.Unsafe(unsafePerformIO)
28 28
29----------------------------------------------------------------------------------- 29-----------------------------------------------------------------------------------
30 30
31infixl 1 # 31infixr 1 #
32a # b = apply a b 32a # b = apply a b
33{-# INLINE (#) #-} 33{-# INLINE (#) #-}
34 34
35a #! b = a # b # id
36{-# INLINE (#!) #-}
37
35----------------------------------------------------------------------------------- 38-----------------------------------------------------------------------------------
36 39
37type TMMM t = t ::> t ::> t ::> Ok 40type TMMM t = t ::> t ::> t ::> Ok
@@ -56,7 +59,7 @@ multiplyAux f st a b = unsafePerformIO $ do
56 when (cols a /= rows b) $ error $ "inconsistent dimensions in matrix product "++ 59 when (cols a /= rows b) $ error $ "inconsistent dimensions in matrix product "++
57 show (rows a,cols a) ++ " x " ++ show (rows b, cols b) 60 show (rows a,cols a) ++ " x " ++ show (rows b, cols b)
58 s <- createMatrix ColumnMajor (rows a) (cols b) 61 s <- createMatrix ColumnMajor (rows a) (cols b)
59 f (isT a) (isT b) # (tt a) # (tt b) # s #| st 62 ((tt a) # (tt b) #! s) (f (isT a) (isT b)) #| st
60 return s 63 return s
61 64
62-- | Matrix product based on BLAS's /dgemm/. 65-- | Matrix product based on BLAS's /dgemm/.
@@ -80,7 +83,7 @@ multiplyI m a b = unsafePerformIO $ do
80 when (cols a /= rows b) $ error $ 83 when (cols a /= rows b) $ error $
81 "inconsistent dimensions in matrix product "++ shSize a ++ " x " ++ shSize b 84 "inconsistent dimensions in matrix product "++ shSize a ++ " x " ++ shSize b
82 s <- createMatrix ColumnMajor (rows a) (cols b) 85 s <- createMatrix ColumnMajor (rows a) (cols b)
83 c_multiplyI m # a # b # s #|"c_multiplyI" 86 (a # b #! s) (c_multiplyI m) #|"c_multiplyI"
84 return s 87 return s
85 88
86multiplyL :: Z -> Matrix Z -> Matrix Z -> Matrix Z 89multiplyL :: Z -> Matrix Z -> Matrix Z -> Matrix Z
@@ -88,7 +91,7 @@ multiplyL m a b = unsafePerformIO $ do
88 when (cols a /= rows b) $ error $ 91 when (cols a /= rows b) $ error $
89 "inconsistent dimensions in matrix product "++ shSize a ++ " x " ++ shSize b 92 "inconsistent dimensions in matrix product "++ shSize a ++ " x " ++ shSize b
90 s <- createMatrix ColumnMajor (rows a) (cols b) 93 s <- createMatrix ColumnMajor (rows a) (cols b)
91 c_multiplyL m # a # b # s #|"c_multiplyL" 94 (a # b #! s) (c_multiplyL m) #|"c_multiplyL"
92 return s 95 return s
93 96
94----------------------------------------------------------------------------- 97-----------------------------------------------------------------------------
@@ -121,7 +124,7 @@ svdAux f st x = unsafePerformIO $ do
121 u <- createMatrix ColumnMajor r r 124 u <- createMatrix ColumnMajor r r
122 s <- createVector (min r c) 125 s <- createVector (min r c)
123 v <- createMatrix ColumnMajor c c 126 v <- createMatrix ColumnMajor c c
124 f # a # u # s # v #| st 127 (a # u # s #! v) f #| st
125 return (u,s,v) 128 return (u,s,v)
126 where 129 where
127 r = rows x 130 r = rows x
@@ -149,7 +152,7 @@ thinSVDAux f st x = unsafePerformIO $ do
149 u <- createMatrix ColumnMajor r q 152 u <- createMatrix ColumnMajor r q
150 s <- createVector q 153 s <- createVector q
151 v <- createMatrix ColumnMajor q c 154 v <- createMatrix ColumnMajor q c
152 f # a # u # s # v #| st 155 (a # u # s #! v) f #| st
153 return (u,s,v) 156 return (u,s,v)
154 where 157 where
155 r = rows x 158 r = rows x
@@ -176,7 +179,7 @@ svCd = svAux zgesdd "svCd"
176svAux f st x = unsafePerformIO $ do 179svAux f st x = unsafePerformIO $ do
177 a <- copy ColumnMajor x 180 a <- copy ColumnMajor x
178 s <- createVector q 181 s <- createVector q
179 g # a # s #| st 182 (a #! s) g #| st
180 return s 183 return s
181 where 184 where
182 r = rows x 185 r = rows x
@@ -197,7 +200,7 @@ rightSVAux f st x = unsafePerformIO $ do
197 a <- copy ColumnMajor x 200 a <- copy ColumnMajor x
198 s <- createVector q 201 s <- createVector q
199 v <- createMatrix ColumnMajor c c 202 v <- createMatrix ColumnMajor c c
200 g # a # s # v #| st 203 (a # s #! v) g #| st
201 return (s,v) 204 return (s,v)
202 where 205 where
203 r = rows x 206 r = rows x
@@ -218,7 +221,7 @@ leftSVAux f st x = unsafePerformIO $ do
218 a <- copy ColumnMajor x 221 a <- copy ColumnMajor x
219 u <- createMatrix ColumnMajor r r 222 u <- createMatrix ColumnMajor r r
220 s <- createVector q 223 s <- createVector q
221 g # a # u # s #| st 224 (a # u #! s) g #| st
222 return (u,s) 225 return (u,s)
223 where 226 where
224 r = rows x 227 r = rows x
@@ -237,7 +240,7 @@ eigAux f st m = unsafePerformIO $ do
237 a <- copy ColumnMajor m 240 a <- copy ColumnMajor m
238 l <- createVector r 241 l <- createVector r
239 v <- createMatrix ColumnMajor r r 242 v <- createMatrix ColumnMajor r r
240 g # a # l # v #| st 243 (a # l #! v) g #| st
241 return (l,v) 244 return (l,v)
242 where 245 where
243 r = rows m 246 r = rows m
@@ -252,7 +255,7 @@ eigC = eigAux zgeev "eigC"
252eigOnlyAux f st m = unsafePerformIO $ do 255eigOnlyAux f st m = unsafePerformIO $ do
253 a <- copy ColumnMajor m 256 a <- copy ColumnMajor m
254 l <- createVector r 257 l <- createVector r
255 g # a # l #| st 258 (a #! l) g #| st
256 return l 259 return l
257 where 260 where
258 r = rows m 261 r = rows m
@@ -277,7 +280,7 @@ eigRaux m = unsafePerformIO $ do
277 a <- copy ColumnMajor m 280 a <- copy ColumnMajor m
278 l <- createVector r 281 l <- createVector r
279 v <- createMatrix ColumnMajor r r 282 v <- createMatrix ColumnMajor r r
280 g # a # l # v #| "eigR" 283 (a # l #! v) g #| "eigR"
281 return (l,v) 284 return (l,v)
282 where 285 where
283 r = rows m 286 r = rows m
@@ -305,7 +308,7 @@ eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR"
305eigSHAux f st m = unsafePerformIO $ do 308eigSHAux f st m = unsafePerformIO $ do
306 l <- createVector r 309 l <- createVector r
307 v <- copy ColumnMajor m 310 v <- copy ColumnMajor m
308 f # l # v #| st 311 (l #! v) f #| st
309 return (l,v) 312 return (l,v)
310 where 313 where
311 r = rows m 314 r = rows m
@@ -356,7 +359,7 @@ linearSolveSQAux g f st a b
356 | n1==n2 && n1==r = unsafePerformIO . g $ do 359 | n1==n2 && n1==r = unsafePerformIO . g $ do
357 a' <- copy ColumnMajor a 360 a' <- copy ColumnMajor a
358 s <- copy ColumnMajor b 361 s <- copy ColumnMajor b
359 f # a' # s #| st 362 (a' #! s) f #| st
360 return s 363 return s
361 | otherwise = error $ st ++ " of nonsquare matrix" 364 | otherwise = error $ st ++ " of nonsquare matrix"
362 where 365 where
@@ -387,7 +390,7 @@ foreign import ccall unsafe "cholSolveC_l" zpotrs :: C ::> C ::> Ok
387linearSolveSQAux2 g f st a b 390linearSolveSQAux2 g f st a b
388 | n1==n2 && n1==r = unsafePerformIO . g $ do 391 | n1==n2 && n1==r = unsafePerformIO . g $ do
389 s <- copy ColumnMajor b 392 s <- copy ColumnMajor b
390 f # a # s #| st 393 (a #! s) f #| st
391 return s 394 return s
392 | otherwise = error $ st ++ " of nonsquare matrix" 395 | otherwise = error $ st ++ " of nonsquare matrix"
393 where 396 where
@@ -415,7 +418,7 @@ linearSolveAux f st a b
415 a' <- copy ColumnMajor a 418 a' <- copy ColumnMajor a
416 r <- createMatrix ColumnMajor (max m n) nrhs 419 r <- createMatrix ColumnMajor (max m n) nrhs
417 setRect 0 0 b r 420 setRect 0 0 b r
418 f # a' # r #| st 421 (a' #! r) f #| st
419 return r 422 return r
420 | otherwise = error $ "different number of rows in linearSolve ("++st++")" 423 | otherwise = error $ "different number of rows in linearSolve ("++st++")"
421 where 424 where
@@ -458,7 +461,7 @@ foreign import ccall unsafe "chol_l_S" dpotrf :: R ::> Ok
458 461
459cholAux f st a = do 462cholAux f st a = do
460 r <- copy ColumnMajor a 463 r <- copy ColumnMajor a
461 f # r #| st 464 (r # id) f #| st
462 return r 465 return r
463 466
464-- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/. 467-- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/.
@@ -495,7 +498,7 @@ qrC = qrAux zgeqr2 "qrC"
495qrAux f st a = unsafePerformIO $ do 498qrAux f st a = unsafePerformIO $ do
496 r <- copy ColumnMajor a 499 r <- copy ColumnMajor a
497 tau <- createVector mn 500 tau <- createVector mn
498 f # tau # r #| st 501 (tau #! r) f #| st
499 return (r,tau) 502 return (r,tau)
500 where 503 where
501 m = rows a 504 m = rows a
@@ -514,7 +517,7 @@ qrgrC = qrgrAux zungqr "qrgrC"
514 517
515qrgrAux f st n (a, tau) = unsafePerformIO $ do 518qrgrAux f st n (a, tau) = unsafePerformIO $ do
516 res <- copy ColumnMajor (subMatrix (0,0) (rows a,n) a) 519 res <- copy ColumnMajor (subMatrix (0,0) (rows a,n) a)
517 f # (subVector 0 n tau') # res #| st 520 ((subVector 0 n tau') #! res) f #| st
518 return res 521 return res
519 where 522 where
520 tau' = vjoin [tau, constantD 0 n] 523 tau' = vjoin [tau, constantD 0 n]
@@ -534,7 +537,7 @@ hessC = hessAux zgehrd "hessC"
534hessAux f st a = unsafePerformIO $ do 537hessAux f st a = unsafePerformIO $ do
535 r <- copy ColumnMajor a 538 r <- copy ColumnMajor a
536 tau <- createVector (mn-1) 539 tau <- createVector (mn-1)
537 f # tau # r #| st 540 (tau #! r) f #| st
538 return (r,tau) 541 return (r,tau)
539 where 542 where
540 m = rows a 543 m = rows a
@@ -556,7 +559,7 @@ schurC = schurAux zgees "schurC"
556schurAux f st a = unsafePerformIO $ do 559schurAux f st a = unsafePerformIO $ do
557 u <- createMatrix ColumnMajor n n 560 u <- createMatrix ColumnMajor n n
558 s <- copy ColumnMajor a 561 s <- copy ColumnMajor a
559 f # u # s #| st 562 (u #! s) f #| st
560 return (u,s) 563 return (u,s)
561 where 564 where
562 n = rows a 565 n = rows a
@@ -576,7 +579,7 @@ luC = luAux zgetrf "luC"
576luAux f st a = unsafePerformIO $ do 579luAux f st a = unsafePerformIO $ do
577 lu <- copy ColumnMajor a 580 lu <- copy ColumnMajor a
578 piv <- createVector (min n m) 581 piv <- createVector (min n m)
579 f # piv # lu #| st 582 (piv #! lu) f #| st
580 return (lu, map (pred.round) (toList piv)) 583 return (lu, map (pred.round) (toList piv))
581 where 584 where
582 n = rows a 585 n = rows a
@@ -598,7 +601,7 @@ lusC a piv b = lusAux zgetrs "lusC" (fmat a) piv b
598lusAux f st a piv b 601lusAux f st a piv b
599 | n1==n2 && n2==n =unsafePerformIO $ do 602 | n1==n2 && n2==n =unsafePerformIO $ do
600 x <- copy ColumnMajor b 603 x <- copy ColumnMajor b
601 f # a # piv' # x #| st 604 (a # piv' #! x) f #| st
602 return x 605 return x
603 | otherwise = error st 606 | otherwise = error st
604 where 607 where
@@ -622,7 +625,7 @@ ldlC = ldlAux zhetrf "ldlC"
622ldlAux f st a = unsafePerformIO $ do 625ldlAux f st a = unsafePerformIO $ do
623 ldl <- copy ColumnMajor a 626 ldl <- copy ColumnMajor a
624 piv <- createVector (rows a) 627 piv <- createVector (rows a)
625 f # piv # ldl #| st 628 (piv #! ldl) f #| st
626 return (ldl, map (pred.round) (toList piv)) 629 return (ldl, map (pred.round) (toList piv))
627 630
628----------------------------------------------------------------------------------- 631-----------------------------------------------------------------------------------
@@ -637,4 +640,3 @@ ldlsR a piv b = lusAux dsytrs "ldlsR" (fmat a) piv b
637-- | Solve a complex linear system from a precomputed LDL decomposition ('ldlC'), using LAPACK's /zsytrs/. 640-- | Solve a complex linear system from a precomputed LDL decomposition ('ldlC'), using LAPACK's /zsytrs/.
638ldlsC :: Matrix (Complex Double) -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double) 641ldlsC :: Matrix (Complex Double) -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double)
639ldlsC a piv b = lusAux zsytrs "ldlsC" (fmat a) piv b 642ldlsC a piv b = lusAux zsytrs "ldlsC" (fmat a) piv b
640
diff --git a/packages/base/src/Internal/Matrix.hs b/packages/base/src/Internal/Matrix.hs
index c47c625..0135288 100644
--- a/packages/base/src/Internal/Matrix.hs
+++ b/packages/base/src/Internal/Matrix.hs
@@ -22,7 +22,7 @@ module Internal.Matrix where
22 22
23import Internal.Vector 23import Internal.Vector
24import Internal.Devel 24import Internal.Devel
25import Internal.Vectorized hiding ((#)) 25import Internal.Vectorized hiding ((#), (#!))
26import Foreign.Marshal.Alloc ( free ) 26import Foreign.Marshal.Alloc ( free )
27import Foreign.Marshal.Array(newArray) 27import Foreign.Marshal.Array(newArray)
28import Foreign.Ptr ( Ptr ) 28import Foreign.Ptr ( Ptr )
@@ -110,15 +110,15 @@ fmat m
110 110
111-- C-Haskell matrix adapters 111-- C-Haskell matrix adapters
112{-# INLINE amatr #-} 112{-# INLINE amatr #-}
113amatr :: Storable a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b 113amatr :: Storable a => Matrix a -> (f -> IO r) -> (CInt -> CInt -> Ptr a -> f) -> IO r
114amatr f x = inlinePerformIO (unsafeWith (xdat x) (return . f r c)) 114amatr x f g = unsafeWith (xdat x) (f . g r c)
115 where 115 where
116 r = fi (rows x) 116 r = fi (rows x)
117 c = fi (cols x) 117 c = fi (cols x)
118 118
119{-# INLINE amat #-} 119{-# INLINE amat #-}
120amat :: Storable a => (CInt -> CInt -> CInt -> CInt -> Ptr a -> b) -> Matrix a -> b 120amat :: Storable a => Matrix a -> (f -> IO r) -> (CInt -> CInt -> CInt -> CInt -> Ptr a -> f) -> IO r
121amat f x = inlinePerformIO (unsafeWith (xdat x) (return . f r c sr sc)) 121amat x f g = unsafeWith (xdat x) (f . g r c sr sc)
122 where 122 where
123 r = fi (rows x) 123 r = fi (rows x)
124 c = fi (cols x) 124 c = fi (cols x)
@@ -135,10 +135,13 @@ instance Storable t => TransArray (Matrix t)
135 applyRaw = amatr 135 applyRaw = amatr
136 {-# INLINE applyRaw #-} 136 {-# INLINE applyRaw #-}
137 137
138infixl 1 # 138infixr 1 #
139a # b = apply a b 139a # b = apply a b
140{-# INLINE (#) #-} 140{-# INLINE (#) #-}
141 141
142a #! b = a # b # id
143{-# INLINE (#!) #-}
144
142-------------------------------------------------------------------------------- 145--------------------------------------------------------------------------------
143 146
144copy ord m = extractR ord m 0 (idxs[0,rows m-1]) 0 (idxs[0,cols m-1]) 147copy ord m = extractR ord m 0 (idxs[0,rows m-1]) 0 (idxs[0,cols m-1])
@@ -426,7 +429,8 @@ extractAux f ord m moder vr modec vc = do
426 let nr = if moder == 0 then fromIntegral $ vr@>1 - vr@>0 + 1 else dim vr 429 let nr = if moder == 0 then fromIntegral $ vr@>1 - vr@>0 + 1 else dim vr
427 nc = if modec == 0 then fromIntegral $ vc@>1 - vc@>0 + 1 else dim vc 430 nc = if modec == 0 then fromIntegral $ vc@>1 - vc@>0 + 1 else dim vc
428 r <- createMatrix ord nr nc 431 r <- createMatrix ord nr nc
429 f moder modec # vr # vc # m # r #|"extract" 432 (vr # vc # m #! r) (f moder modec) #|"extract"
433
430 return r 434 return r
431 435
432type Extr x = CInt -> CInt -> CIdxs (CIdxs (OM x (OM x (IO CInt)))) 436type Extr x = CInt -> CInt -> CIdxs (CIdxs (OM x (OM x (IO CInt))))
@@ -440,7 +444,7 @@ foreign import ccall unsafe "extractL" c_extractL :: Extr Z
440 444
441--------------------------------------------------------------- 445---------------------------------------------------------------
442 446
443setRectAux f i j m r = f (fi i) (fi j) # m # r #|"setRect" 447setRectAux f i j m r = (m #! r) (f (fi i) (fi j)) #|"setRect"
444 448
445type SetRect x = I -> I -> x ::> x::> Ok 449type SetRect x = I -> I -> x ::> x::> Ok
446 450
@@ -455,7 +459,7 @@ foreign import ccall unsafe "setRectL" c_setRectL :: SetRect Z
455 459
456sortG f v = unsafePerformIO $ do 460sortG f v = unsafePerformIO $ do
457 r <- createVector (dim v) 461 r <- createVector (dim v)
458 f # v # r #|"sortG" 462 (v #! r) f #|"sortG"
459 return r 463 return r
460 464
461sortIdxD = sortG c_sort_indexD 465sortIdxD = sortG c_sort_indexD
@@ -482,7 +486,7 @@ foreign import ccall unsafe "sort_valuesL" c_sort_valL :: Z :> Z :> Ok
482 486
483compareG f u v = unsafePerformIO $ do 487compareG f u v = unsafePerformIO $ do
484 r <- createVector (dim v) 488 r <- createVector (dim v)
485 f # u # v # r #|"compareG" 489 (u # v #! r) f #|"compareG"
486 return r 490 return r
487 491
488compareD = compareG c_compareD 492compareD = compareG c_compareD
@@ -499,7 +503,7 @@ foreign import ccall unsafe "compareL" c_compareL :: Z :> Z :> I :> Ok
499 503
500selectG f c u v w = unsafePerformIO $ do 504selectG f c u v w = unsafePerformIO $ do
501 r <- createVector (dim v) 505 r <- createVector (dim v)
502 f # c # u # v # w # r #|"selectG" 506 (c # u # v # w #! r) f #|"selectG"
503 return r 507 return r
504 508
505selectD = selectG c_selectD 509selectD = selectG c_selectD
@@ -522,7 +526,7 @@ foreign import ccall unsafe "chooseL" c_selectL :: Sel Z
522 526
523remapG f i j m = unsafePerformIO $ do 527remapG f i j m = unsafePerformIO $ do
524 r <- createMatrix RowMajor (rows i) (cols i) 528 r <- createMatrix RowMajor (rows i) (cols i)
525 f # i # j # m # r #|"remapG" 529 (i # j # m #! r) f #|"remapG"
526 return r 530 return r
527 531
528remapD = remapG c_remapD 532remapD = remapG c_remapD
@@ -545,7 +549,7 @@ foreign import ccall unsafe "remapL" c_remapL :: Rem Z
545 549
546rowOpAux f c x i1 i2 j1 j2 m = do 550rowOpAux f c x i1 i2 j1 j2 m = do
547 px <- newArray [x] 551 px <- newArray [x]
548 f (fi c) px (fi i1) (fi i2) (fi j1) (fi j2) # m #|"rowOp" 552 (m # id) (f (fi c) px (fi i1) (fi i2) (fi j1) (fi j2)) #|"rowOp"
549 free px 553 free px
550 554
551type RowOp x = CInt -> Ptr x -> CInt -> CInt -> CInt -> CInt -> x ::> Ok 555type RowOp x = CInt -> Ptr x -> CInt -> CInt -> CInt -> CInt -> x ::> Ok
@@ -561,7 +565,7 @@ foreign import ccall unsafe "rowop_mod_int64_t" c_rowOpML :: Z -> RowOp Z
561 565
562-------------------------------------------------------------------------------- 566--------------------------------------------------------------------------------
563 567
564gemmg f v m1 m2 m3 = f # v # m1 # m2 # m3 #|"gemmg" 568gemmg f v m1 m2 m3 = (v # m1 # m2 #! m3) f #|"gemmg"
565 569
566type Tgemm x = x :> x ::> x ::> x ::> Ok 570type Tgemm x = x :> x ::> x ::> x ::> Ok
567 571
@@ -589,10 +593,9 @@ saveMatrix
589saveMatrix name format m = do 593saveMatrix name format m = do
590 cname <- newCString name 594 cname <- newCString name
591 cformat <- newCString format 595 cformat <- newCString format
592 c_saveMatrix cname cformat # m #|"saveMatrix" 596 (m # id) (c_saveMatrix cname cformat) #|"saveMatrix"
593 free cname 597 free cname
594 free cformat 598 free cformat
595 return () 599 return ()
596 600
597-------------------------------------------------------------------------------- 601--------------------------------------------------------------------------------
598
diff --git a/packages/base/src/Internal/Sparse.hs b/packages/base/src/Internal/Sparse.hs
index 1604e7e..1ff3f57 100644
--- a/packages/base/src/Internal/Sparse.hs
+++ b/packages/base/src/Internal/Sparse.hs
@@ -144,13 +144,13 @@ gmXv :: GMatrix -> Vector Double -> Vector Double
144gmXv SparseR { gmCSR = CSR{..}, .. } v = unsafePerformIO $ do 144gmXv SparseR { gmCSR = CSR{..}, .. } v = unsafePerformIO $ do
145 dim v /= nCols ~!~ printf "gmXv (CSR): incorrect sizes: (%d,%d) x %d" nRows nCols (dim v) 145 dim v /= nCols ~!~ printf "gmXv (CSR): incorrect sizes: (%d,%d) x %d" nRows nCols (dim v)
146 r <- createVector nRows 146 r <- createVector nRows
147 c_smXv # csrVals # csrCols # csrRows # v # r #|"CSRXv" 147 (csrVals # csrCols # csrRows # v #! r) c_smXv #|"CSRXv"
148 return r 148 return r
149 149
150gmXv SparseC { gmCSC = CSC{..}, .. } v = unsafePerformIO $ do 150gmXv SparseC { gmCSC = CSC{..}, .. } v = unsafePerformIO $ do
151 dim v /= nCols ~!~ printf "gmXv (CSC): incorrect sizes: (%d,%d) x %d" nRows nCols (dim v) 151 dim v /= nCols ~!~ printf "gmXv (CSC): incorrect sizes: (%d,%d) x %d" nRows nCols (dim v)
152 r <- createVector nRows 152 r <- createVector nRows
153 c_smTXv # cscVals # cscRows # cscCols # v # r #|"CSCXv" 153 (cscVals # cscRows # cscCols # v #! r) c_smTXv #|"CSCXv"
154 return r 154 return r
155 155
156gmXv Diag{..} v 156gmXv Diag{..} v
@@ -211,4 +211,3 @@ instance Transposable GMatrix GMatrix
211 tr (Diag v n m) = Diag v m n 211 tr (Diag v n m) = Diag v m n
212 tr (Dense a n m) = Dense (tr a) m n 212 tr (Dense a n m) = Dense (tr a) m n
213 tr' = tr 213 tr' = tr
214
diff --git a/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs
index cf42961..7f52a9c 100644
--- a/packages/base/src/Internal/Util.hs
+++ b/packages/base/src/Internal/Util.hs
@@ -255,6 +255,7 @@ norm :: Vector Double -> Double
255-- ^ 2-norm of real vector 255-- ^ 2-norm of real vector
256norm = pnorm PNorm2 256norm = pnorm PNorm2
257 257
258-- | p-norm for vectors, operator norm for matrices
258class Normed a 259class Normed a
259 where 260 where
260 norm_0 :: a -> R 261 norm_0 :: a -> R
@@ -319,10 +320,11 @@ instance Normed (Vector (Complex Float))
319 norm_2 = norm_2 . double 320 norm_2 = norm_2 . double
320 norm_Inf = norm_Inf . double 321 norm_Inf = norm_Inf . double
321 322
322 323-- | Frobenius norm (Schatten p-norm with p=2)
323norm_Frob :: (Normed (Vector t), Element t) => Matrix t -> R 324norm_Frob :: (Normed (Vector t), Element t) => Matrix t -> R
324norm_Frob = norm_2 . flatten 325norm_Frob = norm_2 . flatten
325 326
327-- | Sum of singular values (Schatten p-norm with p=1)
326norm_nuclear :: Field t => Matrix t -> R 328norm_nuclear :: Field t => Matrix t -> R
327norm_nuclear = sumElements . singularValues 329norm_nuclear = sumElements . singularValues
328 330
diff --git a/packages/base/src/Internal/Vector.hs b/packages/base/src/Internal/Vector.hs
index b4e235c..c4a310d 100644
--- a/packages/base/src/Internal/Vector.hs
+++ b/packages/base/src/Internal/Vector.hs
@@ -66,9 +66,8 @@ dim = Vector.length
66 66
67-- C-Haskell vector adapter 67-- C-Haskell vector adapter
68{-# INLINE avec #-} 68{-# INLINE avec #-}
69avec :: Storable a => (CInt -> Ptr a -> b) -> Vector a -> b 69avec :: Storable a => Vector a -> (f -> IO r) -> ((CInt -> Ptr a -> f) -> IO r)
70avec f v = inlinePerformIO (unsafeWith v (return . f (fromIntegral (Vector.length v)))) 70avec v f g = unsafeWith v $ \ptr -> f (g (fromIntegral (Vector.length v)) ptr)
71infixl 1 `avec`
72 71
73-- allocates memory for a new vector 72-- allocates memory for a new vector
74createVector :: Storable a => Int -> IO (Vector a) 73createVector :: Storable a => Int -> IO (Vector a)
@@ -199,7 +198,7 @@ takesV ms w | sum ms > dim w = error $ "takesV " ++ show ms ++ " on dim = " ++ (
199 198
200--------------------------------------------------------------- 199---------------------------------------------------------------
201 200
202-- | transforms a complex vector into a real vector with alternating real and imaginary parts 201-- | transforms a complex vector into a real vector with alternating real and imaginary parts
203asReal :: (RealFloat a, Storable a) => Vector (Complex a) -> Vector a 202asReal :: (RealFloat a, Storable a) => Vector (Complex a) -> Vector a
204asReal v = unsafeFromForeignPtr (castForeignPtr fp) (2*i) (2*n) 203asReal v = unsafeFromForeignPtr (castForeignPtr fp) (2*i) (2*n)
205 where (fp,i,n) = unsafeToForeignPtr v 204 where (fp,i,n) = unsafeToForeignPtr v
@@ -244,7 +243,7 @@ zipVectorWith f u v = unsafePerformIO $ do
244{-# INLINE zipVectorWith #-} 243{-# INLINE zipVectorWith #-}
245 244
246-- | unzipWith for Vectors 245-- | unzipWith for Vectors
247unzipVectorWith :: (Storable (a,b), Storable c, Storable d) 246unzipVectorWith :: (Storable (a,b), Storable c, Storable d)
248 => ((a,b) -> (c,d)) -> Vector (a,b) -> (Vector c,Vector d) 247 => ((a,b) -> (c,d)) -> Vector (a,b) -> (Vector c,Vector d)
249unzipVectorWith f u = unsafePerformIO $ do 248unzipVectorWith f u = unsafePerformIO $ do
250 let n = dim u 249 let n = dim u
@@ -255,7 +254,7 @@ unzipVectorWith f u = unsafePerformIO $ do
255 unsafeWith w $ \pw -> do 254 unsafeWith w $ \pw -> do
256 let go (-1) = return () 255 let go (-1) = return ()
257 go !k = do z <- peekElemOff pu k 256 go !k = do z <- peekElemOff pu k
258 let (x,y) = f z 257 let (x,y) = f z
259 pokeElemOff pv k x 258 pokeElemOff pv k x
260 pokeElemOff pw k y 259 pokeElemOff pw k y
261 go (k-1) 260 go (k-1)
@@ -303,11 +302,11 @@ mapVectorM f v = do
303 return w 302 return w
304 where mapVectorM' w' !k !t 303 where mapVectorM' w' !k !t
305 | k == t = do 304 | k == t = do
306 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 305 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
307 y <- f x 306 y <- f x
308 return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y 307 return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
309 | otherwise = do 308 | otherwise = do
310 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 309 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
311 y <- f x 310 y <- f x
312 _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y 311 _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
313 mapVectorM' w' (k+1) t 312 mapVectorM' w' (k+1) t
@@ -322,7 +321,7 @@ mapVectorM_ f v = do
322 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 321 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
323 f x 322 f x
324 | otherwise = do 323 | otherwise = do
325 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 324 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
326 _ <- f x 325 _ <- f x
327 mapVectorM' (k+1) t 326 mapVectorM' (k+1) t
328{-# INLINE mapVectorM_ #-} 327{-# INLINE mapVectorM_ #-}
@@ -336,11 +335,11 @@ mapVectorWithIndexM f v = do
336 return w 335 return w
337 where mapVectorM' w' !k !t 336 where mapVectorM' w' !k !t
338 | k == t = do 337 | k == t = do
339 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 338 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
340 y <- f k x 339 y <- f k x
341 return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y 340 return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
342 | otherwise = do 341 | otherwise = do
343 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 342 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
344 y <- f k x 343 y <- f k x
345 _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y 344 _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
346 mapVectorM' w' (k+1) t 345 mapVectorM' w' (k+1) t
@@ -355,7 +354,7 @@ mapVectorWithIndexM_ f v = do
355 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 354 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
356 f k x 355 f k x
357 | otherwise = do 356 | otherwise = do
358 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 357 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
359 _ <- f k x 358 _ <- f k x
360 mapVectorM' (k+1) t 359 mapVectorM' (k+1) t
361{-# INLINE mapVectorWithIndexM_ #-} 360{-# INLINE mapVectorWithIndexM_ #-}
@@ -454,4 +453,3 @@ unzipVector :: (Storable a, Storable b, Storable (a,b)) => Vector (a,b) -> (Vect
454unzipVector = unzipVectorWith id 453unzipVector = unzipVectorWith id
455 454
456------------------------------------------------------------------- 455-------------------------------------------------------------------
457
diff --git a/packages/base/src/Internal/Vectorized.hs b/packages/base/src/Internal/Vectorized.hs
index 03bcf90..a410bb2 100644
--- a/packages/base/src/Internal/Vectorized.hs
+++ b/packages/base/src/Internal/Vectorized.hs
@@ -27,10 +27,13 @@ import Foreign.C.String
27import System.IO.Unsafe(unsafePerformIO) 27import System.IO.Unsafe(unsafePerformIO)
28import Control.Monad(when) 28import Control.Monad(when)
29 29
30infixl 1 # 30infixr 1 #
31a # b = applyRaw a b 31a # b = applyRaw a b
32{-# INLINE (#) #-} 32{-# INLINE (#) #-}
33 33
34a #! b = a # b # id
35{-# INLINE (#!) #-}
36
34fromei x = fromIntegral (fromEnum x) :: CInt 37fromei x = fromIntegral (fromEnum x) :: CInt
35 38
36data FunCodeV = Sin 39data FunCodeV = Sin
@@ -103,7 +106,7 @@ sumL m = sumg (c_sumL m)
103 106
104sumg f x = unsafePerformIO $ do 107sumg f x = unsafePerformIO $ do
105 r <- createVector 1 108 r <- createVector 1
106 f # x # r #| "sum" 109 (x #! r) f #| "sum"
107 return $ r @> 0 110 return $ r @> 0
108 111
109type TVV t = t :> t :> Ok 112type TVV t = t :> t :> Ok
@@ -139,7 +142,7 @@ prodL = prodg . c_prodL
139 142
140prodg f x = unsafePerformIO $ do 143prodg f x = unsafePerformIO $ do
141 r <- createVector 1 144 r <- createVector 1
142 f # x # r #| "prod" 145 (x #! r) f #| "prod"
143 return $ r @> 0 146 return $ r @> 0
144 147
145 148
@@ -154,24 +157,24 @@ foreign import ccall unsafe "prodL" c_prodL :: Z -> TVV Z
154 157
155toScalarAux fun code v = unsafePerformIO $ do 158toScalarAux fun code v = unsafePerformIO $ do
156 r <- createVector 1 159 r <- createVector 1
157 fun (fromei code) # v # r #|"toScalarAux" 160 (v #! r) (fun (fromei code)) #|"toScalarAux"
158 return (r @> 0) 161 return (r @> 0)
159 162
160vectorMapAux fun code v = unsafePerformIO $ do 163vectorMapAux fun code v = unsafePerformIO $ do
161 r <- createVector (dim v) 164 r <- createVector (dim v)
162 fun (fromei code) # v # r #|"vectorMapAux" 165 (v #! r) (fun (fromei code)) #|"vectorMapAux"
163 return r 166 return r
164 167
165vectorMapValAux fun code val v = unsafePerformIO $ do 168vectorMapValAux fun code val v = unsafePerformIO $ do
166 r <- createVector (dim v) 169 r <- createVector (dim v)
167 pval <- newArray [val] 170 pval <- newArray [val]
168 fun (fromei code) pval # v # r #|"vectorMapValAux" 171 (v #! r) (fun (fromei code) pval) #|"vectorMapValAux"
169 free pval 172 free pval
170 return r 173 return r
171 174
172vectorZipAux fun code u v = unsafePerformIO $ do 175vectorZipAux fun code u v = unsafePerformIO $ do
173 r <- createVector (dim u) 176 r <- createVector (dim u)
174 fun (fromei code) # u # v # r #|"vectorZipAux" 177 (u # v #! r) (fun (fromei code)) #|"vectorZipAux"
175 return r 178 return r
176 179
177--------------------------------------------------------------------- 180---------------------------------------------------------------------
@@ -368,7 +371,7 @@ randomVector :: Seed
368 -> Vector Double 371 -> Vector Double
369randomVector seed dist n = unsafePerformIO $ do 372randomVector seed dist n = unsafePerformIO $ do
370 r <- createVector n 373 r <- createVector n
371 c_random_vector (fi seed) ((fi.fromEnum) dist) # r #|"randomVector" 374 (r # id) (c_random_vector (fi seed) ((fi.fromEnum) dist)) #|"randomVector"
372 return r 375 return r
373 376
374foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> Double :> Ok 377foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> Double :> Ok
@@ -377,7 +380,7 @@ foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> D
377 380
378roundVector v = unsafePerformIO $ do 381roundVector v = unsafePerformIO $ do
379 r <- createVector (dim v) 382 r <- createVector (dim v)
380 c_round_vector # v # r #|"roundVector" 383 (v #! r) c_round_vector #|"roundVector"
381 return r 384 return r
382 385
383foreign import ccall unsafe "round_vector" c_round_vector :: TVV Double 386foreign import ccall unsafe "round_vector" c_round_vector :: TVV Double
@@ -391,7 +394,7 @@ foreign import ccall unsafe "round_vector" c_round_vector :: TVV Double
391range :: Int -> Vector I 394range :: Int -> Vector I
392range n = unsafePerformIO $ do 395range n = unsafePerformIO $ do
393 r <- createVector n 396 r <- createVector n
394 c_range_vector # r #|"range" 397 (r # id) c_range_vector #|"range"
395 return r 398 return r
396 399
397foreign import ccall unsafe "range_vector" c_range_vector :: CInt :> Ok 400foreign import ccall unsafe "range_vector" c_range_vector :: CInt :> Ok
@@ -431,7 +434,7 @@ long2intV = tog c_long2int
431 434
432tog f v = unsafePerformIO $ do 435tog f v = unsafePerformIO $ do
433 r <- createVector (dim v) 436 r <- createVector (dim v)
434 f # v # r #|"tog" 437 (v #! r) f #|"tog"
435 return r 438 return r
436 439
437foreign import ccall unsafe "float2double" c_float2double :: Float :> Double :> Ok 440foreign import ccall unsafe "float2double" c_float2double :: Float :> Double :> Ok
@@ -450,7 +453,7 @@ foreign import ccall unsafe "long2int" c_long2int :: Z :> I :> Ok
450 453
451stepg f v = unsafePerformIO $ do 454stepg f v = unsafePerformIO $ do
452 r <- createVector (dim v) 455 r <- createVector (dim v)
453 f # v # r #|"step" 456 (v #! r) f #|"step"
454 return r 457 return r
455 458
456stepD :: Vector Double -> Vector Double 459stepD :: Vector Double -> Vector Double
@@ -475,7 +478,7 @@ foreign import ccall unsafe "stepL" c_stepL :: TVV Z
475 478
476conjugateAux fun x = unsafePerformIO $ do 479conjugateAux fun x = unsafePerformIO $ do
477 v <- createVector (dim x) 480 v <- createVector (dim x)
478 fun # x # v #|"conjugateAux" 481 (x #! v) fun #|"conjugateAux"
479 return v 482 return v
480 483
481conjugateQ :: Vector (Complex Float) -> Vector (Complex Float) 484conjugateQ :: Vector (Complex Float) -> Vector (Complex Float)
@@ -493,7 +496,7 @@ cloneVector v = do
493 let n = dim v 496 let n = dim v
494 r <- createVector n 497 r <- createVector n
495 let f _ s _ d = copyArray d s n >> return 0 498 let f _ s _ d = copyArray d s n >> return 0
496 f # v # r #|"cloneVector" 499 (v #! r) f #|"cloneVector"
497 return r 500 return r
498 501
499-------------------------------------------------------------------------------- 502--------------------------------------------------------------------------------
@@ -501,7 +504,7 @@ cloneVector v = do
501constantAux fun x n = unsafePerformIO $ do 504constantAux fun x n = unsafePerformIO $ do
502 v <- createVector n 505 v <- createVector n
503 px <- newArray [x] 506 px <- newArray [x]
504 fun px # v #|"constantAux" 507 (v # id) (fun px) #|"constantAux"
505 free px 508 free px
506 return v 509 return v
507 510
@@ -515,4 +518,3 @@ foreign import ccall unsafe "constantI" cconstantI :: TConst CInt
515foreign import ccall unsafe "constantL" cconstantL :: TConst Z 518foreign import ccall unsafe "constantL" cconstantL :: TConst Z
516 519
517---------------------------------------------------------------------- 520----------------------------------------------------------------------
518
diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs
index b5bacf4..626ee52 100644
--- a/packages/base/src/Numeric/LinearAlgebra.hs
+++ b/packages/base/src/Numeric/LinearAlgebra.hs
@@ -117,6 +117,7 @@ module Numeric.LinearAlgebra (
117 svd, 117 svd,
118 thinSVD, 118 thinSVD,
119 compactSVD, 119 compactSVD,
120 compactSVDTol,
120 singularValues, 121 singularValues,
121 leftSV, rightSV, 122 leftSV, rightSV,
122 123
@@ -126,7 +127,7 @@ module Numeric.LinearAlgebra (
126 geigSH, 127 geigSH,
127 128
128 -- * QR 129 -- * QR
129 qr, rq, qrRaw, qrgr, 130 qr, thinQR, rq, thinRQ, qrRaw, qrgr,
130 131
131 -- * Cholesky 132 -- * Cholesky
132 chol, mbChol, 133 chol, mbChol,
diff --git a/packages/base/src/Numeric/LinearAlgebra/Devel.hs b/packages/base/src/Numeric/LinearAlgebra/Devel.hs
index 57a68e7..941b597 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Devel.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Devel.hs
@@ -12,13 +12,6 @@ The library can be easily extended using the tools in this module.
12-------------------------------------------------------------------------------- 12--------------------------------------------------------------------------------
13 13
14module Numeric.LinearAlgebra.Devel( 14module Numeric.LinearAlgebra.Devel(
15 -- * FFI helpers
16 -- | Sample usage, to upload a perspective matrix to a shader.
17 --
18 -- @ glUniformMatrix4fv 0 1 (fromIntegral gl_TRUE) \`appMatrix\` perspective 0.01 100 (pi\/2) (4\/3)
19 -- @
20 module Internal.Foreign,
21
22 -- * FFI tools 15 -- * FFI tools
23 -- | See @examples/devel@ in the repository. 16 -- | See @examples/devel@ in the repository.
24 17
@@ -66,7 +59,6 @@ module Numeric.LinearAlgebra.Devel(
66 59
67) where 60) where
68 61
69import Internal.Foreign
70import Internal.Devel 62import Internal.Devel
71import Internal.ST 63import Internal.ST
72import Internal.Vector 64import Internal.Vector
diff --git a/packages/gsl/src/Numeric/GSL/ODE.hs b/packages/gsl/src/Numeric/GSL/ODE.hs
index 9e52873..c4fab85 100644
--- a/packages/gsl/src/Numeric/GSL/ODE.hs
+++ b/packages/gsl/src/Numeric/GSL/ODE.hs
@@ -17,7 +17,7 @@ A simple example:
17@ 17@
18import Numeric.GSL.ODE 18import Numeric.GSL.ODE
19import Numeric.LinearAlgebra 19import Numeric.LinearAlgebra
20import Numeric.LinearAlgebra.Util(mplot) 20import Graphics.Plot(mplot)
21 21
22xdot t [x,v] = [v, -0.95*x - 0.1*v] 22xdot t [x,v] = [v, -0.95*x - 0.1*v]
23 23
generated by cgit v1.2.3 (git 2.39.1) at 2024-08-29 10:48:14 +0000