4 files changed, 17 insertions, 24 deletions

diff --git a/examples/Real.hs b/examples/Real.hs
index b32a961..02eba9f 100644
--- a/examples/Real.hs
+++ b/examples/Real.hs
@@ -67,13 +67,13 @@ diagl = diag . vector
 zeros :: Int -- ^ rows
       -> Int -- ^ columns
       -> Matrix Double
-zeros r c = reshape c (constant 0 (r*c))
+zeros r c = konst 0 (r,c)
 
 -- | Create a matrix or ones.
 ones :: Int -- ^ rows
      -> Int -- ^ columns
      -> Matrix Double
-ones r c = reshape c (constant 1 (r*c))
+ones r c = konst 1 (r,c)
 
 -- | Concatenation of real vectors.
 infixl 9 #
diff --git a/lib/Data/Packed/Random.hs b/lib/Data/Packed/Random.hs
index b30b299..c34f212 100644
--- a/lib/Data/Packed/Random.hs
+++ b/lib/Data/Packed/Random.hs
@@ -22,10 +22,8 @@ module Data.Packed.Random (
 import Numeric.GSL.Vector
 import Data.Packed
 import Numeric.Container
-import Data.Packed.Internal(constantD)
 import Numeric.LinearAlgebra.Algorithms
 
-constant k v = constantD k v
 
 -- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
 -- Gaussian distribution.
@@ -36,7 +34,7 @@ gaussianSample :: Int -- ^ seed
                -> Matrix Double -- ^ result
 gaussianSample seed n med cov = m where
     c = dim med
-    meds = constant 1 n `outer` med
+    meds = konst 1 n `outer` med
     rs = reshape c $ randomVector seed Gaussian (c * n)
     m = rs `mXm` cholSH cov `add` meds
 
@@ -52,7 +50,7 @@ uniformSample seed n rgs = m where
     cs = zipWith subtract as bs
     d = dim a
     dat = toRows $ reshape n $ randomVector seed Uniform (n*d)
-    am = constant 1 n `outer` a
+    am = konst 1 n `outer` a
     m = fromColumns (zipWith scale cs dat) `add` am
 
 ------------ utilities -------------------------------
@@ -62,7 +60,7 @@ meanCov :: Matrix Double -> (Vector Double, Matrix Double)
 meanCov x = (med,cov) where
     r    = rows x
     k    = 1 / fromIntegral r
-    med  = constant k r `vXm` x
-    meds = constant 1 r `outer` med
+    med  = konst k r `vXm` x
+    meds = konst 1 r `outer` med
     xc   = x `sub` meds
     cov  = flip scale (trans xc `mXm` xc) (recip (fromIntegral (r-1)))
diff --git a/lib/Numeric/Container.hs b/lib/Numeric/Container.hs
index 44e2711..2c9c500 100644
--- a/lib/Numeric/Container.hs
+++ b/lib/Numeric/Container.hs
@@ -71,9 +71,11 @@ class (Element e) => Container c e where
     -- | element by element division
     divide      :: c e -> c e -> c e
     equal       :: c e -> c e -> Bool
-
+    --
     -- | cannot implement instance Functor because of Element class constraint
     cmap        :: (Element a, Element b) => (a -> b) -> c a -> c b
+    -- | constant structure of given size
+    konst    :: e -> IndexOf c -> c e
     --
     -- | indexing function
     atIndex     :: c e -> IndexOf c -> e
@@ -107,8 +109,7 @@ instance Container Vector Float where
     divide = vectorZipF Div
     equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0
     scalar x = fromList [x]
-    --
---instance Container Vector Float where
+    konst = constantD
     cmap = mapVector
     atIndex = (@>)
     minIndex     = round . toScalarF MinIdx
@@ -128,8 +129,7 @@ instance Container Vector Double where
     divide = vectorZipR Div
     equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0
     scalar x = fromList [x]
-    --
---instance Container Vector Double where
+    konst = constantD
     cmap = mapVector
     atIndex = (@>)
     minIndex     = round . toScalarR MinIdx
@@ -149,8 +149,7 @@ instance Container Vector (Complex Double) where
     divide = vectorZipC Div
     equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
     scalar x = fromList [x]
-    --
---instance Container Vector (Complex Double) where
+    konst = constantD
     cmap = mapVector
     atIndex = (@>)
     minIndex     = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
@@ -170,8 +169,7 @@ instance Container Vector (Complex Float) where
     divide = vectorZipQ Div
     equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
     scalar x = fromList [x]
-    --
---instance Container Vector (Complex Float) where
+    konst = constantD
     cmap = mapVector
     atIndex = (@>)
     minIndex     = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
@@ -193,8 +191,7 @@ instance (Container Vector a) => Container Matrix a where
     divide = liftMatrix2 divide
     equal a b = cols a == cols b && flatten a `equal` flatten b
     scalar x = (1><1) [x]
-    --
---instance (Container Vector a) => Container Matrix a where
+    konst v (r,c) = reshape c (konst v (r*c))
     cmap f = liftMatrix (mapVector f)
     atIndex = (@@>)
     minIndex m = let (r,c) = (rows m,cols m)
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index 6e703fa..5bc9ca5 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -83,8 +83,6 @@ import Data.List(foldl1')
 import Data.Array
 import Numeric.Container
 
-constant x = constantD x
-
 -- | Auxiliary typeclass used to define generic computations for both real and complex matrices.
 class (Product t, Container Vector t, Container Matrix t) => Field t where
     svd'         :: Matrix t -> (Matrix t, Vector Double, Matrix t)
@@ -469,7 +467,7 @@ zh k v = fromList $ replicate (k-1) 0 ++ (1:drop k xs)
               where xs = toList v
 
 zt 0 v = v
-zt k v = join [subVector 0 (dim v - k) v, constant 0 k]
+zt k v = join [subVector 0 (dim v - k) v, konst 0 k]
 
 
 unpackQR :: (Field t) => (Matrix t, Vector t) -> (Matrix t, Matrix t)
@@ -623,10 +621,10 @@ luFact (l_u,perm) | r <= c    = (l ,u ,p, s)
     c = cols l_u
     tu = triang r c 0 1
     tl = triang r c 0 0
-    l = takeColumns r (l_u |*| tl) |+| diagRect (constant 1 r) r r
+    l = takeColumns r (l_u |*| tl) |+| diagRect (konst 1 r) r r
     u = l_u |*| tu
     (p,s) = fixPerm r perm
-    l' = (l_u |*| tl) |+| diagRect (constant 1 c) r c
+    l' = (l_u |*| tl) |+| diagRect (konst 1 c) r c
     u' = takeRows c (l_u |*| tu)
     (|+|) = add
     (|*|) = mul