hmatrix 0.9.3.0 → 0.10.0.0

raw patch · 43 files changed

+3256/−1152 lines, 43 filesdep +binarydep +randomdep ~basedep ~vectorsetup-changedPVP ok

version bump matches the API change (PVP)

Dependencies added: binary, random

Dependency ranges changed: base, vector

API changes (from Hackage documentation)

Files

• CHANGES +18/−0
• Setup.lhs +1/−1
• THANKS +7/−0
• configure.hs +29/−20
• examples/Real.hs +2/−2
• examples/monadic.hs +118/−0
• examples/parallel.hs +5/−5
• examples/pca1.hs +1/−1
• examples/pca2.hs +1/−1
• examples/vector.hs +2/−3
• hmatrix.cabal +30/−18
• lib/Data/Packed.hs +8/−7
• lib/Data/Packed/Development.hs +0/−1
• lib/Data/Packed/Internal/Common.hs +20/−34
• lib/Data/Packed/Internal/Matrix.hs +66/−23
• lib/Data/Packed/Internal/Signatures.hs +18/−0
• lib/Data/Packed/Internal/Vector.hs +133/−15
• lib/Data/Packed/Matrix.hs +90/−210
• lib/Data/Packed/Random.hs +11/−11
• lib/Data/Packed/ST.hs +4/−4
• lib/Data/Packed/Vector.hs +43/−51
• lib/Graphics/Plot.hs +28/−67
• lib/Numeric/Chain.hs +140/−0
• lib/Numeric/Container.hs +132/−0
• lib/Numeric/ContainerBoot.hs +583/−0
• lib/Numeric/Conversion.hs +91/−0
• lib/Numeric/GSL/Fitting.hs +1/−1
• lib/Numeric/GSL/Vector.hs +162/−4
• lib/Numeric/GSL/gsl-aux.c +417/−0
• lib/Numeric/IO.hs +160/−0
• lib/Numeric/LinearAlgebra.hs +11/−6
• lib/Numeric/LinearAlgebra/Algorithms.hs +87/−142
• lib/Numeric/LinearAlgebra/Instances.hs +0/−218
• lib/Numeric/LinearAlgebra/Interface.hs +0/−117
• lib/Numeric/LinearAlgebra/LAPACK.hs +20/−9
• lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c +227/−64
• lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h +30/−4
• lib/Numeric/LinearAlgebra/Linear.hs +0/−76
• lib/Numeric/LinearAlgebra/Tests.hs +215/−10
• lib/Numeric/LinearAlgebra/Tests/Instances.hs +25/−11
• lib/Numeric/LinearAlgebra/Tests/Properties.hs +25/−16
• lib/Numeric/Matrix.hs +72/−0
• lib/Numeric/Vector.hs +223/−0

@@ -1,3 +1,21 @@+0.10.0.0+========++- Module reorganization++- Support for Float and Complex Float elements (excluding LAPACK computations)++- Binary instances for Vector and Matrix++- optimiseMult++- mapVectorM, mapVectorWithIndexM, unzipVectorWith, and related functions.++- diagRect admits diagonal vectors of any length without producing an error,+  and takes an additional argument for the off-diagonal elements.++- different signatures in some functions+ 0.9.3.0 ======= 

@@ -1,7 +1,7 @@ #! /usr/bin/env runhaskell  > import Distribution.Simple-> import System(system)+> import System.Process(system)  > main = defaultMainWithHooks autoconfUserHooks {runTests = t} 

@@ -1,6 +1,11 @@ I thank Don Stewart, Henning Thielemann, Bulat Ziganshin, Heinrich Apfelmus, and all the people in the Haskell mailing lists for their help. +I am particularly grateful to Vivian McPhail for his excellent+contributions: improved configure.hs, Binary instances for+Vector and Matrix, support for Float and Complex Float elements,+module reorganization, monadic mapVectorM, and many other improvements.+ - Nico Mahlo discovered a bug in the eigendecomposition wrapper.  - Frederik Eaton discovered a bug in the design of the wrappers.@@ -71,4 +76,6 @@ - Tim Sears reported the zgesdd problem also in intel mac.  - Max Suica simplified the installation on Windows and improved the instructions.++- John Billings reported an incompatibility with QuickCheck>=2.1.1 

@@ -17,7 +17,10 @@  -} -import System+import System.Process+import System.Exit+import System.Environment+import System.Directory(createDirectoryIfMissing) import Data.List(isPrefixOf, intercalate) import Distribution.Simple.LocalBuildInfo import Distribution.Simple.Configure@@ -37,13 +40,14 @@        ]  -- compile a simple program with symbols from GSL and LAPACK with the given libs-testprog buildInfo libs fmks =+testprog bInfo buildInfo libs fmks =     "echo \"#include <gsl/gsl_sf_gamma.h>\nint main(){zgesvd_(); gsl_sf_gamma(5);}\""-                     ++" > /tmp/dummy.c; gcc "+                     ++" > " ++ (buildDir bInfo) ++ "/dummy.c; gcc "                      ++ (join $ ccOptions buildInfo) ++ " "                      ++ (join $ cppOptions buildInfo) ++ " "-                     ++ (join $ map ("-I"++) $ includeDirs buildInfo)-                     ++" /tmp/dummy.c -o /tmp/dummy "+                     ++ (join $ map ("-I"++) $ includeDirs buildInfo) ++ " " +                     ++ (buildDir bInfo) ++ "/dummy.c -o "+                     ++ (buildDir bInfo) ++ "/dummy "                      ++ (join $ map ("-L"++) $ extraLibDirs buildInfo) ++ " "                      ++ (prepend "-l" $ libs) ++ " "                      ++ (prepend "-framework " fmks) ++ " > /dev/null 2> /dev/null"@@ -51,26 +55,28 @@ join = intercalate " " prepend x = unwords . map (x++) . words -check buildInfo libs fmks = (ExitSuccess ==) `fmap` system (testprog buildInfo libs fmks)+check bInfo buildInfo libs fmks = (ExitSuccess ==) `fmap` system (testprog bInfo buildInfo libs fmks)  -- simple test for GSL-gsl buildInfo = "echo \"#include <gsl/gsl_sf_gamma.h>\nint main(){gsl_sf_gamma(5);}\""-           ++" > /tmp/dummy.c; gcc "+gsl bInfo buildInfo = "echo \"#include <gsl/gsl_sf_gamma.h>\nint main(){gsl_sf_gamma(5);}\""+           ++" > " ++ (buildDir bInfo) ++ "/dummy.c; gcc "            ++ (join $ ccOptions buildInfo) ++ " "            ++ (join $ cppOptions buildInfo) ++ " "-           ++ (join $ map ("-I"++) $ includeDirs buildInfo)-           ++ " /tmp/dummy.c -o /tmp/dummy "+           ++ (join $ map ("-I"++) $ includeDirs buildInfo) ++ " " +           ++ (buildDir bInfo) ++ "/dummy.c -o "+           ++ (buildDir bInfo) ++ "/dummy "            ++ (join $ map ("-L"++) $ extraLibDirs buildInfo) ++ " -lgsl -lgslcblas"            ++ " > /dev/null 2> /dev/null"  -- test for gsl >= 1.12-gsl112 buildInfo =+gsl112 bInfo buildInfo =     "echo \"#include <gsl/gsl_sf_exp.h>\nint main(){gsl_sf_exprel_n_CF_e(1,1,0);}\""-           ++" > /tmp/dummy.c; gcc /tmp/dummy.c "+           ++" > " ++ (buildDir bInfo) ++ "/dummy.c; gcc " +           ++ (buildDir bInfo) ++ "/dummy.c "            ++ (join $ ccOptions buildInfo) ++ " "            ++ (join $ cppOptions buildInfo) ++ " "            ++ (join $ map ("-I"++) $ includeDirs buildInfo)-           ++" -o /tmp/dummy "+           ++" -o " ++ (buildDir bInfo) ++ "/dummy "            ++ (join $ map ("-L"++) $ extraLibDirs buildInfo) ++ " -lgsl -lgslcblas"            ++ " > /dev/null 2> /dev/null" @@ -78,11 +84,11 @@ checkCommand c = (ExitSuccess ==) `fmap` system c  -- test different configurations until the first one works-try _ _ _ [] = return Nothing-try i b f (opt:rest) = do-    ok <- check i (b ++ " " ++ opt) f+try _ _ _ _ [] = return Nothing+try l i b f (opt:rest) = do+    ok <- check l i (b ++ " " ++ opt) f     if ok then return (Just opt)-          else try i b f rest+          else try l i b f rest  -- read --configure-option=link:lib1,lib2,lib3,etc linkop = "link:"@@ -110,11 +116,14 @@     let pref = if null (words (base ++ " " ++ auxpref)) then "gsl lapack" else auxpref         fullOpts = map ((pref++" ")++) opts -    r <- try buildInfo base fwks fullOpts+    -- create the build directory (used for tmp files) if necessary+    createDirectoryIfMissing True $ buildDir bInfo+    +    r <- try bInfo buildInfo base fwks fullOpts     case r of         Nothing -> do             putStrLn " FAIL"-            g  <- checkCommand $ gsl buildInfo+            g  <- checkCommand $ gsl bInfo buildInfo             if g                 then putStrLn " *** Sorry, I can't link LAPACK."                 else putStrLn " *** Sorry, I can't link GSL."@@ -124,7 +133,7 @@             writeFile "hmatrix.buildinfo" ("buildable: False\n")         Just ops -> do             putStrLn " OK"-            g <- checkCommand $ gsl112 buildInfo+            g <- checkCommand $ gsl112 bInfo buildInfo             writeFile "hmatrix.buildinfo" $ "extra-libraries: " ++                 ops ++ "\n" ++                 if g

@@ -67,13 +67,13 @@ zeros :: Int -- ^ rows       -> Int -- ^ columns       -> Matrix Double-zeros r c = reshape c (constant 0 (r*c))+zeros r c = konst 0 (r,c)  -- | Create a matrix or ones. ones :: Int -- ^ rows      -> Int -- ^ columns      -> Matrix Double-ones r c = reshape c (constant 1 (r*c))+ones r c = konst 1 (r,c)  -- | Concatenation of real vectors. infixl 9 #

