compatible with vector

7 files changed, 35 insertions, 17 deletions

diff --git a/lib/Numeric/GSL/Internal.hs b/lib/Numeric/GSL/Internal.hs
index 37bcc1b..2ccc3ef 100644
--- a/lib/Numeric/GSL/Internal.hs
+++ b/lib/Numeric/GSL/Internal.hs
@@ -35,12 +35,12 @@ foreign import ccall "wrapper"
 
 aux_vTov :: (Vector Double -> Vector Double) -> TVV
 aux_vTov f n p nr r = g where
-    V {fptr = pr} = f x
+    v = f x
     x = createV (fromIntegral n) copy "aux_vTov"
     copy n' q = do
         copyArray q p (fromIntegral n')
         return 0
-    g = do withForeignPtr pr $ \p' -> copyArray r p' (fromIntegral nr)
+    g = do unsafeWith v $ \p' -> copyArray r p' (fromIntegral nr)
            return 0
 
 foreign import ccall "wrapper"
@@ -51,12 +51,12 @@ foreign import ccall "wrapper"
 
 aux_vTom :: (Vector Double -> Matrix Double) -> TVM
 aux_vTom f n p rr cr r = g where
-    V {fptr = pr} = flatten $ f x
+    v = flatten $ f x
     x = createV (fromIntegral n) copy "aux_vTov"
     copy n' q = do
         copyArray q p (fromIntegral n')
         return 0
-    g = do withForeignPtr pr $ \p' -> copyArray r p' (fromIntegral $ rr*cr)
+    g = do unsafeWith v $ \p' -> copyArray r p' (fromIntegral $ rr*cr)
            return 0
 
 createV n fun msg = unsafePerformIO $ do
diff --git a/lib/Numeric/GSL/Minimization.hs b/lib/Numeric/GSL/Minimization.hs
index 3240ddf..3bbb2b3 100644
--- a/lib/Numeric/GSL/Minimization.hs
+++ b/lib/Numeric/GSL/Minimization.hs
@@ -109,7 +109,7 @@ ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2
 minimizeV method eps maxit szv f xiv = unsafePerformIO $ do
     let n   = dim xiv
     fp <- mkVecfun (iv f)
-    rawpath <- ww2 withVector xiv withVector szv $ \xiv' szv' ->
+    rawpath <- ww2 vec xiv vec szv $ \xiv' szv' ->
                    createMIO maxit (n+3)
                          (c_minimize (fi (fromEnum method)) fp eps (fi maxit) // xiv' // szv')
                          "minimize"
@@ -166,7 +166,7 @@ minimizeVD method eps maxit istep tol f df xiv = unsafePerformIO $ do
         df' = (checkdim1 n . df)
     fp <- mkVecfun (iv f')
     dfp <- mkVecVecfun (aux_vTov df')
-    rawpath <- withVector xiv $ \xiv' ->
+    rawpath <- vec xiv $ \xiv' ->
                     createMIO maxit (n+2)
                          (c_minimizeD (fi (fromEnum method)) fp dfp istep tol eps (fi maxit) // xiv')
                          "minimizeD"
diff --git a/lib/Numeric/GSL/ODE.hs b/lib/Numeric/GSL/ODE.hs
index 9d890fc..f6f11f9 100644
--- a/lib/Numeric/GSL/ODE.hs
+++ b/lib/Numeric/GSL/ODE.hs
@@ -83,8 +83,8 @@ odeSolveV method h epsAbs epsRel f mbjac xiv ts = unsafePerformIO $ do
     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' ->
+    sol <- vec xiv $ \xiv' ->
+            vec (checkTimes ts) $ \ts' ->
              createMIO (dim ts) n
               (ode_c (fi (fromEnum method)) h epsAbs epsRel fp jp // xiv' // ts' )
               "ode"
diff --git a/lib/Numeric/GSL/Root.hs b/lib/Numeric/GSL/Root.hs
index 840e8ee..ddb090d 100644
--- a/lib/Numeric/GSL/Root.hs
+++ b/lib/Numeric/GSL/Root.hs
@@ -79,7 +79,7 @@ rootGen m f xi epsabs maxit = unsafePerformIO $ do
     let xiv = fromList xi
         n   = dim xiv
     fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
-    rawpath <- withVector xiv $ \xiv' ->
+    rawpath <- vec xiv $ \xiv' ->
                    createMIO maxit (2*n+1)
                          (c_root m fp epsabs (fi maxit) // xiv')
                          "root"
@@ -117,7 +117,7 @@ rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do
         n   = dim xiv
     fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
     jp <- mkVecMatfun (aux_vTom (checkdim2 n . fromLists . jac . toList))
-    rawpath <- withVector xiv $ \xiv' ->
+    rawpath <- vec xiv $ \xiv' ->
                    createMIO maxit (2*n+1)
                          (c_rootj m fp jp epsabs (fi maxit) // xiv')
                          "root"
diff --git a/lib/Numeric/LinearAlgebra/Instances.hs b/lib/Numeric/LinearAlgebra/Instances.hs
index 67496f2..b17da5a 100644
--- a/lib/Numeric/LinearAlgebra/Instances.hs
+++ b/lib/Numeric/LinearAlgebra/Instances.hs
@@ -25,7 +25,7 @@ import Data.Packed.Matrix
 import Complex
 import Data.List(transpose,intersperse)
 import Foreign(Storable)
-import Data.Monoid
+-- import Data.Monoid
 import Data.Packed.Internal.Vector
 -- import Control.Parallel.Strategies
 
@@ -182,12 +182,12 @@ instance (Linear Vector a, Floating (Vector a), Fractional (Matrix a)) => Floati
 
 ---------------------------------------------------------------
 
-instance (Storable a, Num (Vector a)) => Monoid (Vector a) where
-    mempty = 0 { idim = 0 }
-    mappend a b = mconcat [a,b]
-    mconcat = j . filter ((>0).dim)
-        where j [] = mempty
-              j l  = join l
+-- instance (Storable a, Num (Vector a)) => Monoid (Vector a) where
+--     mempty = 0 { idim = 0 }
+--     mappend a b = mconcat [a,b]
+--     mconcat = j . filter ((>0).dim)
+--         where j [] = mempty
+--               j l  = join l
 
 ---------------------------------------------------------------
 
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs
index 016b9a1..3f7c847 100644
--- a/lib/Numeric/LinearAlgebra/Tests.hs
+++ b/lib/Numeric/LinearAlgebra/Tests.hs
@@ -180,6 +180,9 @@ runTests n = do
     test (multProp1  . cConsist)
     test (multProp2  . rConsist)
     test (multProp2  . cConsist)
+    putStrLn "------ sub-trans"
+    test (subProp . rM)
+    test (subProp . cM)
     putStrLn "------ lu"
     test (luProp    . rM)
     test (luProp    . cM)
@@ -305,6 +308,7 @@ makeUnitary v | realPart n > 1    = v / scalar n
 -- | Performance measurements.
 runBenchmarks :: IO ()
 runBenchmarks = do
+    subBench
     multBench
     svdBench
     eigBench
@@ -335,6 +339,16 @@ multb n = foldl1' (<>) (replicate (10^6) (ident n :: Matrix Double))
 
 --------------------------------
 
+subBench = do
+    putStrLn ""
+    let g = foldl1' (.) (replicate (10^5) (\v -> subVector 1 (dim v -1) v))
+    time "0.1M subVector   " (g (constant 1 (1+10^5) :: Vector Double) @> 0)
+    let f = foldl1' (.) (replicate (10^5) (fromRows.toRows))
+    time "subVector-join  3" (f (ident  3 :: Matrix Double) @@>(0,0))
+    time "subVector-join 10" (f (ident 10 :: Matrix Double) @@>(0,0))
+
+--------------------------------
+
 multBench = do
     let a = ident 1000 :: Matrix Double
     let b = ident 2000 :: Matrix Double
diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
index 93175f2..6d176fa 100644
--- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -38,6 +38,7 @@ module Numeric.LinearAlgebra.Tests.Properties (
     cholProp,
     expmDiagProp,
     multProp1, multProp2,
+    subProp,
     linearSolveProp, linearSolveProp2
 ) where
 
@@ -243,3 +244,6 @@ linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b)
     where q = min (rows a) (cols a)
           b = a <> x
           wc = rank a == q
+
+subProp m = m == (trans . fromColumns . toRows) m
+