summaryrefslogtreecommitdiff
path: root/lib/Data
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Data')
-rw-r--r--lib/Data/Packed.hs11
-rw-r--r--lib/Data/Packed/Internal/Matrix.hs32
-rw-r--r--lib/Data/Packed/Internal/Vector.hs1
-rw-r--r--lib/Data/Packed/Internal/auxi.c28
-rw-r--r--lib/Data/Packed/Internal/auxi.h3
-rw-r--r--lib/Data/Packed/ST.hs20
-rw-r--r--lib/Data/Packed/Vector.hs9
7 files changed, 35 insertions, 69 deletions
diff --git a/lib/Data/Packed.hs b/lib/Data/Packed.hs
index bb25c26..7d6d200 100644
--- a/lib/Data/Packed.hs
+++ b/lib/Data/Packed.hs
@@ -24,7 +24,7 @@ module Data.Packed (
24import Data.Packed.Vector 24import Data.Packed.Vector
25import Data.Packed.Matrix 25import Data.Packed.Matrix
26import Data.Complex 26import Data.Complex
27import Data.Packed.Internal(fromComplex,toComplex,comp,conj) 27import Data.Packed.Internal(fromComplex,toComplex,conj)
28 28
29-- | conversion utilities 29-- | conversion utilities
30class (Element e) => Container c e where 30class (Element e) => Container c e where
@@ -38,7 +38,7 @@ class (Element e) => Container c e where
38instance Container Vector Double where 38instance Container Vector Double where
39 toComplex = Data.Packed.Internal.toComplex 39 toComplex = Data.Packed.Internal.toComplex
40 fromComplex = Data.Packed.Internal.fromComplex 40 fromComplex = Data.Packed.Internal.fromComplex
41 comp = Data.Packed.Internal.comp 41 comp = internalcomp
42 conj = Data.Packed.Internal.conj 42 conj = Data.Packed.Internal.conj
43 real = id 43 real = id
44 complex = Data.Packed.comp 44 complex = Data.Packed.comp
@@ -56,7 +56,7 @@ instance Container Matrix Double where
56 fromComplex z = (reshape c r, reshape c i) 56 fromComplex z = (reshape c r, reshape c i)
57 where (r,i) = Data.Packed.fromComplex (flatten z) 57 where (r,i) = Data.Packed.fromComplex (flatten z)
58 c = cols z 58 c = cols z
59 comp = liftMatrix Data.Packed.Internal.comp 59 comp = liftMatrix internalcomp
60 conj = liftMatrix Data.Packed.Internal.conj 60 conj = liftMatrix Data.Packed.Internal.conj
61 real = id 61 real = id
62 complex = Data.Packed.comp 62 complex = Data.Packed.comp
@@ -68,3 +68,8 @@ instance Container Matrix (Complex Double) where
68 conj = undefined 68 conj = undefined
69 real = Data.Packed.comp 69 real = Data.Packed.comp
70 complex = id 70 complex = id
71
72
73-- | converts a real vector into a complex representation (with zero imaginary parts)
74internalcomp :: Vector Double -> Vector (Complex Double)
75internalcomp v = Data.Packed.Internal.toComplex (v,constant 0 (dim v))
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 51fb6f8..477b453 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -208,19 +208,16 @@ compat m1 m2 = rows m1 == rows m2 && cols m1 == cols m2
208 208
209-- | Optimized matrix computations are provided for elements in the Element class. 209-- | Optimized matrix computations are provided for elements in the Element class.
210class (Storable a, Floating a) => Element a where 210class (Storable a, Floating a) => Element a where
211 constantD :: a -> Int -> Vector a
212 transdata :: Int -> Vector a -> Int -> Vector a 211 transdata :: Int -> Vector a -> Int -> Vector a
213 subMatrixD :: (Int,Int) -- ^ (r0,c0) starting position 212 subMatrixD :: (Int,Int) -- ^ (r0,c0) starting position
214 -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix 213 -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
215 -> Matrix a -> Matrix a 214 -> Matrix a -> Matrix a
216 215
217instance Element Double where 216instance Element Double where
218 constantD = constantR
219 transdata = transdataR 217 transdata = transdataR
220 subMatrixD = subMatrixR 218 subMatrixD = subMatrixR
221 219
222instance Element (Complex Double) where 220instance Element (Complex Double) where
223 constantD = constantC
224 transdata = transdataC 221 transdata = transdataC
225 subMatrixD = subMatrixC 222 subMatrixD = subMatrixC
226 223
@@ -284,31 +281,6 @@ subMatrix :: Element a
284 -> Matrix a -- ^ result 281 -> Matrix a -- ^ result
285subMatrix = subMatrixD 282subMatrix = subMatrixD
286 283
287------------------------------------------------------------------------
288
289constantAux fun x n = unsafePerformIO $ do
290 v <- createVector n
291 px <- newArray [x]
292 app1 (fun px) vec v "constantAux"
293 free px
294 return v
295
296constantR :: Double -> Int -> Vector Double
297constantR = constantAux cconstantR
298foreign import ccall "auxi.h constantR" cconstantR :: Ptr Double -> TV
299
300constantC :: Complex Double -> Int -> Vector (Complex Double)
301constantC = constantAux cconstantC
302foreign import ccall "auxi.h constantC" cconstantC :: Ptr (Complex Double) -> TCV
303
304{- | creates a vector with a given number of equal components:
305
306@> constant 2 7
3077 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@
308-}
309constant :: Element a => a -> Int -> Vector a
310constant = constantD
311
312-------------------------------------------------------------------------- 284--------------------------------------------------------------------------
313 285
314-- | obtains the complex conjugate of a complex vector 286-- | obtains the complex conjugate of a complex vector
@@ -329,10 +301,6 @@ fromComplex :: Vector (Complex Double) -> (Vector Double, Vector Double)
329fromComplex z = (r,i) where 301fromComplex z = (r,i) where
330 [r,i] = toColumns $ reshape 2 $ asReal z 302 [r,i] = toColumns $ reshape 2 $ asReal z
331 303
332-- | converts a real vector into a complex representation (with zero imaginary parts)
333comp :: Vector Double -> Vector (Complex Double)
334comp v = toComplex (v,constant 0 (dim v))
335
336-- | loads a matrix efficiently from formatted ASCII text file (the number of rows and columns must be known in advance). 304-- | loads a matrix efficiently from formatted ASCII text file (the number of rows and columns must be known in advance).
337fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double) 305fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double)
338fromFile filename (r,c) = do 306fromFile filename (r,c) = do
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
index 2df33e0..f590919 100644
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -99,6 +99,7 @@ n |> l = if length l == n then fromList l else error "|> with wrong size"
99-- | access to Vector elements without range checking 99-- | access to Vector elements without range checking
100at' :: Storable a => Vector a -> Int -> a 100at' :: Storable a => Vector a -> Int -> a
101at' v n = safeRead v $ flip peekElemOff n 101at' v n = safeRead v $ flip peekElemOff n
102{-# INLINE at' #-}
102 103
103-- 104--
104-- turn off bounds checking with -funsafe at configure time. 105-- turn off bounds checking with -funsafe at configure time.
diff --git a/lib/Data/Packed/Internal/auxi.c b/lib/Data/Packed/Internal/auxi.c
index c449b9a..48b05e8 100644
--- a/lib/Data/Packed/Internal/auxi.c
+++ b/lib/Data/Packed/Internal/auxi.c
@@ -1,12 +1,5 @@
1#include "auxi.h" 1#include "auxi.h"
2#include <gsl/gsl_blas.h>
3#include <gsl/gsl_linalg.h>
4#include <gsl/gsl_matrix.h> 2#include <gsl/gsl_matrix.h>
5#include <gsl/gsl_math.h>
6#include <gsl/gsl_errno.h>
7#include <gsl/gsl_complex.h>
8#include <gsl/gsl_complex_math.h>
9#include <gsl/gsl_cblas.h>
10#include <string.h> 3#include <string.h>
11#include <stdio.h> 4#include <stdio.h>
12 5
@@ -92,27 +85,6 @@ int submatrixR(int r1, int r2, int c1, int c2, KRMAT(x),RMAT(r)) {
92} 85}
93 86
94 87
95int constantR(double * pval, RVEC(r)) {
96 DEBUGMSG("constantR")
97 int k;
98 double val = *pval;
99 for(k=0;k<rn;k++) {
100 rp[k]=val;
101 }
102 OK
103}
104
105int constantC(gsl_complex* pval, CVEC(r)) {
106 DEBUGMSG("constantC")
107 int k;
108 gsl_complex val = *pval;
109 for(k=0;k<rn;k++) {
110 rp[k]=val;
111 }
112 OK
113}
114
115
116int conjugate(KCVEC(x),CVEC(t)) { 88int conjugate(KCVEC(x),CVEC(t)) {
117 REQUIRES(xn==tn,BAD_SIZE); 89 REQUIRES(xn==tn,BAD_SIZE);
118 DEBUGMSG("conjugate"); 90 DEBUGMSG("conjugate");
diff --git a/lib/Data/Packed/Internal/auxi.h b/lib/Data/Packed/Internal/auxi.h
index b57fdaf..4698696 100644
--- a/lib/Data/Packed/Internal/auxi.h
+++ b/lib/Data/Packed/Internal/auxi.h
@@ -13,9 +13,6 @@
13int transR(KRMAT(x),RMAT(t)); 13int transR(KRMAT(x),RMAT(t));
14int transC(KCMAT(x),CMAT(t)); 14int transC(KCMAT(x),CMAT(t));
15 15
16int constantR(double *val , RVEC(r));
17int constantC(gsl_complex *val, CVEC(r));
18
19int submatrixR(int r1, int r2, int c1, int c2, KRMAT(x),RMAT(r)); 16int submatrixR(int r1, int r2, int c1, int c2, KRMAT(x),RMAT(r));
20 17
21const char * gsl_strerror (const int gsl_errno); 18const char * gsl_strerror (const int gsl_errno);
diff --git a/lib/Data/Packed/ST.hs b/lib/Data/Packed/ST.hs
index 91b7968..9c9c3b3 100644
--- a/lib/Data/Packed/ST.hs
+++ b/lib/Data/Packed/ST.hs
@@ -1,4 +1,4 @@
1{-# OPTIONS -XTypeOperators -XRank2Types -XFlexibleContexts #-} 1{-# OPTIONS -XTypeOperators -XRank2Types -XFlexibleContexts -XBangPatterns #-}
2 2
3----------------------------------------------------------------------------- 3-----------------------------------------------------------------------------
4-- | 4-- |
@@ -23,8 +23,10 @@ module Data.Packed.ST (
23 STMatrix, newMatrix, thawMatrix, freezeMatrix, runSTMatrix, 23 STMatrix, newMatrix, thawMatrix, freezeMatrix, runSTMatrix,
24 readMatrix, writeMatrix, modifyMatrix, liftSTMatrix, 24 readMatrix, writeMatrix, modifyMatrix, liftSTMatrix,
25 -- * Unsafe functions 25 -- * Unsafe functions
26 newUndefinedVector,
26 unsafeReadVector, unsafeWriteVector, 27 unsafeReadVector, unsafeWriteVector,
27 unsafeThawVector, unsafeFreezeVector, 28 unsafeThawVector, unsafeFreezeVector,
29 newUndefinedMatrix,
28 unsafeReadMatrix, unsafeWriteMatrix, 30 unsafeReadMatrix, unsafeWriteMatrix,
29 unsafeThawMatrix, unsafeFreezeMatrix 31 unsafeThawMatrix, unsafeFreezeMatrix
30) where 32) where
@@ -87,9 +89,17 @@ readVector = safeIndexV unsafeReadVector
87writeVector :: Storable t => STVector s t -> Int -> t -> ST s () 89writeVector :: Storable t => STVector s t -> Int -> t -> ST s ()
88writeVector = safeIndexV unsafeWriteVector 90writeVector = safeIndexV unsafeWriteVector
89 91
92{-# NOINLINE newUndefinedVector #-}
93newUndefinedVector :: Element t => Int -> ST s (STVector s t)
94newUndefinedVector = unsafeIOToST . fmap STVector . createVector
95
90{-# NOINLINE newVector #-} 96{-# NOINLINE newVector #-}
91newVector :: Element t => t -> Int -> ST s (STVector s t) 97newVector :: Element t => t -> Int -> ST s (STVector s t)
92newVector v = unsafeThawVector . constant v 98newVector x n = do
99 v <- newUndefinedVector n
100 let go (-1) = return v
101 go !k = unsafeWriteVector v k x >> go (k-1 :: Int)
102 go (n-1)
93 103
94------------------------------------------------------------------------- 104-------------------------------------------------------------------------
95 105
@@ -153,6 +163,10 @@ readMatrix = safeIndexM unsafeReadMatrix
153writeMatrix :: Storable t => STMatrix s t -> Int -> Int -> t -> ST s () 163writeMatrix :: Storable t => STMatrix s t -> Int -> Int -> t -> ST s ()
154writeMatrix = safeIndexM unsafeWriteMatrix 164writeMatrix = safeIndexM unsafeWriteMatrix
155 165
166{-# NOINLINE newUndefinedMatrix #-}
167newUndefinedMatrix :: Element t => MatrixOrder -> Int -> Int -> ST s (STMatrix s t)
168newUndefinedMatrix order r c = unsafeIOToST $ fmap STMatrix $ createMatrix order r c
169
156{-# NOINLINE newMatrix #-} 170{-# NOINLINE newMatrix #-}
157newMatrix :: Element t => t -> Int -> Int -> ST s (STMatrix s t) 171newMatrix :: Element t => t -> Int -> Int -> ST s (STMatrix s t)
158newMatrix v r c = unsafeThawMatrix . reshape c . constant v $ r*c 172newMatrix v r c = unsafeThawMatrix $ reshape c $ runSTVector $ newVector v (r*c)
diff --git a/lib/Data/Packed/Vector.hs b/lib/Data/Packed/Vector.hs
index 6415c5c..0bbbc34 100644
--- a/lib/Data/Packed/Vector.hs
+++ b/lib/Data/Packed/Vector.hs
@@ -24,6 +24,7 @@ module Data.Packed.Vector (
24 24
25import Data.Packed.Internal 25import Data.Packed.Internal
26import Numeric.GSL.Vector 26import Numeric.GSL.Vector
27import Data.Packed.ST
27 28
28{- | Creates a real vector containing a range of values: 29{- | Creates a real vector containing a range of values:
29 30
@@ -47,3 +48,11 @@ vectorMaxIndex = round . toScalarR MaxIdx
47 48
48vectorMinIndex :: Vector Double -> Int 49vectorMinIndex :: Vector Double -> Int
49vectorMinIndex = round . toScalarR MinIdx 50vectorMinIndex = round . toScalarR MinIdx
51
52{- | creates a vector with a given number of equal components:
53
54@> constant 2 7
557 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@
56-}
57constant :: Element a => a -> Int -> Vector a
58constant x n = runSTVector (newVector x n)