improved constant

9 files changed, 38 insertions, 75 deletions

diff --git a/lib/Data/Packed.hs b/lib/Data/Packed.hs
index bb25c26..7d6d200 100644
--- a/lib/Data/Packed.hs
+++ b/lib/Data/Packed.hs
@@ -24,7 +24,7 @@ module Data.Packed (
 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 @@ class (Element e) => Container c e where
 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 @@ instance Container Matrix Double where
     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 @@ instance Container Matrix (Complex Double) where
     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))
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 51fb6f8..477b453 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -208,19 +208,16 @@ compat m1 m2 = rows m1 == rows m2 && cols m1 == cols m2
 
 -- | Optimized matrix computations are provided for elements in the Element 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
 
 instance Element Double where
-    constantD  = constantR
     transdata = transdataR
     subMatrixD = subMatrixR
 
 instance Element (Complex Double) where
-    constantD  = constantC
     transdata  = transdataC
     subMatrixD = subMatrixC
 
@@ -284,31 +281,6 @@ subMatrix :: Element a
           -> Matrix a -- ^ result
 subMatrix = subMatrixD
 
-------------------------------------------------------------------------
-
-constantAux fun x n = unsafePerformIO $ do
-    v <- createVector n
-    px <- newArray [x]
-    app1 (fun px) vec v "constantAux"
-    free px
-    return v
-
-constantR :: Double -> Int -> Vector Double
-constantR = constantAux cconstantR
-foreign import ccall "auxi.h 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
-
-{- | 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 = constantD
-
 --------------------------------------------------------------------------
 
 -- | obtains the complex conjugate of a complex vector
@@ -329,10 +301,6 @@ fromComplex :: Vector (Complex Double) -> (Vector Double, Vector Double)
 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
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
index 2df33e0..f590919 100644
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -99,6 +99,7 @@ n |> l = if length l == n then fromList l else error "|> with wrong size"
 -- | 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.
diff --git a/lib/Data/Packed/Internal/auxi.c b/lib/Data/Packed/Internal/auxi.c
index c449b9a..48b05e8 100644
--- a/lib/Data/Packed/Internal/auxi.c
+++ b/lib/Data/Packed/Internal/auxi.c
@@ -1,12 +1,5 @@
 #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>
 
@@ -92,27 +85,6 @@ int submatrixR(int r1, int r2, int c1, int c2, KRMAT(x),RMAT(r)) {
 }
 
 
-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 conjugate(KCVEC(x),CVEC(t)) {
     REQUIRES(xn==tn,BAD_SIZE);
     DEBUGMSG("conjugate");
diff --git a/lib/Data/Packed/Internal/auxi.h b/lib/Data/Packed/Internal/auxi.h
index b57fdaf..4698696 100644
--- a/lib/Data/Packed/Internal/auxi.h
+++ b/lib/Data/Packed/Internal/auxi.h
@@ -13,9 +13,6 @@
 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));
 
 const char * gsl_strerror (const int gsl_errno);
diff --git a/lib/Data/Packed/ST.hs b/lib/Data/Packed/ST.hs
index 91b7968..9c9c3b3 100644
--- a/lib/Data/Packed/ST.hs
+++ b/lib/Data/Packed/ST.hs
@@ -1,4 +1,4 @@
-{-# OPTIONS -XTypeOperators -XRank2Types  -XFlexibleContexts #-}
+{-# OPTIONS -XTypeOperators -XRank2Types  -XFlexibleContexts -XBangPatterns #-}
 
 -----------------------------------------------------------------------------
 -- |
@@ -23,8 +23,10 @@ module Data.Packed.ST (
     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 @@ readVector = safeIndexV unsafeReadVector
 writeVector :: Storable t => STVector s t -> Int -> t -> ST s ()
 writeVector = safeIndexV unsafeWriteVector
 
+{-# NOINLINE newUndefinedVector #-}
+newUndefinedVector :: Element t => Int -> ST s (STVector s t)
+newUndefinedVector = unsafeIOToST . fmap STVector . createVector
+
 {-# NOINLINE 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 @@ readMatrix = safeIndexM unsafeReadMatrix
 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)
diff --git a/lib/Data/Packed/Vector.hs b/lib/Data/Packed/Vector.hs
index 6415c5c..0bbbc34 100644
--- a/lib/Data/Packed/Vector.hs
+++ b/lib/Data/Packed/Vector.hs
@@ -24,6 +24,7 @@ module Data.Packed.Vector (
 
 import Data.Packed.Internal
 import Numeric.GSL.Vector
+import Data.Packed.ST
 
 {- | Creates a real vector containing a range of values:
 
@@ -47,3 +48,11 @@ vectorMaxIndex = round . toScalarR MaxIdx
 
 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)
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index f259db5..d978fa4 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -61,7 +61,7 @@ module Numeric.LinearAlgebra.Algorithms (
 ) 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
diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs
index 56945d7..ffa2cb5 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK.hs
+++ b/lib/Numeric/LinearAlgebra/LAPACK.hs
@@ -27,11 +27,8 @@ module Numeric.LinearAlgebra.LAPACK (
     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