diff --git a/CHANGES b/CHANGES
--- a/CHANGES
+++ b/CHANGES
@@ -1,3 +1,17 @@
+0.8.2.0
+=======
+
+- fromRows/fromColumns now automatically expand vectors of dim 1
+  to match the common dimension.
+  fromBlocks also replicates single row/column matrices.
+  Previously all dimensions had to be exactly the same.
+
+- display utilities: dispf, disps, vecdisp
+
+- scalar
+
+- minimizeV, minimizeVD, using Vector instead of lists.
+
 0.8.1.0
 =======
 
diff --git a/README b/README
--- a/README
+++ b/README
@@ -106,3 +106,5 @@
 - Erik de Castro Lopo added buildVector and buildMatrix, which take a
   size parameter(s) and a function that maps vector/matrix indices
   to the values at that position.
+
+- Jean-Francois Tremblay discovered an error in the tutorial.
diff --git a/examples/lie.hs b/examples/lie.hs
--- a/examples/lie.hs
+++ b/examples/lie.hs
@@ -1,9 +1,8 @@
 -- The magic of Lie Algebra
 
 import Numeric.LinearAlgebra
-import Text.Printf(printf)
 
-disp = putStrLn . format "  " (printf "%.5f")
+disp = putStrLn . dispf 5
 
 rot1 :: Double -> Matrix Double
 rot1 a = (3><3)
@@ -50,14 +49,14 @@
 a & b = a <> b - b <> a
 
 infixl 6 |+|
-a |+| b = a + b + a & b */2  + (a-b) & (a & b) */12
+a |+| b = a + b + a&b /2  + (a-b)&(a & b) /12
 
 main = do
    let a =  45*deg
        b =  50*deg
        c = -30*deg
        exact = rot3 a <> rot1 b <> rot2 c
-       lie = a.*g3 |+| b.*g1 |+| c.*g2
+       lie = scalar a * g3 |+| scalar b * g1 |+| scalar c * g2
    putStrLn "position in the tangent space:"
    disp lie
    putStrLn "exponential map back to the group (2 terms):"
diff --git a/examples/parallel.hs b/examples/parallel.hs
--- a/examples/parallel.hs
+++ b/examples/parallel.hs
@@ -12,10 +12,13 @@
 parMul p x y = fromBlocks [ inParallel ( map (x <>) ys ) ]
     where ys = splitColumns p y
 
+-- x <||> y = fromColumns . inParallel . map (x <>) . toColumns $ y
+
 main = do
     n <- (read . head) `fmap` getArgs
     let m = ident n :: Matrix Double
     time $ print $ vectorMax $ takeDiag $ m <> m
+--    time $ print $ vectorMax $ takeDiag $ m <||> m
     time $ print $ vectorMax $ takeDiag $ parMul 2 m m
     time $ print $ vectorMax $ takeDiag $ parMul 4 m m
     time $ print $ vectorMax $ takeDiag $ parMul 8 m m
diff --git a/hmatrix.cabal b/hmatrix.cabal
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -1,12 +1,12 @@
 Name:               hmatrix
-Version:            0.8.1.1
+Version:            0.8.2.0
 License:            GPL
 License-file:       LICENSE
 Author:             Alberto Ruiz
 Maintainer:         Alberto Ruiz <aruiz@um.es>
 Stability:          provisional
-Homepage:           http://www.hmatrix.googlepages.com
-Synopsis:           Linear algebra and numerical computations
+Homepage:           http://code.haskell.org/hmatrix
+Synopsis:           Linear algebra and numerical computation
 Description:        Purely functional interface to basic linear algebra
                     and other numerical computations, internally implemented using
                     GSL, BLAS and LAPACK.
@@ -64,14 +64,14 @@
 
 library
     if flag(splitBase)
-      build-depends:    base >= 3 && < 5, array < 0.3.0
+      build-depends:    base >= 3 && < 5, array
     else
       build-depends:    base < 3
 
