diff --git a/INSTALL b/INSTALL
new file mode 100644
--- /dev/null
+++ b/INSTALL
@@ -0,0 +1,3 @@
+INSTALLATION
+
+$ cabal install hTensor
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,2 @@
+Copyright Alberto Ruiz 2009
+GPL license
diff --git a/README b/README
new file mode 100644
--- /dev/null
+++ b/README
diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,4 @@
+#! /usr/bin/env runhaskell
+
+> import Distribution.Simple
+> main = defaultMain
diff --git a/examples/array.hs b/examples/array.hs
new file mode 100644
--- /dev/null
+++ b/examples/array.hs
@@ -0,0 +1,59 @@
+import Numeric.LinearAlgebra.Array
+import Numeric.LinearAlgebra.Array.Util
+import Control.Applicative
+
+-- 'listArray' specialized for Array Double
+infixl 9 #
+(#) :: [Int] -> [Double] -> Array Double
+(#) = listArray
+
+(<|) :: Name -> [Array Double] -> Array Double
+infixl 8 <|
+n <| ls = index n ls
+
+i = ("i" <|)
+j = ("j" <|)
+k = ("k" <|)
+
+sh x = putStrLn . formatFixed 2 $ x
+
+a = [3,4,2] # [1..]
+b = [2,3] # [5,6]
+c = 7 :: Array Double
+s = [2,2,2,2] # [1..]
+
+t = [3,3,3,3]#[1 ..]!"ijkl"
+
+q = [2,4,3] # (fun <$> r 2 <*> r 4 <*> r 3) !"ijk"
+    where r k = [1..k]
+          fun = \i j k -> i*2*j-k
+
+mk ds f = ds # map f (sequence $ map (enumFromTo 1 . fromIntegral) $ ds)
+
+m = j [i[2,0,0], i[1,0,1], i[0,3,0]] ~> "ij"
+
+main = do
+    putStrLn "8-dimensional array"
+    sh $ (replicate 8 2)#[1::Double ..]!(take 8 ['a'..])
+    ------------------------
+    putStrLn "different display formats"
+    sh $ a!"ijk"
+    printA "%7.3f" a
+    putStrLn . formatScaled 2 $ a!"ijk"
+    sh . noIdx $ a
+    ------------------------
+    putStrLn "array defined using a function"
+    sh q
+    sh $ mk [2,4,3] (\[i,j,k] -> i*2*j-k) !"ijk"
+    ------------------------
+    putStrLn "contraction"
+    sh t
+    sh $ t!"ijkk"
+    ------------------------
+    putStrLn "tensor product"
+    sh $ m
+    sh $ (t !"pqrs" * m!"kr") ~> "pqks"
+    ------------------------
+    putStrLn "automatic conformability"
+    sh $ j[1,2,3] + k[10,20]
+    sh $ k [m, 3*m-1, 7, i[1,2,3]]
diff --git a/examples/exterior.hs b/examples/exterior.hs
new file mode 100644
--- /dev/null
+++ b/examples/exterior.hs
@@ -0,0 +1,41 @@
+import Numeric.LinearAlgebra.Exterior
+import Numeric.LinearAlgebra.Array.Util(formatFixed,asMatrix)
+import Numeric.LinearAlgebra (det,(><))
+
+printAS = print . asMultivector
+
+sh = putStrLn . formatFixed 2
+
+-- 'listTensor' specialized for Tensor Double
+infixl 9 #
+(#) :: [Int] -> [Double] -> Tensor Double
+(#) = listTensor
+
+m = [3,-3]#[ 1,2,5,
+             1,2,8,
+            -2,0,4]
+
+eps = leviCivita 3
+
+a = vector [1,0,0] /\ vector [0,1,0]
+b = vector [2,0,100,0] /\ vector [0,3,0,0] /\ vector [0,0,4,0]
+
+im = eps!"ijb"* m!"pi" * m!"qj" * cov eps!"apq"
+
+main = do
+    putStrLn "exterior product"
+    print a
+    sh a
+    printAS a
+    printAS b
+    putStrLn "\ndeterminant"
+    printAS $ cov eps!"pqr" * m!"pi" * m!"qj" * m!"rk"
+    print $ det (asMatrix m)
+    putStrLn "\ninverse"
+    sh $ im
+    sh $ im!"ik" * m!"kj"
+    putStrLn "\nmeet and join"
+    printAS $ (vector [1,0,1] /\ vector [0,1,0]) \/ (vector [1,1,0] /\ vector [0,0,1])
+    putStrLn "\nEuclidean inner product of r-vectors"
+    printAS $ (vector [1,0,1] /\ vector [0,1,0]) `inner` (vector[1,1,1])
+    print $ vector[3,5] `inner` vector[2,1]
diff --git a/examples/geom.hs b/examples/geom.hs
new file mode 100644
--- /dev/null
+++ b/examples/geom.hs
@@ -0,0 +1,29 @@
+import Numeric.LinearAlgebra.Multivector
+
+o  = e 4
+x = vector[1,0,0,1]
+y = o + e 2
+z = vector[0,0,1] + e 4
+
+p = x /\ y /\ z
+
+p1 = o /\ x /\ y
+p2 = z /\ vector[1,0,1,1] /\ vector[0,1,1,1]
+
+l = o /\ vector[1,1,1,1]
+
+rot = rotor 3 (pi/4) (e 3)
+
+l' = rot * l * rever rot
+
+inh v = v / (v -| e 4) - e 4
+
+main = do
+     print l
+     print $ l \/ p1
+     print $ l \/ p2
+     print $ l \/ p
+     print $ inh $ l \/ p1
+     print $ inh $ l \/ p2
+     print $ inh $ l' \/ p1
+     print $ inh $ l' \/ p2
diff --git a/hTensor.cabal b/hTensor.cabal
new file mode 100644
--- /dev/null
+++ b/hTensor.cabal
@@ -0,0 +1,65 @@
+Name:               hTensor
+Version:            0.1.0
+License:            GPL
+License-file:       LICENSE
+Author:             Alberto Ruiz
+Maintainer:         Alberto Ruiz <aruiz@um.es>
+Stability:          experimental
+Homepage:           http://perception.inf.um.es/tensor
+Synopsis:           Multidimensional arrays and simple tensor computations.
+Description:
+ This is an experimental library for multidimensional arrays,
+ oriented to support simple tensor computations and multilinear
+ algebra.
+ .
+ Array dimensions have an \"identity\" which is preserved
+ in data manipulation. Indices are explicitly selected by name in
+ expressions, and Einstein's summation convention for repeated indices
+ is automatically applied.
+ .
+ The library has a purely functional interface: arrays are immutable,
+ and operations typically work on whole structures which can be assembled
+ and decomposed using simple primitives. Arguments are automatically made conformant
+ by replicating them along extra dimensions appearing in an operation.
+ There is preliminary support for Geometric Algebra.
+ .
+ A tutorial can be found in the website of the project.
+
+Category:           Math
+tested-with:        GHC ==6.10.3
+
+cabal-version:      >=1.2
+build-type:         Simple
+
+extra-source-files: README INSTALL
+
+extra-source-files: examples/array.hs
+                    examples/exterior.hs
+                    examples/geom.hs
+
+flag splitBase
+    description:    Choose the new smaller, split-up base package.
+
+library
+    if flag(splitBase)
+      build-depends:    base >= 3 && < 5
+    else
+      build-depends:    base < 3
+
+    Build-Depends:      haskell98, hmatrix >= 0.5, containers
+
+    hs-source-dirs:     lib
+    Exposed-modules:    Numeric.LinearAlgebra.Array.Simple
+                        Numeric.LinearAlgebra.Array.Util
+                        Numeric.LinearAlgebra.Array
+                        Numeric.LinearAlgebra.Tensor
+                        Numeric.LinearAlgebra.Exterior
+                        Numeric.LinearAlgebra.Multivector
+
+    other-modules:      Numeric.LinearAlgebra.Array.Internal
+
+    ghc-prof-options:   -auto-all
+
+    ghc-options:        -Wall
+                        -fno-warn-missing-signatures
+                        -fno-warn-orphans
diff --git a/lib/Numeric/LinearAlgebra/Array.hs b/lib/Numeric/LinearAlgebra/Array.hs
new file mode 100644
--- /dev/null
+++ b/lib/Numeric/LinearAlgebra/Array.hs
@@ -0,0 +1,77 @@
+{-# LANGUAGE UndecidableInstances #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.LinearAlgebra.Array
+-- Copyright   :  (c) Alberto Ruiz 2009
+-- License     :  GPL
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-- Simple multidimensional array with useful numeric instances.
+--
+-- Contractions only require equal dimension.
+--
+-----------------------------------------------------------------------------
+
+module Numeric.LinearAlgebra.Array (
+    None(..),
+    Array,
+    listArray,
+    scalar,
+    index,
+    (!),(~>),
+    (.*),
+    printA
+) where
+
+import Numeric.LinearAlgebra.Array.Simple
+import Numeric.LinearAlgebra.Array.Util
+import Numeric.LinearAlgebra.Array.Internal(printA)
+import Data.Packed(Vector)
+
+-- | Create an 'Array' from a list of parts (@index = 'newIndex' 'None'@).
+index :: Coord t => Name -> [Array t] -> Array t
+index = newIndex None
+
+
+-- | Element by element product.
+infixl 7 .*
+(.*) :: (Coord a, Compat i) => NArray i a -> NArray i a -> NArray i a
+(.*) = zipArray (*)
+
+instance (Coord t, Compat i) => Eq (NArray i t) where
+    t1 == t2 = sameStructure t1 t2 && coords t1 == coords (reorder (names t1) t2)
+
+instance (Show (NArray i t), Coord t, Compat i) => Num (NArray i t) where
+    (+) = zipArray (+)
+    (*) = (|*|)
+    negate t = scalar (-1) * t
+    fromInteger n = scalar (fromInteger n)
+    abs _ = error "abs for arrays not defined"
+    signum _ = error "signum for arrays not defined"
+
+instance (Coord t, Compat i, Num (NArray i t)) => Fractional (NArray i t) where
+    fromRational = scalar . fromRational
+    (/) = zipArray (/)
+    recip = mapArray recip
+
+instance (Coord t, Compat i, Fractional (NArray i t), Floating (Vector t)) => Floating (NArray i t) where
+    sin   = mapArray sin
+    cos   = mapArray cos
+    tan   = mapArray tan
+    asin  = mapArray asin
+    acos  = mapArray acos
+    atan  = mapArray atan
+    sinh  = mapArray sinh
+    cosh  = mapArray cosh
+    tanh  = mapArray tanh
+    asinh = mapArray asinh
+    acosh = mapArray acosh
+    atanh = mapArray atanh
+    exp   = mapArray exp
+    log   = mapArray log
+    (**)  = zipArray (**)
+    sqrt  = mapArray sqrt
+    pi    = scalar pi
diff --git a/lib/Numeric/LinearAlgebra/Array/Internal.hs b/lib/Numeric/LinearAlgebra/Array/Internal.hs
new file mode 100644
--- /dev/null
+++ b/lib/Numeric/LinearAlgebra/Array/Internal.hs
@@ -0,0 +1,565 @@
+{-# OPTIONS_HADDOCK hide #-}
+{-# LANGUAGE FlexibleInstances, FlexibleContexts, MultiParamTypeClasses #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.Array.Internal
+-- Copyright   :  (c) Alberto Ruiz 2009
+-- License     :  GPL
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-- Multidimensional arrays.
+--
+-- The arrays provided by this library are immutable, built on top of hmatrix
+-- structures.
+-- Operations work on complete structures (indexless), and dimensions have \"names\", 
+-- in order to select the desired contractions in tensor computations.
+--
+-- This module contains auxiliary functions not required by the end user.
+
+-----------------------------------------------------------------------------
+
+module Numeric.LinearAlgebra.Array.Internal (
+    -- * Data structures
+    NArray, Idx(..), Name,
+    rank, names, size, typeOf , dims, coords,
+    Compat(..),
+    -- * Array creation
+    scalar,
+    mkNArray,
+    fromVector, fromMatrix,
+    -- * Array manipulation
+    rename,(!),
+    parts,
+    (|*|),
+    zipArray,
+    mapArray,
+    extract,
+    onIndex,
+    -- * Utilities
+    reorder, (~>),
+    sameStructure,
+    conformable,
+    makeConformant,
+    mapTypes,
+    renameRaw,
+    formatArray, formatFixed, formatScaled, printA,
+    showBases,
+    newIndex,
+    dummyAt, noIdx,
+    basisOf,
+    common,
+    Coord,
+    asMatrix, asVector, asScalar
+) where
+
+import Data.Packed
+import Data.List
+import Numeric.LinearAlgebra(outer,multiply,Field)
+import Control.Applicative
+import Data.Function(on)
+import Text.Printf
+
+-- | Types that can be elements of the multidimensional arrays.
+class (Num (Vector t), Field t) => Coord t
+instance Coord Double
+instance Coord (Complex Double)
+
+-- import Debug.Trace
+-- 
+-- debug s f x = trace (s ++ ": " ++ show (f x)) x
+
+-- | indices are denoted by strings, (frequently single-letter)
+type Name = String
+
+-- | Dimension descriptor.
+data Idx i = Idx { iDim  :: Int
+                 , iName :: Name
+                 , iType :: i
+                 } deriving (Eq)
+
+
+
+-- | A multidimensional array with index type i and elements t.
+data NArray i t = A { dims   :: [Idx i]   -- ^ Get detailed dimension information about the array.
+                    , coords :: Vector t  -- ^ Get the coordinates of an array as a
+                                          -- flattened structure (in the order specified by 'dims').
+                    }
+
+
+-- | development function not intended for the end user
+mkNArray :: [Idx i] -> Vector a -> NArray i a
+mkNArray [] _ = error "array with empty dimensions, use scalar"
+mkNArray dms vec = A dms v where
+    ds = map iDim dms
+    n = product ds
+    v = if dim vec == n && minimum ds > 0
+            then vec
+            else error $ show ds ++ " dimensions and " ++
+                         show (dim vec) ++ " coordinates for mkNArray"
+
+-- | Create a 0-dimensional structure.
+scalar :: Coord t => t -> NArray i t
+scalar x = A [] (fromList [x])
+
+
+-- | 'rename' the indices with single-letter names. Equal indices of compatible type are contracted out.
+infixl 8 !
+(!) :: (Coord t, Compat i)
+       => NArray i t
+       -> String   -- ^ new indices
+       -> NArray i t
+t ! ns = rename t (map return ns)
+
+-- | Rename indices. Equal indices are contracted out.
+rename :: (Coord t, Compat i)
+       => NArray i t
+       -> [Name]     -- ^ new names
+       -> NArray i t
+rename t ns = reorder orig (contract t')
+    where t' = renameRaw t ns
+          orig = nub (names t') \\ common1 t'
+
+
+renameRaw (A d v) l | length l == length d = A d' v
+                    | otherwise = error $ "rename " ++ show d ++ " with " ++ show l
+    where d' = zipWith f d l
+          f i n = i {iName=n}
+
+mapDims f (A d v) = A (map f d) v
+
+mapTypes :: (i1 -> i) -> NArray i1 t -> NArray i t
+mapTypes f = mapDims (\i -> i {iType = f (iType i)})
+
+-- mapNames f = mapDims (\i -> i {iName = f (iName i)})
+
+-- | Index names.
+names :: NArray i t -> [Name]
+names = map iName . dims
+
+-- | Dimension of given index.
+size :: Name -> NArray i t -> Int
+size n t = (iDim . head) (filter ((n==).iName) (dims t))
+
+-- | Type of given index.
+typeOf :: Compat i => Name -> NArray i t -> i
+typeOf n t = (iType . head) (filter ((n==).iName) (dims t))
+
+
+-- | The number of dimensions of a multidimensional array.
+rank :: NArray i t -> Int
+rank = length . dims
+
+----------------------------------------------------------
+
+lastIdx name t = ((d1,d2),m) where
+    (d1,d2) = span (\d -> iName d /= name) (dims t)
+    c = product (map iDim d2)
+    m = reshape c (coords t)
+
+firstIdx name t = (nd,m')
+    where ((d1,d2),m) = lastIdx name t
+          m' = reshape c $ flatten $ trans m
+          nd = d2++d1
+          c = dim (coords t) `div` (iDim $ head d2)
+
+
+-- | Create a list of the substructures at the given level.
+parts :: (Coord t) 
+      => NArray i t
+      -> Name        -- ^ index to expand
+      -> [NArray i t]
+parts a name | name `elem` (names a) = map (reorder orig) (partsRaw a name)
+             | otherwise = error $ "parts: " ++ show name ++ " is not a dimension of "++(show $ names a)
+    where orig = names a \\ [name]
+
+partsRaw a name = map f (toRows m)
+    where (_:ds,m) = firstIdx name a
+          f t = A {dims=ds, coords=t}
+
+tridx [] t = t
+tridx (name:rest) t = A (d:ds) (join ts) where
+    d = case lastIdx name t of
+            ((_,d':_),_) -> d'
+            _ -> error "wrong index sequence to reorder"
+    ps = map (tridx rest) (partsRaw t name)
+    ts = map coords ps
+    ds = dims (head ps)
+
+-- | Change the internal layout of coordinates.
+-- The array, considered as an abstract object, does not change.
+reorder :: (Coord t) => [Name] -> NArray i t -> NArray i t
+reorder ns b | ns == names b = b
+             | sort ns == sort (names b) = tridx ns b
+             | otherwise = error $ "wrong index sequence " ++ show ns
+                                    ++ " to reorder "++(show $ names b)
+
+
+-- | 'reorder' (transpose) the dimensions of the array (with single letter names).
+--
+-- Operations are defined by named indices, so the transposed array is operationally equivalent to the original one.
+infixl 8 ~>
+(~>) :: (Coord t) => NArray i t -> String -> NArray i t
+t ~> ns = reorder (map return ns) t
+
+-------------------------------------------------------------
+
+rawProduct (A d1 v1) (A d2 v2) = A (d1++d2) (flatten (outer v1 v2))
+
+----------------------------------------------------------------------
+
+-- | Apply a function (defined on hmatrix 'Vector's) to all elements of a structure.
+-- Use @mapArray (mapVector f)@ for general functions.
+mapArray :: Coord b => (Vector a -> Vector b) -> NArray i a -> NArray i b
+mapArray f t
+    | null (dims t) = scalar (f (coords t)@>0)
+    | otherwise = mkNArray (dims t) (f (coords t))
+
+liftNA2 f (A d1 v1) (A _d2 v2) = A d1 (f v1 v2)
+
+-- | Class of compatible indices for contractions.
+class (Eq a, Show (Idx a)) => Compat a where
+    compat :: Idx a -> Idx a -> Bool
+
+
+
+contract1 t name1 name2 | ok = foldl1' (liftNA2 (+)) y
+                        | otherwise = error $ "wrong contraction1: "
+                                    ++(show $ dims t)++" "
+                                    ++ name1++" "++name2
+    where ok = (compat <$> getName t name1 <*> getName t name2) == Just True
+          x = map (flip partsRaw name2) (partsRaw t name1)
+          y = map head $ zipWith drop [0..] x
+
+getName t name = d where
+    l = filter ((==name).iName) (dims t)
+    d = if null l
+            then Nothing
+            else Just (head l)
+
+contract1c t n = contract1 renamed n n'
+    where n' = " "++n++" " -- forbid spaces in names...
+          renamed = renameRaw (t) auxnames
+          auxnames = h ++ (n':r)
+          (h,_:r) = break (==n) (names t)
+
+common1 t = [ n1 | (a,n1) <- x , (b,n2) <- x, a>b, n1==n2]
+    where x = zip [0 ::Int ..] (names t)
+
+contract t = foldl' contract1c t (common1 t)
+
+----------------------------------------------------------------------
+
+contract2 t1 t2 n | ok = A (tail ds1 ++ tail ds2) (flatten m)
+                  | otherwise = error $ "wrong contraction2: "++ n ++ " of "++
+                                      (show $ dims t1)++" and "++ (show $ dims t2)
+  where ok = (compat <$> getName t1 n <*> getName t2 n) == Just True
+        (ds1,m1) = firstIdx n t1
+        (ds2,m2) = firstIdx n t2
+        m = (trans m1) `multiply` m2
+
+common2 t1 t2 = [ n1 | n1 <- names t1, n2 <- names t2, n1==n2]
+
+infixl 5 |*|
+-- | Tensor product with automatic contraction of repeated indices, following Einstein summation convention.
+(|*|) :: (Coord t, Compat i)
+      => NArray i t -> NArray i t -> NArray i t
+t1 |*| t2 = r where
+    cs = common2 t1 t2
+    r = case cs of
+        [] -> rawProduct t1 t2
+        n:_ -> reorder orig $ contract (contract2 t1 t2 n)
+    orig = nub (names t1 ++ names t2) \\ cs
+
+-------------------------------------------------------------
+
+-- | Check if two arrays have the same structure.
+sameStructure :: (Eq i) => NArray i t1 -> NArray i t2 -> Bool
+sameStructure a b = sortBy (compare `on` iName) (dims a) == sortBy (compare `on` iName) (dims b)
+
+-------------------------------------------------------------
+
+-- | Apply a function on vectors to conformant arrays. Two arrays are 'conformant' if
+-- the dimensional structure of one of them is contained in the other one. The smaller
+-- structure is replicated along the extra dimensions. The result has the same index
+-- order as the largest structure (or as the first argument, if they are equal).
+zipArray :: (Coord a, Coord b, Compat i)
+   => (Vector a -> Vector b -> Vector c) -- ^ transformation
+   -> NArray i a
+   -> NArray i b
+   -> NArray i c
+zipArray o a b = liftNA2 o a' b' where
+    (a',b') = makeConformantT (a,b)
+
+-------------------------------------------------------
+
+showBases x = f $ concatMap (shbld) x
+    where   shbld (c,[]) = shsign c ++ showc c
+            shbld (c,l) = shsign c ++ g (showc c) ++ "{"++ concatMap show l++"}"
+            shsign c = if c < 0 then " - " else " + "
+            showc c
+                | abs (fromIntegral (round c :: Int) - c) <1E-10  = show (round $ abs c::Int)
+                | otherwise = printf "%.3f" (abs c)
+            f (' ':'+':' ':rs) = rs
+            f (' ':'-':' ':rs) = '-':rs
+            f a = a
+            g "1" = ""
+            g a = a
+
+---------------------------------------------------------
+
+data Rect = Rect { li :: Int, co :: Int, els :: [String] }
+
+rect s = pad r c (Rect r 0 ss)
+    where ss = lines s
+          r  = length ss
+          c  = maximum (map length ss)
+
+pad nr nc (Rect r c ss) = Rect (r+r') (c+c') ss'' where
+    r' = max 0 (nr-r)
+    c' = max 0 (nc-c)
+    ss' = map (padH nc) ss
+    ss'' = replicate r' (replicate nc '-') ++ ss'
+    padH l s = take (l-length s) (" | "++repeat ' ') ++ s
+
+dispH :: Int -> [Rect] -> Rect
+dispH k rs = Rect nr nc nss where
+    nr = maximum (map li rs)
+    nss' = mapTail (\x-> pad nr (co x + k) x) rs
+    nss = foldl1' (zipWith (++)) (map els nss')
+    nc = length (head nss)
+
+dispV :: Int -> [Rect] -> Rect
+dispV k rs = Rect nr nc nss where
+    nc = maximum (map co rs)
+    nss' = mapTail (\x-> pad (li x + k) nc x) rs
+    nss = concatMap els nss'
+    nr = length nss
+
+mapTail f (a:b) = a : map f b
+mapTail _ x     = x
+
+
+
+
+formatAux f x = unlines . addds . els . fmt ms $ x where
+    fmt [] _ = undefined -- cannot happen
+    fmt (g:gs) t
+        | rank t == 0 = rect (f (coords t @> 0))
+        | rank t == 1 =  rect $ unwords $ map f (toList $ coords t)
+        | rank t == 2 =  decor t $ rect $ w1 $ format " " f (reshape (iDim $ last $ dims t) (coords t))
+        | otherwise    = decor t (g ps)
+      where ps = map (fmt gs ) (partsRaw t (head (names t)))
+    ds = showNice (filter ((/='*').head.iName) $ dims x)
+    addds = if null ds then (showRawDims (dims x) :) else (ds:)
+    w1 = unlines . map (' ':) . lines
+    ms = cycle [dispV 1, dispH 2]
+    decor t | odd (rank t) = id
+            | otherwise = decorLeft  (names t!!0) . decorUp (names t!!1)
+
+
+showNice x = unwords . intersperse "x" . map show $ x
+showRawDims = showNice . map iDim . filter ((/="*").iName)
+
+------------------------------------------------------
+
+-- | Show a multidimensional array as a nested 2D table.
+formatArray :: (Coord t, Compat i)
+      => (t -> String) -- ^ format function (eg. printf \"5.2f\")
+      -> NArray i t
+      -> String
+formatArray f t | odd (rank t) = formatAux f (dummyAt 0 t)
+            | otherwise    = formatAux f t
+
+
+decorUp s rec
+    | head s == '*' = rec
+    | otherwise     = dispV 0 [rs,rec]
+  where
+    c = co rec
+    c1 = (c - length s) `div` 2
+    c2 = c - length s - c1
+    rs = rect $ replicate c1 ' ' ++ s ++ replicate c2 ' '
+
+decorLeft s rec
+    | head s == '*' = rec
+    | otherwise     = dispH 0 [rs,rec]
+  where
+    c = li rec
+    r1 = (c - length s+1) `div` 2
+    r2 = c - length s - r1
+    rs = rect $ unlines $ replicate r1 spc ++ s : replicate (r2) spc
+    spc = replicate (length s) ' '
+
+------------------------------------------------------
+
+-- | Print the array as a nested table with the desired format (e.g. %7.2f) (see also 'formatArray', and 'formatScaled').
+printA :: (Coord t, Compat i, PrintfArg t) => String -> NArray i t -> IO ()
+printA f t = putStrLn (formatArray (printf f) t)
+
+
+-- | Show the array as a nested table with autoscaled entries.
+formatScaled :: (Compat i)
+      => Int -- ^ number of of decimal places
+      -> NArray i Double
+      -> String
+formatScaled dec t = unlines (('(':d++")  E"++show o) : m)
+    where ss = formatArray (printf fmt. g) t
+          d:m = lines ss
+          g x = x/10^(o::Int)
+          o = floor $ maximum $ map (logBase 10 . abs) $ toList $ coords t
+          fmt = '%':show (dec+3) ++ '.':show dec ++"f"
+
+-- | Show the array as a nested table with a \"\%.nf\" format. If all entries
+-- are approximate integers the array is shown without the .00.. digits.
+formatFixed :: (Compat i)
+      => Int -- ^ number of of decimal places
+      -> NArray i Double
+      -> String
+formatFixed dec t
+    | isInt t   = formatArray (printf ('%': show (width t) ++".0f")) t
+    | otherwise = formatArray (printf ('%': show (width t+dec+1) ++"."++show dec ++"f")) t
+
+isInt = all lookslikeInt . toList . coords
+lookslikeInt x = show (round x :: Int) ++".0" == shx || "-0.0" == shx
+    where shx = show x
+-- needsSign t = vectorMin (coords t) < 0
+-- width :: Compat i => NArray i Double -> Int
+width = maximum . map (length . (printf "%.0f"::Double->String)) . toList . coords
+-- width t = k + floor (logBase 10 (max 1 $ vectorMax (abs $ coords t))) :: Int
+--      where k | needsSign t = 2
+--              | otherwise = 1
+
+------------------------------------------------------
+
+-- | Create an array from a list of subarrays. (The inverse of 'parts'.)
+newIndex:: (Coord t, Compat i) =>
+     i  -- ^ index type
+     -> Name
+     -> [NArray i t]
+     -> NArray i t
+newIndex i name ts = r where
+    ds = Idx (length ts) name i : (dims (head cts))
+    cts = makeConformant ts
+    r = mkNArray ds (join $ map coords cts)
+
+
+-- | Insert a dummy index of dimension 1 at a given level (for formatting purposes).
+dummyAt :: Int -> NArray i t -> NArray i t
+dummyAt k t = mkNArray d' (coords t) where
+    (d1,d2) = splitAt k (dims t)
+    d' = d1 ++ d : d2
+    d = Idx 1 "*" undefined
+
+-- | Rename indices so that they are not shown in formatted output.
+noIdx :: Compat i => NArray i t -> NArray i t
+noIdx t = renameRaw t (map ('*':) (names t))
+
+-- | Obtain a canonical base for the array.
+basisOf :: Coord t => NArray i t -> [NArray i t]
+basisOf t = map (dims t `mkNArray`) $ toRows (ident . dim . coords $ t)
+
+-------------------------------------------------------------
+
+instance (Coord t, Coord (Complex t), Compat i, Container Vector t) => Container (NArray i) t where
+    toComplex (r,c) = zipArray (curry toComplex) r c
+    fromComplex t = let (r,c) = fromComplex (coords t)
+                     in (mapArray (const r) t, mapArray (const c) t)
+    comp = mapArray comp
+    conj = mapArray conj
+    real = mapArray real
+    complex = mapArray complex
+
+----------------------------------------------------------------------
+
+-- | obtains the common value of a property of a list
+common :: (Eq a) => (b->a) -> [b] -> Maybe a
+common f = commonval . map f where
+    commonval :: (Eq a) => [a] -> Maybe a
+    commonval [] = Nothing
+    commonval [a] = Just a
+    commonval (a:b:xs) = if a==b then commonval (b:xs) else Nothing
+
+------------------------------------------------------------------------
+
+-- | Extract the 'Matrix' corresponding to a two-dimensional array,
+-- in the rows,cols order.
+asMatrix :: (Coord t) => NArray i t -> Matrix t
+asMatrix a | rank a == 2 = reshape c (coords a)
+           | otherwise = error $ "asMatrix requires a rank 2 array."
+    where c = size (last (names a)) a
+
+-- | Extract the 'Vector' corresponding to a one-dimensional array.
+asVector :: (Coord t) => NArray i t -> Vector t
+asVector a | rank a == 1 = coords a
+           | otherwise = error $ "asVector requires a rank 1 array."
+
+-- | Extract the scalar element corresponding to a 0-dimensional array.
+asScalar :: (Coord t) => NArray i t -> t
+asScalar a | rank a == 0 = coords a @>0
+           | otherwise = error $ "asScalar requires a rank 0 array."
+
+------------------------------------------------------------------------
+
+-- | Create a rank-1 array from an hmatrix 'Vector'.
+fromVector :: Compat i => i -> Vector t -> NArray i t
+fromVector i v = mkNArray [Idx (dim v) "1" i ] v
+
+-- | Create a rank-2 array from an hmatrix 'Matrix'.
+fromMatrix :: (Compat i, Coord t) => i -> i -> Matrix t -> NArray i t
+fromMatrix ir ic m = mkNArray [Idx (rows m) "1" ir,
+                               Idx (cols m) "2" ic] (flatten m)
+
+------------------------------------------------------------------------
+
+-- | Select some parts of a tensor, taking into account position and value.
+extract :: (Compat i, Coord t)
+        => (Int -> NArray i t -> Bool)
+        -> Name
+        -> NArray i t
+        -> NArray i t
+extract f name arr = reorder (names arr)
+                   . newIndex (typeOf name arr) name
+                   . map snd . filter (uncurry f)
+                   $ zip [1..] (parts arr name)
+
+-- | Apply a list function to the parts of an array at a given index.
+onIndex :: (Coord a, Coord b, Compat i) =>
+     ([NArray i a] -> [NArray i b])
+     -> Name
+     -> NArray i a
+     -> NArray i b
+onIndex f name t = reorder (names t) $ newIndex (typeOf name t) name (f (parts t name))
+
+------------------------------------------------------------------------
+
+extend alldims (A d v) = reorder (allnames) s where
+    allnames = map iName alldims
+    pref = alldims \\ d
+    n = product (map iDim pref)
+    s = A (pref++d) (join (replicate n v))
+
+-- | Obtains most general structure of a list of dimension specifications
+conformable :: Compat i => [[Idx i]] -> Maybe [Idx i]
+conformable ds | ok        = Just alldims
+               | otherwise = Nothing
+    where alldims = nub (concat ds)
+          allnames = map iName alldims
+          ok = length (allnames) == length (nub allnames)
+
+-- | Converts a list of arrays to a common structure.
+makeConformant :: (Coord t, Compat i) => [NArray i t] -> [NArray i t]
+makeConformant ts =
+    case conformable (map dims ts) of
+        Just alldims -> map (extend alldims) ts
+        Nothing -> error $ "makeConformant with inconsistent dimensions "
+                         ++ show (map dims ts)
+
+-- the same version for tuples with possibly different element types
+makeConformantT (t1,t2) =
+    case conformable [dims t1, dims t2] of
+        Just alldims -> (extend alldims t1, extend alldims t2)
+        Nothing -> error $ "makeConformantT with inconsistent dimensions "
+                         ++ show (dims t1, dims t2)
diff --git a/lib/Numeric/LinearAlgebra/Array/Simple.hs b/lib/Numeric/LinearAlgebra/Array/Simple.hs
new file mode 100644
--- /dev/null
+++ b/lib/Numeric/LinearAlgebra/Array/Simple.hs
@@ -0,0 +1,55 @@
+{-# LANGUAGE FlexibleInstances, FlexibleContexts, TypeSynonymInstances #-}
+{-# OPTIONS_HADDOCK hide #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.Array.Simple
+-- Copyright   :  (c) Alberto Ruiz 2009
+-- License     :  GPL
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-- Simple multidimensional arrays.
+-- Contractions only require equal dimension.
+--
+-----------------------------------------------------------------------------
+
+module Numeric.LinearAlgebra.Array.Simple (
+    None(..),
+    Array,
+    listArray
+) where
+
+import Numeric.LinearAlgebra.Array.Internal
+import Data.Packed
+
+
+instance Show (Idx None) where
+    show (Idx n s _t) = show n ++ ":" ++ s
+
+-- | Unespecified coordinate type. Contractions only
+-- require equal dimension.
+data None = None deriving Eq
+
+
+instance Compat None where
+    compat d1 d2 = iDim d1 == iDim d2
+
+
+-- | Multidimensional array with unespecified coordinate type.
+type Array t = NArray None t
+
+instance (Coord t) => Show (Array t) where
+    show t | null (dims t) = "scalar "++ show (coords t @>0)
+           | otherwise = "listArray "++ show (dims t) ++ " "++ show (toList $ coords t)
+
+-- | Construction of an 'Array' from a list of dimensions and a list of elements in left to right order.
+listArray :: (Coord t)
+    => [Int] -- ^ dimensions
+    -> [t]   -- ^ elements
+    -> Array t
+listArray ds cs = mkNArray dms (product ds |> (cs ++ repeat 0))
+    where dms = zipWith3 Idx ds (map show [1::Int ..]) (repeat None)
+
+
diff --git a/lib/Numeric/LinearAlgebra/Array/Util.hs b/lib/Numeric/LinearAlgebra/Array/Util.hs
new file mode 100644
--- /dev/null
+++ b/lib/Numeric/LinearAlgebra/Array/Util.hs
@@ -0,0 +1,44 @@
+-- {-# LANGUAGE FlexibleInstances, FlexibleContexts #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.Array.Util
+-- Copyright   :  (c) Alberto Ruiz 2009
+-- License     :  GPL
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-- Additional tools for manipulation of multidimensional arrays.
+--
+-----------------------------------------------------------------------------
+
+module Numeric.LinearAlgebra.Array.Util (
+    Coord, Compat(..),
+    NArray, Idx(..), Name,
+    scalar,
+    rank, names, size, typeOf, dims, coords,
+
+    rename, (!),
+
+    parts,
+    newIndex,
+
+    mapArray, zipArray, (|*|),
+
+    extract, onIndex,
+
+    reorder, (~>),
+    formatArray, formatFixed, formatScaled,
+    dummyAt, noIdx,
+    conformable,
+    sameStructure,
+    makeConformant,
+    basisOf,
+    asScalar, asVector, asMatrix,
+    fromVector, fromMatrix,
+    Container(..),
+) where
+
+import Numeric.LinearAlgebra.Array.Internal
+import Data.Packed(Container(..))
diff --git a/lib/Numeric/LinearAlgebra/Exterior.hs b/lib/Numeric/LinearAlgebra/Exterior.hs
new file mode 100644
--- /dev/null
+++ b/lib/Numeric/LinearAlgebra/Exterior.hs
@@ -0,0 +1,128 @@
+{-# LANGUAGE FlexibleContexts #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.LinearAlgebra.Exterior
+-- Copyright   :  (c) Alberto Ruiz 2009
+-- License     :  GPL
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  experimental
+--
+-- Exterior Algebra.
+--
+--
+-----------------------------------------------------------------------------
+
+module Numeric.LinearAlgebra.Exterior (
+    (/\),
+    inner,
+    leviCivita,
+    dual,
+    (\/),
+    module Numeric.LinearAlgebra.Tensor,
+    asMultivector, fromMultivector
+) where
+
+import Numeric.LinearAlgebra.Tensor
+import Numeric.LinearAlgebra.Array.Internal
+import Numeric.LinearAlgebra.Multivector(Multivector,fromTensor,maxDim,grade)
+import qualified Numeric.LinearAlgebra.Multivector as MV
+import Data.List
+
+-- import Debug.Trace
+-- debug x = trace (show x) x
+
+interchanges :: (Ord a) => [a] -> Int
+interchanges ls = sum (map (count ls) ls)
+    where count l p = length $ filter (>p) $ take pel l
+              where Just pel = elemIndex p l
+
+signature :: (Num t, Ord a) => [a] -> t
+signature l | length (nub l) < length l =  0
+            | even (interchanges l)     =  1
+            | otherwise                 = -1
+
+gsym f t = mkNArray (dims t) (coords $ sum ts) where
+    ns = map show [1 .. rank t]
+    t' = cov $ renameRaw t ns
+    per = permutations ns
+    ts  = map (flip renameRaw ns . f . flip reorder t') per
+
+-- symmetrize t = gsym id t
+
+antisymmetrize t = gsym scsig t
+    where scsig x = scalar (signature (names x)) * x
+
+fact n = product [1..n]
+
+wedge a b = antisymmetrize (a*b) * (recip . fromIntegral) (fact (rank a) * fact (rank b))
+
+infixl 5 /\
+-- | The exterior (wedge) product of two tensors. Obtains the union of subspaces.
+--
+--   Implemented as the antisymmetrization of the tensor product.
+(/\) :: (Coord t)
+     => Tensor t
+     -> Tensor t
+     -> Tensor t
+a /\ b = renseq (wedge a' b')
+    where a' = renseq  a
+          b' = renseq' b
+
+-- levi n = antisymmetrize $ product $ zipWith renameRaw ts is
+--     where is = map (return.show) [1 .. n]
+--           ts = map (listTensor [n]) (toLists $ ident n)
+
+levi n = listTensor (replicate n n) $ map signature $ sequence (replicate n [1..n])
+
+-- | The full antisymmetric tensor of rank n (contravariant version).
+leviCivita :: Int -> Tensor Double
+leviCivita = (map levi [0..] !!)
+
+infixl 4 \/
+-- | The \"meet\" operator. Obtains the intersection of subspaces.
+--
+-- @a \\\/ b = dual (dual a \/\\ dual b)@
+(\/) :: Tensor Double -> Tensor Double -> Tensor Double
+a \/ b = dual (dual a /\ dual b)
+
+dual' n t = inner (leviCivita n) t
+
+-- | Inner product of a r-vector with the whole space.
+--
+-- @dual t = inner (leviCivita n) t@
+dual :: Tensor Double -> Tensor Double
+dual t | isScalar t = error $ "cannot deduce dimension for dual of a scalar. Use s * leviCivita n"
+       | otherwise  = dual' n t
+    where n = case common iDim (dims t) of
+                Just x -> x
+                Nothing -> error $ "dual with different dimensions"
+
+-- | Euclidean inner product of multivectors.
+inner :: (Coord t)
+      => Tensor t
+      -> Tensor t
+      -> Tensor t
+inner a b | rank a < rank b = switch (renseq a) * renseq b * k
+          | otherwise       = renseq a * switch (renseq b) * k
+    where k = recip . fromIntegral $ fact $ min (rank a) (rank b)
+
+renseq t = renameRaw t (map show [1..rank t])
+renseq' t = renameRaw t (map ((' ':).show) [1..rank t])
+
+isScalar = null . dims
+
+-- | Extract a compact multivector representation from a full antisymmetric tensor.
+--
+-- asMultivector = Multivector.'fromTensor'.
+--
+-- (We do not check that the tensor is actually antisymmetric.)
+asMultivector :: Tensor Double -> Multivector
+asMultivector = fromTensor
+
+-- | Create an explicit antisymmetric 'Tensor' from the components of a Multivector of a given grade.
+fromMultivector :: Int -> Multivector -> Tensor Double
+fromMultivector k t = sum $ map f (MV.coords $ grade k t) where
+    f (x,es) = scalar x * foldl1' (/\) (map g es)
+    n = maxDim t
+    g i = vector $ replicate (i-1) 0 ++ 1 : replicate (n-i) 0
diff --git a/lib/Numeric/LinearAlgebra/Multivector.hs b/lib/Numeric/LinearAlgebra/Multivector.hs
new file mode 100644
--- /dev/null
+++ b/lib/Numeric/LinearAlgebra/Multivector.hs
@@ -0,0 +1,265 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.LinearAlgebra.Multivector
+-- Copyright   :  (c) Alberto Ruiz 2009
+-- License     :  GPL
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  experimental
+--
+-- A simple implementation of Geometric Algebra.
+--
+-- The Num instance provides the geometric product, and the Fractional
+-- instance provides the inverse of multivectors.
+--
+-- This module provides a simple Euclidean embedding.
+
+-----------------------------------------------------------------------------
+
+module Numeric.LinearAlgebra.Multivector (
+    Multivector, coords,
+    scalar, vector, e, (/\), (-|), (\/), rever, full, rotor,
+    apply,
+    grade, maxGrade, maxDim,
+    fromTensor
+) where
+
+import Numeric.LinearAlgebra(toList,reshape,(<\>),(@>))
+import Numeric.LinearAlgebra.Array.Internal hiding (scalar,coords)
+import Numeric.LinearAlgebra.Tensor hiding (scalar,vector)
+import qualified Numeric.LinearAlgebra.Array.Internal as Array
+import Data.List
+import Control.Monad(filterM)
+import Data.Function(on)
+import qualified Data.Map as Map
+
+powerset = filterM (const [True, False])  -- !!
+
+base :: Int -> [[Int]]
+base k = sortBy (compare `on` length) (powerset [1..k])
+
+base' k = map (\b -> MV [(1,b)]) (base k)
+
+data Multivector = MV { coords :: [(Double,[Int])] } deriving Eq
+
+instance Show Multivector where
+    show = showMV
+
+maxGrade :: Multivector -> Int
+maxGrade (MV l) = maximum . map (length.snd) $ l
+
+grade :: Int -> Multivector -> Multivector
+grade k (MV l) = MV $ filter ((k==).length.snd) l
+
+maxDim :: Multivector -> Int
+maxDim (MV [(_,[])]) = 0
+maxDim (MV l) = maximum . concat . map snd $ l
+
+-- | The reversion operator.
+rever :: Multivector -> Multivector
+rever (MV l) = MV (map r l) where
+    r (c,b) = (c*fromIntegral s ,b)
+        where s = signum (-1)^(k*(k-1)`div`2) :: Int
+              k = length b
+
+-- | Show the non zero coordinates of a multivector in a nicer format.
+showMV :: Multivector -> String
+showMV (MV x) = showBases x
+
+-- | Creates a scalar multivector.
+scalar :: Double -> Multivector
+scalar s = MV [(s,[])]
+
+-- | Creates a grade 1 multivector of from a list of coordinates.
+vector :: [Double] -> Multivector
+vector v = MV $ simplify $ zip v (map (:[]) [1..])
+
+
+-- different product rules
+
+-- reorders the base indices remembering the original position
+r1 :: [Int] -> [(Int,[Int])]
+r1 [] = []
+r1 l = (m,elemIndices m l):(r1 (filter (/=m) l))
+    where m = minimum l
+
+-- geometric product
+r2 :: [(Int, [Int])] -> (Double, [Int])
+r2 = foldl' g (1,[])
+    where g (k,l) (x,ps) = (k*s,l++t)
+              where t = if even (length ps) then [] else [x]
+                    s = product (map f ps')
+                         where f z = if even z then 1 else -1
+                               ps' = zipWith (subtract) ps [0..]
+
+-- exterior product
+r3 :: [(Int, [Int])] -> (Double, [Int])
+r3 = foldl' g (1,[])
+    where g (k,l) (x,ps) = (k*s,l++[x])
+              where s = if length ps > 1 then 0 else if even (head ps) then 1 else -1
+
+
+-- simplification and cleaning of the list of coordinates
+simplify = chop . grp . sortBy (compare `on` snd)
+    where grp [] = []
+          grp [a] = [a]
+          grp ((c1,b1):(c2,b2):rest)
+              | b1 == b2  = grp ( (c1+c2,b1) : rest)
+              | otherwise = (c1,b1): grp ((c2,b2):rest)
+          zero (c,_) = abs c < 1E-8
+          chop = cz . filter (not.zero)
+          cz [] = [(0,[])]
+          cz x  = x
+
+-- sum of multivectors
+gs (MV l1) (MV l2) = MV $ simplify (l1++l2)
+
+-- geometric product
+gp (MV l1) (MV l2) = MV $ simplify [g x y | x<-l1, y <-l2]
+    where g (c1,b1) (c2,b2) = (k*c1*c2,b3) where (k,b3) = gpr b1 b2 --(r2.r1) (b1++b2)
+
+-- exterior product
+ge (MV l1) (MV l2) = MV $ simplify [g x y | x<-l1, y <-l2]
+    where g (c1,b1) (c2,b2) = (k*c1*c2,b3) where (k,b3) = epr b1 b2 -- (r3.r1) (b1++b2)
+
+-- contraction inner product
+gi (MV l1) (MV l2) = sum [g x y | x<-l1, y <-l2]
+    where g (c1,[]) (c2,is) = MV [(c1*c2,is)]
+          g _ (_,[])        = 0
+
+          g (c1,[i]) (c2,[j]) = if i==j then MV [(c1*c2,[])] else 0
+          g (c1,[i]) (c2,j:js) = (g (c1,[i]) (c2,[j]) /\ MV [(1,js)])
+                               - (MV [(c2,[j])] /\ g (c1,[i]) (1,js))
+
+          g (c1,i:is) b = gi (MV [(c1,[i])]) (gi (MV[(1,is)]) (MV [b]))
+
+
+instance Num Multivector where
+    (+) = gs
+    (*) = gp
+    negate (MV l) = MV (map neg l) where neg (k,b) = (-k,b)
+    abs _ = error "abs of multivector not yet defined"
+    signum _ = error "signum of multivector not yet defined"
+    fromInteger x = MV [(fromInteger x,[])]
+
+instance Fractional Multivector where
+    fromRational x = MV [(fromRational x,[])]
+    recip (MV [(x,[])]) = MV [(recip x,[])]
+    recip x = mvrecip x
+
+-- | The k-th basis element.
+e :: Int -> Multivector
+e k = MV [(1,[k])]
+
+-- | The exterior (outer) product.
+(/\) :: Multivector -> Multivector -> Multivector
+infixl 7 /\
+(/\) = ge
+
+
+-- | The contractive inner product.
+(-|) :: Multivector -> Multivector -> Multivector
+infixl 7 -|
+(-|) = gi
+
+-- | The full space of the given dimension. This is the leviCivita simbol, and the basis of the pseudoscalar.
+full :: Int -> Multivector
+full k = MV [(1,[1 .. k])] --product . map e $ [1 .. k]
+
+-- | Intersection of subspaces.
+(\/) :: Multivector -> Multivector -> Multivector
+infixl 7 \/
+(\/) a b = (b -| rever (full k)) -| a
+    where k = max (maxDim a) (maxDim b)
+
+-- check that it is a vector
+normVec v = sqrt x where MV [(x,[])] = v * v
+
+unitary v = v / scalar (normVec v)
+
+-- | The rotor operator, used in a sandwich product.
+rotor :: Int          -- ^ dimension of the space
+      -> Double       -- ^ angle
+      -> Multivector  -- ^ axis
+      -> Multivector  -- ^ result
+rotor k phi axis = scalar (cos (phi/2)) - scalar (sin (phi/2)) * (unitary axis*full k)
+
+
+-- memoization of the rules
+gprules k = Map.fromList [(x, Map.fromList [(y,(r2.r1)(x++y)) | y<-base k] )| x<-base k]
+
+eprules k = Map.fromList [(x, Map.fromList [(y,(r3.r1)(x++y)) | y<-base k] )| x<-base k]
+
+--reasonable limit
+gpr a b = g Map.! a Map.! b
+    where g = gprules 6
+
+epr a b = g Map.! a Map.! b
+    where g = eprules 6
+
+----------------------- tensor expansion -----------------------
+
+expand k = g
+    where g (MV l) = foldl1' (zipWith (+)) $ map f l
+          basepos b = m Map.! b
+          m = baseraw k
+          f (c,b) = en pk (basepos b) c
+          pk = 2^k
+          baseraw q = Map.fromList $ zip (base q) [0..]
+          en n q v = replicate q 0 ++ v : replicate (n-q-1) 0
+
+compact k t = sum $ zipWith (*) (map scalar $ toList (Array.coords t)) (base' k)
+
+gatensor k = listTensor [-pk,-pk,pk] (concat . concat $ gacoords)
+    where pk = 2^k
+          gacoords = [[ f (x * y) | y<-b] | x<-b]
+          b = base' k
+          f = expand k
+
+tmv k x = listTensor [2^k] (expand k x)
+
+
+-- tp a b = comp (g!"ijk" * ta!"i" * tb!"j")
+--     where k = max (maxDim a) (maxDim b)
+--           g = gatensor k
+--           ta = tmv k a
+--           tb = tmv k b
+--           comp = compact k
+
+mat rowidx t = reshape c $ Array.coords t'
+    where c = iDim $ last (dims t')
+          t' = reorder (rowidx: (names t\\[rowidx])) t
+
+-- on the right
+pmat k b = mat "k" $ g!"ijk" * tb!"j"
+    where g = gatensor k
+          tb = tmv k b
+
+divi k a b = compact k $ listTensor [2^k] (toList $ pmat k b <\> Array.coords (tmv k a))
+
+mvrecip b = divi (maxDim b) 1 b
+
+--------------------------------------------------------
+
+-- | Extract a multivector representation from a full antisymmetric tensor.
+--
+-- (We do not check that the tensor is actually antisymmetric.)
+fromTensor :: Tensor Double -> Multivector
+fromTensor t = MV $ filter ((/=0.0).fst) $ zip vals basis
+    where vals = map ((@> 0). Array.coords .foldl' partF t) (map (map pred) basis)
+          r = length (dims t)
+          n = iDim . head . dims $ t
+          partF s i = part s (name,i) where name = iName . head . dims $ s
+          basis = filter (\x-> (x==nub x && x==sort x)) $ sequence $ replicate r [1..n]
+
+part t (name,i) = parts t name !! i
+
+--------------------------------------------------------
+
+-- | Apply a linear transformation, expressed as the image of the element i-th of the basis.
+--
+--  (This is a monadic bind!)
+apply :: (Int -> Multivector) -> Multivector -> Multivector
+apply f t = sum $ map g (coords t) where
+    g (x,[]) = scalar x
+    g (x,es) = scalar x * foldl1' (/\) (map f es)
diff --git a/lib/Numeric/LinearAlgebra/Tensor.hs b/lib/Numeric/LinearAlgebra/Tensor.hs
new file mode 100644
--- /dev/null
+++ b/lib/Numeric/LinearAlgebra/Tensor.hs
@@ -0,0 +1,102 @@
+{-# LANGUAGE FlexibleInstances, FlexibleContexts, TypeSynonymInstances #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.LinearAlgebra.Tensor
+-- Copyright   :  (c) Alberto Ruiz 2009
+-- License     :  GPL
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  experimental
+--
+-- Tensor computations. Indices can only be contracted if they are of different 'Variant' type.
+--
+-----------------------------------------------------------------------------
+
+
+module Numeric.LinearAlgebra.Tensor (
+    -- * The Tensor type
+    Tensor, Variant(..),
+    listTensor,
+    -- * Tensor creation utilities
+    superindex, subindex,
+    vector, covector, transf,
+    -- * Index manipulation
+    switch, cov, contrav, forget,
+    -- * General array operations
+    module Numeric.LinearAlgebra.Array
+) where
+
+import Numeric.LinearAlgebra.Array.Internal
+import Numeric.LinearAlgebra hiding (rank)
+import Numeric.LinearAlgebra.Array
+
+type Tensor t = NArray Variant t
+
+data Variant = Co | Contra deriving (Eq)
+
+instance Compat Variant where
+    compat d1 d2 = iDim d1 == iDim d2 && iType d1 /= iType d2
+
+instance Show (Idx Variant) where
+    show (Idx n s Co)     = show n ++ "_" ++ s
+    show (Idx n s Contra) = show n ++ "^" ++ s
+
+instance (Coord t) => Show (Tensor t) where
+    show t | null (dims t) = show (coords t @>0)
+           | otherwise = "listTensor " ++ show (dims t) ++ " "++ show (toList (coords t))
+
+flipV Co = Contra
+flipV Contra = Co
+
+-- | Creates a tensor from a list of dimensions and a list of coordinates.
+-- A positive dimension means that the index is assumed to be contravariant (vector-like), and
+-- a negative dimension means that the index is assumed to be covariant (like a linear function, or covector). Contractions can only be performed between indices of different type.
+listTensor :: Coord t
+           => [Int] -- ^ dimensions
+           -> [t]   -- ^ coordinates
+           -> Tensor t
+listTensor ds cs = mkNArray dms (product ds' |> (cs ++ repeat 0))
+    where dms = zipWith3 Idx ds' (map show [1::Int ..]) (map f ds)
+          ds' = map abs ds
+          f n | n>0       = Contra
+              | otherwise = Co
+
+-- | Create an 'Tensor' from a list of parts with a contravariant index (@superindex = 'newIndex' 'Contra'@).
+superindex :: Coord t => Name -> [Tensor t] -> Tensor t
+superindex = newIndex Contra
+
+-- | Create an 'Tensor' from a list of parts with a covariant index (@subindex = 'newIndex' 'Co'@).
+subindex :: Coord t => Name -> [Tensor t] -> Tensor t
+subindex   = newIndex Co
+
+
+
+-- | Change the 'Variant' nature of all dimensions to the opposite ones.
+switch :: Tensor t -> Tensor t
+switch = mapTypes flipV
+
+-- | Make all dimensions covariant.
+cov :: NArray i t -> Tensor t
+cov     = mapTypes (const Co)
+
+-- | Make all dimensions contravariant.
+contrav :: NArray i t -> Tensor t
+contrav = mapTypes (const Contra)
+
+-- | Remove the 'Variant' nature of coordinates.
+forget :: NArray i t -> Array t
+forget = mapTypes (const None)
+
+--------------------------------------------------------------
+
+-- | Create a contravariant rank-1 tensor from a list of coordinates.
+vector :: [Double] -> Tensor Double
+vector   = fromVector Contra . fromList
+
+-- | Create a covariant rank-1 tensor from a list of coordinates.
+covector :: [Double] -> Tensor Double
+covector = fromVector Co . fromList
+
+-- | Create a 1-contravariant, 1-covariant rank-2 tensor from list of lists of coordinates.
+transf :: [[Double]] -> Tensor Double
+transf   = fromMatrix Contra Co . fromLists