hmatrix 0.8.2.0 → 0.8.3.1

raw patch · 23 files changed

+546/−300 lines, 23 filesdep ~basePVP ok

version bump matches the API change (PVP)

Dependency ranges changed: base

API changes (from Hackage documentation)

Files

• CHANGES +10/−0
• README +4/−2
• examples/Real.hs +104/−0
• examples/benchmarks.hs +0/−151
• examples/latexmat.hs +0/−11
• examples/ode.hs +34/−0
• examples/pca1.hs +7/−8
• examples/pca2.hs +13/−15
• examples/pinv.hs +1/−1
• hmatrix.cabal +20/−20
• lib/Data/Packed/Internal/Common.hs +17/−6
• lib/Data/Packed/Internal/Matrix.hs +24/−18
• lib/Data/Packed/Internal/Vector.hs +11/−7
• lib/Data/Packed/Matrix.hs +71/−44
• lib/Numeric/GSL.hs +2/−0
• lib/Numeric/GSL/Internal.hs +6/−0
• lib/Numeric/GSL/ODE.hs +111/−0
• lib/Numeric/GSL/gsl-aux.c +88/−1
• lib/Numeric/LinearAlgebra/Algorithms.hs +1/−1
• lib/Numeric/LinearAlgebra/Instances.hs +7/−14
• lib/Numeric/LinearAlgebra/Interface.hs +4/−0
• lib/Numeric/LinearAlgebra/Tests.hs +10/−0
• lib/Numeric/LinearAlgebra/Tests/Properties.hs +1/−1

@@ -1,3 +1,12 @@+0.8.3.0+=======++- odeSolve++- Matrix arithmetic automatically replicates matrix with single row/column++- latexFormat, dispcf+ 0.8.2.0 ======= @@ -56,3 +65,4 @@ - added NFData instances for Matrix and Vector.  - liftVector, liftVector2 replaced by mapVector, zipVector.+

@@ -50,8 +50,8 @@  ACKNOWLEDGEMENTS ----------------------------------------------------- -I thank Don Stewart, Henning Thielemann, Bulat Ziganshin and all the people-in the Haskell mailing lists for their help.+I thank Don Stewart, Henning Thielemann, Bulat Ziganshin, Heinrich Apfelmus,+and all the people in the Haskell mailing lists for their help.  - Nico Mahlo discovered a bug in the eigendecomposition wrapper. @@ -108,3 +108,5 @@   to the values at that position.  - Jean-Francois Tremblay discovered an error in the tutorial.++- Heinrich Apfelmus fixed hmatrix.cabal for OS/X.