-    Build-Depends:      haskell98 < 1.0.1.1,
+    Build-Depends:      haskell98,
                         QuickCheck, HUnit,
                         storable-complex,
-                        process < 1.0.1.2
+                        process
 
     Extensions:         ForeignFunctionInterface,
                         CPP
diff --git a/lib/Data/Packed/Internal/Common.hs b/lib/Data/Packed/Internal/Common.hs
--- a/lib/Data/Packed/Internal/Common.hs
+++ b/lib/Data/Packed/Internal/Common.hs
@@ -18,7 +18,7 @@
   Adapt,
   app1, app2, app3, app4,
   (//), check,
-  partit, common,
+  partit, common, compatdim,
   fi
 ) where
 
@@ -40,6 +40,12 @@
     commonval [] = Nothing
     commonval [a] = Just a
     commonval (a:b:xs) = if a==b then commonval (b:xs) else Nothing
+
+-- | common value with \"adaptable\" 1
+compatdim :: [Int] -> Maybe Int
+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
 
 -- | postfix function application (@flip ($)@)
 (//) :: x -> (x -> y) -> y
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -135,11 +135,16 @@
 toLists :: (Element t) => Matrix t -> [[t]]
 toLists m = partit (cols m) . toList . flatten $ m
 
--- | creates a Matrix from a list of vectors
+-- | Create a matrix from a list of vectors.
+-- All vectors must have the same dimension,
+-- or dimension 1, which is are automatically expanded.
 fromRows :: Element t => [Vector t] -> Matrix t
-fromRows vs = case common dim vs of
+fromRows vs = case compatdim (map dim vs) of
     Nothing -> error "fromRows applied to [] or to vectors with different sizes"
-    Just c  -> reshape c (join vs)
+    Just c  -> reshape c . join . map (adapt c) $ vs
+  where
+    adapt c v | dim v == c = v
+              | otherwise = constantD (v@>0) c
 
 -- | extracts the rows of a matrix as a list of vectors
 toRows :: Element t => Matrix t -> [Vector t]
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -1,4 +1,3 @@
-{-# OPTIONS_GHC -fglasgow-exts #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Data.Packed.Matrix
@@ -29,7 +28,7 @@
     extractRows,
     ident, diag, diagRect, takeDiag,
     liftMatrix, liftMatrix2,
-    format,
+    format, dispf, disps,
     loadMatrix, saveMatrix, fromFile, fileDimensions,
     readMatrix, fromArray2D
 ) where
@@ -40,33 +39,67 @@
 import Data.List(transpose,intersperse)
 import Data.Array
 import System.Process(readProcess)
+import Text.Printf(printf)
 
 -- | creates a matrix from a vertical list of matrices
 joinVert :: Element t => [Matrix t] -> Matrix t
 joinVert ms = case common cols ms of
-    Nothing -> error "joinVert on matrices with different number of columns"
+    Nothing -> error "(impossible) joinVert on matrices with different number of columns"
     Just c  -> reshape c $ join (map flatten ms)
 
 -- | creates a matrix from a horizontal list of matrices
 joinHoriz :: Element t => [Matrix t] -> Matrix t
 joinHoriz ms = trans. joinVert . map trans $ ms
 
-{- | Creates a matrix from blocks given as a list of lists of matrices:
+{- | Creates a matrix from blocks given as a list of lists of matrices.
 
-@\> let a = 'diag' $ 'fromList' [5,7,2]
-\> let b = 'reshape' 4 $ 'constant' (-1) 12
-\> fromBlocks [[a,b],[b,a]]
-(6><7)
- [  5.0,  0.0,  0.0, -1.0, -1.0, -1.0, -1.0
- ,  0.0,  7.0,  0.0, -1.0, -1.0, -1.0, -1.0
- ,  0.0,  0.0,  2.0, -1.0, -1.0, -1.0, -1.0
- , -1.0, -1.0, -1.0, -1.0,  5.0,  0.0,  0.0
- , -1.0, -1.0, -1.0, -1.0,  0.0,  7.0,  0.0
- , -1.0, -1.0, -1.0, -1.0,  0.0,  0.0,  2.0 ]@
+Single row/column components are automatically expanded to match the
+corresponding common row and column:
+
+@\> let disp = putStr . dispf 2
+\> let vector xs = fromList xs :: Vector Double
+\> let diagl = diag . vector
+\> let rowm = asRow . vector
+
+\> disp $ fromBlocks [[ident 5, 7, rowm[10,20]], [3, diagl[1,2,3], 0]]
+
+8x10
+1  0  0  0  0  7  7  7  10  20
+0  1  0  0  0  7  7  7  10  20
+0  0  1  0  0  7  7  7  10  20
+0  0  0  1  0  7  7  7  10  20
+0  0  0  0  1  7  7  7  10  20
+3  3  3  3  3  1  0  0   0   0
+3  3  3  3  3  0  2  0   0   0
+3  3  3  3  3  0  0  3   0   0@
 -}
 fromBlocks :: Element t => [[Matrix t]] -> Matrix t
-fromBlocks = joinVert . map joinHoriz 
+fromBlocks = fromBlocksRaw . adaptBlocks
 
+fromBlocksRaw mms = joinVert . map joinHoriz $ mms
+
+adaptBlocks ms = ms' where
+    bc = case common length ms of
+          Just c -> c
+          Nothing -> error "fromBlocks requires rectangular [[Matrix]]"
+    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)
+
+    g [Just nr,Just nc] m
+                | nr == r && nc == c = m
+                | r == 1 && c == 1 = reshape nc (constant x (nr*nc))
+                | r == 1 = fromRows (replicate nr (flatten m))
+                | otherwise = fromColumns (replicate nc (flatten m))
+      where
+        r = rows m
+        c = cols m
+        x = m@@>(0,0)
+    g _ _ = error "inconsistent dimensions in fromBlocks"
+
+-----------------------------------------------------------
+
 -- | Reverse rows 
 flipud :: Element t => Matrix t -> Matrix t
 flipud m = fromRows . reverse . toRows $ m
@@ -229,7 +262,7 @@
 
 @import Text.Printf(printf)@
 
-@disp = putStrLn . format \"  \" (printf \"%.2f\")@
+@disp = putStr . format \"  \" (printf \"%.2f\")@
 
 -}
 format :: (Element t) => String -> (t -> String) -> Matrix t -> String