@@ -0,0 +1,118 @@+-- monadic computations+-- (contributed by Vivian McPhail)++import Numeric.LinearAlgebra+import Control.Monad.State.Strict+import Control.Monad.Maybe+import Foreign.Storable(Storable)+import System.Random(randomIO)++-------------------------------------------++-- an instance of MonadIO, a monad transformer+type VectorMonadT = StateT Int IO++test1 :: Vector Int -> IO (Vector Int)+test1 = mapVectorM $ \x -> do+    putStr $ (show x) ++ " "+    return (x + 1)++-- we can have an arbitrary monad AND do IO+addInitialM :: Vector Int -> VectorMonadT ()+addInitialM = mapVectorM_ $ \x -> do+    i <- get+    liftIO $ putStr $ (show $ x + i) ++ " "+    put $ x + i++-- sum the values of the even indiced elements+sumEvens :: Vector Int -> Int+sumEvens = foldVectorWithIndex (\x a b -> if x `mod` 2 == 0 then a + b else b) 0++-- sum and print running total of evens+sumEvensAndPrint :: Vector Int -> VectorMonadT ()+sumEvensAndPrint = mapVectorWithIndexM_ $ \ i x -> do+    when (i `mod` 2 == 0) $ do+        v <- get+        put $ v + x+        v' <- get+        liftIO $ putStr $ (show v') ++ " "+++indexPlusSum :: Vector Int -> VectorMonadT ()+indexPlusSum v' = do+    let f i x = do+            s <- get+            let inc = x+s+            liftIO $ putStr $ show (i,inc) ++ " "+            put inc+            return inc+    v <- mapVectorWithIndexM f v'+    liftIO $ do+        putStrLn ""+        putStrLn $ show v++-------------------------------------------++-- short circuit+monoStep :: Double -> MaybeT (State Double) ()+monoStep d = do+    dp <- get+    when (d < dp) (fail "negative difference")+    put d+{-# INLINE monoStep #-}++isMonotoneIncreasing :: Vector Double -> Bool+isMonotoneIncreasing v =+    let res = evalState (runMaybeT $ (mapVectorM_ monoStep v)) (v @> 0)+     in case res of+        Nothing -> False+        Just _  -> True+++-------------------------------------------++-- | apply a test to successive elements of a vector, evaluates to true iff test passes for all pairs+successive_ :: Storable a => (a -> a -> Bool) -> Vector a -> Bool+successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ step (subVector 1 (dim v - 1) v))) (v @> 0)+   where step e = do+                  ep <- lift $ get+                  if t e ep+                     then lift $ put e+                     else (fail "successive_ test failed")++-- | operate on successive elements of a vector and return the resulting vector, whose length 1 less than that of the input+successive :: (Storable a, Storable b) => (a -> a -> b) -> Vector a -> Vector b+successive f v = evalState (mapVectorM step (subVector 1 (dim v - 1) v)) (v @> 0)+   where step e = do+                  ep <- get+                  put e+                  return $ f ep e++-------------------------------------------++v :: Vector Int+v = 10 |> [0..]++w = fromList ([1..10]++[10,9..1]) :: Vector Double+++main = do+    v' <- test1 v+    putStrLn ""+    putStrLn $ show v'+    evalStateT (addInitialM v) 0+    putStrLn ""+    putStrLn $ show (sumEvens v)+    evalStateT (sumEvensAndPrint v) 0+    putStrLn ""+    evalStateT (indexPlusSum v) 0+    putStrLn "-----------------------"+    mapVectorM_ print v+    print =<< (mapVectorM (const randomIO) v :: IO (Vector Double))+    print =<< (mapVectorM (\a -> fmap (+a) randomIO) (5|>[0,100..1000]) :: IO (Vector Double))+    putStrLn "-----------------------"+    print $ isMonotoneIncreasing w+    print $ isMonotoneIncreasing (subVector 0 7 w)+    print $ successive_ (>) v+    print $ successive_ (>) w+    print $ successive (+) v

@@ -1,6 +1,6 @@ -- $ runhaskell parallel.hs 2000 -import System(getArgs)+import System.Environment(getArgs) import Numeric.LinearAlgebra import Control.Parallel.Strategies import System.Time@@ -15,10 +15,10 @@ main = do     n <- (read . head) `fmap` getArgs     let m = ident n :: Matrix Double-    time $ print $ vectorMax $ takeDiag $ m <> 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 $ print $ maxElement $ takeDiag $ m <> m+    time $ print $ maxElement $ takeDiag $ parMul 2 m m+    time $ print $ maxElement $ takeDiag $ parMul 4 m m+    time $ print $ maxElement $ takeDiag $ parMul 8 m m  time act = do     t0 <- getClockTime

@@ -2,7 +2,7 @@  import Numeric.LinearAlgebra import System.Directory(doesFileExist)-import System(system)+import System.Process(system) import Control.Monad(when)  type Vec = Vector Double

@@ -3,7 +3,7 @@ import Numeric.LinearAlgebra import Graphics.Plot import System.Directory(doesFileExist)-import System(system)+import System.Process(system) import Control.Monad(when)  type Vec = Vector Double

@@ -14,7 +14,7 @@ fromVector v = unsafeFromForeignPtr p i n where     (p,i,n) = V.unsafeToForeignPtr v -toVector :: H.Vector t -> V.Vector t+toVector :: Storable t => H.Vector t -> V.Vector t toVector v = V.unsafeFromForeignPtr p i n where     (p,i,n) = unsafeToForeignPtr v @@ -22,11 +22,10 @@  v = V.slice 5 10 (V.fromList [1 .. 10::Double] V.++ V.replicate 10 7) -w = subVector 2 3 (linspace 10 (0,2))+w = subVector 2 3 (linspace 5 (0,1)) :: Vector Double  main = do     print v     print $ fromVector v     print w     print $ toVector w-

@@ -1,5 +1,5 @@ Name:               hmatrix-Version:            0.9.3.0+Version:            0.10.0.0 License:            GPL License-file:       LICENSE Author:             Alberto Ruiz@@ -11,11 +11,19 @@                     and other numerical computations, internally implemented using                     GSL, BLAS and LAPACK.                     .-                    See also hmatrix-special and hmatrix-glpk.+                    The Linear Algebra API is organized as follows:+                    .+                    - "Data.Packed": structure manipulation+                    .+                    - "Numeric.Container": simple numeric functions+                    .+                    - "Numeric.LinearAlgebra.Algorithms": matrix computations+                    .+                    - "Numeric.LinearAlgebra": everything + instances of standard Haskell numeric classes Category:           Math tested-with:        GHC ==6.10.4, GHC ==6.12.1 -cabal-version:      >=1.2+cabal-version:      >=1.6  build-type:         Custom @@ -45,6 +53,7 @@                     examples/devel/ej2/functions.c                     examples/Real.hs                     examples/vector.hs+                    examples/monadic.hs  extra-source-files: lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h,                     lib/Numeric/LinearAlgebra/LAPACK/clapack.h@@ -74,7 +83,8 @@     Build-Depends:      base >= 4 && < 5,                         array,                         storable-complex,-                        process+                        process,+                        binary      Extensions:         ForeignFunctionInterface,                         CPP@@ -91,35 +101,38 @@                         Numeric.GSL.Root,                         Numeric.GSL.Fitting,                         Numeric.GSL.ODE,-                        Numeric.GSL.Vector,                         Numeric.GSL,+                        Numeric.Container,                         Numeric.LinearAlgebra,                         Numeric.LinearAlgebra.LAPACK,-                        Numeric.LinearAlgebra.Interface,                         Numeric.LinearAlgebra.Algorithms,                         Graphics.Plot,-                     -- Data.Packed.Convert,                         Data.Packed.ST,-                        Data.Packed.Development,-                        Data.Packed.Random+                        Data.Packed.Development     other-modules:      Data.Packed.Internal,                         Data.Packed.Internal.Common,                         Data.Packed.Internal.Signatures,                         Data.Packed.Internal.Vector,                         Data.Packed.Internal.Matrix,-                        Numeric.LinearAlgebra.Linear,-                        Numeric.LinearAlgebra.Instances,-                        Numeric.GSL.Internal+                        Data.Packed.Random,+                        Numeric.GSL.Internal,+                        Numeric.GSL.Vector,+                        Numeric.Conversion,+                        Numeric.ContainerBoot,+                        Numeric.IO,+                        Numeric.Chain,+                        Numeric.Vector,+                        Numeric.Matrix      C-sources:          lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c,                         lib/Numeric/GSL/gsl-aux.c      if flag(vector)-       Build-Depends:   vector+       Build-Depends:   vector >= 0.7        cpp-options:     -DVECTOR      if flag(tests)-       Build-Depends:   QuickCheck, HUnit+       Build-Depends:   QuickCheck, HUnit, random        exposed-modules: Numeric.LinearAlgebra.Tests        other-modules:   Numeric.LinearAlgebra.Tests.Instances,                         Numeric.LinearAlgebra.Tests.Properties@@ -160,7 +173,6 @@     extra-libraries:     extra-lib-dirs: -    source-repository head-        type:     darcs-        location: http://code.haskell.org/hmatrix-+source-repository head+    type:     darcs+    location: http://code.haskell.org/hmatrix

@@ -1,14 +1,14 @@ ----------------------------------------------------------------------------- {- | Module      :  Data.Packed-Copyright   :  (c) Alberto Ruiz 2006-7+Copyright   :  (c) Alberto Ruiz 2006-2010 License     :  GPL-style  Maintainer  :  Alberto Ruiz (aruiz at um dot es) Stability   :  provisional Portability :  uses ffi -The Vector and Matrix types and some utilities.+Types for dense 'Vector' and 'Matrix' of 'Storable' elements.  -} -----------------------------------------------------------------------------@@ -16,12 +16,13 @@ module Data.Packed (     module Data.Packed.Vector,     module Data.Packed.Matrix,-    module Data.Packed.Random,-    module Data.Complex+--    module Numeric.Conversion,+--    module Data.Packed.Random,+--    module Data.Complex ) where  import Data.Packed.Vector import Data.Packed.Matrix-import Data.Packed.Random-import Data.Complex-+--import Data.Packed.Random+--import Data.Complex+--import Numeric.Conversion

@@ -17,7 +17,6 @@  module Data.Packed.Development (     createVector, createMatrix,-    Adapt,     vec, mat,     app1, app2, app3, app4,     app5, app6, app7, app8, app9, app10,

@@ -82,41 +82,27 @@  type Adapt f t r = t -> ((f -> r) -> IO()) -> IO() -app1 :: f-     -> Adapt f t (IO CInt)-     -> t-     -> String-     -> IO()--app2 :: f-     -> Adapt f t1 r-     -> t1-     -> Adapt r t2 (IO CInt)-     -> t2-     -> String-     -> IO()--app3 :: f-     -> Adapt f t1 r1-     -> t1-     -> Adapt r1 t2 r2-     -> t2-     -> Adapt r2 t3 (IO CInt)-     -> t3-     -> String-     -> IO()+type Adapt1 f t1 = Adapt f t1 (IO CInt) -> t1 -> String -> IO()+type Adapt2 f t1 r1 t2 = Adapt f t1 r1 -> t1 -> Adapt1 r1 t2+type Adapt3 f t1 r1 t2 r2 t3 = Adapt f t1 r1 -> t1 -> Adapt2 r1 t2 r2 t3+type Adapt4 f t1 r1 t2 r2 t3 r3 t4 = Adapt f t1 r1 -> t1 -> Adapt3 r1 t2 r2 t3 r3 t4+type Adapt5 f t1 r1 t2 r2 t3 r3 t4 r4 t5 = Adapt f t1 r1 -> t1 -> Adapt4 r1 t2 r2 t3 r3 t4 r4 t5+type Adapt6 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 = Adapt f t1 r1 -> t1 -> Adapt5 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6+type Adapt7 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 = Adapt f t1 r1 -> t1 -> Adapt6 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7+type Adapt8 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 = Adapt f t1 r1 -> t1 -> Adapt7 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8+type Adapt9 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 = Adapt f t1 r1 -> t1 -> Adapt8 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9+type Adapt10 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 r9 t10 = Adapt f t1 r1 -> t1 -> Adapt9 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 r9 t10 -app4 :: f-     -> Adapt f t1 r1-     -> t1-     -> Adapt r1 t2 r2-     -> t2-     -> Adapt r2 t3 r3-     -> t3-     -> Adapt r3 t4 (IO CInt)-     -> t4-     -> String-     -> IO()+app1 :: f -> Adapt1 f t1+app2 :: f -> Adapt2 f t1 r1 t2+app3 :: f -> Adapt3 f t1 r1 t2 r2 t3+app4 :: f -> Adapt4 f t1 r1 t2 r2 t3 r3 t4+app5 :: f -> Adapt5 f t1 r1 t2 r2 t3 r3 t4 r4 t5+app6 :: f -> Adapt6 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6+app7 :: f -> Adapt7 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7+app8 :: f -> Adapt8 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8+app9 :: f -> Adapt9 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9+app10 :: f -> Adapt10 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 r9 t10  app1 f w1 o1 s = w1 o1 $ \a1 -> f // a1 // check s app2 f w1 o1 w2 o2 s = ww2 w1 o1 w2 o2 $ \a1 a2 -> f // a1 // a2 // check s

@@ -29,7 +29,6 @@     liftMatrix, liftMatrix2,     (@@>),     saveMatrix,-    fromComplexV, toComplexV, conjV,     singleton ) where @@ -76,7 +75,11 @@  data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq) --- | Matrix representation suitable for GSL and LAPACK computations.+{- | Matrix representation suitable for GSL and LAPACK computations.++The elements are stored in a continuous memory array.++-} data Matrix t = MC { irows :: {-# UNPACK #-} !Int                    , icols :: {-# UNPACK #-} !Int                    , cdat :: {-# UNPACK #-} !(Vector t) }@@ -222,13 +225,13 @@  , 9.0, 10.0, 11.0, 12.0 ]@  -}-reshape :: Element t => Int -> Vector t -> Matrix t+reshape :: Storable t => Int -> Vector t -> Matrix t reshape c v = matrixFromVector RowMajor c v  singleton x = reshape 1 (fromList [x])  -- | application of a vector function on the flattened matrix elements-liftMatrix :: (Element a, Element b) => (Vector a -> Vector b) -> Matrix a -> Matrix b+liftMatrix :: (Storable a, Storable b) => (Vector a -> Vector b) -> Matrix a -> Matrix b liftMatrix f MC { icols = c, cdat = d } = matrixFromVector RowMajor    c (f d) liftMatrix f MF { icols = c, fdat = d } = matrixFromVector ColumnMajor c (f d) @@ -246,21 +249,37 @@  ------------------------------------------------------------------ --- | Auxiliary class.-class (Storable a, Floating a) => Element a where+{- | Supported matrix elements.++    This class provides optimized internal+    operations for selected element types.+    It provides unoptimised defaults for any 'Storable' type,+    so you can create instances simply as:+    @instance Element Foo@.+-}+class (Storable a) => Element a where     subMatrixD :: (Int,Int) -- ^ (r0,c0) starting position                 -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix                -> Matrix a -> Matrix a     subMatrixD = subMatrix'     transdata :: Int -> Vector a -> Int -> Vector a-    transdata = transdata'+    transdata = transdataP -- transdata'     constantD  :: a -> Int -> Vector a-    constantD = constant'+    constantD = constantP -- constant' ++instance Element Float where+    transdata  = transdataAux ctransF+    constantD  = constantAux cconstantF+ instance Element Double where     transdata  = transdataAux ctransR     constantD  = constantAux cconstantR +instance Element (Complex Float) where+    transdata  = transdataAux ctransQ+    constantD  = constantAux cconstantQ+ instance Element (Complex Double) where     transdata  = transdataAux ctransC     constantD  = constantAux cconstantC@@ -308,8 +327,27 @@         r2 = dim d `div` c2         noneed = r1 == 1 || c1 == 1 +transdataP :: Storable a => Int -> Vector a -> Int -> Vector a+transdataP c1 d c2 =+    if noneed+       then d+       else unsafePerformIO $ do+          v <- createVector (dim d)+          unsafeWith d $ \pd ->+              unsafeWith v $ \pv ->+                  ctransP (fi r1) (fi c1) (castPtr pd) (fi sz) (fi r2) (fi c2) (castPtr pv) (fi sz) // check "transdataP"+          return v+   where r1 = dim d `div` c1+         r2 = dim d `div` c2+         sz = sizeOf (d @> 0)+         noneed = r1 == 1 || c1 == 1++foreign import ccall "transF" ctransF :: TFMFM foreign import ccall "transR" ctransR :: TMM+foreign import ccall "transQ" ctransQ :: TQMQM foreign import ccall "transC" ctransC :: TCMCM+foreign import ccall "transP" ctransP :: CInt -> CInt -> Ptr () -> CInt -> CInt -> CInt -> Ptr () -> CInt -> IO CInt+ ----------------------------------------------------------------------  constant' v n = unsafePerformIO $ do@@ -329,13 +367,33 @@     free px     return v +constantF :: Float -> Int -> Vector Float+constantF = constantAux cconstantF+foreign import ccall "constantF" cconstantF :: Ptr Float -> TF+ constantR :: Double -> Int -> Vector Double constantR = constantAux cconstantR foreign import ccall "constantR" cconstantR :: Ptr Double -> TV +constantQ :: Complex Float -> Int -> Vector (Complex Float)+constantQ = constantAux cconstantQ+foreign import ccall "constantQ" cconstantQ :: Ptr (Complex Float) -> TQV+ constantC :: Complex Double -> Int -> Vector (Complex Double) constantC = constantAux cconstantC foreign import ccall "constantC" cconstantC :: Ptr (Complex Double) -> TCV++constantP :: Storable a => a -> Int -> Vector a+constantP a n = unsafePerformIO $ do+    let sz = sizeOf a+    v <- createVector n+    unsafeWith v $ \p -> do+       alloca $ \k -> do+                      poke k a+                      cconstantP (castPtr k) (fi n) (castPtr p) (fi sz) // check "constantP"+    return v+foreign import ccall "constantP" cconstantP :: Ptr () -> CInt -> Ptr () -> CInt -> IO CInt+ ----------------------------------------------------------------------  -- | Extracts a submatrix from a matrix.@@ -364,21 +422,6 @@  subMatrix' (r0,c0) (rt,ct) (MC _r c v) = MC rt ct $ subMatrix'' (r0,c0) (rt,ct) c v subMatrix' (r0,c0) (rt,ct) m = trans $ subMatrix' (c0,r0) (ct,rt) (trans m)-------------------------------------------------------------------------------- | obtains the complex conjugate of a complex vector-conjV :: Vector (Complex Double) -> Vector (Complex Double)-conjV = mapVector conjugate---- | creates a complex vector from vectors with real and imaginary parts-toComplexV :: (Vector Double, Vector Double) ->  Vector (Complex Double)-toComplexV (r,i) = asComplex $ flatten $ fromColumns [r,i]---- | the inverse of 'toComplex'-fromComplexV :: Vector (Complex Double) -> (Vector Double, Vector Double)-fromComplexV z = (r,i) where-    [r,i] = toColumns $ reshape 2 $ asReal z  -------------------------------------------------------------------------- 

@@ -18,11 +18,21 @@ import Data.Complex import Foreign.C.Types +type PF = Ptr Float                             -- type PD = Ptr Double                            --+type PQ = Ptr (Complex Float)                   -- type PC = Ptr (Complex Double)                  --+type TF = CInt -> PF -> IO CInt                 --+type TFF = CInt -> PF -> TF                     --+type TFV = CInt -> PF -> TV                     --+type TVF = CInt -> PD -> TF                     --+type TFFF = CInt -> PF -> TFF                   -- type TV = CInt -> PD -> IO CInt                 -- type TVV = CInt -> PD -> TV                     -- type TVVV = CInt -> PD -> TVV                   --+type TFM = CInt -> CInt -> PF -> IO CInt        --+type TFMFM =  CInt -> CInt -> PF -> TFM         --+type TFMFMFM =  CInt -> CInt -> PF -> TFMFM     -- type TM = CInt -> CInt -> PD -> IO CInt         -- type TMM =  CInt -> CInt -> PD -> TM            -- type TVMM = CInt -> PD -> TMM                   --@@ -47,6 +57,14 @@ type TCV = CInt -> PC -> IO CInt                -- type TCVCV = CInt -> PC -> TCV                  -- type TCVCVCV = CInt -> PC -> TCVCV              --+type TCVV = CInt -> PC -> TV                    --+type TQV = CInt -> PQ -> IO CInt                --+type TQVQV = CInt -> PQ -> TQV                  --+type TQVQVQV = CInt -> PQ -> TQVQV              --+type TQVF = CInt -> PQ -> TF                    --+type TQM = CInt -> CInt -> PQ -> IO CInt        --+type TQMQM = CInt -> CInt -> PQ -> TQM          --+type TQMQMQM = CInt -> CInt -> PQ -> TQMQM      -- type TCMCV = CInt -> CInt -> PC -> TCV          -- type TVCV = CInt -> PD -> TCV                   -- type TCVM = CInt -> PC -> TM                    --

@@ -1,4 +1,4 @@-{-# LANGUAGE MagicHash, CPP, UnboxedTuples, BangPatterns #-}+{-# LANGUAGE MagicHash, CPP, UnboxedTuples, BangPatterns, FlexibleContexts #-} ----------------------------------------------------------------------------- -- | -- Module      :  Data.Packed.Internal.Vector@@ -17,10 +17,11 @@     Vector, dim,     fromList, toList, (|>),     join, (@>), safe, at, at', subVector, takesV,-    mapVector, zipVector,-    foldVector, foldVectorG, foldLoop,+    mapVector, zipVectorWith, unzipVectorWith,+    mapVectorM, mapVectorM_, mapVectorWithIndexM, mapVectorWithIndexM_,+    foldVector, foldVectorG, foldLoop, foldVectorWithIndex,     createVector, vec,-    asComplex, asReal,+    asComplex, asReal, float2DoubleV, double2FloatV,     fwriteVector, freadVector, fprintfVector, fscanfVector,     cloneVector,     unsafeToForeignPtr,@@ -44,7 +45,7 @@  import GHC.Base #if __GLASGOW_HASKELL__ < 612-import GHC.IOBase+import GHC.IOBase hiding (liftIO) #endif  #ifdef VECTOR@@ -70,11 +71,11 @@       , fptr :: {-# UNPACK #-} !(ForeignPtr t)   -- ^ foreign pointer to the memory block       } -unsafeToForeignPtr :: Vector a -> (ForeignPtr a, Int, Int)+unsafeToForeignPtr :: Storable a => Vector a -> (ForeignPtr a, Int, Int) unsafeToForeignPtr v = (fptr v, ioff v, idim v)  -- | Same convention as in Roman Leshchinskiy's vector package.-unsafeFromForeignPtr :: ForeignPtr a -> Int -> Int -> Vector a+unsafeFromForeignPtr :: Storable a => ForeignPtr a -> Int -> Int -> Vector a unsafeFromForeignPtr fp i n | n > 0 = V {ioff = i, idim = n, fptr = fp}                             | otherwise = error "unsafeFromForeignPtr with dim < 1" @@ -264,17 +265,33 @@ ---------------------------------------------------------------  -- | transforms a complex vector into a real vector with alternating real and imaginary parts -asReal :: Vector (Complex Double) -> Vector Double---asReal v = V { ioff = 2*ioff v, idim = 2*dim v, fptr =  castForeignPtr (fptr v) }+asReal :: (RealFloat a, Storable a) => Vector (Complex a) -> Vector a asReal v = unsafeFromForeignPtr (castForeignPtr fp) (2*i) (2*n)     where (fp,i,n) = unsafeToForeignPtr v  -- | transforms a real vector into a complex vector with alternating real and imaginary parts-asComplex :: Vector Double -> Vector (Complex Double)---asComplex v = V { ioff = ioff v `div` 2, idim = dim v `div` 2, fptr =  castForeignPtr (fptr v) }+asComplex :: (RealFloat a, Storable a) => Vector a -> Vector (Complex a) asComplex v = unsafeFromForeignPtr (castForeignPtr fp) (i `div` 2) (n `div` 2)     where (fp,i,n) = unsafeToForeignPtr v +---------------------------------------------------------------++float2DoubleV :: Vector Float -> Vector Double+float2DoubleV v = unsafePerformIO $ do+    r <- createVector (dim v)+    app2 c_float2double vec v vec r "float2double"+    return r++double2FloatV :: Vector Double -> Vector Float+double2FloatV v = unsafePerformIO $ do+    r <- createVector (dim v)+    app2 c_double2float vec v vec r "double2float2"+    return r+++foreign import ccall "float2double" c_float2double:: TFV+foreign import ccall "double2float" c_double2float:: TVF+ ----------------------------------------------------------------  cloneVector :: Storable t => Vector t -> IO (Vector t)@@ -302,8 +319,8 @@ {-# INLINE mapVector #-}  -- | zipWith for Vectors-zipVector :: (Storable a, Storable b, Storable c) => (a-> b -> c) -> Vector a -> Vector b -> Vector c-zipVector f u v = unsafePerformIO $ do+zipVectorWith :: (Storable a, Storable b, Storable c) => (a-> b -> c) -> Vector a -> Vector b -> Vector c+zipVectorWith f u v = unsafePerformIO $ do     let n = min (dim u) (dim v)     w <- createVector n     unsafeWith u $ \pu ->@@ -316,16 +333,47 @@                                go (k-1)                 go (n -1)     return w-{-# INLINE zipVector #-}+{-# INLINE zipVectorWith #-} +-- | unzipWith for Vectors+unzipVectorWith :: (Storable (a,b), Storable c, Storable d) +                   => ((a,b) -> (c,d)) -> Vector (a,b) -> (Vector c,Vector d)+unzipVectorWith f u = unsafePerformIO $ do+      let n = dim u+      v <- createVector n+      w <- createVector n+      unsafeWith u $ \pu ->+          unsafeWith v $ \pv ->+              unsafeWith w $ \pw -> do+                  let go (-1) = return ()+                      go !k   = do z <- peekElemOff pu k+                                   let (x,y) = f z +                                   pokeElemOff      pv k x+                                   pokeElemOff      pw k y+                                   go (k-1)+                  go (n-1)+      return (v,w)+{-# INLINE unzipVectorWith #-}++foldVector :: Storable a => (a -> b -> b) -> b -> Vector a -> b foldVector f x v = unsafePerformIO $-    unsafeWith (v::Vector Double) $ \p -> do+    unsafeWith v $ \p -> do         let go (-1) s = return s             go !k !s = do y <- peekElemOff p k                           go (k-1::Int) (f y s)         go (dim v -1) x {-# INLINE foldVector #-} +-- the zero-indexed index is passed to the folding function+foldVectorWithIndex :: Storable a => (Int -> a -> b -> b) -> b -> Vector a -> b+foldVectorWithIndex f x v = unsafePerformIO $+    unsafeWith v $ \p -> do+        let go (-1) s = return s+            go !k !s = do y <- peekElemOff p k+                          go (k-1::Int) (f k y s)+        go (dim v -1) x+{-# INLINE foldVectorWithIndex #-}+ foldLoop f s0 d = go (d - 1) s0      where        go 0 s = f (0::Int) s@@ -338,6 +386,75 @@  ------------------------------------------------------------------- +-- | monadic map over Vectors+--    the monad @m@ must be strict+mapVectorM :: (Storable a, Storable b, Monad m) => (a -> m b) -> Vector a -> m (Vector b)+mapVectorM f v = do+    w <- return $! unsafePerformIO $! createVector (dim v)+    mapVectorM' w 0 (dim v -1)+    return w+    where mapVectorM' w' !k !t+              | k == t               = do+                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k +                                       y <- f x+                                       return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y+              | otherwise            = do+                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k +                                       y <- f x+                                       _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y+                                       mapVectorM' w' (k+1) t+{-# INLINE mapVectorM #-}++-- | monadic map over Vectors+mapVectorM_ :: (Storable a, Monad m) => (a -> m ()) -> Vector a -> m ()+mapVectorM_ f v = do+    mapVectorM' 0 (dim v -1)+    where mapVectorM' !k !t+              | k == t            = do+                                    x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k+                                    f x+              | otherwise         = do+                                    x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k +                                    _ <- f x+                                    mapVectorM' (k+1) t+{-# INLINE mapVectorM_ #-}++-- | monadic map over Vectors with the zero-indexed index passed to the mapping function+--    the monad @m@ must be strict+mapVectorWithIndexM :: (Storable a, Storable b, Monad m) => (Int -> a -> m b) -> Vector a -> m (Vector b)+mapVectorWithIndexM f v = do+    w <- return $! unsafePerformIO $! createVector (dim v)+    mapVectorM' w 0 (dim v -1)+    return w+    where mapVectorM' w' !k !t+              | k == t               = do+                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k +                                       y <- f k x+                                       return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y+              | otherwise            = do+                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k +                                       y <- f k x+                                       _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y+                                       mapVectorM' w' (k+1) t+{-# INLINE mapVectorWithIndexM #-}++-- | monadic map over Vectors with the zero-indexed index passed to the mapping function+mapVectorWithIndexM_ :: (Storable a, Monad m) => (Int -> a -> m ()) -> Vector a -> m ()+mapVectorWithIndexM_ f v = do+    mapVectorM' 0 (dim v -1)+    where mapVectorM' !k !t+              | k == t            = do+                                    x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k+                                    f k x+              | otherwise         = do+                                    x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k +                                    _ <- f k x+                                    mapVectorM' (k+1) t+{-# INLINE mapVectorWithIndexM_ #-}++-------------------------------------------------------------------++ -- | Loads a vector from an ASCII file (the number of elements must be known in advance). fscanfVector :: FilePath -> Int -> IO (Vector Double) fscanfVector filename n = do@@ -379,3 +496,4 @@     free charname  foreign import ccall "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV+

@@ -1,8 +1,13 @@-{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, FlexibleContexts #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FunctionalDependencies #-}+ ----------------------------------------------------------------------------- -- | -- Module      :  Data.Packed.Matrix--- Copyright   :  (c) Alberto Ruiz 2007+-- Copyright   :  (c) Alberto Ruiz 2007-10 -- License     :  GPL-style -- -- Maintainer  :  Alberto Ruiz <aruiz@um.es>@@ -11,11 +16,14 @@ -- -- A Matrix representation suitable for numerical computations using LAPACK and GSL. --+-- This module provides basic functions for manipulation of structure.+ -----------------------------------------------------------------------------  module Data.Packed.Matrix (-    Element, Container(..),-    Matrix,rows,cols,+    Matrix,+    Element,+    rows,cols,     (><),     trans,     reshape, flatten,@@ -28,22 +36,67 @@     flipud, fliprl,     subMatrix, takeRows, dropRows, takeColumns, dropColumns,     extractRows,-    ident, diag, diagRect, takeDiag,-    liftMatrix, liftMatrix2, liftMatrix2Auto,-    dispf, disps, dispcf, latexFormat, format,-    loadMatrix, saveMatrix, fromFile, fileDimensions,-    readMatrix, fromArray2D+    diagRect, takeDiag,+    liftMatrix, liftMatrix2, liftMatrix2Auto,fromArray2D ) where  import Data.Packed.Internal import qualified Data.Packed.ST as ST-import Data.Packed.Vector-import Data.Array-import System.Process(readProcess)-import Text.Printf(printf) import Data.List(transpose,intersperse)-import Data.Complex+import Data.Array ++import Data.Binary+import Foreign.Storable+import Control.Monad(replicateM)+--import Control.Arrow((***))+--import GHC.Float(double2Float,float2Double)+++-------------------------------------------------------------------++instance (Binary a, Element a, Storable a) => Binary (Matrix a) where+    put m = do+            let r = rows m+            let c = cols m+            put r+            put c+            mapM_ (\i -> mapM_ (\j -> put $ m @@> (i,j)) [0..(c-1)]) [0..(r-1)]+    get = do+          r <- get+          c <- get+          xs <- replicateM r $ replicateM c get+          return $ fromLists xs++-------------------------------------------------------------------++instance (Show a, Element a) => (Show (Matrix a)) where+    show m = (sizes++) . dsp . map (map show) . toLists $ m+        where sizes = "("++show (rows m)++"><"++show (cols m)++")\n"++dsp as = (++" ]") . (" ["++) . init . drop 2 . unlines . map (" , "++) . map unwords' $ transpose mtp+    where+        mt = transpose as+        longs = map (maximum . map length) mt+        mtp = zipWith (\a b -> map (pad a) b) longs mt+        pad n str = replicate (n - length str) ' ' ++ str+        unwords' = concat . intersperse ", "++------------------------------------------------------------------++instance (Element a, Read a) => Read (Matrix a) where+    readsPrec _ s = [((rs><cs) . read $ listnums, rest)]+        where (thing,rest) = breakAt ']' s+              (dims,listnums) = breakAt ')' thing+              cs = read . init . fst. breakAt ')' . snd . breakAt '<' $ dims+              rs = read . snd . breakAt '(' .init . fst . breakAt '>' $ dims+++breakAt c l = (a++[c],tail b) where+    (a,b) = break (==c) l++------------------------------------------------------------------+ -- | creates a matrix from a vertical list of matrices joinVert :: Element t => [Matrix t] -> Matrix t joinVert ms = case common cols ms of@@ -92,7 +145,7 @@      g [Just nr,Just nc] m                 | nr == r && nc == c = m-                | r == 1 && c == 1 = reshape nc (constant x (nr*nc))+                | r == 1 && c == 1 = reshape nc (constantD x (nr*nc))                 | r == 1 = fromRows (replicate nr (flatten m))                 | otherwise = fromColumns (replicate nc (flatten m))       where@@ -113,28 +166,19 @@  ------------------------------------------------------------ --- | Creates a square matrix with a given diagonal.-diag :: Element a => Vector a -> Matrix a-diag v = ST.runSTMatrix $ do-    let d = dim v-    m <- ST.newMatrix 0 d d-    mapM_ (\k -> ST.writeMatrix m k k (v@>k)) [0..d-1]-    return m--{- | creates a rectangular diagonal matrix+{- | creates a rectangular diagonal matrix: -@> diagRect (constant 5 3) 3 4 :: Matrix Double-(3><4)- [ 5.0, 0.0, 0.0, 0.0- , 0.0, 5.0, 0.0, 0.0- , 0.0, 0.0, 5.0, 0.0 ]@+@> diagRect 7 (fromList [10,20,30]) 4 5 :: Matrix Double+(4><5)+ [ 10.0,  7.0,  7.0, 7.0, 7.0+ ,  7.0, 20.0,  7.0, 7.0, 7.0+ ,  7.0,  7.0, 30.0, 7.0, 7.0+ ,  7.0,  7.0,  7.0, 7.0, 7.0 ]@ -}-diagRect :: (Element t, Num t) => Vector t -> Int -> Int -> Matrix t-diagRect v r c-    | dim v < min r c = error "diagRect called with dim v < min r c"-    | otherwise = ST.runSTMatrix $ do-        m <- ST.newMatrix 0 r c-        let d = min r c+diagRect :: (Storable t) => t -> Vector t -> Int -> Int -> Matrix t+diagRect z v r c = ST.runSTMatrix $ do+        m <- ST.newMatrix z r c+        let d = min r c `min` (dim v)         mapM_ (\k -> ST.writeMatrix m k k (v@>k)) [0..d-1]         return m @@ -142,10 +186,6 @@ takeDiag :: (Element t) => Matrix t -> Vector t takeDiag m = fromList [flatten m `at` (k*cols m+k) | k <- [0 .. min (rows m) (cols m) -1]] --- | creates the identity matrix of given dimension-ident :: Element a => Int -> Matrix a-ident n = diag (constant 1 n)- ------------------------------------------------------------  {- | An easy way to create a matrix:@@ -169,7 +209,7 @@  , 4.0, 5.0, 6.0 ]@  -}-(><) :: (Element a) => Int -> Int -> [a] -> Matrix a+(><) :: (Storable a) => Int -> Int -> [a] -> Matrix a r >< c = f where     f l | dim v == r*c = matrixFromVector RowMajor c v         | otherwise    = error $ "inconsistent list size = "@@ -205,22 +245,28 @@ fromLists = fromRows . map fromList  -- | creates a 1-row matrix from a vector-asRow :: Element a => Vector a -> Matrix a+asRow :: Storable a => Vector a -> Matrix a asRow v = reshape (dim v) v  -- | creates a 1-column matrix from a vector-asColumn :: Element a => Vector a -> Matrix a+asColumn :: Storable a => Vector a -> Matrix a asColumn v = reshape 1 v  + {- | creates a Matrix of the specified size using the supplied function to-     to map the row/column position to the value at that row/column position.+     to map the row\/column position to the value at that row\/column position.  @> buildMatrix 3 4 (\ (r,c) -> fromIntegral r * fromIntegral c) (3><4)  [ 0.0, 0.0, 0.0, 0.0, 0.0  , 0.0, 1.0, 2.0, 3.0, 4.0  , 0.0, 2.0, 4.0, 6.0, 8.0]@++Hilbert matrix of order N:++@hilb n = buildMatrix n n (\(i,j)->1/(fromIntegral i + fromIntegral j +1))@+ -} buildMatrix :: Element a => Int -> Int -> ((Int, Int) -> a) -> Matrix a buildMatrix rc cc f =@@ -229,135 +275,13 @@  ----------------------------------------------------- -fromArray2D :: (Element e) => Array (Int, Int) e -> Matrix e+fromArray2D :: (Storable e) => Array (Int, Int) e -> Matrix e fromArray2D m = (r><c) (elems m)     where ((r0,c0),(r1,c1)) = bounds m           r = r1-r0+1           c = c1-c0+1  ----------------------------------------------------------------------- display utilities---{- | Creates a string from a matrix given a separator and a function to show each entry. Using-this function the user can easily define any desired display function:--@import Text.Printf(printf)@--@disp = putStr . format \"  \" (printf \"%.2f\")@---}-format :: (Element t) => String -> (t -> String) -> Matrix t -> String-format sep f m = table sep . map (map f) . toLists $ m--{- | Show a matrix with \"autoscaling\" and a given number of decimal places.--@disp = putStr . disps 2--\> disp $ 120 * (3><4) [1..]-3x4  E3- 0.12  0.24  0.36  0.48- 0.60  0.72  0.84  0.96- 1.08  1.20  1.32  1.44-@--}-disps :: Int -> Matrix Double -> String-disps d x = sdims x ++ "  " ++ formatScaled d x--{- | Show a matrix with a given number of decimal places.--@disp = putStr . dispf 3--\> disp (1/3 + ident 4)-4x4-1.333  0.333  0.333  0.333-0.333  1.333  0.333  0.333-0.333  0.333  1.333  0.333-0.333  0.333  0.333  1.333-@--}-dispf :: Int -> Matrix Double -> String-dispf d x = sdims x ++ "\n" ++ formatFixed (if isInt x then 0 else d) x--sdims x = show (rows x) ++ "x" ++ show (cols x)--formatFixed d x = format "  " (printf ("%."++show d++"f")) $ x--isInt = all lookslikeInt . toList . flatten--formatScaled dec t = "E"++show o++"\n" ++ ss-    where ss = format " " (printf fmt. g) t-          g x | o >= 0    = x/10^(o::Int)-              | otherwise = x*10^(-o)-          o = floor $ maximum $ map (logBase 10 . abs) $ toList $ flatten t-          fmt = '%':show (dec+3) ++ '.':show dec ++"f"---- | Tool to display matrices with latex syntax.-latexFormat :: String -- ^ type of braces: \"matrix\", \"bmatrix\", \"pmatrix\", etc.-            -> String -- ^ Formatted matrix, with elements separated by spaces and newlines-            -> String-latexFormat del tab = "\\begin{"++del++"}\n" ++ f tab ++ "\\end{"++del++"}"-    where f = unlines . intersperse "\\\\" . map unwords . map (intersperse " & " . words) . tail . lines---- | Pretty print a complex number with at most n decimal digits.-showComplex :: Int -> Complex Double -> String-showComplex d (a:+b)-    | isZero a && isZero b = "0"-    | isZero b = sa-    | isZero a && isOne b = s2++"i"-    | isZero a = sb++"i"-    | isOne b = sa++s3++"i"-    | otherwise = sa++s1++sb++"i"-  where-    sa = shcr d a-    sb = shcr d b-    s1 = if b<0 then "" else "+"-    s2 = if b<0 then "-" else ""-    s3 = if b<0 then "-" else "+"--shcr d a | lookslikeInt a = printf "%.0f" a-         | otherwise      = printf ("%."++show d++"f") a---lookslikeInt x = show (round x :: Int) ++".0" == shx || "-0.0" == shx-   where shx = show x--isZero x = show x `elem` ["0.0","-0.0"]-isOne  x = show x `elem` ["1.0","-1.0"]---- | Pretty print a complex matrix with at most n decimal digits.-dispcf :: Int -> Matrix (Complex Double) -> String-dispcf d m = sdims m ++ "\n" ++ format "  " (showComplex d) m-------------------------------------------------------------------------- | reads a matrix from a string containing a table of numbers.-readMatrix :: String -> Matrix Double-readMatrix = fromLists . map (map read). map words . filter (not.null) . lines--{- |  obtains the number of rows and columns in an ASCII data file-      (provisionally using unix's wc).--}-fileDimensions :: FilePath -> IO (Int,Int)-fileDimensions fname = do-    wcres <- readProcess "wc" ["-w",fname] ""-    contents <- readFile fname-    let tot = read . head . words $ wcres-        c   = length . head . dropWhile null . map words . lines $ contents-    if tot > 0-        then return (tot `div` c, c)-        else return (0,0)---- | Loads a matrix from an ASCII file formatted as a 2D table.-loadMatrix :: FilePath -> IO (Matrix Double)-loadMatrix file = fromFile file =<< fileDimensions file---- | Loads a matrix from an ASCII file (the number of rows and columns must be known in advance).-fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double)-fromFile filename (r,c) = reshape c `fmap` fscanfVector filename (r*c)-- -- | rearranges the rows of a matrix according to the order given in a list of integers.  extractRows :: Element t => [Int] -> Matrix t -> Matrix t extractRows l m = fromRows $ extract (toRows $ m) l@@ -424,47 +348,3 @@     cs = replicate qc c ++ if rc > 0 then [rc] else []  ----------------------------------------------------------------------- | conversion utilities-class (Element e) => Container c e where-    toComplex   :: RealFloat e => (c e, c e) -> c (Complex e)-    fromComplex :: RealFloat e => c (Complex e) -> (c e, c e)-    comp        :: RealFloat e => c e -> c (Complex e)-    conj        :: RealFloat e => c (Complex e) -> c (Complex e)-    real        :: c Double -> c e-    complex     :: c e -> c (Complex Double)--instance Container Vector Double where-    toComplex = toComplexV-    fromComplex = fromComplexV-    comp v = toComplex (v,constant 0 (dim v))-    conj = conjV-    real = id-    complex = comp--instance Container Vector (Complex Double) where-    toComplex = undefined -- can't match-    fromComplex = undefined-    comp = undefined-    conj = undefined-    real = comp-    complex = id--instance Container Matrix Double where-    toComplex = uncurry $ liftMatrix2 $ curry toComplex-    fromComplex z = (reshape c r, reshape c i)-        where (r,i) = fromComplex (flatten z)-              c = cols z-    comp = liftMatrix comp-    conj = liftMatrix conj-    real = id-    complex = comp--instance Container Matrix (Complex Double) where-    toComplex = undefined-    fromComplex = undefined-    comp = undefined-    conj = undefined-    real = comp-    complex = id-

@@ -20,11 +20,11 @@ ) where  import Numeric.GSL.Vector-import Data.Packed.Matrix-import Data.Packed.Vector+import Data.Packed+import Numeric.ContainerBoot import Numeric.LinearAlgebra.Algorithms-import Numeric.LinearAlgebra.Interface + -- | Obtains a matrix whose rows are pseudorandom samples from a multivariate -- Gaussian distribution. gaussianSample :: Int -- ^ seed@@ -34,9 +34,9 @@                -> Matrix Double -- ^ result gaussianSample seed n med cov = m where     c = dim med-    meds = constant 1 n `outer` med+    meds = konst 1 n `outer` med     rs = reshape c $ randomVector seed Gaussian (c * n)-    m = rs <> cholSH cov + meds+    m = rs `mXm` cholSH cov `add` meds  -- | Obtains a matrix whose rows are pseudorandom samples from a multivariate -- uniform distribution.@@ -50,8 +50,8 @@     cs = zipWith subtract as bs     d = dim a     dat = toRows $ reshape n $ randomVector seed Uniform (n*d)-    am = constant 1 n `outer` a-    m = fromColumns (zipWith scale cs dat) + am+    am = konst 1 n `outer` a+    m = fromColumns (zipWith scale cs dat) `add` am  ------------ utilities ------------------------------- @@ -60,7 +60,7 @@ meanCov x = (med,cov) where     r    = rows x     k    = 1 / fromIntegral r-    med  = constant k r <> x-    meds = constant 1 r `outer` med-    xc   = x - meds-    cov  = (trans xc <> xc) / fromIntegral (r-1)+    med  = konst k r `vXm` x+    meds = konst 1 r `outer` med+    xc   = x `sub` meds+    cov  = flip scale (trans xc `mXm` xc) (recip (fromIntegral (r-1)))

@@ -90,11 +90,11 @@ writeVector = safeIndexV unsafeWriteVector  {-# NOINLINE newUndefinedVector #-}-newUndefinedVector :: Element t => Int -> ST s (STVector s t)+newUndefinedVector :: Storable t => Int -> ST s (STVector s t) newUndefinedVector = unsafeIOToST . fmap STVector . createVector  {-# INLINE newVector #-}-newVector :: Element t => t -> Int -> ST s (STVector s t)+newVector :: Storable t => t -> Int -> ST s (STVector s t) newVector x n = do     v <- newUndefinedVector n     let go (-1) = return v@@ -164,9 +164,9 @@ writeMatrix = safeIndexM unsafeWriteMatrix  {-# NOINLINE newUndefinedMatrix #-}-newUndefinedMatrix :: Element t => MatrixOrder -> Int -> Int -> ST s (STMatrix s t)+newUndefinedMatrix :: Storable 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 :: Storable t => t -> Int -> Int -> ST s (STMatrix s t) newMatrix v r c = unsafeThawMatrix $ reshape c $ runSTVector $ newVector v (r*c)

@@ -1,3 +1,4 @@+{-# LANGUAGE FlexibleContexts #-} ----------------------------------------------------------------------------- -- | -- Module      :  Data.Packed.Vector@@ -10,6 +11,8 @@ -- -- 1D arrays suitable for numeric computations using external libraries. --+-- This module provides basic functions for manipulation of structure.+-- -----------------------------------------------------------------------------  module Data.Packed.Vector (@@ -17,53 +20,47 @@     fromList, (|>), toList, buildVector,     dim, (@>),     subVector, takesV, join,-    constant, linspace,-    vecdisp,-    vectorMax, vectorMin, vectorMaxIndex, vectorMinIndex,-    mapVector, zipVector,-    fscanfVector, fprintfVector, freadVector, fwriteVector,-    foldLoop, foldVector, foldVectorG+    mapVector, zipVector, zipVectorWith, unzipVector, unzipVectorWith,+    mapVectorM, mapVectorM_, mapVectorWithIndexM, mapVectorWithIndexM_,+    foldLoop, foldVector, foldVectorG, foldVectorWithIndex ) where -import Data.Packed.Internal-import Numeric.GSL.Vector--- import Data.Packed.ST--{- | Creates a real vector containing a range of values:--@\> linspace 5 (-3,7)-5 |> [-3.0,-0.5,2.0,4.5,7.0]@--Logarithmic spacing can be defined as follows:+import Data.Packed.Internal.Vector+import Data.Binary+import Foreign.Storable+import Control.Monad(replicateM) -@logspace n (a,b) = 10 ** linspace n (a,b)@--}-linspace :: Int -> (Double, Double) -> Vector Double-linspace n (a,b) = add a $ scale s  $ fromList [0 .. fromIntegral n-1]-    where scale = vectorMapValR Scale-          add   = vectorMapValR AddConstant-          s = (b-a)/fromIntegral (n-1)+------------------------------------------------------------------- -vectorMax :: Vector Double -> Double-vectorMax = toScalarR Max+-- a 64K cache, with a Double taking 13 bytes in Bytestring,+-- implies a chunk size of 5041+chunk :: Int+chunk = 5000 -vectorMin :: Vector Double -> Double-vectorMin = toScalarR Min+chunks :: Int -> [Int]+chunks d = let c = d `div` chunk+               m = d `mod` chunk+           in if m /= 0 then reverse (m:(replicate c chunk)) else (replicate c chunk)   -vectorMaxIndex :: Vector Double -> Int-vectorMaxIndex = round . toScalarR MaxIdx+putVector v = do+              let d = dim v+              mapM_ (\i -> put $ v @> i) [0..(d-1)] -vectorMinIndex :: Vector Double -> Int-vectorMinIndex = round . toScalarR MinIdx+getVector d = do+              xs <- replicateM d get+              return $! fromList xs -{- | creates a vector with a given number of equal components:+instance (Binary a, Storable a) => Binary (Vector a) where+    put v = do+            let d = dim v+            put d+            mapM_ putVector $! takesV (chunks d) v+    get = do+          d <- get+          vs <- mapM getVector $ chunks d+          return $! join vs -@> 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)-constant = constantD -- about 2x faster+-------------------------------------------------------------------  {- | creates a Vector of the specified length using the supplied function to      to map the index to the value at that index.@@ -72,22 +69,17 @@ 4 |> [0.0,1.0,2.0,3.0]@  -}-buildVector :: Element a => Int -> (Int -> a) -> Vector a+buildVector :: Storable a => Int -> (Int -> a) -> Vector a buildVector len f =     fromList $ map f [0 .. (len - 1)]  -{- | Show a vector using a function for showing matrices.+-- | zip for Vectors+zipVector :: (Storable a, Storable b, Storable (a,b)) => Vector a -> Vector b -> Vector (a,b)+zipVector = zipVectorWith (,) -@disp = putStr . vecdisp ('dispf' 2)+-- | unzip for Vectors+unzipVector :: (Storable a, Storable b, Storable (a,b)) => Vector (a,b) -> (Vector a,Vector b)+unzipVector = unzipVectorWith id -\> disp ('linspace' 10 (0,1))-10 |> 0.00  0.11  0.22  0.33  0.44  0.56  0.67  0.78  0.89  1.00-@--}-vecdisp :: (Element t) => (Matrix t -> String) -> Vector t -> String-vecdisp f v-    = ((show (dim v) ++ " |> ") ++) . (++"\n")-    . unwords . lines .  tail . dropWhile (not . (`elem` " \n"))-    . f . trans . reshape 1-    $ v+-------------------------------------------------------------------

@@ -3,15 +3,13 @@ -- Module      :  Graphics.Plot -- Copyright   :  (c) Alberto Ruiz 2005-8 -- License     :  GPL-style--- +-- -- Maintainer  :  Alberto Ruiz (aruiz at um dot es) -- Stability   :  provisional -- Portability :  uses gnuplot and ImageMagick ----- Very basic (and provisional) drawing tools using gnuplot and imageMagick.--- --- This module is deprecated. It will be replaced by improved drawing tools based--- on the Gnuplot package by Henning Thielemann.+-- This module is deprecated. It can be replaced by improved drawing tools+-- available in the plot\\plot-gtk packages by Vivian McPhail or Gnuplot by Henning Thielemann. -----------------------------------------------------------------------------  module Graphics.Plot(@@ -20,7 +18,7 @@      plot, parametricPlot,  -    splot, mesh, mesh', meshdom,+    splot, mesh, meshdom,      matrixToPGM, imshow, @@ -28,40 +26,18 @@  ) where -import Data.Packed-import Numeric.LinearAlgebra(outer)+import Numeric.Container import Data.List(intersperse) import System.Process (system) -size = dim---- | Loads a real matrix from a formatted ASCII text file ---fromFile :: FilePath -> IO Matrix---fromFile filename = readFile filename >>= return . readMatrix read---- | Saves a real matrix to a formatted ascii text file-toFile' :: FilePath -> Matrix Double -> IO ()-toFile' filename matrix = writeFile filename (unlines . map unwords. map (map show) . toLists $ matrix)---------------------------------------------------------------------------- -- | From vectors x and y, it generates a pair of matrices to be used as x and y arguments for matrix functions. meshdom :: Vector Double -> Vector Double -> (Matrix Double , Matrix Double)-meshdom r1 r2 = (outer r1 (constant 1 (size r2)), outer (constant 1 (size r1)) r2)--gnuplotX :: String -> IO ()-gnuplotX command = do { _ <- system cmdstr; return()} where-    cmdstr = "echo \""++command++"\" | gnuplot -persist"--datafollows = "\\\"-\\\""--prep = (++"e\n\n") . unlines . map (unwords . (map show))+meshdom r1 r2 = (outer r1 (constant 1 (dim r2)), outer (constant 1 (dim r1)) r2)   {- | Draws a 3D surface representation of a real matrix. -> > mesh (hilb 20)+> > mesh $ build (10,10) (\\i j -> i + (j-5)^2)  In certain versions you can interactively rotate the graphic using the mouse. @@ -71,15 +47,6 @@     command = "splot "++datafollows++" matrix with lines\n"     dat = prep $ toLists $ m -mesh' :: Matrix Double -> IO ()-mesh' m = do-    writeFile "splot-gnu-command" "splot \"splot-tmp.txt\" matrix with lines; pause -1"; -    toFile' "splot-tmp.txt" m-    putStr "Press [Return] to close the graphic and continue... "-    _ <- system "gnuplot -persist splot-gnu-command"-    _ <- system "rm splot-tmp.txt splot-gnu-command"-    return ()- {- | Draws the surface represented by the function f in the desired ranges and number of points, internally using 'mesh'.  > > let f x y = cos (x + y) @@ -87,11 +54,15 @@  -} splot :: (Matrix Double->Matrix Double->Matrix Double) -> (Double,Double) -> (Double,Double) -> Int -> IO () -splot f rx ry n = mesh' z where+splot f rx ry n = mesh z where     (x,y) = meshdom (linspace n rx) (linspace n ry)     z = f x y -{- | plots several vectors against the first one -}+{- | plots several vectors against the first one ++> > let t = linspace 100 (-3,3) in mplot [t, sin t, exp (-t^2)]++-} mplot :: [Vector Double] -> IO () mplot m = gnuplotX (commands++dats) where     commands = if length m == 1 then command1 else commandmore@@ -103,26 +74,6 @@     dats = concat (replicate (length m-1) dat)  -{--mplot' m = do-    writeFile "plot-gnu-command" (commands++endcmd)-    toFile "plot-tmp.txt" (fromColumns m)-    putStr "Press [Return] to close the graphic and continue... "-    system "gnuplot plot-gnu-command"-    system "rm plot-tmp.txt plot-gnu-command"-    return ()- where-    commands = if length m == 1 then command1 else commandmore-    command1 = "plot \"plot-tmp.txt\" with lines\n"-    commandmore = "plot " ++ plots ++ "\n"-    plots = concat $ intersperse ", " (map cmd [2 .. length m])-    cmd k = "\"plot-tmp.txt\" using 1:"++show k++" with lines"-    endcmd = "pause -1"--}---- apply several functions to one object-mapf fs x = map ($ x) fs- {- | Draws a list of functions over a desired range and with a desired number of points   > > plot [sin, cos, sin.(3*)] (0,2*pi) 1000@@ -130,7 +81,8 @@ -} plot :: [Vector Double->Vector Double] -> (Double,Double) -> Int -> IO () plot fs rx n = mplot (x: mapf fs x)-    where x = linspace n rx  +    where x = linspace n rx+          mapf gs y = map ($ y) gs  {- | Draws a parametric curve. For instance, to draw a spiral we can do something like: @@ -150,12 +102,12 @@     r = rows m     header = "P2 "++show c++" "++show r++" "++show (round maxgray :: Int)++"\n"     maxgray = 255.0-    maxval = vectorMax $ flatten $ m-    minval = vectorMin $ flatten $ m-    scale = if (maxval == minval) +    maxval = maxElement m+    minval = minElement m+    scale' = if (maxval == minval)          then 0.0         else maxgray / (maxval - minval)-    f x = show ( round ( scale *(x - minval) ) :: Int )+    f x = show ( round ( scale' *(x - minval) ) :: Int )     ll = map (map f) (toLists m)  -- | imshow shows a representation of a matrix as a gray level image using ImageMagick's display.@@ -165,6 +117,15 @@     return ()  ----------------------------------------------------++gnuplotX :: String -> IO ()+gnuplotX command = do { _ <- system cmdstr; return()} where+    cmdstr = "echo \""++command++"\" | gnuplot -persist"++datafollows = "\\\"-\\\""++prep = (++"e\n\n") . unlines . map (unwords . (map show))+  gnuplotpdf :: String -> String -> [([[Double]], String)] -> IO () gnuplotpdf title command ds = gnuplot (prelude ++ command ++" "++ draw) >> postproc where

@@ -0,0 +1,140 @@+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.Chain+-- Copyright   :  (c) Vivian McPhail 2010+-- License     :  GPL-style+--+-- Maintainer  :  Vivian McPhail <haskell.vivian.mcphail <at> gmail.com>+-- Stability   :  provisional+-- Portability :  portable+--+-- optimisation of association order for chains of matrix multiplication+--+-----------------------------------------------------------------------------++module Numeric.Chain (+                      optimiseMult,+                     ) where++import Data.Maybe++import Data.Packed.Matrix+import Numeric.ContainerBoot++import qualified Data.Array.IArray as A++-----------------------------------------------------------------------------+{- | +     Provide optimal association order for a chain of matrix multiplications +     and apply the multiplications.++     The algorithm is the well-known O(n\^3) dynamic programming algorithm+     that builds a pyramid of optimal associations.++> m1, m2, m3, m4 :: Matrix Double+> m1 = (10><15) [1..]+> m2 = (15><20) [1..]+> m3 = (20><5) [1..]+> m4 = (5><10) [1..]++> >>> optimiseMult [m1,m2,m3,m4]++will perform @((m1 `multiply` (m2 `multiply` m3)) `multiply` m4)@++The naive left-to-right multiplication would take @4500@ scalar multiplications+whereas the optimised version performs @2750@ scalar multiplications.  The complexity+in this case is 32 (= 4^3/2) * (2 comparisons, 3 scalar multiplications, 3 scalar additions,+5 lookups, 2 updates) + a constant (= three table allocations)+-}+optimiseMult :: Product t => [Matrix t] -> Matrix t+optimiseMult = chain++-----------------------------------------------------------------------------++type Matrices a = A.Array Int (Matrix a)+type Sizes      = A.Array Int (Int,Int)+type Cost       = A.Array Int (A.Array Int (Maybe Int))+type Indexes    = A.Array Int (A.Array Int (Maybe ((Int,Int),(Int,Int))))++update :: A.Array Int (A.Array Int a) -> (Int,Int) -> a -> A.Array Int (A.Array Int a)+update a (r,c) e = a A.// [(r,(a A.! r) A.// [(c,e)])]++newWorkSpaceCost :: Int -> A.Array Int (A.Array Int (Maybe Int))+newWorkSpaceCost n = A.array (1,n) $ map (\i -> (i, subArray i)) [1..n]+   where subArray i = A.listArray (1,i) (repeat Nothing)++newWorkSpaceIndexes :: Int -> A.Array Int (A.Array Int (Maybe ((Int,Int),(Int,Int))))+newWorkSpaceIndexes n = A.array (1,n) $ map (\i -> (i, subArray i)) [1..n]+   where subArray i = A.listArray (1,i) (repeat Nothing)++matricesToSizes :: [Matrix a] -> Sizes+matricesToSizes ms = A.listArray (1,length ms) $ map (\m -> (rows m,cols m)) ms++chain :: Product a => [Matrix a] -> Matrix a+chain []  = error "chain: zero matrices to multiply"+chain [m] = m+chain [ml,mr] = ml `multiply` mr+chain ms = let ln = length ms+               ma = A.listArray (1,ln) ms+               mz = matricesToSizes ms+               i = chain_cost mz+           in chain_paren (ln,ln) i ma++chain_cost :: Sizes -> Indexes+chain_cost mz = let (_,u) = A.bounds mz+                    cost = newWorkSpaceCost u+                    ixes = newWorkSpaceIndexes u+                    (_,_,i) =  foldl chain_cost' (mz,cost,ixes) (order u)+                in i++chain_cost' :: (Sizes,Cost,Indexes) -> (Int,Int) -> (Sizes,Cost,Indexes)+chain_cost' sci@(mz,cost,ixes) (r,c) +    | c == 1                     = let cost' = update cost (r,c) (Just 0)+                                       ixes' = update ixes (r,c) (Just ((r,c),(r,c)))+                                       in (mz,cost',ixes')+    | otherwise                  = minimum_cost sci (r,c)++minimum_cost :: (Sizes,Cost,Indexes) -> (Int,Int) -> (Sizes,Cost,Indexes)+minimum_cost sci fu = foldl (smaller_cost fu) sci (fulcrum_order fu)++smaller_cost :: (Int,Int) -> (Sizes,Cost,Indexes) -> ((Int,Int),(Int,Int)) -> (Sizes,Cost,Indexes)+smaller_cost (r,c) (mz,cost,ixes) ix@((lr,lc),(rr,rc)) = let op_cost = (fromJust ((cost A.! lr) A.! lc))+                                                                       + (fromJust ((cost A.! rr) A.! rc))+                                                                       + ((fst $ mz A.! (lr-lc+1))+                                                                          *(snd $ mz A.! lc)+                                                                          *(snd $ mz A.! rr))+                                                             cost' = (cost A.! r) A.! c+                                                         in case cost' of+                                                                       Nothing -> let cost'' = update cost (r,c) (Just op_cost)+                                                                                      ixes'' = update ixes (r,c) (Just ix)+                                                                                  in (mz,cost'',ixes'')+                                                                       Just ct -> if op_cost < ct then+                                                                                  let cost'' = update cost (r,c) (Just op_cost)+                                                                                      ixes'' = update ixes (r,c) (Just ix)+                                                                                  in (mz,cost'',ixes'')+                                                                                  else (mz,cost,ixes)+                                                                         ++fulcrum_order (r,c) = let fs' = zip (repeat r) [1..(c-1)]+                      in map (partner (r,c)) fs'++partner (r,c) (a,b) = (((r-b),(c-b)),(a,b))++order 0 = []+order n = (order (n-1)) ++ (zip (repeat n) [1..n])++chain_paren :: Product a => (Int,Int) -> Indexes -> Matrices a -> Matrix a+chain_paren (r,c) ixes ma = let ((lr,lc),(rr,rc)) = fromJust $ (ixes A.! r) A.! c+                            in if lr == rr && lc == rc then (ma A.! lr)+                               else (chain_paren (lr,lc) ixes ma) `multiply` (chain_paren (rr,rc) ixes ma) ++--------------------------------------------------------------------------++{- TESTS -}++-- optimal association is ((m1*(m2*m3))*m4)+m1, m2, m3, m4 :: Matrix Double+m1 = (10><15) [1..]+m2 = (15><20) [1..]+m3 = (20><5) [1..]+m4 = (5><10) [1..]

@@ -0,0 +1,132 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE UndecidableInstances #-}++-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.Container+-- Copyright   :  (c) Alberto Ruiz 2010+-- License     :  GPL-style+--+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>+-- Stability   :  provisional+-- Portability :  portable+--+-- Basic numeric operations on 'Vector' and 'Matrix', including conversion routines.+--+-- The 'Container' class is used to define optimized generic functions which work+-- on 'Vector' and 'Matrix' with real or complex elements.+--+-- Some of these functions are also available in the instances of the standard+-- numeric Haskell classes provided by "Numeric.LinearAlgebra".+--+-----------------------------------------------------------------------------++module Numeric.Container (+    -- * Basic functions+    module Data.Packed,+    constant, linspace,+    diag, ident,+    ctrans,+    -- * Generic operations+    Container(..),+    -- * Matrix product+    Product(..),+    optimiseMult,+    mXm,mXv,vXm,(<.>),(<>),(<\>),+    outer, kronecker,+    -- * Random numbers+    RandDist(..),+    randomVector,+    gaussianSample,+    uniformSample,+    meanCov,+    -- * Element conversion+    Convert(..),+    Complexable(),+    RealElement(),++    RealOf, ComplexOf, SingleOf, DoubleOf,++    IndexOf,+    module Data.Complex,+    -- * Input / Output+    dispf, disps, dispcf, vecdisp, latexFormat, format,+    loadMatrix, saveMatrix, fromFile, fileDimensions,+    readMatrix,+    fscanfVector, fprintfVector, freadVector, fwriteVector,+    -- * Experimental+    build', konst',+    -- * Deprecated+    (.*),(*/),(<|>),(<->),+    vectorMax,vectorMin,+    vectorMaxIndex, vectorMinIndex+) where++import Data.Packed+import Data.Packed.Internal(constantD)+import Numeric.ContainerBoot+import Numeric.Chain+import Numeric.IO+import Data.Complex+import Numeric.LinearAlgebra.Algorithms(Field,linearSolveSVD)+import Data.Packed.Random++------------------------------------------------------------------++{- | 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)+constant = constantD-- about 2x faster++{- | Creates a real vector containing a range of values:++@\> linspace 5 (-3,7)+5 |> [-3.0,-0.5,2.0,4.5,7.0]@++Logarithmic spacing can be defined as follows:++@logspace n (a,b) = 10 ** linspace n (a,b)@+-}+linspace :: (Enum e, Container Vector e) => Int -> (e, e) -> Vector e+linspace n (a,b) = addConstant a $ scale s $ fromList [0 .. fromIntegral n-1]+    where s = (b-a)/fromIntegral (n-1)++-- | Dot product: @u \<.\> v = dot u v@+(<.>) :: Product t => Vector t -> Vector t -> t+infixl 7 <.>+(<.>) = dot++++--------------------------------------------------------++class Mul a b c | a b -> c where+ infixl 7 <>+ -- | Matrix-matrix, matrix-vector, and vector-matrix products.+ (<>)  :: Product t => a t -> b t -> c t++instance Mul Matrix Matrix Matrix where+    (<>) = mXm++instance Mul Matrix Vector Vector where+    (<>) m v = flatten $ m <> (asColumn v)++instance Mul Vector Matrix Vector where+    (<>) v m = flatten $ (asRow v) <> m++--------------------------------------------------------++-- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD).+(<\>) :: (Field a) => Matrix a -> Vector a -> Vector a+infixl 7 <\>+m <\> v = flatten (linearSolveSVD m (reshape 1 v))++--------------------------------------------------------

@@ -0,0 +1,583 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE UndecidableInstances #-}++-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.ContainerBoot+-- Copyright   :  (c) Alberto Ruiz 2010+-- License     :  GPL-style+--+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>+-- Stability   :  provisional+-- Portability :  portable+--+-- Module to avoid cyclyc dependencies.+--+-----------------------------------------------------------------------------++module Numeric.ContainerBoot (+    -- * Basic functions+    ident, diag, ctrans,+    -- * Generic operations+    Container(..),+    -- * Matrix product and related functions+    Product(..),+    mXm,mXv,vXm,+    outer, kronecker,+    -- * Element conversion+    Convert(..),+    Complexable(),+    RealElement(),++    RealOf, ComplexOf, SingleOf, DoubleOf,++    IndexOf,+    module Data.Complex,+    -- * Experimental+    build', konst',+    -- * Deprecated+    (.*),(*/),(<|>),(<->),+    vectorMax,vectorMin,+    vectorMaxIndex, vectorMinIndex+) where++import Data.Packed+import Numeric.Conversion+import Data.Packed.Internal+import Numeric.GSL.Vector++import Data.Complex+import Control.Monad(ap)++import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ)++import System.IO.Unsafe++-------------------------------------------------------------------++type family IndexOf c++type instance IndexOf Vector = Int+type instance IndexOf Matrix = (Int,Int)++type family ArgOf c a++type instance ArgOf Vector a = a -> a+type instance ArgOf Matrix a = a -> a -> a++-------------------------------------------------------------------++-- | Basic element-by-element functions for numeric containers+class (Complexable c, Fractional e, Element e) => Container c e where+    -- | create a structure with a single element+    scalar      :: e -> c e+    -- | complex conjugate+    conj        :: c e -> c e+    scale       :: e -> c e -> c e+    -- | scale the element by element reciprocal of the object:+    --+    -- @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@+    scaleRecip  :: e -> c e -> c e+    addConstant :: e -> c e -> c e+    add         :: c e -> c e -> c e+    sub         :: c e -> c e -> c e+    -- | element by element multiplication+    mul         :: c e -> c e -> c e+    -- | element by element division+    divide      :: c e -> c e -> c e+    equal       :: c e -> c e -> Bool+    --+    -- element by element inverse tangent+    arctan2     :: c e -> c e -> c e+    --+    -- | cannot implement instance Functor because of Element class constraint+    cmap        :: (Element a, Element b) => (a -> b) -> c a -> c b+    -- | constant structure of given size+    konst       :: e -> IndexOf c -> c e+    -- | create a structure using a function+    --+    -- Hilbert matrix of order N:+    --+    -- @hilb n = build (n,n) (\\i j -> 1/(i+j+1))@+    build       :: IndexOf c -> (ArgOf c e) -> c e+    --build       :: BoundsOf f -> f -> (ContainerOf f) e+    --+    -- | indexing function+    atIndex     :: c e -> IndexOf c -> e+    -- | index of min element+    minIndex    :: c e -> IndexOf c+    -- | index of max element+    maxIndex    :: c e -> IndexOf c+    -- | value of min element+    minElement  :: c e -> e+    -- | value of max element+    maxElement  :: c e -> e+    -- the C functions sumX/prodX are twice as fast as using foldVector+    -- | the sum of elements (faster than using @fold@)+    sumElements :: c e -> e+    -- | the product of elements (faster than using @fold@)+    prodElements :: c e -> e++--------------------------------------------------------------------------++instance Container Vector Float where+    scale = vectorMapValF Scale+    scaleRecip = vectorMapValF Recip+    addConstant = vectorMapValF AddConstant+    add = vectorZipF Add+    sub = vectorZipF Sub+    mul = vectorZipF Mul+    divide = vectorZipF Div+    equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0+    arctan2 = vectorZipF ATan2+    scalar x = fromList [x]+    konst = constantD+    build = buildV+    conj = id+    cmap = mapVector+    atIndex = (@>)+    minIndex     = round . toScalarF MinIdx+    maxIndex     = round . toScalarF MaxIdx+    minElement  = toScalarF Min+    maxElement  = toScalarF Max+    sumElements  = sumF+    prodElements = prodF++instance Container Vector Double where+    scale = vectorMapValR Scale+    scaleRecip = vectorMapValR Recip+    addConstant = vectorMapValR AddConstant+    add = vectorZipR Add+    sub = vectorZipR Sub+    mul = vectorZipR Mul+    divide = vectorZipR Div+    equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0+    arctan2 = vectorZipR ATan2+    scalar x = fromList [x]+    konst = constantD+    build = buildV+    conj = id+    cmap = mapVector+    atIndex = (@>)+    minIndex     = round . toScalarR MinIdx+    maxIndex     = round . toScalarR MaxIdx+    minElement  = toScalarR Min+    maxElement  = toScalarR Max+    sumElements  = sumR+    prodElements = prodR++instance Container Vector (Complex Double) where+    scale = vectorMapValC Scale+    scaleRecip = vectorMapValC Recip+    addConstant = vectorMapValC AddConstant+    add = vectorZipC Add+    sub = vectorZipC Sub+    mul = vectorZipC Mul+    divide = vectorZipC Div+    equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0+    arctan2 = vectorZipC ATan2+    scalar x = fromList [x]+    konst = constantD+    build = buildV+    conj = conjugateC+    cmap = mapVector+    atIndex = (@>)+    minIndex     = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)+    maxIndex     = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)+    minElement  = ap (@>) minIndex+    maxElement  = ap (@>) maxIndex+    sumElements  = sumC+    prodElements = prodC++instance Container Vector (Complex Float) where+    scale = vectorMapValQ Scale+    scaleRecip = vectorMapValQ Recip+    addConstant = vectorMapValQ AddConstant+    add = vectorZipQ Add+    sub = vectorZipQ Sub+    mul = vectorZipQ Mul+    divide = vectorZipQ Div+    equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0+    arctan2 = vectorZipQ ATan2+    scalar x = fromList [x]+    konst = constantD+    build = buildV+    conj = conjugateQ+    cmap = mapVector+    atIndex = (@>)+    minIndex     = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)+    maxIndex     = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)+    minElement  = ap (@>) minIndex+    maxElement  = ap (@>) maxIndex+    sumElements  = sumQ+    prodElements = prodQ++---------------------------------------------------------------++instance (Container Vector a) => Container Matrix a where+    scale x = liftMatrix (scale x)+    scaleRecip x = liftMatrix (scaleRecip x)+    addConstant x = liftMatrix (addConstant x)+    add = liftMatrix2 add+    sub = liftMatrix2 sub+    mul = liftMatrix2 mul+    divide = liftMatrix2 divide+    equal a b = cols a == cols b && flatten a `equal` flatten b+    arctan2 = liftMatrix2 arctan2+    scalar x = (1><1) [x]+    konst v (r,c) = reshape c (konst v (r*c))+    build = buildM+    conj = liftMatrix conj+    cmap f = liftMatrix (mapVector f)+    atIndex = (@@>)+    minIndex m = let (r,c) = (rows m,cols m)+                     i = (minIndex $ flatten m)+                 in (i `div` c,i `mod` c)+    maxIndex m = let (r,c) = (rows m,cols m)+                     i = (maxIndex $ flatten m)+                 in (i `div` c,i `mod` c)+    minElement = ap (@@>) minIndex+    maxElement = ap (@@>) maxIndex+    sumElements = sumElements . flatten+    prodElements = prodElements . flatten++----------------------------------------------------++-- | Matrix product and related functions+class Element e => Product e where+    -- | matrix product+    multiply :: Matrix e -> Matrix e -> Matrix e+    -- | dot (inner) product+    dot        :: Vector e -> Vector e -> e+    -- | sum of absolute value of elements (differs in complex case from @norm1@)+    absSum     :: Vector e -> RealOf e+    -- | sum of absolute value of elements+    norm1      :: Vector e -> RealOf e+    -- | euclidean norm+    norm2      :: Vector e -> RealOf e+    -- | element of maximum magnitude+    normInf    :: Vector e -> RealOf e++instance Product Float where+    norm2      = toScalarF Norm2+    absSum     = toScalarF AbsSum+    dot        = dotF+    norm1      = toScalarF AbsSum+    normInf    = maxElement . vectorMapF Abs+    multiply = multiplyF++instance Product Double where+    norm2      = toScalarR Norm2+    absSum     = toScalarR AbsSum+    dot        = dotR+    norm1      = toScalarR AbsSum+    normInf    = maxElement . vectorMapR Abs+    multiply = multiplyR++instance Product (Complex Float) where+    norm2      = toScalarQ Norm2+    absSum     = toScalarQ AbsSum+    dot        = dotQ+    norm1      = sumElements . fst . fromComplex . vectorMapQ Abs+    normInf    = maxElement . fst . fromComplex . vectorMapQ Abs+    multiply = multiplyQ++instance Product (Complex Double) where+    norm2      = toScalarC Norm2+    absSum     = toScalarC AbsSum+    dot        = dotC+    norm1      = sumElements . fst . fromComplex . vectorMapC Abs+    normInf    = maxElement . fst . fromComplex . vectorMapC Abs+    multiply = multiplyC++----------------------------------------------------------++-- synonym for matrix product+mXm :: Product t => Matrix t -> Matrix t -> Matrix t+mXm = multiply++-- matrix - vector product+mXv :: Product t => Matrix t -> Vector t -> Vector t+mXv m v = flatten $ m `mXm` (asColumn v)++-- vector - matrix product+vXm :: Product t => Vector t -> Matrix t -> Vector t+vXm v m = flatten $ (asRow v) `mXm` m++{- | Outer product of two vectors.++@\> 'fromList' [1,2,3] \`outer\` 'fromList' [5,2,3]+(3><3)+ [  5.0, 2.0, 3.0+ , 10.0, 4.0, 6.0+ , 15.0, 6.0, 9.0 ]@+-}+outer :: (Product t) => Vector t -> Vector t -> Matrix t+outer u v = asColumn u `multiply` asRow v++{- | Kronecker product of two matrices.++@m1=(2><3)+ [ 1.0,  2.0, 0.0+ , 0.0, -1.0, 3.0 ]+m2=(4><3)+ [  1.0,  2.0,  3.0+ ,  4.0,  5.0,  6.0+ ,  7.0,  8.0,  9.0+ , 10.0, 11.0, 12.0 ]@++@\> kronecker m1 m2+(8><9)+ [  1.0,  2.0,  3.0,   2.0,   4.0,   6.0,  0.0,  0.0,  0.0+ ,  4.0,  5.0,  6.0,   8.0,  10.0,  12.0,  0.0,  0.0,  0.0+ ,  7.0,  8.0,  9.0,  14.0,  16.0,  18.0,  0.0,  0.0,  0.0+ , 10.0, 11.0, 12.0,  20.0,  22.0,  24.0,  0.0,  0.0,  0.0+ ,  0.0,  0.0,  0.0,  -1.0,  -2.0,  -3.0,  3.0,  6.0,  9.0+ ,  0.0,  0.0,  0.0,  -4.0,  -5.0,  -6.0, 12.0, 15.0, 18.0+ ,  0.0,  0.0,  0.0,  -7.0,  -8.0,  -9.0, 21.0, 24.0, 27.0+ ,  0.0,  0.0,  0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@+-}+kronecker :: (Product t) => Matrix t -> Matrix t -> Matrix t+kronecker a b = fromBlocks+              . splitEvery (cols a)+              . map (reshape (cols b))+              . toRows+              $ flatten a `outer` flatten b++-------------------------------------------------------------------+++class Convert t where+    real    :: Container c t => c (RealOf t) -> c t+    complex :: Container c t => c t -> c (ComplexOf t)+    single  :: Container c t => c t -> c (SingleOf t)+    double  :: Container c t => c t -> c (DoubleOf t)+    toComplex   :: (Container c t, RealElement t) => (c t, c t) -> c (Complex t)+    fromComplex :: (Container c t, RealElement t) => c (Complex t) -> (c t, c t)+++instance Convert Double where+    real = id+    complex = comp'+    single = single'+    double = id+    toComplex = toComplex'+    fromComplex = fromComplex'++instance Convert Float where+    real = id+    complex = comp'+    single = id+    double = double'+    toComplex = toComplex'+    fromComplex = fromComplex'++instance Convert (Complex Double) where+    real = comp'+    complex = id+    single = single'+    double = id+    toComplex = toComplex'+    fromComplex = fromComplex'++instance Convert (Complex Float) where+    real = comp'+    complex = id+    single = id+    double = double'+    toComplex = toComplex'+    fromComplex = fromComplex'++-------------------------------------------------------------------++type family RealOf x++type instance RealOf Double = Double+type instance RealOf (Complex Double) = Double++type instance RealOf Float = Float+type instance RealOf (Complex Float) = Float++type family ComplexOf x++type instance ComplexOf Double = Complex Double+type instance ComplexOf (Complex Double) = Complex Double++type instance ComplexOf Float = Complex Float+type instance ComplexOf (Complex Float) = Complex Float++type family SingleOf x++type instance SingleOf Double = Float+type instance SingleOf Float  = Float++type instance SingleOf (Complex a) = Complex (SingleOf a)++type family DoubleOf x++type instance DoubleOf Double = Double+type instance DoubleOf Float  = Double++type instance DoubleOf (Complex a) = Complex (DoubleOf a)++type family ElementOf c++type instance ElementOf (Vector a) = a+type instance ElementOf (Matrix a) = a++------------------------------------------------------------++conjugateAux fun x = unsafePerformIO $ do+    v <- createVector (dim x)+    app2 fun vec x vec v "conjugateAux"+    return v++conjugateQ :: Vector (Complex Float) -> Vector (Complex Float)+conjugateQ = conjugateAux c_conjugateQ+foreign import ccall "conjugateQ" c_conjugateQ :: TQVQV++conjugateC :: Vector (Complex Double) -> Vector (Complex Double)+conjugateC = conjugateAux c_conjugateC+foreign import ccall "conjugateC" c_conjugateC :: TCVCV++----------------------------------------------------++{-# DEPRECATED (.*) "use scale a x or scalar a * x" #-}++-- -- | @x .* a = scale x a@+-- (.*) :: (Linear c a) => a -> c a -> c a+infixl 7 .*+a .* x = scale a x++----------------------------------------------------++{-# DEPRECATED (*/) "use scale (recip a) x or x / scalar a" #-}++-- -- | @a *\/ x = scale (recip x) a@+-- (*/) :: (Linear c a) => c a -> a -> c a+infixl 7 */+v */ x = scale (recip x) v+++------------------------------------------------++{-# DEPRECATED (<|>) "define operator a & b = fromBlocks[[a,b]] and use asRow/asColumn to join vectors" #-}+{-# DEPRECATED (<->) "define operator a // b = fromBlocks[[a],[b]] and use asRow/asColumn to join vectors" #-}++class Joinable a b where+    joinH :: Element t => a t -> b t -> Matrix t+    joinV :: Element t => a t -> b t -> Matrix t++instance Joinable Matrix Matrix where+    joinH m1 m2 = fromBlocks [[m1,m2]]+    joinV m1 m2 = fromBlocks [[m1],[m2]]++instance Joinable Matrix Vector where+    joinH m v = joinH m (asColumn v)+    joinV m v = joinV m (asRow v)++instance Joinable Vector Matrix where+    joinH v m = joinH (asColumn v) m+    joinV v m = joinV (asRow v) m++infixl 4 <|>+infixl 3 <->++{-- - | Horizontal concatenation of matrices and vectors:++@> (ident 3 \<-\> 3 * ident 3) \<|\> fromList [1..6.0]+(6><4)+ [ 1.0, 0.0, 0.0, 1.0+ , 0.0, 1.0, 0.0, 2.0+ , 0.0, 0.0, 1.0, 3.0+ , 3.0, 0.0, 0.0, 4.0+ , 0.0, 3.0, 0.0, 5.0+ , 0.0, 0.0, 3.0, 6.0 ]@+-}+-- (<|>) :: (Element t, Joinable a b) => a t -> b t -> Matrix t+a <|> b = joinH a b++-- -- | Vertical concatenation of matrices and vectors.+-- (<->) :: (Element t, Joinable a b) => a t -> b t -> Matrix t+a <-> b = joinV a b++-------------------------------------------------------------------++{-# DEPRECATED vectorMin "use minElement" #-}+vectorMin :: (Container Vector t, Element t) => Vector t -> t+vectorMin = minElement++{-# DEPRECATED vectorMax "use maxElement" #-}+vectorMax :: (Container Vector t, Element t) => Vector t -> t+vectorMax = maxElement+++{-# DEPRECATED vectorMaxIndex "use minIndex" #-}+vectorMaxIndex :: Vector Double -> Int+vectorMaxIndex = round . toScalarR MaxIdx++{-# DEPRECATED vectorMinIndex "use maxIndex" #-}+vectorMinIndex :: Vector Double -> Int+vectorMinIndex = round . toScalarR MinIdx++-----------------------------------------------------++class Build f where+    build' :: BoundsOf f -> f -> ContainerOf f++type family BoundsOf x++type instance BoundsOf (a->a) = Int+type instance BoundsOf (a->a->a) = (Int,Int)++type family ContainerOf x++type instance ContainerOf (a->a) = Vector a+type instance ContainerOf (a->a->a) = Matrix a++instance (Element a, Num a) => Build (a->a) where+    build' = buildV++instance (Element a, Num a) => Build (a->a->a) where+    build' = buildM++buildM (rc,cc) f = fromLists [ [f r c | c <- cs] | r <- rs ]+    where rs = map fromIntegral [0 .. (rc-1)]+          cs = map fromIntegral [0 .. (cc-1)]++buildV n f = fromList [f k | k <- ks]+    where ks = map fromIntegral [0 .. (n-1)]++----------------------------------------------------+-- experimental++class Konst s where+    konst' :: Element e => e -> s -> ContainerOf' s e++type family ContainerOf' x y++type instance ContainerOf' Int a = Vector a+type instance ContainerOf' (Int,Int) a = Matrix a++instance Konst Int where+    konst' = constantD++instance Konst (Int,Int) where+    konst' k (r,c) = reshape c $ konst' k (r*c)++--------------------------------------------------------+-- | conjugate transpose+ctrans :: (Container Vector e, Element e) => Matrix e -> Matrix e+ctrans = liftMatrix conj . trans++-- | Creates a square matrix with a given diagonal.+diag :: (Num a, Element a) => Vector a -> Matrix a+diag v = diagRect 0 v n n where n = dim v++-- | creates the identity matrix of given dimension+ident :: (Num a, Element a) => Int -> Matrix a+ident n = diag (constantD 1 n)

@@ -0,0 +1,91 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE UndecidableInstances #-}++-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.Conversion+-- Copyright   :  (c) Alberto Ruiz 2010+-- License     :  GPL-style+--+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>+-- Stability   :  provisional+-- Portability :  portable+--+-- Conversion routines+--+-----------------------------------------------------------------------------++module Numeric.Conversion (+    Complexable(..), RealElement,+    module Data.Complex+) where++import Data.Packed.Internal.Vector+import Data.Packed.Internal.Matrix+import Data.Complex+import Control.Arrow((***))++-------------------------------------------------------------------++-- | Supported single-double precision type pairs+class (Element s, Element d) => Precision s d | s -> d, d -> s where+    double2FloatG :: Vector d -> Vector s+    float2DoubleG :: Vector s -> Vector d++instance Precision Float Double where+    double2FloatG = double2FloatV+    float2DoubleG = float2DoubleV++instance Precision (Complex Float) (Complex Double) where+    double2FloatG = asComplex . double2FloatV . asReal+    float2DoubleG = asComplex . float2DoubleV . asReal++-- | Supported real types+class (Element t, Element (Complex t), RealFloat t+--       , RealOf t ~ t, RealOf (Complex t) ~ t+       )+    => RealElement t++instance RealElement Double+instance RealElement Float+++-- | Structures that may contain complex numbers+class Complexable c where+    toComplex'   :: (RealElement e) => (c e, c e) -> c (Complex e)+    fromComplex' :: (RealElement e) => c (Complex e) -> (c e, c e)+    comp'        :: (RealElement e) => c e -> c (Complex e)+    single'      :: Precision a b => c b -> c a+    double'      :: Precision a b => c a -> c b+++instance Complexable Vector where+    toComplex' = toComplexV+    fromComplex' = fromComplexV+    comp' v = toComplex' (v,constantD 0 (dim v))+    single' = double2FloatG+    double' = float2DoubleG+++-- | creates a complex vector from vectors with real and imaginary parts+toComplexV :: (RealElement a) => (Vector a, Vector a) ->  Vector (Complex a)+toComplexV (r,i) = asComplex $ flatten $ fromColumns [r,i]++-- | the inverse of 'toComplex'+fromComplexV :: (RealElement a) => Vector (Complex a) -> (Vector a, Vector a)+fromComplexV z = (r,i) where+    [r,i] = toColumns $ reshape 2 $ asReal z+++instance Complexable Matrix where+    toComplex' = uncurry $ liftMatrix2 $ curry toComplex'+    fromComplex' z = (reshape c *** reshape c) . fromComplex' . flatten $ z+        where c = cols z+    comp' = liftMatrix comp'+    single' = liftMatrix single'+    double' = liftMatrix double'+

@@ -109,7 +109,7 @@     sol = toList vsol     c = max 1 (chi/sqrt (fromIntegral dof))     dof = length dat - (rows cov)-    chi = pnorm PNorm2 (fromList $ cost (resMs model) dat sol)+    chi = norm2 (fromList $ cost (resMs model) dat sol)     js = fromLists $ jacobian (resDs deriv) dat sol     cov = inv $ trans js <> js     errs = toList $ scalar c * sqrt (takeDiag cov)

@@ -14,10 +14,13 @@ -----------------------------------------------------------------------------  module Numeric.GSL.Vector (-    FunCodeS(..), toScalarR,-    FunCodeV(..), vectorMapR, vectorMapC,-    FunCodeSV(..), vectorMapValR, vectorMapValC,-    FunCodeVV(..), vectorZipR, vectorZipC,+    sumF, sumR, sumQ, sumC,+    prodF, prodR, prodQ, prodC,+    dotF, dotR, dotQ, dotC,+    FunCodeS(..), toScalarR, toScalarF, toScalarC, toScalarQ,+    FunCodeV(..), vectorMapR, vectorMapC, vectorMapF, vectorMapQ,+    FunCodeSV(..), vectorMapValR, vectorMapValC, vectorMapValF, vectorMapValQ,+    FunCodeVV(..), vectorZipR, vectorZipC, vectorZipF, vectorZipQ,     RandDist(..), randomVector ) where @@ -76,6 +79,107 @@  ------------------------------------------------------------------ +-- | sum of elements+sumF :: Vector Float -> Float+sumF x = unsafePerformIO $ do+           r <- createVector 1+           app2 c_sumF vec x vec r "sumF"+           return $ r @> 0++-- | sum of elements+sumR :: Vector Double -> Double+sumR x = unsafePerformIO $ do+           r <- createVector 1+           app2 c_sumR vec x vec r "sumR"+           return $ r @> 0++-- | sum of elements+sumQ :: Vector (Complex Float) -> Complex Float+sumQ x = unsafePerformIO $ do+           r <- createVector 1+           app2 c_sumQ vec x vec r "sumQ"+           return $ r @> 0++-- | sum of elements+sumC :: Vector (Complex Double) -> Complex Double+sumC x = unsafePerformIO $ do+           r <- createVector 1+           app2 c_sumC vec x vec r "sumC"+           return $ r @> 0++foreign import ccall safe "gsl-aux.h sumF" c_sumF :: TFF+foreign import ccall safe "gsl-aux.h sumR" c_sumR :: TVV+foreign import ccall safe "gsl-aux.h sumQ" c_sumQ :: TQVQV+foreign import ccall safe "gsl-aux.h sumC" c_sumC :: TCVCV++-- | product of elements+prodF :: Vector Float -> Float+prodF x = unsafePerformIO $ do+           r <- createVector 1+           app2 c_prodF vec x vec r "prodF"+           return $ r @> 0++-- | product of elements+prodR :: Vector Double -> Double+prodR x = unsafePerformIO $ do+           r <- createVector 1+           app2 c_prodR vec x vec r "prodR"+           return $ r @> 0++-- | product of elements+prodQ :: Vector (Complex Float) -> Complex Float+prodQ x = unsafePerformIO $ do+           r <- createVector 1+           app2 c_prodQ vec x vec r "prodQ"+           return $ r @> 0++-- | product of elements+prodC :: Vector (Complex Double) -> Complex Double+prodC x = unsafePerformIO $ do+           r <- createVector 1+           app2 c_prodC vec x vec r "prodC"+           return $ r @> 0++foreign import ccall safe "gsl-aux.h prodF" c_prodF :: TFF+foreign import ccall safe "gsl-aux.h prodR" c_prodR :: TVV+foreign import ccall safe "gsl-aux.h prodQ" c_prodQ :: TQVQV+foreign import ccall safe "gsl-aux.h prodC" c_prodC :: TCVCV++-- | dot product+dotF :: Vector Float -> Vector Float -> Float+dotF x y = unsafePerformIO $ do+           r <- createVector 1+           app3 c_dotF vec x vec y vec r "dotF"+           return $ r @> 0++-- | dot product+dotR :: Vector Double -> Vector Double -> Double+dotR x y = unsafePerformIO $ do+           r <- createVector 1+           app3 c_dotR vec x vec y vec r "dotR"+           return $ r @> 0++-- | dot product+dotQ :: Vector (Complex Float) -> Vector (Complex Float) -> Complex Float+dotQ x y = unsafePerformIO $ do+           r <- createVector 1+           app3 c_dotQ vec x vec y vec r "dotQ"+           return $ r @> 0++-- | dot product+dotC :: Vector (Complex Double) -> Vector (Complex Double) -> Complex Double+dotC x y = unsafePerformIO $ do+           r <- createVector 1+           app3 c_dotC vec x vec y vec r "dotC"+           return $ r @> 0++foreign import ccall safe "gsl-aux.h dotF" c_dotF :: TFFF+foreign import ccall safe "gsl-aux.h dotR" c_dotR :: TVVV+foreign import ccall safe "gsl-aux.h dotQ" c_dotQ :: TQVQVQV+foreign import ccall safe "gsl-aux.h dotC" c_dotC :: TCVCVCV++------------------------------------------------------------------+ toScalarAux fun code v = unsafePerformIO $ do     r <- createVector 1     app2 (fun (fromei code)) vec v vec r "toScalarAux"@@ -106,6 +210,24 @@  foreign import ccall safe "gsl-aux.h toScalarR" c_toScalarR :: CInt -> TVV +-- | obtains different functions of a vector: norm1, norm2, max, min, posmax, posmin, etc.+toScalarF :: FunCodeS -> Vector Float -> Float+toScalarF oper =  toScalarAux c_toScalarF (fromei oper)++foreign import ccall safe "gsl-aux.h toScalarF" c_toScalarF :: CInt -> TFF++-- | obtains different functions of a vector: only norm1, norm2+toScalarC :: FunCodeS -> Vector (Complex Double) -> Double+toScalarC oper =  toScalarAux c_toScalarC (fromei oper)++foreign import ccall safe "gsl-aux.h toScalarC" c_toScalarC :: CInt -> TCVV++-- | obtains different functions of a vector: only norm1, norm2+toScalarQ :: FunCodeS -> Vector (Complex Float) -> Float+toScalarQ oper =  toScalarAux c_toScalarQ (fromei oper)++foreign import ccall safe "gsl-aux.h toScalarQ" c_toScalarQ :: CInt -> TQVF+ ------------------------------------------------------------------  -- | map of real vectors with given function@@ -120,6 +242,18 @@  foreign import ccall safe "gsl-aux.h mapC" c_vectorMapC :: CInt -> TCVCV +-- | map of real vectors with given function+vectorMapF :: FunCodeV -> Vector Float -> Vector Float+vectorMapF = vectorMapAux c_vectorMapF++foreign import ccall safe "gsl-aux.h mapF" c_vectorMapF :: CInt -> TFF++-- | map of real vectors with given function+vectorMapQ :: FunCodeV -> Vector (Complex Float) -> Vector (Complex Float)+vectorMapQ = vectorMapAux c_vectorMapQ++foreign import ccall safe "gsl-aux.h mapQ" c_vectorMapQ :: CInt -> TQVQV+ -------------------------------------------------------------------  -- | map of real vectors with given function@@ -134,6 +268,18 @@  foreign import ccall safe "gsl-aux.h mapValC" c_vectorMapValC :: CInt -> Ptr (Complex Double) -> TCVCV +-- | map of real vectors with given function+vectorMapValF :: FunCodeSV -> Float -> Vector Float -> Vector Float+vectorMapValF oper = vectorMapValAux c_vectorMapValF (fromei oper)++foreign import ccall safe "gsl-aux.h mapValF" c_vectorMapValF :: CInt -> Ptr Float -> TFF++-- | map of complex vectors with given function+vectorMapValQ :: FunCodeSV -> Complex Float -> Vector (Complex Float) -> Vector (Complex Float)+vectorMapValQ oper = vectorMapValAux c_vectorMapValQ (fromei oper)++foreign import ccall safe "gsl-aux.h mapValQ" c_vectorMapValQ :: CInt -> Ptr (Complex Float) -> TQVQV+ -------------------------------------------------------------------  -- | elementwise operation on real vectors@@ -147,6 +293,18 @@ vectorZipC = vectorZipAux c_vectorZipC  foreign import ccall safe "gsl-aux.h zipC" c_vectorZipC :: CInt -> TCVCVCV++-- | elementwise operation on real vectors+vectorZipF :: FunCodeVV -> Vector Float -> Vector Float -> Vector Float+vectorZipF = vectorZipAux c_vectorZipF++foreign import ccall safe "gsl-aux.h zipF" c_vectorZipF :: CInt -> TFFF++-- | elementwise operation on complex vectors+vectorZipQ :: FunCodeVV -> Vector (Complex Float) -> Vector (Complex Float) -> Vector (Complex Float)+vectorZipQ = vectorZipAux c_vectorZipQ++foreign import ccall safe "gsl-aux.h zipQ" c_vectorZipQ :: CInt -> TQVQVQV  ----------------------------------------------------------------------- 

@@ -10,6 +10,16 @@ #define KCVEC(A) int A##n, const gsl_complex*A##p #define KCMAT(A) int A##r, int A##c, const gsl_complex* A##p +#define FVEC(A) int A##n, float*A##p+#define FMAT(A) int A##r, int A##c, float* A##p+#define KFVEC(A) int A##n, const float*A##p+#define KFMAT(A) int A##r, int A##c, const float* A##p++#define QVEC(A) int A##n, gsl_complex_float*A##p+#define QMAT(A) int A##r, int A##c, gsl_complex_float* A##p+#define KQVEC(A) int A##n, const gsl_complex_float*A##p+#define KQMAT(A) int A##r, int A##c, const gsl_complex_float* A##p+ #include <gsl/gsl_blas.h> #include <gsl/gsl_math.h> #include <gsl/gsl_errno.h>@@ -64,12 +74,24 @@ #define KCVVIEW(A) gsl_vector_complex_const_view A = gsl_vector_complex_const_view_array((double*)A##p,A##n) #define KCMVIEW(A) gsl_matrix_complex_const_view A = gsl_matrix_complex_const_view_array((double*)A##p,A##r,A##c) +#define FVVIEW(A) gsl_vector_float_view A = gsl_vector_float_view_array(A##p,A##n)+#define FMVIEW(A) gsl_matrix_float_view A = gsl_matrix_float_view_array(A##p,A##r,A##c)+#define QVVIEW(A) gsl_vector_complex_float_view A = gsl_vector_float_complex_view_array((float*)A##p,A##n)+#define QMVIEW(A) gsl_matrix_complex_float_view A = gsl_matrix_float_complex_view_array((float*)A##p,A##r,A##c)+#define KFVVIEW(A) gsl_vector_float_const_view A = gsl_vector_float_const_view_array(A##p,A##n)+#define KFMVIEW(A) gsl_matrix_float_const_view A = gsl_matrix_float_const_view_array(A##p,A##r,A##c)+#define KQVVIEW(A) gsl_vector_complex_float_const_view A = gsl_vector_complex_float_const_view_array((float*)A##p,A##n)+#define KQMVIEW(A) gsl_matrix_complex_float_const_view A = gsl_matrix_complex_float_const_view_array((float*)A##p,A##r,A##c)+ #define V(a) (&a.vector) #define M(a) (&a.matrix)  #define GCVEC(A) int A##n, gsl_complex*A##p #define KGCVEC(A) int A##n, const gsl_complex*A##p +#define GQVEC(A) int A##n, gsl_complex_float*A##p+#define KGQVEC(A) int A##n, const gsl_complex_float*A##p+ #define BAD_SIZE 2000 #define BAD_CODE 2001 #define MEM      2002@@ -81,6 +103,154 @@ }  +int sumF(KFVEC(x),FVEC(r)) {+    DEBUGMSG("sumF");+    REQUIRES(rn==1,BAD_SIZE);+    int i;+    float res = 0;+    for (i = 0; i < xn; i++) res += xp[i];+    rp[0] = res;+    OK+}+    +int sumR(KRVEC(x),RVEC(r)) {+    DEBUGMSG("sumR");+    REQUIRES(rn==1,BAD_SIZE);+    int i;+    double res = 0;+    for (i = 0; i < xn; i++) res += xp[i];+    rp[0] = res;+    OK+}+    +int sumQ(KQVEC(x),QVEC(r)) {+    DEBUGMSG("sumQ");+    REQUIRES(rn==1,BAD_SIZE);+    int i;+    gsl_complex_float res;+    res.dat[0] = 0;+    res.dat[1] = 0;+    for (i = 0; i < xn; i++) {+      res.dat[0] += xp[i].dat[0];+      res.dat[1] += xp[i].dat[1];+    }+    rp[0] = res;+    OK+}+    +int sumC(KCVEC(x),CVEC(r)) {+    DEBUGMSG("sumC");+    REQUIRES(rn==1,BAD_SIZE);+    int i;+    gsl_complex res;+    res.dat[0] = 0;+    res.dat[1] = 0;+    for (i = 0; i < xn; i++)  {+      res.dat[0] += xp[i].dat[0];+      res.dat[1] += xp[i].dat[1];+    }+    rp[0] = res;+    OK+}++int prodF(KFVEC(x),FVEC(r)) {+    DEBUGMSG("prodF");+    REQUIRES(rn==1,BAD_SIZE);+    int i;+    float res = 1;+    for (i = 0; i < xn; i++) res *= xp[i];+    rp[0] = res;+    OK+}+    +int prodR(KRVEC(x),RVEC(r)) {+    DEBUGMSG("prodR");+    REQUIRES(rn==1,BAD_SIZE);+    int i;+    double res = 1;+    for (i = 0; i < xn; i++) res *= xp[i];+    rp[0] = res;+    OK+}+    +int prodQ(KQVEC(x),QVEC(r)) {+    DEBUGMSG("prodQ");+    REQUIRES(rn==1,BAD_SIZE);+    int i;+    gsl_complex_float res;+    float temp;+    res.dat[0] = 1;+    res.dat[1] = 0;+    for (i = 0; i < xn; i++) {+      temp       = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1];+      res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0];+      res.dat[0] = temp;+    }+    rp[0] = res;+    OK+}+    +int prodC(KCVEC(x),CVEC(r)) {+    DEBUGMSG("prodC");+    REQUIRES(rn==1,BAD_SIZE);+    int i;+    gsl_complex res;+    double temp;+    res.dat[0] = 1;+    res.dat[1] = 0;+    for (i = 0; i < xn; i++)  {+      temp       = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1];+      res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0];+      res.dat[0] = temp;+    }+    rp[0] = res;+    OK+}++int dotF(KFVEC(x), KFVEC(y), FVEC(r)) {+    DEBUGMSG("dotF");+    REQUIRES(xn==yn,BAD_SIZE); +    REQUIRES(rn==1,BAD_SIZE);+    DEBUGMSG("dotF");+    KFVVIEW(x);+    KFVVIEW(y);+    gsl_blas_sdot(V(x),V(y),rp);+    OK+}+    +int dotR(KRVEC(x), KRVEC(y), RVEC(r)) {+    DEBUGMSG("dotR");+    REQUIRES(xn==yn,BAD_SIZE); +    REQUIRES(rn==1,BAD_SIZE);+    DEBUGMSG("dotR");+    KDVVIEW(x);+    KDVVIEW(y);+    gsl_blas_ddot(V(x),V(y),rp);+    OK+}+    +int dotQ(KQVEC(x), KQVEC(y), QVEC(r)) {+    DEBUGMSG("dotQ");+    REQUIRES(xn==yn,BAD_SIZE); +    REQUIRES(rn==1,BAD_SIZE);+    DEBUGMSG("dotQ");+    KQVVIEW(x);+    KQVVIEW(y);+    gsl_blas_cdotu(V(x),V(y),rp);+    OK+}+    +int dotC(KCVEC(x), KCVEC(y), CVEC(r)) {+    DEBUGMSG("dotC");+    REQUIRES(xn==yn,BAD_SIZE); +    REQUIRES(rn==1,BAD_SIZE);+    DEBUGMSG("dotC");+    KCVVIEW(x);+    KCVVIEW(y);+    gsl_blas_zdotu(V(x),V(y),rp);+    OK+}+     int toScalarR(int code, KRVEC(x), RVEC(r)) {      REQUIRES(rn==1,BAD_SIZE);     DEBUGMSG("toScalarR");@@ -99,7 +269,54 @@     OK } +int toScalarF(int code, KFVEC(x), FVEC(r)) { +    REQUIRES(rn==1,BAD_SIZE);+    DEBUGMSG("toScalarF");+    KFVVIEW(x);+    float res;+    switch(code) {+        case 0: { res = gsl_blas_snrm2(V(x)); break; } +        case 1: { res = gsl_blas_sasum(V(x));  break; }+        case 2: { res = gsl_vector_float_max_index(V(x));  break; }+        case 3: { res = gsl_vector_float_max(V(x));  break; }+        case 4: { res = gsl_vector_float_min_index(V(x)); break; }+        case 5: { res = gsl_vector_float_min(V(x)); break; }+        default: ERROR(BAD_CODE);+    }+    rp[0] = res;+    OK+} ++int toScalarC(int code, KCVEC(x), RVEC(r)) { +    REQUIRES(rn==1,BAD_SIZE);+    DEBUGMSG("toScalarC");+    KCVVIEW(x);+    double res;+    switch(code) {+        case 0: { res = gsl_blas_dznrm2(V(x)); break; } +        case 1: { res = gsl_blas_dzasum(V(x));  break; }+        default: ERROR(BAD_CODE);+    }+    rp[0] = res;+    OK+}++int toScalarQ(int code, KQVEC(x), FVEC(r)) { +    REQUIRES(rn==1,BAD_SIZE);+    DEBUGMSG("toScalarQ");+    KQVVIEW(x);+    float res;+    switch(code) {+        case 0: { res = gsl_blas_scnrm2(V(x)); break; } +        case 1: { res = gsl_blas_scasum(V(x));  break; }+        default: ERROR(BAD_CODE);+    }+    rp[0] = res;+    OK+}++ inline double sign(double x) {     if(x>0) {         return +1.0;@@ -110,6 +327,16 @@     } } +inline float float_sign(float x) {+    if(x>0) {+        return +1.0;+    } else if (x<0) {+        return -1.0;+    } else {+        return 0.0;+    }+}+ inline gsl_complex complex_abs(gsl_complex z) {     gsl_complex r;     r.dat[0] = gsl_complex_abs(z);@@ -159,7 +386,33 @@     } } +int mapF(int code, KFVEC(x), FVEC(r)) {+    int k;+    REQUIRES(xn == rn,BAD_SIZE);+    DEBUGMSG("mapF");+    switch (code) {+        OP(0,sin)+        OP(1,cos)+        OP(2,tan)+        OP(3,fabs)+        OP(4,asin)+        OP(5,acos)+        OP(6,atan) /* atan2 mediante vectorZip */+        OP(7,sinh)+        OP(8,cosh)+        OP(9,tanh)+        OP(10,gsl_asinh)+        OP(11,gsl_acosh)+        OP(12,gsl_atanh)+        OP(13,exp)+        OP(14,log)+        OP(15,sign)+        OP(16,sqrt)+        default: ERROR(BAD_CODE);+    }+} + int mapCAux(int code, KGCVEC(x), GCVEC(r)) {     int k;     REQUIRES(xn == rn,BAD_SIZE);@@ -194,6 +447,83 @@ }  +gsl_complex_float complex_float_math_fun(gsl_complex (*cf)(gsl_complex), gsl_complex_float a)+{+  gsl_complex c;+  gsl_complex r;++  gsl_complex_float float_r;++  c.dat[0] = a.dat[0];+  c.dat[1] = a.dat[1];++  r = (*cf)(c);++  float_r.dat[0] = r.dat[0];+  float_r.dat[1] = r.dat[1];++  return float_r;+}++gsl_complex_float complex_float_math_op(gsl_complex (*cf)(gsl_complex,gsl_complex), +					gsl_complex_float a,gsl_complex_float b)+{+  gsl_complex c1;+  gsl_complex c2;+  gsl_complex r;++  gsl_complex_float float_r;++  c1.dat[0] = a.dat[0];+  c1.dat[1] = a.dat[1];++  c2.dat[0] = b.dat[0];+  c2.dat[1] = b.dat[1];++  r = (*cf)(c1,c2);++  float_r.dat[0] = r.dat[0];+  float_r.dat[1] = r.dat[1];++  return float_r;+}++#define OPC(C,F) case C: { for(k=0;k<xn;k++) rp[k] = complex_float_math_fun(&F,xp[k]); OK }+#define OPCA(C,F,A,B) case C: { for(k=0;k<xn;k++) rp[k] = complex_float_math_op(&F,A,B); OK }+int mapQAux(int code, KGQVEC(x), GQVEC(r)) {+    int k;+    REQUIRES(xn == rn,BAD_SIZE);+    DEBUGMSG("mapQ");+    switch (code) {+        OPC(0,gsl_complex_sin)+        OPC(1,gsl_complex_cos)+        OPC(2,gsl_complex_tan)+        OPC(3,complex_abs)+        OPC(4,gsl_complex_arcsin)+        OPC(5,gsl_complex_arccos)+        OPC(6,gsl_complex_arctan)+        OPC(7,gsl_complex_sinh)+        OPC(8,gsl_complex_cosh)+        OPC(9,gsl_complex_tanh)+        OPC(10,gsl_complex_arcsinh)+        OPC(11,gsl_complex_arccosh)+        OPC(12,gsl_complex_arctanh)+        OPC(13,gsl_complex_exp)+        OPC(14,gsl_complex_log)+        OPC(15,complex_signum)+        OPC(16,gsl_complex_sqrt)++        // gsl_complex_arg+        // gsl_complex_abs+        default: ERROR(BAD_CODE);+    }+}++int mapQ(int code, KQVEC(x), QVEC(r)) {+    return mapQAux(code, xn, (gsl_complex_float*)xp, rn, (gsl_complex_float*)rp);+}++ int mapValR(int code, double* pval, KRVEC(x), RVEC(r)) {     int k;     double val = *pval;@@ -210,6 +540,22 @@     } } +int mapValF(int code, float* pval, KFVEC(x), FVEC(r)) {+    int k;+    float val = *pval;+    REQUIRES(xn == rn,BAD_SIZE);+    DEBUGMSG("mapValF");+    switch (code) {+        OPV(0,val*xp[k])+        OPV(1,val/xp[k])+        OPV(2,val+xp[k])+        OPV(3,val-xp[k])+        OPV(4,pow(val,xp[k]))+        OPV(5,pow(xp[k],val))+        default: ERROR(BAD_CODE);+    }+}+ int mapValCAux(int code, gsl_complex* pval, KGCVEC(x), GCVEC(r)) {     int k;     gsl_complex val = *pval;@@ -231,6 +577,27 @@ }  +int mapValQAux(int code, gsl_complex_float* pval, KQVEC(x), GQVEC(r)) {+    int k;+    gsl_complex_float val = *pval;+    REQUIRES(xn == rn,BAD_SIZE);+    DEBUGMSG("mapValQ");+    switch (code) {+        OPCA(0,gsl_complex_mul,val,xp[k])+	OPCA(1,gsl_complex_div,val,xp[k])+	OPCA(2,gsl_complex_add,val,xp[k])+	OPCA(3,gsl_complex_sub,val,xp[k])+	OPCA(4,gsl_complex_pow,val,xp[k])+	OPCA(5,gsl_complex_pow,xp[k],val)+        default: ERROR(BAD_CODE);+    }+}++int mapValQ(int code, gsl_complex_float* val, KQVEC(x), QVEC(r)) {+    return mapValQAux(code, val, xn, (gsl_complex_float*)xp, rn, (gsl_complex_float*)rp);+}++ #define OPZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = E(ap[k],bp[k]); OK } #define OPZV(C,msg,E) case C: {DEBUGMSG(msg) res = E(V(r),V(b)); CHECK(res,res); OK } int zipR(int code, KRVEC(a), KRVEC(b), RVEC(r)) {@@ -255,6 +622,28 @@ }  +int zipF(int code, KFVEC(a), KFVEC(b), FVEC(r)) {+    REQUIRES(an == bn && an == rn, BAD_SIZE);+    int k;+    switch(code) {+        OPZE(4,"zipF Pow",pow)+        OPZE(5,"zipF ATan2",atan2)+    }+    KFVVIEW(a);+    KFVVIEW(b);+    FVVIEW(r);+    gsl_vector_float_memcpy(V(r),V(a));+    int res;+    switch(code) {+        OPZV(0,"zipF Add",gsl_vector_float_add)+        OPZV(1,"zipF Sub",gsl_vector_float_sub)+        OPZV(2,"zipF Mul",gsl_vector_float_mul)+        OPZV(3,"zipF Div",gsl_vector_float_div)+        default: ERROR(BAD_CODE);+    }+}++ int zipCAux(int code, KGCVEC(a), KGCVEC(b), GCVEC(r)) {     REQUIRES(an == bn && an == rn, BAD_SIZE);     int k;@@ -279,6 +668,34 @@  int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r)) {     return zipCAux(code, an, (gsl_complex*)ap, bn, (gsl_complex*)bp, rn, (gsl_complex*)rp);+}+++#define OPCZE(C,msg,E) case C: {DEBUGMSG(msg) for(k=0;k<an;k++) rp[k] = complex_float_math_op(&E,ap[k],bp[k]); OK }+int zipQAux(int code, KGQVEC(a), KGQVEC(b), GQVEC(r)) {+    REQUIRES(an == bn && an == rn, BAD_SIZE);+    int k;+    switch(code) {+        OPCZE(0,"zipQ Add",gsl_complex_add)+        OPCZE(1,"zipQ Sub",gsl_complex_sub)+        OPCZE(2,"zipQ Mul",gsl_complex_mul)+        OPCZE(3,"zipQ Div",gsl_complex_div)+        OPCZE(4,"zipQ Pow",gsl_complex_pow)+        //OPZE(5,"zipR ATan2",atan2)+    }+    //KCVVIEW(a);+    //KCVVIEW(b);+    //CVVIEW(r);+    //gsl_vector_memcpy(V(r),V(a));+    //int res;+    switch(code) {+        default: ERROR(BAD_CODE);+    }+}+++int zipQ(int code, KQVEC(a), KQVEC(b), QVEC(r)) {+    return zipQAux(code, an, (gsl_complex_float*)ap, bn, (gsl_complex_float*)bp, rn, (gsl_complex_float*)rp); }  

@@ -0,0 +1,160 @@+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.IO+-- Copyright   :  (c) Alberto Ruiz 2010+-- License     :  GPL+--+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>+-- Stability   :  provisional+-- Portability :  portable+--+-- Display, formatting and IO functions for numeric 'Vector' and 'Matrix'+--+-----------------------------------------------------------------------------++module Numeric.IO (+    dispf, disps, dispcf, vecdisp, latexFormat, format,+    loadMatrix, saveMatrix, fromFile, fileDimensions,+    readMatrix, fromArray2D,+    fscanfVector, fprintfVector, freadVector, fwriteVector+) where++import Data.Packed+import Data.Packed.Internal+import System.Process(readProcess)+import Text.Printf(printf)+import Data.List(intersperse)+import Data.Complex++{- | Creates a string from a matrix given a separator and a function to show each entry. Using+this function the user can easily define any desired display function:++@import Text.Printf(printf)@++@disp = putStr . format \"  \" (printf \"%.2f\")@++-}+format :: (Element t) => String -> (t -> String) -> Matrix t -> String+format sep f m = table sep . map (map f) . toLists $ m++{- | Show a matrix with \"autoscaling\" and a given number of decimal places.++@disp = putStr . disps 2++\> disp $ 120 * (3><4) [1..]+3x4  E3+ 0.12  0.24  0.36  0.48+ 0.60  0.72  0.84  0.96+ 1.08  1.20  1.32  1.44+@+-}+disps :: Int -> Matrix Double -> String+disps d x = sdims x ++ "  " ++ formatScaled d x++{- | Show a matrix with a given number of decimal places.++@disp = putStr . dispf 3++\> disp (1/3 + ident 4)+4x4+1.333  0.333  0.333  0.333+0.333  1.333  0.333  0.333+0.333  0.333  1.333  0.333+0.333  0.333  0.333  1.333+@+-}+dispf :: Int -> Matrix Double -> String+dispf d x = sdims x ++ "\n" ++ formatFixed (if isInt x then 0 else d) x++sdims x = show (rows x) ++ "x" ++ show (cols x)++formatFixed d x = format "  " (printf ("%."++show d++"f")) $ x++isInt = all lookslikeInt . toList . flatten++formatScaled dec t = "E"++show o++"\n" ++ ss+    where ss = format " " (printf fmt. g) t+          g x | o >= 0    = x/10^(o::Int)+              | otherwise = x*10^(-o)+          o = floor $ maximum $ map (logBase 10 . abs) $ toList $ flatten t+          fmt = '%':show (dec+3) ++ '.':show dec ++"f"++{- | Show a vector using a function for showing matrices.++@disp = putStr . vecdisp ('dispf' 2)++\> disp ('linspace' 10 (0,1))+10 |> 0.00  0.11  0.22  0.33  0.44  0.56  0.67  0.78  0.89  1.00+@+-}+vecdisp :: (Element t) => (Matrix t -> String) -> Vector t -> String+vecdisp f v+    = ((show (dim v) ++ " |> ") ++) . (++"\n")+    . unwords . lines .  tail . dropWhile (not . (`elem` " \n"))+    . f . trans . reshape 1+    $ v++-- | Tool to display matrices with latex syntax.+latexFormat :: String -- ^ type of braces: \"matrix\", \"bmatrix\", \"pmatrix\", etc.+            -> String -- ^ Formatted matrix, with elements separated by spaces and newlines+            -> String+latexFormat del tab = "\\begin{"++del++"}\n" ++ f tab ++ "\\end{"++del++"}"+    where f = unlines . intersperse "\\\\" . map unwords . map (intersperse " & " . words) . tail . lines++-- | Pretty print a complex number with at most n decimal digits.+showComplex :: Int -> Complex Double -> String+showComplex d (a:+b)+    | isZero a && isZero b = "0"+    | isZero b = sa+    | isZero a && isOne b = s2++"i"+    | isZero a = sb++"i"+    | isOne b = sa++s3++"i"+    | otherwise = sa++s1++sb++"i"+  where+    sa = shcr d a+    sb = shcr d b+    s1 = if b<0 then "" else "+"+    s2 = if b<0 then "-" else ""+    s3 = if b<0 then "-" else "+"++shcr d a | lookslikeInt a = printf "%.0f" a+         | otherwise      = printf ("%."++show d++"f") a+++lookslikeInt x = show (round x :: Int) ++".0" == shx || "-0.0" == shx+   where shx = show x++isZero x = show x `elem` ["0.0","-0.0"]+isOne  x = show x `elem` ["1.0","-1.0"]++-- | Pretty print a complex matrix with at most n decimal digits.+dispcf :: Int -> Matrix (Complex Double) -> String+dispcf d m = sdims m ++ "\n" ++ format "  " (showComplex d) m++--------------------------------------------------------------------++-- | reads a matrix from a string containing a table of numbers.+readMatrix :: String -> Matrix Double+readMatrix = fromLists . map (map read). map words . filter (not.null) . lines++{- |  obtains the number of rows and columns in an ASCII data file+      (provisionally using unix's wc).+-}+fileDimensions :: FilePath -> IO (Int,Int)+fileDimensions fname = do+    wcres <- readProcess "wc" ["-w",fname] ""+    contents <- readFile fname+    let tot = read . head . words $ wcres+        c   = length . head . dropWhile null . map words . lines $ contents+    if tot > 0+        then return (tot `div` c, c)+        else return (0,0)++-- | Loads a matrix from an ASCII file formatted as a 2D table.+loadMatrix :: FilePath -> IO (Matrix Double)+loadMatrix file = fromFile file =<< fileDimensions file++-- | Loads a matrix from an ASCII file (the number of rows and columns must be known in advance).+fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double)+fromFile filename (r,c) = reshape c `fmap` fscanfVector filename (r*c)+

@@ -1,7 +1,7 @@ ----------------------------------------------------------------------------- {- | Module      :  Numeric.LinearAlgebra-Copyright   :  (c) Alberto Ruiz 2006-9+Copyright   :  (c) Alberto Ruiz 2006-10 License     :  GPL-style  Maintainer  :  Alberto Ruiz (aruiz at um dot es)@@ -10,14 +10,19 @@  This module reexports all normally required functions for Linear Algebra applications. +It also provides instances of standard classes 'Show', 'Read', 'Eq',+'Num', 'Fractional', and 'Floating' for 'Vector' and 'Matrix'.+In arithmetic operations one-component vectors and matrices automatically+expand to match the dimensions of the other operand.+ -} ----------------------------------------------------------------------------- module Numeric.LinearAlgebra (-    module Data.Packed,-    module Numeric.LinearAlgebra.Algorithms,-    module Numeric.LinearAlgebra.Interface+    module Numeric.Container,+    module Numeric.LinearAlgebra.Algorithms ) where -import Data.Packed+import Numeric.Container import Numeric.LinearAlgebra.Algorithms-import Numeric.LinearAlgebra.Interface+import Numeric.Matrix()+import Numeric.Vector()

@@ -1,5 +1,8 @@-{-# OPTIONS_GHC -XFlexibleContexts -XFlexibleInstances #-}+{-# LANGUAGE FlexibleContexts, FlexibleInstances #-} {-# LANGUAGE CPP #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE TypeFamilies #-} ----------------------------------------------------------------------------- {- | Module      :  Numeric.LinearAlgebra.Algorithms@@ -10,7 +13,7 @@ Stability   :  provisional Portability :  uses ffi -Generic interface for the most common functions. Using it we can write higher level algorithms and testing properties for both real and complex matrices.+High level generic interface to common matrix computations.  Specific functions for particular base types can also be explicitly imported from "Numeric.LinearAlgebra.LAPACK".@@ -21,9 +24,6 @@ module Numeric.LinearAlgebra.Algorithms ( -- * Supported types     Field(),--- * Products-    multiply, dot,-    outer, kronecker, -- * Linear Systems     linearSolve,     luSolve,@@ -63,10 +63,9 @@     nullspaceSVD, -- * Norms     Normed(..), NormType(..),+    relativeError, -- * Misc-    ctrans,-    eps, i,-    Linear(..),+    eps, peps, i, -- * Util     haussholder,     unpackQR, unpackHess,@@ -77,18 +76,23 @@   import Data.Packed.Internal hiding ((//))-import Data.Packed.Vector import Data.Packed.Matrix-import Data.Complex-import Numeric.GSL.Vector import Numeric.LinearAlgebra.LAPACK as LAPACK-import Numeric.LinearAlgebra.Linear import Data.List(foldl1') import Data.Array+import Numeric.ContainerBoot hiding ((.*),(*/))  --- | Auxiliary typeclass used to define generic computations for both real and complex matrices.-class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where+{- | Class used to define generic linear algebra computations for both real and complex matrices. Only double precision is supported in this version (we can+transform single precision objects using 'single' and 'double').++-}+class (Product t,+       Convert t,+       Container Vector t,+       Container Matrix t,+       Normed Matrix t,+       Normed Vector t) => Field t where     svd'         :: Matrix t -> (Matrix t, Vector Double, Matrix t)     thinSVD'     :: Matrix t -> (Matrix t, Vector Double, Matrix t)     sv'          :: Matrix t -> Vector Double@@ -107,8 +111,6 @@     qr'          :: Matrix t -> (Matrix t, Matrix t)     hess'        :: Matrix t -> (Matrix t, Matrix t)     schur'       :: Matrix t -> (Matrix t, Matrix t)-    ctrans'      :: Matrix t -> Matrix t-    multiply'    :: Matrix t -> Matrix t -> Matrix t   instance Field Double where@@ -121,7 +123,6 @@     cholSolve' = cholSolveR     linearSolveLS' = linearSolveLSR     linearSolveSVD' = linearSolveSVDR Nothing-    ctrans' = trans     eig' = eigR     eigSH'' = eigS     eigOnly = eigOnlyR@@ -131,7 +132,6 @@     qr' = unpackQR . qrR     hess' = unpackHess hessR     schur' = schurR-    multiply' = multiplyR  instance Field (Complex Double) where #ifdef NOZGESDD@@ -148,7 +148,6 @@     cholSolve' = cholSolveC     linearSolveLS' = linearSolveLSC     linearSolveSVD' = linearSolveSVDC Nothing-    ctrans' = conj . trans     eig' = eigC     eigOnly = eigOnlyC     eigSH'' = eigH@@ -158,7 +157,6 @@     qr' = unpackQR . qrC     hess' = unpackHess hessC     schur' = schurC-    multiply' = multiplyC  -------------------------------------------------------------- @@ -190,7 +188,7 @@ fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t) fullSVD m = (u,d,v) where     (u,s,v) = svd m-    d = diagRect s r c+    d = diagRect 0 s r c     r = rows m     c = cols m @@ -217,7 +215,7 @@ {-# DEPRECATED full "use fullSVD instead" #-} full svdFun m = (u, d ,v) where     (u,s,v) = svdFun m-    d = diagRect s r c+    d = diagRect 0 s r c     r = rows m     c = cols m @@ -326,14 +324,7 @@ schur       :: Field t => Matrix t -> (Matrix t, Matrix t) schur = schur' --- | Generic conjugate transpose.-ctrans :: Field t => Matrix t -> Matrix t-ctrans = ctrans' --- | Matrix product.-multiply :: Field t => Matrix t -> Matrix t -> Matrix t-multiply = {-# SCC "multiply" #-} multiply'- -- | Similar to 'cholSH', but instead of an error (e.g., caused by a matrix not positive definite) it returns 'Nothing'. mbCholSH :: Field t => Matrix t -> Maybe (Matrix t) mbCholSH = {-# SCC "mbCholSH" #-} mbCholSH'@@ -398,74 +389,15 @@ eps :: Double eps =  2.22044604925031e-16 --- | The imaginary unit: @i = 0.0 :+ 1.0@-i :: Complex Double-i = 0:+1 ---- matrix product-mXm :: (Num t, Field t) => Matrix t -> Matrix t -> Matrix t-mXm = multiply---- matrix - vector product-mXv :: (Num t, Field t) => Matrix t -> Vector t -> Vector t-mXv m v = flatten $ m `mXm` (asColumn v)---- vector - matrix product-vXm :: (Num t, Field t) => Vector t -> Matrix t -> Vector t-vXm v m = flatten $ (asRow v) `mXm` m+-- | 1 + 0.5*peps == 1,  1 + 0.6*peps /= 1+peps :: RealFloat x => x+peps = x where x = 2.0**(fromIntegral $ 1-floatDigits x)  ------------------------------------------------------------------------------norm2 :: Vector Double -> Double-norm2 = toScalarR Norm2--norm1 :: Vector Double -> Double-norm1 = toScalarR AbsSum--data NormType = Infinity | PNorm1 | PNorm2 -- PNorm Int--pnormRV PNorm2 = norm2-pnormRV PNorm1 = norm1-pnormRV Infinity = vectorMax . vectorMapR Abs---pnormRV _ = error "pnormRV not yet defined"--pnormCV PNorm2 = norm2 . asReal-pnormCV PNorm1 = norm1 . mapVector magnitude-pnormCV Infinity = vectorMax . mapVector magnitude---pnormCV _ = error "pnormCV not yet defined"--pnormRM PNorm2 m = singularValues m @> 0-pnormRM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (vectorMapR Abs) m-pnormRM Infinity m = vectorMax $ liftMatrix (vectorMapR Abs) m `mXv` constant 1 (cols m)---pnormRM _ _ = error "p norm not yet defined"--pnormCM PNorm2 m = singularValues m @> 0-pnormCM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (mapVector magnitude) m-pnormCM Infinity m = vectorMax $ liftMatrix (mapVector magnitude) m `mXv` constant 1 (cols m)---pnormCM _ _ = error "p norm not yet defined"---- | Objects which have a p-norm.--- Using it you can define convenient shortcuts:------ @norm2 x = pnorm PNorm2 x@------ @frobenius m = norm2 . flatten $ m@-class Normed t where-    pnorm :: NormType -> t -> Double--instance Normed (Vector Double) where-    pnorm = pnormRV--instance Normed (Vector (Complex Double)) where-    pnorm = pnormCV--instance Normed (Matrix Double) where-    pnorm = pnormRM--instance Normed (Matrix (Complex Double)) where-    pnorm = pnormCM+-- | The imaginary unit: @i = 0.0 :+ 1.0@+i :: Complex Double+i = 0:+1  ----------------------------------------------------------------------- @@ -543,7 +475,7 @@               where xs = toList v  zt 0 v = v-zt k v = join [subVector 0 (dim v - k) v, constant 0 k]+zt k v = join [subVector 0 (dim v - k) v, konst 0 k]   unpackQR :: (Field t) => (Matrix t, Vector t) -> (Matrix t, Matrix t)@@ -602,8 +534,8 @@ -- -- @logm = matFunc log@ ---matFunc :: Field t => (Complex Double -> Complex Double) -> Matrix t -> Matrix (Complex Double)-matFunc f m = case diagonalize (complex m) of+matFunc :: (Complex Double -> Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)+matFunc f m = case diagonalize m of     Just (l,v) -> v `mXm` diag (mapVector f l) `mXm` inv v     Nothing -> error "Sorry, matFunc requires a diagonalizable matrix"  @@ -660,11 +592,11 @@  [ 2.0, 2.25  , 0.0,  2.0 ]@ -}-sqrtm :: Field t => Matrix t -> Matrix t+sqrtm ::  Field t => Matrix t -> Matrix t sqrtm = sqrtmInv  sqrtmInv x = fst $ fixedPoint $ iterate f (x, ident (rows x))-    where fixedPoint (a:b:rest) | pnorm PNorm1 (fst a |-| fst b) < eps   = a+    where fixedPoint (a:b:rest) | pnorm PNorm1 (fst a |-| fst b) < peps   = a                                 | otherwise = fixedPoint (b:rest)           fixedPoint _ = error "fixedpoint with impossible inputs"           f (y,z) = (0.5 .* (y |+| inv z),@@ -697,59 +629,72 @@     c = cols l_u     tu = triang r c 0 1     tl = triang r c 0 0-    l = takeColumns r (l_u |*| tl) |+| diagRect (constant 1 r) r r+    l = takeColumns r (l_u |*| tl) |+| diagRect 0 (konst 1 r) r r     u = l_u |*| tu     (p,s) = fixPerm r perm-    l' = (l_u |*| tl) |+| diagRect (constant 1 c) r c+    l' = (l_u |*| tl) |+| diagRect 0 (konst 1 c) r c     u' = takeRows c (l_u |*| tu)     (|+|) = add     (|*|) = mul ---------------------------------------------------+--------------------------------------------------------------------------- --- | Euclidean inner product.-dot :: (Field t) => Vector t -> Vector t -> t-dot u v = multiply r c  @@> (0,0)-    where r = asRow u-          c = asColumn v+data NormType = Infinity | PNorm1 | PNorm2 | Frobenius +class (RealFloat (RealOf t)) => Normed c t where+    pnorm :: NormType -> c t -> RealOf t -{- | Outer product of two vectors.+instance Normed Vector Double where+    pnorm PNorm1    = norm1+    pnorm PNorm2    = norm2+    pnorm Infinity  = normInf+    pnorm Frobenius = norm2 -@\> 'fromList' [1,2,3] \`outer\` 'fromList' [5,2,3]-(3><3)- [  5.0, 2.0, 3.0- , 10.0, 4.0, 6.0- , 15.0, 6.0, 9.0 ]@--}-outer :: (Field t) => Vector t -> Vector t -> Matrix t-outer u v = asColumn u `multiply` asRow v+instance Normed Vector (Complex Double) where+    pnorm PNorm1    = norm1+    pnorm PNorm2    = norm2+    pnorm Infinity  = normInf+    pnorm Frobenius = pnorm PNorm2 -{- | Kronecker product of two matrices.+instance Normed Vector Float where+    pnorm PNorm1    = norm1+    pnorm PNorm2    = norm2+    pnorm Infinity  = normInf+    pnorm Frobenius = pnorm PNorm2 -@m1=(2><3)- [ 1.0,  2.0, 0.0- , 0.0, -1.0, 3.0 ]-m2=(4><3)- [  1.0,  2.0,  3.0- ,  4.0,  5.0,  6.0- ,  7.0,  8.0,  9.0- , 10.0, 11.0, 12.0 ]@+instance Normed Vector (Complex Float) where+    pnorm PNorm1    = norm1+    pnorm PNorm2    = norm2+    pnorm Infinity  = normInf+    pnorm Frobenius = pnorm PNorm2 -@\> kronecker m1 m2-(8><9)- [  1.0,  2.0,  3.0,   2.0,   4.0,   6.0,  0.0,  0.0,  0.0- ,  4.0,  5.0,  6.0,   8.0,  10.0,  12.0,  0.0,  0.0,  0.0- ,  7.0,  8.0,  9.0,  14.0,  16.0,  18.0,  0.0,  0.0,  0.0- , 10.0, 11.0, 12.0,  20.0,  22.0,  24.0,  0.0,  0.0,  0.0- ,  0.0,  0.0,  0.0,  -1.0,  -2.0,  -3.0,  3.0,  6.0,  9.0- ,  0.0,  0.0,  0.0,  -4.0,  -5.0,  -6.0, 12.0, 15.0, 18.0- ,  0.0,  0.0,  0.0,  -7.0,  -8.0,  -9.0, 21.0, 24.0, 27.0- ,  0.0,  0.0,  0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@--}-kronecker :: (Field t) => Matrix t -> Matrix t -> Matrix t-kronecker a b = fromBlocks-              . splitEvery (cols a)-              . map (reshape (cols b))-              . toRows-              $ flatten a `outer` flatten b++instance Normed Matrix Double where+    pnorm PNorm1    = maximum . map (pnorm PNorm1) . toColumns+    pnorm PNorm2    = (@>0) . singularValues+    pnorm Infinity  = pnorm PNorm1 . trans+    pnorm Frobenius = pnorm PNorm2 . flatten++instance Normed Matrix (Complex Double) where+    pnorm PNorm1    = maximum . map (pnorm PNorm1) . toColumns+    pnorm PNorm2    = (@>0) . singularValues+    pnorm Infinity  = pnorm PNorm1 . trans+    pnorm Frobenius = pnorm PNorm2 . flatten++instance Normed Matrix Float where+    pnorm PNorm1    = maximum . map (pnorm PNorm1) . toColumns+    pnorm PNorm2    = realToFrac . (@>0) . singularValues . double+    pnorm Infinity  = pnorm PNorm1 . trans+    pnorm Frobenius = pnorm PNorm2 . flatten++instance Normed Matrix (Complex Float) where+    pnorm PNorm1    = maximum . map (pnorm PNorm1) . toColumns+    pnorm PNorm2    = realToFrac . (@>0) . singularValues . double+    pnorm Infinity  = pnorm PNorm1 . trans+    pnorm Frobenius = pnorm PNorm2 . flatten++-- | Approximate number of common digits in the maximum element.+relativeError :: (Normed c t, Container c t) => c t -> c t -> Int+relativeError x y = dig (norm (x `sub` y) / norm x)+    where norm = pnorm Infinity+          dig r = round $ -logBase 10 (realToFrac r :: Double)

@@ -1,218 +0,0 @@-{-# LANGUAGE UndecidableInstances, FlexibleInstances #-}-------------------------------------------------------------------------------{- |-Module      :  Numeric.LinearAlgebra.Instances-Copyright   :  (c) Alberto Ruiz 2006-License     :  GPL-style--Maintainer  :  Alberto Ruiz (aruiz at um dot es)-Stability   :  provisional-Portability :  portable--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.---}--------------------------------------------------------------------------------module Numeric.LinearAlgebra.Instances(-) where--import Numeric.LinearAlgebra.Linear-import Numeric.GSL.Vector-import Data.Packed.Matrix-import Data.Complex-import Data.List(transpose,intersperse)-import Data.Packed.Internal.Vector--#ifndef VECTOR-import Foreign(Storable)-#endif----------------------------------------------------------------------instance (Show a, Element a) => (Show (Matrix a)) where-    show m = (sizes++) . dsp . map (map show) . toLists $ m-        where sizes = "("++show (rows m)++"><"++show (cols m)++")\n"--dsp as = (++" ]") . (" ["++) . init . drop 2 . unlines . map (" , "++) . map unwords' $ transpose mtp-    where-        mt = transpose as-        longs = map (maximum . map length) mt-        mtp = zipWith (\a b -> map (pad a) b) longs mt-        pad n str = replicate (n - length str) ' ' ++ str-        unwords' = concat . intersperse ", "--#ifndef VECTOR--instance (Show a, Storable a) => (Show (Vector a)) where-    show v = (show (dim v))++" |> " ++ show (toList v)--#endif----------------------------------------------------------------------instance (Element a, Read a) => Read (Matrix a) where-    readsPrec _ s = [((rs><cs) . read $ listnums, rest)]-        where (thing,rest) = breakAt ']' s-              (dims,listnums) = breakAt ')' thing-              cs = read . init . fst. breakAt ')' . snd . breakAt '<' $ dims-              rs = read . snd . breakAt '(' .init . fst . breakAt '>' $ dims--#ifdef VECTOR--instance (Element a, Read a) => Read (Vector a) where-    readsPrec _ s = [(fromList . read $ listnums, rest)]-        where (thing,trest) = breakAt ']' s-              (dims,listnums) = breakAt ' ' (dropWhile (==' ') thing)-              rest = drop 31 trest-#else--instance (Element a, Read a) => Read (Vector a) where-    readsPrec _ s = [((d |>) . read $ listnums, rest)]-        where (thing,rest) = breakAt ']' s-              (dims,listnums) = breakAt '>' thing-              d = read . init . fst . breakAt '|' $ dims--#endif--breakAt c l = (a++[c],tail b) where-    (a,b) = break (==c) l----------------------------------------------------------------------adaptScalar f1 f2 f3 x y-    | dim x == 1 = f1   (x@>0) y-    | dim y == 1 = f3 x (y@>0)-    | otherwise = f2 x y--#ifndef VECTOR--instance Linear Vector a => Eq (Vector a) where-    (==) = equal--#endif--instance Num (Vector Double) where-    (+) = adaptScalar addConstant add (flip addConstant)-    negate = scale (-1)-    (*) = adaptScalar scale mul (flip scale)-    signum = vectorMapR Sign-    abs = vectorMapR Abs-    fromInteger = fromList . return . fromInteger--instance Num (Vector (Complex Double)) where-    (+) = adaptScalar addConstant add (flip addConstant)-    negate = scale (-1)-    (*) = adaptScalar scale mul (flip scale)-    signum = vectorMapC Sign-    abs = vectorMapC Abs-    fromInteger = fromList . return . fromInteger--instance Linear Matrix a => Eq (Matrix a) where-    (==) = equal--instance (Linear Matrix a, Num (Vector a)) => Num (Matrix a) where-    (+) = liftMatrix2Auto (+)-    (-) = liftMatrix2Auto (-)-    negate = liftMatrix negate-    (*) = liftMatrix2Auto (*)-    signum = liftMatrix signum-    abs = liftMatrix abs-    fromInteger = (1><1) . return . fromInteger-------------------------------------------------------instance (Linear Vector a, Num (Vector a)) => Fractional (Vector a) where-    fromRational n = fromList [fromRational n]-    (/) = adaptScalar f divide g where-        r `f` v = scaleRecip r v-        v `g` r = scale (recip r) v-----------------------------------------------------------instance (Linear Vector a, Fractional (Vector a), Num (Matrix a)) => Fractional (Matrix a) where-    fromRational n = (1><1) [fromRational n]-    (/) = liftMatrix2Auto (/)-------------------------------------------------------------instance Floating (Vector Double) where-    sin   = vectorMapR Sin-    cos   = vectorMapR Cos-    tan   = vectorMapR Tan-    asin  = vectorMapR ASin-    acos  = vectorMapR ACos-    atan  = vectorMapR ATan-    sinh  = vectorMapR Sinh-    cosh  = vectorMapR Cosh-    tanh  = vectorMapR Tanh-    asinh = vectorMapR ASinh-    acosh = vectorMapR ACosh-    atanh = vectorMapR ATanh-    exp   = vectorMapR Exp-    log   = vectorMapR Log-    sqrt  = vectorMapR Sqrt-    (**)  = adaptScalar (vectorMapValR PowSV) (vectorZipR Pow) (flip (vectorMapValR PowVS))-    pi    = fromList [pi]-----------------------------------------------------------------instance Floating (Vector (Complex Double)) where-    sin   = vectorMapC Sin-    cos   = vectorMapC Cos-    tan   = vectorMapC Tan-    asin  = vectorMapC ASin-    acos  = vectorMapC ACos-    atan  = vectorMapC ATan-    sinh  = vectorMapC Sinh-    cosh  = vectorMapC Cosh-    tanh  = vectorMapC Tanh-    asinh = vectorMapC ASinh-    acosh = vectorMapC ACosh-    atanh = vectorMapC ATanh-    exp   = vectorMapC Exp-    log   = vectorMapC Log-    sqrt  = vectorMapC Sqrt-    (**)  = adaptScalar (vectorMapValC PowSV) (vectorZipC Pow) (flip (vectorMapValC PowVS))-    pi    = fromList [pi]---------------------------------------------------------------instance (Linear Vector a, Floating (Vector a), Fractional (Matrix a)) => Floating (Matrix a) where-    sin   = liftMatrix sin-    cos   = liftMatrix cos-    tan   = liftMatrix tan-    asin  = liftMatrix asin-    acos  = liftMatrix acos-    atan  = liftMatrix atan-    sinh  = liftMatrix sinh-    cosh  = liftMatrix cosh-    tanh  = liftMatrix tanh-    asinh = liftMatrix asinh-    acosh = liftMatrix acosh-    atanh = liftMatrix atanh-    exp   = liftMatrix exp-    log   = liftMatrix log-    (**)  = liftMatrix2Auto (**)-    sqrt  = liftMatrix sqrt-    pi    = (1><1) [pi]--------------------------------------------------------------------- instance (Storable a, Num (Vector a)) => Monoid (Vector a) where---     mempty = 0 { idim = 0 }---     mappend a b = mconcat [a,b]---     mconcat = j . filter ((>0).dim)---         where j [] = mempty---               j l  = join l--------------------------------------------------------------------- instance (NFData a, Storable a) => NFData (Vector a) where---     rnf = rnf . (@>0)------ instance (NFData a, Element a) => NFData (Matrix a) where---     rnf = rnf . flatten-

@@ -1,117 +0,0 @@-{-# OPTIONS_GHC -fglasgow-exts #-}-------------------------------------------------------------------------------{- |-Module      :  Numeric.LinearAlgebra.Interface-Copyright   :  (c) Alberto Ruiz 2007-License     :  GPL-style--Maintainer  :  Alberto Ruiz (aruiz at um dot es)-Stability   :  provisional-Portability :  portable--Some useful operators, and 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.----}--------------------------------------------------------------------------------module Numeric.LinearAlgebra.Interface(-    (<>),(<.>),-    (<\>),-    (.*),(*/),-    (<|>),(<->),-) where--import Numeric.LinearAlgebra.Instances()-import Data.Packed.Vector-import Data.Packed.Matrix-import Numeric.LinearAlgebra.Algorithms--class Mul a b c | a b -> c where- infixl 7 <>- -- | Matrix-matrix, matrix-vector, and vector-matrix products.- (<>) :: Field t => a t -> b t -> c t--instance Mul Matrix Matrix Matrix where-    (<>) = multiply--instance Mul Matrix Vector Vector where-    (<>) m v = flatten $ m <> (asColumn v)--instance Mul Vector Matrix Vector where-    (<>) v m = flatten $ (asRow v) <> m--------------------------------------------------------- | Dot product: @u \<.\> v = dot u v@-(<.>) :: (Field t) => Vector t -> Vector t -> t-infixl 7 <.>-(<.>) = dot--------------------------------------------------------{-# DEPRECATED (.*) "use scale a x or scalar a * x" #-}---- -- | @x .* a = scale x a@--- (.*) :: (Linear c a) => a -> c a -> c a-infixl 7 .*-a .* x = scale a x--------------------------------------------------------{-# DEPRECATED (*/) "use scale (recip a) x or x / scalar a" #-}---- -- | @a *\/ x = scale (recip x) a@--- (*/) :: (Linear c a) => c a -> a -> c a-infixl 7 */-v */ x = scale (recip x) v---- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD).-(<\>) :: (Field a) => Matrix a -> Vector a -> Vector a-infixl 7 <\>-m <\> v = flatten (linearSolveSVD m (reshape 1 v))----------------------------------------------------{-# DEPRECATED (<|>) "define operator a & b = fromBlocks[[a,b]] and use asRow/asColumn to join vectors" #-}-{-# DEPRECATED (<->) "define operator a // b = fromBlocks[[a],[b]] and use asRow/asColumn to join vectors" #-}--class Joinable a b where-    joinH :: Element t => a t -> b t -> Matrix t-    joinV :: Element t => a t -> b t -> Matrix t--instance Joinable Matrix Matrix where-    joinH m1 m2 = fromBlocks [[m1,m2]]-    joinV m1 m2 = fromBlocks [[m1],[m2]]--instance Joinable Matrix Vector where-    joinH m v = joinH m (asColumn v)-    joinV m v = joinV m (asRow v)--instance Joinable Vector Matrix where-    joinH v m = joinH (asColumn v) m-    joinV v m = joinV (asRow v) m--infixl 4 <|>-infixl 3 <->--{-- - | Horizontal concatenation of matrices and vectors:--@> (ident 3 \<-\> 3 * ident 3) \<|\> fromList [1..6.0]-(6><4)- [ 1.0, 0.0, 0.0, 1.0- , 0.0, 1.0, 0.0, 2.0- , 0.0, 0.0, 1.0, 3.0- , 3.0, 0.0, 0.0, 4.0- , 0.0, 3.0, 0.0, 5.0- , 0.0, 0.0, 3.0, 6.0 ]@--}--- (<|>) :: (Element t, Joinable a b) => a t -> b t -> Matrix t-a <|> b = joinH a b---- -- | Vertical concatenation of matrices and vectors.--- (<->) :: (Element t, Joinable a b) => a t -> b t -> Matrix t-a <-> b = joinV a b-

@@ -14,7 +14,7 @@  module Numeric.LinearAlgebra.LAPACK (     -- * Matrix product-    multiplyR, multiplyC,+    multiplyR, multiplyC, multiplyF, multiplyQ,     -- * Linear systems     linearSolveR, linearSolveC,     lusR, lusC,@@ -43,7 +43,8 @@  import Data.Packed.Internal import Data.Packed.Matrix-import Data.Complex+--import Data.Complex+import Numeric.Conversion import Numeric.GSL.Vector(vectorMapValR, FunCodeSV(Scale)) import Foreign import Foreign.C.Types (CInt)@@ -51,8 +52,10 @@  ----------------------------------------------------------------------------------- -foreign import ccall "LAPACK/lapack-aux.h multiplyR" dgemmc :: CInt -> CInt -> TMMM-foreign import ccall "LAPACK/lapack-aux.h multiplyC" zgemmc :: CInt -> CInt -> TCMCMCM+foreign import ccall "multiplyR" dgemmc :: CInt -> CInt -> TMMM+foreign import ccall "multiplyC" zgemmc :: CInt -> CInt -> TCMCMCM+foreign import ccall "multiplyF" sgemmc :: CInt -> CInt -> TFMFMFM+foreign import ccall "multiplyQ" cgemmc :: CInt -> CInt -> TQMQMQM  isT MF{} = 0 isT MC{} = 1@@ -69,12 +72,20 @@  -- | Matrix product based on BLAS's /dgemm/. multiplyR :: Matrix Double -> Matrix Double -> Matrix Double-multiplyR a b = multiplyAux dgemmc "dgemmc" a b+multiplyR a b = {-# SCC "multiplyR" #-} multiplyAux dgemmc "dgemmc" a b  -- | Matrix product based on BLAS's /zgemm/. multiplyC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) multiplyC a b = multiplyAux zgemmc "zgemmc" a b +-- | Matrix product based on BLAS's /sgemm/.+multiplyF :: Matrix Float -> Matrix Float -> Matrix Float+multiplyF a b = multiplyAux sgemmc "sgemmc" a b++-- | Matrix product based on BLAS's /cgemm/.+multiplyQ :: Matrix (Complex Float) -> Matrix (Complex Float) -> Matrix (Complex Float)+multiplyQ a b = multiplyAux cgemmc "cgemmc" a b+ ----------------------------------------------------------------------------- foreign import ccall "svd_l_R" dgesvd :: TMMVM foreign import ccall "svd_l_C" zgesvd :: TCMCMVCM@@ -248,14 +259,14 @@   where r = rows m         g ra ca pa = dgeev ra ca pa 0 0 nullPtr -fixeig1 s = toComplex (subVector 0 r (asReal s), subVector r r (asReal s))+fixeig1 s = toComplex' (subVector 0 r (asReal s), subVector r r (asReal s))     where r = dim s  fixeig  []  _ =  []-fixeig [_] [v] = [comp v]+fixeig [_] [v] = [comp' v] fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs)-    | r1 == r2 && i1 == (-i2) = toComplex (v1,v2) : toComplex (v1,scale (-1) v2) : fixeig r vs-    | otherwise = comp v1 : fixeig ((r2:+i2):r) (v2:vs)+    | r1 == r2 && i1 == (-i2) = toComplex' (v1,v2) : toComplex' (v1,scale (-1) v2) : fixeig r vs+    | otherwise = comp' v1 : fixeig ((r2:+i2):r) (v2:vs)   where scale = vectorMapValR Scale fixeig _ _ = error "fixeig with impossible inputs" 

@@ -11,15 +11,25 @@  #define MIN(A,B) ((A)<(B)?(A):(B)) #define MAX(A,B) ((A)>(B)?(A):(B))- ++// #define DBGL+ #ifdef DBGL-#define DEBUGMSG(M) printf("LAPACK Wrapper "M"\n: "); size_t t0 = time(NULL);-#define OK MACRO(printf("%ld s\n",time(0)-t0); return 0;);+#define DEBUGMSG(M) printf("\nLAPACK "M"\n"); #else #define DEBUGMSG(M)-#define OK return 0; #endif +#define OK return 0;++// #ifdef DBGL+// #define DEBUGMSG(M) printf("LAPACK Wrapper "M"\n: "); size_t t0 = time(NULL);+// #define OK MACRO(printf("%ld s\n",time(0)-t0); return 0;);+// #else+// #define DEBUGMSG(M)+// #define OK return 0;+// #endif+ #define TRACEMAT(M) {int q; printf(" %d x %d: ",M##r,M##c); \                      for(q=0;q<M##r*M##c;q++) printf("%.1f ",M##p[q]); printf("\n");} @@ -177,9 +187,9 @@             ldvt = q;         }     }DEBUGMSG("svd_l_C");-    double *B = (double*)malloc(2*m*n*sizeof(double));+    doublecomplex *B = (doublecomplex*)malloc(m*n*sizeof(doublecomplex));     CHECK(!B,MEM);-    memcpy(B,ap,m*n*2*sizeof(double));+    memcpy(B,ap,m*n*sizeof(doublecomplex));      double *rwork = (double*) malloc(5*q*sizeof(double));     CHECK(!rwork,MEM);@@ -188,21 +198,21 @@     // ask for optimal lwork     doublecomplex ans;     zgesvd_ (jobu,jobvt,-             &m,&n,(doublecomplex*)B,&m,+             &m,&n,B,&m,              sp,-             (doublecomplex*)up,&m,-             (doublecomplex*)vp,&ldvt,+             up,&m,+             vp,&ldvt,              &ans, &lwork,              rwork,              &res);     lwork = ceil(ans.r);-    doublecomplex * work = (doublecomplex*)malloc(lwork*2*sizeof(double));+    doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));     CHECK(!work,MEM);     zgesvd_ (jobu,jobvt,-             &m,&n,(doublecomplex*)B,&m,+             &m,&n,B,&m,              sp,-             (doublecomplex*)up,&m,-             (doublecomplex*)vp,&ldvt,+             up,&m,+             vp,&ldvt,              work, &lwork,              rwork,              &res);@@ -257,12 +267,12 @@     integer res;     // ask for optimal lwk     doublecomplex ans;-    zgesdd_ (jobz,&m,&n,B,&m,sp,(doublecomplex*)up,&m,(doublecomplex*)vp,&ldvt,&ans,&lwk,rwk,iwk,&res);+    zgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,rwk,iwk,&res);     lwk = ans.r;     //printf("lwk = %ld\n",lwk);     doublecomplex * workv = (doublecomplex*)malloc(lwk*sizeof(doublecomplex));     CHECK(!workv,MEM);-    zgesdd_ (jobz,&m,&n,B,&m,sp,(doublecomplex*)up,&m,(doublecomplex*)vp,&ldvt,workv,&lwk,rwk,iwk,&res);+    zgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,workv,&lwk,rwk,iwk,&res);     //printf("res = %ld\n",res);     CHECK(res,res);     free(workv); // printf("freed workv\n");@@ -293,10 +303,10 @@     doublecomplex ans;     //printf("ask zgeev\n");     zgeev_  (&jobvl,&jobvr,-             &n,(doublecomplex*)B,&n,-             (doublecomplex*)sp,-             (doublecomplex*)up,&n,-             (doublecomplex*)vp,&n,+             &n,B,&n,+             sp,+             up,&n,+             vp,&n,              &ans, &lwork,              rwork,              &res);@@ -306,10 +316,10 @@     CHECK(!work,MEM);     //printf("zgeev\n");     zgeev_  (&jobvl,&jobvr,-             &n,(doublecomplex*)B,&n,-             (doublecomplex*)sp,-             (doublecomplex*)up,&n,-             (doublecomplex*)vp,&n,+             &n,B,&n,+             sp,+             up,&n,+             vp,&n,              work, &lwork,              rwork,              &res);@@ -342,7 +352,7 @@     //printf("ask dgeev\n");     dgeev_  (&jobvl,&jobvr,              &n,B,&n,-             sp, sp+n,+             (double*)sp, (double*)sp+n,              up,&n,              vp,&n,              &ans, &lwork,@@ -354,7 +364,7 @@     //printf("dgeev\n");     dgeev_  (&jobvl,&jobvr,              &n,B,&n,-             sp, sp+n,+             (double*)sp, (double*)sp+n,              up,&n,              vp,&n,              work, &lwork,@@ -419,7 +429,7 @@     doublecomplex ans;     //printf("ask zheev\n");     zheev_  (&jobz,&uplo,-             &n,(doublecomplex*)vp,&n,+             &n,vp,&n,              sp,              &ans, &lwork,              rwork,@@ -429,7 +439,7 @@     doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));     CHECK(!work,MEM);     zheev_  (&jobz,&uplo,-             &n,(doublecomplex*)vp,&n,+             &n,vp,&n,              sp,              work, &lwork,              rwork,@@ -473,15 +483,15 @@     integer nhrs = bc;     REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);     DEBUGMSG("linearSolveC_l");-    double*AC = (double*)malloc(2*n*n*sizeof(double));-    memcpy(AC,ap,2*n*n*sizeof(double));-    memcpy(xp,bp,2*n*nhrs*sizeof(double));+    doublecomplex*AC = (doublecomplex*)malloc(n*n*sizeof(doublecomplex));+    memcpy(AC,ap,n*n*sizeof(doublecomplex));+    memcpy(xp,bp,n*nhrs*sizeof(doublecomplex));     integer * ipiv = (integer*)malloc(n*sizeof(integer));     integer res;     zgesv_  (&n,&nhrs,-             (doublecomplex*)AC, &n,+             AC, &n,              ipiv,-             (doublecomplex*)xp, &n,+             xp, &n,              &res);     if(res>0) {         return SINGULAR;@@ -517,12 +527,12 @@     integer nhrs = bc;     REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);     DEBUGMSG("cholSolveC_l");-    memcpy(xp,bp,2*n*nhrs*sizeof(double));+    memcpy(xp,bp,n*nhrs*sizeof(doublecomplex));     integer res;     zpotrs_  ("U",              &n,&nhrs,              (doublecomplex*)ap, &n,-             (doublecomplex*)xp, &n,+             xp, &n,              &res);     CHECK(res,res);     OK@@ -581,31 +591,30 @@     integer ldb = xr;     REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE);     DEBUGMSG("linearSolveLSC_l");-    double*AC = (double*)malloc(2*m*n*sizeof(double));-    memcpy(AC,ap,2*m*n*sizeof(double));-    memcpy(AC,ap,2*m*n*sizeof(double));+    doublecomplex*AC = (doublecomplex*)malloc(m*n*sizeof(doublecomplex));+    memcpy(AC,ap,m*n*sizeof(doublecomplex));     if (m>=n) {-        memcpy(xp,bp,2*m*nrhs*sizeof(double));+        memcpy(xp,bp,m*nrhs*sizeof(doublecomplex));     } else {         int k;         for(k = 0; k<nrhs; k++) {-            memcpy(xp+2*ldb*k,bp+2*m*k,m*2*sizeof(double));+            memcpy(xp+ldb*k,bp+m*k,m*sizeof(doublecomplex));         }     }     integer res;     integer lwork = -1;     doublecomplex ans;     zgels_  ("N",&m,&n,&nrhs,-             (doublecomplex*)AC,&m,-             (doublecomplex*)xp,&ldb,+             AC,&m,+             xp,&ldb,              &ans,&lwork,              &res);     lwork = ceil(ans.r);     //printf("ans = %d\n",lwork);     doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));     zgels_  ("N",&m,&n,&nrhs,-             (doublecomplex*)AC,&m,-             (doublecomplex*)xp,&ldb,+             AC,&m,+             xp,&ldb,              work,&lwork,              &res);     if(res>0) {@@ -685,16 +694,16 @@     integer ldb = xr;     REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE);     DEBUGMSG("linearSolveSVDC_l");-    double*AC = (double*)malloc(2*m*n*sizeof(double));+    doublecomplex*AC = (doublecomplex*)malloc(m*n*sizeof(doublecomplex));     double*S = (double*)malloc(MIN(m,n)*sizeof(double));     double*RWORK = (double*)malloc(5*MIN(m,n)*sizeof(double));-    memcpy(AC,ap,2*m*n*sizeof(double));+    memcpy(AC,ap,m*n*sizeof(doublecomplex));     if (m>=n) {-        memcpy(xp,bp,2*m*nrhs*sizeof(double));+        memcpy(xp,bp,m*nrhs*sizeof(doublecomplex));     } else {         int k;         for(k = 0; k<nrhs; k++) {-            memcpy(xp+2*ldb*k,bp+2*m*k,m*2*sizeof(double));+            memcpy(xp+ldb*k,bp+m*k,m*sizeof(doublecomplex));         }     }     integer res;@@ -702,8 +711,8 @@     integer rank;     doublecomplex ans;     zgelss_  (&m,&n,&nrhs,-             (doublecomplex*)AC,&m,-             (doublecomplex*)xp,&ldb,+             AC,&m,+             xp,&ldb,              S,              &rcond,&rank,              &ans,&lwork,@@ -713,8 +722,8 @@     //printf("ans = %d\n",lwork);     doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));     zgelss_  (&m,&n,&nrhs,-             (doublecomplex*)AC,&m,-             (doublecomplex*)xp,&ldb,+             AC,&m,+             xp,&ldb,              S,              &rcond,&rank,              work,&lwork,@@ -740,14 +749,14 @@     memcpy(lp,ap,n*n*sizeof(doublecomplex));     char uplo = 'U';     integer res;-    zpotrf_ (&uplo,&n,(doublecomplex*)lp,&n,&res);+    zpotrf_ (&uplo,&n,lp,&n,&res);     CHECK(res>0,NODEFPOS);     CHECK(res,res);     doublecomplex zero = {0.,0.};     int r,c;     for (r=0; r<lr-1; r++) {         for(c=r+1; c<lc; c++) {-            ((doublecomplex*)lp)[r*lc+c] = zero;+            lp[r*lc+c] = zero;         }     }     OK@@ -800,7 +809,7 @@     CHECK(!WORK,MEM);     memcpy(rp,ap,m*n*sizeof(doublecomplex));     integer res;-    zgeqr2_ (&m,&n,(doublecomplex*)rp,&m,(doublecomplex*)taup,WORK,&res);+    zgeqr2_ (&m,&n,rp,&m,taup,WORK,&res);     CHECK(res,res);     free(WORK);     OK@@ -838,7 +847,7 @@     memcpy(rp,ap,m*n*sizeof(doublecomplex));     integer res;     integer one = 1;-    zgehrd_ (&n,&one,&n,(doublecomplex*)rp,&n,(doublecomplex*)taup,WORK,&lwork,&res);+    zgehrd_ (&n,&one,&n,rp,&n,taup,WORK,&lwork,&res);     CHECK(res,res);     free(WORK);     OK@@ -894,8 +903,8 @@     double *RWORK = (double*)malloc(n*sizeof(double));     integer res;     integer sdim;-    zgees_ ("V","N",NULL,&n,(doublecomplex*)sp,&n,&sdim,W,-                            (doublecomplex*)up,&n,+    zgees_ ("V","N",NULL,&n,sp,&n,&sdim,W,+                            up,&n,                             WORK,&lwork,RWORK,BWORK,&res);     if(res>0) {         return NOCONVER;@@ -940,7 +949,7 @@     integer* auxipiv = (integer*)malloc(mn*sizeof(integer));     memcpy(rp,ap,m*n*sizeof(doublecomplex));     integer res;-    zgetrf_ (&m,&n,(doublecomplex*)rp,&m,auxipiv,&res);+    zgetrf_ (&m,&n,rp,&m,auxipiv,&res);     if(res>0) {         res = 0; // fixme     }@@ -990,7 +999,7 @@     }     integer res;     memcpy(xp,bp,mrhs*nrhs*sizeof(doublecomplex));-    zgetrs_ ("N",&n,&nrhs,(doublecomplex*)ap,&m,auxipiv,(doublecomplex*)xp,&mrhs,&res);+    zgetrs_ ("N",&n,&nrhs,(doublecomplex*)ap,&m,auxipiv,xp,&mrhs,&res);     CHECK(res,res);     free(auxipiv);     OK@@ -1004,6 +1013,7 @@  int multiplyR(int ta, int tb, KDMAT(a),KDMAT(b),DMAT(r)) {     //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);+    DEBUGMSG("dgemm_");     integer m = ta?ac:ar;     integer n = tb?br:bc;     integer k = ta?ar:ac;@@ -1022,6 +1032,7 @@  int multiplyC(int ta, int tb, KCMAT(a),KCMAT(b),CMAT(r)) {     //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);+    DEBUGMSG("zgemm_");     integer m = ta?ac:ar;     integer n = tb?br:bc;     integer k = ta?ar:ac;@@ -1031,14 +1042,67 @@     doublecomplex alpha = {1,0};     doublecomplex beta = {0,0};     zgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,-           (doublecomplex*)ap,&lda,-           (doublecomplex*)bp,&ldb,&beta,-           (doublecomplex*)rp,&ldc);+           ap,&lda,+           bp,&ldb,&beta,+           rp,&ldc);     OK } +void sgemm_(char *, char *, integer *, integer *, integer *,+            float *, const float *, integer *, const float *,+           integer *, float *, float *, integer *);++int multiplyF(int ta, int tb, KFMAT(a),KFMAT(b),FMAT(r)) {+    //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);+    DEBUGMSG("sgemm_");+    integer m = ta?ac:ar;+    integer n = tb?br:bc;+    integer k = ta?ar:ac;+    integer lda = ar;+    integer ldb = br;+    integer ldc = rr;+    float alpha = 1;+    float beta = 0;+    sgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,ap,&lda,bp,&ldb,&beta,rp,&ldc);+    OK+}++void cgemm_(char *, char *, integer *, integer *, integer *,+           complex *, const complex *, integer *, const complex *,+           integer *, complex *, complex *, integer *);++int multiplyQ(int ta, int tb, KQMAT(a),KQMAT(b),QMAT(r)) {+    //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);+    DEBUGMSG("cgemm_");+    integer m = ta?ac:ar;+    integer n = tb?br:bc;+    integer k = ta?ar:ac;+    integer lda = ar;+    integer ldb = br;+    integer ldc = rr;+    complex alpha = {1,0};+    complex beta = {0,0};+    cgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,+           ap,&lda,+           bp,&ldb,&beta,+           rp,&ldc);+    OK+}+ //////////////////// transpose ///////////////////////// +int transF(KFMAT(x),FMAT(t)) {+    REQUIRES(xr==tc && xc==tr,BAD_SIZE);+    DEBUGMSG("transF");+    int i,j;+    for (i=0; i<tr; i++) {+        for (j=0; j<tc; j++) {+        tp[i*tc+j] = xp[j*xc+i];+        }+    }+    OK+}+ int transR(KDMAT(x),DMAT(t)) {     REQUIRES(xr==tc && xc==tr,BAD_SIZE);     DEBUGMSG("transR");@@ -1051,20 +1115,55 @@     OK } +int transQ(KQMAT(x),QMAT(t)) {+    REQUIRES(xr==tc && xc==tr,BAD_SIZE);+    DEBUGMSG("transQ");+    int i,j;+    for (i=0; i<tr; i++) {+        for (j=0; j<tc; j++) {+        tp[i*tc+j] = xp[j*xc+i];+        }+    }+    OK+}+ int transC(KCMAT(x),CMAT(t)) {     REQUIRES(xr==tc && xc==tr,BAD_SIZE);     DEBUGMSG("transC");     int i,j;     for (i=0; i<tr; i++) {         for (j=0; j<tc; j++) {-        ((doublecomplex*)tp)[i*tc+j] = ((doublecomplex*)xp)[j*xc+i];+        tp[i*tc+j] = xp[j*xc+i];         }     }     OK } +int transP(KPMAT(x), PMAT(t)) {+    REQUIRES(xr==tc && xc==tr,BAD_SIZE);+    REQUIRES(xs==ts,NOCONVER);+    DEBUGMSG("transP");+    int i,j;+    for (i=0; i<tr; i++) {+        for (j=0; j<tc; j++) {+	  memcpy(tp+(i*tc+j)*xs,xp +(j*xc+i)*xs,xs);+        }+    }+    OK+}+ //////////////////// constant ///////////////////////// +int constantF(float * pval, FVEC(r)) {+    DEBUGMSG("constantF")+    int k;+    double val = *pval;+    for(k=0;k<rn;k++) {+        rp[k]=val;+    }+    OK+}+ int constantR(double * pval, DVEC(r)) {     DEBUGMSG("constantR")     int k;@@ -1075,12 +1174,76 @@     OK } +int constantQ(complex* pval, QVEC(r)) {+    DEBUGMSG("constantQ")+    int k;+    complex val = *pval;+    for(k=0;k<rn;k++) {+        rp[k]=val;+    }+    OK+}+ int constantC(doublecomplex* pval, CVEC(r)) {     DEBUGMSG("constantC")     int k;     doublecomplex val = *pval;     for(k=0;k<rn;k++) {-        ((doublecomplex*)rp)[k]=val;+        rp[k]=val;     }     OK }++int constantP(void* pval, PVEC(r)) {+    DEBUGMSG("constantP")+    int k;+    for(k=0;k<rn;k++) {+      memcpy(rp+k*rs,pval,rs);+    }+    OK+}++//////////////////// float-double conversion /////////////////////////++int float2double(FVEC(x),DVEC(y)) {+    DEBUGMSG("float2double")+    int k;+    for(k=0;k<xn;k++) {+        yp[k]=xp[k];+    }+    OK+}++int double2float(DVEC(x),FVEC(y)) {+    DEBUGMSG("double2float")+    int k;+    for(k=0;k<xn;k++) {+        yp[k]=xp[k];+    }+    OK+}++//////////////////// conjugate /////////////////////////++int conjugateQ(KQVEC(x),QVEC(t)) {+    REQUIRES(xn==tn,BAD_SIZE);+    DEBUGMSG("conjugateQ");+    int k;+    for(k=0;k<xn;k++) {+        tp[k].r =  xp[k].r;+        tp[k].i = -xp[k].i;+    }+    OK+}++int conjugateC(KCVEC(x),CVEC(t)) {+    REQUIRES(xn==tn,BAD_SIZE);+    DEBUGMSG("conjugateC");+    int k;+    for(k=0;k<xn;k++) {+        tp[k].r =  xp[k].r;+        tp[k].i = -xp[k].i;+    }+    OK+}+

@@ -40,26 +40,52 @@  /********************************************************/ +#define FVEC(A) int A##n, float*A##p #define DVEC(A) int A##n, double*A##p-#define CVEC(A) int A##n, double*A##p+#define QVEC(A) int A##n, complex*A##p+#define CVEC(A) int A##n, doublecomplex*A##p+#define PVEC(A) int A##n, void* A##p, int A##s+#define FMAT(A) int A##r, int A##c, float* A##p #define DMAT(A) int A##r, int A##c, double* A##p-#define CMAT(A) int A##r, int A##c, double* A##p+#define QMAT(A) int A##r, int A##c, complex* A##p+#define CMAT(A) int A##r, int A##c, doublecomplex* A##p+#define PMAT(A) int A##r, int A##c, void* A##p, int A##s +#define KFVEC(A) int A##n, const float*A##p #define KDVEC(A) int A##n, const double*A##p-#define KCVEC(A) int A##n, const double*A##p+#define KQVEC(A) int A##n, const complex*A##p+#define KCVEC(A) int A##n, const doublecomplex*A##p+#define KPVEC(A) int A##n, const void* A##p, int A##s+#define KFMAT(A) int A##r, int A##c, const float* A##p #define KDMAT(A) int A##r, int A##c, const double* A##p-#define KCMAT(A) int A##r, int A##c, const double* A##p+#define KQMAT(A) int A##r, int A##c, const complex* A##p+#define KCMAT(A) int A##r, int A##c, const doublecomplex* A##p+#define KPMAT(A) int A##r, int A##c, const void* A##p, int A##s  /********************************************************/ +int multiplyF(int ta, int tb, KFMAT(a),KFMAT(b),FMAT(r)); int multiplyR(int ta, int tb, KDMAT(a),KDMAT(b),DMAT(r)); int multiplyC(int ta, int tb, KCMAT(a),KCMAT(b),CMAT(r));+int multiplyQ(int ta, int tb, KQMAT(a),KQMAT(b),QMAT(r)); +int transF(KFMAT(x),FMAT(t)); int transR(KDMAT(x),DMAT(t));+int transQ(KQMAT(x),QMAT(t)); int transC(KCMAT(x),CMAT(t));+int transP(KPMAT(x),PMAT(t)); +int constantF(float * pval, FVEC(r)); int constantR(double * pval, DVEC(r));+int constantQ(complex* pval, QVEC(r)); int constantC(doublecomplex* pval, CVEC(r));+int constantP(void* pval, PVEC(r));++int float2double(FVEC(x),DVEC(y));+int double2float(DVEC(x),FVEC(y));++int conjugateQ(KQVEC(x),QVEC(t));+int conjugateC(KCVEC(x),CVEC(t));  int svd_l_R(KDMAT(x),DMAT(u),DVEC(s),DMAT(v)); int svd_l_Rdd(KDMAT(x),DMAT(u),DVEC(s),DMAT(v));

@@ -1,76 +0,0 @@-{-# LANGUAGE UndecidableInstances, MultiParamTypeClasses, FlexibleInstances #-}-------------------------------------------------------------------------------{- |-Module      :  Numeric.LinearAlgebra.Linear-Copyright   :  (c) Alberto Ruiz 2006-7-License     :  GPL-style--Maintainer  :  Alberto Ruiz (aruiz at um dot es)-Stability   :  provisional-Portability :  uses ffi--Basic optimized operations on vectors and matrices.---}--------------------------------------------------------------------------------module Numeric.LinearAlgebra.Linear (-    Linear(..)-) where--import Data.Packed.Vector-import Data.Packed.Matrix-import Data.Complex-import Numeric.GSL.Vector---- | Basic element-by-element functions.-class (Container c e) => Linear c e where-    -- | create a structure with a single element-    scalar      :: e -> c e-    scale       :: e -> c e -> c e-    -- | scale the element by element reciprocal of the object:-    ---    -- @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@-    scaleRecip  :: e -> c e -> c e-    addConstant :: e -> c e -> c e-    add         :: c e -> c e -> c e-    sub         :: c e -> c e -> c e-    -- | element by element multiplication-    mul         :: c e -> c e -> c e-    -- | element by element division-    divide      :: c e -> c e -> c e-    equal       :: c e -> c e -> Bool---instance Linear Vector Double where-    scale = vectorMapValR Scale-    scaleRecip = vectorMapValR Recip-    addConstant = vectorMapValR AddConstant-    add = vectorZipR Add-    sub = vectorZipR Sub-    mul = vectorZipR Mul-    divide = vectorZipR Div-    equal u v = dim u == dim v && vectorMax (vectorMapR Abs (sub u v)) == 0.0-    scalar x = fromList [x]--instance Linear Vector (Complex Double) where-    scale = vectorMapValC Scale-    scaleRecip = vectorMapValC Recip-    addConstant = vectorMapValC AddConstant-    add = vectorZipC Add-    sub = vectorZipC Sub-    mul = vectorZipC Mul-    divide = vectorZipC Div-    equal u v = dim u == dim v && vectorMax (mapVector magnitude (sub u v)) == 0.0-    scalar x = fromList [x]--instance (Linear Vector a, Container Matrix a) => (Linear Matrix a) where-    scale x = liftMatrix (scale x)-    scaleRecip x = liftMatrix (scaleRecip x)-    addConstant x = liftMatrix (addConstant x)-    add = liftMatrix2 add-    sub = liftMatrix2 sub-    mul = liftMatrix2 mul-    divide = liftMatrix2 divide-    equal a b = cols a == cols b && flatten a `equal` flatten b-    scalar x = (1><1) [x]

@@ -21,11 +21,12 @@ --, runBigTests ) where +import Data.Packed.Random import Numeric.LinearAlgebra import Numeric.LinearAlgebra.LAPACK import Numeric.LinearAlgebra.Tests.Instances import Numeric.LinearAlgebra.Tests.Properties-import Test.HUnit hiding ((~:),test,Testable)+import Test.HUnit hiding ((~:),test,Testable,State) import System.Info import Data.List(foldl1') import Numeric.GSL@@ -33,9 +34,14 @@ import qualified Prelude import System.CPUTime import Text.Printf+import Data.Packed.Development(unsafeFromForeignPtr,unsafeToForeignPtr)+import Control.Arrow((***))+import Debug.Trace  #include "Tests/quickCheckCompat.h" +debug x = trace (show x) x+ a ^ b = a Prelude.^ (b :: Int)  utest str b = TestCase $ assertBool str b@@ -44,6 +50,8 @@  feye n = flipud (ident n) :: Matrix Double +-----------------------------------------------------------+ detTest1 = det m == 26         && det mc == 38 :+ (-3)         && det (feye 2) == -1@@ -164,7 +172,7 @@     sol = fst $ fitModel 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0]      ok1 = and (zipWith f sols [5,0.1,1]) where f (x,d) r = abs (x-r)<2*d-    ok2 = pnorm PNorm2 (fromList (map fst sols) - fromList sol) < 1E-5+    ok2 = norm2 (fromList (map fst sols) - fromList sol) < 1E-5  ----------------------------------------------------- @@ -201,6 +209,155 @@     where fun n = foldl1' (<>) (map rot angles)               where angles = toList $ linspace n (0,1) +---------------------------------------------------------------------+-- vector <= 0.6.0.2 bug discovered by Patrick Perry+-- http://trac.haskell.org/vector/ticket/31++offsetTest = y == y' where+    x = fromList [0..3 :: Double]+    y = subVector 1 3 x+    (f,o,n) = unsafeToForeignPtr y+    y' = unsafeFromForeignPtr f o n++---------------------------------------------------------------------++normsVTest = TestList [+    utest "normv2CD" $ norm2PropC v+  , utest "normv2CF" $ norm2PropC (single v)+  , utest "normv2D"  $ norm2PropR x+  , utest "normv2F"  $ norm2PropR (single x)++  , utest "normv1CD" $ norm1 v          == 8+  , utest "normv1CF" $ norm1 (single v) == 8+  , utest "normv1D"  $ norm1 x          == 6+  , utest "normv1F"  $ norm1 (single x) == 6++  , utest "normvInfCD" $ normInf v          == 5+  , utest "normvInfCF" $ normInf (single v) == 5+  , utest "normvInfD"  $ normInf x          == 3+  , utest "normvInfF"  $ normInf (single x) == 3++ ] where v = fromList [1,-2,3:+4] :: Vector (Complex Double)+         x = fromList [1,2,-3] :: Vector Double+         norm2PropR a = norm2 a =~= sqrt (dot a a)+         norm2PropC a = norm2 a =~= realPart (sqrt (dot a (conj a)))+         a =~= b = fromList [a] |~| fromList [b]++normsMTest = TestList [+    utest "norm2mCD" $ pnorm PNorm2 v          =~= 8.86164970498005+  , utest "norm2mCF" $ pnorm PNorm2 (single v) =~= 8.86164970498005+  , utest "norm2mD"  $ pnorm PNorm2 x          =~= 5.96667765076216+  , utest "norm2mF"  $ pnorm PNorm2 (single x) =~= 5.96667765076216++  , utest "norm1mCD" $ pnorm PNorm1 v          == 9+  , utest "norm1mCF" $ pnorm PNorm1 (single v) == 9+  , utest "norm1mD"  $ pnorm PNorm1 x          == 7+  , utest "norm1mF"  $ pnorm PNorm1 (single x) == 7++  , utest "normmInfCD" $ pnorm Infinity v          == 12+  , utest "normmInfCF" $ pnorm Infinity (single v) == 12+  , utest "normmInfD"  $ pnorm Infinity x          == 8+  , utest "normmInfF"  $ pnorm Infinity (single x) == 8++  , utest "normmFroCD" $ pnorm Frobenius v          =~= 8.88819441731559+  , utest "normmFroCF" $ pnorm Frobenius (single v) =~~= 8.88819441731559+  , utest "normmFroD"  $ pnorm Frobenius x          =~= 6.24499799839840+  , utest "normmFroF"  $ pnorm Frobenius (single x) =~~= 6.24499799839840++ ] where v = (2><2) [1,-2*i,3:+4,7] :: Matrix (Complex Double)+         x = (2><2) [1,2,-3,5] :: Matrix Double+         a =~= b = fromList [a] :~10~: fromList [b]+         a =~~= b = fromList [a] :~5~: fromList [b]++---------------------------------------------------------------------++sumprodTest = TestList [+    utest "sumCD" $ sumElements z            == 6+  , utest "sumCF" $ sumElements (single z)   == 6+  , utest "sumD"  $ sumElements v            == 6+  , utest "sumF"  $ sumElements (single v)   == 6++  , utest "prodCD" $ prodProp z+  , utest "prodCF" $ prodProp (single z)+  , utest "prodD"  $ prodProp v+  , utest "prodF"  $ prodProp (single v)+ ] where v = fromList [1,2,3] :: Vector Double+         z = fromList [1,2-i,3+i]+         prodProp x = prodElements x == product (toList x)++---------------------------------------------------------------------++chainTest = utest "chain" $ foldl1' (<>) ms |~| optimiseMult ms where+    ms = [ diag (fromList [1,2,3 :: Double])+         , konst 3 (3,5)+         , (5><10) [1 .. ]+         , konst 5 (10,2)+         ]++---------------------------------------------------------------------++conjuTest m = mapVector conjugate (flatten (trans m)) == flatten (ctrans m)++---------------------------------------------------------------------++newtype State s a = State { runState :: s -> (a,s) }++instance Monad (State s) where+    return a = State $ \s -> (a,s)+    m >>= f = State $ \s -> let (a,s') = runState m s+                            in runState (f a) s'++state_get :: State s s+state_get = State $ \s -> (s,s)++state_put :: s -> State s ()+state_put s = State $ \_ -> ((),s)++evalState :: State s a -> s -> a+evalState m s = let (a,s') = runState m s+                in seq s' a++newtype MaybeT m a = MaybeT { runMaybeT :: m (Maybe a) }++instance Monad m => Monad (MaybeT m) where+    return a = MaybeT $ return $ Just a+    m >>= f  = MaybeT $ do+                        res <- runMaybeT m+                        case res of+                                 Nothing -> return Nothing+                                 Just r  -> runMaybeT (f r)+    fail _   = MaybeT $ return Nothing++lift_maybe m = MaybeT $ do+                        res <- m+                        return $ Just res++-- | apply a test to successive elements of a vector, evaluates to true iff test passes for all pairs+--successive_ :: Storable a => (a -> a -> Bool) -> Vector a -> Bool+successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ step (subVector 1 (dim v - 1) v))) (v @> 0)+   where step e = do+                  ep <- lift_maybe $ state_get+                  if t e ep+                     then lift_maybe $ state_put e+                     else (fail "successive_ test failed")++-- | operate on successive elements of a vector and return the resulting vector, whose length 1 less than that of the input+--successive :: (Storable a, Storable b) => (a -> a -> b) -> Vector a -> Vector b+successive f v = evalState (mapVectorM step (subVector 1 (dim v - 1) v)) (v @> 0)+   where step e = do+                  ep <- state_get+                  state_put e+                  return $ f ep e+++succTest = utest "successive" $+       successive_ (>) (fromList [1 :: Double,2,3,4]) == True+    && successive_ (>) (fromList [1 :: Double,3,2,4]) == False+    && successive (+) (fromList [1..10 :: Double]) == 9 |> [3,5,7,9,11,13,15,17,19]++---------------------------------------------------------------------++ -- | All tests must pass with a maximum dimension of about 20 --  (some tests may fail with bigger sizes due to precision loss). runTests :: Int  -- ^ maximum dimension@@ -208,14 +365,22 @@ runTests n = do     setErrorHandlerOff     let test p = qCheck n p-    putStrLn "------ mult"-    test (multProp1  . rConsist)-    test (multProp1  . cConsist)-    test (multProp2  . rConsist)-    test (multProp2  . cConsist)+    putStrLn "------ mult Double"+    test (multProp1 10 . rConsist)+    test (multProp1 10 . cConsist)+    test (multProp2 10 . rConsist)+    test (multProp2 10 . cConsist)+    putStrLn "------ mult Float"+    test (multProp1  6 . (single *** single) . rConsist)+    test (multProp1  6 . (single *** single) . cConsist)+    test (multProp2  6 . (single *** single) . rConsist)+    test (multProp2  6 . (single *** single) . cConsist)     putStrLn "------ sub-trans"     test (subProp . rM)     test (subProp . cM)+    putStrLn "------ ctrans"+    test (conjuTest . cM)+    test (conjuTest . zM)     putStrLn "------ lu"     test (luProp    . rM)     test (luProp    . cM)@@ -286,6 +451,9 @@     test (qrProp     . cM)     test (rqProp     . rM)     test (rqProp     . cM)+    test (rqProp1     . cM)+    test (rqProp2     . cM)+    test (rqProp3     . cM)     putStrLn "------ hess"     test (hessProp   . rSq)     test (hessProp   . cSq)@@ -296,21 +464,31 @@     test (cholProp   . rPosDef)     test (cholProp   . cPosDef)     putStrLn "------ expm"-    test (expmDiagProp . rSqWC)+    test (expmDiagProp . complex. rSqWC)     test (expmDiagProp . cSqWC)     putStrLn "------ fft"     test (\v -> ifft (fft v) |~| v)-    putStrLn "------ vector operations"+    putStrLn "------ vector operations - Double"     test (\u -> sin u ^ 2 + cos u ^ 2 |~| (1::RM))     test $ (\u -> sin u ^ 2 + cos u ^ 2 |~| (1::CM)) . liftMatrix makeUnitary     test (\u -> sin u ** 2 + cos u ** 2 |~| (1::RM))     test (\u -> cos u * tan u |~| sin (u::RM))     test $ (\u -> cos u * tan u |~| sin (u::CM)) . liftMatrix makeUnitary+    putStrLn "------ vector operations - Float"+    test (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::FM))+    test $ (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::ZM)) . liftMatrix makeUnitary+    test (\u -> sin u ** 2 + cos u ** 2 |~~| (1::FM))+    test (\u -> cos u * tan u |~~| sin (u::FM))+    test $ (\u -> cos u * tan u |~~| sin (u::ZM)) . liftMatrix makeUnitary     putStrLn "------ read . show"     test (\m -> (m::RM) == read (show m))     test (\m -> (m::CM) == read (show m))     test (\m -> toRows (m::RM) == read (show (toRows m)))     test (\m -> toRows (m::CM) == read (show (toRows m)))+    test (\m -> (m::FM) == read (show m))+    test (\m -> (m::ZM) == read (show m))+    test (\m -> toRows (m::FM) == read (show (toRows m)))+    test (\m -> toRows (m::ZM) == read (show (toRows m)))     putStrLn "------ some unit tests"     _ <- runTestTT $ TestList         [ utest "1E5 rots" rotTest@@ -332,7 +510,7 @@         , utest "randomGaussian" randomTestGaussian         , utest "randomUniform" randomTestUniform         , utest "buildVector/Matrix" $-                        comp (10 |> [0::Double ..]) == buildVector 10 fromIntegral+                        complex (10 |> [0::Double ..]) == buildVector 10 fromIntegral                      && ident 5 == buildMatrix 5 5 (\(r,c) -> if r==c then 1::Double else 0)         , utest "rank" $  rank ((2><3)[1,0,0,1,6*eps,0]) == 1                        && rank ((2><3)[1,0,0,1,7*eps,0]) == 2@@ -340,9 +518,20 @@         , odeTest         , fittingTest         , mbCholTest+        , utest "offset" offsetTest+        , normsVTest+        , normsMTest+        , sumprodTest+        , chainTest+        , succTest         ]     return () ++-- single precision approximate equality+infixl 4 |~~|+a |~~| b = a :~6~: b+ makeUnitary v | realPart n > 1    = v / scalar n               | otherwise = v     where n = sqrt (conj v <.> v)@@ -356,6 +545,7 @@ -- | Performance measurements. runBenchmarks :: IO () runBenchmarks = do+  --cholBench     solveBench     subBench     multBench@@ -455,3 +645,18 @@     solveBenchN 500     solveBenchN 1000     -- solveBenchN 1500++--------------------------------++cholBenchN n = do+    let x = uniformSample 777 (2*n) (replicate n (-1,1))+        a = trans x <> x+    a `seq` putStrLn ""+    time ("chol " ++ show n) (chol a)++cholBench = do+    cholBenchN 1200+    cholBenchN 600+    cholBenchN 300+--    cholBenchN 150+--    cholBenchN 50

@@ -1,4 +1,4 @@-{-# LANGUAGE FlexibleContexts, UndecidableInstances, CPP #-}+{-# LANGUAGE FlexibleContexts, UndecidableInstances, CPP, FlexibleInstances #-} {-# OPTIONS_GHC -fno-warn-unused-imports #-} ----------------------------------------------------------------------------- {- |@@ -22,17 +22,16 @@     SqWC(..),   rSqWC, cSqWC,     PosDef(..), rPosDef, cPosDef,     Consistent(..), rConsist, cConsist,-    RM,CM, rM,cM+    RM,CM, rM,cM,+    FM,ZM, fM,zM ) where --+import System.Random  import Numeric.LinearAlgebra import Control.Monad(replicateM) #include "quickCheckCompat.h" - #if MIN_VERSION_QuickCheck(2,0,0) shrinkListElementwise :: (Arbitrary a) => [a] -> [[a]] shrinkListElementwise []     = []@@ -43,7 +42,8 @@ shrinkPair (a,b) = [ (a,x) | x <- shrink b ] ++ [ (x,b) | x <- shrink a ] #endif -+#if MIN_VERSION_QuickCheck(2,1,1)+#else instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where     arbitrary = do         re <- arbitrary@@ -58,6 +58,8 @@     coarbitrary = undefined  #endif +#endif+ chooseDim = sized $ \m -> choose (1,max 1 m)  instance (Field a, Arbitrary a) => Arbitrary (Vector a) where @@ -68,7 +70,7 @@ #if MIN_VERSION_QuickCheck(2,0,0)     -- shrink any one of the components     shrink = map fromList . shrinkListElementwise . toList-                              + #else     coarbitrary = undefined #endif@@ -133,10 +135,14 @@     coarbitrary = undefined #endif +class (Field a, Arbitrary a, Element (RealOf a), Random (RealOf a)) => ArbitraryField a+instance ArbitraryField Double+instance ArbitraryField (Complex Double) + -- a well-conditioned general matrix (the singular values are between 1 and 100) newtype (WC a) = WC (Matrix a) deriving Show-instance (Field a, Arbitrary a) => Arbitrary (WC a) where+instance (ArbitraryField a) => Arbitrary (WC a) where     arbitrary = do         m <- arbitrary         let (u,_,v) = svd m@@ -144,7 +150,7 @@             c = cols m             n = min r c         sv' <- replicateM n (choose (1,100))-        let s = diagRect (fromList sv') r c+        let s = diagRect 0 (fromList sv') r c         return $ WC (u <> real s <> trans v)  #if MIN_VERSION_QuickCheck(2,0,0)@@ -155,7 +161,7 @@  -- a well-conditioned square matrix (the singular values are between 1 and 100) newtype (SqWC a) = SqWC (Matrix a) deriving Show-instance (Field a, Arbitrary a) => Arbitrary (SqWC a) where+instance (ArbitraryField a) => Arbitrary (SqWC a) where     arbitrary = do         Sq m <- arbitrary         let (u,_,v) = svd m@@ -172,7 +178,8 @@  -- a positive definite square matrix (the eigenvalues are between 0 and 100) newtype (PosDef a) = PosDef (Matrix a) deriving Show-instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (PosDef a) where+instance (ArbitraryField a, Num (Vector a)) +    => Arbitrary (PosDef a) where     arbitrary = do         Her m <- arbitrary         let (_,v) = eigSH m@@ -209,10 +216,16 @@  type RM = Matrix Double type CM = Matrix (Complex Double)+type FM = Matrix Float+type ZM = Matrix (Complex Float) + rM m = m :: RM cM m = m :: CM+fM m = m :: FM+zM m = m :: ZM + rHer (Her m) = m :: RM cHer (Her m) = m :: CM @@ -233,3 +246,4 @@  rConsist (Consistent (a,b)) = (a,b::RM) cConsist (Consistent (a,b)) = (a,b::CM)+

@@ -32,7 +32,7 @@     svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4,     svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7,     eigProp, eigSHProp, eigProp2, eigSHProp2,-    qrProp, rqProp,+    qrProp, rqProp, rqProp1, rqProp2, rqProp3,     hessProp,     schurProp1, schurProp2,     cholProp,@@ -42,24 +42,27 @@     linearSolveProp, linearSolveProp2 ) where -import Numeric.LinearAlgebra+import Numeric.LinearAlgebra --hiding (real,complex) import Numeric.LinearAlgebra.LAPACK import Debug.Trace #include "quickCheckCompat.h"  +--real x = real'' x+--complex x = complex'' x+ debug x = trace (show x) x  -- relative error-dist :: (Normed t, Num t) => t -> t -> Double-dist a b = r+dist :: (Normed c t, Num (c t)) => c t -> c t -> Double+dist a b = realToFrac r     where norm = pnorm Infinity           na = norm a           nb = norm b           nab = norm (a-b)           mx = max na nb           mn = min na nb-          r = if mn < eps+          r = if mn < peps                 then mx                 else nab/mx @@ -68,7 +71,7 @@ --a |~| b = dist a b < 10^^(-10)  data Aprox a = (:~) a Int-(~:) :: (Normed a, Num a) => Aprox a -> a -> Bool+-- (~:) :: (Normed a, Num a) => Aprox a -> a -> Bool a :~n~: b = dist a b < 10^^(-n)  ------------------------------------------------------@@ -135,7 +138,7 @@  svdProp1a svdfun m = m |~| u <> real d <> trans v && unitary u && unitary v where     (u,s,v) = svdfun m-    d = diagRect s (rows m) (cols m)+    d = diagRect 0 s (rows m) (cols m)  svdProp1b svdfun m = unitary u && unitary v where     (u,_,v) = svdfun m@@ -207,16 +210,22 @@ qrProp m = q <> r |~| m && unitary q && upperTriang r     where (q,r) = qr m -rqProp m = r <> q |~| m && unitary q && utr+rqProp m = r <> q |~| m && unitary q && upperTriang' r     where (r,q) = rq m-          upptr f c = buildMatrix f c $ \(r',c') -> if r'-t > c' then 0 else 1-              where t = f-c-          utr = upptr (rows r) (cols r) * r |~| r -upperTriang' m = rows m == 1 || down |~| z-    where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m))-          z = constant 0 (dim down)+rqProp1 m = r <> q |~| m+    where (r,q) = rq m +rqProp2 m = unitary q+    where (r,q) = rq m++rqProp3 m = upperTriang' r+    where (r,q) = rq m++upperTriang' r = upptr (rows r) (cols r) * r |~| r+    where upptr f c = buildMatrix f c $ \(r',c') -> if r'-t > c' then 0 else 1+              where t = f-c+ hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h     where (p,h) = hess m @@ -237,9 +246,9 @@ mulH a b = fromLists [[ doth ai bj | bj <- toColumns b] | ai <- toRows a ]     where doth u v = sum $ zipWith (*) (toList u) (toList v) -multProp1 (a,b) = a <> b |~| mulH a b+multProp1 p (a,b) = (a <> b) :~p~: (mulH a b) -multProp2 (a,b) = ctrans (a <> b) |~| ctrans b <> ctrans a+multProp2 p (a,b) = (ctrans (a <> b)) :~p~: (ctrans b <> ctrans a)  linearSolveProp f m = f m m |~| ident (rows m) 

@@ -0,0 +1,72 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FunctionalDependencies #-}++-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.Matrix+-- Copyright   :  (c) Alberto Ruiz 2010+-- License     :  GPL-style+--+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>+-- Stability   :  provisional+-- Portability :  portable+--+-- Provides instances of standard classes 'Show', 'Read', 'Eq',+-- 'Num', 'Fractional', and 'Floating' for 'Matrix'.+--+-- In arithmetic operations one-component+-- vectors and matrices automatically expand to match the dimensions of the other operand.++-----------------------------------------------------------------------------++module Numeric.Matrix (+                      ) where++-------------------------------------------------------------------++import Numeric.Container++-------------------------------------------------------------------++instance Container Matrix a => Eq (Matrix a) where+    (==) = equal++instance (Container Matrix a, Num (Vector a)) => Num (Matrix a) where+    (+) = liftMatrix2Auto (+)+    (-) = liftMatrix2Auto (-)+    negate = liftMatrix negate+    (*) = liftMatrix2Auto (*)+    signum = liftMatrix signum+    abs = liftMatrix abs+    fromInteger = (1><1) . return . fromInteger++---------------------------------------------------++instance (Container Vector a, Fractional (Vector a), Num (Matrix a)) => Fractional (Matrix a) where+    fromRational n = (1><1) [fromRational n]+    (/) = liftMatrix2Auto (/)++---------------------------------------------------------++instance (Floating a, Container Vector a, Floating (Vector a), Fractional (Matrix a)) => Floating (Matrix a) where+    sin   = liftMatrix sin+    cos   = liftMatrix cos+    tan   = liftMatrix tan+    asin  = liftMatrix asin+    acos  = liftMatrix acos+    atan  = liftMatrix atan+    sinh  = liftMatrix sinh+    cosh  = liftMatrix cosh+    tanh  = liftMatrix tanh+    asinh = liftMatrix asinh+    acosh = liftMatrix acosh+    atanh = liftMatrix atanh+    exp   = liftMatrix exp+    log   = liftMatrix log+    (**)  = liftMatrix2Auto (**)+    sqrt  = liftMatrix sqrt+    pi    = (1><1) [pi]

@@ -0,0 +1,223 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.Vector+-- Copyright   :  (c) Alberto Ruiz 2010+-- License     :  GPL-style+--+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>+-- Stability   :  provisional+-- Portability :  portable+--+-- Provides instances of standard classes 'Show', 'Read', 'Eq',+-- 'Num', 'Fractional',  and 'Floating' for 'Vector'.+-- +-----------------------------------------------------------------------------++module Numeric.Vector (+                      ) where++import Numeric.GSL.Vector+import Numeric.Container++-------------------------------------------------------------------++#ifndef VECTOR+import Foreign(Storable)+#endif++------------------------------------------------------------------++#ifndef VECTOR++instance (Show a, Storable a) => (Show (Vector a)) where+    show v = (show (dim v))++" |> " ++ show (toList v)++#endif++#ifdef VECTOR++instance (Element a, Read a) => Read (Vector a) where+    readsPrec _ s = [(fromList . read $ listnums, rest)]+        where (thing,trest) = breakAt ']' s+              (dims,listnums) = breakAt ' ' (dropWhile (==' ') thing)+              rest = drop 31 trest+#else++instance (Element a, Read a) => Read (Vector a) where+    readsPrec _ s = [((d |>) . read $ listnums, rest)]+        where (thing,rest) = breakAt ']' s+              (dims,listnums) = breakAt '>' thing+              d = read . init . fst . breakAt '|' $ dims++#endif++breakAt c l = (a++[c],tail b) where+    (a,b) = break (==c) l+++------------------------------------------------------------------++adaptScalar f1 f2 f3 x y+    | dim x == 1 = f1   (x@>0) y+    | dim y == 1 = f3 x (y@>0)+    | otherwise = f2 x y++------------------------------------------------------------------++#ifndef VECTOR++instance Container Vector a => Eq (Vector a) where+    (==) = equal++#endif++instance Num (Vector Float) where+    (+) = adaptScalar addConstant add (flip addConstant)+    negate = scale (-1)+    (*) = adaptScalar scale mul (flip scale)+    signum = vectorMapF Sign+    abs = vectorMapF Abs+    fromInteger = fromList . return . fromInteger++instance Num (Vector Double) where+    (+) = adaptScalar addConstant add (flip addConstant)+    negate = scale (-1)+    (*) = adaptScalar scale mul (flip scale)+    signum = vectorMapR Sign+    abs = vectorMapR Abs+    fromInteger = fromList . return . fromInteger++instance Num (Vector (Complex Double)) where+    (+) = adaptScalar addConstant add (flip addConstant)+    negate = scale (-1)+    (*) = adaptScalar scale mul (flip scale)+    signum = vectorMapC Sign+    abs = vectorMapC Abs+    fromInteger = fromList . return . fromInteger++instance Num (Vector (Complex Float)) where+    (+) = adaptScalar addConstant add (flip addConstant)+    negate = scale (-1)+    (*) = adaptScalar scale mul (flip scale)+    signum = vectorMapQ Sign+    abs = vectorMapQ Abs+    fromInteger = fromList . return . fromInteger++---------------------------------------------------++instance (Container Vector a, Num (Vector a)) => Fractional (Vector a) where+    fromRational n = fromList [fromRational n]+    (/) = adaptScalar f divide g where+        r `f` v = scaleRecip r v+        v `g` r = scale (recip r) v++-------------------------------------------------------++instance Floating (Vector Float) where+    sin   = vectorMapF Sin+    cos   = vectorMapF Cos+    tan   = vectorMapF Tan+    asin  = vectorMapF ASin+    acos  = vectorMapF ACos+    atan  = vectorMapF ATan+    sinh  = vectorMapF Sinh+    cosh  = vectorMapF Cosh+    tanh  = vectorMapF Tanh+    asinh = vectorMapF ASinh+    acosh = vectorMapF ACosh+    atanh = vectorMapF ATanh+    exp   = vectorMapF Exp+    log   = vectorMapF Log+    sqrt  = vectorMapF Sqrt+    (**)  = adaptScalar (vectorMapValF PowSV) (vectorZipF Pow) (flip (vectorMapValF PowVS))+    pi    = fromList [pi]++-------------------------------------------------------------++instance Floating (Vector Double) where+    sin   = vectorMapR Sin+    cos   = vectorMapR Cos+    tan   = vectorMapR Tan+    asin  = vectorMapR ASin+    acos  = vectorMapR ACos+    atan  = vectorMapR ATan+    sinh  = vectorMapR Sinh+    cosh  = vectorMapR Cosh+    tanh  = vectorMapR Tanh+    asinh = vectorMapR ASinh+    acosh = vectorMapR ACosh+    atanh = vectorMapR ATanh+    exp   = vectorMapR Exp+    log   = vectorMapR Log+    sqrt  = vectorMapR Sqrt+    (**)  = adaptScalar (vectorMapValR PowSV) (vectorZipR Pow) (flip (vectorMapValR PowVS))+    pi    = fromList [pi]++-------------------------------------------------------------++instance Floating (Vector (Complex Double)) where+    sin   = vectorMapC Sin+    cos   = vectorMapC Cos+    tan   = vectorMapC Tan+    asin  = vectorMapC ASin+    acos  = vectorMapC ACos+    atan  = vectorMapC ATan+    sinh  = vectorMapC Sinh+    cosh  = vectorMapC Cosh+    tanh  = vectorMapC Tanh+    asinh = vectorMapC ASinh+    acosh = vectorMapC ACosh+    atanh = vectorMapC ATanh+    exp   = vectorMapC Exp+    log   = vectorMapC Log+    sqrt  = vectorMapC Sqrt+    (**)  = adaptScalar (vectorMapValC PowSV) (vectorZipC Pow) (flip (vectorMapValC PowVS))+    pi    = fromList [pi]++-----------------------------------------------------------++instance Floating (Vector (Complex Float)) where+    sin   = vectorMapQ Sin+    cos   = vectorMapQ Cos+    tan   = vectorMapQ Tan+    asin  = vectorMapQ ASin+    acos  = vectorMapQ ACos+    atan  = vectorMapQ ATan+    sinh  = vectorMapQ Sinh+    cosh  = vectorMapQ Cosh+    tanh  = vectorMapQ Tanh+    asinh = vectorMapQ ASinh+    acosh = vectorMapQ ACosh+    atanh = vectorMapQ ATanh+    exp   = vectorMapQ Exp+    log   = vectorMapQ Log+    sqrt  = vectorMapQ Sqrt+    (**)  = adaptScalar (vectorMapValQ PowSV) (vectorZipQ Pow) (flip (vectorMapValQ PowVS))+    pi    = fromList [pi]++-----------------------------------------------------------+++-- instance (Storable a, Num (Vector a)) => Monoid (Vector a) where+--     mempty = 0 { idim = 0 }+--     mappend a b = mconcat [a,b]+--     mconcat = j . filter ((>0).dim)+--         where j [] = mempty+--               j l  = join l++---------------------------------------------------------------++-- instance (NFData a, Storable a) => NFData (Vector a) where+--     rnf = rnf . (@>0)+--+-- instance (NFData a, Element a) => NFData (Matrix a) where+--     rnf = rnf . flatten+++--------------------------------------------------------------------------+