@@ -0,0 +1,104 @@++-- Alternative interface and utilities for creation of real arrays, useful to work in interactive mode.++module Real(+    module Numeric.LinearAlgebra,+    (<>), (*>), (<*), (<\>), (\>),+    vector,+    eye,+    zeros, ones,+    diagl,+    row,+    col,+    (#),(&), (//), blocks,+    rand+) where++import Numeric.LinearAlgebra hiding ((<>), (<|>), (<->), (<\>), (.*), (*/))+import System.Random(randomIO)++infixl 7 <>+-- | Matrix product ('multiply')+(<>) :: Field t => Matrix t -> Matrix t -> Matrix t+(<>) = multiply++infixl 7 *>+-- | matrix x vector+(*>) :: Field t => Matrix t -> Vector t -> Vector t+m *> v = flatten $ m <> (asColumn v)++infixl 7 <*+-- | vector x matrix+(<*) :: Field t => Vector t -> Matrix t -> Vector t+v <* m = flatten $ (asRow v) <> m+++-- | Least squares solution of a linear system for several right-hand sides, similar to the \\ operator of Matlab\/Octave. (\<\\\>) = 'linearSolveSVD'.+(<\>) :: (Field a) => Matrix a -> Matrix a -> Matrix a+infixl 7 <\>+(<\>) = linearSolveSVD++-- | Least squares solution of a linear system for a single right-hand side. See '(\<\\\>)'.+(\>) :: (Field a) => Matrix a -> Vector a -> Vector a+infixl 7 \>+m \> v = flatten (m <\> reshape 1 v)++-- | Pseudorandom matrix with uniform elements between 0 and 1.+rand :: Int -- ^ rows+     -> Int -- ^ columns+     -> IO (Matrix Double)+rand r c = do+    seed <- randomIO+    return (reshape c $ randomVector seed Uniform (r*c))++-- | Real identity matrix.+eye :: Int -> Matrix Double+eye = ident++-- | Create a real vector from a list.+vector :: [Double] -> Vector Double+vector = fromList++-- | Create a real diagonal matrix from a list.+diagl :: [Double] -> Matrix Double+diagl = diag . vector++-- | Create a matrix or zeros.+zeros :: Int -- ^ rows+      -> Int -- ^ columns+      -> Matrix Double+zeros r c = reshape c (constant 0 (r*c))++-- | Create a matrix or ones.+ones :: Int -- ^ rows+     -> Int -- ^ columns+     -> Matrix Double+ones r c = reshape c (constant 1 (r*c))++-- | Concatenation of real vectors.+infixl 9 #+(#) :: Vector Double -> Vector Double -> Vector Double+a # b = join [a,b]++-- | Horizontal concatenation of real matrices.+infixl 8 &+(&) :: Matrix Double -> Matrix Double -> Matrix Double+a & b = fromBlocks [[a,b]]++-- | Vertical concatenation of real matrices.+infixl 7 //+(//) :: Matrix Double -> Matrix Double -> Matrix Double+a // b = fromBlocks [[a],[b]]++-- | Real block matrix from a rectangular list of lists.+blocks :: [[Matrix Double]] -> Matrix Double+blocks = fromBlocks++-- | A real matrix with a single row, create from a list of elements.+row :: [Double] -> Matrix Double+row = asRow . vector++-- | A real matrix with a single column, created from a list of elements.+col :: [Double] -> Matrix Double+col = asColumn . vector+

@@ -1,151 +0,0 @@-{-# LANGUAGE BangPatterns #-}---- $ ghc --make -O2 benchmarks.hs---import Numeric.LinearAlgebra-import System.Time-import System.CPUTime-import Text.Printf-import Data.List(foldl1')---time act = do-    t0 <- getCPUTime-    act-    t1 <- getCPUTime-    printf "%.3f s CPU\n" $ (fromIntegral (t1 - t0) / (10^12 :: Double)) :: IO ()------------------------------------------------------------------------------------main = sequence_ [bench1,-                  bench2,-                  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--sumVH v = go (d - 1) 0-     where-       d = dim v-       go :: Int -> Double -> Double-       go 0 s = s + (v @> 0)-       go !j !s = go (j - 1) (s + (v @> j))--innerH u v = go (d - 1) 0-     where-       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-    putStrLn "-------------------------------------------------------"-    putStrLn "Multiplication of 1M different 3x3 matrices:"---    putStrLn "from [[]]"---    time $ print $ manymult (10^6) rot'---    putStrLn "from (3><3) []"-    time $ print $ manymult (10^6) rot-    print $ cos (10^6/2)---rot' :: Double -> Matrix Double-rot' a = matrix [[ c,0,s],-                 [ 0,1,0],-                 [-s,0,c]]-    where c = cos a-          s = sin a-          matrix = fromLists--rot :: Double -> Matrix Double-rot a = (3><3) [ c,0,s-               , 0,1,0-               ,-s,0,c ]-    where c = cos a-          s = sin a--manymult n r = foldl1' (<>) (map r angles)-    where angles = toList $ linspace n (0,1)-          -- angles = map (k*) [0..n']-          -- n' = fromIntegral n - 1-          -- k  = recip n'------------------------------------------------------------------------------------bench4 = do-    putStrLn "-------------------------------------------------------"-    putStrLn "1000x1000 inverse"-    let a = ident 1000 :: Matrix Double-    let b = 2*a-    print $ vectorMax $ flatten (a+b) -- evaluate it-    time $ print $ vectorMax $ flatten $ linearSolve a b------------------------------------------------------------------------------------op1 a b = a <> trans b-op2 a b = a  + trans b--timep = time . print . vectorMax . flatten--bench5 n d = do-    putStrLn "-------------------------------------------------------"-    putStrLn "transpose in add"-    let ms = replicate n ((ident d :: Matrix Double))-    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 "many constants"-    time $ print $ sum $ map ((@>0). flip constant sz) [1..n]

@@ -1,11 +0,0 @@-import Numeric.LinearAlgebra-import Text.Printf--disp w l fmt m = unlines $ map (++l) $ lines $ format w (printf fmt) m--latex fmt m = "\\begin{bmatrix}\n" ++ disp " & " " \\\\" fmt m ++ "\\end{bmatrix}"--main = do-    let m = (3><4) [1..12::Double]-    putStrLn $ disp " | " "" "%.2f" m-    putStrLn $ latex "%.3f" m

@@ -0,0 +1,34 @@+import Numeric.GSL.ODE+import Numeric.LinearAlgebra+import Graphics.Plot++vanderpol mu = do+    let xdot mu t [x,v] = [v, -x + mu * v * (1-x^2)]+        ts = linspace 1000 (0,50)+        sol = toColumns $ odeSolve (xdot mu) [1,0] ts+    mplot (ts : sol)+    mplot sol+++harmonic w d = do+    let xdot w d t [x,v] = [v, a*x + b*v] where a = -w^2; b = -2*d*w+        ts = linspace 100 (0,20)+        sol = odeSolve (xdot w d) [1,0] ts+    mplot (ts : toColumns sol)+++kepler v a = mplot (take 2 $ toColumns sol) where+    xdot t [x,y,vx,vy] = [vx,vy,x*k,y*k]+        where g=1+              k=(-g)*(x*x+y*y)**(-1.5)+    ts = linspace 100 (0,30)+    sol = odeSolve xdot [4, 0, v * cos (a*degree), v * sin (a*degree)] ts+    degree = pi/180+++main = do+    vanderpol 2+    harmonic 1 0+    harmonic 1 0.1+    kepler 0.3 60+    kepler 0.4 70

@@ -8,16 +8,15 @@ type Vec = Vector Double type Mat = Matrix Double -sumColumns m = constant 1 (rows m) <> m --- Vec with the mean value of the columns of a Mat-mean x = sumColumns x / fromIntegral (rows x)+-- Vector with the mean value of the columns of a matrix+mean a = constant (recip . fromIntegral . rows $ a) (rows a) <> a --- covariance Mat of a list of observations as rows of a Mat-cov x = (trans xc <> xc) / fromIntegral (rows x -1)-    where xc = center x-          center m = m - constant 1 (rows m) `outer` mean m+-- covariance matrix of a list of observations stored as rows+cov x = (trans xc <> xc) / fromIntegral (rows x - 1)+    where xc = x - asRow (mean x) + -- creates the compression and decompression functions from the desired number of components pca :: Int -> Mat -> (Vec -> Vec , Vec -> Vec) pca n dataSet = (encode,decode)@@ -38,7 +37,7 @@         system("wget -nv http://dis.um.es/~alberto/material/sp/mnist.txt.gz")         system("gunzip mnist.txt.gz")         return ()-    m <- fromFile "mnist.txt" (5000,785)+    m <- loadMatrix "mnist.txt" -- fromFile "mnist.txt" (5000,785)     let xs = takeColumns (cols m -1) m -- the last column is the digit type (class label)     let x = toRows xs !! 4  -- an arbitrary test Vec     let (pe,pd) = pca 10 xs

@@ -9,33 +9,31 @@ type Vec = Vector Double type Mat = Matrix Double -sumColumns m = constant 1 (rows m) <> m+-- Vector with the mean value of the columns of a matrix+mean a = constant (recip . fromIntegral . rows $ a) (rows a) <> a --- Vector with the mean value of the columns of a Mat-mean x = sumColumns x / fromIntegral (rows x)+-- covariance matrix of a list of observations stored as rows+cov x = (trans xc <> xc) / fromIntegral (rows x - 1)+    where xc = x - asRow (mean x) --- covariance matrix of a list of observations as rows of a matrix-cov x = (trans xc <> xc) / fromIntegral (rows x -1) -    where xc = center x-          center m = m - constant 1 (rows m) `outer` mean m  type Stat = (Vec, [Double], Mat)--- 1st and 2nd order statistics of a dataset (mean, eigenvalues and eigenvectors of cov)  +-- 1st and 2nd order statistics of a dataset (mean, eigenvalues and eigenvectors of cov) stat :: Mat -> Stat-stat x = (m, toList s, trans v) where   +stat x = (m, toList s, trans v) where     m = mean x-    (s,v) = eigSH' (cov x)   +    (s,v) = eigSH' (cov x) --- creates the compression and decompression functions from the desired reconstruction +-- creates the compression and decompression functions from the desired reconstruction -- quality and the statistics of a data set pca :: Double -> Stat -> (Vec -> Vec , Vec -> Vec) pca prec (m,s,v) = (encode,decode)       where         encode x = vp <> (x - m)     decode x = x <> vp + m-    vp = takeRows n v    +    vp = takeRows n v     n = 1 + (length $ fst $ span (< (prec'*sum s)) $ cumSum s)-    cumSum = tail . scanl (+) 0.0     +    cumSum = tail . scanl (+) 0.0     prec' = if prec <=0.0 || prec >= 1.0                 then error "the precision in pca must be 0<prec<1"                 else prec@@ -49,7 +47,7 @@     let (pe,pd) = pca prec st     let y = pe x     print $ dim y-    shdigit (pd y)    +    shdigit (pd y)  main = do     ok <- doesFileExist ("mnist.txt")@@ -58,7 +56,7 @@         system("wget -nv http://dis.um.es/~alberto/material/sp/mnist.txt.gz")         system("gunzip mnist.txt.gz")         return ()-    m <- fromFile "mnist.txt" (5000,785)+    m <- loadMatrix "mnist.txt"     let xs = takeColumns (cols m -1) m     let x = toRows xs !! 4  -- an arbitrary test vector     shdigit x

@@ -3,7 +3,7 @@ import Text.Printf(printf)  expand :: Int -> Vector Double -> Matrix Double-expand n x = fromColumns $ constant 1 (dim x): map (x^) [1 .. n]+expand n x = fromColumns $ map (x^) [0 .. n]  polynomialModel :: Vector Double -> Vector Double -> Int                 -> (Vector Double -> Vector Double)

@@ -1,5 +1,5 @@ Name:               hmatrix-Version:            0.8.2.0+Version:            0.8.3.1 License:            GPL License-file:       LICENSE Author:             Alberto Ruiz@@ -14,10 +14,11 @@ tested-with:        GHC ==6.10.4, GHC ==6.12.1  cabal-version:      >=1.2+ build-type:         Custom-extra-source-files: lib/Numeric/LinearAlgebra/Tests/quickCheckCompat.h -extra-source-files: configure configure.hs README INSTALL CHANGES+extra-source-files: lib/Numeric/LinearAlgebra/Tests/quickCheckCompat.h+                    configure configure.hs README INSTALL CHANGES extra-tmp-files:    hmatrix.buildinfo  extra-source-files: examples/tests.hs@@ -25,6 +26,7 @@                     examples/integrate.hs                     examples/minimize.hs                     examples/root.hs+                    examples/ode.hs                     examples/pca1.hs                     examples/pca2.hs                     examples/pinv.hs@@ -33,12 +35,11 @@                     examples/kalman.hs                     examples/parallel.hs                     examples/plot.hs-                    examples/latexmat.hs                     examples/inplace.hs-                    examples/benchmarks.hs                     examples/error.hs                     examples/devel/wrappers.hs                     examples/devel/functions.c+                    examples/Real.hs  extra-source-files: lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h,                     lib/Numeric/LinearAlgebra/LAPACK/clapack.h@@ -47,28 +48,19 @@                     lib/Numeric/GSL/Special/autoall.sh,                     lib/Numeric/GSL/Special/replace.hs -flag splitBase-    description:    Choose the new smaller, split-up base package.- flag mkl     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  library-    if flag(splitBase)-      build-depends:    base >= 3 && < 5, array-    else-      build-depends:    base < 3 -    Build-Depends:      haskell98,+    Build-Depends:      base >= 3 && < 5,+                        array,+                        haskell98,                         QuickCheck, HUnit,                         storable-complex,                         process@@ -86,6 +78,7 @@                         Numeric.GSL.Polynomials,                         Numeric.GSL.Minimization,                         Numeric.GSL.Root,+                        Numeric.GSL.ODE,                         Numeric.GSL.Vector,                         Numeric.GSL.Special,                         Numeric.GSL.Special.Gamma,@@ -127,7 +120,7 @@                         Numeric.LinearAlgebra.Tests,                         Data.Packed.Convert,                         Data.Packed.ST,-                        Data.Packed.Development+                        Data.Packed.Development,                         Data.Packed.Random     other-modules:      Data.Packed.Internal,                         Data.Packed.Internal.Common,@@ -159,12 +152,19 @@       else         extra-libraries:   gsl mkl_lapack mkl_intel mkl_sequential mkl_core -    if flag(accelerate)-        frameworks: Accelerate+    if os(OSX)+        extra-lib-dirs: /opt/local/lib/+        include-dirs: /opt/local/include         extra-libraries: gsl+        frameworks: Accelerate  -- The extra-libraries required for GSL and LAPACK -- should now be automatically detected by configure(.hs)      extra-libraries:     extra-lib-dirs:++    source-repository head+        type:     darcs+        location: http://code.haskell.org/hmatrix+

@@ -18,8 +18,9 @@   Adapt,   app1, app2, app3, app4,   (//), check,-  partit, common, compatdim,-  fi+  splitEvery, common, compatdim,+  fi,+  table ) where  import Foreign@@ -27,11 +28,12 @@ import Foreign.C.String(peekCString) import Foreign.C.Types import Foreign.Storable.Complex()+import Data.List(transpose,intersperse) --- | @partit 3 [1..9] == [[1,2,3],[4,5,6],[7,8,9]]@-partit :: Int -> [a] -> [[a]]-partit _ [] = []-partit n l  = take n l : partit n (drop n l)+-- | @splitEvery 3 [1..9] == [[1,2,3],[4,5,6],[7,8,9]]@+splitEvery :: Int -> [a] -> [[a]]+splitEvery _ [] = []+splitEvery k l = take k l : splitEvery k (drop k l)  -- | obtains the common value of a property of a list common :: (Eq a) => (b->a) -> [b] -> Maybe a@@ -46,6 +48,15 @@ compatdim [] = Nothing compatdim [a] = Just a compatdim (a:b:xs) = if a==b || a==1 || b==1 then compatdim (max a b:xs) else Nothing++-- | Formatting tool+table :: String -> [[String]] -> String+table sep as = unlines . 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 sep  -- | postfix function application (@flip ($)@) (//) :: x -> (x -> y) -> y

@@ -16,7 +16,7 @@ -- #hide  module Data.Packed.Internal.Matrix(-    Matrix(..),+    Matrix(..), rows, cols,     MatrixOrder(..), orderOf,     createMatrix, withMatrix, mat,     cmat, fmat,@@ -77,17 +77,23 @@ data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq)  -- | Matrix representation suitable for GSL and LAPACK computations.-data Matrix t = MC { rows :: {-# UNPACK #-} !Int-                   , cols :: {-# UNPACK #-} !Int+data Matrix t = MC { irows :: {-# UNPACK #-} !Int+                   , icols :: {-# UNPACK #-} !Int                    , cdat :: {-# UNPACK #-} !(Vector t) } -              | MF { rows :: {-# UNPACK #-} !Int-                   , cols :: {-# UNPACK #-} !Int+              | MF { irows :: {-# UNPACK #-} !Int+                   , icols :: {-# UNPACK #-} !Int                    , fdat :: {-# UNPACK #-} !(Vector t) }  -- MC: preferred by C, fdat may require a transposition -- MF: preferred by LAPACK, cdat may require a transposition +rows :: Matrix t -> Int+rows = irows++cols :: Matrix t -> Int+cols = icols+ xdat MC {cdat = d } = d xdat MF {fdat = d } = d @@ -97,16 +103,16 @@  -- | Matrix transpose. trans :: Matrix t -> Matrix t-trans MC {rows = r, cols = c, cdat = d } = MF {rows = c, cols = r, fdat = d }-trans MF {rows = r, cols = c, fdat = d } = MC {rows = c, cols = r, cdat = d }+trans MC {irows = r, icols = c, cdat = d } = MF {irows = c, icols = r, fdat = d }+trans MF {irows = r, icols = c, fdat = d } = MC {irows = c, icols = r, cdat = d }  cmat :: (Element t) => Matrix t -> Matrix t cmat m@MC{} = m-cmat MF {rows = r, cols = c, fdat = d } = MC {rows = r, cols = c, cdat = transdata r d c}+cmat MF {irows = r, icols = c, fdat = d } = MC {irows = r, icols = c, cdat = transdata r d c}  fmat :: (Element t) => Matrix t -> Matrix t fmat m@MF{} = m-fmat MC {rows = r, cols = c, cdat = d } = MF {rows = r, cols = c, fdat = transdata c d r}+fmat MC {irows = r, icols = c, cdat = d } = MF {irows = r, icols = c, fdat = transdata c d r}  -- C-Haskell matrix adapter mat :: Adapt (CInt -> CInt -> Ptr t -> r) (Matrix t) r@@ -133,7 +139,7 @@  -- | the inverse of 'Data.Packed.Matrix.fromLists' toLists :: (Element t) => Matrix t -> [[t]]-toLists m = partit (cols m) . toList . flatten $ m+toLists m = splitEvery (cols m) . toList . flatten $ m  -- | Create a matrix from a list of vectors. -- All vectors must have the same dimension,@@ -170,13 +176,13 @@ --    | i<0 || i>=r || j<0 || j>=c = error "matrix indexing out of range" --    | otherwise   = cdat m `at` (i*c+j) -MC {rows = r, cols = c, cdat = v} @@> (i,j)+MC {irows = r, icols = c, cdat = v} @@> (i,j)     | safe      = if i<0 || i>=r || j<0 || j>=c                     then error "matrix indexing out of range"                     else v `at` (i*c+j)     | otherwise = v `at` (i*c+j) -MF {rows = r, cols = c, fdat = v} @@> (i,j)+MF {irows = r, icols = c, fdat = v} @@> (i,j)     | safe      = if i<0 || i>=r || j<0 || j>=c                     then error "matrix indexing out of range"                     else v `at` (j*r+i)@@ -184,18 +190,18 @@ {-# 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)+atM' MC {icols = c, cdat = v} i j = v `at'` (i*c+j)+atM' MF {irows = r, fdat = v} i j = v `at'` (j*r+i) {-# INLINE atM' #-}  ------------------------------------------------------------------ -matrixFromVector RowMajor c v = MC { rows = r, cols = c, cdat = v }+matrixFromVector RowMajor c v = MC { irows = r, icols = c, cdat = v }     where (d,m) = dim v `divMod` c           r | m==0 = d             | otherwise = error "matrixFromVector" -matrixFromVector ColumnMajor c v = MF { rows = r, cols = c, fdat = v }+matrixFromVector ColumnMajor c v = MF { irows = r, icols = c, fdat = v }     where (d,m) = dim v `divMod` c           r | m==0 = d             | otherwise = error "matrixFromVector"@@ -223,8 +229,8 @@  -- | application of a vector function on the flattened matrix elements liftMatrix :: (Element a, Element b) => (Vector a -> Vector b) -> Matrix a -> Matrix b-liftMatrix f MC { cols = c, cdat = d } = matrixFromVector RowMajor    c (f d)-liftMatrix f MF { cols = c, fdat = d } = matrixFromVector ColumnMajor c (f d)+liftMatrix f MC { icols = c, cdat = d } = matrixFromVector RowMajor    c (f d)+liftMatrix f MF { icols = c, fdat = d } = matrixFromVector ColumnMajor c (f d)  -- | application of a vector function on the flattened matrices elements liftMatrix2 :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t

@@ -15,7 +15,7 @@ -- #hide  module Data.Packed.Internal.Vector (-    Vector(..),+    Vector(..), dim,     fromList, toList, (|>),     join, (@>), safe, at, at', subVector,     mapVector, zipVector,@@ -47,10 +47,14 @@  -- | A one-dimensional array of objects stored in a contiguous memory block. data Vector t =-    V { dim  :: {-# UNPACK #-} !Int               -- ^ number of elements+    V { idim  :: {-# UNPACK #-} !Int              -- ^ number of elements       , fptr :: {-# UNPACK #-} !(ForeignPtr t)    -- ^ foreign pointer to the memory block       } +-- | Number of elements+dim :: Vector t -> Int+dim = idim+ -- C-Haskell vector adapter vec :: Adapt (CInt -> Ptr t -> r) (Vector t) r vec = withVector@@ -157,7 +161,7 @@                         -> Int       -- ^ number of elements to extract                         -> Vector t  -- ^ source                         -> Vector t  -- ^ result-subVector k l (v@V {dim=n})+subVector k l (v@V {idim=n})     | k<0 || k >= n || k+l > n || l < 0 = error "subVector out of range"     | otherwise = unsafePerformIO $ do         r <- createVector l@@ -192,23 +196,23 @@         joiner as tot ptr     return r   where joiner [] _ _ = return ()-        joiner (V {dim = n, fptr = b} : cs) _ p = do+        joiner (V {idim = n, fptr = b} : cs) _ p = do             withForeignPtr b $ \pb -> copyArray p pb n             joiner cs 0 (advancePtr p n)   -- | transforms a complex vector into a real vector with alternating real and imaginary parts  asReal :: Vector (Complex Double) -> Vector Double-asReal v = V { dim = 2*dim v, fptr =  castForeignPtr (fptr v) }+asReal v = V { idim = 2*dim v, fptr =  castForeignPtr (fptr v) }  -- | transforms a real vector into a complex vector with alternating real and imaginary parts asComplex :: Vector Double -> Vector (Complex Double)-asComplex v = V { dim = dim v `div` 2, fptr =  castForeignPtr (fptr v) }+asComplex v = V { idim = dim v `div` 2, fptr =  castForeignPtr (fptr v) }  ----------------------------------------------------------------  cloneVector :: Storable t => Vector t -> IO (Vector t)-cloneVector (v@V {dim=n}) = do+cloneVector (v@V {idim=n}) = do         r <- createVector n         let f _ s _ d =  copyArray d s n >> return 0         app2 f vec v vec r "cloneVector"

@@ -27,8 +27,8 @@     subMatrix, takeRows, dropRows, takeColumns, dropColumns,     extractRows,     ident, diag, diagRect, takeDiag,-    liftMatrix, liftMatrix2,-    format, dispf, disps,+    liftMatrix, liftMatrix2, liftMatrix2Auto,+    dispf, disps, dispcf, latexFormat, format,     loadMatrix, saveMatrix, fromFile, fileDimensions,     readMatrix, fromArray2D ) where@@ -36,10 +36,11 @@ import Data.Packed.Internal import qualified Data.Packed.ST as ST import Data.Packed.Vector-import Data.List(transpose,intersperse) import Data.Array import System.Process(readProcess) import Text.Printf(printf)+import Data.List(transpose,intersperse)+import Data.Complex  -- | creates a matrix from a vertical list of matrices joinVert :: Element t => [Matrix t] -> Matrix t@@ -85,7 +86,7 @@     rs = map (compatdim . map rows) ms     cs = map (compatdim . map cols) (transpose ms)     szs = sequence [rs,cs]-    ms' = partit bc $ zipWith g szs (concat ms)+    ms' = splitEvery bc $ zipWith g szs (concat ms)      g [Just nr,Just nc] m                 | nr == r && nc == c = m@@ -232,30 +233,10 @@           r = r1-r0+1           c = c1-c0+1 --------------------------------------------------------{---- shows a Double with n digits after the decimal point    -shf :: (RealFloat a) => Int -> a -> String     -shf dec n | abs n < 1e-10 = "0."-          | abs (n - (fromIntegral.round $ n)) < 1e-10 = show (round n) ++"."-          | otherwise = showGFloat (Just dec) n ""    --- shows a Complex Double as a pair, with n digits after the decimal point    -shfc n z@ (a:+b) -    | magnitude z <1e-10 = "0."-    | abs b < 1e-10 = shf n a-    | abs a < 1e-10 = shf n b ++"i"-    | b > 0         = shf n a ++"+"++shf n b ++"i"-    | otherwise     = shf n a ++shf n b ++"i"          --}+-------------------------------------------------------------------+-- display utilities -dsp' :: String -> [[String]] -> String-dsp' sep as = unlines . 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 sep  {- | 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:@@ -266,20 +247,7 @@  -} format :: (Element t) => String -> (t -> String) -> Matrix t -> String-format sep f m = dsp' sep . map (map f) . toLists $ m--{--disp m f = putStrLn $ "matrix ("++show (rows m) ++"x"++ show (cols m) ++")\n"++format " | " f m--dispR :: Int -> Matrix Double -> IO ()-dispR d m = disp m (shf d)--dispC :: Int -> Matrix (Complex Double) -> IO ()-dispC d m = disp m (shfc d)--}------------------------------------------------------------------------ display utilities+format sep f m = table sep . map (map f) . toLists $ m  {- | Show a matrix with \"autoscaling\" and a given number of decimal places. @@ -314,9 +282,7 @@  formatFixed d x = format "  " (printf ("%."++show d++"f")) $ x -isInt = all lookslikeInt . toList . flatten where-    lookslikeInt x = show (round x :: Int) ++".0" == shx || "-0.0" == shx-        where shx = show x+isInt = all lookslikeInt . toList . flatten  formatScaled dec t = "E"++show o++"\n" ++ ss     where ss = format " " (printf fmt. g) t@@ -340,7 +306,43 @@     . 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.@@ -385,4 +387,29 @@  -} repmat :: (Element t) => Matrix t -> Int -> Int -> Matrix t-repmat m r c = fromBlocks $ partit c $ replicate (r*c) m+repmat m r c = fromBlocks $ splitEvery c $ replicate (r*c) m++-- | A version of 'liftMatrix2' which automatically adapt matrices with a single row or column to match the dimensions of the other matrix.+liftMatrix2Auto :: (Element t, Element a, Element b)+             => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t+liftMatrix2Auto f m1 m2 | compat' m1 m2 = lM f m1 m2+                     | rows m1 == rows m2 && cols m2 == 1 = lM f m1 (repCols (cols m1) m2)+                     | rows m1 == rows m2 && cols m1 == 1 = lM f (repCols (cols m2) m1) m2+                     | cols m1 == cols m2 && rows m2 == 1 = lM f m1 (repRows (rows m1) m2)+                     | cols m1 == cols m2 && cols m1 == 1 = lM f (repRows (rows m2) m1) m2+                     | rows m1 == 1 && cols m2 == 1 = lM f (repRows (rows m2) m1) (repCols (cols m1) m2)+                     | cols m1 == 1 && rows m2 == 1 = lM f (repCols (cols m2) m1) (repRows (rows m1) m2)+                     | otherwise    = error $ "nonconformable matrices in liftMatrix2Auto: " ++ show (size m1) ++ ", " ++ show (size m2)++size m = (rows m, cols m)++lM f m1 m2 = reshape (max (cols m1) (cols m2)) (f (flatten m1) (flatten m2))++repRows n x = fromRows (replicate n (flatten x))+repCols n x = fromColumns (replicate n (flatten x))++compat' :: Matrix a -> Matrix b -> Bool+compat' m1 m2 = rows m1 == 1 && cols m1 == 1+             || rows m2 == 1 && cols m2 == 1+             || rows m1 == rows m2 && cols m1 == cols m2+

@@ -19,6 +19,7 @@ , module Numeric.GSL.Polynomials , module Numeric.GSL.Minimization , module Numeric.GSL.Root+, module Numeric.GSL.ODE , module Numeric.GSL.Special , module Complex , setErrorHandlerOff@@ -31,6 +32,7 @@ import Numeric.GSL.Polynomials import Numeric.GSL.Minimization import Numeric.GSL.Root+import Numeric.GSL.ODE import Complex  

@@ -30,6 +30,9 @@ foreign import ccall "wrapper"     mkVecVecfun :: TVV -> IO (FunPtr TVV) +foreign import ccall "wrapper"+    mkDoubleVecVecfun :: (Double -> TVV) -> IO (FunPtr (Double -> TVV))+ aux_vTov :: (Vector Double -> Vector Double) -> TVV aux_vTov f n p nr r = g where     V {fptr = pr} = f x@@ -42,6 +45,9 @@  foreign import ccall "wrapper"     mkVecMatfun :: TVM -> IO (FunPtr TVM)++foreign import ccall "wrapper"+    mkDoubleVecMatfun :: (Double -> TVM) -> IO (FunPtr (Double -> TVM))  aux_vTom :: (Vector Double -> Matrix Double) -> TVM aux_vTom f n p rr cr r = g where

@@ -0,0 +1,111 @@+{- |+Module      :  Numeric.GSL.ODE+Copyright   :  (c) Alberto Ruiz 2010+License     :  GPL++Maintainer  :  Alberto Ruiz (aruiz at um dot es)+Stability   :  provisional+Portability :  uses ffi++Solution of ordinary differential equation (ODE) initial value problems.++<http://www.gnu.org/software/gsl/manual/html_node/Ordinary-Differential-Equations.html>++A simple example:++@import Numeric.GSL+import Numeric.LinearAlgebra+import Graphics.Plot++xdot t [x,v] = [v, -0.95*x - 0.1*v]++ts = linspace 100 (0,20)++sol = odeSolve xdot [10,0] ts++main = mplot (ts : toColumns sol)@++-}+-----------------------------------------------------------------------------++module Numeric.GSL.ODE (+    odeSolve, odeSolveV, ODEMethod(..)+) where++import Data.Packed.Internal+import Data.Packed.Matrix+import Foreign+import Foreign.C.Types(CInt)+import Numeric.GSL.Internal++-------------------------------------------------------------------------++-- | Stepping functions+data ODEMethod = RK2 -- ^ Embedded Runge-Kutta (2, 3) method.+               | RK4 -- ^ 4th order (classical) Runge-Kutta. The error estimate is obtained by halving the step-size. For more efficient estimate of the error, use 'RKf45'.+               | RKf45 -- ^ Embedded Runge-Kutta-Fehlberg (4, 5) method. This method is a good general-purpose integrator.+               | RKck -- ^ Embedded Runge-Kutta Cash-Karp (4, 5) method.+               | RK8pd -- ^ Embedded Runge-Kutta Prince-Dormand (8,9) method.+               | RK2imp -- ^ Implicit 2nd order Runge-Kutta at Gaussian points.+               | RK4imp -- ^ Implicit 4th order Runge-Kutta at Gaussian points.+               | BSimp -- ^ Implicit Bulirsch-Stoer method of Bader and Deuflhard. This algorithm requires the Jacobian.+               | Gear1 -- ^ M=1 implicit Gear method.+               | Gear2 -- ^ M=2 implicit Gear method.+               deriving (Enum,Eq,Show,Bounded)++-- | A version of 'odeSolveV' with reasonable default parameters and system of equations defined using lists.+odeSolve+    :: (Double -> [Double] -> [Double])        -- ^ xdot(t,x)+    -> [Double]        -- ^ initial conditions+    -> Vector Double   -- ^ desired solution times+    -> Matrix Double   -- ^ solution+odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) Nothing (fromList xi) ts+    where hi = (ts@>1 - ts@>0)/100+          epsAbs = 1.49012e-08+          epsRel = 1.49012e-08+          l2v f = \t -> fromList  . f t . toList+          l2m f = \t -> fromLists . f t . toList++-- | Evolution of the system with adaptive step-size control.+odeSolveV+    :: ODEMethod+    -> Double -- ^ initial step size+    -> Double -- ^ absolute tolerance for the state vector+    -> Double -- ^ relative tolerance for the state vector+    -> (Double -> Vector Double -> Vector Double)   -- ^ xdot(t,x)+    -> Maybe (Double -> Vector Double -> Matrix Double)   -- ^ optional jacobian+    -> Vector Double     -- ^ initial conditions+    -> Vector Double     -- ^ desired solution times+    -> Matrix Double     -- ^ solution+odeSolveV method h epsAbs epsRel f mbjac xiv ts = unsafePerformIO $ do+    let n   = dim xiv+    fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t))+    jp <- case mbjac of+        Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t))+        Nothing  -> return nullFunPtr+    sol <- withVector xiv $ \xiv' ->+            withVector (checkTimes ts) $ \ts' ->+             createMIO (dim ts) n+              (ode_c (fi (fromEnum method)) h epsAbs epsRel fp jp // xiv' // ts' )+              "ode"+    freeHaskellFunPtr fp+    return sol++foreign import ccall "ode"+    ode_c :: CInt -> Double -> Double -> Double -> FunPtr (Double -> TVV) -> FunPtr (Double -> TVM) -> TVVM++-------------------------------------------------------++checkdim1 n v+    | dim v == n = v+    | otherwise = error $ "Error: "++ show n+                        ++ " components expected in the result of the function supplied to odeSolve"++checkdim2 n m+    | rows m == n && cols m == n = m+    | otherwise = error $ "Error: "++ show n ++ "x" ++ show n+                        ++ " Jacobian expected in odeSolve"++checkTimes ts | dim ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts+              | otherwise = error "odeSolve requires increasing times"+    where ts' = toList ts

@@ -22,6 +22,7 @@ #include <gsl/gsl_complex_math.h> #include <gsl/gsl_rng.h> #include <gsl/gsl_randist.h>+#include <gsl/gsl_odeiv.h> #include <string.h> #include <stdio.h> @@ -507,7 +508,7 @@   void df_aux_min(const gsl_vector * x, void * pars, gsl_vector * g) {-    Tfdf * fdf = ((Tfdf*) pars);  +    Tfdf * fdf = ((Tfdf*) pars);     double* p = (double*)calloc(x->size,sizeof(double));     double* q = (double*)calloc(g->size,sizeof(double));     int k;@@ -797,3 +798,89 @@     } } #undef RAN++//////////////////////////////////////////////////////+//                        ODE++typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode;++int odefunc (double t, const double y[], double f[], void *params) { +    Tode * P = (Tode*) params;+    (P->f)(t,P->n,y,P->n,f);+    return GSL_SUCCESS;+}++int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) {+     Tode * P = ((Tode*) params);+     (P->j)(t,P->n,y,P->n,P->n,dfdy);+     int j;+     for (j=0; j< P->n; j++)+        dfdt[j] = 0.0;+     return GSL_SUCCESS;+}+++int ode(int method, double h, double eps_abs, double eps_rel,+        int f(double, int, const double*, int, double*),+        int jac(double, int, const double*, int, int, double*),+        KRVEC(xi), KRVEC(ts), RMAT(sol)) {++    const gsl_odeiv_step_type * T;++    switch(method) {+        case 0 : {T = gsl_odeiv_step_rk2; break; }+        case 1 : {T = gsl_odeiv_step_rk4; break; }+        case 2 : {T = gsl_odeiv_step_rkf45; break; }+        case 3 : {T = gsl_odeiv_step_rkck; break; }+        case 4 : {T = gsl_odeiv_step_rk8pd; break; }+        case 5 : {T = gsl_odeiv_step_rk2imp; break; }+        case 6 : {T = gsl_odeiv_step_rk4imp; break; }+        case 7 : {T = gsl_odeiv_step_bsimp; break; }+        case 8 : {T = gsl_odeiv_step_gear1; break; }+        case 9 : {T = gsl_odeiv_step_gear2; break; }+        default: ERROR(BAD_CODE);+    }+++    gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, xin);+    gsl_odeiv_control * c = gsl_odeiv_control_y_new (eps_abs, eps_rel);+    gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (xin);++    Tode P;+    P.f = f;+    P.j = jac;+    P.n = xin;++    gsl_odeiv_system sys = {odefunc, odejac, xin, &P};++    double t = tsp[0];++    double* y = (double*)calloc(xin,sizeof(double));+    int i,j;+    for(i=0; i< xin; i++) {+        y[i] = xip[i];+        solp[i] = xip[i];+    }++       for (i = 1; i < tsn ; i++)+         {+           double ti = tsp[i];+           while (t < ti)+             {+               gsl_odeiv_evolve_apply (e, c, s,+                                       &sys,+                                       &t, ti, &h,+                                       y);+               // if (h < hmin) h = hmin;+             }+           for(j=0; j<xin; j++) {+               solp[i*xin + j] = y[j];+           }+         }++    free(y);+    gsl_odeiv_evolve_free (e);+    gsl_odeiv_control_free (c);+    gsl_odeiv_step_free (s);+    return 0;+}

@@ -729,7 +729,7 @@ -} kronecker :: (Field t) => Matrix t -> Matrix t -> Matrix t kronecker a b = fromBlocks-              . partit (cols a)+              . splitEvery (cols a)               . map (reshape (cols b))               . toRows               $ flatten a `outer` flatten b

@@ -71,15 +71,7 @@     | dim y == 1 = f3 x (y@>0)     | otherwise = f2 x y -liftMatrix2' :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t-liftMatrix2' f m1 m2 | compat' m1 m2 = reshape (max (cols m1) (cols m2)) (f (flatten m1) (flatten m2))-                     | otherwise    = error "nonconformant matrices in liftMatrix2'" -compat' :: Matrix a -> Matrix b -> Bool-compat' m1 m2 = rows m1 == 1 && cols m1 == 1-             || rows m2 == 1 && cols m2 == 1-             || rows m1 == rows m2 && cols m1 == cols m2- instance Linear Vector a => Eq (Vector a) where     (==) = equal @@ -103,10 +95,10 @@     (==) = equal  instance (Linear Matrix a, Num (Vector a)) => Num (Matrix a) where-    (+) = liftMatrix2' (+)-    (-) = liftMatrix2' (-)+    (+) = liftMatrix2Auto (+)+    (-) = liftMatrix2Auto (-)     negate = liftMatrix negate-    (*) = liftMatrix2' (*)+    (*) = liftMatrix2Auto (*)     signum = liftMatrix signum     abs = liftMatrix abs     fromInteger = (1><1) . return . fromInteger@@ -123,7 +115,7 @@  instance (Linear Vector a, Fractional (Vector a), Num (Matrix a)) => Fractional (Matrix a) where     fromRational n = (1><1) [fromRational n]-    (/) = liftMatrix2' (/)+    (/) = liftMatrix2Auto (/)  --------------------------------------------------------- @@ -184,14 +176,14 @@     atanh = liftMatrix atanh     exp   = liftMatrix exp     log   = liftMatrix log-    (**)  = liftMatrix2' (**)+    (**)  = liftMatrix2Auto (**)     sqrt  = liftMatrix sqrt     pi    = (1><1) [pi]  ---------------------------------------------------------------  instance (Storable a, Num (Vector a)) => Monoid (Vector a) where-    mempty = 0 { dim = 0 }+    mempty = 0 { idim = 0 }     mappend a b = mconcat [a,b]     mconcat = j . filter ((>0).dim)         where j [] = mempty@@ -204,3 +196,4 @@ -- -- instance (NFData a, Element a) => NFData (Matrix a) where --     rnf = rnf . flatten+

@@ -72,6 +72,9 @@  ------------------------------------------------ +{-# 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@@ -108,3 +111,4 @@ -- | Vertical concatenation of matrices and vectors. (<->) :: (Element t, Joinable a b) => a t -> b t -> Matrix t a <-> b = joinV a b+

@@ -135,6 +135,14 @@  --------------------------------------------------------------------- +odeTest = utest "ode" (last (toLists sol) ~~ [-1.7588880332411019, 8.364348908711941e-2])+    where sol = odeSolveV RK8pd 1E-6 1E-6 0 (l2v $ vanderpol 10) Nothing (fromList [1,0]) ts+          ts = linspace 101 (0,100)+          l2v f = \t -> fromList  . f t . toList+          vanderpol mu _t [x,y] = [y, -x + mu * y * (1-x^2) ]++---------------------------------------------------------------------+ randomTestGaussian = c :~1~: snd (meanCov dat) where     a = (3><3) [1,2,3,                 2,4,0,@@ -280,6 +288,7 @@         , utest "rank" $  rank ((2><3)[1,0,0,1,6*eps,0]) == 1                        && rank ((2><3)[1,0,0,1,7*eps,0]) == 2         , utest "block" $ fromBlocks [[ident 3,0],[0,ident 4]] == (ident 7 :: CM)+        , odeTest         ]     return () @@ -366,3 +375,4 @@     time "full svd         3000x500" (fv $ svd a)     time "singular values 1000x1000" (singularValues b)     time "full svd        1000x1000" (fv $ svd b)+

@@ -148,7 +148,7 @@            && orthonormal u && orthonormal v            && (dim s == r || r == 0 && dim s == 1)     where (u,s,v) = compactSVD m-          m = m' <-> m'+          m = fromBlocks [[m'],[m']]           r = rank m'  svdProp5a m = and (map (s1|~|) [s2,s3,s4,s5,s6]) where