@@ -244,6 +277,71 @@
 dispC :: Int -> Matrix (Complex Double) -> IO ()
 dispC d m = disp m (shfc d)
 -}
+
+-------------------------------------------------------------------
+-- display utilities
+
+{- | Show a matrix with \"autoscaling\" and a given number of decimal places.
+
+@disp = putStr . disps 2
+
+\> disp $ 120 * (3><4) [1..]
+3x4  E3
+ 0.12  0.24  0.36  0.48
+ 0.60  0.72  0.84  0.96
+ 1.08  1.20  1.32  1.44
+@
+-}
+disps :: Int -> Matrix Double -> String
+disps d x = sdims x ++ "  " ++ formatScaled d x
+
+{- | Show a matrix with a given number of decimal places.
+
+@disp = putStr . dispf 3
+
+\> disp (1/3 + ident 4)
+4x4
+1.333  0.333  0.333  0.333
+0.333  1.333  0.333  0.333
+0.333  0.333  1.333  0.333
+0.333  0.333  0.333  1.333
+@
+-}
+dispf :: Int -> Matrix Double -> String
+dispf d x = sdims x ++ "\n" ++ formatFixed (if isInt x then 0 else d) x
+
+sdims x = show (rows x) ++ "x" ++ show (cols x)
+
+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
+
+formatScaled dec t = "E"++show o++"\n" ++ ss
+    where ss = format " " (printf fmt. g) t
+          g x | o >= 0    = x/10^(o::Int)
+              | otherwise = x*10^(-o)
+          o = floor $ maximum $ map (logBase 10 . abs) $ toList $ flatten t
+          fmt = '%':show (dec+3) ++ '.':show dec ++"f"
+
+{- | Show a vector using a function for showing matrices.
+
+@disp = putStr . vecdisp (dispf 2)
+
+\> disp (linspace 10 (0,1))
+10 |> 0.00  0.11  0.22  0.33  0.44  0.56  0.67  0.78  0.89  1.00
+@
+-}
+vecdisp :: (Element t) => (Matrix t -> String) -> Vector t -> String
+vecdisp f v
+    = ((show (dim v) ++ " |> ") ++) . (++"\n")
+    . unwords . lines .  tail . dropWhile (not . (`elem` " \n"))
+    . f . trans . reshape 1
+    $ v
+
+
+--------------------------------------------------------------------
 
 -- | reads a matrix from a string containing a table of numbers.
 readMatrix :: String -> Matrix Double
