packages feed

hmatrix 0.5.0.1 → 0.5.1.1

raw patch · 28 files changed

+660/−645 lines, 28 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

+ Data.Packed.ST: newUndefinedMatrix :: (Element t) => MatrixOrder -> Int -> Int -> ST s (STMatrix s t)
+ Data.Packed.ST: newUndefinedVector :: (Element t) => Int -> ST s (STVector s t)
+ Data.Packed.Vector: foldLoop :: (Int -> t -> t) -> t -> Int -> t
+ Data.Packed.Vector: foldVector :: (Double -> b -> b) -> b -> Vector Double -> b
+ Data.Packed.Vector: foldVectorG :: (Storable a) => (Int -> (Int -> a) -> t -> t) -> t -> Vector a -> t
+ Data.Packed.Vector: mapVector :: (Storable a, Storable b) => (a -> b) -> Vector a -> Vector b
+ Data.Packed.Vector: zipVector :: (Storable a, Storable b, Storable c) => (a -> b -> c) -> Vector a -> Vector b -> Vector c
+ Numeric.GSL.Minimization: minimizeVectorBFGS2 :: Double -> Double -> Double -> Int -> ([Double] -> Double) -> ([Double] -> [Double]) -> [Double] -> ([Double], Matrix Double)

Files

README view
@@ -103,3 +103,11 @@   avoid the wrong NaNs produced by foreign functions.  - Reiner Pope added support for luSolve, based on (d|z)getrs.++- Simon Beaumont reported the need of QuickCheck<2 and the invalid+  asm("finit") on ppc. He also contributed the configuration options+  for the accelerate framework on OS X.++- Daniel Schüssler added compatibility with QuickCheck 2 as well+  as QuickCheck 1 using the C preprocessor. He also added some+  implementations for the new "shrink" method of class Arbitrary.
examples/benchmarks.hs view
@@ -18,19 +18,48 @@  -------------------------------------------------------------------------------- -main = sequence_ [bench1,bench2,bench3,bench4,bench5 1000000 3]+main = sequence_ [bench1,+                  bench2,+                  bench4,+                  bench5 1000000 3, bench5 100000 50,+                  bench6 100 (100000::Double), bench6 100000 (100::Double), bench6 10000 (1000::Double)]  w :: Vector Double w = constant 1 5000000 w2 = 1 * w +v = flatten $ ident 500 :: Vector Double++ bench1 = do+    time $ print$ vectorMax (w+w2) -- evaluate it     putStrLn "Sum of a vector with 5M doubles:"     print$ vectorMax (w+w2) -- evaluate it     time $ printf "         BLAS: %.2f: " $ sumVB w     time $ printf "BLAS only dot: %.2f: " $ w <.> w2     time $ printf "      Haskell: %.2f: " $ sumVH w     time $ printf "       innerH: %.2f: " $ innerH w w2+    time $ printf "foldVector: %.2f: " $ sumVector w+    let getPos k s = if k `mod` 500 < 200 && w@>k > 0 then k:s else s+    putStrLn "foldLoop for element selection:"+    time $ print $ (`divMod` 500) $ maximum $ foldLoop getPos [] (dim w)+    putStrLn "constant 5M:"+    time $ print $ constant (1::Double) 5000001 @> 7+    time $ print $ constant           i 5000001 @> 7+    time $ print $ conj (constant i 5000001) @> 7+    putStrLn "zips C vs H:"+    time $ print $ (w / w2) @> 7+    time $ print $ (zipVector (/) w w2) @> 7+    putStrLn "folds C/BLAS vs H:"+    let t = constant (1::Double) 5000002+    print $ t @> 7+    time $ print $ foldVector max (t@>0) t+    time $ print $ vectorMax t+    time $ print $ sqrt $ foldVector (\v s -> v*v+s) 0 t+    time $ print $ pnorm PNorm2 t+    putStrLn "scale C/BLAS vs H:"+    time $ print $ mapVector (*2) t @> 7+    time $ print $ (2 * t) @> 7  sumVB v = constant 1 (dim v) <.> v @@ -43,11 +72,15 @@  innerH u v = go (d - 1) 0      where-       d = dim u+       d = min (dim u) (dim v)        go :: Int -> Double -> Double        go 0 s = s + (u @> 0) * (v @> 0)        go !j !s = go (j - 1) (s + (u @> j) * (v @> j)) ++-- sumVector = foldVectorG (\k v s -> v k + s) 0.0+sumVector = foldVector (+) 0.0+ --------------------------------------------------------------------------------  bench2 = do@@ -83,39 +116,6 @@  -------------------------------------------------------------------------------- -bench3 = do-    putStrLn "-------------------------------------------------------"-    putStrLn "foldVector"-    let v = flatten $ ident 500 :: Vector Double-    print $ vectorMax v  -- evaluate it--    putStrLn "sum, dim=5M:"-    -- time $ print $ foldLoop (\k s -> w@>k + s) 0.0 (dim w)-    time $ print $ sumVector w--    putStrLn "sum, dim=0.25M:"-    --time $ print $ foldLoop (\k s -> v@>k + s) 0.0 (dim v)-    time $ print $ sumVector v--    let getPos k s = if k `mod` 500 < 200 && v@>k > 0 then k:s else s-    putStrLn "foldLoop for element selection, dim=0.25M:"-    time $ print $ (`divMod` 500) $ maximum $ foldLoop getPos [] (dim v)--foldLoop f s d = go (d - 1) s-     where-       go 0 s = f (0::Int) s-       go !j !s = go (j - 1) (f j s)--foldVector f s v = foldLoop g s (dim v)-    where g !k !s = f k (v@>) s-          {-# INLINE g #-} -- Thanks Ryan Ingram (http://permalink.gmane.org/gmane.comp.lang.haskell.cafe/46479)--sumVector = foldVector (\k v s -> v k + s) 0.0---- foldVector is slower if used in two places unless we use the above INLINE--- this does not happen with foldLoop---------------------------------------------------------------------------------- bench4 = do     putStrLn "-------------------------------------------------------"     putStrLn "1000x1000 inverse"@@ -127,19 +127,25 @@ --------------------------------------------------------------------------------  op1 a b = a <> trans b--op2 a b = a + trans b+op2 a b = a  + trans b  timep = time . print . vectorMax . flatten  bench5 n d = do     putStrLn "-------------------------------------------------------"-    putStrLn "transpose in multiply"+    putStrLn "transpose in add"     let ms = replicate n ((ident d :: Matrix Double))-    let mz = replicate n (diag (constant (0::Double) d))+    timep $ foldl1' (+)  ms+    timep $ foldl1' op2  ms+    putStrLn "-------------------------------------------------------"+    putStrLn "transpose in multiply"+     timep $ foldl1' (<>) ms     timep $ foldl1' op1  ms++--------------------------------------------------------------------------------++bench6 sz n = do     putStrLn "-------------------------------------------------------"-    putStrLn "transpose in add"-    timep $ foldl1' (+)  ms-    timep $ foldl1' op2  ms+    putStrLn "many constants"+    time $ print $ sum $ map ((@>0). flip constant sz) [1..n]
− examples/experiments/Static.hs 
@@ -1,115 +0,0 @@-{-# OPTIONS_GHC -fglasgow-exts -fth -fallow-overlapping-instances -fallow-undecidable-instances #-}--module Static where--import Language.Haskell.TH-import Numeric.LinearAlgebra-import Foreign-import Language.Haskell.TH.Syntax--instance Lift Double where-  lift x = return (LitE (RationalL (toRational x)))--instance Lift (Vector a ) where-    lift v = [e| v |]--instance Lift (Matrix a) where-    lift m = [e| m |]--tdim :: Int -> ExpQ-tdim 0 = [| Z |]-tdim n = [| S $(tdim (n-1)) |]---data Z = Z deriving Show-data S a = S a deriving Show--class Dim a--instance Dim Z-instance Dim a => Dim (S a)--class Sum a b c | a b -> c -- , a c -> b, b c -> a--instance Sum Z a a-instance Sum a Z a-instance Sum a b c => Sum a (S b) (S c)--newtype SVec d t = SVec (Vector t) deriving Show-newtype SMat r c t = SMat (Matrix t) deriving Show--createl :: d -> [Double] -> SVec d Double-createl d l = SVec (fromList l)--createv :: Storable t => d -> Vector t -> SVec d t-createv d v = SVec v--vec'' v = [|createv ($(tdim (dim v))) v|]--vec' :: [Double] -> ExpQ-vec' d = [| createl ($(tdim (length d))) d |]---createm :: (Dim r, Dim c) => r -> c -> (Matrix Double) -> SMat r c Double-createm _ _ m = SMat m--createml :: (Dim r, Dim c) => r -> c -> Int -> Int -> [Double] -> SMat r c Double-createml _ _ r c l = SMat ((r><c) l)--mat :: Int -> Int -> [Double] -> ExpQ-mat r c l = [| createml ($(tdim r)) ($(tdim c)) r c l |]--vec :: [Double] -> ExpQ-vec d = mat (length d) 1 d-----mat' :: Matrix Double -> ExpQ---mat' m = [| createm ($(tdim (rows m))) ($(tdim (cols m))) m |]--covec :: [Double] -> ExpQ-covec d = mat 1 (length d) d--scalar :: SMat (S Z) (S Z) Double -> Double-scalar (SMat m) = flatten m @> 0--v = fromList [1..5] :: Vector Double-l = [1,1.5..5::Double]--k = [11..30::Int]--rawv (SVec v) = v-raw (SMat m) = m--liftStatic :: (Matrix a -> Matrix b -> Matrix c) -> SMat dr dc a -> SMat dr dc b -> SMat dr dc c-liftStatic f a b = SMat (f (raw a) (raw b))--a |+| b = liftStatic (+) a b--prod :: SMat r k Double -> SMat k c Double -> SMat r c Double-prod a b = SMat (raw a <> raw b)--strans :: SMat r c Double -> SMat c r Double-strans = SMat . trans . raw--sdot a b = scalar (prod a b)--jv :: (Field t, Sum r1 r2 r3) => SMat r1 c t -> SMat r2 c t -> SMat r3 c t-jv a b = SMat ((raw a) <-> (raw b))---- curiously, we cannot easily fold jv because the matrics are not of the same type.--jh a b = strans (jv (strans a) (strans b))---homog :: (Field t) => SMat r c t -> SMat (S r) c t-homog m = SMat (raw m <-> constant 1 (cols (raw m)))--inhomog :: (Linear Vector t) => SMat (S (S r)) c t -> SMat r c t-inhomog (SMat m) = SMat (sm <> d) where-    sm = takeRows r' m-    d = diag $ 1 / (flatten $ dropRows r' m)-    r' = rows m -1---ht t vs = inhomog (t `prod` homog vs)-
− examples/experiments/useStatic.hs 
@@ -1,36 +0,0 @@-{-# OPTIONS -fno-monomorphism-restriction #-}--import Static-import Numeric.LinearAlgebra---x = $(vec [1,2])--y = $(vec [5,7])--z a = vec [a,a]--w = $(vec [1,2,3])--cx = $(covec [1,2,3])---t3 = $(tdim 3)--crm33 = createml t3 t3 3 3--rot a = crm33 [a,0,0,0,a,0,0,0,1]----q = x |+| y |+| $(z 5)--m = $(mat 2 3 [1..6])--n = $(mat 3 5 [1..15])--infixl 7 <*>-(<*>) = prod--r1 = m <*> n-r2 = strans (strans n <*> strans m)----r' = prod n m
examples/minimize.hs view
@@ -13,6 +13,9 @@ -- the conjugate gradient method minimizeCG = minimizeConjugateGradient 1E-2 1E-4 1E-3 30 +-- the BFGS2 method+minimizeBFGS2 = minimizeVectorBFGS2 1E-2 1E-2 1E-3 30+ -- a minimization algorithm which does not require the gradient minimizeS f xi = minimizeNMSimplex f xi (replicate (length xi) 1) 1E-2 100 @@ -24,20 +27,27 @@     (a,_:b) = splitAt n v  main = do-    -- conjugate gradient with true gradient-    let (s,p) = minimizeCG f df [5,7]+    putStrLn "BFGS2 with true gradient"+    let (s,p) = minimizeBFGS2 f df [5,7]     print s -- solution     disp p -- evolution of the algorithm     let [x,y] = drop 2 (toColumns p)     mplot [x,y]  -- path from the starting point to the solution -    -- conjugate gradient with estimated gradient+    putStrLn "conjugate gradient with true gradient"+    let (s,p) = minimizeCG f df [5,7]+    print s+    disp p+    let [x,y] = drop 2 (toColumns p)+    mplot [x,y]++    putStrLn "conjugate gradient with estimated gradient"     let (s,p) = minimizeCG f (gradient f) [5,7]     print s     disp p     mplot $ drop 2 (toColumns p) -    -- without gradient, using the NM Simplex method+    putStrLn "without gradient, using the NM Simplex method"     let (s,p) = minimizeS f [5,7]     print s     disp p
hmatrix.cabal view
@@ -1,5 +1,5 @@ Name:               hmatrix-Version:            0.5.0.1+Version:            0.5.1.1 License:            GPL License-file:       LICENSE Author:             Alberto Ruiz@@ -11,11 +11,13 @@                     and other numerical computations, internally implemented using                     GSL, BLAS and LAPACK. Category:           Math-tested-with:        GHC ==6.10.0+tested-with:        GHC ==6.10.2  cabal-version:      >=1.2 build-type:         Simple+extra-source-files: lib/Numeric/LinearAlgebra/Tests/quickCheckCompat.h + flag splitBase     description:    Choose the new smaller, split-up base package. @@ -23,6 +25,10 @@     description:    Link with Intel's MKL optimized libraries.     default:        False +flag accelerate+    description:    Use the accelerate framework for LAPACK/BLAS on OS X+    default:        False+ flag unsafe     description: Compile the library with bound checking disabled.     default: False@@ -93,8 +99,7 @@                         Numeric.GSL.Special.Internal,                         Numeric.LinearAlgebra.Tests.Instances,                         Numeric.LinearAlgebra.Tests.Properties-    C-sources:          lib/Data/Packed/Internal/auxi.c,-                        lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c,+    C-sources:          lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c,                         lib/Numeric/GSL/gsl-aux.c      ghc-prof-options:   -auto-all@@ -104,6 +109,8 @@     if flag(unsafe)         cpp-options: -DUNSAFE +    if impl(ghc < 6.10.2)+        cpp-options: -DFINIT      if flag(mkl)       if arch(x86_64)@@ -111,6 +118,10 @@       else         extra-libraries:   gsl mkl_lapack mkl_intel mkl_sequential mkl_core     else+      if flag(accelerate)+        frameworks: Accelerate+        extra-libraries: gsl+      else          extra-libraries: gsl lapack
lib/Data/Packed.hs view
@@ -24,7 +24,7 @@ import Data.Packed.Vector import Data.Packed.Matrix import Data.Complex-import Data.Packed.Internal(fromComplex,toComplex,comp,conj)+import Data.Packed.Internal(fromComplex,toComplex,conj)  -- | conversion utilities class (Element e) => Container c e where@@ -38,7 +38,7 @@ instance Container Vector Double where     toComplex = Data.Packed.Internal.toComplex     fromComplex = Data.Packed.Internal.fromComplex-    comp = Data.Packed.Internal.comp+    comp = internalcomp     conj = Data.Packed.Internal.conj     real = id     complex = Data.Packed.comp@@ -56,7 +56,7 @@     fromComplex z = (reshape c r, reshape c i)         where (r,i) = Data.Packed.fromComplex (flatten z)               c = cols z-    comp = liftMatrix Data.Packed.Internal.comp+    comp = liftMatrix internalcomp     conj = liftMatrix Data.Packed.Internal.conj     real = id     complex = Data.Packed.comp@@ -68,3 +68,8 @@     conj = undefined     real = Data.Packed.comp     complex = id+++-- | converts a real vector into a complex representation (with zero imaginary parts)+internalcomp :: Vector Double -> Vector (Complex Double)+internalcomp v = Data.Packed.Internal.toComplex (v,constant 0 (dim v))
lib/Data/Packed/Internal/Common.hs view
@@ -1,4 +1,4 @@-{-# OPTIONS_GHC -fglasgow-exts #-}+{-# LANGUAGE CPP #-} ----------------------------------------------------------------------------- -- | -- Module      :  Data.Packed.Internal.Common@@ -82,12 +82,14 @@   -- | clear the fpu-foreign import ccall "auxi.h asm_finit" finit :: IO ()+foreign import ccall "asm_finit" finit :: IO ()  -- | check the error code check :: String -> IO CInt -> IO () check msg f = do+#if FINIT     finit+#endif     err <- f     when (err/=0) $ if err > 1024                       then (error (msg++": "++errorCode err)) -- our errors@@ -101,9 +103,6 @@ foreign import ccall "auxi.h gsl_strerror" gsl_strerror :: CInt -> IO (Ptr CChar)  ------------------------------------------------------ ugly, but my haddock version doesn't understand--- yet infix type constructors---------------------------------------------------- ---------- signatures of the C functions --------- -------------------------------------------------- type PD = Ptr Double                            --@@ -141,8 +140,3 @@ type TMCVM = CInt -> CInt -> PD -> TCVM         -- type TMMCVM = CInt -> CInt -> PD -> TMCVM       -- ----------------------------------------------------type TauxMul a = CInt -> CInt -> CInt -> Ptr a-               -> CInt -> CInt -> CInt -> Ptr a-               -> CInt -> CInt -> Ptr a-               -> IO CInt
lib/Data/Packed/Internal/Matrix.hs view
@@ -1,5 +1,5 @@ {-# OPTIONS_GHC -fglasgow-exts #-}-{-# LANGUAGE CPP #-}+{-# LANGUAGE CPP, BangPatterns #-} ----------------------------------------------------------------------------- -- | -- Module      :  Data.Packed.Internal.Matrix@@ -71,6 +71,12 @@ -- MC: preferred by C, fdat may require a transposition -- MF: preferred by LAPACK, cdat may require a transposition +xdat MC {cdat = d } = d+xdat MF {fdat = d } = d++orderOf MF{} = ColumnMajor+orderOf MC{} = RowMajor+ -- | Matrix transpose. trans :: Matrix t -> Matrix t trans MC {rows = r, cols = c, cdat = d } = MF {rows = c, cols = r, fdat = d }@@ -84,16 +90,10 @@  mat = withMatrix -withMatrix MC {rows = r, cols = c, cdat = d } f =-    withForeignPtr (fptr d) $ \p -> do-        let m g = do-            g (fi r) (fi c) p-        f m--withMatrix MF {rows = r, cols = c, fdat = d } f =-    withForeignPtr (fptr d) $ \p -> do+withMatrix a f =+    withForeignPtr (fptr (xdat a)) $ \p -> do         let m g = do-            g (fi r) (fi c) p+            g (fi (rows a)) (fi (cols a)) p         f m  {- | Creates a vector by concatenation of rows@@ -156,6 +156,11 @@     | otherwise = v `at` (j*r+i) {-# INLINE (@@>) #-} +-- | Unsafe matrix access without range checking+atM' MC {cols = c, cdat = v} i j = v `at'` (i*c+j)+atM' MF {rows = r, fdat = v} i j = v `at'` (j*r+i)+{-# INLINE atM' #-}+ ------------------------------------------------------------------  matrixFromVector RowMajor c v = MC { rows = r, cols = c, cdat = v }@@ -204,45 +209,55 @@ compat :: Matrix a -> Matrix b -> Bool compat m1 m2 = rows m1 == rows m2 && cols m1 == cols m2 -----------------------------------------------------------------+------------------------------------------------------------------ --- | Optimized matrix computations are provided for elements in the Element class.+-- | Auxiliary class. class (Storable a, Floating a) => Element a where-    constantD :: a -> Int -> Vector a-    transdata :: Int -> Vector a -> Int -> Vector a     subMatrixD :: (Int,Int) -- ^ (r0,c0) starting position                 -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix                -> Matrix a -> Matrix a-    diagD :: Vector a -> Matrix a+    transdata :: Int -> Vector a -> Int -> Vector a+    constantD  :: a -> Int -> Vector a  instance Element Double where-    constantD  = constantR-    transdata = transdataR-    subMatrixD = subMatrixR-    diagD      = diagR+    subMatrixD = subMatrix'+    transdata  = transdataAux ctransR     -- transdata'+    constantD  = constantAux cconstantR   -- constant'  instance Element (Complex Double) where-    constantD  = constantC-    transdata  = transdataC-    subMatrixD = subMatrixC-    diagD      = diagC+    subMatrixD = subMatrix'+    transdata  = transdataAux ctransC     -- transdata'+    constantD  = constantAux cconstantC   -- constant' -------------------------------------------------------------------+------------------------------------------------------------------- -(>|<) :: (Element a) => Int -> Int -> [a] -> Matrix a-r >|< c = f where-    f l | dim v == r*c = matrixFromVector ColumnMajor c v-        | otherwise    = error $ "inconsistent list size = "-                                 ++show (dim v) ++" in ("++show r++"><"++show c++")"-        where v = fromList l+transdata' :: Storable a => Int -> Vector a -> Int -> Vector a+transdata' c1 v c2 =+    if noneed+        then v+        else unsafePerformIO $ do+                w <- createVector (r2*c2)+                withForeignPtr (fptr v) $ \p ->+                    withForeignPtr (fptr w) $ \q -> do+                        let go (-1) _ = return ()+                            go !i (-1) = go (i-1) (c1-1)+                            go !i !j = do x <- peekElemOff p (i*c1+j)+                                          pokeElemOff      q (j*c2+i) x+                                          go i (j-1)+                        go (r1-1) (c1-1)+                return w+  where r1 = dim v `div` c1+        r2 = dim v `div` c2+        noneed = r1 == 1 || c1 == 1 --------------------------------------------------------------------+-- {-# SPECIALIZE transdata' :: Int -> Vector Double -> Int ->  Vector Double #-}+-- {-# SPECIALIZE transdata' :: Int -> Vector (Complex Double) -> Int -> Vector (Complex Double) #-} -transdataR :: Int -> Vector Double -> Int -> Vector Double-transdataR = transdataAux ctransR+-- I don't know how to specialize...+-- The above pragmas only seem to work on top level defs+-- Fortunately everything seems to work using the above class -transdataC :: Int -> Vector (Complex Double) -> Int -> Vector (Complex Double)-transdataC = transdataAux ctransC+-- C versions, still a little faster:  transdataAux fun c1 d c2 =     if noneed@@ -252,64 +267,24 @@             withForeignPtr (fptr d) $ \pd ->                 withForeignPtr (fptr v) $ \pv ->                     fun (fi r1) (fi c1) pd (fi r2) (fi c2) pv // check "transdataAux"-            -- putStrLn $ "---> transdataAux" ++ show (toList d) ++ show (toList v)             return v   where r1 = dim d `div` c1         r2 = dim d `div` c2         noneed = r1 == 1 || c1 == 1 -foreign import ccall "auxi.h transR" ctransR :: TMM-foreign import ccall "auxi.h transC" ctransC :: TCMCM-+foreign import ccall "transR" ctransR :: TMM+foreign import ccall "transC" ctransC :: TCMCM ---------------------------------------------------------------------- --- | extraction of a submatrix from a real matrix-subMatrixR :: (Int,Int) -> (Int,Int) -> Matrix Double -> Matrix Double-subMatrixR (r0,c0) (rt,ct) x' = unsafePerformIO $ do-    r <- createMatrix RowMajor rt ct-    let x = cmat x'-    app2 (c_submatrixR (fi r0) (fi $ r0+rt-1) (fi c0) (fi $ c0+ct-1)) mat x mat r "subMatrixR"-    return r-foreign import ccall "auxi.h submatrixR" c_submatrixR :: CInt -> CInt -> CInt -> CInt -> TMM---- | extraction of a submatrix from a complex matrix-subMatrixC :: (Int,Int) -> (Int,Int) -> Matrix (Complex Double) -> Matrix (Complex Double)-subMatrixC (r0,c0) (rt,ct) x =-    reshape ct . asComplex . flatten .-    subMatrixR (r0,2*c0) (rt,2*ct) .-    reshape (2*cols x) . asReal . flatten $ x---- | Extracts a submatrix from a matrix.-subMatrix :: Element a-          => (Int,Int) -- ^ (r0,c0) starting position -          -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix-          -> Matrix a -- ^ input matrix-          -> Matrix a -- ^ result-subMatrix = subMatrixD--------------------------------------------------------------------------diagAux fun msg (v@V {dim = n}) = unsafePerformIO $ do-    m <- createMatrix RowMajor n n-    app2 fun vec v mat m msg-    return m---- | diagonal matrix from a real vector-diagR :: Vector Double -> Matrix Double-diagR = diagAux c_diagR "diagR"-foreign import ccall "auxi.h diagR" c_diagR :: TVM---- | diagonal matrix from a real vector-diagC :: Vector (Complex Double) -> Matrix (Complex Double)-diagC = diagAux c_diagC "diagC"-foreign import ccall "auxi.h diagC" c_diagC :: TCVCM---- | creates a square matrix with the given diagonal-diag :: Element a => Vector a -> Matrix a-diag = diagD+constant' v n = unsafePerformIO $ do+    w <- createVector n+    withForeignPtr (fptr w) $ \p -> do+        let go (-1) = return ()+            go !k = pokeElemOff p k v >> go (k-1)+        go (n-1)+    return w -------------------------------------------------------------------------+-- C versions  constantAux fun x n = unsafePerformIO $ do     v <- createVector n@@ -320,30 +295,45 @@  constantR :: Double -> Int -> Vector Double constantR = constantAux cconstantR-foreign import ccall "auxi.h constantR" cconstantR :: Ptr Double -> TV+foreign import ccall "constantR" cconstantR :: Ptr Double -> TV  constantC :: Complex Double -> Int -> Vector (Complex Double) constantC = constantAux cconstantC-foreign import ccall "auxi.h constantC" cconstantC :: Ptr (Complex Double) -> TCV+foreign import ccall "constantC" cconstantC :: Ptr (Complex Double) -> TCV+---------------------------------------------------------------------- -{- | creates a vector with a given number of equal components:+-- | Extracts a submatrix from a matrix.+subMatrix :: Element a+          => (Int,Int) -- ^ (r0,c0) starting position +          -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix+          -> Matrix a -- ^ input matrix+          -> Matrix a -- ^ result+subMatrix (r0,c0) (rt,ct) m+    | 0 <= r0 && 0 < rt && r0+rt <= (rows m) &&+      0 <= c0 && 0 < ct && c0+ct <= (cols m) = subMatrixD (r0,c0) (rt,ct) m+    | otherwise = error $ "wrong subMatrix "+++                          show ((r0,c0),(rt,ct))++" of "++show(rows m)++"x"++ show (cols m) -@> constant 2 7-7 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@--}-constant :: Element a => a -> Int -> Vector a-constant = constantD+subMatrix'' (r0,c0) (rt,ct) c v = unsafePerformIO $ do+    w <- createVector (rt*ct)+    withForeignPtr (fptr v) $ \p ->+        withForeignPtr (fptr w) $ \q -> do+            let go (-1) _ = return ()+                go !i (-1) = go (i-1) (ct-1)+                go !i !j = do x <- peekElemOff p ((i+r0)*c+j+c0)+                              pokeElemOff      q (i*ct+j) x+                              go i (j-1)+            go (rt-1) (ct-1)+    return w +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 conj :: Vector (Complex Double) -> Vector (Complex Double)-conj v = unsafePerformIO $ do-    r <- createVector (dim v)-    app2 cconjugate vec v vec r "cconjugate"-    return r-foreign import ccall "auxi.h conjugate" cconjugate :: TCVCV-+conj = mapVector conjugate  -- | creates a complex vector from vectors with real and imaginary parts toComplex :: (Vector Double, Vector Double) ->  Vector (Complex Double)@@ -354,10 +344,6 @@ fromComplex z = (r,i) where     [r,i] = toColumns $ reshape 2 $ asReal z --- | converts a real vector into a complex representation (with zero imaginary parts)-comp :: Vector Double -> Vector (Complex Double)-comp v = toComplex (v,constant 0 (dim v))- -- | loads a matrix efficiently from formatted ASCII text file (the number of rows and columns must be known in advance). fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double) fromFile filename (r,c) = do@@ -366,4 +352,4 @@     app1 (c_gslReadMatrix charname) mat res "gslReadMatrix"     --free charname  -- TO DO: free the auxiliary CString     return res-foreign import ccall "auxi.h matrix_fscanf" c_gslReadMatrix:: Ptr CChar -> TM+foreign import ccall "matrix_fscanf" c_gslReadMatrix:: Ptr CChar -> TM
lib/Data/Packed/Internal/Vector.hs view
@@ -1,4 +1,4 @@-{-# LANGUAGE MagicHash, CPP, UnboxedTuples #-}+{-# LANGUAGE MagicHash, CPP, UnboxedTuples, BangPatterns #-} ----------------------------------------------------------------------------- -- | -- Module      :  Data.Packed.Internal.Vector@@ -99,6 +99,7 @@ -- | access to Vector elements without range checking at' :: Storable a => Vector a -> Int -> a at' v n = safeRead v $ flip peekElemOff n+{-# INLINE at' #-}  -- -- turn off bounds checking with -funsafe at configure time.@@ -179,13 +180,9 @@  ---------------------------------------------------------------- --- | map on Vectors-liftVector :: (Storable a, Storable b) => (a-> b) -> Vector a -> Vector b-liftVector  f = fromList . map f . toList+liftVector f x = mapVector f x --- | zipWith for Vectors-liftVector2 :: (Storable a, Storable b, Storable c) => (a-> b -> c) -> Vector a -> Vector b -> Vector c-liftVector2 f u v = fromList $ zipWith f (toList u) (toList v)+liftVector2 f u v = zipVector f u v  ----------------------------------------------------------------- @@ -195,3 +192,54 @@         let f _ s _ d =  copyArray d s n >> return 0         app2 f vec v vec r "cloneVector"         return r++------------------------------------------------------------------++-- | map on Vectors+mapVector :: (Storable a, Storable b) => (a-> b) -> Vector a -> Vector b+mapVector f v = unsafePerformIO $ do+    w <- createVector (dim v)+    withForeignPtr (fptr v) $ \p ->+        withForeignPtr (fptr w) $ \q -> do+            let go (-1) = return ()+                go !k = do x <- peekElemOff p k+                           pokeElemOff      q k (f x)+                           go (k-1)+            go (dim v -1)+    return w+{-# 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+    let n = min (dim u) (dim v)+    w <- createVector n+    withForeignPtr (fptr u) $ \pu ->+        withForeignPtr (fptr v) $ \pv ->+            withForeignPtr (fptr w) $ \pw -> do+                let go (-1) = return ()+                    go !k = do x <- peekElemOff pu k+                               y <- peekElemOff pv k+                               pokeElemOff      pw k (f x y)+                               go (k-1)+                go (n -1)+    return w+{-# INLINE zipVector #-}++foldVector f x v = unsafePerformIO $+    withForeignPtr (fptr (v::Vector Double)) $ \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 #-}++foldLoop f s0 d = go (d - 1) s0+     where+       go 0 s = f (0::Int) s+       go !j !s = go (j - 1) (f j s)++foldVectorG f s0 v = foldLoop g s0 (dim v)+    where g !k !s = f k (at' v) s+          {-# INLINE g #-} -- Thanks to Ryan Ingram (http://permalink.gmane.org/gmane.comp.lang.haskell.cafe/46479)+{-# INLINE foldVectorG #-}
− lib/Data/Packed/Internal/auxi.c 
@@ -1,157 +0,0 @@-#include "auxi.h"-#include <gsl/gsl_blas.h>-#include <gsl/gsl_linalg.h>-#include <gsl/gsl_matrix.h>-#include <gsl/gsl_math.h>-#include <gsl/gsl_errno.h>-#include <gsl/gsl_complex.h>-#include <gsl/gsl_complex_math.h>-#include <gsl/gsl_cblas.h>-#include <string.h>-#include <stdio.h>--#define MACRO(B) do {B} while (0)-#define ERROR(CODE) MACRO(return CODE;)-#define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);})-#define OK return 0;--#define MIN(A,B) ((A)<(B)?(A):(B))-#define MAX(A,B) ((A)>(B)?(A):(B))--#ifdef DBG-#define DEBUGMSG(M) printf("*** calling aux C function: %s\n",M);-#else-#define DEBUGMSG(M)-#endif--#define CHECK(RES,CODE) MACRO(if(RES) return CODE;)--#ifdef DBG-#define DEBUGMAT(MSG,X) printf(MSG" = \n"); gsl_matrix_fprintf(stdout,X,"%f"); printf("\n");-#else-#define DEBUGMAT(MSG,X)-#endif--#ifdef DBG-#define DEBUGVEC(MSG,X) printf(MSG" = \n"); gsl_vector_fprintf(stdout,X,"%f"); printf("\n");-#else-#define DEBUGVEC(MSG,X)-#endif--#define DVVIEW(A) gsl_vector_view A = gsl_vector_view_array(A##p,A##n)-#define DMVIEW(A) gsl_matrix_view A = gsl_matrix_view_array(A##p,A##r,A##c)-#define CVVIEW(A) gsl_vector_complex_view A = gsl_vector_complex_view_array((double*)A##p,A##n)-#define CMVIEW(A) gsl_matrix_complex_view A = gsl_matrix_complex_view_array((double*)A##p,A##r,A##c)-#define KDVVIEW(A) gsl_vector_const_view A = gsl_vector_const_view_array(A##p,A##n)-#define KDMVIEW(A) gsl_matrix_const_view A = gsl_matrix_const_view_array(A##p,A##r,A##c)-#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 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 BAD_SIZE 2000-#define BAD_CODE 2001-#define MEM      2002-#define BAD_FILE 2003--int transR(KRMAT(x),RMAT(t)) {-    REQUIRES(xr==tc && xc==tr,BAD_SIZE);-    DEBUGMSG("transR");-    KDMVIEW(x);-    DMVIEW(t);-    int res = gsl_matrix_transpose_memcpy(M(t),M(x));-    CHECK(res,res);-    OK-}--int transC(KCMAT(x),CMAT(t)) {-    REQUIRES(xr==tc && xc==tr,BAD_SIZE);-    DEBUGMSG("transC");-    KCMVIEW(x);-    CMVIEW(t);-    int res = gsl_matrix_complex_transpose_memcpy(M(t),M(x));-    CHECK(res,res);-    OK-}---int submatrixR(int r1, int r2, int c1, int c2, KRMAT(x),RMAT(r)) {-    REQUIRES(0<=r1 && r1<=r2 && r2<xr && 0<=c1 && c1<=c2 && c2<xc &&-            rr==r2-r1+1 && rc==c2-c1+1,BAD_SIZE);-    DEBUGMSG("submatrixR");-    KDMVIEW(x);-    DMVIEW(r);-    gsl_matrix_const_view S = gsl_matrix_const_submatrix(M(x),r1,c1,rr,rc);-    int res = gsl_matrix_memcpy(M(r),M(S));-    CHECK(res,res);-    OK-}---int constantR(double * pval, RVEC(r)) {-    DEBUGMSG("constantR")-    int k;-    double val = *pval;-    for(k=0;k<rn;k++) {-        rp[k]=val;-    }-    OK-}--int constantC(gsl_complex* pval, CVEC(r)) {-    DEBUGMSG("constantC")-    int k;-    gsl_complex val = *pval;-    for(k=0;k<rn;k++) {-        rp[k]=val;-    }-    OK-}---int diagR(KRVEC(d),RMAT(r)) {-    REQUIRES(dn==rr && rr==rc,BAD_SIZE);-    DEBUGMSG("diagR");-    int i,j;-    for (i=0;i<rr;i++) {-        for(j=0;j<rc;j++) {-            rp[i*rc+j] = i==j?dp[i]:0.;-        }-    }-    OK-}--int diagC(KCVEC(d),CMAT(r)) {-    REQUIRES(dn==rr && rr==rc,BAD_SIZE);-    DEBUGMSG("diagC");-    int i,j;-    gsl_complex zero;-    GSL_SET_COMPLEX(&zero,0.,0.);-    for (i=0;i<rr;i++) {-        for(j=0;j<rc;j++) {-            rp[i*rc+j] = i==j?dp[i]:zero;-        }-    }-    OK-}--int conjugate(KCVEC(x),CVEC(t)) {-    REQUIRES(xn==tn,BAD_SIZE);-    DEBUGMSG("conjugate");-    int k;-    for (k=0; k<xn; k++) {-        tp[k].dat[0] =   xp[k].dat[0];-        tp[k].dat[1] = - xp[k].dat[1];-    }-    OK-}--//----------------------------------------void asm_finit() {-    asm("finit");-}-//---------------------------------------
− lib/Data/Packed/Internal/auxi.h 
@@ -1,30 +0,0 @@-#include <gsl/gsl_complex.h>--#define RVEC(A) int A##n, double*A##p-#define RMAT(A) int A##r, int A##c, double* A##p-#define KRVEC(A) int A##n, const double*A##p-#define KRMAT(A) int A##r, int A##c, const double* A##p--#define CVEC(A) int A##n, gsl_complex*A##p-#define CMAT(A) int A##r, int A##c, gsl_complex* A##p-#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--int transR(KRMAT(x),RMAT(t));-int transC(KCMAT(x),CMAT(t));--int constantR(double *val     , RVEC(r));-int constantC(gsl_complex *val, CVEC(r));--int submatrixR(int r1, int r2, int c1, int c2, KRMAT(x),RMAT(r));--int diagR(KRVEC(d),RMAT(r));-int diagC(KCVEC(d),CMAT(r));--const char * gsl_strerror (const int gsl_errno);--int matrix_fscanf(char*filename, RMAT(a));--int conjugate(KCVEC(x),CVEC(t));--void asm_finit();
lib/Data/Packed/Matrix.hs view
@@ -33,6 +33,7 @@ ) where  import Data.Packed.Internal+import qualified Data.Packed.ST as ST import Data.Packed.Vector import Data.List(transpose,intersperse) import Data.Array@@ -73,21 +74,30 @@  ------------------------------------------------------------ +-- | 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 -@> diagRect (constant 5 3) 3 4+@> 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 :: (Element t, Num t) => Vector t -> Int -> Int -> Matrix t-diagRect s r c-    | dim s < min r c = error "diagRect"-    | r == c    = diag s-    | r < c     = trans $ diagRect s c r-    | otherwise = joinVert  [diag s , zeros (r-c,c)]-    where zeros (r',c') = reshape c' $ constantD 0 (r'*c')+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+        mapM_ (\k -> ST.writeMatrix m k k (v@>k)) [0..d-1]+        return m  -- | extracts the diagonal from a rectangular matrix takeDiag :: (Element t) => Matrix t -> Vector t
lib/Data/Packed/ST.hs view
@@ -1,4 +1,4 @@-{-# OPTIONS -XTypeOperators -XRank2Types  -XFlexibleContexts #-}+{-# OPTIONS -XTypeOperators -XRank2Types  -XFlexibleContexts -XBangPatterns #-}  ----------------------------------------------------------------------------- -- |@@ -23,8 +23,10 @@     STMatrix, newMatrix, thawMatrix, freezeMatrix, runSTMatrix,     readMatrix, writeMatrix, modifyMatrix, liftSTMatrix,     -- * Unsafe functions+    newUndefinedVector,     unsafeReadVector, unsafeWriteVector,     unsafeThawVector, unsafeFreezeVector,+    newUndefinedMatrix,     unsafeReadMatrix, unsafeWriteMatrix,     unsafeThawMatrix, unsafeFreezeMatrix ) where@@ -87,9 +89,17 @@ writeVector :: Storable t => STVector s t -> Int -> t -> ST s () writeVector = safeIndexV unsafeWriteVector -{-# NOINLINE newVector #-}+{-# NOINLINE newUndefinedVector #-}+newUndefinedVector :: Element 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 v = unsafeThawVector . constant v+newVector x n = do+    v <- newUndefinedVector n+    let go (-1) = return v+        go !k = unsafeWriteVector v k x >> go (k-1 :: Int)+    go (n-1)  ------------------------------------------------------------------------- @@ -153,6 +163,10 @@ writeMatrix :: Storable t => STMatrix s t -> Int -> Int -> t -> ST s () writeMatrix = safeIndexM unsafeWriteMatrix +{-# NOINLINE newUndefinedMatrix #-}+newUndefinedMatrix :: Element t => MatrixOrder -> Int -> Int -> ST s (STMatrix s t)+newUndefinedMatrix order r c = unsafeIOToST $ fmap STMatrix $ createMatrix order r c+ {-# NOINLINE newMatrix #-} newMatrix :: Element t => t -> Int -> Int -> ST s (STMatrix s t)-newMatrix v r c = unsafeThawMatrix . reshape c . constant v $ r*c+newMatrix v r c = unsafeThawMatrix $ reshape c $ runSTVector $ newVector v (r*c)
lib/Data/Packed/Vector.hs view
@@ -19,11 +19,13 @@     subVector, join,     constant, linspace,     vectorMax, vectorMin, vectorMaxIndex, vectorMinIndex,-    liftVector, liftVector2+    liftVector, liftVector2,+    foldLoop, foldVector, foldVectorG, mapVector, zipVector ) where  import Data.Packed.Internal import Numeric.GSL.Vector+-- import Data.Packed.ST  {- | Creates a real vector containing a range of values: @@ -47,3 +49,12 @@  vectorMinIndex :: Vector Double -> Int vectorMinIndex = round . toScalarR MinIdx++{- | creates a vector with a given number of equal components:++@> constant 2 7+7 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@+-}+constant :: Element a => a -> Int -> Vector a+-- constant x n = runSTVector (newVector x n)+constant = constantD -- about 2x faster
lib/Graphics/Plot.hs view
@@ -1,7 +1,7 @@ ----------------------------------------------------------------------------- -- | -- Module      :  Graphics.Plot--- Copyright   :  (c) Alberto Ruiz 2005+-- Copyright   :  (c) Alberto Ruiz 2005-8 -- License     :  GPL-style --  -- Maintainer  :  Alberto Ruiz (aruiz at um dot es)@@ -28,15 +28,11 @@  ) where -import Data.Packed.Vector-import Data.Packed.Matrix+import Data.Packed import Numeric.LinearAlgebra(outer)-import Numeric.GSL.Vector(FunCodeS(Max,Min),toScalarR) import Data.List(intersperse) import System-import Foreign hiding (rotate) - size = dim  -- | Loads a real matrix from a formatted ASCII text file @@ -154,8 +150,8 @@     r = rows m     header = "P2 "++show c++" "++show r++" "++show (round maxgray :: Int)++"\n"     maxgray = 255.0-    maxval = toScalarR Max $ flatten $ m-    minval = toScalarR Min $ flatten $ m+    maxval = vectorMax $ flatten $ m+    minval = vectorMin $ flatten $ m     scale = if (maxval == minval)          then 0.0         else maxgray / (maxval - minval)
lib/Numeric/GSL/Minimization.hs view
@@ -2,14 +2,14 @@ ----------------------------------------------------------------------------- {- | Module      :  Numeric.GSL.Minimization-Copyright   :  (c) Alberto Ruiz 2006+Copyright   :  (c) Alberto Ruiz 2006-9 License     :  GPL-style  Maintainer  :  Alberto Ruiz (aruiz at um dot es) Stability   :  provisional Portability :  uses ffi -Minimization of a multidimensional function Minimization of a multidimensional function using some of the algorithms described in:+Minimization of a multidimensional function using some of the algorithms described in:  <http://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Minimization.html> @@ -17,6 +17,7 @@ ----------------------------------------------------------------------------- module Numeric.GSL.Minimization (     minimizeConjugateGradient,+    minimizeVectorBFGS2,     minimizeNMSimplex ) where @@ -132,7 +133,7 @@  @'Graphics.Plot.mplot' $ drop 2 ('toColumns' p)@ --}     +-} minimizeConjugateGradient ::        Double        -- ^ initial step size     -> Double        -- ^ minimization parameter   @@ -142,7 +143,27 @@     -> ([Double] -> [Double])      -- ^ gradient     -> [Double]             -- ^ starting point     -> ([Double], Matrix Double)        -- ^ solution vector, and the optimization trajectory followed by the algorithm-minimizeConjugateGradient istep minimpar tol maxit f df xi = unsafePerformIO $ do+minimizeConjugateGradient = minimizeWithDeriv 0++{- | Taken from the GSL manual:++The vector Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm. This is a quasi-Newton method which builds up an approximation to the second derivatives of the function f using the difference between successive gradient vectors. By combining the first and second derivatives the algorithm is able to take Newton-type steps towards the function minimum, assuming quadratic behavior in that region.++The bfgs2 version of this minimizer is the most efficient version available, and is a faithful implementation of the line minimization scheme described in Fletcher's Practical Methods of Optimization, Algorithms 2.6.2 and 2.6.4. It supercedes the original bfgs routine and requires substantially fewer function and gradient evaluations. The user-supplied tolerance tol corresponds to the parameter \sigma used by Fletcher. A value of 0.1 is recommended for typical use (larger values correspond to less accurate line searches).+-}+minimizeVectorBFGS2 ::+       Double        -- ^ initial step size+    -> Double        -- ^ minimization parameter tol+    -> Double        -- ^ desired precision of the solution (gradient test)+    -> Int           -- ^ maximum number of iterations allowed+    -> ([Double] -> Double) -- ^ function to minimize+    -> ([Double] -> [Double])      -- ^ gradient+    -> [Double]             -- ^ starting point+    -> ([Double], Matrix Double)        -- ^ solution vector, and the optimization trajectory followed by the algorithm+minimizeVectorBFGS2 = minimizeWithDeriv 1+++minimizeWithDeriv method istep minimpar tol maxit f df xi = unsafePerformIO $ do     let xiv = fromList xi         n = dim xiv         f' = f . toList@@ -151,7 +172,7 @@     dfp <- mkVecVecfun (aux_vTov df')     rawpath <- withVector xiv $ \xiv' ->                     createMIO maxit (n+2)-                         (c_minimizeConjugateGradient fp dfp istep minimpar tol (fi maxit) // xiv')+                         (c_minimizeWithDeriv method fp dfp istep minimpar tol (fi maxit) // xiv')                          "minimizeDerivV"     let it = round (rawpath @@> (maxit-1,0))         path = takeRows it rawpath@@ -160,9 +181,8 @@     freeHaskellFunPtr dfp     return (sol,path) - foreign import ccall "gsl-aux.h minimizeWithDeriv"-    c_minimizeConjugateGradient :: FunPtr (CInt -> Ptr Double -> Double)+    c_minimizeWithDeriv :: CInt -> FunPtr (CInt -> Ptr Double -> Double)                                 -> FunPtr (CInt -> Ptr Double -> Ptr Double -> IO ())                                 -> Double -> Double -> Double -> CInt                                 -> TVM
lib/Numeric/GSL/gsl-aux.c view
@@ -441,7 +441,7 @@ }  // conjugate gradient-int minimizeWithDeriv(double f(int, double*), void df(int, double*, double*), +int minimizeWithDeriv(int method, double f(int, double*), void df(int, double*, double*),                        double initstep, double minimpar, double tolgrad, int maxit,                        KRVEC(xi), RMAT(sol)) {     REQUIRES(solr == maxit && solc == 2+xin,BAD_SIZE);@@ -463,7 +463,11 @@     // Starting point     KDVVIEW(xi);     // conjugate gradient fr-    T = gsl_multimin_fdfminimizer_conjugate_fr;+    switch(method) {+        case 0 : {T = gsl_multimin_fdfminimizer_conjugate_fr; break; }+        case 1 : {T = gsl_multimin_fdfminimizer_vector_bfgs2; break; }+        default: ERROR(BAD_CODE);+    }     s = gsl_multimin_fdfminimizer_alloc (T, my_func.n);     gsl_multimin_fdfminimizer_set (s, &my_func, V(xi), initstep, minimpar);     do {
lib/Numeric/GSL/gsl-aux.h view
@@ -39,7 +39,7 @@ int minimize(double f(int, double*), double tolsize, int maxit,                   KRVEC(xi), KRVEC(sz), RMAT(sol)); -int minimizeWithDeriv(double f(int, double*), void df(int, double*, double*),+int minimizeWithDeriv(int method, double f(int, double*), void df(int, double*, double*),                       double initstep, double minimpar, double tolgrad, int maxit,                        KRVEC(xi), RMAT(sol));
lib/Numeric/LinearAlgebra.hs view
@@ -1,7 +1,7 @@ ----------------------------------------------------------------------------- {- | Module      :  Numeric.LinearAlgebra-Copyright   :  (c) Alberto Ruiz 2006-7+Copyright   :  (c) Alberto Ruiz 2006-9 License     :  GPL-style  Maintainer  :  Alberto Ruiz (aruiz at um dot es)@@ -9,6 +9,8 @@ Portability :  uses ffi  Basic matrix computations implemented by BLAS, LAPACK and GSL.++This is module reexports the most comon functions (including "Numeric.LinearAlgebra.Instances").  -} -----------------------------------------------------------------------------
lib/Numeric/LinearAlgebra/Algorithms.hs view
@@ -3,7 +3,7 @@ ----------------------------------------------------------------------------- {- | Module      :  Numeric.LinearAlgebra.Algorithms-Copyright   :  (c) Alberto Ruiz 2006-7+Copyright   :  (c) Alberto Ruiz 2006-9 License     :  GPL-style  Maintainer  :  Alberto Ruiz (aruiz at um dot es)@@ -20,27 +20,33 @@ -----------------------------------------------------------------------------  module Numeric.LinearAlgebra.Algorithms (--- * Linear Systems+-- * Supported types+    Field(),+-- * Products     multiply, dot,+    outer, kronecker,+-- * Linear Systems     linearSolve,+    luSolve,+    linearSolveSVD,     inv, pinv,-    pinvTol, det, rank, rcond,+    det, rank, rcond, -- * Matrix factorizations -- ** Singular value decomposition     svd,     full, economy, --thin, -- ** Eigensystems-    eig, eigSH,+    eig, eigSH, eigSH', -- ** QR     qr, -- ** Cholesky-    chol,+    chol, cholSH, -- ** Hessenberg     hess, -- ** Schur     schur, -- ** LU-    lu, luPacked, luSolve,+    lu, luPacked, -- * Matrix functions     expm,     sqrtm,@@ -53,15 +59,14 @@ -- * Misc     ctrans,     eps, i,-    outer, kronecker, -- * Util     haussholder,     unpackQR, unpackHess,-    Field(linearSolveSVD,eigSH',cholSH)+    pinvTol ) where  -import Data.Packed.Internal hiding (fromComplex, toComplex, comp, conj, (//))+import Data.Packed.Internal hiding (fromComplex, toComplex, conj, (//)) import Data.Packed import Numeric.GSL.Vector import Numeric.LinearAlgebra.LAPACK as LAPACK@@ -71,82 +76,120 @@ import Data.Array  - -- | 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-    -- | Singular value decomposition using lapack's dgesvd or zgesvd.-    svd         :: Matrix t -> (Matrix t, Vector Double, Matrix t)-    -- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'.-    luPacked    :: Matrix t -> (Matrix t, [Int])-    -- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization-    --   obtained by 'luPacked'.-    luSolve     :: (Matrix t, [Int]) -> Matrix t -> Matrix t-    -- | Solution of a general linear system (for several right-hand sides) using lapacks' dgesv or zgesv.-    -- It is similar to 'luSolve' . 'luPacked', but @linearSolve@ raises an error if called on a singular system.-    --  See also other versions of linearSolve in "Numeric.LinearAlgebra.LAPACK".-    linearSolve :: Matrix t -> Matrix t -> Matrix t-    linearSolveSVD :: Matrix t -> Matrix t -> Matrix t-    -- | Eigenvalues and eigenvectors of a general square matrix using lapack's dgeev or zgeev.-    ---    -- If @(s,v) = eig m@ then @m \<> v == v \<> diag s@-    eig         :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double))-    -- | Similar to eigSH without checking that the input matrix is hermitian or symmetric.-    eigSH'      :: Matrix t -> (Vector Double, Matrix t)-    -- | Similar to chol without checking that the input matrix is hermitian or symmetric.-    cholSH      :: Matrix t -> Matrix t-    -- | QR factorization using lapack's dgeqr2 or zgeqr2.-    ---    -- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular.-    qr          :: Matrix t -> (Matrix t, Matrix t)-    -- | Hessenberg factorization using lapack's dgehrd or zgehrd.-    ---    -- If @(p,h) = hess m@ then @m == p \<> h \<> ctrans p@, where p is unitary-    -- and h is in upper Hessenberg form.-    hess        :: Matrix t -> (Matrix t, Matrix t)-    -- | Schur factorization using lapack's dgees or zgees.-    ---    -- If @(u,s) = schur m@ then @m == u \<> s \<> ctrans u@, where u is unitary-    -- and s is a Shur matrix. A complex Schur matrix is upper triangular. A real Schur matrix is-    -- upper triangular in 2x2 blocks.-    ---    -- \"Anything that the Jordan decomposition can do, the Schur decomposition-    -- can do better!\" (Van Loan)-    schur       :: Matrix t -> (Matrix t, Matrix t)-    -- | Conjugate transpose.-    ctrans :: Matrix t -> Matrix t-    -- | Matrix product.-    multiply :: Matrix t -> Matrix t -> Matrix t+    svd'         :: Matrix t -> (Matrix t, Vector Double, Matrix t)+    luPacked'    :: Matrix t -> (Matrix t, [Int])+    luSolve'     :: (Matrix t, [Int]) -> Matrix t -> Matrix t+    linearSolve' :: Matrix t -> Matrix t -> Matrix t+    linearSolveSVD' :: Matrix t -> Matrix t -> Matrix t+    eig'         :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double))+    eigSH''      :: Matrix t -> (Vector Double, Matrix t)+    cholSH'      :: Matrix t -> Matrix t+    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  +-- | Singular value decomposition using lapack's dgesvd or zgesvd.+svd         :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)+svd = svd'++-- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'.+luPacked    :: Field t => Matrix t -> (Matrix t, [Int])+luPacked = luPacked'++-- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization+--   obtained by 'luPacked'.+luSolve     :: Field t => (Matrix t, [Int]) -> Matrix t -> Matrix t+luSolve = luSolve'++-- | Solution of a general linear system (for several right-hand sides) using lapacks' dgesv or zgesv.+-- It is similar to 'luSolve' . 'luPacked', but @linearSolve@ raises an error if called on a singular system.+--  See also other versions of linearSolve in "Numeric.LinearAlgebra.LAPACK".+linearSolve :: Field t => Matrix t -> Matrix t -> Matrix t+linearSolve = linearSolve'+linearSolveSVD :: Field t => Matrix t -> Matrix t -> Matrix t+linearSolveSVD = linearSolveSVD'++-- | Eigenvalues and eigenvectors of a general square matrix using lapack's dgeev or zgeev.+--+-- If @(s,v) = eig m@ then @m \<> v == v \<> diag s@+eig         :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double))+eig = eig'++-- | Similar to 'eigSH' without checking that the input matrix is hermitian or symmetric.+eigSH'      :: Field t => Matrix t -> (Vector Double, Matrix t)+eigSH' = eigSH''++-- | Similar to 'chol' without checking that the input matrix is hermitian or symmetric.+cholSH      :: Field t => Matrix t -> Matrix t+cholSH = cholSH'++-- | QR factorization using lapack's dgeqr2 or zgeqr2.+--+-- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular.+qr          :: Field t => Matrix t -> (Matrix t, Matrix t)+qr = qr'++-- | Hessenberg factorization using lapack's dgehrd or zgehrd.+--+-- If @(p,h) = hess m@ then @m == p \<> h \<> ctrans p@, where p is unitary+-- and h is in upper Hessenberg form.+hess        :: Field t => Matrix t -> (Matrix t, Matrix t)+hess = hess'++-- | Schur factorization using lapack's dgees or zgees.+--+-- If @(u,s) = schur m@ then @m == u \<> s \<> ctrans u@, where u is unitary+-- and s is a Shur matrix. A complex Schur matrix is upper triangular. A real Schur matrix is+-- upper triangular in 2x2 blocks.+--+-- \"Anything that the Jordan decomposition can do, the Schur decomposition+-- can do better!\" (Van Loan)+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 = multiply'++ instance Field Double where-    svd = svdR-    luPacked = luR-    luSolve (l_u,perm) = lusR l_u perm-    linearSolve = linearSolveR                 -- (luSolve . luPacked) ??-    linearSolveSVD = linearSolveSVDR Nothing-    ctrans = trans-    eig = eigR-    eigSH' = eigS-    cholSH = cholS-    qr = unpackQR . qrR-    hess = unpackHess hessR-    schur = schurR-    multiply = multiplyR+    svd' = svdR+    luPacked' = luR+    luSolve' (l_u,perm) = lusR l_u perm+    linearSolve' = linearSolveR                 -- (luSolve . luPacked) ??+    linearSolveSVD' = linearSolveSVDR Nothing+    ctrans' = trans+    eig' = eigR+    eigSH'' = eigS+    cholSH' = cholS+    qr' = unpackQR . qrR+    hess' = unpackHess hessR+    schur' = schurR+    multiply' = multiplyR  instance Field (Complex Double) where-    svd = svdC-    luPacked = luC-    luSolve (l_u,perm) = lusC l_u perm-    linearSolve = linearSolveC-    linearSolveSVD = linearSolveSVDC Nothing-    ctrans = conj . trans-    eig = eigC-    eigSH' = eigH-    cholSH = cholH-    qr = unpackQR . qrC-    hess = unpackHess hessC-    schur = schurC-    multiply = multiplyC+    svd' = svdC+    luPacked' = luC+    luSolve' (l_u,perm) = lusC l_u perm+    linearSolve' = linearSolveC+    linearSolveSVD' = linearSolveSVDC Nothing+    ctrans' = conj . trans+    eig' = eigC+    eigSH'' = eigH+    cholSH' = cholH+    qr' = unpackQR . qrC+    hess' = unpackHess hessC+    schur' = schurC+    multiply' = multiplyC   -- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix using lapack's dsyev or zheev.@@ -193,8 +236,8 @@ -- If @(u,d,v) = full svd m@ then @m == u \<> d \<> trans v@. full :: Element t       => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Matrix Double, Matrix t)-full svd' m = (u, d ,v) where-    (u,s,v) = svd' m+full svdFun m = (u, d ,v) where+    (u,s,v) = svdFun m     d = diagRect s r c     r = rows m     c = cols m@@ -204,8 +247,8 @@ -- If @(u,s,v) = economy svd m@ then @m == u \<> diag s \<> trans v@. economy :: Element t          => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Vector Double, Matrix t)-economy svd' m = (u', subVector 0 d s, v') where-    (u,s,v) = svd' m+economy svdFun m = (u', subVector 0 d s, v') where+    (u,s,v) = svdFun m     sl@(g:_) = toList s     s' = fromList . filter (>tol) $ sl     t = 1
lib/Numeric/LinearAlgebra/LAPACK.hs view
@@ -27,11 +27,8 @@     schurR, schurC ) where -import Data.Packed.Internal-import Data.Packed.Internal.Vector-import Data.Packed.Internal.Matrix-import Data.Packed.Vector-import Data.Packed.Matrix+import Data.Packed.Internal hiding (toComplex)+import Data.Packed import Numeric.GSL.Vector(vectorMapValR, FunCodeSV(Scale)) import Complex import Foreign
lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c view
@@ -34,6 +34,13 @@ #define NODEFPOS 2006 #define NOSPRTD  2007 +//---------------------------------------+void asm_finit() {+#ifdef i386+    asm("finit");+#endif+}+//---------------------------------------  //////////////////// real svd //////////////////////////////////// @@ -899,5 +906,53 @@            (doublecomplex*)ap,&lda,            (doublecomplex*)bp,&ldb,&beta,            (doublecomplex*)rp,&ldc);+    OK+}++//////////////////// transpose /////////////////////////++int transR(KDMAT(x),DMAT(t)) {+    REQUIRES(xr==tc && xc==tr,BAD_SIZE);+    DEBUGMSG("transR");+    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];+        }+    }+    OK+}++//////////////////// constant /////////////////////////++int constantR(double * pval, DVEC(r)) {+    DEBUGMSG("constantR")+    int k;+    double 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;+    }     OK }
lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h view
@@ -55,6 +55,12 @@ 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 transR(KDMAT(x),DMAT(t));+int transC(KCMAT(x),CMAT(t));++int constantR(double * pval, DVEC(r));+int constantC(doublecomplex* pval, CVEC(r));+ int svd_l_R(KDMAT(x),DMAT(u),DVEC(s),DMAT(v)); int svd_l_Rdd(KDMAT(x),DMAT(u),DVEC(s),DMAT(v)); int svd_l_C(KCMAT(a),CMAT(u),DVEC(s),CMAT(v));
lib/Numeric/LinearAlgebra/Tests.hs view
@@ -1,7 +1,9 @@+{-# LANGUAGE CPP #-}+{-# OPTIONS_GHC -fno-warn-unused-imports -fno-warn-incomplete-patterns #-} ----------------------------------------------------------------------------- {- | Module      :  Numeric.LinearAlgebra.Tests-Copyright   :  (c) Alberto Ruiz 2007+Copyright   :  (c) Alberto Ruiz 2007-9 License     :  GPL-style  Maintainer  :  Alberto Ruiz (aruiz at um dot es)@@ -22,18 +24,16 @@ import Numeric.LinearAlgebra import Numeric.LinearAlgebra.Tests.Instances import Numeric.LinearAlgebra.Tests.Properties-import Test.QuickCheck hiding (test)-import Test.HUnit hiding ((~:),test)+import Test.HUnit hiding ((~:),test,Testable) import System.Info import Data.List(foldl1') import Numeric.GSL hiding (sin,cos,exp,choose) import Prelude hiding ((^)) import qualified Prelude+#include "Tests/quickCheckCompat.h"  a ^ b = a Prelude.^ (b :: Int) -qCheck n = check defaultConfig {configSize = const n}- utest str b = TestCase $ assertBool str b  feye n = flipud (ident n) :: Matrix Double@@ -105,7 +105,17 @@  , 2.718281828459045  , 2.718281828459045 ] +--------------------------------------------------------------------- +minimizationTest = TestList [ utest "minimization conj grad" (minim1 f df [5,7] ~~ [1,2])+                            , utest "minimization bg2"       (minim2 f df [5,7] ~~ [1,2])+                            ]+    where f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30+          df [x,y] = [20*(x-1), 40*(y-2)]+          a ~~ b = fromList a |~| fromList b+          minim1 g dg ini = fst $ minimizeConjugateGradient 1E-2 1E-4 1E-3 30 g dg ini+          minim2 g dg ini = fst $ minimizeVectorBFGS2 1E-2 1E-2 1E-3 30 g dg ini+ ---------------------------------------------------------------------  rot :: Double -> Matrix Double@@ -182,9 +192,10 @@     test (\v -> ifft (fft v) |~| v)     putStrLn "------ vector operations"     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))+    test $ (\u -> cos u * tan u |~| sin (u::CM)) . liftMatrix makeUnitary     putStrLn "------ read . show"     test (\m -> (m::RM) == read (show m))     test (\m -> (m::CM) == read (show m))@@ -205,8 +216,13 @@         , exponentialTest         , utest "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5^3) < 1E-8)         , utest "polySolve" (polySolveProp [1,2,3,4])+        , minimizationTest         ]     return ()++makeUnitary v | realPart n > 1    = v */ n+              | otherwise = v+    where n = sqrt (conj v <.> v)  -- -- | Some additional tests on big matrices. They take a few minutes. -- runBigTests :: IO ()
lib/Numeric/LinearAlgebra/Tests/Instances.hs view
@@ -1,4 +1,5 @@-{-# LANGUAGE FlexibleContexts, UndecidableInstances #-}+{-# LANGUAGE FlexibleContexts, UndecidableInstances, CPP #-}+{-# OPTIONS_GHC -fno-warn-unused-imports #-} ----------------------------------------------------------------------------- {- | Module      :  Numeric.LinearAlgebra.Tests.Instances@@ -24,33 +25,72 @@     RM,CM, rM,cM ) where +++ import Numeric.LinearAlgebra-import Test.QuickCheck import Control.Monad(replicateM)+#include "quickCheckCompat.h" ++#if MIN_VERSION_QuickCheck(2,0,0)+shrinkListElementwise :: (Arbitrary a) => [a] -> [[a]]+shrinkListElementwise []     = []+shrinkListElementwise (x:xs) = [ y:xs | y  <- shrink x                 ]+                            ++ [ x:ys | ys <- shrinkListElementwise xs ]++shrinkPair :: (Arbitrary a, Arbitrary b) => (a,b) -> [(a,b)]+shrinkPair (a,b) = [ (a,x) | x <- shrink b ] ++ [ (x,b) | x <- shrink a ]+#endif++ instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where     arbitrary = do         re <- arbitrary         im <- arbitrary         return (re :+ im)-    coarbitrary = undefined +#if MIN_VERSION_QuickCheck(2,0,0)+    shrink (re :+ im) = +        [ u :+ v | (u,v) <- shrinkPair (re,im) ]+#else+    -- this has been moved to the 'Coarbitrary' class in QuickCheck 2+    coarbitrary = undefined +#endif+ chooseDim = sized $ \m -> choose (1,max 1 m)  instance (Field a, Arbitrary a) => Arbitrary (Vector a) where -   arbitrary = do m <- chooseDim-                  l <- vector m-                  return $ fromList l-   coarbitrary = undefined+    arbitrary = do m <- chooseDim+                   l <- vector m+                   return $ fromList l +#if MIN_VERSION_QuickCheck(2,0,0)+    -- shrink any one of the components+    shrink = map fromList . shrinkListElementwise . toList+                              +#else+    coarbitrary = undefined+#endif+ instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where      arbitrary = do         m <- chooseDim         n <- chooseDim         l <- vector (m*n)         return $ (m><n) l++#if MIN_VERSION_QuickCheck(2,0,0)+    -- shrink any one of the components+    shrink a = map ((rows a) >< (cols a))+               . shrinkListElementwise+               . concat . toLists +                     $ a+#else     coarbitrary = undefined+#endif + -- a square matrix newtype (Sq a) = Sq (Matrix a) deriving Show instance (Element a, Arbitrary a) => Arbitrary (Sq a) where@@ -58,8 +98,14 @@         n <- chooseDim         l <- vector (n*n)         return $ Sq $ (n><n) l++#if MIN_VERSION_QuickCheck(2,0,0)+    shrink (Sq a) = [ Sq b | b <- shrink a ]+#else     coarbitrary = undefined+#endif + -- a unitary matrix newtype (Rot a) = Rot (Matrix a) deriving Show instance (Field a, Arbitrary a) => Arbitrary (Rot a) where@@ -67,8 +113,13 @@         Sq m <- arbitrary         let (q,_) = qr m         return (Rot q)++#if MIN_VERSION_QuickCheck(2,0,0)+#else     coarbitrary = undefined+#endif + -- a complex hermitian or real symmetric matrix newtype (Her a) = Her (Matrix a) deriving Show instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Her a) where@@ -76,8 +127,13 @@         Sq m <- arbitrary         let m' = m/2         return $ Her (m' + ctrans m')++#if MIN_VERSION_QuickCheck(2,0,0)+#else     coarbitrary = undefined+#endif + -- 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@@ -90,8 +146,13 @@         sv <- replicateM n (choose (1,100))         let s = diagRect (fromList sv) r c         return $ WC (u <> real s <> trans v)++#if MIN_VERSION_QuickCheck(2,0,0)+#else     coarbitrary = undefined+#endif + -- 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@@ -102,8 +163,13 @@         sv <- replicateM n (choose (1,100))         let s = diag (fromList sv)         return $ SqWC (u <> real s <> trans v)++#if MIN_VERSION_QuickCheck(2,0,0)+#else     coarbitrary = undefined+#endif + -- 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@@ -115,8 +181,13 @@         let s = diag (fromList l)             p = v <> real s <> ctrans v         return $ PosDef (0.5 .* p + 0.5 .* ctrans p)++#if MIN_VERSION_QuickCheck(2,0,0)+#else     coarbitrary = undefined+#endif + -- a pair of matrices that can be multiplied newtype (Consistent a) = Consistent (Matrix a, Matrix a) deriving Show instance (Field a, Arbitrary a) => Arbitrary (Consistent a) where@@ -127,7 +198,13 @@         la <- vector (n*k)         lb <- vector (k*m)         return $ Consistent ((n><k) la, (k><m) lb)++#if MIN_VERSION_QuickCheck(2,0,0)+    shrink (Consistent (x,y)) = [ Consistent (u,v) | (u,v) <- shrinkPair (x,y) ]+#else     coarbitrary = undefined+#endif+   type RM = Matrix Double
lib/Numeric/LinearAlgebra/Tests/Properties.hs view
@@ -1,4 +1,5 @@-{-# OPTIONS #-}+{-# LANGUAGE CPP #-}+{-# OPTIONS_GHC -fno-warn-unused-imports #-} ----------------------------------------------------------------------------- {- | Module      :  Numeric.LinearAlgebra.Tests.Properties@@ -40,7 +41,7 @@ ) where  import Numeric.LinearAlgebra-import Test.QuickCheck+#include "quickCheckCompat.h" -- import Debug.Trace  -- debug x = trace (show x) x
+ lib/Numeric/LinearAlgebra/Tests/quickCheckCompat.h view
@@ -0,0 +1,33 @@+#ifndef MIN_VERSION_QuickCheck+#define MIN_VERSION_QuickCheck(A,B,C) 1+#endif++#if MIN_VERSION_QuickCheck(2,0,0)+import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector+                      ,sized,classify,Testable,Property++                      ,quickCheckWith,maxSize,stdArgs,shrink)++#else+import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector+                      ,sized,classify,Testable,Property++                      ,check,configSize,defaultConfig,trivial)+#endif++++#if MIN_VERSION_QuickCheck(2,0,0)+trivial :: Testable a => Bool -> a -> Property+trivial = (`classify` "trivial")+#else+#endif+++-- define qCheck, which used to be in Tests.hs+#if MIN_VERSION_QuickCheck(2,0,0)+qCheck n = quickCheckWith stdArgs {maxSize = n}+#else+qCheck n = check defaultConfig {configSize = const n}+#endif+