packages feed

hmatrix 0.3.0.0 → 0.4.0.0

raw patch · 16 files changed

+661/−84 lines, 16 files

Files

README view
@@ -198,8 +198,8 @@  ACKNOWLEDGEMENTS ----------------------------------------------------- -I thank Henning Thielemann and all the people in the Haskell mailing-lists for their help.+I thank Don Stewart, Henning Thielemann, Bulat Ziganshin and all the people+in the Haskell mailing lists for their help.  - Nico Mahlo discovered a bug in the eigendecomposition wrapper. @@ -221,3 +221,10 @@   which are wrong for this architecture.  - Jason Schroeder reported an error in the documentation.++- Bulat Ziganshin gave invaluable help for the ST monad interface to+  in-place modifications.++- Don Stewart fixed the implementation of the internal data structures+  to achieve excellent, C-like performance in Haskell functions which+  explicitly work with the elements of vectors and matrices.
+ examples/benchmarks.hs view
@@ -0,0 +1,103 @@+{-# OPTIONS -fbang-patterns #-}++-- compile as:+-- ghc --make -O2 -optc-O2 -fvia-C benchmarks.hs+-- ghc --make -O benchmarks.hs++import Numeric.LinearAlgebra+import System.Time+import System.CPUTime+import Text.Printf+import Data.List(foldl1')+++time act = do+    t0 <- getCPUTime+    act+    t1 <- getCPUTime+    printf "%.3f s CPU\n" $ (fromIntegral (t1 - t0) / (10^12 :: Double)) :: IO ()++--------------------------------------------------------------------------------++main = sequence_ [bench1,bench2,bench3]++w :: Vector Double+w = constant 1 30000000++bench1 = do+    putStrLn "Sum of a vector with 30M doubles:"+    print$ vectorMax w -- evaluate it+    time $ printf "     BLAS: %.2f: " $ sumVB w+    time $ printf "  Haskell: %.2f: " $ sumVH w+    time $ printf "     BLAS: %.2f: " $ sumVB w+    time $ printf "  Haskell: %.2f: " $ sumVH w++sumVB v = constant 1 (dim v) <.> v++sumVH v = go (d - 1) 0+     where+       d = dim v+       go :: Int -> Double -> Double+       go 0 s = s + (v @> 0)+       go !j !s = go (j - 1) (s + (v @> j))++--------------------------------------------------------------------------------++bench2 = do+    putStrLn "-------------------------------------------------------"+    putStrLn "Multiplication of 1M different 3x3 matrices:"+--     putStrLn "from [[]]"+--     time $ print $ fun (10^6) rot'+--     putStrLn "from []"+    time $ print $ manymult (10^6) rot+    print $ cos (10^6/2)+++rot' :: Double -> Matrix Double+rot' a = matrix [[ c,0,s],+                 [ 0,1,0],+                 [-s,0,c]]+    where c = cos a+          s = sin a+          matrix = fromLists++rot :: Double -> Matrix Double+rot a = (3><3) [ c,0,s+               , 0,1,0+               ,-s,0,c ]+    where c = cos a+          s = sin a++manymult n r = foldl1' (<>) (map r angles)+    where angles = toList $ linspace n (0,1)+          -- angles = map (k*) [0..n']+          -- n' = fromIntegral n - 1+          -- k  = recip n'++--------------------------------------------------------------------------------++bench3 = do+    putStrLn "-------------------------------------------------------"+    putStrLn "foldVector"+    let v = flatten $ ident 500 :: Vector Double+    print $ vectorMax v  -- evaluate it+    let getPos k s = if k `mod` 500 < 200 && v@>k > 0 then k:s else s+    putStrLn "indices extraction, dim=0.25M:"+    time $ print $ (`divMod` 500) $ maximum $ foldLoop getPos [] (dim v)+    putStrLn "sum, dim=30M:"+    --time $ print $ foldLoop (\k s -> w@>k + s) 0.0 (dim w)+    time $ print $ foldVector (\k v s -> v k + s) 0.0 w+    putStrLn "sum, dim=0.25M:"+    --time $ print $ foldVector (\k v s -> v k + s) 0.0 v+    time $ print $ foldLoop (\k s -> v@>k + s) 0.0 (dim v)++-- foldVector is slower if it is used in two places. (!?)+-- this does not happen with foldLoop++foldLoop f s d = go (d - 1) s+     where+       go 0 s = f (0::Int) s+       go !j !s = go (j - 1) (f j s)++foldVector f s v = foldLoop g s (dim v)+    where g !k !s = f k (v@>) s
examples/experiments/parallel.hs view
@@ -5,24 +5,32 @@  inParallel = parMap rwhnf id -parMul x y = fromBlocks [inParallel[x <> y1, x <> y2]]-    where p = cols y `div` 2-          (y1,y2) = splitColumnsAt p y+parMul p x y = fromBlocks [ inParallel ( map (x <>) ys ) ]+    where ys = splitColumns p y + main = do     n <- (read . head) `fmap` getArgs     let m = ident n :: Matrix Double-    time $ print $ vectorMax $ takeDiag $ parMul m m     time $ print $ vectorMax $ takeDiag $ m <> m--a = (2><3) [1..6::Double]-b = (3><4) [1..12::Double]--splitRowsAt p m    = (takeRows p m, dropRows p m)-splitColumnsAt p m = (takeColumns p m, dropColumns p m)+    time $ print $ vectorMax $ takeDiag $ parMul 2 m m+    time $ print $ vectorMax $ takeDiag $ parMul 4 m m+    time $ print $ vectorMax $ takeDiag $ parMul 8 m m  time act = do     t0 <- getClockTime     act     t1 <- getClockTime     print $ tdSec $ normalizeTimeDiff $ diffClockTimes t1 t0++splitColumns n m = splitColumns' (f n (cols m)) m+    where+    splitColumns' [] m = []+    splitColumns' ((a,b):rest) m = subMatrix (0,a) (rows m, b-a+1) m : splitColumns' rest m++    f :: Int -> Int -> [(Int,Int)]+    f n c = zip ks (map pred $ tail ks)+        where ks = map round $ toList $ linspace (fromIntegral n+1) (0,fromIntegral c)++splitRowsAt p m    = (takeRows p m, dropRows p m)+splitColumnsAt p m = (takeColumns p m, dropColumns p m)
+ examples/inplace.hs view
@@ -0,0 +1,152 @@+-- some tests of the interface for pure+-- computations with inplace updates++import Numeric.LinearAlgebra+import Data.Packed.ST+import Data.Packed.Convert++import Data.Array.Unboxed+import Data.Array.ST+import Control.Monad.ST+import Control.Monad++main = sequence_[+    print test1,+    print test2,+    print test3,+    print test4,+    test5,+    test6,+    print test7,+    test8,+    test0]++-- helper functions+vector l = fromList l :: Vector Double+norm v = pnorm PNorm2 v++-- hmatrix vector and matrix+v = vector [1..10]+m = (5><10) [1..50::Double]++----------------------------------------------------------------------++-- vector creation by in-place updates on a copy of the argument+test1 = fun v++fun :: Element t => Vector t -> Vector t+fun x = runSTVector $ do+    a <- thawVector x+    mapM_ (flip (modifyVector a) (+57)) [0 .. dim x `div` 2 - 1]+    return a++-- another example: creation of an antidiagonal matrix from a list+test2 = antiDiag 5 8 [1..] :: Matrix Double++antiDiag :: (Element b) => Int -> Int -> [b] -> Matrix b+antiDiag r c l = runSTMatrix $ do+    m <- newMatrix 0 r c+    let d = min r c - 1+    sequence_ $ zipWith (\i v -> writeMatrix m i (c-1-i) v) [0..d] l+    return m++-- using vector or matrix functions on mutable objects requires freezing:+test3 = g1 v++g1 x = runST $ do+    a <- thawVector x+    writeVector a (dim x -1) 0+    b <- freezeVector a+    return (norm b)++--  another possibility:+test4 = g2 v++g2 x = runST $ do+    a <- thawVector x+    writeVector a (dim x -1) 0+    t <- liftSTVector norm a+    return t++--------------------------------------------------------------++-- haskell arrays+hv = listArray (0,9) [1..10::Double]+hm = listArray ((0,0),(4,9)) [1..50::Double]++++-- conversion from standard Haskell arrays+test5 = do+    print $ norm (vectorFromArray hv)+    print $ norm v+    print $ rcond (matrixFromArray hm)+    print $ rcond m+++-- conversion to mutable ST arrays+test6 = do+    let y = clearColumn m 1+    print y+    print (matrixFromArray y)++clearColumn x c = runSTUArray $ do+    a <- mArrayFromMatrix x+    forM_ [0..rows x-1] $ \i->+        writeArray a (i,c) (0::Double)+    return a++-- hmatrix functions applied to mutable ST arrays+test7 = unitary (listArray (1,4) [3,5,7,2] :: UArray Int Double)++unitary v = runSTUArray $ do+    a <- thaw v+    n <- norm `fmap` vectorFromMArray a+    b <- mapArray (/n) a+    return b++-------------------------------------------------++-- (just to check that they are not affected)+test0 = do+    print v+    print m+    --print hv+    --print hm++-------------------------------------------------++histogram n ds = runSTVector $ do+    h <- newVector (0::Double) n -- number of bins+    let inc k = modifyVector h k (+1)+    mapM_ inc ds+    return h++-- check that newVector is really called with a fresh new array+histoCheck ds = runSTVector $ do+    h <- newVector (0::Double) 15 -- > constant for this test+    let inc k = modifyVector h k (+1)+    mapM_ inc ds+    return h++hc = fromList [1 .. 15::Double]++-- check that thawVector creates a new array+histoCheck2 ds = runSTVector $ do+    h <- thawVector hc+    let inc k = modifyVector h k (+1)+    mapM_ inc ds+    return h++test8 = do+    let ds = [0..14]+    print $ histogram 15 ds+    print $ histogram 15 ds+    print $ histogram 15 ds+    print $ histoCheck ds+    print $ histoCheck ds+    print $ histoCheck ds+    print $ histoCheck2 ds+    print $ histoCheck2 ds+    print $ histoCheck2 ds+    putStrLn "----------------------"
hmatrix.cabal view
@@ -1,19 +1,19 @@ Name:               hmatrix-Version:            0.3.0.0+Version:            0.4.0.0 License:            GPL License-file:       LICENSE Author:             Alberto Ruiz Maintainer:         Alberto Ruiz <aruiz@um.es> Stability:          provisional-Homepage:           http://alberrto.googlepages.com/gslhaskell+Homepage:           http://www.hmatrix.googlepages.com Synopsis:           Linear algebra and numerical computations Description:        A purely functional interface to basic linear algebra computations                     and other numerical routines, internally implemented using                     GSL, BLAS and LAPACK.                     .-                    More information: <http://alberrto.googlepages.com/gslhaskell>+                    More information: <http://www.hmatrix.googlepages.com> Category:           Numerical, Math-tested-with:        GHC ==6.8.2+tested-with:        GHC ==6.8.3  cabal-version:      >=1.2 build-type:         Simple@@ -25,6 +25,14 @@     description:    Link with Intel's MKL optimized libraries.     default:        False +flag gsl+    description:    Link with GSL unoptimized blas.+    default:        False++flag unsafe+    description: Compile the library with bound checking disabled.+    default: False+ library     if flag(splitBase)       build-depends:    base >= 3, array, QuickCheck, HUnit, storable-complex@@ -36,19 +44,19 @@         ghc-options:    -fvia-C      Build-Depends:      haskell98-    Extensions:         ForeignFunctionInterface+    Extensions:         ForeignFunctionInterface,+                        CPP      hs-source-dirs:     lib     Exposed-modules:    Data.Packed,                         Data.Packed.Vector,                         Data.Packed.Matrix,-                        Numeric.GSL.Vector,-                        Numeric.GSL.Matrix,                         Numeric.GSL.Differentiation,                         Numeric.GSL.Integration,                         Numeric.GSL.Fourier,                         Numeric.GSL.Polynomials,                         Numeric.GSL.Minimization,+                        Numeric.GSL.Vector,                         Numeric.GSL.Special,                         Numeric.GSL.Special.Gamma,                         Numeric.GSL.Special.Erf,@@ -85,12 +93,15 @@                         Numeric.LinearAlgebra.Algorithms,                         Graphics.Plot,                         Numeric.LinearAlgebra.Tests+                        Data.Packed.Convert+                        Data.Packed.ST     other-modules:      Data.Packed.Internal,                         Data.Packed.Internal.Common,                         Data.Packed.Internal.Vector,                         Data.Packed.Internal.Matrix,                         Numeric.GSL.Special.Internal,-                        Numeric.LinearAlgebra.Tests.Instances+                        Numeric.GSL.Matrix,+                        Numeric.LinearAlgebra.Tests.Instances,                         Numeric.LinearAlgebra.Tests.Properties     C-sources:          lib/Data/Packed/Internal/auxi.c,                         lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c,@@ -101,7 +112,14 @@       else         extra-libraries:   gsl mkl_lapack mkl_intel mkl_sequential mkl_core     else-      extra-libraries:  gsl blas lapack+      if flag(gsl)+        extra-libraries:  gsl gslcblas lapack+      else+        extra-libraries:  gsl blas lapack      cc-options:         -O4     ghc-prof-options:   -auto-all++    if flag(unsafe)+        cpp-options: -DUNSAFE+
+ lib/Data/Packed/Convert.hs view
@@ -0,0 +1,102 @@+{-# OPTIONS -XTypeOperators -XRank2Types  -XFlexibleContexts #-}++-----------------------------------------------------------------------------+-- |+-- Module      :  Data.Packed.Convert+-- Copyright   :  (c) Alberto Ruiz 2008+-- License     :  GPL-style+--+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>+-- Stability   :  provisional+-- Portability :  portable+--+-- Conversion of Vectors and Matrices to and from the standard Haskell arrays.+-- (provisional)+--+-----------------------------------------------------------------------------++module Data.Packed.Convert (+    arrayFromVector, vectorFromArray,+    mArrayFromVector, vectorFromMArray,+    vectorToStorableArray, storableArrayToVector,+    arrayFromMatrix, matrixFromArray,+    mArrayFromMatrix, matrixFromMArray,+--    matrixToStorableArray, storableArrayToMatrix+) where++import Data.Packed.Internal+import Data.Array.Storable+import Foreign+import Control.Monad.ST+import Data.Array.ST+import Data.Array.IArray+import Data.Array.Unboxed++-- | Creates a StorableArray indexed from 0 to dim -1.+--   (Memory is efficiently copied, so you can then freely modify the obtained array)+vectorToStorableArray :: Storable t => Vector t -> IO (StorableArray Int t)+vectorToStorableArray v = do+    r <- cloneVector v+    unsafeForeignPtrToStorableArray (fptr r) (0,dim r -1)++-- | Creates a Vector from a StorableArray.+--   (Memory is efficiently copied, so posterior changes in the array will not affect the result)+storableArrayToVector :: Storable t => StorableArray Int t -> IO (Vector t)+storableArrayToVector s = do+    (a,b) <- getBounds s+    let n = (b-a+1)+    r <- createVector n+    withStorableArray s $ \p -> do+        let f _ d =  copyArray d p n >> return 0+        app1 f vec r "storableArrayToVector"+    return r+++unsafeVectorToStorableArray :: Storable t => Vector t -> IO (StorableArray Int t)+unsafeVectorToStorableArray v = unsafeForeignPtrToStorableArray (fptr v) (0,dim v -1)++--unsafeStorableArrayToVector :: Storable t => StorableArray Int t -> IO (Vector t)+--unsafeStorableArrayToVector s = undefined -- the foreign ptr of Storable Array is not available?++-----------------------------------------------------------------+-- provisional, we need Unboxed arrays for Complex Double+++unsafeFreeze' :: (MArray a e m, Ix i) => a i e -> m (Array i e)+unsafeFreeze' = unsafeFreeze++-- | creates an immutable Array from an hmatrix Vector (to do: unboxed)+arrayFromVector :: (Storable t) => Vector t -> Array Int t+arrayFromVector x = runSTArray (mArrayFromVector x)++-- | creates a mutable array from an hmatrix Vector (to do: unboxed)+mArrayFromVector :: (MArray b t (ST s), Storable t) => Vector t -> ST s (b Int t)+mArrayFromVector v = unsafeThaw =<< unsafeIOToST ( unsafeFreeze' =<< (vectorToStorableArray $ v))+++-- (creates an hmatrix Vector from an immutable array (to do: unboxed))+vectorFromArray :: (Storable t) => Array Int t -> Vector t+vectorFromArray a = unsafePerformIO $ storableArrayToVector =<< unsafeThaw a++-- | creates a mutable Array from an hmatrix Vector for manipulation with runSTUArray (to do: unboxed)+vectorFromMArray :: (Storable t, MArray a t (ST s)) => a Int t -> ST s (Vector t)+vectorFromMArray x = fmap vectorFromArray (unsafeFreeze' x)++--------------------------------------------------------------------+-- provisional++matrixFromArray :: UArray (Int, Int) Double -> Matrix Double+matrixFromArray m = reshape c . fromList . elems $ m+    where ((r1,c1),(r2,c2)) = bounds m+          r = r2-r1+1+          c = c2-c1+1++arrayFromMatrix :: Matrix Double -> UArray (Int, Int) Double+arrayFromMatrix m = listArray ((0,0),(rows m -1, cols m -1)) (toList $ flatten m)+++mArrayFromMatrix :: (MArray b Double m) => Matrix Double -> m (b (Int, Int) Double)+mArrayFromMatrix = unsafeThaw . arrayFromMatrix++matrixFromMArray :: (MArray a Double (ST s)) => a (Int,Int) Double -> ST s (Matrix Double)+matrixFromMArray x = fmap matrixFromArray (unsafeFreeze x)
lib/Data/Packed/Internal.hs view
lib/Data/Packed/Internal/Matrix.hs view
@@ -1,4 +1,5 @@ {-# OPTIONS_GHC -fglasgow-exts #-}+{-# LANGUAGE CPP #-} ----------------------------------------------------------------------------- -- | -- Module      :  Data.Packed.Internal.Matrix@@ -61,9 +62,14 @@ data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq)  -- | Matrix representation suitable for GSL and LAPACK computations.-data Matrix t = MC { rows :: Int, cols :: Int, cdat :: Vector t }-              | MF { rows :: Int, cols :: Int, fdat :: Vector t }+data Matrix t = MC { rows :: {-# UNPACK #-} !Int+                   , cols :: {-# UNPACK #-} !Int+                   , cdat :: {-# UNPACK #-} !(Vector t) } +              | MF { rows :: {-# UNPACK #-} !Int+                   , cols :: {-# UNPACK #-} !Int+                   , fdat :: {-# UNPACK #-} !(Vector t) }+ -- MC: preferred by C, fdat may require a transposition -- MF: preferred by LAPACK, cdat may require a transposition @@ -100,7 +106,6 @@ flatten :: Element t => Matrix t -> Vector t flatten = cdat . cmat - type Mt t s = Int -> Int -> Ptr t -> s -- not yet admitted by my haddock version -- infixr 6 ::>@@ -133,7 +138,6 @@ toColumns :: Element t => Matrix t -> [Vector t] toColumns m = toRows . trans $ m - -- | Reads a matrix position. (@@>) :: Storable t => Matrix t -> (Int,Int) -> t infixl 9 @@>@@ -142,12 +146,17 @@ --    | otherwise   = cdat m `at` (i*c+j)  MC {rows = r, cols = c, cdat = v} @@> (i,j)-    | i<0 || i>=r || j<0 || j>=c = error "matrix indexing out of range"-    | otherwise                  = v `at` (i*c+j)+    | safe      = if i<0 || i>=r || j<0 || j>=c+                    then error "matrix indexing out of range"+                    else v `at` (i*c+j)+    | otherwise = v `at` (i*c+j)  MF {rows = r, cols = c, fdat = v} @@> (i,j)-    | i<0 || i>=r || j<0 || j>=c = error "matrix indexing out of range"-    | otherwise                  = v `at` (j*r+i)+    | safe      = if i<0 || i>=r || j<0 || j>=c+                    then error "matrix indexing out of range"+                    else v `at` (j*r+i)+    | otherwise = v `at` (j*r+i)+{-# INLINE (@@>) #-}  ------------------------------------------------------------------ @@ -254,8 +263,8 @@         r2 = dim d `div` c2         noneed = r1 == 1 || c1 == 1 -foreign import ccall unsafe "auxi.h transR" ctransR :: TMM-foreign import ccall unsafe "auxi.h transC" ctransC :: TCMCM+foreign import ccall "auxi.h transR" ctransR :: TMM+foreign import ccall "auxi.h transC" ctransC :: TCMCM  ------------------------------------------------------------------ @@ -275,10 +284,10 @@     return r  multiplyR = multiplyAux cmultiplyR-foreign import ccall unsafe "auxi.h multiplyR" cmultiplyR :: TauxMul Double+foreign import ccall "auxi.h multiplyR" cmultiplyR :: TauxMul Double  multiplyC = multiplyAux cmultiplyC-foreign import ccall unsafe "auxi.h multiplyC" cmultiplyC :: TauxMul (Complex Double)+foreign import ccall "auxi.h multiplyC" cmultiplyC :: TauxMul (Complex Double)  -- | matrix product multiply :: (Element a) => Matrix a -> Matrix a -> Matrix a
lib/Data/Packed/Internal/Vector.hs view
@@ -1,4 +1,4 @@-{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}+{-# LANGUAGE MagicHash, CPP, UnboxedTuples #-} ----------------------------------------------------------------------------- -- | -- Module      :  Data.Packed.Internal.Vector@@ -21,10 +21,20 @@ import Complex import Control.Monad(when) +#if __GLASGOW_HASKELL__ >= 605+import GHC.ForeignPtr           (mallocPlainForeignPtrBytes)+#else+import Foreign.ForeignPtr       (mallocForeignPtrBytes)+#endif++import GHC.Base+import GHC.IOBase+ -- | A one-dimensional array of objects stored in a contiguous memory block.-data Vector t = V { dim  :: Int              -- ^ number of elements-                  , fptr :: ForeignPtr t     -- ^ foreign pointer to the memory block-                  }+data Vector t =+    V { dim  :: {-# UNPACK #-} !Int               -- ^ number of elements+      , fptr :: {-# UNPACK #-} !(ForeignPtr t)    -- ^ foreign pointer to the memory block+      }  vec = withVector @@ -37,8 +47,20 @@ createVector :: Storable a => Int -> IO (Vector a) createVector n = do     when (n <= 0) $ error ("trying to createVector of dim "++show n)-    fp <- mallocForeignPtrArray n+    fp <- doMalloc undefined     return $ V n fp+  where+    --+    -- Use the much cheaper Haskell heap allocated storage+    -- for foreign pointer space we control+    --+    doMalloc :: Storable b => b -> IO (ForeignPtr b)+    doMalloc dummy = do+#if __GLASGOW_HASKELL__ >= 605+        mallocPlainForeignPtrBytes (n * sizeOf dummy)+#else+        mallocForeignPtrBytes      (n * sizeOf dummy)+#endif  {- | creates a Vector from a list: @@ -53,8 +75,13 @@     app1 f vec v "fromList"     return v -safeRead v = unsafePerformIO . withForeignPtr (fptr v)+safeRead v = inlinePerformIO . withForeignPtr (fptr v)+{-# INLINE safeRead #-} +inlinePerformIO :: IO a -> a+inlinePerformIO (IO m) = case m realWorld# of (# _, r #) -> r+{-# INLINE inlinePerformIO #-}+ {- | extracts the Vector elements to a list  @> toList (linspace 5 (1,10))@@ -73,10 +100,25 @@ at' :: Storable a => Vector a -> Int -> a at' v n = safeRead v $ flip peekElemOff n +--+-- turn off bounds checking with -funsafe at configure time.+-- ghc will optimise away the salways true case at compile time.+--+#if defined(UNSAFE)+safe :: Bool+safe = False+#else+safe = True+#endif+ -- | access to Vector elements with range checking. at :: Storable a => Vector a -> Int -> a-at v n | n >= 0 && n < dim v = at' v n-       | otherwise          = error "vector index out of range"+at v n+    | safe      = if n >= 0 && n < dim v+                    then at' v n+                    else error "vector index out of range"+    | otherwise = at' v n+{-# INLINE at #-}  {- | takes a number of consecutive elements from a Vector @@ -146,3 +188,10 @@ liftVector2 f u v = fromList $ zipWith f (toList u) (toList v)  -----------------------------------------------------------------++cloneVector :: Storable t => Vector t -> IO (Vector t)+cloneVector (v@V {dim=n}) = do+        r <- createVector n+        let f _ s _ d =  copyArray d s n >> return 0+        app2 f vec v vec r "cloneVector"+        return r
lib/Data/Packed/Matrix.hs view
+ lib/Data/Packed/ST.hs view
@@ -0,0 +1,160 @@+{-# OPTIONS -XTypeOperators -XRank2Types  -XFlexibleContexts #-}++-----------------------------------------------------------------------------+-- |+-- Module      :  Data.Packed.ST+-- Copyright   :  (c) Alberto Ruiz 2008+-- License     :  GPL-style+--+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>+-- Stability   :  provisional+-- Portability :  portable+--+-- In-place manipulation inside the ST monad.+-- See examples/inplace.hs in the distribution.+--+-----------------------------------------------------------------------------++module Data.Packed.ST (+    -- * Mutable Vectors+    STVector, newVector, thawVector, freezeVector, runSTVector,+    readVector, writeVector, modifyVector, liftSTVector,+    -- * Mutable Matrices+    STMatrix, newMatrix, thawMatrix, freezeMatrix, runSTMatrix,+    readMatrix, writeMatrix, modifyMatrix, liftSTMatrix,+    -- * Unsafe functions+    unsafeReadVector, unsafeWriteVector,+    unsafeThawVector, unsafeFreezeVector,+    unsafeReadMatrix, unsafeWriteMatrix,+    unsafeThawMatrix, unsafeFreezeMatrix+) where++import Data.Packed.Internal+import Data.Array.Storable+import Control.Monad.ST+import Data.Array.ST+import Foreign++{-# INLINE ioReadV #-}+ioReadV :: Storable t => Vector t -> Int -> IO t+ioReadV v k = withForeignPtr (fptr v) $ \s -> peekElemOff s k++{-# INLINE ioWriteV #-}+ioWriteV :: Storable t => Vector t -> Int -> t -> IO ()+ioWriteV v k x = withForeignPtr (fptr v) $ \s -> pokeElemOff s k x++newtype STVector s t = STVector (Vector t)++thawVector :: Storable t => Vector t -> ST s (STVector s t)+thawVector = unsafeIOToST . fmap STVector . cloneVector++unsafeThawVector :: Storable t => Vector t -> ST s (STVector s t)+unsafeThawVector = unsafeIOToST . return . STVector++runSTVector :: Storable t => (forall s . ST s (STVector s t)) -> Vector t+runSTVector st = runST (st >>= unsafeFreezeVector)++{-# INLINE unsafeReadVector #-}+unsafeReadVector :: Storable t => STVector s t -> Int -> ST s t+unsafeReadVector   (STVector x) = unsafeIOToST . ioReadV x++{-# INLINE unsafeWriteVector #-}+unsafeWriteVector :: Storable t => STVector s t -> Int -> t -> ST s ()+unsafeWriteVector  (STVector x) k = unsafeIOToST . ioWriteV x k++{-# INLINE modifyVector #-}+modifyVector :: (Storable t) => STVector s t -> Int -> (t -> t) -> ST s ()+modifyVector x k f = readVector x k >>= return . f >>= unsafeWriteVector x k++liftSTVector :: (Storable t) => (Vector t -> a) -> STVector s1 t -> ST s2 a+liftSTVector f (STVector x) = unsafeIOToST . fmap f . cloneVector $ x++freezeVector :: (Storable t) => STVector s1 t -> ST s2 (Vector t)+freezeVector v = liftSTVector id v++unsafeFreezeVector :: (Storable t) => STVector s1 t -> ST s2 (Vector t)+unsafeFreezeVector (STVector x) = unsafeIOToST . return $ x++{-# INLINE safeIndexV #-}+safeIndexV f (STVector v) k+    | k < 0 || k>= dim v = error $ "out of range error in vector (dim="+                                   ++show (dim v)++", pos="++show k++")"+    | otherwise = f (STVector v) k++{-# INLINE readVector #-}+readVector :: Storable t => STVector s t -> Int -> ST s t+readVector = safeIndexV unsafeReadVector++{-# INLINE writeVector #-}+writeVector :: Storable t => STVector s t -> Int -> t -> ST s ()+writeVector = safeIndexV unsafeWriteVector++{-# NOINLINE newVector #-}+newVector :: Element t => t -> Int -> ST s (STVector s t)+newVector v = unsafeThawVector . constant v++-------------------------------------------------------------------------++{-# INLINE ioReadM #-}+ioReadM :: Storable t => Matrix t -> Int -> Int -> IO t+ioReadM (MC nr nc cv) r c = ioReadV cv (r*nc+c)+ioReadM (MF nr nc fv) r c = ioReadV fv (c*nr+r)++{-# INLINE ioWriteM #-}+ioWriteM :: Storable t => Matrix t -> Int -> Int -> t -> IO ()+ioWriteM (MC nr nc cv) r c val = ioWriteV cv (r*nc+c) val+ioWriteM (MF nr nc fv) r c val = ioWriteV fv (c*nr+r) val++newtype STMatrix s t = STMatrix (Matrix t)++thawMatrix :: Storable t => Matrix t -> ST s (STMatrix s t)+thawMatrix = unsafeIOToST . fmap STMatrix . cloneMatrix++unsafeThawMatrix :: Storable t => Matrix t -> ST s (STMatrix s t)+unsafeThawMatrix = unsafeIOToST . return . STMatrix++runSTMatrix :: Storable t => (forall s . ST s (STMatrix s t)) -> Matrix t+runSTMatrix st = runST (st >>= unsafeFreezeMatrix)++{-# INLINE unsafeReadMatrix #-}+unsafeReadMatrix :: Storable t => STMatrix s t -> Int -> Int -> ST s t+unsafeReadMatrix   (STMatrix x) r = unsafeIOToST . ioReadM x r++{-# INLINE unsafeWriteMatrix #-}+unsafeWriteMatrix :: Storable t => STMatrix s t -> Int -> Int -> t -> ST s ()+unsafeWriteMatrix  (STMatrix x) r c = unsafeIOToST . ioWriteM x r c++{-# INLINE modifyMatrix #-}+modifyMatrix :: (Storable t) => STMatrix s t -> Int -> Int -> (t -> t) -> ST s ()+modifyMatrix x r c f = readMatrix x r c >>= return . f >>= unsafeWriteMatrix x r c++liftSTMatrix :: (Storable t) => (Matrix t -> a) -> STMatrix s1 t -> ST s2 a+liftSTMatrix f (STMatrix x) = unsafeIOToST . fmap f . cloneMatrix $ x++unsafeFreezeMatrix :: (Storable t) => STMatrix s1 t -> ST s2 (Matrix t)+unsafeFreezeMatrix (STMatrix x) = unsafeIOToST . return $ x++freezeMatrix :: (Storable t) => STMatrix s1 t -> ST s2 (Matrix t)+freezeMatrix m = liftSTMatrix id m++cloneMatrix (MC r c d) = cloneVector d >>= return . MC r c+cloneMatrix (MF r c d) = cloneVector d >>= return . MF r c++{-# INLINE safeIndexM #-}+safeIndexM f (STMatrix m) r c+    | r<0 || r>=rows m ||+      c<0 || c>=cols m = error $ "out of range error in matrix (size="+                                 ++show (rows m,cols m)++", pos="++show (r,c)++")"+    | otherwise = f (STMatrix m) r c++{-# INLINE readMatrix #-}+readMatrix :: Storable t => STMatrix s t -> Int -> Int -> ST s t+readMatrix = safeIndexM unsafeReadMatrix++{-# INLINE writeMatrix #-}+writeMatrix :: Storable t => STMatrix s t -> Int -> Int -> t -> ST s ()+writeMatrix = safeIndexM unsafeWriteMatrix++{-# NOINLINE newMatrix #-}+newMatrix :: Element t => t -> Int -> Int -> ST s (STMatrix s t)+newMatrix v r c = unsafeThawMatrix . reshape c . constant v $ r*c
lib/Numeric/GSL/Fourier.hs view
@@ -35,7 +35,7 @@  {- | Fast 1D Fourier transform of a 'Vector' @(@'Complex' 'Double'@)@ using /gsl_fft_complex_forward/. It uses the same scaling conventions as GNU Octave. -@> fft ('GSL.Matrix.fromList' [1,2,3,4])+@> fft ('fromList' [1,2,3,4]) vector (4) [10.0 :+ 0.0,(-2.0) :+ 2.0,(-2.0) :+ 0.0,(-2.0) :+ (-2.0)]@  -}
lib/Numeric/GSL/Matrix.hs view
@@ -8,7 +8,7 @@ -- Stability   :  provisional -- Portability :  portable (uses FFI) ----- A few linear algebra computations based on the Numeric.GSL (<http://www.gnu.org/software/Numeric.GSL>).+-- A few linear algebra computations based on GSL. -- ----------------------------------------------------------------------------- -- #hide@@ -92,23 +92,6 @@  {- | Singular value decomposition of a real matrix, using /gsl_linalg_SV_decomp_mod/: -@\> let (u,s,v) = svdg $ 'fromLists' [[1,2,3],[-4,1,7]]-\-\> u-0.310 -0.951-0.951  0.310-\-\> s-8.497 2.792-\-\> v--0.411 -0.785- 0.185 -0.570- 0.893 -0.243-\-\> u \<\> 'diag' s \<\> 'trans' v- 1. 2. 3.--4. 1. 7.@  -} svdg :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)@@ -128,20 +111,6 @@  {- | QR decomposition of a real matrix using /gsl_linalg_QR_decomp/ and /gsl_linalg_QR_unpack/. -@\> let (q,r) = qr $ 'fromLists' [[1,3,5,7],[2,0,-2,4]]-\-\> q--0.447 -0.894--0.894  0.447-\-\> r--2.236 -1.342 -0.447 -6.708-    0. -2.683 -5.367 -4.472-\-\> q \<\> r-1.000 3.000  5.000 7.000-2.000    0. -2.000 4.000@- -} qr :: Matrix Double -> (Matrix Double, Matrix Double) qr = qr' . cmat@@ -299,7 +268,7 @@   2.+3.i   -7.   0.       1.    2.  -3.       1.  -1.i  2.i-\ +\ -- CPP \> l \<\> u  2.+3.i   -7.   0.      1.    2.  -3.
lib/Numeric/GSL/Minimization.hs view
@@ -31,14 +31,14 @@ {- | 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:  @minimize f xi = minimizeNMSimplex f xi (replicate (length xi) 1) 1e-2 100-\ +\  --  f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30-\ +\  -- main = do     let (s,p) = minimize f [5,7]     print s     print p-\ +\  -- \> main [0.9920430849306285,1.9969168063253164] 0. 512.500    1.082 6.500    5.@@ -104,14 +104,14 @@  @minimize = minimizeConjugateGradient 1E-2 1E-4 1E-3 30 f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30-\ +\  -- df [x,y] = [20*(x-1), 40*(y-2)]-\  +\   -- main = do     let (s,p) = minimize f df [5,7]     print s     print p-\ +\  -- \> main [1.0,2.0]  0. 687.848 4.996 6.991
lib/Numeric/LinearAlgebra/Algorithms.hs view
@@ -300,12 +300,12 @@ @\> let m = 'fromLists' [[1,0,    0]                     ,[0,1,    0]                     ,[0,0,1e-10]]-\ +\  -- \> 'pinv' m  1. 0.           0. 0. 1.           0. 0. 0. 10000000000.-\ +\  -- \> pinvTol 1E8 m 1. 0. 0. 0. 1. 0.
lib/Numeric/LinearAlgebra/Instances.hs view
@@ -9,7 +9,7 @@ Stability   :  provisional Portability :  portable -This module exports Show, Eq, Num, Fractional, and Floating instances for Vector and Matrix.+This module exports Show, Read, Eq, Num, Fractional, and Floating instances for Vector and Matrix.  In the context of the standard numeric operators, one-component vectors and matrices automatically expand to match the dimensions of the other operand.