diff --git a/lib/Data/Packed/Random.hs b/lib/Data/Packed/Random.hs
--- a/lib/Data/Packed/Random.hs
+++ b/lib/Data/Packed/Random.hs
@@ -21,6 +21,7 @@
 
 import Numeric.GSL.Vector
 import Data.Packed
+import Numeric.LinearAlgebra.Linear
 import Numeric.LinearAlgebra.Algorithms
 import Numeric.LinearAlgebra.Instances()
 import Numeric.LinearAlgebra.Interface
@@ -54,7 +55,7 @@
     d = dim a
     dat = toRows $ reshape n $ randomVector seed Uniform (n*d)
     am = constant 1 n `outer` a
-    m = fromColumns (zipWith (.*) cs dat) + am
+    m = fromColumns (zipWith scale cs dat) + am
 
 ------------ utilities -------------------------------
 
diff --git a/lib/Data/Packed/Vector.hs b/lib/Data/Packed/Vector.hs
--- a/lib/Data/Packed/Vector.hs
+++ b/lib/Data/Packed/Vector.hs
@@ -18,6 +18,7 @@
     dim, (@>),
     subVector, join,
     constant, linspace,
+    vecdisp,
     vectorMax, vectorMin, vectorMaxIndex, vectorMinIndex,
     mapVector, zipVector,
     fscanfVector, fprintfVector, freadVector, fwriteVector,
@@ -32,6 +33,10 @@
 
 @\> linspace 5 (-3,7)
 5 |> [-3.0,-0.5,2.0,4.5,7.0]@
+
+Logarithmic spacing can be defined as follows:
+
+@logspace n (a,b) = 10 ** linspace n (a,b)@
 -}
 linspace :: Int -> (Double, Double) -> Vector Double
 linspace n (a,b) = add a $ scale s  $ fromList [0 .. fromIntegral n-1]
@@ -70,3 +75,19 @@
 buildVector :: Element a => Int -> (Int -> a) -> Vector a
 buildVector len f =
     fromList $ map f [0 .. (len - 1)]
