From 2044be5a83d77536daed0f0f34d2baa6aa548dd1 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Wed, 5 Nov 2008 14:01:11 +0000 Subject: improved constant --- lib/Data/Packed.hs | 11 ++++++++--- lib/Data/Packed/Internal/Matrix.hs | 32 -------------------------------- lib/Data/Packed/Internal/Vector.hs | 1 + lib/Data/Packed/Internal/auxi.c | 28 ---------------------------- lib/Data/Packed/Internal/auxi.h | 3 --- lib/Data/Packed/ST.hs | 20 +++++++++++++++++--- lib/Data/Packed/Vector.hs | 9 +++++++++ 7 files changed, 35 insertions(+), 69 deletions(-) (limited to 'lib/Data') 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 ( import Data.Packed.Vector import Data.Packed.Matrix import Data.Complex -import Data.Packed.Internal(fromComplex,toComplex,comp,conj) +import Data.Packed.Internal(fromComplex,toComplex,conj) -- | conversion utilities class (Element e) => Container c e where @@ -38,7 +38,7 @@ class (Element e) => Container c e where instance Container Vector Double where toComplex = Data.Packed.Internal.toComplex fromComplex = Data.Packed.Internal.fromComplex - comp = Data.Packed.Internal.comp + comp = internalcomp conj = Data.Packed.Internal.conj real = id complex = Data.Packed.comp @@ -56,7 +56,7 @@ instance Container Matrix Double where fromComplex z = (reshape c r, reshape c i) where (r,i) = Data.Packed.fromComplex (flatten z) c = cols z - comp = liftMatrix Data.Packed.Internal.comp + comp = liftMatrix internalcomp conj = liftMatrix Data.Packed.Internal.conj real = id complex = Data.Packed.comp @@ -68,3 +68,8 @@ instance Container Matrix (Complex Double) where conj = undefined real = Data.Packed.comp complex = id + + +-- | converts a real vector into a complex representation (with zero imaginary parts) +internalcomp :: Vector Double -> Vector (Complex Double) +internalcomp 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 -- | Optimized matrix computations are provided for elements in the Element class. class (Storable a, Floating a) => Element a where - constantD :: a -> Int -> Vector a transdata :: Int -> Vector a -> Int -> Vector a subMatrixD :: (Int,Int) -- ^ (r0,c0) starting position -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix -> Matrix a -> Matrix a instance Element Double where - constantD = constantR transdata = transdataR subMatrixD = subMatrixR instance Element (Complex Double) where - constantD = constantC transdata = transdataC subMatrixD = subMatrixC @@ -284,31 +281,6 @@ subMatrix :: Element a -> Matrix a -- ^ result subMatrix = subMatrixD ------------------------------------------------------------------------- - -constantAux fun x n = unsafePerformIO $ do - v <- createVector n - px <- newArray [x] - app1 (fun px) vec v "constantAux" - free px - return v - -constantR :: Double -> Int -> Vector Double -constantR = constantAux cconstantR -foreign import ccall "auxi.h constantR" cconstantR :: Ptr Double -> TV - -constantC :: Complex Double -> Int -> Vector (Complex Double) -constantC = constantAux cconstantC -foreign import ccall "auxi.h constantC" cconstantC :: Ptr (Complex Double) -> TCV - -{- | creates a vector with a given number of equal components: - -@> constant 2 7 -7 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@ --} -constant :: Element a => a -> Int -> Vector a -constant = constantD - -------------------------------------------------------------------------- -- | obtains the complex conjugate of a complex vector @@ -329,10 +301,6 @@ fromComplex :: Vector (Complex Double) -> (Vector Double, Vector Double) fromComplex z = (r,i) where [r,i] = toColumns $ reshape 2 $ asReal z --- | converts a real vector into a complex representation (with zero imaginary parts) -comp :: Vector Double -> Vector (Complex Double) -comp v = toComplex (v,constant 0 (dim v)) - -- | loads a matrix efficiently from formatted ASCII text file (the number of rows and columns must be known in advance). fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double) fromFile 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" -- | access to Vector elements without range checking at' :: Storable a => Vector a -> Int -> a at' v n = safeRead v $ flip peekElemOff n +{-# INLINE at' #-} -- -- 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 @@ #include "auxi.h" -#include -#include #include -#include -#include -#include -#include -#include #include #include @@ -92,27 +85,6 @@ int submatrixR(int r1, int r2, int c1, int c2, KRMAT(x),RMAT(r)) { } -int constantR(double * pval, RVEC(r)) { - DEBUGMSG("constantR") - int k; - double val = *pval; - for(k=0;k STVector s t -> Int -> t -> ST s () writeVector = safeIndexV unsafeWriteVector +{-# NOINLINE newUndefinedVector #-} +newUndefinedVector :: Element t => Int -> ST s (STVector s t) +newUndefinedVector = unsafeIOToST . fmap STVector . createVector + {-# NOINLINE newVector #-} newVector :: Element t => t -> Int -> ST s (STVector s t) -newVector v = unsafeThawVector . constant v +newVector x n = do + v <- newUndefinedVector n + let go (-1) = return v + go !k = unsafeWriteVector v k x >> go (k-1 :: Int) + go (n-1) ------------------------------------------------------------------------- @@ -153,6 +163,10 @@ readMatrix = safeIndexM unsafeReadMatrix writeMatrix :: Storable t => STMatrix s t -> Int -> Int -> t -> ST s () writeMatrix = safeIndexM unsafeWriteMatrix +{-# NOINLINE newUndefinedMatrix #-} +newUndefinedMatrix :: Element t => MatrixOrder -> Int -> Int -> ST s (STMatrix s t) +newUndefinedMatrix order r c = unsafeIOToST $ fmap STMatrix $ createMatrix order r c + {-# NOINLINE newMatrix #-} newMatrix :: Element t => t -> Int -> Int -> ST s (STMatrix s t) -newMatrix v r c = unsafeThawMatrix . reshape c . constant v $ r*c +newMatrix 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 ( import Data.Packed.Internal import Numeric.GSL.Vector +import Data.Packed.ST {- | Creates a real vector containing a range of values: @@ -47,3 +48,11 @@ vectorMaxIndex = round . toScalarR MaxIdx vectorMinIndex :: Vector Double -> Int vectorMinIndex = round . toScalarR MinIdx + +{- | creates a vector with a given number of equal components: + +@> constant 2 7 +7 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@ +-} +constant :: Element a => a -> Int -> Vector a +constant x n = runSTVector (newVector x n) -- cgit v1.2.3