From 1f8f53ed0f1374631882e7d6c816b05d029edb4e Mon Sep 17 00:00:00 2001
From: Alberto Ruiz <aruiz@um.es>
Date: Fri, 13 Jun 2014 10:12:41 +0200
Subject: move modules (I)

---
 packages/base/src/Numeric/HMatrix.hs               | 619 ---------------------
 packages/base/src/Numeric/LinearAlgebra/Static.hs  | 422 --------------
 .../src/Numeric/LinearAlgebra/Static/Internal.hs   | 422 ++++++++++++++
 .../base/src/Numeric/LinearAlgebra/tmpStatic.hs    | 619 +++++++++++++++++++++
 4 files changed, 1041 insertions(+), 1041 deletions(-)
 delete mode 100644 packages/base/src/Numeric/HMatrix.hs
 delete mode 100644 packages/base/src/Numeric/LinearAlgebra/Static.hs
 create mode 100644 packages/base/src/Numeric/LinearAlgebra/Static/Internal.hs
 create mode 100644 packages/base/src/Numeric/LinearAlgebra/tmpStatic.hs

diff --git a/packages/base/src/Numeric/HMatrix.hs b/packages/base/src/Numeric/HMatrix.hs
deleted file mode 100644
index c88918f..0000000
--- a/packages/base/src/Numeric/HMatrix.hs
+++ /dev/null
@@ -1,619 +0,0 @@
-{-# LANGUAGE DataKinds #-}
-{-# LANGUAGE KindSignatures #-}
-{-# LANGUAGE GeneralizedNewtypeDeriving #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE FunctionalDependencies #-}
-{-# LANGUAGE FlexibleContexts #-}
-{-# LANGUAGE ScopedTypeVariables #-}
-{-# LANGUAGE EmptyDataDecls #-}
-{-# LANGUAGE Rank2Types #-}
-{-# LANGUAGE FlexibleInstances #-}
-{-# LANGUAGE TypeOperators #-}
-{-# LANGUAGE ViewPatterns #-}
-{-# LANGUAGE GADTs #-}
-{-# LANGUAGE OverlappingInstances #-}
-{-# LANGUAGE TypeFamilies #-}
-
-
-{- |
-Module      :  Numeric.HMatrix
-Copyright   :  (c) Alberto Ruiz 2006-14
-License     :  BSD3
-Stability   :  experimental
-
-Experimental interface with statically checked dimensions.
-
--}
-
-module Numeric.HMatrix(
-    -- * Vector
-    ℝ, R,
-    vec2, vec3, vec4, (&), (#), split, headTail,
-    vector,
-    linspace, range, dim,
-    -- * Matrix
-    L, Sq, build,
-    row, col, (¦),(——), splitRows, splitCols,
-    unrow, uncol,
-    tr,
-    eye,
-    diag,
-    blockAt,
-    matrix,
-    -- * Complex
-    C, M, Her, her, 𝑖,
-    -- * Products
-    (<>),(#>),(<·>),
-    -- * Linear Systems
-    linSolve, (<\>),
-    -- * Factorizations
-    svd, withCompactSVD, svdTall, svdFlat, Eigen(..),
-    withNullspace, qr,
-    -- * Misc
-    mean,
-    Disp(..), Domain(..),
-    withVector, withMatrix,
-    toRows, toColumns,
-    Sized(..), Diag(..), Sym, sym, mTm, unSym
-) where
-
-
-import GHC.TypeLits
-import Numeric.LinearAlgebra.HMatrix hiding (
-    (<>),(#>),(<·>),Konst(..),diag, disp,(¦),(——),row,col,vector,matrix,linspace,toRows,toColumns,
-    (<\>),fromList,takeDiag,svd,eig,eigSH,eigSH',eigenvalues,eigenvaluesSH,eigenvaluesSH',build,
-    qr)
-import qualified Numeric.LinearAlgebra.HMatrix as LA
-import Data.Proxy(Proxy)
-import Numeric.LinearAlgebra.Static
-import Control.Arrow((***))
-
-
-
-
-
-ud1 :: R n -> Vector ℝ
-ud1 (R (Dim v)) = v
-
-
-infixl 4 &
-(&) :: forall n . (KnownNat n, 1 <= n)
-    => R n -> ℝ -> R (n+1)
-u & x = u # (konst x :: R 1)
-
-infixl 4 #
-(#) :: forall n m . (KnownNat n, KnownNat m)
-    => R n -> R m -> R (n+m)
-(R u) # (R v) = R (vconcat u v)
-
-
-
-vec2 :: ℝ -> ℝ -> R 2
-vec2 a b = R (gvec2 a b)
-
-vec3 :: ℝ -> ℝ -> ℝ -> R 3
-vec3 a b c = R (gvec3 a b c)
-
-
-vec4 :: ℝ -> ℝ -> ℝ -> ℝ -> R 4
-vec4 a b c d = R (gvec4 a b c d)
-
-vector :: KnownNat n => [ℝ] -> R n
-vector = fromList
-
-matrix :: (KnownNat m, KnownNat n) => [ℝ] -> L m n
-matrix = fromList
-
-linspace :: forall n . KnownNat n => (ℝ,ℝ) -> R n
-linspace (a,b) = mkR (LA.linspace d (a,b))
-  where
-    d = fromIntegral . natVal $ (undefined :: Proxy n)
-
-range :: forall n . KnownNat n => R n
-range = mkR (LA.linspace d (1,fromIntegral d))
-  where
-    d = fromIntegral . natVal $ (undefined :: Proxy n)
-
-dim :: forall n . KnownNat n => R n
-dim = mkR (scalar d)
-  where
-    d = fromIntegral . natVal $ (undefined :: Proxy n)
-
-
---------------------------------------------------------------------------------
-
-
-ud2 :: L m n -> Matrix ℝ
-ud2 (L (Dim (Dim x))) = x
-
-
---------------------------------------------------------------------------------
-
-diag :: KnownNat n => R n -> Sq n
-diag = diagR 0
-
-eye :: KnownNat n => Sq n
-eye = diag 1
-
---------------------------------------------------------------------------------
-
-blockAt :: forall m n . (KnownNat m, KnownNat n) =>  ℝ -> Int -> Int -> Matrix Double -> L m n
-blockAt x r c a = mkL res
-  where
-    z = scalar x
-    z1 = LA.konst x (r,c)
-    z2 = LA.konst x (max 0 (m'-(ra+r)), max 0 (n'-(ca+c)))
-    ra = min (rows a) . max 0 $ m'-r
-    ca = min (cols a) . max 0 $ n'-c
-    sa = subMatrix (0,0) (ra, ca) a
-    m' = fromIntegral . natVal $ (undefined :: Proxy m)
-    n' = fromIntegral . natVal $ (undefined :: Proxy n)
-    res = fromBlocks [[z1,z,z],[z,sa,z],[z,z,z2]]
-
-
-
-
-
---------------------------------------------------------------------------------
-
-
-row :: R n -> L 1 n
-row = mkL . asRow . ud1
-
---col :: R n -> L n 1
-col v = tr . row $ v
-
-unrow :: L 1 n -> R n
-unrow = mkR . head . LA.toRows . ud2
-
---uncol :: L n 1 -> R n
-uncol v = unrow . tr $ v
-
-
-infixl 2 ——
-(——) :: (KnownNat r1, KnownNat r2, KnownNat c) => L r1 c -> L r2 c -> L (r1+r2) c
-a —— b = mkL (extract a LA.—— extract b)
-
-
-infixl 3 ¦
--- (¦) :: (KnownNat r, KnownNat c1, KnownNat c2) => L r c1 -> L r c2 -> L r (c1+c2)
-a ¦ b = tr (tr a —— tr b)
-
-
-type Sq n  = L n n
---type CSq n = CL n n
-
-type GL = (KnownNat n, KnownNat m) => L m n
-type GSq = KnownNat n => Sq n
-
-isKonst :: forall m n . (KnownNat m, KnownNat n) => L m n -> Maybe (ℝ,(Int,Int))
-isKonst (unwrap -> x)
-    | singleM x = Just (x `atIndex` (0,0), (m',n'))
-    | otherwise = Nothing
-  where
-    m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
-    n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
-
-
-isKonstC :: forall m n . (KnownNat m, KnownNat n) => M m n -> Maybe (ℂ,(Int,Int))
-isKonstC (unwrap -> x)
-    | singleM x = Just (x `atIndex` (0,0), (m',n'))
-    | otherwise = Nothing
-  where
-    m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
-    n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
-
-
-
-infixr 8 <>
-(<>) :: forall m k n. (KnownNat m, KnownNat k, KnownNat n) => L m k -> L k n -> L m n
-(<>) = mulR
-
-
-infixr 8 #>
-(#>) :: (KnownNat m, KnownNat n) => L m n -> R n -> R m
-(#>) = appR
-
-
-infixr 8 <·>
-(<·>) :: R n -> R n -> ℝ
-(<·>) = dotR
-
---------------------------------------------------------------------------------
-
-class Diag m d | m -> d
-  where
-    takeDiag :: m -> d
-
-
-instance forall n . (KnownNat n) => Diag (L n n) (R n)
-  where
-    takeDiag m = mkR (LA.takeDiag (extract m))
-
-
-instance forall m n . (KnownNat m, KnownNat n, m <= n+1) => Diag (L m n) (R m)
-  where
-    takeDiag m = mkR (LA.takeDiag (extract m))
-
-
-instance forall m n . (KnownNat m, KnownNat n, n <= m+1) => Diag (L m n) (R n)
-  where
-    takeDiag m = mkR (LA.takeDiag (extract m))
-
-
---------------------------------------------------------------------------------
-
-linSolve :: (KnownNat m, KnownNat n) => L m m -> L m n -> Maybe (L m n)
-linSolve (extract -> a) (extract -> b) = fmap mkL (LA.linearSolve a b)
-
-(<\>) :: (KnownNat m, KnownNat n, KnownNat r) => L m n -> L m r -> L n r
-(extract -> a) <\> (extract -> b) = mkL (a LA.<\> b)
-
-svd :: (KnownNat m, KnownNat n) => L m n -> (L m m, R n, L n n)
-svd (extract -> m) = (mkL u, mkR s', mkL v)
-  where
-    (u,s,v) = LA.svd m
-    s' = vjoin [s, z]
-    z = LA.konst 0 (max 0 (cols m - size s))
-
-
-svdTall :: (KnownNat m, KnownNat n, n <= m) => L m n -> (L m n, R n, L n n)
-svdTall (extract -> m) = (mkL u, mkR s, mkL v)
-  where
-    (u,s,v) = LA.thinSVD m
-
-
-svdFlat :: (KnownNat m, KnownNat n, m <= n) => L m n -> (L m m, R m, L n m)
-svdFlat (extract -> m) = (mkL u, mkR s, mkL v)
-  where
-    (u,s,v) = LA.thinSVD m
-
---------------------------------------------------------------------------------
-
-class Eigen m l v | m -> l, m -> v
-  where
-    eigensystem :: m -> (l,v)
-    eigenvalues :: m -> l
-
-newtype Sym n = Sym (Sq n) deriving Show
-
-
-sym :: KnownNat n => Sq n -> Sym n
-sym m = Sym $ (m + tr m)/2
-
-mTm :: (KnownNat m, KnownNat n) => L m n -> Sym n
-mTm x = Sym (tr x <> x)
-
-unSym :: Sym n -> Sq n
-unSym (Sym x) = x
-
-
-𝑖 :: Sized ℂ s c => s
-𝑖 = konst iC
-
-newtype Her n = Her (M n n)
-
-her :: KnownNat n => M n n -> Her n
-her m = Her $ (m + LA.tr m)/2
-
-
-
-instance KnownNat n => Eigen (Sym n) (R n) (L n n)
-  where
-    eigenvalues (Sym (extract -> m)) =  mkR . LA.eigenvaluesSH' $ m
-    eigensystem (Sym (extract -> m)) = (mkR l, mkL v)
-      where
-        (l,v) = LA.eigSH' m
-
-instance KnownNat n => Eigen (Sq n) (C n) (M n n)
-  where
-    eigenvalues (extract -> m) = mkC . LA.eigenvalues $ m
-    eigensystem (extract -> m) = (mkC l, mkM v)
-      where
-        (l,v) = LA.eig m
-
---------------------------------------------------------------------------------
-
-withNullspace
-    :: forall m n z . (KnownNat m, KnownNat n)
-    => L m n
-    -> (forall k . (KnownNat k) => L n k -> z)
-    -> z
-withNullspace (LA.nullspace . extract -> a) f =
-    case someNatVal $ fromIntegral $ cols a of
-       Nothing -> error "static/dynamic mismatch"
-       Just (SomeNat (_ :: Proxy k)) -> f (mkL a :: L n k)
-
-
-withCompactSVD
-    :: forall m n z . (KnownNat m, KnownNat n)
-    => L m n
-    -> (forall k . (KnownNat k) => (L m k, R k, L n k) -> z)
-    -> z
-withCompactSVD (LA.compactSVD . extract -> (u,s,v)) f =
-    case someNatVal $ fromIntegral $ size s of
-       Nothing -> error "static/dynamic mismatch"
-       Just (SomeNat (_ :: Proxy k)) -> f (mkL u :: L m k, mkR s :: R k, mkL v :: L n k)
-
---------------------------------------------------------------------------------
-
-qr :: (KnownNat m, KnownNat n) => L m n -> (L m m, L m n)
-qr (extract -> x) = (mkL q, mkL r)
-  where
-    (q,r) = LA.qr x
-
--- use qrRaw?
-
---------------------------------------------------------------------------------
-
-split :: forall p n . (KnownNat p, KnownNat n, p<=n) => R n -> (R p, R (n-p))
-split (extract -> v) = ( mkR (subVector 0 p' v) ,
-                         mkR (subVector p' (size v - p') v) )
-  where
-    p' = fromIntegral . natVal $ (undefined :: Proxy p) :: Int
-
-
-headTail :: (KnownNat n, 1<=n) => R n -> (ℝ, R (n-1))
-headTail = ((!0) . extract *** id) . split
-
-
-splitRows :: forall p m n . (KnownNat p, KnownNat m, KnownNat n, p<=m) => L m n -> (L p n, L (m-p) n)
-splitRows (extract -> x) = ( mkL (takeRows p' x) ,
-                             mkL (dropRows p' x) )
-  where
-    p' = fromIntegral . natVal $ (undefined :: Proxy p) :: Int
-
-splitCols :: forall p m n. (KnownNat p, KnownNat m, KnownNat n, KnownNat (n-p), p<=n) => L m n -> (L m p, L m (n-p))
-splitCols = (tr *** tr) . splitRows . tr
-
-
-toRows :: forall m n . (KnownNat m, KnownNat n) => L m n -> [R n]
-toRows (LA.toRows . extract -> vs) = map mkR vs
-
-
-toColumns :: forall m n . (KnownNat m, KnownNat n) => L m n -> [R m]
-toColumns (LA.toColumns . extract -> vs) = map mkR vs
-
-
---------------------------------------------------------------------------------
-
-build
-  :: forall m n . (KnownNat n, KnownNat m)
-    => (ℝ -> ℝ -> ℝ)
-    -> L m n
-build f = mkL $ LA.build (m',n') f
-  where
-    m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
-    n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
-
---------------------------------------------------------------------------------
-
-withVector
-    :: forall z
-     . Vector ℝ
-    -> (forall n . (KnownNat n) => R n -> z)
-    -> z
-withVector v f =
-    case someNatVal $ fromIntegral $ size v of
-       Nothing -> error "static/dynamic mismatch"
-       Just (SomeNat (_ :: Proxy m)) -> f (mkR v :: R m)
-
-
-withMatrix
-    :: forall z
-     . Matrix ℝ
-    -> (forall m n . (KnownNat m, KnownNat n) => L m n -> z)
-    -> z
-withMatrix a f =
-    case someNatVal $ fromIntegral $ rows a of
-       Nothing -> error "static/dynamic mismatch"
-       Just (SomeNat (_ :: Proxy m)) ->
-           case someNatVal $ fromIntegral $ cols a of
-               Nothing -> error "static/dynamic mismatch"
-               Just (SomeNat (_ :: Proxy n)) ->
-                  f (mkL a :: L m n)
-
---------------------------------------------------------------------------------
-
-class Domain field vec mat | mat -> vec field, vec -> mat field, field -> mat vec
-  where
-    mul :: forall m k n. (KnownNat m, KnownNat k, KnownNat n) => mat m k -> mat k n -> mat m n
-    app :: forall m n . (KnownNat m, KnownNat n) => mat m n -> vec n -> vec m
-    dot :: forall n . (KnownNat n) => vec n -> vec n -> field
-    cross :: vec 3 -> vec 3 -> vec 3
-    diagR ::  forall m n k . (KnownNat m, KnownNat n, KnownNat k) => field -> vec k -> mat m n
-
-
-instance Domain ℝ R L
-  where
-    mul = mulR
-    app = appR
-    dot = dotR
-    cross = crossR
-    diagR = diagRectR
-
-instance Domain ℂ C M
-  where
-    mul = mulC
-    app = appC
-    dot = dotC
-    cross = crossC
-    diagR = diagRectC
-
---------------------------------------------------------------------------------
-
-mulR :: forall m k n. (KnownNat m, KnownNat k, KnownNat n) => L m k -> L k n -> L m n
-
-mulR (isKonst -> Just (a,(_,k))) (isKonst -> Just (b,_)) = konst (a * b * fromIntegral k)
-
-mulR (isDiag -> Just (0,a,_)) (isDiag -> Just (0,b,_)) = diagR 0 (mkR v :: R k)
-  where
-    v = a' * b'
-    n = min (size a) (size b)
-    a' = subVector 0 n a
-    b' = subVector 0 n b
-
-mulR (isDiag -> Just (0,a,_)) (extract -> b) = mkL (asColumn a * takeRows (size a) b)
-
-mulR (extract -> a) (isDiag -> Just (0,b,_)) = mkL (takeColumns (size b) a * asRow b)
-
-mulR a b = mkL (extract a LA.<> extract b)
-
-
-appR :: (KnownNat m, KnownNat n) => L m n -> R n -> R m
-appR (isDiag -> Just (0, w, _)) v = mkR (w * subVector 0 (size w) (extract v))
-appR m v = mkR (extract m LA.#> extract v)
-
-
-dotR :: R n -> R n -> ℝ
-dotR (ud1 -> u) (ud1 -> v)
-    | singleV u || singleV v = sumElements (u * v)
-    | otherwise = udot u v
-
-
-crossR :: R 3 -> R 3 -> R 3
-crossR (extract -> x) (extract -> y) = vec3 z1 z2 z3
-  where
-    z1 = x!1*y!2-x!2*y!1
-    z2 = x!2*y!0-x!0*y!2
-    z3 = x!0*y!1-x!1*y!0
-
---------------------------------------------------------------------------------
-
-mulC :: forall m k n. (KnownNat m, KnownNat k, KnownNat n) => M m k -> M k n -> M m n
-
-mulC (isKonstC -> Just (a,(_,k))) (isKonstC -> Just (b,_)) = konst (a * b * fromIntegral k)
-
-mulC (isDiagC -> Just (0,a,_)) (isDiagC -> Just (0,b,_)) = diagR 0 (mkC v :: C k)
-  where
-    v = a' * b'
-    n = min (size a) (size b)
-    a' = subVector 0 n a
-    b' = subVector 0 n b
-
-mulC (isDiagC -> Just (0,a,_)) (extract -> b) = mkM (asColumn a * takeRows (size a) b)
-
-mulC (extract -> a) (isDiagC -> Just (0,b,_)) = mkM (takeColumns (size b) a * asRow b)
-
-mulC a b = mkM (extract a LA.<> extract b)
-
-
-appC :: (KnownNat m, KnownNat n) => M m n -> C n -> C m
-appC (isDiagC -> Just (0, w, _)) v = mkC (w * subVector 0 (size w) (extract v))
-appC m v = mkC (extract m LA.#> extract v)
-
-
-dotC :: KnownNat n => C n -> C n -> ℂ
-dotC (unwrap -> u) (unwrap -> v)
-    | singleV u || singleV v = sumElements (conj u * v)
-    | otherwise = u LA.<·> v
-
-
-crossC :: C 3 -> C 3 -> C 3
-crossC (extract -> x) (extract -> y) = mkC (LA.fromList [z1, z2, z3])
-  where
-    z1 = x!1*y!2-x!2*y!1
-    z2 = x!2*y!0-x!0*y!2
-    z3 = x!0*y!1-x!1*y!0
-
---------------------------------------------------------------------------------
-
-diagRectR :: forall m n k . (KnownNat m, KnownNat n, KnownNat k) => ℝ -> R k -> L m n
-diagRectR x v = mkL (asRow (vjoin [scalar x, ev, zeros]))
-  where
-    ev = extract v
-    zeros = LA.konst x (max 0 ((min m' n') - size ev))
-    m' = fromIntegral . natVal $ (undefined :: Proxy m)
-    n' = fromIntegral . natVal $ (undefined :: Proxy n)
-
-
-diagRectC :: forall m n k . (KnownNat m, KnownNat n, KnownNat k) => ℂ -> C k -> M m n
-diagRectC x v = mkM (asRow (vjoin [scalar x, ev, zeros]))
-  where
-    ev = extract v
-    zeros = LA.konst x (max 0 ((min m' n') - size ev))
-    m' = fromIntegral . natVal $ (undefined :: Proxy m)
-    n' = fromIntegral . natVal $ (undefined :: Proxy n)
-
---------------------------------------------------------------------------------
-
-mean :: (KnownNat n, 1<=n) => R n -> ℝ
-mean v = v <·> (1/dim)
-
-test :: (Bool, IO ())
-test = (ok,info)
-  where
-    ok =   extract (eye :: Sq 5) == ident 5
-           && (unwrap .unSym) (mTm sm :: Sym 3) == tr ((3><3)[1..]) LA.<> (3><3)[1..]
-           && unwrap (tm :: L 3 5) == LA.matrix 5 [1..15]
-           && thingS == thingD
-           && precS == precD
-           && withVector (LA.vector [1..15]) sumV == sumElements (LA.fromList [1..15])
-
-    info = do
-        print $ u
-        print $ v
-        print (eye :: Sq 3)
-        print $ ((u & 5) + 1) <·> v
-        print (tm :: L 2 5)
-        print (tm <> sm :: L 2 3)
-        print thingS
-        print thingD
-        print precS
-        print precD
-        print $ withVector (LA.vector [1..15]) sumV
-        splittest
-
-    sumV w = w <·> konst 1
-
-    u = vec2 3 5
-
-    𝕧 x = vector [x] :: R 1
-
-    v = 𝕧 2 & 4 & 7
-
-    tm :: GL
-    tm = lmat 0 [1..]
-
-    lmat :: forall m n . (KnownNat m, KnownNat n) => ℝ -> [ℝ] -> L m n
-    lmat z xs = mkL . reshape n' . LA.fromList . take (m'*n') $ xs ++ repeat z
-      where
-        m' = fromIntegral . natVal $ (undefined :: Proxy m)
-        n' = fromIntegral . natVal $ (undefined :: Proxy n)
-
-    sm :: GSq
-    sm = lmat 0 [1..]
-
-    thingS = (u & 1) <·> tr q #> q #> v
-      where
-        q = tm :: L 10 3
-
-    thingD = vjoin [ud1 u, 1] LA.<·> tr m LA.#> m LA.#> ud1 v
-      where
-        m = LA.matrix 3 [1..30]
-
-    precS = (1::Double) + (2::Double) * ((1 :: R 3) * (u & 6)) <·> konst 2 #> v
-    precD = 1 + 2 * vjoin[ud1 u, 6] LA.<·> LA.konst 2 (size (ud1 u) +1, size (ud1 v)) LA.#> ud1 v
-
-
-splittest
-    = do
-    let v = range :: R 7
-        a = snd (split v) :: R 4
-    print $ a
-    print $ snd . headTail . snd . headTail $ v
-    print $ first (vec3 1 2 3)
-    print $ second (vec3 1 2 3)
-    print $ third (vec3 1 2 3)
-    print $ (snd $ splitRows eye :: L 4 6)
- where
-    first v = fst . headTail $ v
-    second v = first . snd . headTail $ v
-    third v = first . snd . headTail . snd . headTail $ v
-
-
-instance (KnownNat n', KnownNat m') => Testable (L n' m')
-  where
-    checkT _ = test
-
-
diff --git a/packages/base/src/Numeric/LinearAlgebra/Static.hs b/packages/base/src/Numeric/LinearAlgebra/Static.hs
deleted file mode 100644
index daf8d80..0000000
--- a/packages/base/src/Numeric/LinearAlgebra/Static.hs
+++ /dev/null
@@ -1,422 +0,0 @@
-{-# LANGUAGE DataKinds #-}
-{-# LANGUAGE KindSignatures #-}
-{-# LANGUAGE GeneralizedNewtypeDeriving #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE FunctionalDependencies #-}
-{-# LANGUAGE FlexibleContexts #-}
-{-# LANGUAGE ScopedTypeVariables #-}
-{-# LANGUAGE EmptyDataDecls #-}
-{-# LANGUAGE Rank2Types #-}
-{-# LANGUAGE FlexibleInstances #-}
-{-# LANGUAGE TypeOperators #-}
-{-# LANGUAGE ViewPatterns #-}
-{-# LANGUAGE GADTs #-}
-
-
-{- |
-Module      :  Numeric.LinearAlgebra.Static
-Copyright   :  (c) Alberto Ruiz 2006-14
-License     :  BSD3
-Stability   :  provisional
-
--}
-
-module Numeric.LinearAlgebra.Static where
-
-
-import GHC.TypeLits
-import Numeric.LinearAlgebra.HMatrix as LA
-import Data.Packed as D
-import Data.Packed.ST
-import Data.Proxy(Proxy)
-import Foreign.Storable(Storable)
-import Text.Printf
-
---------------------------------------------------------------------------------
-
-newtype Dim (n :: Nat) t = Dim t
-  deriving Show
-
-lift1F
-  :: (c t -> c t)
-  -> Dim n (c t) -> Dim n (c t)
-lift1F f (Dim v) = Dim (f v)
-
-lift2F
-  :: (c t -> c t -> c t)
-  -> Dim n (c t) -> Dim n (c t) -> Dim n (c t)
-lift2F f (Dim u) (Dim v) = Dim (f u v)
-
---------------------------------------------------------------------------------
-
-newtype R n = R (Dim n (Vector ℝ))
-  deriving (Num,Fractional)
-
-newtype C n = C (Dim n (Vector ℂ))
-  deriving (Num,Fractional)
-
-newtype L m n = L (Dim m (Dim n (Matrix ℝ)))
-
-newtype M m n = M (Dim m (Dim n (Matrix  ℂ)))
-
-
-mkR :: Vector ℝ -> R n
-mkR = R . Dim
-
-mkC :: Vector ℂ -> C n
-mkC = C . Dim
-
-mkL :: Matrix ℝ -> L m n
-mkL x = L (Dim (Dim x))
-
-mkM :: Matrix ℂ -> M m n
-mkM x = M (Dim (Dim x))
-
---------------------------------------------------------------------------------
-
-type V n t = Dim n (Vector t)
-
-ud :: Dim n (Vector t) -> Vector t
-ud (Dim v) = v
-
-mkV :: forall (n :: Nat) t . t -> Dim n t
-mkV = Dim 
-
-
-vconcat :: forall n m t . (KnownNat n, KnownNat m, Numeric t)
-    => V n t -> V m t -> V (n+m) t
-(ud -> u) `vconcat` (ud -> v) = mkV (vjoin [u', v'])
-  where
-    du = fromIntegral . natVal $ (undefined :: Proxy n)
-    dv = fromIntegral . natVal $ (undefined :: Proxy m)
-    u' | du > 1 && size u == 1 = LA.konst (u D.@> 0) du
-       | otherwise = u
-    v' | dv > 1 && size v == 1 = LA.konst (v D.@> 0) dv
-       | otherwise = v
-
-
-gvec2 :: Storable t => t -> t -> V 2 t
-gvec2 a b = mkV $ runSTVector $ do
-    v <- newUndefinedVector 2
-    writeVector v 0 a
-    writeVector v 1 b
-    return v
-
-gvec3 :: Storable t => t -> t -> t -> V 3 t
-gvec3 a b c = mkV $ runSTVector $ do
-    v <- newUndefinedVector 3
-    writeVector v 0 a
-    writeVector v 1 b
-    writeVector v 2 c
-    return v
-
-
-gvec4 :: Storable t => t -> t -> t -> t -> V 4 t
-gvec4 a b c d = mkV $ runSTVector $ do
-    v <- newUndefinedVector 4
-    writeVector v 0 a
-    writeVector v 1 b
-    writeVector v 2 c
-    writeVector v 3 d
-    return v
-
-
-gvect :: forall n t . (Show t, KnownNat n, Numeric t) => String -> [t] -> V n t
-gvect st xs'
-    | ok = mkV v
-    | not (null rest) && null (tail rest) = abort (show xs')
-    | not (null rest) = abort (init (show (xs++take 1 rest))++", ... ]")
-    | otherwise = abort (show xs)
-  where
-    (xs,rest) = splitAt d xs'
-    ok = size v == d && null rest
-    v = LA.fromList xs
-    d = fromIntegral . natVal $ (undefined :: Proxy n)
-    abort info = error $ st++" "++show d++" can't be created from elements "++info
-
-
---------------------------------------------------------------------------------
-
-type GM m n t = Dim m (Dim n (Matrix t))
-
-
-gmat :: forall m n t . (Show t, KnownNat m, KnownNat n, Numeric t) => String -> [t] -> GM m n t
-gmat st xs'
-    | ok = Dim (Dim x)
-    | not (null rest) && null (tail rest) = abort (show xs')
-    | not (null rest) = abort (init (show (xs++take 1 rest))++", ... ]")
-    | otherwise = abort (show xs)
-  where
-    (xs,rest) = splitAt (m'*n') xs'
-    v = LA.fromList xs
-    x = reshape n' v
-    ok = rem (size v) n' == 0 && size x == (m',n') && null rest
-    m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
-    n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
-    abort info = error $ st ++" "++show m' ++ " " ++ show n'++" can't be created from elements " ++ info
-
---------------------------------------------------------------------------------
-
-class Num t => Sized t s d | s -> t, s -> d
-  where
-    konst     ::  t  -> s
-    unwrap    ::  s  -> d
-    fromList  :: [t] -> s
-    extract   ::  s  -> d
-
-singleV v = size v == 1
-singleM m = rows m == 1 && cols m == 1
-
-
-instance forall n. KnownNat n => Sized ℂ (C n) (Vector ℂ)
-  where
-    konst x = mkC (LA.scalar x)
-    unwrap (C (Dim v)) = v
-    fromList xs = C (gvect "C" xs)
-    extract (unwrap -> v)
-      | singleV v = LA.konst (v!0) d
-      | otherwise = v
-     where
-       d = fromIntegral . natVal $ (undefined :: Proxy n)
-
-
-instance forall n. KnownNat n => Sized ℝ (R n) (Vector ℝ)
-  where
-    konst x = mkR (LA.scalar x)
-    unwrap (R (Dim v)) = v
-    fromList xs = R (gvect "R" xs)
-    extract (unwrap -> v)
-      | singleV v = LA.konst (v!0) d
-      | otherwise = v
-     where
-       d = fromIntegral . natVal $ (undefined :: Proxy n)
-
-
-
-instance forall m n . (KnownNat m, KnownNat n) => Sized ℝ (L m n) (Matrix ℝ)
-  where
-    konst x = mkL (LA.scalar x)
-    fromList xs = L (gmat "L" xs)
-    unwrap (L (Dim (Dim m))) = m
-    extract (isDiag -> Just (z,y,(m',n'))) = diagRect z y m' n'
-    extract (unwrap -> a)
-        | singleM a = LA.konst (a `atIndex` (0,0)) (m',n')
-        | otherwise = a
-      where
-        m' = fromIntegral . natVal $ (undefined :: Proxy m)
-        n' = fromIntegral . natVal $ (undefined :: Proxy n)
-
-
-instance forall m n . (KnownNat m, KnownNat n) => Sized ℂ (M m n) (Matrix ℂ)
-  where
-    konst x = mkM (LA.scalar x)
-    fromList xs = M (gmat "M" xs)
-    unwrap (M (Dim (Dim m))) = m
-    extract (isDiagC -> Just (z,y,(m',n'))) = diagRect z y m' n'
-    extract (unwrap -> a)
-        | singleM a = LA.konst (a `atIndex` (0,0)) (m',n')
-        | otherwise = a
-      where
-        m' = fromIntegral . natVal $ (undefined :: Proxy m)
-        n' = fromIntegral . natVal $ (undefined :: Proxy n)
-
---------------------------------------------------------------------------------
-
-instance (KnownNat n, KnownNat m) => Transposable (L m n) (L n m)
-  where
-    tr a@(isDiag -> Just _) = mkL (extract a)
-    tr (extract -> a) = mkL (tr a)
-
-instance (KnownNat n, KnownNat m) => Transposable (M m n) (M n m)
-  where
-    tr a@(isDiagC -> Just _) = mkM (extract a)
-    tr (extract -> a) = mkM (tr a)
-
---------------------------------------------------------------------------------
-
-isDiag :: forall m n . (KnownNat m, KnownNat n) => L m n -> Maybe (ℝ, Vector ℝ, (Int,Int))
-isDiag (L x) = isDiagg x
-
-isDiagC :: forall m n . (KnownNat m, KnownNat n) => M m n -> Maybe (ℂ, Vector ℂ, (Int,Int))
-isDiagC (M x) = isDiagg x
-
-
-isDiagg :: forall m n t . (Numeric t, KnownNat m, KnownNat n) => GM m n t -> Maybe (t, Vector t, (Int,Int))
-isDiagg (Dim (Dim x))
-    | singleM x = Nothing
-    | rows x == 1 && m' > 1 || cols x == 1 && n' > 1 = Just (z,yz,(m',n'))
-    | otherwise = Nothing
-  where
-    m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
-    n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
-    v = flatten x
-    z = v `atIndex` 0
-    y = subVector 1 (size v-1) v
-    ny = size y
-    zeros = LA.konst 0 (max 0 (min m' n' - ny))
-    yz = vjoin [y,zeros]
-
---------------------------------------------------------------------------------
-
-instance forall n . KnownNat n => Show (R n)
-  where
-    show (R (Dim v))
-      | singleV v = "("++show (v!0)++" :: R "++show d++")"
-      | otherwise   = "(vector"++ drop 8 (show v)++" :: R "++show d++")"
-      where
-        d = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
-
-instance forall n . KnownNat n => Show (C n)
-  where
-    show (C (Dim v))
-      | singleV v = "("++show (v!0)++" :: C "++show d++")"
-      | otherwise   = "(vector"++ drop 8 (show v)++" :: C "++show d++")"
-      where
-        d = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
-
-instance forall m n . (KnownNat m, KnownNat n) => Show (L m n)
-  where
-    show (isDiag -> Just (z,y,(m',n'))) = printf "(diag %s %s :: L %d %d)" (show z) (drop 9 $ show y) m' n'
-    show (L (Dim (Dim x)))
-       | singleM x = printf "(%s :: L %d %d)" (show (x `atIndex` (0,0))) m' n'
-       | otherwise = "(matrix"++ dropWhile (/='\n') (show x)++" :: L "++show m'++" "++show n'++")"
-      where
-        m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
-        n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
-
-instance forall m n . (KnownNat m, KnownNat n) => Show (M m n)
-  where
-    show (isDiagC -> Just (z,y,(m',n'))) = printf "(diag %s %s :: M %d %d)" (show z) (drop 9 $ show y) m' n'
-    show (M (Dim (Dim x)))
-       | singleM x = printf "(%s :: M %d %d)" (show (x `atIndex` (0,0))) m' n'
-       | otherwise = "(matrix"++ dropWhile (/='\n') (show x)++" :: M "++show m'++" "++show n'++")"
-      where
-        m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
-        n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
-
---------------------------------------------------------------------------------
-
-instance forall n t . (Num (Vector t), Numeric t )=> Num (Dim n (Vector t))
-  where
-    (+) = lift2F (+)
-    (*) = lift2F (*)
-    (-) = lift2F (-)
-    abs = lift1F abs
-    signum = lift1F signum
-    negate = lift1F negate
-    fromInteger x = Dim (fromInteger x)
-
-instance (Num (Vector t), Num (Matrix t), Numeric t) => Fractional (Dim n (Vector t))
-  where
-    fromRational x = Dim (fromRational x)
-    (/) = lift2F (/)
-
-
-instance (Num (Matrix t), Numeric t) => Num (Dim m (Dim n (Matrix t)))
-  where
-    (+) = (lift2F . lift2F) (+)
-    (*) = (lift2F . lift2F) (*)
-    (-) = (lift2F . lift2F) (-)
-    abs = (lift1F . lift1F) abs
-    signum = (lift1F . lift1F) signum
-    negate = (lift1F . lift1F) negate
-    fromInteger x = Dim (Dim (fromInteger x))
-
-instance (Num (Vector t), Num (Matrix t), Numeric t) => Fractional (Dim m (Dim n (Matrix t)))
-  where
-    fromRational x = Dim (Dim (fromRational x))
-    (/) = (lift2F.lift2F) (/)
-
---------------------------------------------------------------------------------
-
-
-adaptDiag f a@(isDiag -> Just _) b | isFull b = f (mkL (extract a)) b
-adaptDiag f a b@(isDiag -> Just _) | isFull a = f a (mkL (extract b))
-adaptDiag f a b = f a b
-
-isFull m = isDiag m == Nothing && not (singleM (unwrap m))
-
-
-lift1L f (L v) = L (f v)
-lift2L f (L a) (L b) = L (f a b)
-lift2LD f = adaptDiag (lift2L f)
-
-
-instance (KnownNat n, KnownNat m) =>  Num (L n m)
-  where
-    (+) = lift2LD (+)
-    (*) = lift2LD (*)
-    (-) = lift2LD (-)
-    abs = lift1L abs
-    signum = lift1L signum
-    negate = lift1L negate
-    fromInteger = L . Dim . Dim . fromInteger
-
-instance (KnownNat n, KnownNat m) => Fractional (L n m)
-  where
-    fromRational = L . Dim . Dim . fromRational
-    (/) = lift2LD (/)
-
---------------------------------------------------------------------------------
-
-adaptDiagC f a@(isDiagC -> Just _) b | isFullC b = f (mkM (extract a)) b
-adaptDiagC f a b@(isDiagC -> Just _) | isFullC a = f a (mkM (extract b))
-adaptDiagC f a b = f a b
-
-isFullC m = isDiagC m == Nothing && not (singleM (unwrap m))
-
-lift1M f (M v) = M (f v)
-lift2M f (M a) (M b) = M (f a b)
-lift2MD f = adaptDiagC (lift2M f)
-
-instance (KnownNat n, KnownNat m) =>  Num (M n m)
-  where
-    (+) = lift2MD (+)
-    (*) = lift2MD (*)
-    (-) = lift2MD (-)
-    abs = lift1M abs
-    signum = lift1M signum
-    negate = lift1M negate
-    fromInteger = M . Dim . Dim . fromInteger
-
-instance (KnownNat n, KnownNat m) => Fractional (M n m)
-  where
-    fromRational = M . Dim . Dim . fromRational
-    (/) = lift2MD (/)
-
---------------------------------------------------------------------------------
-
-
-class Disp t
-  where
-    disp :: Int -> t -> IO ()
-
-
-instance (KnownNat m, KnownNat n) => Disp (L m n)
-  where
-    disp n x = do
-        let a = extract x
-        let su = LA.dispf n a
-        printf "L %d %d" (rows a) (cols a) >> putStr (dropWhile (/='\n') $ su)
-
-instance (KnownNat m, KnownNat n) => Disp (M m n)
-  where
-    disp n x = do
-        let a = extract x
-        let su = LA.dispcf n a
-        printf "M %d %d" (rows a) (cols a) >> putStr (dropWhile (/='\n') $ su)
-
-
-instance KnownNat n => Disp (R n)
-  where
-    disp n v = do
-        let su = LA.dispf n (asRow $ extract v)
-        putStr "R " >> putStr (tail . dropWhile (/='x') $ su)
-
-instance KnownNat n => Disp (C n)
-  where
-    disp n v = do
-        let su = LA.dispcf n (asRow $ extract v)
-        putStr "C " >> putStr (tail . dropWhile (/='x') $ su)
-
---------------------------------------------------------------------------------
diff --git a/packages/base/src/Numeric/LinearAlgebra/Static/Internal.hs b/packages/base/src/Numeric/LinearAlgebra/Static/Internal.hs
new file mode 100644
index 0000000..c9641d5
--- /dev/null
+++ b/packages/base/src/Numeric/LinearAlgebra/Static/Internal.hs
@@ -0,0 +1,422 @@
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE KindSignatures #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FunctionalDependencies #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE EmptyDataDecls #-}
+{-# LANGUAGE Rank2Types #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE ViewPatterns #-}
+{-# LANGUAGE GADTs #-}
+
+
+{- |
+Module      :  Numeric.LinearAlgebra.Static.Internal
+Copyright   :  (c) Alberto Ruiz 2006-14
+License     :  BSD3
+Stability   :  provisional
+
+-}
+
+module Numeric.LinearAlgebra.Static.Internal where
+
+
+import GHC.TypeLits
+import Numeric.LinearAlgebra.HMatrix as LA
+import Data.Packed as D
+import Data.Packed.ST
+import Data.Proxy(Proxy)
+import Foreign.Storable(Storable)
+import Text.Printf
+
+--------------------------------------------------------------------------------
+
+newtype Dim (n :: Nat) t = Dim t
+  deriving Show
+
+lift1F
+  :: (c t -> c t)
+  -> Dim n (c t) -> Dim n (c t)
+lift1F f (Dim v) = Dim (f v)
+
+lift2F
+  :: (c t -> c t -> c t)
+  -> Dim n (c t) -> Dim n (c t) -> Dim n (c t)
+lift2F f (Dim u) (Dim v) = Dim (f u v)
+
+--------------------------------------------------------------------------------
+
+newtype R n = R (Dim n (Vector ℝ))
+  deriving (Num,Fractional)
+
+newtype C n = C (Dim n (Vector ℂ))
+  deriving (Num,Fractional)
+
+newtype L m n = L (Dim m (Dim n (Matrix ℝ)))
+
+newtype M m n = M (Dim m (Dim n (Matrix  ℂ)))
+
+
+mkR :: Vector ℝ -> R n
+mkR = R . Dim
+
+mkC :: Vector ℂ -> C n
+mkC = C . Dim
+
+mkL :: Matrix ℝ -> L m n
+mkL x = L (Dim (Dim x))
+
+mkM :: Matrix ℂ -> M m n
+mkM x = M (Dim (Dim x))
+
+--------------------------------------------------------------------------------
+
+type V n t = Dim n (Vector t)
+
+ud :: Dim n (Vector t) -> Vector t
+ud (Dim v) = v
+
+mkV :: forall (n :: Nat) t . t -> Dim n t
+mkV = Dim 
+
+
+vconcat :: forall n m t . (KnownNat n, KnownNat m, Numeric t)
+    => V n t -> V m t -> V (n+m) t
+(ud -> u) `vconcat` (ud -> v) = mkV (vjoin [u', v'])
+  where
+    du = fromIntegral . natVal $ (undefined :: Proxy n)
+    dv = fromIntegral . natVal $ (undefined :: Proxy m)
+    u' | du > 1 && size u == 1 = LA.konst (u D.@> 0) du
+       | otherwise = u
+    v' | dv > 1 && size v == 1 = LA.konst (v D.@> 0) dv
+       | otherwise = v
+
+
+gvec2 :: Storable t => t -> t -> V 2 t
+gvec2 a b = mkV $ runSTVector $ do
+    v <- newUndefinedVector 2
+    writeVector v 0 a
+    writeVector v 1 b
+    return v
+
+gvec3 :: Storable t => t -> t -> t -> V 3 t
+gvec3 a b c = mkV $ runSTVector $ do
+    v <- newUndefinedVector 3
+    writeVector v 0 a
+    writeVector v 1 b
+    writeVector v 2 c
+    return v
+
+
+gvec4 :: Storable t => t -> t -> t -> t -> V 4 t
+gvec4 a b c d = mkV $ runSTVector $ do
+    v <- newUndefinedVector 4
+    writeVector v 0 a
+    writeVector v 1 b
+    writeVector v 2 c
+    writeVector v 3 d
+    return v
+
+
+gvect :: forall n t . (Show t, KnownNat n, Numeric t) => String -> [t] -> V n t
+gvect st xs'
+    | ok = mkV v
+    | not (null rest) && null (tail rest) = abort (show xs')
+    | not (null rest) = abort (init (show (xs++take 1 rest))++", ... ]")
+    | otherwise = abort (show xs)
+  where
+    (xs,rest) = splitAt d xs'
+    ok = size v == d && null rest
+    v = LA.fromList xs
+    d = fromIntegral . natVal $ (undefined :: Proxy n)
+    abort info = error $ st++" "++show d++" can't be created from elements "++info
+
+
+--------------------------------------------------------------------------------
+
+type GM m n t = Dim m (Dim n (Matrix t))
+
+
+gmat :: forall m n t . (Show t, KnownNat m, KnownNat n, Numeric t) => String -> [t] -> GM m n t
+gmat st xs'
+    | ok = Dim (Dim x)
+    | not (null rest) && null (tail rest) = abort (show xs')
+    | not (null rest) = abort (init (show (xs++take 1 rest))++", ... ]")
+    | otherwise = abort (show xs)
+  where
+    (xs,rest) = splitAt (m'*n') xs'
+    v = LA.fromList xs
+    x = reshape n' v
+    ok = rem (size v) n' == 0 && size x == (m',n') && null rest
+    m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
+    n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
+    abort info = error $ st ++" "++show m' ++ " " ++ show n'++" can't be created from elements " ++ info
+
+--------------------------------------------------------------------------------
+
+class Num t => Sized t s d | s -> t, s -> d
+  where
+    konst     ::  t  -> s
+    unwrap    ::  s  -> d
+    fromList  :: [t] -> s
+    extract   ::  s  -> d
+
+singleV v = size v == 1
+singleM m = rows m == 1 && cols m == 1
+
+
+instance forall n. KnownNat n => Sized ℂ (C n) (Vector ℂ)
+  where
+    konst x = mkC (LA.scalar x)
+    unwrap (C (Dim v)) = v
+    fromList xs = C (gvect "C" xs)
+    extract (unwrap -> v)
+      | singleV v = LA.konst (v!0) d
+      | otherwise = v
+     where
+       d = fromIntegral . natVal $ (undefined :: Proxy n)
+
+
+instance forall n. KnownNat n => Sized ℝ (R n) (Vector ℝ)
+  where
+    konst x = mkR (LA.scalar x)
+    unwrap (R (Dim v)) = v
+    fromList xs = R (gvect "R" xs)
+    extract (unwrap -> v)
+      | singleV v = LA.konst (v!0) d
+      | otherwise = v
+     where
+       d = fromIntegral . natVal $ (undefined :: Proxy n)
+
+
+
+instance forall m n . (KnownNat m, KnownNat n) => Sized ℝ (L m n) (Matrix ℝ)
+  where
+    konst x = mkL (LA.scalar x)
+    fromList xs = L (gmat "L" xs)
+    unwrap (L (Dim (Dim m))) = m
+    extract (isDiag -> Just (z,y,(m',n'))) = diagRect z y m' n'
+    extract (unwrap -> a)
+        | singleM a = LA.konst (a `atIndex` (0,0)) (m',n')
+        | otherwise = a
+      where
+        m' = fromIntegral . natVal $ (undefined :: Proxy m)
+        n' = fromIntegral . natVal $ (undefined :: Proxy n)
+
+
+instance forall m n . (KnownNat m, KnownNat n) => Sized ℂ (M m n) (Matrix ℂ)
+  where
+    konst x = mkM (LA.scalar x)
+    fromList xs = M (gmat "M" xs)
+    unwrap (M (Dim (Dim m))) = m
+    extract (isDiagC -> Just (z,y,(m',n'))) = diagRect z y m' n'
+    extract (unwrap -> a)
+        | singleM a = LA.konst (a `atIndex` (0,0)) (m',n')
+        | otherwise = a
+      where
+        m' = fromIntegral . natVal $ (undefined :: Proxy m)
+        n' = fromIntegral . natVal $ (undefined :: Proxy n)
+
+--------------------------------------------------------------------------------
+
+instance (KnownNat n, KnownNat m) => Transposable (L m n) (L n m)
+  where
+    tr a@(isDiag -> Just _) = mkL (extract a)
+    tr (extract -> a) = mkL (tr a)
+
+instance (KnownNat n, KnownNat m) => Transposable (M m n) (M n m)
+  where
+    tr a@(isDiagC -> Just _) = mkM (extract a)
+    tr (extract -> a) = mkM (tr a)
+
+--------------------------------------------------------------------------------
+
+isDiag :: forall m n . (KnownNat m, KnownNat n) => L m n -> Maybe (ℝ, Vector ℝ, (Int,Int))
+isDiag (L x) = isDiagg x
+
+isDiagC :: forall m n . (KnownNat m, KnownNat n) => M m n -> Maybe (ℂ, Vector ℂ, (Int,Int))
+isDiagC (M x) = isDiagg x
+
+
+isDiagg :: forall m n t . (Numeric t, KnownNat m, KnownNat n) => GM m n t -> Maybe (t, Vector t, (Int,Int))
+isDiagg (Dim (Dim x))
+    | singleM x = Nothing
+    | rows x == 1 && m' > 1 || cols x == 1 && n' > 1 = Just (z,yz,(m',n'))
+    | otherwise = Nothing
+  where
+    m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
+    n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
+    v = flatten x
+    z = v `atIndex` 0
+    y = subVector 1 (size v-1) v
+    ny = size y
+    zeros = LA.konst 0 (max 0 (min m' n' - ny))
+    yz = vjoin [y,zeros]
+
+--------------------------------------------------------------------------------
+
+instance forall n . KnownNat n => Show (R n)
+  where
+    show (R (Dim v))
+      | singleV v = "("++show (v!0)++" :: R "++show d++")"
+      | otherwise   = "(vector"++ drop 8 (show v)++" :: R "++show d++")"
+      where
+        d = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
+
+instance forall n . KnownNat n => Show (C n)
+  where
+    show (C (Dim v))
+      | singleV v = "("++show (v!0)++" :: C "++show d++")"
+      | otherwise   = "(vector"++ drop 8 (show v)++" :: C "++show d++")"
+      where
+        d = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
+
+instance forall m n . (KnownNat m, KnownNat n) => Show (L m n)
+  where
+    show (isDiag -> Just (z,y,(m',n'))) = printf "(diag %s %s :: L %d %d)" (show z) (drop 9 $ show y) m' n'
+    show (L (Dim (Dim x)))
+       | singleM x = printf "(%s :: L %d %d)" (show (x `atIndex` (0,0))) m' n'
+       | otherwise = "(matrix"++ dropWhile (/='\n') (show x)++" :: L "++show m'++" "++show n'++")"
+      where
+        m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
+        n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
+
+instance forall m n . (KnownNat m, KnownNat n) => Show (M m n)
+  where
+    show (isDiagC -> Just (z,y,(m',n'))) = printf "(diag %s %s :: M %d %d)" (show z) (drop 9 $ show y) m' n'
+    show (M (Dim (Dim x)))
+       | singleM x = printf "(%s :: M %d %d)" (show (x `atIndex` (0,0))) m' n'
+       | otherwise = "(matrix"++ dropWhile (/='\n') (show x)++" :: M "++show m'++" "++show n'++")"
+      where
+        m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
+        n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
+
+--------------------------------------------------------------------------------
+
+instance forall n t . (Num (Vector t), Numeric t )=> Num (Dim n (Vector t))
+  where
+    (+) = lift2F (+)
+    (*) = lift2F (*)
+    (-) = lift2F (-)
+    abs = lift1F abs
+    signum = lift1F signum
+    negate = lift1F negate
+    fromInteger x = Dim (fromInteger x)
+
+instance (Num (Vector t), Num (Matrix t), Numeric t) => Fractional (Dim n (Vector t))
+  where
+    fromRational x = Dim (fromRational x)
+    (/) = lift2F (/)
+
+
+instance (Num (Matrix t), Numeric t) => Num (Dim m (Dim n (Matrix t)))
+  where
+    (+) = (lift2F . lift2F) (+)
+    (*) = (lift2F . lift2F) (*)
+    (-) = (lift2F . lift2F) (-)
+    abs = (lift1F . lift1F) abs
+    signum = (lift1F . lift1F) signum
+    negate = (lift1F . lift1F) negate
+    fromInteger x = Dim (Dim (fromInteger x))
+
+instance (Num (Vector t), Num (Matrix t), Numeric t) => Fractional (Dim m (Dim n (Matrix t)))
+  where
+    fromRational x = Dim (Dim (fromRational x))
+    (/) = (lift2F.lift2F) (/)
+
+--------------------------------------------------------------------------------
+
+
+adaptDiag f a@(isDiag -> Just _) b | isFull b = f (mkL (extract a)) b
+adaptDiag f a b@(isDiag -> Just _) | isFull a = f a (mkL (extract b))
+adaptDiag f a b = f a b
+
+isFull m = isDiag m == Nothing && not (singleM (unwrap m))
+
+
+lift1L f (L v) = L (f v)
+lift2L f (L a) (L b) = L (f a b)
+lift2LD f = adaptDiag (lift2L f)
+
+
+instance (KnownNat n, KnownNat m) =>  Num (L n m)
+  where
+    (+) = lift2LD (+)
+    (*) = lift2LD (*)
+    (-) = lift2LD (-)
+    abs = lift1L abs
+    signum = lift1L signum
+    negate = lift1L negate
+    fromInteger = L . Dim . Dim . fromInteger
+
+instance (KnownNat n, KnownNat m) => Fractional (L n m)
+  where
+    fromRational = L . Dim . Dim . fromRational
+    (/) = lift2LD (/)
+
+--------------------------------------------------------------------------------
+
+adaptDiagC f a@(isDiagC -> Just _) b | isFullC b = f (mkM (extract a)) b
+adaptDiagC f a b@(isDiagC -> Just _) | isFullC a = f a (mkM (extract b))
+adaptDiagC f a b = f a b
+
+isFullC m = isDiagC m == Nothing && not (singleM (unwrap m))
+
+lift1M f (M v) = M (f v)
+lift2M f (M a) (M b) = M (f a b)
+lift2MD f = adaptDiagC (lift2M f)
+
+instance (KnownNat n, KnownNat m) =>  Num (M n m)
+  where
+    (+) = lift2MD (+)
+    (*) = lift2MD (*)
+    (-) = lift2MD (-)
+    abs = lift1M abs
+    signum = lift1M signum
+    negate = lift1M negate
+    fromInteger = M . Dim . Dim . fromInteger
+
+instance (KnownNat n, KnownNat m) => Fractional (M n m)
+  where
+    fromRational = M . Dim . Dim . fromRational
+    (/) = lift2MD (/)
+
+--------------------------------------------------------------------------------
+
+
+class Disp t
+  where
+    disp :: Int -> t -> IO ()
+
+
+instance (KnownNat m, KnownNat n) => Disp (L m n)
+  where
+    disp n x = do
+        let a = extract x
+        let su = LA.dispf n a
+        printf "L %d %d" (rows a) (cols a) >> putStr (dropWhile (/='\n') $ su)
+
+instance (KnownNat m, KnownNat n) => Disp (M m n)
+  where
+    disp n x = do
+        let a = extract x
+        let su = LA.dispcf n a
+        printf "M %d %d" (rows a) (cols a) >> putStr (dropWhile (/='\n') $ su)
+
+
+instance KnownNat n => Disp (R n)
+  where
+    disp n v = do
+        let su = LA.dispf n (asRow $ extract v)
+        putStr "R " >> putStr (tail . dropWhile (/='x') $ su)
+
+instance KnownNat n => Disp (C n)
+  where
+    disp n v = do
+        let su = LA.dispcf n (asRow $ extract v)
+        putStr "C " >> putStr (tail . dropWhile (/='x') $ su)
+
+--------------------------------------------------------------------------------
diff --git a/packages/base/src/Numeric/LinearAlgebra/tmpStatic.hs b/packages/base/src/Numeric/LinearAlgebra/tmpStatic.hs
new file mode 100644
index 0000000..4258d6b
--- /dev/null
+++ b/packages/base/src/Numeric/LinearAlgebra/tmpStatic.hs
@@ -0,0 +1,619 @@
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE KindSignatures #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FunctionalDependencies #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE EmptyDataDecls #-}
+{-# LANGUAGE Rank2Types #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE ViewPatterns #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE OverlappingInstances #-}
+{-# LANGUAGE TypeFamilies #-}
+
+
+{- |
+Module      :  Numeric.LinearAlgebra.Static
+Copyright   :  (c) Alberto Ruiz 2014
+License     :  BSD3
+Stability   :  experimental
+
+Experimental interface with statically checked dimensions.
+
+-}
+
+module Numeric.LinearAlgebra.Static(
+    -- * Vector
+    ℝ, R,
+    vec2, vec3, vec4, (&), (#), split, headTail,
+    vector,
+    linspace, range, dim,
+    -- * Matrix
+    L, Sq, build,
+    row, col, (¦),(——), splitRows, splitCols,
+    unrow, uncol,
+    tr,
+    eye,
+    diag,
+    blockAt,
+    matrix,
+    -- * Complex
+    C, M, Her, her, 𝑖,
+    -- * Products
+    (<>),(#>),(<·>),
+    -- * Linear Systems
+    linSolve, (<\>),
+    -- * Factorizations
+    svd, withCompactSVD, svdTall, svdFlat, Eigen(..),
+    withNullspace, qr,
+    -- * Misc
+    mean,
+    Disp(..), Domain(..),
+    withVector, withMatrix,
+    toRows, toColumns,
+    Sized(..), Diag(..), Sym, sym, mTm, unSym
+) where
+
+
+import GHC.TypeLits
+import Numeric.LinearAlgebra.HMatrix hiding (
+    (<>),(#>),(<·>),Konst(..),diag, disp,(¦),(——),row,col,vector,matrix,linspace,toRows,toColumns,
+    (<\>),fromList,takeDiag,svd,eig,eigSH,eigSH',eigenvalues,eigenvaluesSH,eigenvaluesSH',build,
+    qr)
+import qualified Numeric.LinearAlgebra.HMatrix as LA
+import Data.Proxy(Proxy)
+import Numeric.LinearAlgebra.Static.Internal
+import Control.Arrow((***))
+
+
+
+
+
+ud1 :: R n -> Vector ℝ
+ud1 (R (Dim v)) = v
+
+
+infixl 4 &
+(&) :: forall n . (KnownNat n, 1 <= n)
+    => R n -> ℝ -> R (n+1)
+u & x = u # (konst x :: R 1)
+
+infixl 4 #
+(#) :: forall n m . (KnownNat n, KnownNat m)
+    => R n -> R m -> R (n+m)
+(R u) # (R v) = R (vconcat u v)
+
+
+
+vec2 :: ℝ -> ℝ -> R 2
+vec2 a b = R (gvec2 a b)
+
+vec3 :: ℝ -> ℝ -> ℝ -> R 3
+vec3 a b c = R (gvec3 a b c)
+
+
+vec4 :: ℝ -> ℝ -> ℝ -> ℝ -> R 4
+vec4 a b c d = R (gvec4 a b c d)
+
+vector :: KnownNat n => [ℝ] -> R n
+vector = fromList
+
+matrix :: (KnownNat m, KnownNat n) => [ℝ] -> L m n
+matrix = fromList
+
+linspace :: forall n . KnownNat n => (ℝ,ℝ) -> R n
+linspace (a,b) = mkR (LA.linspace d (a,b))
+  where
+    d = fromIntegral . natVal $ (undefined :: Proxy n)
+
+range :: forall n . KnownNat n => R n
+range = mkR (LA.linspace d (1,fromIntegral d))
+  where
+    d = fromIntegral . natVal $ (undefined :: Proxy n)
+
+dim :: forall n . KnownNat n => R n
+dim = mkR (scalar d)
+  where
+    d = fromIntegral . natVal $ (undefined :: Proxy n)
+
+
+--------------------------------------------------------------------------------
+
+
+ud2 :: L m n -> Matrix ℝ
+ud2 (L (Dim (Dim x))) = x
+
+
+--------------------------------------------------------------------------------
+
+diag :: KnownNat n => R n -> Sq n
+diag = diagR 0
+
+eye :: KnownNat n => Sq n
+eye = diag 1
+
+--------------------------------------------------------------------------------
+
+blockAt :: forall m n . (KnownNat m, KnownNat n) =>  ℝ -> Int -> Int -> Matrix Double -> L m n
+blockAt x r c a = mkL res
+  where
+    z = scalar x
+    z1 = LA.konst x (r,c)
+    z2 = LA.konst x (max 0 (m'-(ra+r)), max 0 (n'-(ca+c)))
+    ra = min (rows a) . max 0 $ m'-r
+    ca = min (cols a) . max 0 $ n'-c
+    sa = subMatrix (0,0) (ra, ca) a
+    m' = fromIntegral . natVal $ (undefined :: Proxy m)
+    n' = fromIntegral . natVal $ (undefined :: Proxy n)
+    res = fromBlocks [[z1,z,z],[z,sa,z],[z,z,z2]]
+
+
+
+
+
+--------------------------------------------------------------------------------
+
+
+row :: R n -> L 1 n
+row = mkL . asRow . ud1
+
+--col :: R n -> L n 1
+col v = tr . row $ v
+
+unrow :: L 1 n -> R n
+unrow = mkR . head . LA.toRows . ud2
+
+--uncol :: L n 1 -> R n
+uncol v = unrow . tr $ v
+
+
+infixl 2 ——
+(——) :: (KnownNat r1, KnownNat r2, KnownNat c) => L r1 c -> L r2 c -> L (r1+r2) c
+a —— b = mkL (extract a LA.—— extract b)
+
+
+infixl 3 ¦
+-- (¦) :: (KnownNat r, KnownNat c1, KnownNat c2) => L r c1 -> L r c2 -> L r (c1+c2)
+a ¦ b = tr (tr a —— tr b)
+
+
+type Sq n  = L n n
+--type CSq n = CL n n
+
+type GL = (KnownNat n, KnownNat m) => L m n
+type GSq = KnownNat n => Sq n
+
+isKonst :: forall m n . (KnownNat m, KnownNat n) => L m n -> Maybe (ℝ,(Int,Int))
+isKonst (unwrap -> x)
+    | singleM x = Just (x `atIndex` (0,0), (m',n'))
+    | otherwise = Nothing
+  where
+    m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
+    n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
+
+
+isKonstC :: forall m n . (KnownNat m, KnownNat n) => M m n -> Maybe (ℂ,(Int,Int))
+isKonstC (unwrap -> x)
+    | singleM x = Just (x `atIndex` (0,0), (m',n'))
+    | otherwise = Nothing
+  where
+    m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
+    n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
+
+
+
+infixr 8 <>
+(<>) :: forall m k n. (KnownNat m, KnownNat k, KnownNat n) => L m k -> L k n -> L m n
+(<>) = mulR
+
+
+infixr 8 #>
+(#>) :: (KnownNat m, KnownNat n) => L m n -> R n -> R m
+(#>) = appR
+
+
+infixr 8 <·>
+(<·>) :: R n -> R n -> ℝ
+(<·>) = dotR
+
+--------------------------------------------------------------------------------
+
+class Diag m d | m -> d
+  where
+    takeDiag :: m -> d
+
+
+instance forall n . (KnownNat n) => Diag (L n n) (R n)
+  where
+    takeDiag m = mkR (LA.takeDiag (extract m))
+
+
+instance forall m n . (KnownNat m, KnownNat n, m <= n+1) => Diag (L m n) (R m)
+  where
+    takeDiag m = mkR (LA.takeDiag (extract m))
+
+
+instance forall m n . (KnownNat m, KnownNat n, n <= m+1) => Diag (L m n) (R n)
+  where
+    takeDiag m = mkR (LA.takeDiag (extract m))
+
+
+--------------------------------------------------------------------------------
+
+linSolve :: (KnownNat m, KnownNat n) => L m m -> L m n -> Maybe (L m n)
+linSolve (extract -> a) (extract -> b) = fmap mkL (LA.linearSolve a b)
+
+(<\>) :: (KnownNat m, KnownNat n, KnownNat r) => L m n -> L m r -> L n r
+(extract -> a) <\> (extract -> b) = mkL (a LA.<\> b)
+
+svd :: (KnownNat m, KnownNat n) => L m n -> (L m m, R n, L n n)
+svd (extract -> m) = (mkL u, mkR s', mkL v)
+  where
+    (u,s,v) = LA.svd m
+    s' = vjoin [s, z]
+    z = LA.konst 0 (max 0 (cols m - size s))
+
+
+svdTall :: (KnownNat m, KnownNat n, n <= m) => L m n -> (L m n, R n, L n n)
+svdTall (extract -> m) = (mkL u, mkR s, mkL v)
+  where
+    (u,s,v) = LA.thinSVD m
+
+
+svdFlat :: (KnownNat m, KnownNat n, m <= n) => L m n -> (L m m, R m, L n m)
+svdFlat (extract -> m) = (mkL u, mkR s, mkL v)
+  where
+    (u,s,v) = LA.thinSVD m
+
+--------------------------------------------------------------------------------
+
+class Eigen m l v | m -> l, m -> v
+  where
+    eigensystem :: m -> (l,v)
+    eigenvalues :: m -> l
+
+newtype Sym n = Sym (Sq n) deriving Show
+
+
+sym :: KnownNat n => Sq n -> Sym n
+sym m = Sym $ (m + tr m)/2
+
+mTm :: (KnownNat m, KnownNat n) => L m n -> Sym n
+mTm x = Sym (tr x <> x)
+
+unSym :: Sym n -> Sq n
+unSym (Sym x) = x
+
+
+𝑖 :: Sized ℂ s c => s
+𝑖 = konst iC
+
+newtype Her n = Her (M n n)
+
+her :: KnownNat n => M n n -> Her n
+her m = Her $ (m + LA.tr m)/2
+
+
+
+instance KnownNat n => Eigen (Sym n) (R n) (L n n)
+  where
+    eigenvalues (Sym (extract -> m)) =  mkR . LA.eigenvaluesSH' $ m
+    eigensystem (Sym (extract -> m)) = (mkR l, mkL v)
+      where
+        (l,v) = LA.eigSH' m
+
+instance KnownNat n => Eigen (Sq n) (C n) (M n n)
+  where
+    eigenvalues (extract -> m) = mkC . LA.eigenvalues $ m
+    eigensystem (extract -> m) = (mkC l, mkM v)
+      where
+        (l,v) = LA.eig m
+
+--------------------------------------------------------------------------------
+
+withNullspace
+    :: forall m n z . (KnownNat m, KnownNat n)
+    => L m n
+    -> (forall k . (KnownNat k) => L n k -> z)
+    -> z
+withNullspace (LA.nullspace . extract -> a) f =
+    case someNatVal $ fromIntegral $ cols a of
+       Nothing -> error "static/dynamic mismatch"
+       Just (SomeNat (_ :: Proxy k)) -> f (mkL a :: L n k)
+
+
+withCompactSVD
+    :: forall m n z . (KnownNat m, KnownNat n)
+    => L m n
+    -> (forall k . (KnownNat k) => (L m k, R k, L n k) -> z)
+    -> z
+withCompactSVD (LA.compactSVD . extract -> (u,s,v)) f =
+    case someNatVal $ fromIntegral $ size s of
+       Nothing -> error "static/dynamic mismatch"
+       Just (SomeNat (_ :: Proxy k)) -> f (mkL u :: L m k, mkR s :: R k, mkL v :: L n k)
+
+--------------------------------------------------------------------------------
+
+qr :: (KnownNat m, KnownNat n) => L m n -> (L m m, L m n)
+qr (extract -> x) = (mkL q, mkL r)
+  where
+    (q,r) = LA.qr x
+
+-- use qrRaw?
+
+--------------------------------------------------------------------------------
+
+split :: forall p n . (KnownNat p, KnownNat n, p<=n) => R n -> (R p, R (n-p))
+split (extract -> v) = ( mkR (subVector 0 p' v) ,
+                         mkR (subVector p' (size v - p') v) )
+  where
+    p' = fromIntegral . natVal $ (undefined :: Proxy p) :: Int
+
+
+headTail :: (KnownNat n, 1<=n) => R n -> (ℝ, R (n-1))
+headTail = ((!0) . extract *** id) . split
+
+
+splitRows :: forall p m n . (KnownNat p, KnownNat m, KnownNat n, p<=m) => L m n -> (L p n, L (m-p) n)
+splitRows (extract -> x) = ( mkL (takeRows p' x) ,
+                             mkL (dropRows p' x) )
+  where
+    p' = fromIntegral . natVal $ (undefined :: Proxy p) :: Int
+
+splitCols :: forall p m n. (KnownNat p, KnownNat m, KnownNat n, KnownNat (n-p), p<=n) => L m n -> (L m p, L m (n-p))
+splitCols = (tr *** tr) . splitRows . tr
+
+
+toRows :: forall m n . (KnownNat m, KnownNat n) => L m n -> [R n]
+toRows (LA.toRows . extract -> vs) = map mkR vs
+
+
+toColumns :: forall m n . (KnownNat m, KnownNat n) => L m n -> [R m]
+toColumns (LA.toColumns . extract -> vs) = map mkR vs
+
+
+--------------------------------------------------------------------------------
+
+build
+  :: forall m n . (KnownNat n, KnownNat m)
+    => (ℝ -> ℝ -> ℝ)
+    -> L m n
+build f = mkL $ LA.build (m',n') f
+  where
+    m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
+    n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
+
+--------------------------------------------------------------------------------
+
+withVector
+    :: forall z
+     . Vector ℝ
+    -> (forall n . (KnownNat n) => R n -> z)
+    -> z
+withVector v f =
+    case someNatVal $ fromIntegral $ size v of
+       Nothing -> error "static/dynamic mismatch"
+       Just (SomeNat (_ :: Proxy m)) -> f (mkR v :: R m)
+
+
+withMatrix
+    :: forall z
+     . Matrix ℝ
+    -> (forall m n . (KnownNat m, KnownNat n) => L m n -> z)
+    -> z
+withMatrix a f =
+    case someNatVal $ fromIntegral $ rows a of
+       Nothing -> error "static/dynamic mismatch"
+       Just (SomeNat (_ :: Proxy m)) ->
+           case someNatVal $ fromIntegral $ cols a of
+               Nothing -> error "static/dynamic mismatch"
+               Just (SomeNat (_ :: Proxy n)) ->
+                  f (mkL a :: L m n)
+
+--------------------------------------------------------------------------------
+
+class Domain field vec mat | mat -> vec field, vec -> mat field, field -> mat vec
+  where
+    mul :: forall m k n. (KnownNat m, KnownNat k, KnownNat n) => mat m k -> mat k n -> mat m n
+    app :: forall m n . (KnownNat m, KnownNat n) => mat m n -> vec n -> vec m
+    dot :: forall n . (KnownNat n) => vec n -> vec n -> field
+    cross :: vec 3 -> vec 3 -> vec 3
+    diagR ::  forall m n k . (KnownNat m, KnownNat n, KnownNat k) => field -> vec k -> mat m n
+
+
+instance Domain ℝ R L
+  where
+    mul = mulR
+    app = appR
+    dot = dotR
+    cross = crossR
+    diagR = diagRectR
+
+instance Domain ℂ C M
+  where
+    mul = mulC
+    app = appC
+    dot = dotC
+    cross = crossC
+    diagR = diagRectC
+
+--------------------------------------------------------------------------------
+
+mulR :: forall m k n. (KnownNat m, KnownNat k, KnownNat n) => L m k -> L k n -> L m n
+
+mulR (isKonst -> Just (a,(_,k))) (isKonst -> Just (b,_)) = konst (a * b * fromIntegral k)
+
+mulR (isDiag -> Just (0,a,_)) (isDiag -> Just (0,b,_)) = diagR 0 (mkR v :: R k)
+  where
+    v = a' * b'
+    n = min (size a) (size b)
+    a' = subVector 0 n a
+    b' = subVector 0 n b
+
+mulR (isDiag -> Just (0,a,_)) (extract -> b) = mkL (asColumn a * takeRows (size a) b)
+
+mulR (extract -> a) (isDiag -> Just (0,b,_)) = mkL (takeColumns (size b) a * asRow b)
+
+mulR a b = mkL (extract a LA.<> extract b)
+
+
+appR :: (KnownNat m, KnownNat n) => L m n -> R n -> R m
+appR (isDiag -> Just (0, w, _)) v = mkR (w * subVector 0 (size w) (extract v))
+appR m v = mkR (extract m LA.#> extract v)
+
+
+dotR :: R n -> R n -> ℝ
+dotR (ud1 -> u) (ud1 -> v)
+    | singleV u || singleV v = sumElements (u * v)
+    | otherwise = udot u v
+
+
+crossR :: R 3 -> R 3 -> R 3
+crossR (extract -> x) (extract -> y) = vec3 z1 z2 z3
+  where
+    z1 = x!1*y!2-x!2*y!1
+    z2 = x!2*y!0-x!0*y!2
+    z3 = x!0*y!1-x!1*y!0
+
+--------------------------------------------------------------------------------
+
+mulC :: forall m k n. (KnownNat m, KnownNat k, KnownNat n) => M m k -> M k n -> M m n
+
+mulC (isKonstC -> Just (a,(_,k))) (isKonstC -> Just (b,_)) = konst (a * b * fromIntegral k)
+
+mulC (isDiagC -> Just (0,a,_)) (isDiagC -> Just (0,b,_)) = diagR 0 (mkC v :: C k)
+  where
+    v = a' * b'
+    n = min (size a) (size b)
+    a' = subVector 0 n a
+    b' = subVector 0 n b
+
+mulC (isDiagC -> Just (0,a,_)) (extract -> b) = mkM (asColumn a * takeRows (size a) b)
+
+mulC (extract -> a) (isDiagC -> Just (0,b,_)) = mkM (takeColumns (size b) a * asRow b)
+
+mulC a b = mkM (extract a LA.<> extract b)
+
+
+appC :: (KnownNat m, KnownNat n) => M m n -> C n -> C m
+appC (isDiagC -> Just (0, w, _)) v = mkC (w * subVector 0 (size w) (extract v))
+appC m v = mkC (extract m LA.#> extract v)
+
+
+dotC :: KnownNat n => C n -> C n -> ℂ
+dotC (unwrap -> u) (unwrap -> v)
+    | singleV u || singleV v = sumElements (conj u * v)
+    | otherwise = u LA.<·> v
+
+
+crossC :: C 3 -> C 3 -> C 3
+crossC (extract -> x) (extract -> y) = mkC (LA.fromList [z1, z2, z3])
+  where
+    z1 = x!1*y!2-x!2*y!1
+    z2 = x!2*y!0-x!0*y!2
+    z3 = x!0*y!1-x!1*y!0
+
+--------------------------------------------------------------------------------
+
+diagRectR :: forall m n k . (KnownNat m, KnownNat n, KnownNat k) => ℝ -> R k -> L m n
+diagRectR x v = mkL (asRow (vjoin [scalar x, ev, zeros]))
+  where
+    ev = extract v
+    zeros = LA.konst x (max 0 ((min m' n') - size ev))
+    m' = fromIntegral . natVal $ (undefined :: Proxy m)
+    n' = fromIntegral . natVal $ (undefined :: Proxy n)
+
+
+diagRectC :: forall m n k . (KnownNat m, KnownNat n, KnownNat k) => ℂ -> C k -> M m n
+diagRectC x v = mkM (asRow (vjoin [scalar x, ev, zeros]))
+  where
+    ev = extract v
+    zeros = LA.konst x (max 0 ((min m' n') - size ev))
+    m' = fromIntegral . natVal $ (undefined :: Proxy m)
+    n' = fromIntegral . natVal $ (undefined :: Proxy n)
+
+--------------------------------------------------------------------------------
+
+mean :: (KnownNat n, 1<=n) => R n -> ℝ
+mean v = v <·> (1/dim)
+
+test :: (Bool, IO ())
+test = (ok,info)
+  where
+    ok =   extract (eye :: Sq 5) == ident 5
+           && (unwrap .unSym) (mTm sm :: Sym 3) == tr ((3><3)[1..]) LA.<> (3><3)[1..]
+           && unwrap (tm :: L 3 5) == LA.matrix 5 [1..15]
+           && thingS == thingD
+           && precS == precD
+           && withVector (LA.vector [1..15]) sumV == sumElements (LA.fromList [1..15])
+
+    info = do
+        print $ u
+        print $ v
+        print (eye :: Sq 3)
+        print $ ((u & 5) + 1) <·> v
+        print (tm :: L 2 5)
+        print (tm <> sm :: L 2 3)
+        print thingS
+        print thingD
+        print precS
+        print precD
+        print $ withVector (LA.vector [1..15]) sumV
+        splittest
+
+    sumV w = w <·> konst 1
+
+    u = vec2 3 5
+
+    𝕧 x = vector [x] :: R 1
+
+    v = 𝕧 2 & 4 & 7
+
+    tm :: GL
+    tm = lmat 0 [1..]
+
+    lmat :: forall m n . (KnownNat m, KnownNat n) => ℝ -> [ℝ] -> L m n
+    lmat z xs = mkL . reshape n' . LA.fromList . take (m'*n') $ xs ++ repeat z
+      where
+        m' = fromIntegral . natVal $ (undefined :: Proxy m)
+        n' = fromIntegral . natVal $ (undefined :: Proxy n)
+
+    sm :: GSq
+    sm = lmat 0 [1..]
+
+    thingS = (u & 1) <·> tr q #> q #> v
+      where
+        q = tm :: L 10 3
+
+    thingD = vjoin [ud1 u, 1] LA.<·> tr m LA.#> m LA.#> ud1 v
+      where
+        m = LA.matrix 3 [1..30]
+
+    precS = (1::Double) + (2::Double) * ((1 :: R 3) * (u & 6)) <·> konst 2 #> v
+    precD = 1 + 2 * vjoin[ud1 u, 6] LA.<·> LA.konst 2 (size (ud1 u) +1, size (ud1 v)) LA.#> ud1 v
+
+
+splittest
+    = do
+    let v = range :: R 7
+        a = snd (split v) :: R 4
+    print $ a
+    print $ snd . headTail . snd . headTail $ v
+    print $ first (vec3 1 2 3)
+    print $ second (vec3 1 2 3)
+    print $ third (vec3 1 2 3)
+    print $ (snd $ splitRows eye :: L 4 6)
+ where
+    first v = fst . headTail $ v
+    second v = first . snd . headTail $ v
+    third v = first . snd . headTail . snd . headTail $ v
+
+
+instance (KnownNat n', KnownNat m') => Testable (L n' m')
+  where
+    checkT _ = test
+
+
-- 
cgit v1.2.3