+
+
+{- | Show a vector using a function for showing matrices.
+
+@disp = putStr . vecdisp ('dispf' 2)
+
+\> disp ('linspace' 10 (0,1))
+10 |> 0.00  0.11  0.22  0.33  0.44  0.56  0.67  0.78  0.89  1.00
+@
+-}
+vecdisp :: (Element t) => (Matrix t -> String) -> Vector t -> String
+vecdisp f v
+    = ((show (dim v) ++ " |> ") ++) . (++"\n")
+    . unwords . lines .  tail . dropWhile (not . (`elem` " \n"))
+    . f . trans . reshape 1
+    $ v
diff --git a/lib/Numeric/GSL/Minimization.hs b/lib/Numeric/GSL/Minimization.hs
--- a/lib/Numeric/GSL/Minimization.hs
+++ b/lib/Numeric/GSL/Minimization.hs
@@ -51,8 +51,9 @@
 
 -----------------------------------------------------------------------------
 module Numeric.GSL.Minimization (
-    minimize, MinimizeMethod(..),
-    minimizeD, MinimizeMethodD(..),
+    minimize, minimizeV, MinimizeMethod(..),
+    minimizeD, minimizeVD, MinimizeMethodD(..),
+
     minimizeNMSimplex,
     minimizeConjugateGradient,
     minimizeVectorBFGS2
@@ -82,7 +83,7 @@
                     | NMSimplex2
                     deriving (Enum,Eq,Show,Bounded)
 
--- | Minimization without derivatives.
+-- | Minimization without derivatives
 minimize :: MinimizeMethod
          -> Double              -- ^ desired precision of the solution (size test)
          -> Int                 -- ^ maximum number of iterations allowed
@@ -91,8 +92,40 @@
          -> [Double]            -- ^ starting point
          -> ([Double], Matrix Double) -- ^ solution vector and optimization path
 
-minimize method = minimizeGen (fi (fromEnum method))
+-- | Minimization without derivatives (vector version)
+minimizeV :: MinimizeMethod
+         -> Double              -- ^ desired precision of the solution (size test)
+         -> Int                 -- ^ maximum number of iterations allowed
+         -> Vector Double       -- ^ sizes of the initial search box
+         -> (Vector Double -> Double) -- ^ function to minimize
+         -> Vector Double            -- ^ starting point
+         -> (Vector Double, Matrix Double) -- ^ solution vector and optimization path
 
+minimize method eps maxit sz f xi = v2l $ minimizeV method eps maxit (fromList sz) (f.toList) (fromList xi)
+    where v2l (v,m) = (toList v, m)
+
+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' ->
+                   createMIO maxit (n+3)
+                         (c_minimize (fi (fromEnum method)) fp eps (fi maxit) // xiv' // szv')
+                         "minimize"
+    let it = round (rawpath @@> (maxit-1,0))
+        path = takeRows it rawpath
+        sol = cdat $ dropColumns 3 $ dropRows (it-1) path
+    freeHaskellFunPtr fp
+    return (sol, path)
+
+
+foreign import ccall "gsl-aux.h minimize"
+    c_minimize:: CInt -> FunPtr (CInt -> Ptr Double -> Double) -> Double -> CInt -> TVVM
+
+----------------------------------------------------------------------------------
+
+
 data MinimizeMethodD = ConjugateFR
                      | ConjugatePR
                      | VectorBFGS
@@ -111,49 +144,35 @@
     -> [Double]               -- ^ starting point
     -> ([Double], Matrix Double) -- ^ solution vector and optimization path
 
-minimizeD method = minimizeDGen (fi (fromEnum method))
-
--------------------------------------------------------------------------
-
-ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2
-
-minimizeGen method eps maxit sz f xi = unsafePerformIO $ do
-    let xiv = fromList xi
-        szv = fromList sz
-        n   = dim xiv
-    fp <- mkVecfun (iv (f.toList))
-    rawpath <- ww2 withVector xiv withVector szv $ \xiv' szv' ->
-                   createMIO maxit (n+3)
-                         (c_minimize method fp eps (fi maxit) // xiv' // szv')
-                         "minimize"
-    let it = round (rawpath @@> (maxit-1,0))
-        path = takeRows it rawpath
-        [sol] = toLists $ dropRows (it-1) path
-    freeHaskellFunPtr fp
-    return (drop 3 sol, path)
-
-
-foreign import ccall "gsl-aux.h minimize"
-    c_minimize:: CInt -> FunPtr (CInt -> Ptr Double -> Double) -> Double -> CInt -> TVVM
-
-----------------------------------------------------------------------------------
+-- | Minimization with derivatives (vector version)
+minimizeVD :: MinimizeMethodD
+    -> Double                 -- ^ desired precision of the solution (gradient test)
+    -> Int                    -- ^ maximum number of iterations allowed
+    -> Double                 -- ^ size of the first trial step
+    -> Double                 -- ^ tol (precise meaning depends on method)
+    -> (Vector Double -> Double)   -- ^ function to minimize
+    -> (Vector Double -> Vector Double) -- ^ gradient
+    -> Vector Double               -- ^ starting point
+    -> (Vector Double, Matrix Double) -- ^ solution vector and optimization path
 
+minimizeD method eps maxit istep tol f df xi = v2l $ minimizeVD
+          method eps maxit istep tol (f.toList) (fromList.df.toList) (fromList xi)
+    where v2l (v,m) = (toList v, m)
 
 
-minimizeDGen method eps maxit istep tol f df xi = unsafePerformIO $ do
-    let xiv = fromList xi
-        n = dim xiv
-        f' = f . toList
-        df' = (checkdim1 n .fromList . df . toList)
+minimizeVD method eps maxit istep tol f df xiv = unsafePerformIO $ do
+    let n = dim xiv
+        f' = f
+        df' = (checkdim1 n . df)
     fp <- mkVecfun (iv f')
     dfp <- mkVecVecfun (aux_vTov df')
     rawpath <- withVector xiv $ \xiv' ->
                     createMIO maxit (n+2)
-                         (c_minimizeD method fp dfp istep tol eps (fi maxit) // xiv')
+                         (c_minimizeD (fi (fromEnum method)) fp dfp istep tol eps (fi maxit) // xiv')
                          "minimizeD"
     let it = round (rawpath @@> (maxit-1,0))
         path = takeRows it rawpath
-        sol = toList $ cdat $ dropColumns 2 $ dropRows (it-1) path
+        sol = cdat $ dropColumns 2 $ dropRows (it-1) path
     freeHaskellFunPtr fp
     freeHaskellFunPtr dfp
     return (sol,path)
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -147,6 +147,14 @@
 
 --------------------------------------------------------------
 
+square m = rows m == cols m
+
+vertical m = rows m >= cols m
+
+exactHermitian m = m `equal` ctrans m
+
+--------------------------------------------------------------
+
 -- | Full singular value decomposition.
 svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
 svd = svd'
@@ -179,7 +187,6 @@
     u' = takeColumns d u
     v' = takeColumns d v
 
-vertical m = rows m >= cols m
 
 -- | Singular values and all right singular vectors.
 rightSV :: Field t => Matrix t -> (Vector Double, Matrix t)
@@ -255,12 +262,12 @@
 --
 -- If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@
 eigSH :: Field t => Matrix t -> (Vector Double, Matrix t)
-eigSH m | m `equal` ctrans m = eigSH' m
+eigSH m | exactHermitian m = eigSH' m
         | otherwise = error "eigSH requires complex hermitian or real symmetric matrix"
 
 -- | Eigenvalues of a complex hermitian or real symmetric matrix.
 eigenvaluesSH :: Field t => Matrix t -> Vector Double
-eigenvaluesSH m | m `equal` ctrans m = eigenvaluesSH' m
+eigenvaluesSH m | exactHermitian m = eigenvaluesSH' m
                 | otherwise = error "eigenvaluesSH requires complex hermitian or real symmetric matrix"
 
 --------------------------------------------------------------
@@ -317,14 +324,12 @@
 --
 -- If @c = chol m@ then @m == ctrans c \<> c@.
 chol :: Field t => Matrix t ->  Matrix t
-chol m | m `equal` ctrans m = cholSH m
+chol m | exactHermitian m = cholSH m
        | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix"
 
 
 
 
-square m = rows m == cols m
-
 -- | Determinant of a square matrix.
 det :: Field t => Matrix t -> t
 det m | square m = s * (product $ toList $ takeDiag $ lup)
@@ -569,7 +574,7 @@
                     then Just (l,v)
                     else Nothing
     where n = rows m
-          (l,v) = if m `equal` ctrans m
+          (l,v) = if exactHermitian m
                     then let (l',v') = eigSH m in (real l', v')
                     else eig m
 
diff --git a/lib/Numeric/LinearAlgebra/Instances.hs b/lib/Numeric/LinearAlgebra/Instances.hs
--- a/lib/Numeric/LinearAlgebra/Instances.hs
+++ b/lib/Numeric/LinearAlgebra/Instances.hs
@@ -80,8 +80,8 @@
              || rows m2 == 1 && cols m2 == 1
              || rows m1 == rows m2 && cols m1 == cols m2
 
-instance (Eq a, Element a) => Eq (Vector a) where
-    a == b = dim a == dim b && toList a == toList b
+instance Linear Vector a => Eq (Vector a) where
+    (==) = equal
 
 instance Num (Vector Double) where
     (+) = adaptScalar addConstant add (flip addConstant)
@@ -99,8 +99,8 @@
     abs = vectorMapC Abs
     fromInteger = fromList . return . fromInteger
 
-instance (Eq a, Element a) => Eq (Matrix a) where
-    a == b = cols a == cols b && flatten a == flatten b
+instance Linear Matrix a => Eq (Matrix a) where
+    (==) = equal
 
 instance (Linear Matrix a, Num (Vector a)) => Num (Matrix a) where
     (+) = liftMatrix2' (+)
diff --git a/lib/Numeric/LinearAlgebra/Interface.hs b/lib/Numeric/LinearAlgebra/Interface.hs
--- a/lib/Numeric/LinearAlgebra/Interface.hs
+++ b/lib/Numeric/LinearAlgebra/Interface.hs
@@ -49,12 +49,16 @@
 
 ----------------------------------------------------
 
+{-# DEPRECATED (.*) "use scale a x or scalar a * x" #-}
+
 -- | @x .* a = scale x a@
 (.*) :: (Linear c a) => a -> c a -> c a
 infixl 7 .*
 a .* x = scale a x
 
 ----------------------------------------------------
+
+{-# DEPRECATED (*/) "use scale (recip a) x or x / scalar a" #-}
 
 -- | @a *\/ x = scale (recip x) a@
 (*/) :: (Linear c a) => c a -> a -> c a
diff --git a/lib/Numeric/LinearAlgebra/Linear.hs b/lib/Numeric/LinearAlgebra/Linear.hs
--- a/lib/Numeric/LinearAlgebra/Linear.hs
+++ b/lib/Numeric/LinearAlgebra/Linear.hs
@@ -23,7 +23,13 @@
 
 -- | A generic interface for vectors and matrices to a few element-by-element functions in Numeric.GSL.Vector.
 class (Container c e) => Linear c e where
+    -- | create a structure with a single element
+    scalar      :: e -> c e
     scale       :: e -> c e -> c e
+    -- | scale the element by element reciprocal of the object:
+    --
+    -- @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@
+    scaleRecip  :: e -> c e -> c e
     addConstant :: e -> c e -> c e
     add         :: c e -> c e -> c e
     sub         :: c e -> c e -> c e
@@ -31,11 +37,9 @@
     mul         :: c e -> c e -> c e
     -- | element by element division
     divide      :: c e -> c e -> c e
-    -- | scale the element by element reciprocal of the object: @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@
-    scaleRecip  :: e -> c e -> c e
     equal       :: c e -> c e -> Bool
---  numequal    :: Double -> c e -> c e -> Bool
 
+
 instance Linear Vector Double where
     scale = vectorMapValR Scale
     scaleRecip = vectorMapValR Recip
@@ -45,6 +49,7 @@
     mul = vectorZipR Mul
     divide = vectorZipR Div
     equal u v = dim u == dim v && vectorMax (vectorMapR Abs (sub u v)) == 0.0
+    scalar x = fromList [x]
 
 instance Linear Vector (Complex Double) where
     scale = vectorMapValC Scale
@@ -55,6 +60,7 @@
     mul = vectorZipC Mul
     divide = vectorZipC Div
     equal u v = dim u == dim v && vectorMax (mapVector magnitude (sub u v)) == 0.0
+    scalar x = fromList [x]
 
 instance (Linear Vector a, Container Matrix a) => (Linear Matrix a) where
     scale x = liftMatrix (scale x)
@@ -65,3 +71,4 @@
     mul = liftMatrix2 mul
     divide = liftMatrix2 divide
     equal a b = cols a == cols b && flatten a `equal` flatten b
+    scalar x = (1><1) [x]
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs
--- a/lib/Numeric/LinearAlgebra/Tests.hs
+++ b/lib/Numeric/LinearAlgebra/Tests.hs
@@ -262,8 +262,8 @@
         , utest "expm1" (expmTest1)
         , utest "expm2" (expmTest2)
         , utest "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM)
-        , utest "arith2" $ (((1+i) .* ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( (140*i-51).*1 :: CM)
-        , utest "arith3" $ exp (i.*ones(10,10)*pi) + 1 |~| 0
+        , utest "arith2" $ ((scalar (1+i) * ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( scalar (140*i-51) :: CM)
+        , utest "arith3" $ exp (scalar i * ones(10,10)*pi) + 1 |~| 0
         , utest "<\\>"   $ (3><2) [2,0,0,3,1,1::Double] <\> 3|>[4,9,5] |~| 2|>[2,3]
         , utest "gamma" (gamma 5 == 24.0)
         , besselTest
@@ -279,10 +279,11 @@
                      && ident 5 == buildMatrix 5 5 (\(r,c) -> if r==c then 1::Double else 0)
         , 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)
         ]
     return ()
 
-makeUnitary v | realPart n > 1    = v */ n
+makeUnitary v | realPart n > 1    = v / scalar n
               | otherwise = v
     where n = sqrt (conj v <.> v)
 
@@ -321,6 +322,8 @@
         where c = cos a
               s = sin a
 
+multb n = foldl1' (<>) (replicate (10^6) (ident n :: Matrix Double))
+
 --------------------------------
 
 multBench = do
@@ -328,6 +331,15 @@
     let b = ident 2000 :: Matrix Double
     a `seq` b `seq` putStrLn ""
     time "product of 1M different 3x3 matrices" (manymult (10^6))
+    putStrLn ""
+    time "product of 1M constant  1x1 matrices" (multb 1)
+    time "product of 1M constant  3x3 matrices" (multb 3)
+    --time "product of 1M constant  5x5 matrices" (multb 5)
+    time "product of 1M const.  10x10 matrices" (multb 10)
+    --time "product of 1M const.  15x15 matrices" (multb 15)
+    time "product of 1M const.  20x20 matrices" (multb 20)
+    --time "product of 1M const.  25x25 matrices" (multb 25)
+    putStrLn ""
     time "product (1000 x 1000)<>(1000 x 1000)" (a<>a)
     time "product (2000 x 2000)<>(2000 x 2000)" (b<>b)
 
diff --git a/lib/Numeric/LinearAlgebra/Tests/Instances.hs b/lib/Numeric/LinearAlgebra/Tests/Instances.hs
--- a/lib/Numeric/LinearAlgebra/Tests/Instances.hs
+++ b/lib/Numeric/LinearAlgebra/Tests/Instances.hs
@@ -180,7 +180,7 @@
         l <- replicateM n (choose (0,100))
         let s = diag (fromList l)
             p = v <> real s <> ctrans v
-        return $ PosDef (0.5 .* p + 0.5 .* ctrans p)
+        return $ PosDef (0.5 * p + 0.5 * ctrans p)
 
 #if MIN_VERSION_QuickCheck(2,0,0)
 #else