diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,1 @@
+version 0.1.0.0
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c) 2013, Carter Tazio Schonwald and Wellposed Limited
+
+
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+    * Neither the name of Carter Tazio Schonwald nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,94 @@
+[![Wellposed](http://www.wellposed.com/mini.png)](http://www.wellposed.com)™
+
+# Currently in Pre alpha release engineering, so nearly ready for experimental consumption
+(though please dont publicize yet)
+
+# About  Numerical-Core
+This is the core Package for Numerical Haskell, a project by Carter Schonwald aka
+Wellposed Ltd, and (soon I hope!) other contributors.
+
+Numerical-Core is an open source component of the [Wellposed](http://www.wellposed.com)® Numerical Haskell software suite.
+
+##Build Status
+
+[![Build Status](https://secure.travis-ci.org/wellposed/numerical.png?branch=master)](http://travis-ci.org/wellposed/numerical)
+
+
+#note
+
+this library is **pre alpha release** so not all examples / codes may work as expected.
+That said, the current api should be enough to prototype and typecheck algorithms.
+
+
+
+# Numerical Haskell
+Numerical Haskell is an effort to bring great numerical computation and data analysis
+tooling to haskell, and be the best possible platform for sophisticated efforts in those same domains
+
+## What array Formats are Supported
+
+The initial
+
+## What convention is used for indexing?
+
+When you have an index tuple, just think  ``x,y,z``  to keep track of the meaning.
+Indexing tuples are written as statically sized lists, eg ``x:*y:*z:*Nil``.
+This follows the tradition of x,y,z axes used in plotting. Note well: the underlying memory
+order can be row OR column major or other!
+
+
+All the computations on these static sized lists get specialized away into
+nonrecursive computations at their use sites. So in this special scenario, lists aren't a problem!
+
+
+
+# Contributing
+Great! Theres so many awesome ways you could help out. Look at CONTRIBUTING.md for more details.
+Right now theres a lot of low hanging fruit in improving test coverage,
+and soon there'll be many opportunities on the performance tuning and numerical
+algorithms/tooling areas.
+
+## bug reports
+see bug.md for how to file a bug report
+
+
+# Performance FAQ
+1. How do I use Numerical haskell to write fast code thats outstandingly high level !?
+    * The leading cause of poor performance in numerical routines (aside from poor choice
+    in algorithms) is bad memory locality,
+    which has but a single easy cure: ** block recursive algorithms **
+    * Yes, you heard me, in compiled languages recursion is pretty cheap outside of the inner
+    most loops! It also is a fantastic tool for facilitating good memory locality!
+    * I'm totally serious, try out the benchmarks for the various versions of the same routines we
+    provide!
+2. But, what about fusion?
+    * Because of certain aspects of the numerical haskell design, we can't *automagically* use
+    the fusion optimization facilities of the underlying array representations such as Vector.
+
+# Community
+Many member of the Numerical Haskell community can be found on `#numerical-haskell` on freenode IRC
+There is also  the [numericalhaskell mailing list](https://groups.google.com/forum/#!forum/numericalhaskell)
+
+# Support
+The community provides some basic support through the IRC channel, Mailing list,
+and the relevant project [issue trackers](http://github.com/wellposed).
+
+If your support needs can't be resolved though those channels, please do not
+hesistate to contact Wellposed (aka Carter) to find out more about our support and
+professional services options.
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/numerical.cabal b/numerical.cabal
new file mode 100644
--- /dev/null
+++ b/numerical.cabal
@@ -0,0 +1,156 @@
+cabal-version:       2.2
+-- Initial numerics-types.cabal generated by cabal init.  For further
+-- documentation, see http://haskell.org/cabal/users-guide/
+
+-- The name of the package.
+name:                numerical
+
+-- The package version.  See the Haskell package versioning policy (PVP)
+-- for standards guiding when and how versions should be incremented.
+-- http://www.haskell.org/haskellwiki/Package_versioning_policy
+-- PVP summary:      +-+------- breaking API changes
+--                   | | +----- non-breaking API additions
+--                   | | | +--- code changes with no API change
+version:             0.0.0.0
+
+-- A short (one-line) description of the package.
+synopsis:           core package for Numerical Haskell project
+
+-- A longer description of the package.
+description:  the core package for Numerical Haskell. Still immature and incomplete
+
+-- The license under which the package is released.
+license:             BSD-2-Clause
+
+-- The file containing the license text.
+license-file:        LICENSE
+
+-- The package author(s).
+author:              Carter Tazio Schonwald
+
+-- An email address to which users can send suggestions, bug reports, and
+-- patches.
+maintainer:          carter at wellposed dot com
+
+-- A copyright notice.
+ copyright:     Carter Schonwald
+
+category:            Math
+
+build-type:          Simple
+tested-with:    GHC == 8.2.2 , GHC == 8.4.3, GHC == 8.6.2
+
+
+
+
+-- Constraint on the version of Cabal needed to build this package.
+
+
+
+extra-source-files:
+  README.md
+  CHANGELOG.md
+
+
+source-repository head
+  type: git
+  location: http://github.com/wellposed/numerical.git
+
+
+library
+    -- Modules exported by the library.
+    exposed-modules:
+        Numerical.Matrix.Basic
+        Numerical.Array.Layout
+        Numerical.Array.Layout.Base
+        Numerical.Array.Layout.Dense
+        Numerical.Array.Layout.Sparse
+        --Numerical.Array.Layout.Dense.Builder
+        Numerical.Array.Layout.Builder
+        Numerical.Array
+        Numerical.World
+        Numerical.InternalUtils
+        Control.NumericalMonad.State.Strict
+        Control.NumericalApplicative.Backwards
+        --Numerical.Array.Operations
+        Numerical.Array.Shape
+        Numerical.Array.Storage
+        Numerical.Nat
+        Numerical.Array.Mutable
+        Numerical.Array.Pure
+        Numerical.Array.Range
+        Numerical.Array.Locality
+        Numerical.Array.Address
+        Numerical.Data.Vector.Pair
+        Numerical.Data.Vector.HPair
+
+    ghc-options: -Wall
+    -- -ddump-simpl -ddump-to-file
+    --  --ghc-option=-ddump-simpl --ghc-option=-ddump-to-file  --ghc-options=-dsuppress-all
+    -- -O2
+    -- Modules included in this library but not exported.
+    -- other-modules:
+
+    if impl(ghc >= 8.0) && impl(ghc < 8.2)
+      ghc-options: -Wno-redundant-constraints
+
+
+    -- LANGUAGE extensions used by modules in this package.
+    other-extensions:
+        PolyKinds
+        BangPatterns
+        DataKinds
+        TypeFamilies
+        DeriveDataTypeable
+        TypeOperators
+        FlexibleInstances
+        FlexibleContexts
+        ScopedTypeVariables
+
+
+    -- Other library packages from which modules are imported.
+    build-depends:       base >= 4.10  && < 5
+
+                        ,primitive >= 0.5 && < 0.8
+                        ,vector >= 0.11 && < 0.13
+                        ,tagged >= 0.7 && < 0.9
+                        ,monad-ste >= 0.1 && < 0.2
+                        ,transformers >= 0.4 && < 0.6
+                        ,transformers-compat >= 0.4 && < 0.6
+                        ,ghc-prim >=0.2 && <0.6
+                        ,vector-algorithms >= 0.6.0.1 && < 0.9
+                        -- ,pqueue >= 1.2 && < 1.3
+                         -- , quickcheck >=
+
+                        -- what version constraints?
+
+
+
+    -- Directories containing source files.
+    hs-source-dirs:     src
+
+    -- Base language which the package is written in.
+    default-language:    Haskell2010
+
+Test-suite testsuite
+   default-language: Haskell2010
+   type: exitcode-stdio-1.0
+   build-depends:    base
+        ,hspec >=2.2 && <2.5
+        ,hspec-expectations
+        ,HUnit  >= 1.2.5 && < 1.7
+        ,primitive
+        ,vector
+        ,tagged
+        ,transformers
+        ,ghc-prim
+        ,numerical
+   ghc-options: -threaded
+   hs-source-dirs: tests
+   main-is: Main.hs
+   other-modules:
+        -- Only modules which are part of the test suite and not the library
+        -- should be listed here. These modules are excluded from the coverage
+        -- report because their coverage will be 100%.
+        NumericalUnit.Shape
+        NumericalUnit.Layout
diff --git a/src/Control/NumericalApplicative/Backwards.hs b/src/Control/NumericalApplicative/Backwards.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/NumericalApplicative/Backwards.hs
@@ -0,0 +1,58 @@
+
+
+module Control.NumericalApplicative.Backwards where
+
+import Prelude hiding (foldr, foldr1, foldl, foldl1)
+import qualified Control.Applicative as A
+import Data.Foldable as F
+import Data.Traversable as T
+
+-- | The same functor, but with an 'Applicative' instance that performs
+-- actions in the reverse order.
+newtype Backwards f a = Backwards { forwards :: f a }
+
+-- | Derived instance.
+instance (Functor f) => Functor (Backwards f) where
+    fmap f (Backwards a) = Backwards (fmap f a)
+    {-# INLINE fmap #-}
+
+-- | Apply @f@-actions in the reverse order.
+instance (A.Applicative f) => A.Applicative (Backwards f) where
+    pure a = Backwards (A.pure a)
+    {-# INLINE pure #-}
+    Backwards f <*> Backwards a = Backwards (a <**> f)
+    {-# INLINE (<*>) #-}
+
+
+-- | Try alternatives in the same order as @f@.
+instance (A.Alternative f) => A.Alternative (Backwards f) where
+    empty = Backwards A.empty
+    Backwards x <|> Backwards y = Backwards (x A.<|> y)
+
+-- | Derived instance.
+instance (Foldable f) => Foldable (Backwards f) where
+    foldMap f (Backwards t) = foldMap f t
+    foldr f z (Backwards t) = foldr f z t
+    foldl f z (Backwards t) = foldl f z t
+    foldr1 f (Backwards t) = foldl1 f t
+    foldl1 f (Backwards t) = foldr1 f t
+
+-- | Derived instance.
+instance (Traversable f) => Traversable (Backwards f) where
+    traverse f (Backwards t) = fmap Backwards (traverse f t)
+    sequenceA (Backwards t) = fmap Backwards (sequenceA t)
+    mapM f = A.unwrapMonad . T.traverse (A.WrapMonad . f)
+    sequence = T.mapM id
+    {-#INLINE traverse #-}
+    {-#INLINE sequenceA #-}
+    {-#INLINE mapM #-}
+    {-#INLINE sequence #-}
+
+(<**>) :: A.Applicative f => f a -> f (a -> b) -> f b
+(<**>) = liftA2 (flip ($))
+{-# INLINE (<**>) #-}
+
+liftA2 :: A.Applicative f => (a -> b -> c) -> f a -> f b -> f c
+liftA2 f a b = f `fmap` a A.<*> b
+{-# INLINE liftA2 #-}
+
diff --git a/src/Control/NumericalMonad/State/Strict.hs b/src/Control/NumericalMonad/State/Strict.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/NumericalMonad/State/Strict.hs
@@ -0,0 +1,246 @@
+
+module Control.NumericalMonad.State.Strict where
+
+--import Data.Functor.Identity
+import Control.Monad.Trans.Class
+import Control.Monad.IO.Class
+import Control.Applicative
+import Control.Monad
+import Control.Monad.Fix
+
+
+{-
+
+This module is a private copy of the Strict State Monad by Ross Patterson,
+patched to unconditionally inline.
+
+Its only purpose is to ensure that certain generic routines in
+Numerical.Array.Shape will compositionally unconditionally inline in their use sites
+
+ONLY use if writing generic code in your inner loops
+
+-}
+
+
+import Data.Foldable (Foldable(foldMap))
+import Data.Traversable (Traversable(traverse))
+
+-- | Identity functor and monad.
+newtype Identity a = Identity { runIdentity :: a }
+
+-- ---------------------------------------------------------------------------
+-- Identity instances for Functor and Monad
+
+instance Functor Identity where
+    fmap f m = Identity (f (runIdentity m))
+    {-# INLINE fmap #-}
+
+instance Foldable Identity where
+    foldMap f (Identity x) = f x
+    {-# INLINE foldMap #-}
+
+instance Traversable Identity where
+    traverse f (Identity x) = Identity <$> f x
+    {-# INLINE traverse #-}
+
+instance Applicative Identity where
+    pure a = Identity a
+    {-# INLINE pure #-}
+    Identity f <*> Identity x = Identity (f x)
+    {-# INLINE (<*>) #-}
+
+instance Monad Identity where
+    return a = Identity a
+    {-# INLINE return #-}
+    m >>= k  = k (runIdentity m)
+    {-# INLINE (>>=)#-}
+instance MonadFix Identity where
+    mfix f = Identity (fix (runIdentity . f))
+    {-# INLINE mfix  #-}
+
+
+-- ---------------------------------------------------------------------------
+-- | A state monad parameterized by the type @s@ of the state to carry.
+--
+-- The 'return' function leaves the state unchanged, while @>>=@ uses
+-- the final state of the first computation as the initial state of
+-- the second.
+type State s = StateT s Identity
+
+-- | Construct a state monad computation from a function.
+-- (The inverse of 'runState'.)
+state :: Monad m
+      => (s -> (a, s))  -- ^pure state transformer
+      -> StateT s m a   -- ^equivalent state-passing computation
+state  = \f -> StateT (return . f)
+{-# INLINE state #-}
+
+-- | Unwrap a state monad computation as a function.
+-- (The inverse of 'state'.)
+runState :: State s a   -- ^state-passing computation to execute
+         -> s           -- ^initial state
+         -> (a, s)      -- ^return value and final state
+runState  = \ m ->  runIdentity . runStateT m
+{-# INLINE runState#-}
+
+
+-- | Evaluate a state computation with the given initial state
+-- and return the final value, discarding the final state.
+--
+-- * @'evalState' m s = 'fst' ('runState' m s)@
+evalState :: State s a  -- ^state-passing computation to execute
+          -> s          -- ^initial value
+          -> a          -- ^return value of the state computation
+evalState  = \m s -> fst (runState m s)
+{-# INLINE evalState #-}
+
+-- | Evaluate a state computation with the given initial state
+-- and return the final state, discarding the final value.
+--
+-- * @'execState' m s = 'snd' ('runState' m s)@
+execState :: State s a  -- ^state-passing computation to execute
+          -> s          -- ^initial value
+          -> s          -- ^final state
+execState  =  \m s -> snd (runState m s)
+{-# INLINE execState#-}
+
+-- | Map both the return value and final state of a computation using
+-- the given function.
+--
+-- * @'runState' ('mapState' f m) = f . 'runState' m@
+mapState :: ((a, s) -> (b, s)) -> State s a -> State s b
+mapState   = \ f ->  mapStateT (Identity . f . runIdentity)
+{-# INLINE mapState #-}
+
+-- | @'withState' f m@ executes action @m@ on a state modified by
+-- applying @f@.
+--
+-- * @'withState' f m = 'modify' f >> m@
+withState :: (s -> s) -> State s a -> State s a
+withState = \f st -> withStateT f st
+{-# INLINE withState #-}
+-- ---------------------------------------------------------------------------
+-- | A state transformer monad parameterized by:
+--
+--   * @s@ - The state.
+--
+--   * @m@ - The inner monad.
+--
+-- The 'return' function leaves the state unchanged, while @>>=@ uses
+-- the final state of the first computation as the initial state of
+-- the second.
+newtype StateT s m a = StateT { runStateT :: s -> m (a,s) }
+
+-- | Evaluate a state computation with the given initial state
+-- and return the final value, discarding the final state.
+--
+-- * @'evalStateT' m s = 'liftM' 'fst' ('runStateT' m s)@
+evalStateT :: (Monad m) => StateT s m a -> s -> m a
+evalStateT  = \ m s -> do
+    (a, _) <- runStateT m s
+    return a
+{-# INLINE evalStateT #-}
+
+-- | Evaluate a state computation with the given initial state
+-- and return the final state, discarding the final value.
+--
+-- * @'execStateT' m s = 'liftM' 'snd' ('runStateT' m s)@
+execStateT :: (Monad m) => StateT s m a -> s -> m s
+execStateT  = \ m s -> do
+    (_, s') <- runStateT m s
+    return s'
+{-# INLINE  execStateT #-}
+
+
+-- | Map both the return value and final state of a computation using
+-- the given function.
+--
+-- * @'runStateT' ('mapStateT' f m) = f . 'runStateT' m@
+mapStateT :: (m (a, s) -> n (b, s)) -> StateT s m a -> StateT s n b
+mapStateT  = \ f m ->  StateT $ f . runStateT m
+
+-- | @'withStateT' f m@ executes action @m@ on a state modified by
+-- applying @f@.
+--
+-- * @'withStateT' f m = 'modify' f >> m@
+withStateT :: (s -> s) -> StateT s m a -> StateT s m a
+withStateT  = \ f m -> StateT $ runStateT m . f
+
+instance (Functor m) => Functor (StateT s m) where
+    fmap = \ f m  ->  StateT $ \ s ->
+        fmap (\ (a, s') -> (f a, s')) $ runStateT m s
+    {-# INLINE fmap  #-}
+
+instance (Functor m, Monad m) => Applicative (StateT s m) where
+    pure = \ a ->return a
+    (<*>) = \ a b ->  ap a b
+
+instance (Functor m, MonadPlus m) => Alternative (StateT s m) where
+    empty = mzero
+    {-# INLINE empty  #-}
+    (<|>) = \ a b -> mplus a b
+    {-#INLINE (<|>)#-}
+
+instance (Monad m) => Monad (StateT s m) where
+    {-# INLINE return #-}
+    return  = \ a ->  state $ \s -> (a, s)
+    {-# INLINE (>>=)#-}
+    (>>=)   = \m k ->  StateT $ \s -> do
+        (a, s') <- runStateT m s
+        runStateT (k a) s'
+    fail str = StateT $ \_ -> fail str
+
+instance (MonadPlus m) => MonadPlus (StateT s m) where
+    mzero       = StateT $ \_ -> mzero
+    {-# INLINE mzero #-}
+    mplus = \ m n -> StateT $ \s -> runStateT m s `mplus` runStateT n s
+    {-# INLINE mplus #-}
+instance (MonadFix m) => MonadFix (StateT s m) where
+    mfix  = \ f -> StateT $ \s -> mfix $ \ ~(a, _) -> runStateT (f a) s
+    {-# INLINE mfix #-}
+
+instance MonadTrans (StateT s) where
+    {-#INLINE lift #-}
+    lift   = \ m ->  StateT $ \s -> do
+        a <- m
+        return (a, s)
+
+instance (MonadIO m) => MonadIO (StateT s m) where
+    liftIO = lift . liftIO
+
+-- | Fetch the current value of the state within the monad.
+get :: (Monad m) => StateT s m s
+get = state $ \s -> (s, s)
+{-# INLINE get #-}
+
+-- | @'put' s@ sets the state within the monad to @s@.
+put :: (Monad m) => s -> StateT s m ()
+put  = \s ->  state $ \_ -> ((), s)
+{-# INLINE put #-}
+
+-- | @'modify' f@ is an action that updates the state to the result of
+-- applying @f@ to the current state.
+--
+-- * @'modify' f = 'get' >>= ('put' . f)@
+modify :: (Monad m) => (s -> s) -> StateT s m ()
+modify = \f ->  state $ \s -> ((), f s)
+{-# INLINE modify #-}
+
+-- | A variant of 'modify' in which the computation is strict in the
+-- new state.
+--
+-- * @'modify'' f = 'get' >>= (('$!') 'put' . f)@
+modify' :: (Monad m) => (s -> s) -> StateT s m ()
+modify' f = do
+    s <- get
+    put $! f s
+{-# INLINE  modify' #-}
+
+-- | Get a specific component of the state, using a projection function
+-- supplied.
+--
+-- * @'gets' f = 'liftM' f 'get'@
+gets :: (Monad m) => (s -> a) -> StateT s m a
+gets = \ f ->  state $ \s -> (f s, s)
+{-# INLINE gets #-}
+
diff --git a/src/Numerical/Array.hs b/src/Numerical/Array.hs
new file mode 100644
--- /dev/null
+++ b/src/Numerical/Array.hs
@@ -0,0 +1,106 @@
+{-# LANGUAGE PolyKinds   #-}
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE TypeFamilies #-}
+
+module Numerical.Array where
+
+--import Numerical.Array.Shape
+
+--data  MArray  world rep lay (view:: Locality) sh elem
+
+
+
+
+
+
+{-
+lets do just IO and not ST for now?
+or bite the primstate bullet now?
+-}
+
+--data family Array world rep lay (view:: Locality) sh elm
+
+
+{-
+only Row and Column Major have  dense formats that are unique across ALL
+possible ranks, not so simple for others. Make different data instances per formats
+
+
+
+
+
+-}
+--class Array  where
+
+-----------
+-- | for now locality is a closed type, may change going forward
+-- also should it be in a different module like shape or layout?
+
+
+
+{-
+theres several points in the design space of array apis that are nice, but none quite right
+
+Vector is probably the closest
+pros:
+  which has nice pure vs mutable apis
+  simple interface
+cons:
+  its really designed for Int indexing
+  assumes every pure Vector is internally derived from an imperative one
+    (this is reflected  in where the thaw/freeze)
+
+so there needs to be an Array class, an MArray class,
+
+and the Thawing / Freezing needs to be in a seperate PhasedArray class!
+why? Because we can't assume that pure/mutable arrays are the fundamental data type!
+
+-}
+
+
+{-
+maybe do
+     data Locality = Contiguous | Strided
+
+For now lets assume that the concrete (rather than delayed) arrays
+have a regular structure when strided. (rather than nonuniform gaps)
+-}
+{-
+rep = storable, unboxed, boxed, delay, etc
+
+lay = row major, column major, morton z, morton w (flipped n),
+  --- this  ignores symmetry  and hermitian being properties as well as packed layouts
+  --- also need to have a good sparse story
+    --- as currently done, most don't really make sense for != rank-2 arrays,
+
+-- rowMajor is a foldR, columnMajor is a foldL  over the shape ices
+
+-- Repa and accelerate use a Snoc List so that Row major fuses well for row major
+
+sh= rank / shape, ie matrix or vector, or some  higher tensor thingy
+lets borrow from  repa/ accelea
+
+
+mode= need to have a notion of runnable worlds,
+based on "backend" chosen, eg CBlasish, DPH, Repa, LLVM, Free (get the shallow/ deep ast)
+
+view =
+    Origin, Slice, and Diced, I might make this a fixed universe for now
+    Lets not distinguish whether a contiguous array is the original or derived for now
+    doesn't seem to be a meaningful difference and would make type inference crap / not bijective
+    Note that this does mean that accidental space leaks may happen
+
+    that suggests (but not for now) having a notion of origin / derived
+    that would allow elimiting space leaks. But lets not do that for now
+
+-}
+
+{-
+
+uncheckedReshape :: Array wld rep lay
+
+
+-}
+
+
diff --git a/src/Numerical/Array/Address.hs b/src/Numerical/Array/Address.hs
new file mode 100644
--- /dev/null
+++ b/src/Numerical/Array/Address.hs
@@ -0,0 +1,119 @@
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE DeriveDataTypeable #-}
+{-# LANGUAGE DeriveGeneric #-}
+{-# LANGUAGE TypeFamilies  #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+module Numerical.Array.Address(
+  Address(..)
+  ,SparseAddress(..)
+    ) where
+
+import Data.Data
+import Control.Monad (liftM)
+import qualified Foreign.Storable  as Store
+import qualified Data.Vector.Unboxed as UV
+import qualified Data.Vector.Generic as GV
+import qualified Data.Vector.Generic.Mutable as GMV
+import GHC.Generics
+
+-- | 'Address' is the type used for addressing into the underlying memory buffers
+-- of numerical arrays, Used for Dense Rank n arrays, and 1dim sparse arrays.
+newtype Address = Address  Int
+  deriving (Eq,Ord,Show,Read,Typeable,Generic,Data,Store.Storable)
+
+-- | 'LogicalAddress' is
+-- possibly dead code
+newtype LogicalAddress = LogicalAddress Int
+  deriving (Eq,Ord,Show,Read,Typeable,Generic,Data,Store.Storable)
+-- todo, add unboxed for
+
+
+-- | this m
+--newtype LogicalExtent
+-- sparse address seems to be dead atm
+data SparseAddress = SparseAddress {
+        outerIndex  :: {-# UNPACK #-} !Int
+        ,innerIndex :: {-# UNPACK #-} !Int }
+      deriving (Eq,Show,Data,Generic,Typeable)
+
+{-
+At some point decouple logical and physical address
+Logical Address should always be Int64 -- maybe even MORE?!
+physical address should be native IntPtr (aka Int)
+
+-}
+
+
+
+instance Num Address where
+    {-# INLINE (+) #-}
+    (+) (Address a) (Address b) = Address (a+b)
+    {-# INLINE (-) #-}
+    (-) (Address a) (Address b) =  Address (a-b)
+
+    (*) _ _ = error "you cant  multiply Addresses"
+    negate _ = error "you cant Apply Negate to An Address"
+    signum _ = error "error you cant take signum of an Address"
+    abs _ = error "error you cant take abs of an Address"
+    fromInteger _ = error "you cant use Integer Literals or fromInteger to form an Address"
+
+{-
+note that
+-}
+
+{-
+note that i don't think these vector instances ever matter
+-}
+
+newtype instance UV.MVector s Address  = MV_Address (UV.MVector s Int)
+newtype instance UV.Vector Address  = V_Address  (UV.Vector    Int)
+
+instance UV.Unbox Address where
+
+
+
+instance  GMV.MVector UV.MVector Address where
+  {-# INLINE basicLength #-}
+  {-# INLINE basicUnsafeSlice #-}
+  {-# INLINE basicOverlaps #-}
+  {-# INLINE basicUnsafeNew #-}
+  {-# INLINE basicUnsafeReplicate #-}
+  {-# INLINE basicUnsafeRead #-}
+  {-# INLINE basicUnsafeWrite #-}
+  {-# INLINE basicClear #-}
+  {-# INLINE basicSet #-}
+  {-# INLINE basicUnsafeCopy #-}
+  {-# INLINE basicUnsafeGrow #-}
+  {-# INLINE basicInitialize #-}
+  basicInitialize = \ (MV_Address mva) -> GMV.basicInitialize mva
+  basicLength (MV_Address v) = GMV.basicLength v
+  basicUnsafeSlice i n (MV_Address v) = MV_Address $ GMV.basicUnsafeSlice i n v
+  basicOverlaps (MV_Address v1) (MV_Address v2) = GMV.basicOverlaps v1 v2
+  basicUnsafeNew n = MV_Address `liftM` GMV.basicUnsafeNew n
+  basicUnsafeReplicate n (Address a) = MV_Address `liftM` GMV.basicUnsafeReplicate n a
+  basicUnsafeRead (MV_Address v) i = Address `liftM` GMV.basicUnsafeRead v i
+  basicUnsafeWrite (MV_Address v) i (Address a) = GMV.basicUnsafeWrite v i a
+  basicClear (MV_Address v) = GMV.basicClear v
+  basicSet (MV_Address v) (Address a) = GMV.basicSet v a
+  basicUnsafeCopy (MV_Address v1) (MV_Address v2) = GMV.basicUnsafeCopy v1 v2
+  basicUnsafeMove (MV_Address v1) (MV_Address v2) = GMV.basicUnsafeMove v1 v2
+  basicUnsafeGrow (MV_Address v) n = MV_Address `liftM` GMV.basicUnsafeGrow v n
+
+instance  GV.Vector UV.Vector Address where
+  {-# INLINE basicUnsafeFreeze #-}
+  {-# INLINE basicUnsafeThaw #-}
+  {-# INLINE basicLength #-}
+  {-# INLINE basicUnsafeSlice #-}
+  {-# INLINE basicUnsafeIndexM #-}
+  {-# INLINE elemseq #-}
+  basicUnsafeFreeze (MV_Address v) = V_Address `liftM` GV.basicUnsafeFreeze v
+  basicUnsafeThaw (V_Address v) = MV_Address`liftM` GV.basicUnsafeThaw v
+  basicLength (V_Address v) = GV.basicLength v
+  basicUnsafeSlice i n (V_Address v) = V_Address $ GV.basicUnsafeSlice i n v
+  basicUnsafeIndexM (V_Address v) i
+                = Address `liftM` GV.basicUnsafeIndexM v i
+  basicUnsafeCopy (MV_Address mv) (V_Address v)
+                = GV.basicUnsafeCopy mv v
+  elemseq _ (Address a) z =   GV.elemseq (undefined :: UV.Vector a) a z
+
diff --git a/src/Numerical/Array/Layout.hs b/src/Numerical/Array/Layout.hs
new file mode 100644
--- /dev/null
+++ b/src/Numerical/Array/Layout.hs
@@ -0,0 +1,18 @@
+
+
+module Numerical.Array.Layout(
+  module Numerical.Array.Layout.Base
+  ,module Numerical.Array.Layout.Dense
+  ,module Numerical.Array.Layout.Sparse
+  ,module Numerical.Array.Address
+  ) where
+
+
+import  Numerical.Array.Layout.Base
+import  Numerical.Array.Layout.Dense
+import  Numerical.Array.Layout.Sparse
+import Numerical.Array.Address
+
+
+
+
diff --git a/src/Numerical/Array/Layout/Base.hs b/src/Numerical/Array/Layout/Base.hs
new file mode 100644
--- /dev/null
+++ b/src/Numerical/Array/Layout/Base.hs
@@ -0,0 +1,554 @@
+{- |  Comments for this modules
+
+
+-}
+
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE DeriveDataTypeable #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE NoImplicitPrelude #-}
+{-# LANGUAGE FunctionalDependencies #-}
+{-# LANGUAGE CPP #-}
+
+{-# LANGUAGE StandaloneDeriving #-}
+
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE UndecidableInstances #-}
+
+#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 707
+ {-# LANGUAGE AutoDeriveTypeable #-}
+#endif
+
+
+module Numerical.Array.Layout.Base(
+  Layout(..)
+  ,DenseLayout(..)
+  ,RectilinearLayout(..)
+  ,LayoutAddress
+  ,LayoutLogicalFormat
+  ,Transposed
+  ,FormatStorageRep
+  ,RectOrientationForm
+  ,RectDownRankForm
+  ,InnerContigForm
+  ,Format
+  ,TaggedShape(..)
+  ,GDSlice(..) --- right? right?
+  ,SMajorOrientation(..)
+  ,MajorOrientation(..)
+  ,majorCompareRightToLeft
+  ,majorCompareLeftToRight
+  ,shapeCompareRightToLeft
+  ,shapeCompareLeftToRight
+  -- * All the various helper types
+  ,module Numerical.Array.Storage
+  ,module Numerical.Array.Locality
+  ,module Numerical.Array.Shape
+  ,module Numerical.Array.Range
+  ,module Numerical.Array.Address
+) where
+
+
+import Data.Data
+
+import Numerical.Nat
+import Numerical.Array.Address
+import Numerical.Array.Locality
+import Numerical.Array.Shape
+import Numerical.Array.Storage
+import Numerical.Array.Range
+
+--import Data.Typeable
+#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ < 709
+import  qualified Control.Applicative as A
+import Prelude hiding (foldr,foldr1,foldl1,foldl,map)
+import  qualified  Data.Foldable as F
+#elif __GLASGOW_HASKELL__ >= 709
+import  qualified Control.Applicative as A
+import  qualified  Data.Foldable as F
+#endif
+
+#if MIN_VERSION_base(4,8,0)
+import Prelude hiding (foldl)
+#endif
+
+{-
+NB: may need to add some specialization for low rank indexing,
+theres 4 choices:
+a) INLINE EVERYTHING
+b) rewrite rules that take low rank indexing code into specialized versions thereof
+c) wait till ghc 7.8.2 to resolve https://ghc.haskell.org/trac/ghc/ticket/8848
+    and use SPECIALIZE
+d) benchmark and then decide
+
+for now I choose (a), and defer benchmarking variations till everything works :)
+
+
+a related concern is the interplay of inlining and specialization
+https://ghc.haskell.org/trac/ghc/ticket/5928
+
+-}
+
+
+
+
+
+-- either we need to break ties, or the ties have been broken
+majorCompareLeftToRight :: Ordering -> Ordering -> Ordering
+majorCompareLeftToRight EQ new = new
+majorCompareLeftToRight a _ = a
+
+
+majorCompareRightToLeft :: Ordering -> Ordering -> Ordering
+majorCompareRightToLeft new EQ = new
+majorCompareRightToLeft _ b = b
+
+{-# INLINE shapeCompareLeftToRight #-}
+shapeCompareLeftToRight :: (F.Foldable (Shape r),A.Applicative (Shape r), Ord a)
+    => Shape r a -> Shape r a -> Ordering
+shapeCompareLeftToRight =   \  ls rs -> foldl majorCompareLeftToRight EQ  $ map2 compare ls rs
+
+{-# INLINE shapeCompareRightToLeft #-}
+shapeCompareRightToLeft :: ((F.Foldable (Shape r)),A.Applicative (Shape r), Ord a)
+   => Shape r a -> Shape r a -> Ordering
+shapeCompareRightToLeft =   \  ls rs -> foldl majorCompareRightToLeft EQ  $ map2 compare ls rs
+
+
+
+
+-- | this is kinda a hack
+newtype TaggedShape (form :: *) (rank::Nat) = TaggedShape {unTagShape:: Shape rank Int }
+instance Eq (Shape rank Int)=> Eq (TaggedShape f rank) where
+  (==) l r =  (==) (unTagShape l) (unTagShape r )
+
+instance Show (Shape rank Int) => Show (TaggedShape f rank) where
+  show (TaggedShape ix) =  "TaggedShape (" ++ show ix ++ " )"
+
+instance forall form  rank . (Eq (Shape rank Int),Layout form rank)
+  => Ord (TaggedShape form rank) where
+  compare left right = basicCompareIndex (Proxy:: Proxy form ) (unTagShape left) (unTagShape right)
+
+
+-- | Generalized Dense Slice Projection notation,
+-- not sure if it should be defined in this module or elsewhere
+-- This provides a type safe interface for the classical
+-- general array slice notation.
+-- That said, its only useful for dense array formats,
+-- at least in general. For formats that aren't "rectilinear dense",
+-- this COULD be used as a description format for traversing
+-- over various rectilinear subsets of points though?
+data GDSlice (from :: Nat) (to :: Nat) :: *  where
+  GDNil :: GDSlice 'Z 'Z
+  GDPick :: Int -> !(GDSlice from to) -> GDSlice ('S from) to
+  GDRange :: (Int,Int,Int) {- this is a nonempty interval or error -} -> !(GDSlice from to) -> GDSlice ('S from) ('S to)
+  GDAll :: !(GDSlice from to) -> GDSlice ('S from) ('S to)
+
+{-
+TODO: for things that
+
+-}
+
+
+instance Show (GDSlice 'Z 'Z) where
+  show _ = "GDNil"
+
+instance (Show (GDSlice (f) ('S t)),Show (GDSlice f t))=> Show (GDSlice ('S f) ('S t)) where
+  show (tup `GDRange` rest) = show tup ++ " `GDRange` (" ++ show rest ++ ")"
+  show (GDAll rest) =  "GDAll " ++ show rest
+  show (ix `GDPick` rest) = show ix ++" `GDPick` " ++ show rest
+
+
+instance Show (GDSlice f 'Z)=> Show (GDSlice ('S f) 'Z) where
+  show (ix `GDPick` rest) = show ix ++" `GDPick` " ++ show rest
+--instance Show (GDSlice f t)  where
+--  func =
+
+{-
+In some (moderately precise sense)
+
+-}
+
+
+-- GDRange (from,step,to)
+  -- GDAll is just sugar for a special case of GDRange, but maybe its worthwhile sugar?
+
+--computeSlicePlan:: GDSlice from to -> Shape from Int -> Shape from (Either Int (AffineRange Int))
+--computeSlicePlan GDNil  Nil = Nil
+--computeSlicePlan  ( ix `GDPick` gdRest )
+--                  (bd:* shpRest)| ix < bd   && ix >= 0 = Left ix :* computeSlicePlan gdRest shpRest
+--                      | otherwise = error
+--                          $ "bad indices for computeSlicePlan " ++ show (ix,bd)
+--computeSlicePlan ( (strt,step,end) `GDRange` grest) (bd:* shprest)
+
+
+
+data family Format  lay (contiguity:: Locality)  (rank :: Nat) rep
+
+deriving instance Typeable Format
+
+type family FormatStorageRep ( a:: * ) :: *
+
+type instance FormatStorageRep (Format lay ctg rnk rep)= rep
+
+type family  Transposed (form :: *) :: *
+
+type family  LayoutAddress (form :: *) :: *
+
+-- TODO / FIXME remove the basic* prefix  from all the operations
+-- this was done originally because
+
+
+-- TODO : should this be pushed into the type class?
+-- TODO : should this be pushed into the type class?
+-- | every format has a "logical" sibling, that represents the address translation
+-- when the underlying buffer layer is contiguous and packed. So it could be claimed
+-- that  any type that obeys @a~'LayoutLogicalFormat' a@ is one that an be a legal
+-- instance of LayoutBuilder?
+type family LayoutLogicalFormat (form :: *) :: *
+
+-- | the 'Layout' type class
+class Layout form  (rank :: Nat) | form -> rank  where
+
+    -- | 'basicLogicalShape' gives the extent of the format
+    basicLogicalShape :: form -> Shape rank Int
+
+    -- | 'basicLogicalForm' converts a given format into its "contiguous" analogue
+    -- this is useful for supporting various address translation manipulation tricks
+    -- efficiently. Note that any valid  simple format should strive to ensure this is an O(1) operation.
+    -- though certain composite 'Layout' instances may provide a slower implementation.
+    basicLogicalForm :: (logicalForm ~ LayoutLogicalFormat form ) => form -> logicalForm
+
+
+    -- | 'transposedLayout' transposes the format data type
+    transposedLayout :: (form ~ Transposed transform,transform~Transposed form)=> form  -> transform
+
+    -- | 'basicCompareIndex' lets you compare where two (presumably inbounds)
+    -- 'Index' values are in a formats ordering. The logical 'Shape' of the array
+    -- is not needed
+    basicCompareIndex :: p form-> Shape rank Int ->Shape rank Int -> Ordering
+
+    -- | the (possibly empty) min and max of the valid addresses for a given format.
+    -- @minAddress = fmap _RangeMin . rangedFormatAddress@
+    -- and @maxAddress = fmap _RangeMax . rangedFormatAddress@
+    -- FIXME : This also is a terrible name
+    basicAddressRange ::  (address ~ LayoutAddress form)=> form -> Maybe (Range address)
+    -- FIX ME! this name is crap, i dont like it
+
+    -- | 'basicToAddress' takes an Index, and tries to translate it to an address if its in bounds
+    --
+    basicToAddress :: (address ~ LayoutAddress form)=>
+        form  -> Index rank  -> Maybe  address
+
+    -- | 'basicToIndex' takes an address, and always successfully translates it to
+    -- a valid index. Behavior of invalid addresses constructed by a library user
+    -- is unspecified.
+    basicToIndex ::(address ~ LayoutAddress form)=>
+        form -> address -> Index rank
+
+    -- | 'basicNextAddress' takes an address, and tries to compute the next valid
+    -- address, or returns Nothing if there is no subsequent valid address.
+    basicNextAddress :: (address ~ LayoutAddress form)=>
+        form  -> address -> Maybe  address
+
+    -- |  @'basicNextIndex' form ix mbeAddress@  computes the next valid index after
+    -- @ix@ if it exists. It takes a @'Maybe' address@ as a hint for where to do the search for the successor.
+    -- If the index is in bounds and not the last index, it returns both the index and the associated address.
+    basicNextIndex :: (address ~ LayoutAddress form)=>
+          form  -> Index rank -> Maybe address  -> Maybe ( Index rank, address)
+
+
+    basicAddressPopCount :: (address ~ LayoutAddress form)=>
+        form -> Range address -> Int
+
+    -- | This operation is REALLY unsafe
+    -- This should ONLY be used on Formats that are directly
+    -- paired with a Buffer or Mutable Buffer (ie a Vector)
+    --  This operation being in this class is also kinda a hack
+    -- but lets leave it here for now
+    basicAddressAsInt :: (address ~ LayoutAddress form)=>
+        form ->  address -> Int
+    basicAddressAsInt =
+       \ _ _ ->
+        error "called basicAddressAsInt on a Layout thats not meant for this world"
+
+    -- | The semantics of @`basicAffineAddressShift` form addr step@ is that
+    -- when  step > 0, its equivalent to iteratively computing 'basicNextAddress' @step@ times.
+    -- However, the step size can be negative, which means it can
+    basicAffineAddressShift :: (address ~ LayoutAddress form) =>
+        form -> address -> Int -> Maybe address
+
+
+    {-# MINIMAL basicToAddress, basicToIndex, basicNextAddress,basicNextIndex
+          ,basicAddressRange,basicLogicalShape,basicCompareIndex
+          , transposedLayout, basicAddressPopCount,basicLogicalForm, basicAffineAddressShift #-}
+
+
+{- |
+these names aren't ideal, but lets punt on bikeshedding till theres >= 2 serious
+users
+-}
+data MajorOrientation = Rowed | Columned | BlockedColumn | BlockedRow
+  deriving(Data,Typeable)
+
+data SMajorOrientation (o :: MajorOrientation) where
+    SRowed :: SMajorOrientation 'Rowed
+    SColumned :: SMajorOrientation 'Columned
+    SBlockedRow :: SMajorOrientation 'BlockedRow
+    SBlockedColumn :: SMajorOrientation 'BlockedColumn
+
+
+-- |  Every instance of 'RectilinearLayout' needs to have a corresponding
+-- 'RectOrientationForm', 'RectDownRankForm', and 'InnerContigForm'
+type family RectOrientationForm form :: MajorOrientation
+
+type family RectDownRankForm   form :: *
+
+type family InnerContigForm form :: *
+
+{- | 'RectilinearLayout' is the type class that supports the modle widely
+  usable class of slicing operations in Numerical.
+  for every instance @'RectilinearLayout' format rank orientation@, a corresponding
+  @'RectOrientationForm' form @, @'RectDownRankForm' form@
+  and @'InnerContigForm' form@ type family instance should be defined
+
+  The purpose of 'RectilinearLayout' class is to provide
+
+-}
+class Layout form rank =>
+  RectilinearLayout form (rank :: Nat) (oriented :: MajorOrientation) | form -> rank oriented where
+
+    -- | 'formRectOrientation' provides a runtime mechanism for reflecting
+    -- the orientation of the format
+    formRectOrientation :: p form -> SMajorOrientation oriented
+
+{-
+is array layout always static?
+for now lets say yes, cause you can always just existential up the class
+
+-}
+
+    -- | For  @'rectlinearShape' form==shp@, we always have that
+    -- @'basicLogicalShape' form  `weaklyDominates` shp@.
+    -- when 'strictlyDominates' holds, that implies that the underlying array format
+    -- is a rectilinear layout whose "elements" are tiles of a fixed size array format.
+    -- For this initial release and initial set of applicable rectilinear array formats,
+    -- the following is always true @'basicLogicalShape' form  == basicLogicalShape' form @
+    -- Should be @O(1)@ always. Or more precisely @O(rank)@
+    rectlinearShape :: form -> Index rank
+
+    unconsOuter:: ('S down ~ rank)=> p form -> Shape rank a -> (a, Shape down a)
+    consOuter ::  ('S down ~ rank)=> p form -> a -> Shape down a -> Shape rank a
+
+    -- | @'majorAxisSlice' fm (x,y)@ requires that y-x>=1, ie that more than
+    -- one sub range wrt the major axis be selected, so that the logical
+    -- rank of the selected array stays the same. This operation also preserves
+    -- memory locality as applicable.
+    -- @O(1)@ / @O(rank)@
+    majorAxisSlice :: form -> (Int,Int)-> form
+     -- should this be -> Maybe form?
+
+
+    -- | @'majorAxixProject' form x@ picks a "row" with respect to the outer most
+    -- dimension of the array format. This will
+    -- @O(1)@ or @O(rank)@
+    majorAxisProject :: (RectilinearLayout downForm subRank oriented,
+     rank ~ ('S subRank) , downForm~ RectDownRankForm form) => form -> Int -> downForm
+
+    -- | this is the nonstrided subset of general array slice notation.
+    --  Invoke as @'rectilinearSlice'  form  leastCorner greatestCorner@,
+    -- where the least and greatest corners of the sub array are determined
+    -- by the 'strictlyDominates' partial order on the bounds of the sub array.
+    -- For Dense array formats, this should be @O(1)@ or more precisely @O(rank)@
+    -- For the basic Sparse array formats thus far the complexity should be
+    -- @O(size of outermost dimension)@, which could be computed by
+    --  @fst . unconsOuter [form] . rectilinearShape $ form@
+    rectlinearSlice :: (RectilinearLayout icForm rank oriented,icForm~InnerContigForm form )=>form -> Index rank -> Index rank -> icForm -- FIXME, need the range infos????? (icfFOrm, adddress,address)
+
+
+{- | 'DenseLayout' only has instances for Dense array formats.
+this class will need some sprucing up for the beta,
+but its ok for now. NB that 'DenseLayout' is really strictly meant to be used
+for optimization purposes, and not meant as a default api
+-}
+class Layout form rank =>  DenseLayout form  (rank :: Nat) | form -> rank  where
+
+
+
+    basicToDenseAddress :: form  -> Index rank  ->   Address
+
+    basicToDenseIndex :: form -> Address -> Index rank
+
+
+
+    basicNextDenseAddress :: form  -> Address ->  Address
+    basicNextDenseAddress =  \form shp -> snd
+      (basicNextDenseIndex form  $ basicToDenseIndex form  shp )
+    {-# INLINE basicNextDenseAddress #-}
+
+    basicNextDenseIndex :: form  -> Index rank ->(Index rank ,Address)
+    basicNextDenseIndex  = \form shp -> (\ addr ->( basicToDenseIndex form addr, addr) ) $!
+       basicNextDenseAddress form  $ basicToDenseAddress form  shp
+    {-# INLINE  basicNextDenseIndex #-}
+
+
+#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 707
+    {-# MINIMAL  basicToDenseIndex, basicToDenseAddress,
+     (basicNextDenseIndex | basicNextDenseAddress)   #-}
+#endif
+
+{-
+*Numerical.Array.Layout> basicToAddress (FormColumn (2 :* 3 :* 7 :* Nil)) (0:* 2 :* 2 :* Nil)
+Address 16
+*Numerical.Array.Layout> basicToAddress (FormColumn (2 :* 3 :* 7 :* Nil)) (1:* 0 :* 0 :* Nil)
+Address 1
+*Numerical.Array.Layout> basicToAddress (FormColumn (2 :* 3 :* 7 :* Nil)) (0:* 0 :* 0 :* Nil)
+Address 0
+*Numerical.Array.Layout> basicToAddress (FormColumn (2 :* 3 :* 7 :* Nil)) (0:* 1 :* 0 :* Nil)
+Address 2
+*Numerical.Array.Layout> basicToAddress (FormColumn (2 :* 3 :* 7 :* Nil)) (0:* 0 :* 1 :* Nil)
+
+
+
+-}
+
+
+--data Elem ls el  where
+--    Point :: Elem '[] el
+--    (:#) :: a -> Elem ls el -> Elem (a ': ls) el
+
+
+{-
+    One important invariant about all layouts at all ranks is that for
+    any given ints x < y, that the array index for inr
+
+     toIndex shapedLayout (pure x :: Shape rank Int) is strictly less than
+     toIndex shapedLayout (pure y :: Shape rank Int).
+
+     more generally
+
+     for rank k tuples,
+      xi = x_1 :* ... :* x_k *: Nil  and
+      yj = y_1 :* ... :* x_k *: Nil
+      such that forall \ell, x_\ell  < y_\ell
+    we have that
+       toIndex shapedLayout xi <  toIndex  shapedLayout yj
+
+
+this actually relates to the notion of partial ordering over vectors in convex
+geometry!
+
+
+so roughly: we have layouts that are dense
+we have layouts that can be used as tiles (and are dense)
+
+and we have layouts which can can't be tiled, but can have elements which are tiled
+
+So we have
+
+PrimitiveLayouts
+
+Static Layouts
+
+General Layouts (which are a Top level layout over a static layout)
+
+the Layout class tries to abstract over all three cases
+(NB: this only makes sense when the "rank" for the inner
+and outer layouts have the same rank!)
+
+-}
+
+
+{- Sized is used as a sort of hack to make it easy to express
+   the staticly sized layouts. NB, one trade off is that its only
+   possible to express  "cube" shaped blocks, but on the other
+   hand blocking sizes are expressible for every single rank!
+-}
+--data Sized :: * -> * where
+    --(:@) :: Nat -> a -> Sized a
+
+
+{-
+
+per se I don't need the StaticLay, PrimLay, Lay constructors, BUT
+I really do like how it makes things a teeny bit simpler.. though I may remove them
+-}
+
+
+
+--class SimpleDenseLayout lay (rank :: Nat) where
+--  type SimpleDenseTranpose lay
+--  toIndexSimpleDense :: Shaped rank Int lay -> Shape rank Int -> Int
+
+
+--class PrimLayout lay (rank :: Nat) where
+--    type TranposedPrim lay
+--    toIndexPrim :: Shaped rank Int (PrimLay lay) -> Shape rank Int -> Int
+--    fromIndexPrim :: Shaped rank Int (PrimLay lay) -> Int -> Shape rank Int
+
+
+{-
+for now we will not deal with nested formats, but this will
+be a breaking change i plan for later
+-}
+
+{-
+what is the law for the Layout class?
+forall valid formms
+toIndex sd  (fromIndex sd ix)==ix
+fromIndex sd (toIndex sd shp)==shp
+-}
+
+{-
+if   tup1 is strictly less than tup2 (pointwise),
+  then any lawful Layout will asign tup1 an index strictly less than that
+  asigned to tup2
+
+  transposedLayout . transposedLayout == id
+
+
+
+i treat coordinates as being in x:* y :* z :* Nil, which is Fortran style idexing
+
+in row major we'd have for x:* y :* Nil that X is the inner dimension, and y the outter,
+by contrast, in column major, y would be the inner most, and x the outter most.
+
+
+
+
+-}
+
+
+{- In some respects, the Layout type class is a multidimensional
+analogue of the Enum type class in Haskell Prelude,
+for Dense / Dense Structured matrix formats
+but
+    a) requires a witness value, the "Form"
+    b) needs to handle multivariate structures
+    c) has to deal with structure matrices, like triangular, symmetric, etc
+    e) I think every layout should have pure 0 be a valid index, at least for "Dense"
+    arrays
+    f) transposedLayout . transposedLayout == id
+    g)
+
+  Form needs to carry the shape / extent of the matrix
+
+-}
+{-
+
+-}
+
+--data View = Origin | Slice
+{-
+i'm really really hoping to not need a View parameter,
+but the nature of the addressing logic needs to change when its a slice
+vs a deep copy (for certain classes of arrays that I wish to support very easily)
+
+I will be likely adding this the moment benchmarks validate the distinction
+
+on the
+-}
diff --git a/src/Numerical/Array/Layout/Builder.hs b/src/Numerical/Array/Layout/Builder.hs
new file mode 100644
--- /dev/null
+++ b/src/Numerical/Array/Layout/Builder.hs
@@ -0,0 +1,310 @@
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE DeriveDataTypeable #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FunctionalDependencies #-}
+{-#  LANGUAGE GADTs #-}
+{-# LANGUAGE RankNTypes  #-}
+{-# LANGUAGE ScopedTypeVariables#-}
+-- {-# LANGUAGE DeriveFunctor #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE UndecidableInstances #-}
+{-# LANGUAGE NoImplicitPrelude #-}
+
+module Numerical.Array.Layout.Builder where
+
+import Control.Monad.Primitive ( PrimMonad, PrimState )
+import qualified Data.Vector.Generic as VG
+import qualified Data.Vector.Generic.Mutable as VGM
+import Numerical.Array.Layout.Base
+import  Numerical.Array.Layout.Dense as Dense
+--import Numerical.Array.Layou.Sparse as Sparse
+--import Numerical.Data.Vector.Pair
+import Control.Monad.ST (runST)
+import Data.Typeable
+import qualified  Data.Foldable as F
+import   Data.Traversable as T
+import   Control.Applicative as A
+
+import Numerical.Data.Vector.Pair
+import Numerical.Array.Layout.Sparse
+import Data.Vector.Algorithms.Intro as IntroSort
+import Data.List (group)
+import Numerical.InternalUtils
+import Prelude hiding (error)
+
+
+data BatchInit   v = BatchInit    { batchInitSize :: !Int
+             ,batchInitKV :: !(Either [v]  (IntFun v))    }
+            deriving (Typeable)
+
+materializeBatchMV :: (PrimMonad m, VGM.MVector mv a)  => BatchInit a  -> m (mv (PrimState m) a)
+materializeBatchMV  (BatchInit size (Left ls )) =
+         do
+            newMV <- VGM.new size
+            _ <- Prelude.mapM (\(ix ,val )-> VGM.unsafeWrite newMV  ix val ) (zip [0..] $ take size ls)
+            return newMV
+materializeBatchMV  (BatchInit size (Right (IntFun f) )) =
+         do
+            newMV <- VGM.new size
+            _ <- Prelude.mapM (\ix -> do v <- (f ix) ; VGM.unsafeWrite newMV  ix  v ) $ take size  [0..]
+            return newMV
+
+--- not sure if this is EVER useful
+newtype AnyMV mv e = AMV (forall s . mv s e )
+
+
+instance (Show  a)=> Show (BatchInit a) where
+  show (BatchInit size (Left ls) )  | size > 100 =  "(BatchInit " ++show size  ++
+                                          "-- only showing the first 100 elements\n"
+                                          ++ "(Left "++(show $ take 100 ls ) ++ "))\n"
+                                    | otherwise ="(BatchInit " ++show size  ++
+                                     " (Left "++(show  ls ) ++ "))\n"
+  show (BatchInit size (Right (IntFun f)) ) | size > 100 =  "(BatchInit " ++show size  ++
+                                          "-- only showing the first 100 elements\n"
+                                          ++ "(Left "++(show $ runST (Prelude.mapM f [0..100]) ) ++ "))\n"
+                             | otherwise ="(BatchInit " ++show size
+                                          ++ "(Left "++(show $ runST (Prelude.mapM f [0,1..size -1]) ) ++ "))\n"
+
+
+newtype IntFun a = IntFun  (forall m. (PrimMonad m)=>  Int -> m a )
+-- This may change substantially in a future release, but for now
+-- acts like
+            deriving (Typeable)
+
+instance  Functor IntFun  where
+  fmap f (IntFun g) = IntFun (\x->   g x  >>= (\ y -> return (f y))  )
+  {-# INLINE fmap #-}
+
+instance Functor BatchInit  where
+  {-# INLINE fmap  #-}
+  fmap = \f bival ->
+              case  bival of
+                (BatchInit size (Left ls))->
+                       BatchInit size (Left (Prelude.map   f  ls  ))
+                (BatchInit size (Right gfun))->
+                      BatchInit size (Right  $ fmap  f gfun  )
+
+
+
+-- batchInit size should be Word rather than Int for size, but Vector is lame like that
+
+
+{-
+ChoiceT from monad lib is tempting
+as is one of the ListT done right
+Bundle from Vector 0.11 and Stream from 0.10 are both alluring too
+
+but all of them make things complicated,
+punt for now
+
+
+ALso: I may want/need to distinguish sparse vs dense builders
+and put them into different classes, punting that for now
+-}
+
+
+fromListBI :: [a] -> BatchInit a
+fromListBI ls = BatchInit (length ls) (Left ls)
+
+fromVectorBI :: VG.Vector v e => v e -> BatchInit e
+fromVectorBI v =  BatchInit size
+      (Right
+        (IntFun $
+          \i -> if i >= size
+              then error  $ " out of bounds index on IntFun of size: " ++ show i
+              else return $ v VG.! i
+            ))
+  where
+    size = VG.length v
+
+fromMVectorBI :: (VGM.MVector mv e ) => AnyMV mv e -> BatchInit e
+fromMVectorBI (AMV v) =  BatchInit size
+      (Right
+        (IntFun $
+          \i -> if i >= size
+              then error  $ " out of bounds index on IntFun of size: " ++ show i
+              else  v `VGM.read` i
+            ))
+  where
+    size = VGM.length v
+
+
+
+class Layout form (rank::Nat) => LayoutBuilder form (rank::Nat) | form -> rank where
+
+  buildFormatM :: (store~FormatStorageRep form,Buffer store Int ,Buffer store a,PrimMonad m)=>
+         Index rank  -> proxy form -> a
+         -> Maybe (BatchInit  (Index rank ,a))
+         ->m (form, BufferMut store (PrimState m) a )
+
+
+buildFormatPure:: forall store form rank proxy m  a. (LayoutBuilder form (rank::Nat)
+  ,store~FormatStorageRep form,Buffer store Int  ,Buffer store  a, Monad m ) =>
+     Index rank -> proxy form -> a  -> Maybe (BatchInit  (Index rank ,a))
+                                              ->m (form, BufferPure store  a )
+buildFormatPure shape prox defaultValue builder =
+  do  res@(!_,!_)<-return $! theComputation
+      return res
+  where
+        theComputation :: (form,BufferPure store   a )
+        !theComputation = runST $
+            do  (form,buf) <- buildFormatM shape prox defaultValue builder
+                pureBuff <- VG.unsafeFreeze buf
+                return (form, pureBuff)
+{-
+this is a funky api for both dense and sparse arrays general builder format.
+
+given the target shape, logical dimensions,a default value (only used for dense arrays)
+and the list of manifest values (mostly only used for sparse), build the format
+descriptor and the suitably initialized and sized values buffer
+
+this api is only meant for internal use for building new array values
+
+
+TODO: compare using a catenable priority heap vs just doing fast sorting.
+-}
+
+
+{-
+the dense instances ignore the builder structure, which does suggest that maybe
+there shoudl be a dense builder layout class and a sparse layout class separately
+-}
+
+instance LayoutBuilder (Format  Direct 'Contiguous ('S 'Z) rep) ('S 'Z) where
+
+   buildFormatM (size:* _) _ defaultValue _ =
+      do
+        buf<-  VGM.replicate size defaultValue
+        return (FormatDirectContiguous  size,buf)
+
+
+-- really wish I didn't have to write the foldable and traversable constraints
+-- seems like a code smell?!
+instance (F.Foldable (Shape r),T.Traversable (Shape r) ,A.Applicative (Shape r))
+  => LayoutBuilder (Format  Row 'Contiguous r rep) r  where
+
+   buildFormatM ix  _ defaultValue _ =
+      do
+        buf<-  VGM.replicate (F.foldl' (*) 0   ix) defaultValue
+        return (FormatRowContiguous   ix,buf)
+
+instance (F.Foldable (Shape r),T.Traversable (Shape r) ,A.Applicative (Shape r))
+  =>  LayoutBuilder (Format  Column 'Contiguous r rep) r  where
+
+   buildFormatM ix  _ defaultValue _ =
+      do
+        buf<-  VGM.replicate (F.foldl' (*) 0   ix) defaultValue
+        return (FormatColumnContiguous   ix,buf)
+
+isStrictlyMonotonicV ::(VG.Vector v e)=> (e -> e->Ordering)-> v e -> Maybe Int
+isStrictlyMonotonicV cmp v = go  0 (VG.length v)
+  where
+    go !i !len  | i+1 >= len   = Nothing
+              |  (v VG.! i) `lt` (v VG.! (i+1))= go (i+1) len
+             | otherwise = Just i
+
+    lt a b = case cmp a b  of
+                  LT -> True
+                  _ -> False
+
+
+instance (Buffer rep Int)=>LayoutBuilder (Format DirectSparse 'Contiguous ('S 'Z) rep ) ('S 'Z) where
+
+
+  buildFormatM (size:* _) _ _ Nothing  = do
+      mvI <- VGM.new 0
+      vI <- VG.unsafeFreeze mvI
+      mvV <- VGM.new 0
+      return $!  (FormatDirectSparseContiguous size 0 vI, mvV)
+
+  buildFormatM (size:* _) _ _ (Just builder)= do
+    -- need to use let so type inference doesnt totally barf
+    mvt@(MVPair (MVLeaf ix) (MVLeaf val)) <- materializeBatchMV $ fmap  ( \((ix:*_),v)-> (ix,v)) builder
+    -- if i swap to using this  instead of  ix <- mat.. ; val <- mat..
+    --i get CRAZY type errors
+    -- could this be a ghc bug?
+    --ix <- materializeBatchMV $ fmap fst builtTup
+    --val <- materializeBatchMV $ fmap snd builtTup
+    _<- IntroSort.sortBy  (\x y -> compare (fst x) (fst y)) mvt
+                                                              -- this lets me sort a pair of arrays!
+    vIx <- VG.unsafeFreeze ix
+    optFail  <- return $ isStrictlyMonotonicV   compare vIx
+    --_hoelly
+    case optFail of
+      Nothing -> return (FormatDirectSparseContiguous size 0 vIx, val)
+      Just ixWrong ->  error $ "DirectSparse Index duplication at index "++ show (vIx VG.! ixWrong)
+
+
+instance (Buffer rep Int) => LayoutBuilder (Format CompressedSparseRow 'Contiguous ('S ('S 'Z)) rep ) ('S ('S 'Z)) where
+
+  buildFormatM (x:* y :* _) _   _ Nothing= do
+    mvi <-  VGM.new 0
+    vi <-  VG.unsafeFreeze  mvi
+    mvval <- VGM.new 0
+    return $
+      (FormatContiguousCompressedSparseRow
+              (FormatContiguousCompressedSparseInternal y x  vi vi), mvval )
+
+  buildFormatM (x:* y :* _) proxyFormat  _ (Just builder) = do
+    mvtup@(MVPair (MVPair (MVLeaf mvectYs) (MVLeaf mvectXs)) (MVLeaf mvectVals))<-
+          materializeBatchMV  $ fmap (\((xix:* yix :* _),val)-> ((yix,xix),val) ) builder
+    _ <-  IntroSort.sortBy (\((y1,x1),_) ((y2,x2),_) ->  basicCompareIndex  proxyFormat (x1:*y1 :* Nil) (x2:*y2:* Nil)  )
+                  mvtup
+    vectXs <- unsafeBufferFreeze mvectXs
+    vectYs <- unsafeBufferFreeze mvectYs
+    --- predicate check here wrt monotonicity of
+    --- compute runlength partial sums of where ys go
+
+    -- need to actually check
+    yRunsVect <- return $
+          VG.replicate (y+1) (0::Int) VG.//  computeStarts  (computeRunLengths vectYs) 0 y
+    --_ <- (error "computeRUnCount") vectYs yRunsMVect
+    --yRunsVect <- unsafeBufferFreeze yRunsMVect
+    let xyVect =         (VPair (VLeaf vectXs) (VLeaf vectYs))
+    optFail <- return $
+      isStrictlyMonotonicV (\(x1,y1) (x2,y2)->basicCompareIndex proxyFormat (x1:*y1:*Nil) (x2:*y2:*Nil))
+        xyVect
+    case optFail of
+      Nothing ->  return $
+        (FormatContiguousCompressedSparseRow
+              (FormatContiguousCompressedSparseInternal y x  vectXs yRunsVect), mvectVals )
+      Just i ->
+        error  $ "illegal duplication in CSR builder (x,y) coordinates  "
+            ++ show (xyVect VG.! i) ++ " and " ++ show (xyVect VG.! (i+1))
+            ++ "starting at position "  ++ show i
+
+
+computeRunLengths :: (VG.Vector v e, Eq e)=> v e -> [(e,Int)]
+computeRunLengths =  \y ->   fmap   (\x ->(head x,length x)) $ group $ VG.toList y
+
+
+
+{-# SPECIALIZE INLINE  computeStarts :: [(Int,Int)]->Int->Int ->[(Int,Int)] #-}
+
+computeStarts:: (Enum a, Ord a, Num b )=>[(a,b)]-> a -> a -> [(a,b)]
+computeStarts [] start end | start <= end  = fmap (\x -> (x ,0)) [start..end]
+                          |  otherwise = error "bad start end arguments to computeStarts"
+computeStarts ls start end | start <= end  = go start 0 ls
+                            | otherwise =  error "bad start end arguments to computeStarts"
+  where
+    --go :: a ->b->[(a,b)]-> [(a,b)]
+    go !posNext preSum [] | posNext <= end = fmap (\x -> (x,preSum)) [posNext .. end]
+                      | otherwise  = error "impossible go computeStarts "
+    go !posNext !preSum gls@((posAt,atSum):rest)
+            | posNext < posAt= (posNext,preSum):  go (succ posNext) preSum gls
+            | posNext == posAt = (posNext,preSum) : go (succ posNext) (preSum + atSum) rest
+            | otherwise = error "bad position in prefix stream for computeStarts go, literally unpossible "
+
+
+
+--computeStarts :: (Eq a, Num a)=> [(a,Int)]->Int -> [(a,Int)]
+--computeStarts [] len = map (\x -> (x ,0)) [0..len]
+--computeStarts ls len = go 0 0 ls
+--   where
+--    go preSum place [] |  place > len = []
+--                        | place == len = [(place,preSum)]
+--                        | otherwise = map
diff --git a/src/Numerical/Array/Layout/Dense.hs b/src/Numerical/Array/Layout/Dense.hs
new file mode 100644
--- /dev/null
+++ b/src/Numerical/Array/Layout/Dense.hs
@@ -0,0 +1,914 @@
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE DeriveDataTypeable #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE NoImplicitPrelude #-}
+{-# LANGUAGE FunctionalDependencies #-}
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE StandaloneDeriving #-}
+#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 707
+ {-# LANGUAGE AutoDeriveTypeable #-}
+#endif
+
+module Numerical.Array.Layout.Dense(
+  DenseLayout(..)
+  ,Locality(..)
+  ,Format(..)
+  ,Row
+  ,Column
+  ,Direct
+  ,module Numerical.Array.Layout.Base
+   ) where
+
+
+
+import Numerical.Nat
+import Control.Applicative
+import Numerical.Array.Locality
+import Numerical.Array.Layout.Base
+import Numerical.Array.Shape as S
+import Data.Data(Data,Typeable)
+
+
+--import Data.Traversable (Traversable)
+
+import Control.NumericalMonad.State.Strict
+
+import qualified Data.Foldable as F
+import Data.Traversable
+
+import Prelude hiding (foldr,foldl,map,scanl,scanr,scanl1,scanr1)
+
+
+data Direct
+
+data Row
+
+
+data Column
+
+
+{-
+one important gotcha about shape is that for many formats,
+the Shape is the (fmap (+1)) of the largestIndex,
+often, but perhaps not always.
+
+-}
+
+{-
+
+need to figure out how to support symmetric and hermitian and triangular
+and banded matrices
+
+-}
+
+
+--class Layout form rank => DenseLayout  form  (rank :: Nat) | form -> rank  where
+  {-
+  empty class instances for all the dense Layouts
+  -}
+
+type instance LayoutLogicalFormat (Format Direct  cont ('S 'Z ) rep  )
+     = Format Direct 'Contiguous ('S 'Z) rep
+
+-- | @'Format' 'Direct' 'Contiguous' ('S' 'Z')@ is a 1dim array 'Layout' with unit stride
+data instance Format  Direct 'Contiguous ('S 'Z) rep  =
+    FormatDirectContiguous {
+        logicalShapeDirectContiguous :: {-#UNPACK#-} !Int }
+    deriving (Show,Eq,Data)
+
+-- | @'Format' 'Direct' 'Strided'  ('S' 'Z')@ is a 1dim array 'Layout' with a regular stride >= 1
+data instance Format  Direct 'Strided ('S 'Z) rep  =
+    FormatDirectStrided {
+        logicalShapeDirectStrided :: {-#UNPACK#-}!Int
+        ,logicalStrideDirectStrided:: {-#UNPACK#-}!Int}
+    --deriving (Show,Eq,Data)
+
+type instance LayoutLogicalFormat (Format Row  cont n rep  )
+     = Format Row 'Contiguous n rep
+
+-- | @'Format'  'Row'  'Contiguous' n@ is a rank n Array
+data instance  Format  Row  'Contiguous n rep   =
+    FormatRowContiguous {
+        boundsFormRow :: !(Shape n Int)}
+    --deriving (Show,Eq,Data)
+
+data instance  Format  Row  'Strided n rep  =
+    FormatRowStrided
+        {boundsFormRowStrided:: !(Shape n Int)
+        ,strideFormRowStrided:: !(Shape n Int)}
+    --deriving (Show,Eq,Data)
+
+data instance  Format  Row  'InnerContiguous n rep  =
+    FormatRowInnerContiguous {
+        boundsFormRowInnerContig :: !(Shape n Int)
+        ,strideFormRowInnerContig:: !(Shape n Int)}
+    --deriving (Show,Eq,Data)
+
+type instance LayoutLogicalFormat (Format Column  cont n rep  )
+     = Format Column 'Contiguous n rep
+
+data instance  Format  Column 'Contiguous n  rep =
+    FormatColumnContiguous {
+      boundsColumnContig :: !(Shape n Int)}
+    --deriving (Show,Eq,Data)
+--deriving instance (Data (Shape n Int),Typeable n,Typeable rep) =>
+
+data instance  Format Column 'InnerContiguous n rep  =
+    FormatColumnInnerContiguous {
+        boundsColumnInnerContig :: !(Shape n Int)
+        ,strideFormColumnInnerContig:: !(Shape n Int)
+      }
+
+deriving instance Show (Shape n Int) => Show (Format Column 'InnerContiguous n rep)
+deriving instance (Data (Shape n Int),Typeable n,Typeable rep) =>Data (Format Column 'InnerContiguous n rep)
+    --deriving (Show,Eq,Data)
+
+data instance  Format Column 'Strided n rep  =
+    FormatColumnStrided {
+      boundsColumnStrided :: !(Shape n Int)
+      ,strideFormColumnStrided:: !(Shape n Int)}
+deriving instance Show (Shape n Int) => Show (Format Column 'Strided n rep)
+--deriving instance (Eq (Shape n Int)) => Eq (Format Column Strided n rep)
+deriving instance (Data (Shape n Int),Typeable n,Typeable rep) => Data (Format Column 'Strided n rep)
+    --deriving (Show,Eq,Data)
+
+
+type instance Transposed (Format Direct 'Contiguous ('S 'Z) rep) =
+    Format Direct 'Contiguous ('S 'Z) rep
+type instance Transposed (Format Direct 'Strided ('S 'Z) rep ) =
+   Format Direct 'Strided ('S 'Z) rep
+
+type instance  Transposed (Format Row  'Contiguous rank rep) =
+  Format Column 'Contiguous rank rep
+type instance Transposed (Format Row  'InnerContiguous rank rep) =
+    Format Column  'InnerContiguous rank rep
+type instance  Transposed (Format Row  'Strided rank rep) =
+    Format Column  'Strided rank rep
+
+type instance Transposed (Format Column 'Contiguous rank rep)=
+    Format Row 'Contiguous rank rep
+type instance Transposed (Format Column 'InnerContiguous rank rep)=
+    Format Row  'InnerContiguous rank rep
+type instance  Transposed (Format Column  'Strided rank rep)=
+    Format Row  'Strided rank rep
+
+
+{-
+a bunch of routines used to give various Layout operations for
+array Formats that have  DenseLayout instance
+not exported or for human use
+-}
+{-# INLINE basicAddressRangeGeneric #-}
+basicAddressRangeGeneric ::
+  (Functor (Shape rank),Applicative (Shape rank),F.Foldable (Shape rank),
+            DenseLayout form rank, Address~LayoutAddress form)=> form -> Maybe (Range Address)
+basicAddressRangeGeneric = \ form ->
+  if  (fmap (flip (-) 1)$ basicLogicalShape form) `strictlyDominates`  pure 0
+    then Just $!
+       Range  (basicToDenseAddress form  $! pure 0)
+              (basicToDenseAddress form $!
+                  fmap (flip (-) 1) $! basicLogicalShape form)
+    else Nothing
+
+{-# INLINE basicToAddressDenseGeneric #-}
+basicToAddressDenseGeneric :: (Functor (Shape rank),Applicative (Shape rank),F.Foldable (Shape rank),
+    DenseLayout form rank,Address~LayoutAddress form) => form -> Shape rank Int -> Maybe Address
+basicToAddressDenseGeneric = \ form ix ->
+  if (fmap (flip (-) 1)$ basicLogicalShape form) `weaklyDominates`  ix
+    && ix `weaklyDominates` pure 0
+    then Just $ basicToDenseAddress form ix
+    else Nothing
+{-# INLINE basicToIndexDenseGeneric #-}
+basicToIndexDenseGeneric ::
+  (Functor (Shape rank),F.Foldable (Shape rank),
+    DenseLayout form rank,Address~LayoutAddress form) =>  form -> Address -> Shape rank Int
+basicToIndexDenseGeneric = \form addr ->
+  basicToDenseIndex form addr
+
+{-# INLINE basicNextAddressDenseGeneric #-}
+basicNextAddressDenseGeneric ::
+  (Functor (Shape rank),F.Foldable (Shape rank),
+    DenseLayout form rank,Address~LayoutAddress form) =>  form -> Address-> Maybe Address
+basicNextAddressDenseGeneric = \ form addy ->
+  case  basicAddressRange form of
+    Just  (Range lo hi ) ->  if addy >= lo && addy < hi
+        then Just $! basicNextDenseAddress form addy
+        else Nothing
+    Nothing -> Nothing
+
+{-# INLINE basicNextIndexDenseGeneric #-}
+basicNextIndexDenseGeneric :: (Functor (Shape rank),F.Foldable (Shape rank),Applicative (Shape rank),
+    DenseLayout form rank,Address~LayoutAddress form)  =>
+    form -> Shape rank Int -> Maybe Address ->Maybe (Shape rank Int,Address)
+basicNextIndexDenseGeneric = \form ix _  ->
+  if (fmap (flip (-) 1)$ basicLogicalShape form) `strictlyDominates`  ix
+      && ix `weaklyDominates` pure 0
+    then
+      Just $! basicNextDenseIndex form ix
+    else
+      Nothing
+
+
+
+{- | note that basicAffineAddressShiftGeneric may be suboptimal,
+need to investigate what the core looks like
+also TODO needs tests
+-}
+{-# INLINE basicAffineAddressShiftDenseGeneric #-}
+basicAffineAddressShiftDenseGeneric :: (DenseLayout form rank
+  ,DenseLayout (LayoutLogicalFormat form) rank
+  ,Address~ LayoutAddress (LayoutLogicalFormat form))
+  => form -> Address -> Int -> Maybe Address
+basicAffineAddressShiftDenseGeneric form  = \ addy shift ->
+  let newForm = basicLogicalForm form in
+   do
+    nativeIndex <- return $ basicToDenseIndex form addy
+    popBaseAddress <- return $   basicToDenseAddress newForm nativeIndex
+    rng <- basicAddressRange newForm
+    candidateAddress <- return $ popBaseAddress + Address shift
+    if (getConst $ rangeMin ( Const) rng) <= candidateAddress
+        && candidateAddress  <= (getConst $ rangeMax ( Const) rng)
+      then return $ basicToDenseAddress  form $ basicToDenseIndex newForm candidateAddress
+      else Nothing
+
+
+
+
+-----
+-----
+-----
+
+
+type instance LayoutAddress (Format Direct 'Contiguous ('S 'Z) rep) = Address
+type instance LayoutLogicalFormat  (Format Direct 'Contiguous ('S 'Z) rep) = Format Direct 'Contiguous ('S 'Z) rep
+instance Layout (Format Direct 'Contiguous ('S 'Z) rep)  ('S 'Z)  where
+
+
+    {-# INLINE basicLogicalShape #-}
+    basicLogicalShape = \ x -> (logicalShapeDirectContiguous x) :* Nil
+
+    basicLogicalForm = id
+
+    transposedLayout = id
+
+    {-# INLINE basicCompareIndex #-}
+    basicCompareIndex = \ _  (l:* _) (r:* _) -> compare l r
+
+    basicAddressRange =  basicAddressRangeGeneric
+
+    basicToAddress = basicToAddressDenseGeneric
+
+    basicToIndex = basicToIndexDenseGeneric
+
+
+    basicNextAddress = basicNextAddressDenseGeneric
+
+    basicNextIndex = basicNextIndexDenseGeneric
+
+
+    basicAddressPopCount = \ _   (Range (Address lo) (Address hi )) ->
+      if  hi >= lo then hi - lo
+        else error $ "for basicAddressPopCount requires address obey hi >= lo, given: "
+          ++ show hi ++ " "  ++ show lo
+      -- FIX me, add the range error checking
+      -- in the style of the Sparse instances
+
+
+    basicAddressAsInt = \ _ (Address a) -> a
+
+    basicAffineAddressShift = basicAffineAddressShiftDenseGeneric
+
+    {-# INLINE basicAffineAddressShift #-}
+    {-# INLINE basicAddressRange #-}
+    {-# INLINE basicToAddress #-}
+    {-# INLINE basicToIndex #-}
+    {-# INLINE basicNextAddress #-}
+    {-# INLINE basicNextIndex #-}
+    {-# INLINE basicAddressPopCount #-}
+
+type instance LayoutAddress (Format Direct 'Strided ('S 'Z) rep) = Address
+
+instance  Layout (Format Direct 'Strided ('S 'Z) rep)  ('S 'Z)  where
+
+    {-# INLINE basicLogicalShape #-}
+    basicLogicalShape = \x -> (logicalShapeDirectStrided x) :* Nil
+
+    transposedLayout = id
+
+    basicLogicalForm =  (\  (n :* Nil ) ->  FormatDirectContiguous n) . basicLogicalShape
+
+    {-# INLINE basicCompareIndex #-}
+    basicCompareIndex = \ _  (l:* _) (r:* _) -> compare l r
+
+    basicAddressRange = basicAddressRangeGeneric
+
+    basicToAddress = basicToAddressDenseGeneric
+
+    basicToIndex = basicToIndexDenseGeneric
+
+    basicNextAddress =  basicNextAddressDenseGeneric
+
+    basicNextIndex =  basicNextIndexDenseGeneric
+
+    basicAddressPopCount = \form@(FormatDirectStrided size _ ) (Range loA hiA)->
+      let newForm = (FormatDirectContiguous size)
+        in
+          basicAddressPopCount  newForm
+            (Range (basicToDenseAddress newForm $ basicToDenseIndex form loA)
+                 (basicToDenseAddress newForm $ basicToDenseIndex form hiA) )
+
+    basicAddressAsInt = \ _ (Address a) -> a
+
+    basicAffineAddressShift = basicAffineAddressShiftDenseGeneric
+
+    {-# INLINE basicAffineAddressShift #-}
+    {-# INLINE basicAddressRange #-}
+    {-# INLINE basicToAddress #-}
+    {-# INLINE basicToIndex #-}
+    {-# INLINE basicNextAddress #-}
+    {-# INLINE basicNextIndex #-}
+    {-# INLINE basicAddressPopCount #-}
+
+
+-- one type family instance for all the rows
+type instance LayoutAddress (Format Row locality    rank rep) = Address
+
+instance   (Applicative (Shape rank), Traversable (Shape rank))
+    =>  Layout (Format Row  'Contiguous rank rep) rank where
+
+    transposedLayout = \(FormatRowContiguous shp) -> FormatColumnContiguous $ reverseShape shp
+
+    {-# INLINE basicLogicalShape #-}
+    basicLogicalShape =  boundsFormRow
+
+    basicLogicalForm = id
+
+    {-# INLINE basicCompareIndex #-}
+    basicCompareIndex = \ _  ls rs -> foldl majorCompareLeftToRight EQ  $ S.map2 compare ls rs
+
+    basicAddressPopCount = \ _   (Range (Address lo) (Address hi )) -> hi - lo
+      -- FIX me, add the range error checking
+      -- in the style of the Sparse instances
+    basicAddressRange = basicAddressRangeGeneric
+
+    basicToAddress = basicToAddressDenseGeneric
+
+    basicToIndex = basicToIndexDenseGeneric
+
+    basicNextAddress =  basicNextAddressDenseGeneric
+
+    basicNextIndex =  basicNextIndexDenseGeneric
+
+    basicAddressAsInt = \ _ (Address a) -> a
+
+    basicAffineAddressShift = basicAffineAddressShiftDenseGeneric
+
+    {-# INLINE basicAffineAddressShift #-}
+    {-# INLINE basicAddressRange #-}
+    {-# INLINE basicToAddress #-}
+    {-# INLINE basicToIndex #-}
+    {-# INLINE basicNextAddress #-}
+    {-# INLINE basicNextIndex #-}
+    {-# INLINE basicAddressPopCount #-}
+
+
+
+instance   (Applicative (Shape rank), Traversable (Shape rank))
+  =>  Layout (Format Row  'InnerContiguous rank rep)  rank  where
+
+    {-# INLINE basicLogicalShape  #-}
+    basicLogicalShape = boundsFormRowInnerContig
+
+    basicLogicalForm form = FormatRowContiguous $ basicLogicalShape form
+
+    transposedLayout = \(FormatRowInnerContiguous shp stride) ->
+        FormatColumnInnerContiguous  (reverseShape shp)  (reverseShape stride)
+
+    {-# INLINE basicCompareIndex #-}
+    basicCompareIndex = \ _  ls rs ->
+      foldl majorCompareLeftToRight EQ  $ S.map2 compare ls rs
+
+    basicAddressRange = basicAddressRangeGeneric
+
+    basicToAddress = basicToAddressDenseGeneric
+
+    basicToIndex = basicToIndexDenseGeneric
+
+    basicNextAddress =  basicNextAddressDenseGeneric
+
+    basicNextIndex =  basicNextIndexDenseGeneric
+
+    basicAddressPopCount = \form@(FormatRowInnerContiguous size _) (Range loA hiA)->
+      let newForm = (FormatRowContiguous size)
+        in
+          basicAddressPopCount  newForm
+            (Range (basicToDenseAddress newForm $ basicToDenseIndex form loA)
+                 (basicToDenseAddress newForm $ basicToDenseIndex form hiA) )
+
+    basicAddressAsInt = \ _ (Address a) -> a
+
+    basicAffineAddressShift = basicAffineAddressShiftDenseGeneric
+
+    {-# INLINE basicAffineAddressShift #-}
+    {-# INLINE basicAddressRange #-}
+    {-# INLINE basicToAddress #-}
+    {-# INLINE basicToIndex #-}
+    {-# INLINE basicNextAddress #-}
+    {-# INLINE basicNextIndex #-}
+    {-# INLINE basicAddressPopCount #-}
+
+
+
+instance  (Applicative (Shape rank),Traversable (Shape rank))
+  =>  Layout (Format Row 'Strided rank rep) rank  where
+
+    {-# INLINE basicLogicalShape  #-}
+    basicLogicalShape =  boundsFormRowStrided
+
+    basicLogicalForm form = FormatRowContiguous $ basicLogicalShape form
+
+    transposedLayout = \(FormatRowStrided shp stride) ->
+        FormatColumnStrided  (reverseShape shp)  (reverseShape stride)
+
+    {-# INLINE basicCompareIndex #-}
+    basicCompareIndex = \ _  ls rs ->
+        foldl majorCompareLeftToRight EQ  $ S.map2 compare ls rs
+
+    basicAddressRange = basicAddressRangeGeneric
+
+    basicToAddress = basicToAddressDenseGeneric
+
+    basicToIndex = basicToIndexDenseGeneric
+
+    basicNextAddress =  basicNextAddressDenseGeneric
+
+    basicNextIndex =  basicNextIndexDenseGeneric
+
+    basicAddressPopCount = \form@(FormatRowStrided size _) (Range loA hiA)->
+      let newForm = (FormatRowContiguous size)
+        in
+          basicAddressPopCount  newForm
+            (Range (basicToDenseAddress newForm $ basicToDenseIndex form loA)
+                 (basicToDenseAddress newForm $ basicToDenseIndex form hiA) )
+
+    basicAddressAsInt = \ _ (Address a) -> a
+
+    basicAffineAddressShift = basicAffineAddressShiftDenseGeneric
+
+    {-# INLINE basicAffineAddressShift #-}
+    {-# INLINE basicAddressRange #-}
+    {-# INLINE basicToAddress #-}
+    {-# INLINE basicToIndex #-}
+    {-# INLINE basicNextAddress #-}
+    {-# INLINE basicNextIndex #-}
+    {-# INLINE basicAddressPopCount #-}
+
+
+type instance LayoutAddress (Format Column locality    rank rep) = Address
+instance  (Applicative (Shape rank), Traversable (Shape rank))
+  =>  Layout (Format Column 'Contiguous rank rep)  rank where
+
+    {-# INLINE basicLogicalShape  #-}
+    basicLogicalShape =  boundsColumnContig
+
+    basicLogicalForm = id
+
+    transposedLayout = \(FormatColumnContiguous shp)-> FormatRowContiguous $ reverseShape shp
+
+    {-# INLINE basicCompareIndex #-}
+    basicCompareIndex = \ _  ls rs -> foldr majorCompareRightToLeft EQ  $ S.map2 compare ls rs
+
+    basicAddressPopCount = \ _   (Range (Address lo) (Address hi )) ->
+        if hi >= lo then hi - lo
+            else  error  $ "for basicAddressPopCount, require address hi >= lo, given: "
+              ++ show hi ++ " " ++ show lo
+      -- FIX me, add the range error checking
+      -- in the style of the Sparse instances
+    basicAddressRange = basicAddressRangeGeneric
+
+    basicToAddress = basicToAddressDenseGeneric
+
+    basicToIndex = basicToIndexDenseGeneric
+
+    basicNextAddress =  basicNextAddressDenseGeneric
+
+    basicNextIndex =  basicNextIndexDenseGeneric
+
+    basicAddressAsInt = \ _ (Address a) -> a
+
+
+    basicAffineAddressShift = basicAffineAddressShiftDenseGeneric
+
+    {-# INLINE basicAffineAddressShift #-}
+    {-# INLINE basicAddressRange #-}
+    {-# INLINE basicToAddress #-}
+    {-# INLINE basicToIndex #-}
+    {-# INLINE basicNextAddress #-}
+    {-# INLINE basicNextIndex #-}
+    {-# INLINE basicAddressPopCount #-}
+
+
+instance  (Applicative (Shape rank), Traversable (Shape rank))
+  => Layout (Format Column 'InnerContiguous rank rep) rank  where
+
+
+    {-# INLINE basicLogicalShape  #-}
+    basicLogicalShape =  boundsColumnInnerContig
+
+    basicLogicalForm form = FormatColumnContiguous $ basicLogicalShape form
+
+    transposedLayout = \(FormatColumnInnerContiguous shp stride)->
+         FormatRowInnerContiguous (reverseShape shp) (reverseShape stride)
+
+    {-# INLINE basicCompareIndex #-}
+    basicCompareIndex = \ _  ls rs -> foldr majorCompareRightToLeft EQ  $ S.map2 compare ls rs
+
+    basicAddressRange = basicAddressRangeGeneric
+
+    basicToAddress = basicToAddressDenseGeneric
+
+    basicToIndex = basicToIndexDenseGeneric
+
+    basicNextAddress=  basicNextAddressDenseGeneric
+
+    basicNextIndex=  basicNextIndexDenseGeneric
+
+    basicAddressPopCount = \form@(FormatColumnInnerContiguous size _) (Range loA hiA)->
+      let newForm = (FormatColumnContiguous size)
+        in
+          basicAddressPopCount  newForm
+            (Range (basicToDenseAddress newForm $ basicToDenseIndex form loA)
+                 (basicToDenseAddress newForm $ basicToDenseIndex form hiA) )
+
+    basicAddressAsInt = \ _ (Address a) -> a
+ -- strideRow :: Shape rank Int,
+
+    basicAffineAddressShift = basicAffineAddressShiftDenseGeneric
+
+    {-# INLINE basicAffineAddressShift #-}
+    {-# INLINE basicAddressRange #-}
+    {-# INLINE basicToAddress #-}
+    {-# INLINE basicToIndex #-}
+    {-# INLINE basicNextAddress #-}
+    {-# INLINE basicNextIndex #-}
+    {-# INLINE basicAddressPopCount #-}
+
+instance   (Applicative (Shape rank), Traversable (Shape rank))
+  => Layout (Format Column 'Strided rank rep) rank where
+
+    {-# INLINE basicLogicalShape  #-}
+    basicLogicalShape = boundsColumnStrided
+
+    basicLogicalForm form = FormatColumnContiguous $ basicLogicalShape form
+
+    transposedLayout = \(FormatColumnStrided shp stride)->
+         FormatRowStrided (reverseShape shp) (reverseShape stride)
+
+    {-# INLINE basicCompareIndex #-}
+    basicCompareIndex = \ _  ls rs -> foldr majorCompareRightToLeft EQ $ S.map2 compare ls rs
+
+    basicAddressRange = basicAddressRangeGeneric
+
+    basicToAddress = basicToAddressDenseGeneric
+
+    basicToIndex = basicToIndexDenseGeneric
+
+    basicNextAddress =  basicNextAddressDenseGeneric
+
+    basicNextIndex =  basicNextIndexDenseGeneric
+
+    basicAddressPopCount = \form@(FormatColumnStrided size _) (Range loA hiA)->
+      let newForm = (FormatColumnContiguous size)
+        in
+          basicAddressPopCount  newForm
+            (Range (basicToDenseAddress newForm $ basicToDenseIndex form loA)
+                 (basicToDenseAddress newForm $ basicToDenseIndex form hiA) )
+
+    basicAddressAsInt = \ _ (Address a) -> a
+
+    basicAffineAddressShift = basicAffineAddressShiftDenseGeneric
+
+    {-# INLINE basicAffineAddressShift #-}
+    {-# INLINE basicAddressRange #-}
+    {-# INLINE basicToAddress #-}
+    {-# INLINE basicToIndex #-}
+    {-# INLINE basicNextAddress #-}
+    {-# INLINE basicNextIndex #-}
+    {-# INLINE basicAddressPopCount #-}
+
+----------------------
+----------------------
+-----
+-----
+----------------------
+----------------------
+
+
+
+---
+---
+---
+
+{-
+these are factored out versions of the
+various shared computations in both Row and Column Major
+rank n Array format computations
+
+-}
+
+{-# INLINE computeStrideShape #-}
+computeStrideShape ::
+     ((Int -> State Int Int) -> Shape n Int  -> State Int (Shape n Int )) -> Shape n Int -> Shape n Int
+computeStrideShape = \trvse shp  ->
+    flip evalState 1 $
+                      flip  trvse shp  $
+                      -- basically accumulating the product of the
+                      -- dimensions
+                          \ val ->
+                               do accum <- get ;
+                                  put $! (val * accum) ;
+                                  return accum;
+
+
+
+
+-----
+-----
+-----
+
+instance DenseLayout (Format Direct 'Contiguous ('S 'Z) rep)  ('S 'Z)  where
+
+
+    --maxDenseAddress = \ (FormatDirectContiguous ix) -> Address (ix -1)
+
+
+    {-#INLINE basicToDenseAddress #-}
+    basicToDenseAddress   = \ (FormatDirectContiguous _) (j :* _ ) -> Address j
+
+    --basicNextIndex=  undefined -- \ _ x ->  Just $! x + 1
+    --note its unchecked!
+    {-# INLINE basicToDenseIndex #-}
+    basicToDenseIndex =  \ (FormatDirectContiguous _) (Address ix)  -> (ix ) :* Nil
+
+    {-# INLINE basicNextDenseAddress #-}
+    basicNextDenseAddress = \ _ addr -> addr + 1
+
+
+
+
+
+instance DenseLayout (Format Direct 'Strided ('S 'Z) rep)  ('S 'Z)  where
+
+
+
+
+    {-#INLINE basicToDenseAddress #-}
+    basicToDenseAddress   = \ (FormatDirectStrided _ strid) (j :* Nil )->  Address (strid * j)
+
+    {-# INLINE basicNextDenseAddress #-}
+    basicNextDenseAddress = \ (FormatDirectStrided _ strid) addr ->  addr + Address strid
+
+    {-# INLINE basicNextDenseIndex #-}
+    basicNextDenseIndex =  \ form  (i:* Nil ) ->  (\ix -> (ix,basicToDenseAddress form ix)) $! (i + 1 :* Nil )
+
+
+    {-# INLINE basicToDenseIndex #-}
+    basicToDenseIndex = \ (FormatDirectStrided _ stride) (Address ix)  -> (ix `div` stride ) :* Nil
+
+
+-----
+-----
+-----
+
+
+
+-- strideRow :: Shape rank Int,
+instance   (Applicative (Shape rank),F.Foldable (Shape rank), Traversable (Shape rank))  =>
+    DenseLayout (Format Row  'Contiguous rank rep) rank where
+
+{-
+TODO  AUDIT
+
+-}
+    {-# INLINE basicToDenseAddress #-}
+    --basicToAddress = \rs tup -> let !strider =takePrefix $! S.scanr (*) 1 (boundsFormRow rs)
+    basicToDenseAddress = \rs tup ->
+          let !strider =  computeStrideShape traverse (boundsFormRow rs)
+                  in Address $! S.foldl'  (+) 0 $! map2 (*) strider tup
+
+    {-# INLINE basicNextDenseAddress #-}
+    basicNextDenseAddress = \_ addr -> addr + 1
+
+    {-# INLINE basicToDenseIndex #-}
+    basicToDenseIndex  =   \ rs (Address ix) ->
+        let !striderShape  = computeStrideShape traverse (boundsFormRow rs)
+
+            in
+               flip evalState ix $
+                  flip (S.backwards traverse)  striderShape $
+                  -- want to start from largest stride (which is on the right)
+                      \ currentStride ->
+                             do remainderIx <- get ;
+                                let (!qt,!rm)= quotRem remainderIx currentStride
+                                put  $! rm
+                                return  qt;
+
+
+
+
+-----
+-----
+
+-- strideRow :: Shape rank Int,
+instance   (Applicative (Shape rank),F.Foldable (Shape rank), Traversable (Shape rank))
+  => DenseLayout (Format Row  'InnerContiguous rank rep) rank  where
+
+
+    {-# INLINE basicToDenseAddress #-}
+    basicToDenseAddress = \rs tup ->
+                       Address $! S.foldl'  (+) 0 $!
+                         map2 (*) (strideFormRowInnerContig rs ) tup
+
+    {-# INLINE basicNextDenseIndex #-}
+    basicNextDenseIndex = \ form@(FormatRowInnerContiguous shape _) ix ->
+        --S.map snd $!
+      (\index -> (index,basicToDenseAddress form  index)) $!
+        flip evalState 1 $
+           for   ((,) <$> ix <*> shape) $
+              \(ixv ,shpv   )->
+                  do  carry <-get
+                      let (newCarry,modVal)=divMod (carry + ixv) shpv
+                      put $! newCarry
+                      return modVal
+
+
+    {-# INLINE basicToDenseIndex #-}
+    basicToDenseIndex  =   \ rs (Address ix) ->   flip evalState ix $
+                          flip ( S.backwards traverse)  (strideFormRowInnerContig rs ) $
+                              \ currentStride ->
+                                     do remainderIx <- get ;
+                                        let (!qt,!rm)= quotRem remainderIx currentStride
+                                        put $! rm
+                                        return  qt;
+
+
+
+---
+---
+-- strideRow :: Shape rank Int,
+
+instance  (Applicative (Shape rank),F.Foldable (Shape rank), Traversable (Shape rank))
+  => DenseLayout (Format Row 'Strided rank rep) rank  where
+
+
+
+    {-# INLINE basicToDenseAddress #-}
+    basicToDenseAddress = \rs tup ->   Address $!
+          S.foldl'  (+) 0 $! map2 (*) (strideFormRowStrided rs ) tup
+
+    {-# INLINE basicNextDenseIndex #-}
+    basicNextDenseIndex = \ form@(FormatRowStrided shape _) ix ->
+      (\index -> (index,basicToDenseAddress form index)) $!
+        flip evalState 1 $
+           for  ((,) <$> ix <*> shape) $
+              \(ixv ,shpv   )->
+                  do  carry <-get
+                      let (newCarry,modVal)=divMod (carry + ixv) shpv
+                      put $! newCarry
+                      return modVal
+
+
+    {-# INLINE basicToDenseIndex #-}
+    basicToDenseIndex  =   \ rs (Address ix) ->   flip evalState ix $
+                          flip (S.backwards traverse ) (strideFormRowStrided rs ) $
+                              \ currentStride ->
+                                     do remainderIx <- get ;
+                                        let (!qt,!rm)= quotRem remainderIx currentStride
+                                        put $!  rm
+                                        return  qt;
+
+
+
+
+-----
+-----
+-----
+
+
+ -- strideRow :: Shape rank Int,
+instance  (Applicative (Shape rank),F.Foldable (Shape rank), Traversable (Shape rank))
+  => DenseLayout (Format Column  'Contiguous rank rep)  rank where
+
+
+
+    {-# INLINE basicToDenseAddress #-}
+    basicToDenseAddress = \rs tup ->
+          let !strider = computeStrideShape  (S.backwards traverse) (boundsColumnContig rs)
+                                in Address $! S.foldl'  (+) 0 $! map2 (*) strider tup
+
+    {-# INLINE basicNextDenseAddress #-}
+    basicNextDenseAddress = \_ addr -> addr + 1
+
+    {-# INLINE basicToDenseIndex #-}
+    basicToDenseIndex  =   \ rs (Address ix) ->
+            let !striderShape  =  computeStrideShape  (S.backwards traverse) (boundsColumnContig rs)
+                in
+                   flip evalState ix $
+                        for  striderShape $
+                              \ currentStride ->
+                                     do remainderIx <- get ;
+                                        let (!qt,!rm)= quotRem remainderIx currentStride
+                                        put $!  rm
+                                        return  qt;
+
+
+
+
+
+ -- strideRow :: Shape rank Int,
+instance  (Applicative (Shape rank),F.Foldable (Shape rank), Traversable (Shape rank))
+  => DenseLayout (Format Column  'InnerContiguous rank rep) rank  where
+
+
+    {-# INLINE basicToDenseAddress #-}
+    basicToDenseAddress    =   \ form tup -> let !strider =   strideFormColumnInnerContig form
+                                in Address $! foldl' (+) 0  $! map2 (*) strider tup
+    {-# INLINE basicNextDenseIndex #-}
+    basicNextDenseIndex = \ form@(FormatColumnInnerContiguous shape _) ix ->
+        --S.map snd $!
+      (\index -> (index,basicToDenseAddress form index)) $!
+        flip evalState 1 $
+           flip (S.backwards traverse)  ((,) <$> ix <*> shape) $
+              \(ixv ,shpv   )->
+                  do  carry <-get
+                      let (newCarry,modVal)=divMod (carry + ixv) shpv
+                      put $! newCarry
+                      return modVal
+
+
+    {-# INLINE basicToDenseIndex #-}
+    basicToDenseIndex  =   \ rs (Address ix) ->   flip evalState ix $
+                          flip S.traverse  (strideFormColumnInnerContig rs ) $
+                              \ currentStride ->
+                                     do remainderIx <- get ;
+                                        let (!qt,!rm)= quotRem remainderIx currentStride
+                                        put $! rm
+                                        return  qt;
+
+
+
+
+instance   (Applicative (Shape rank),F.Foldable (Shape rank), Traversable (Shape rank))
+  => DenseLayout (Format Column  'Strided rank rep) rank where
+
+    {-# INLINE basicToDenseAddress #-}
+    basicToDenseAddress    =   \ form tup -> let !strider =   strideFormColumnStrided form
+                                in Address $! foldl' (+) 0  $! map2 (*) strider tup
+
+    {-# INLINE basicNextDenseIndex #-}
+    basicNextDenseIndex = \ form@(FormatColumnStrided shape _) ix ->
+        --S.map snd $!
+      (\index -> (index,basicToDenseAddress form index)) $!
+        flip evalState 1 $
+           flip (S.backwards traverse)  ((,) <$> ix <*> shape) $
+              \(ixv ,shpv   )->
+                  do  carry <-get
+                      let (newCarry,modVal)=divMod (carry + ixv) shpv
+                      put $! newCarry
+                      return modVal
+
+
+    {-# INLINE basicToDenseIndex #-}
+    basicToDenseIndex  =   \ rs (Address ix) ->   flip evalState ix $
+                          flip S.traverse  (strideFormColumnStrided rs ) $
+                              \ currentStride ->
+                                     do remainderIx <- get ;
+                                        let (!qt,!rm)= quotRem remainderIx currentStride
+                                        put $! rm
+                                        return  qt;
+
+
+
+
+
+
+{-
+*Numerical.Array.Layout> basicToAddress (FormColumn (2 :* 3 :* 7 :* Nil)) (0:* 2 :* 2 :* Nil)
+Address 16
+*Numerical.Array.Layout> basicToAddress (FormColumn (2 :* 3 :* 7 :* Nil)) (1:* 0 :* 0 :* Nil)
+Address 1
+*Numerical.Array.Layout> basicToAddress (FormColumn (2 :* 3 :* 7 :* Nil)) (0:* 0 :* 0 :* Nil)
+Address 0
+*Numerical.Array.Layout> basicToAddress (FormColumn (2 :* 3 :* 7 :* Nil)) (0:* 1 :* 0 :* Nil)
+Address 2
+*Numerical.Array.Layout> basicToAddress (FormColumn (2 :* 3 :* 7 :* Nil)) (0:* 0 :* 1 :* Nil)
+-}
+
+
diff --git a/src/Numerical/Array/Layout/Sparse.hs b/src/Numerical/Array/Layout/Sparse.hs
new file mode 100644
--- /dev/null
+++ b/src/Numerical/Array/Layout/Sparse.hs
@@ -0,0 +1,913 @@
+{-
+the following (currently 5) sparse formats will live here
+
+
+DirectSparse 1dim
+
+
+
+one subtlety and a seemingly subtle point will be
+that contiguous / inner contiguous sparse arrays
+in  2dim  (and  1dim) will have an ``inner dimension" shift int.
+This is so that slices can  be zero copy on *BOTH* the array of values,
+and the Format indexing array machinery.
+
+Note that in the 2dim case, it still wont quite be zero copy, because the
+offsets into the inner dimension lookup table (not quite the right word)
+will have to change when a general slice is used rather than a slice
+that acts only on the outermost dimension.
+-}
+
+
+
+-- {-# LANGUAGE PolyKinds   #-}
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FunctionalDependencies #-}
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE StandaloneDeriving#-}
+{-# LANGUAGE FlexibleInstances  #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE NoImplicitPrelude #-}
+{-# LANGUAGE UndecidableInstances #-}
+
+#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 707
+ {-# LANGUAGE AutoDeriveTypeable #-}
+#endif
+module Numerical.Array.Layout.Sparse(
+  Layout(..)
+  ,DirectSparse
+  ,CSR
+  ,CSC
+  ,CompressedSparseRow
+  ,CompressedSparseColumn --  FIX ME, re add column support later
+  ,Format(FormatDirectSparseContiguous
+      ,FormatContiguousCompressedSparseRow
+      ,FormatInnerContiguousCompressedSparseRow
+      ,FormatContiguousCompressedSparseColumn
+      ,FormatInnerContiguousCompressedSparseColumn)
+  ,ContiguousCompressedSparseMatrix(..)
+  ,InnerContiguousCompressedSparseMatrix(..)
+  ,module Numerical.Array.Layout.Base
+  ) where
+
+import Data.Data
+import Data.Bits (unsafeShiftR)
+import Control.Applicative
+import Numerical.Array.Layout.Base
+--import Numerical.Array.Shape
+import Numerical.InternalUtils
+import qualified  Data.Vector.Generic as V
+import Prelude hiding (error )
+
+
+data CompressedSparseRow
+  deriving (Typeable)
+
+type CSR = CompressedSparseRow
+
+data CompressedSparseColumn
+    deriving (Typeable)
+
+type CSC = CompressedSparseColumn
+
+data DirectSparse
+    deriving (Typeable)
+
+
+
+data instance Format DirectSparse 'Contiguous ('S 'Z) rep =
+    FormatDirectSparseContiguous {
+      _logicalShapeDirectSparse:: {-# UNPACK#-} !Int
+      ,_logicalBaseIndexShiftDirectSparse::{-# UNPACK#-} !Int
+      ,_indexTableDirectSparse :: ! (BufferPure rep Int )  }
+
+
+deriving instance Show  (BufferPure rep Int )  => Show (Format DirectSparse 'Contiguous ('S 'Z) rep)
+
+
+
+{-
+for some listings of the design space of Sparse matrices
+as found in other tools,
+see < https://software.intel.com/en-us/mkl_11.2_ref >
+and then navigate to the section  "Sparse Matrix Storage Formats"  within
+"BLAS and Sparse BLAS Routines"
+
+<  http://netlib.org/linalg/html_templates/node90.html > is also pretty readable
+
+theres a subtle detail about the invariants of contiguous vs inner inner contiguous
+for CSR and CSC
+when I do an inner contiguous / contiguous slice / projection,
+what "address shifts" do i need to track to make sure the slices
+are zero copy as much as possible
+
+just slicing on the outer dimension doesn't need any row shifts,
+but a generalized (a,b) ... (a+x,b+y) selection when a,b!=0 does need a inner
+dim shift,
+
+NOTE that translating the inner dimension table's addresses to the corresponding
+value buffer's address can require a shift!
+This will happen when doing a MajorAxis (outer dimension) slice
+the picks out a Suffix of the CSR matrix's rows
+
+
+note that there are 2 formulations of CSR (/ CSC) formats
+
+a) 3 array: value, column index,  and  row start vectors
+
+b) 4 array: value, column index, rowstart, and row end vectors
+
+lets use choice a) for contiguous vectors, and choice b) for
+inner contiguous vectors.
+
+In both cases we need to enrich the type with a "buffer shift"
+to handle correctly doing lookups on submatrices picked out
+by either a major axis slice
+
+-}
+
+
+--deriving instance (Show (Shape (S (S Z)) Int), Show (BufferPure rep Int) )
+    -- => Show (Format CompressedSparseRow Contiguous (S (S Z)) rep)
+
+--deriving instance  (Eq (Shape (S (S Z)) Int), Eq (BufferPure rep Int) )
+    -- => Eq (Format CompressedSparseRow Contiguous (S (S Z)) rep)
+
+--deriving instance (Data (Shape (S (S Z)) Int), Data (BufferPure rep Int) )
+  --- => Data (Format CompressedSparseRow Contiguous (S (S Z)) rep)
+
+--deriving instance  (Typeable (Shape (S (S Z)) Int ), Typeable (BufferPure rep Int) )
+ -- => Typeable (Format CompressedSparseRow Contiguous (S (S Z)) rep)
+    --deriving (Eq,Data,Typeable)
+
+
+{-
+NOTE!!!!!
+_logicalBaseIndexShiftDirectSparse (and friends)
+are so that major axis slices can still use the same buffer,
+(needed for both Contiguous and InnerContiguous cases).
+So When looking up the Address for a value based upon its
+Inner dimension, we need to *SUBTRACT* that shift
+to get the correct offset index into the current SLICE.
+
+NB: THIS IS A TERRRIBLE EXPLANATION, FIXMEEEEE
+
+Phrased differently, This address shift is the *Discrepancy/Difference*
+between the size of the elided prefix of the Vector and the starting
+position of the manifest entries.
+
+(Q: does this ever ever matter, or can i punt that to vector, and only
+need this )
+
+
+This is kinda a good argument for not punting the Slicing on the raw buffers to
+Vector, because it generally makes this a bit more subtle to think about
+and someone IS going to implement something wrong this way!
+
+
+Another subtle and potentially confusing point is distinguishing between
+Affine shifts in the Index Space vs the Address space.
+
+Only the outer dimension lookup table shift is needed in the Contiguous
+2dim case, but the 2dim InnerContiguous case is a bit more confusing
+because of the potential for a slice along the inner dimension
+
+Rank 1 sparse  (like Direct sparse) is only Contiguous,
+and either a) doesn't need a shift, or b) only needs an index shift
+commensurate matching the leading implicit index of a Major Axis Slice
+
+
+theres a BIG corner case in most standard CSR / CSC formats which is
+underspecified in most docs about CSC and CSR formats.
+Consider Without loss of generality, CSR format
+  1) how are empty rows modeled/signaled?
+  2) if the last row is empty, how is that signaled?
+
+2) The last row is signaled to be be empty by having
+  the last entry of _outerDim2InnerDim buffer be set to >=
+  length of _innerDimIndex buffer (ie >= 1 + largest index of _innerDimIndex)
+1)
+
+note that the outer index table has 1+#rows length, with the last one being the
+length of the array
+
+-}
+
+data ContiguousCompressedSparseMatrix rep =
+    FormatContiguousCompressedSparseInternal {
+     -- does this need the index space shift for outer range slices???
+      _outerDimContiguousSparseFormat ::  {-# UNPACK #-} !Int
+      ,_innerDimContiguousSparseFormat ::  {-# UNPACK #-} !Int
+      ,_innerDimIndexContiguousSparseFormat :: !(BufferPure rep Int)
+      ,_outerDim2InnerDimContiguousSparseFormat:: ! (BufferPure rep Int )
+  }
+  deriving (Typeable)
+
+deriving instance (Show (BufferPure rep Int))=> Show (ContiguousCompressedSparseMatrix rep)
+
+{-
+  outerDim innerDim  innerTable  outer2InnerStart
+-}
+
+
+
+{-
+for Row major Compressed Sparse (CSR)
+the X dim (columns) are the inner dimension, and Y dim (rows) are the outer dim
+-}
+
+
+
+data  InnerContiguousCompressedSparseMatrix rep =
+   FormatInnerContiguousCompressedSparseInternal {
+      _outerDimInnerContiguousSparseFormat ::    {-# UNPACK #-} !Int
+      ,_innerDimInnerContiguousSparseFormat ::  {-# UNPACK #-} !Int
+      ,_innerDimIndexShiftInnerContiguousSparseFormat:: {-# UNPACK #-} !Int
+
+      ,_innerDimIndexInnerContiguousSparseFormat :: !(BufferPure rep Int)
+      ,_outerDim2InnerDimStartInnerContiguousSparseFormat:: ! (BufferPure rep Int )
+      ,_outerDim2InnerDimEndInnerContiguousSparseFormat:: ! (BufferPure rep Int )
+         }
+     deriving Typeable
+
+deriving instance (Show (BufferPure rep Int))=> Show (InnerContiguousCompressedSparseMatrix rep)
+
+
+newtype instance Format CompressedSparseRow 'Contiguous ('S ('S 'Z)) rep =
+    FormatContiguousCompressedSparseRow {
+      _getFormatContiguousCSR :: (ContiguousCompressedSparseMatrix rep) }
+
+deriving instance Show (ContiguousCompressedSparseMatrix rep)
+    => Show (Format CompressedSparseRow 'Contiguous ('S ('S 'Z)) rep)
+
+newtype instance Format CompressedSparseColumn 'Contiguous ('S ('S 'Z)) rep =
+    FormatContiguousCompressedSparseColumn {
+      _getFormatContiguousCSC :: (ContiguousCompressedSparseMatrix rep) }
+
+deriving instance Show (ContiguousCompressedSparseMatrix rep)
+    => Show (Format CompressedSparseColumn 'Contiguous ('S ('S 'Z)) rep)
+
+newtype instance Format CompressedSparseRow 'InnerContiguous ('S ('S 'Z)) rep =
+    FormatInnerContiguousCompressedSparseRow {
+      _getFormatInnerContiguousCSR :: (InnerContiguousCompressedSparseMatrix rep )
+  }
+deriving instance  Show (InnerContiguousCompressedSparseMatrix rep )
+    =>  Show (Format CompressedSparseRow 'InnerContiguous ('S ('S 'Z)) rep)
+
+newtype instance Format CompressedSparseColumn 'InnerContiguous ('S ('S 'Z)) rep =
+    FormatInnerContiguousCompressedSparseColumn {
+      _getFormatInnerContiguousCSC :: (InnerContiguousCompressedSparseMatrix rep )
+  }
+
+deriving instance Show (InnerContiguousCompressedSparseMatrix rep )
+  =>  Show (Format CompressedSparseColumn 'InnerContiguous ('S ('S 'Z)) rep)
+
+      --deriving (Show,Eq,Data)
+
+{-
+  FormatInnerContiguous rowsize columnsize
+
+-}
+--newtype instance Format CompressedSparseColumn Contiguous (S (S Z)) rep =
+--    FormatContiguousCompressedSparseColumn {
+--      _getFormatContiguousCSC ::  (ContiguousCompressedSparseMatrix rep)
+--  }
+    --deriving (Show,Eq,Data)
+
+--newtype  instance Format CompressedSparseColumn InnerContiguous (S (S Z)) rep =
+--    FormatInnerContiguousCompressedSparseColumn {
+--     _getFormatInnerContiguousCSC :: (InnerContiguousCompressedSparseMatrix rep)
+--  }
+--    --deriving (Show,Eq,Data)
+
+--CSR and CSC go here, and their version of lookups and next address and next index
+
+
+
+
+
+
+--  Offset binary search --- cribbed with permission from
+-- edward kmett's structured lib
+
+{-
+todo: theres some neat micro optimizations that are
+possible If I know how indexed structures are paged aligned and what not
+eg, when binary search, check both the first and last slot of a page I land on.
+Also on >= Nehalem, pages are "paired" so if you land on the lower page, the
+upper page is always loaded, etc etc. Not doing these for now.
+
+
+also should compare against search strategies defined in
+the vector-algorithms package, namely the
+galloping ones
+-}
+
+
+
+
+{-
+-- Assuming @l <= h@. Returns @h@ if the predicate is never @True@ over @[l..h)@
+-- requires p be a "monotonic" predicate  (FFFFFTTTTT)
+-}
+bsearchUp :: (Int -> Bool) -> Int -> Int -> Int
+bsearchUp p = go where
+  go l h
+    | l == h    = l
+    | p m       = go l m
+    | otherwise = go (m+1) h
+    where hml = h - l
+          m = l + unsafeShiftR hml 1 + unsafeShiftR hml 6
+{-# INLINE bsearchUp #-}
+{-
+ Assuming @l <= h@. Returns @l@ if the predicate is never @True@ over @(l..h]@
+  assumes predicate p is monotonic decreasing TTTTTFFFFF
+  -}
+bsearchDown :: (Int -> Bool) -> Int -> Int -> Int
+bsearchDown p = go where
+  go l h
+    | l == h    = l
+    | p (m+1)       = go (m+1) h
+    | otherwise = go l m
+    where hml = h - l
+          m = l + unsafeShiftR hml 1 + unsafeShiftR hml 6
+{-# INLINE bsearchDown #-}
+
+{-
+-- Assuming @l <= h@. Returns @h@ if the predicate is never @True@ over @[l..h)@
+-- requires p be a "monotonic" predicate  (FFFFFTTTTT)
+-- does a linear scan on the first constant number of slots
+(for now 97 because i had to pick a number thats ~ log MaxInt)
+and then falls back to binary search.
+Meant to have O(1) average case, O(log n) worst case
+-}
+basicHybridSearchUp :: (Int -> Bool ) -> Int -> Int -> Int
+basicHybridSearchUp  p = goCaseMe where
+  goCaseMe l h  | (h-l <= magicConstant) || p magicConstant
+                  {- either the range is short, OR
+                    we know match happens in the first magicConstant size subrange
+                  -}
+                    = linearSearchUp p l (min h magicConstant)
+                | otherwise = bsearchUp p magicConstant h
+{-# INLINE  basicHybridSearchUp #-}
+
+
+basicHybridSearchDown :: (Int -> Bool)-> Int -> Int -> Int
+basicHybridSearchDown  p = goCaseMe where
+  goCaseMe l h  | (h-l <= magicConstant)  || p (h- magicConstant)
+                {-  either the range is short, OR
+                    we know match happens in the first magicConstant size subrange
+                 -}
+                    = linearSearchDown p  (max l (h - magicConstant)) h
+                | otherwise = bsearchDown p l (h - magicConstant)
+{-# INLINE basicHybridSearchDown #-}
+
+{-
+i chose 97 because it seemed like a number thats ~ log MaxInt always (within 4x)
+And is a range that should stay in L1 cache sizes for most purposes
+-}
+magicConstant :: Int
+magicConstant = 97
+
+
+-- Assuming @l <= h@. Returns @h@ if the predicate is never @True@ over @[l..h)@
+linearSearchUp :: (Int -> Bool)-> Int -> Int -> Int
+linearSearchUp p = go where
+  go l h
+    | l ==h = l
+    | p l = l
+    | otherwise = go (l+1) h
+{-#INLINE linearSearchUp #-}
+
+-- Assuming @l <= h@. Returns @l@ if the predicate is never @True@ over @(l..h]@
+linearSearchDown :: (Int -> Bool)-> Int -> Int -> Int
+linearSearchDown p = go where
+  go l h
+    | l ==h = l
+    | p h = h
+    | otherwise = go l (h-1)
+{-#INLINE linearSearchDown #-}
+
+
+
+
+--
+-- now assumed each key is unique and ordered
+--
+-- Assuming @l <= h@. Returns @h@ if the predicate is never @True@ over @[l..h)@
+
+-- should at some point try out a ternary search scheme to have even better
+-- cache behavior (and benchmark of course)
+
+searchOrd :: (Int -> Ordering) -> Int -> Int -> Int
+searchOrd  p = go where
+  go l h
+    | l == h    = l
+    | otherwise = case p m of
+                  LT -> go (m+1) h
+                  ---  entry is less than target, go up!
+                  EQ -> m
+                  -- we're there! Finish early
+                  GT -> go l m
+                  -- entry is greater than target, go down!
+    where hml = h - l
+          m = l + unsafeShiftR hml 1 + unsafeShiftR hml 6
+{-# INLINE searchOrd #-}
+
+lookupExact :: (Ord k, V.Vector vec k) => vec k -> k -> Maybe Int
+lookupExact ks key
+  | j <- searchOrd (\i -> compare (ks V.! i)  key) 0 (V.length ks - 1)
+  , ks V.! j == key = Just $! j
+  | otherwise = Nothing
+{-# INLINE lookupExact #-}
+
+lookupExactRange :: (Ord k, V.Vector vec k) => vec k -> k -> Int -> Int -> Maybe Int
+lookupExactRange  ks key lo hi
+  | j <- searchOrd (\i -> compare (ks V.! i)  key) lo hi
+  , ks V.! j == key = Just $! j
+  | otherwise = Nothing
+{-# INLINE lookupExactRange  #-}
+
+--lookupLUB ::  (Ord k, V.Vector vec k) => vec k -> k -> Maybe Int
+--lookupLUB  ks key
+--  | j <- search  (\i -> compare (ks V.! i)  key) 0 (V.length ks - 1)
+--  , ks V.! j <= key = Just $! j
+--  | otherwise = Nothing
+--{-# INLINE lookupLUB  #-}
+
+type instance  Transposed (Format DirectSparse 'Contiguous ('S 'Z) rep )=
+   (Format DirectSparse 'Contiguous ('S 'Z) rep )
+
+
+
+
+type instance LayoutAddress (Format DirectSparse 'Contiguous ('S 'Z) rep) =  Address
+
+
+instance V.Vector (BufferPure rep) Int
+  => Layout  (Format DirectSparse 'Contiguous ('S 'Z) rep ) ('S 'Z) where
+
+  transposedLayout  = id
+  -- {-# INLINE transposedLayout #-}
+
+  basicLogicalShape = \ form -> _logicalShapeDirectSparse form  :* Nil
+  -- {-# INLINE basicLogicalShape #-}
+
+  basicCompareIndex = \ _ (a:* Nil) (b :* Nil) -> compare a b
+  -- {-# INLINE basicCompareIndex #-}
+
+  basicAddressRange = \form ->
+    case (minAddress form , maxAddress form ) of
+      (Just least, Just greatest) -> Just (Range least greatest )
+      _ -> Nothing
+
+    where
+        minAddress =
+          \ (FormatDirectSparseContiguous _ _   lookupTable) ->
+              if  V.length lookupTable >0 then  Just $! Address 0 else Nothing
+
+        maxAddress =
+          \ (FormatDirectSparseContiguous _ _   lookupTable) ->
+            if (V.length lookupTable >0 )
+               then Just $! Address (V.length lookupTable - 1 )
+               else Nothing
+
+-- TODO, double check that im doing shift correctly
+  {-# INLINE basicToAddress #-}
+  basicToAddress =
+      \ (FormatDirectSparseContiguous shape  indexshift lookupTable) (ix:*_) ->
+         if  not (ix < shape && ix > 0 ) then  Nothing
+          else  fmap Address  $! lookupExact lookupTable (ix + indexshift)
+
+  {-# INLINE basicToIndex #-}
+  basicToIndex =
+    \ (FormatDirectSparseContiguous _ shift lut) (Address addr) ->
+        ((lut V.! addr ) - shift) :* Nil
+  {-# INLINE basicAddressAsInt #-}
+  basicAddressAsInt = \ _ (Address a) -> a
+
+  {-# INLINE basicNextAddress #-}
+  basicNextAddress =
+    \ (FormatDirectSparseContiguous _ _ lut) (Address addr) ->
+      if  addr >= (V.length lut) then Nothing else Just  (Address (addr+1))
+
+  -- {-# INLINE basicAddressPopCount #-}
+  basicAddressPopCount = \ form (Range loadr@(Address lo) hiadr@(Address hi)) ->
+    if not ( lo <= hi ) then
+      error $! "basicAddressPopCount was passed a bad Address Range " ++ show loadr ++" " ++ show hiadr
+      else
+        case  basicAddressRange form of
+          Nothing -> 0
+          Just (Range (Address loBound) (Address  hiBound)) ->
+            if not $ (loBound<= lo ) && (hi <= hiBound)
+              then error $!
+               "basicAddressPopCount was passed a bad Address Range: "
+                ++show lo++" "++ show hi++"\nwith format Address range"
+                ++ show loBound ++ " " ++ show hiBound
+              else hi - lo
+
+
+{-
+    i've said it before, i'll say it again, scanning forward in the index space
+    for sparse structures is really weird, :)
+
+    NOTE: also need to remember to do those index space shifts for
+    1dim direct sparse, and test them thoroughly
+-}
+  -- {-# INLINE basicNextIndex #-}
+  basicNextIndex =
+    \form@(FormatDirectSparseContiguous size shift lut) (ix:*Nil) mebeAddress ->
+      if  ix >= size || ix >= (lut V.! (V.length lut -1) - shift ) then Nothing
+            -- if ix is out of bounds or the last element, we're done!
+      else
+        let
+            resAddr = Address $! bsearchUp  (\lix-> ix < ((lut V.! lix)-shift) )
+                        0 (V.length lut )
+        in case mebeAddress of
+          Nothing ->  resAddr `seq` (Just (basicToIndex form resAddr ,  resAddr))
+                -- Q: do i want the Index part of the tuple to be strict or not?
+                -- leaving it lazy for now
+                -- TODO / FIX / AUDIT ME / NOT SURE
+              -- this is the fall back binary search based lookup
+
+          Just (Address adr)->
+          -- make sure the address hint is in bounds and
+          -- is <= the current position
+              if adr >0 && adr < (V.length lut -1) && ix >=((lut V.! adr )-shift)
+              then
+                -- by construction we know theres at least one applicable index
+                -- thats
+                let !nextAddr = Address $!
+                                basicHybridSearchUp
+                                  (\lix-> ix <  ((lut V.! lix)-shift ) )
+                                  adr (V.length lut -1)
+                  in  Just (basicToIndex form nextAddr ,  nextAddr)
+              else
+                resAddr `seq` (Just (basicToIndex form resAddr ,  resAddr))
+
+
+------------
+------------
+
+type instance Transposed (Format CompressedSparseRow 'Contiguous ('S ('S 'Z)) rep )=
+    (Format CompressedSparseColumn 'Contiguous ('S ('S 'Z)) rep )
+
+
+type instance LayoutAddress (Format CompressedSparseRow 'Contiguous ('S ('S 'Z)) rep ) = SparseAddress
+
+instance  (V.Vector (BufferPure rep) Int )
+  => Layout  (Format CompressedSparseRow 'Contiguous ('S ('S 'Z)) rep ) ('S ('S 'Z)) where
+
+  transposedLayout  = \(FormatContiguousCompressedSparseRow repFormat) ->
+                          (FormatContiguousCompressedSparseColumn  repFormat)
+  {-# INLINE transposedLayout #-}
+
+
+  basicLogicalShape = \ form -> (_innerDimContiguousSparseFormat $ _getFormatContiguousCSR  form ) :*
+         ( _outerDimContiguousSparseFormat $ _getFormatContiguousCSR form ):* Nil
+          --   x_ix :* y_ix
+  {-# INLINE basicLogicalShape #-}
+
+
+  basicCompareIndex = \ _ as  bs -> shapeCompareRightToLeft as bs
+  {-# INLINE basicCompareIndex #-}
+
+
+  {-# INLINE basicAddressPopCount #-}
+  basicAddressPopCount = \ form (Range (SparseAddress _ lo) (SparseAddress _ hi)) ->
+    if not ( lo <= hi ) then
+      error $! "basicAddressPopCount was passed a bad Address Range " ++ show lo ++" " ++ show hi
+      else
+        case  basicAddressRange form of
+          Nothing -> 0
+          Just (Range (SparseAddress _ loBound) (SparseAddress _ hiBound)) ->
+            if not $ (loBound<= lo ) && (hi <= hiBound)
+              then error $!
+               "basicAddressPopCount was passed a bad SparseAddress Range: "
+                ++show lo++" "++ show hi++"\nwith format SparseAddress range"
+                ++ show loBound ++ " " ++ show hiBound
+              else hi - lo
+
+   -- {-# INLINE rangedFormatAddress #-}
+  basicAddressRange = \ form ->
+    case (minAddress form,maxAddress form) of
+      (Just least, Just greatest)-> Just (Range least greatest)
+      _ -> Nothing
+
+    where
+      {-
+      probably should deduplicate min/maxAddress
+      -}
+      minAddress =
+            \(FormatContiguousCompressedSparseRow
+                (FormatContiguousCompressedSparseInternal  y_row_range x_col_range
+                                                            columnIndex rowStartIndex)) ->
+                    if  y_row_range < 1  || x_col_range < 1|| (V.length columnIndex  < 1)
+                      then Nothing
+                      else
+                      -- the value buffer has the invariant the the end points
+                      -- of the buffer MUST be valid  in bounds values if length buffer > 0
+                    --SparseAddress $! 0 $! 0
+
+                    -- hoisted where into if branch as let so lets could be strict
+                        let
+                          !addrShift = columnIndex V.! 0
+
+                          -- for now assuming candidateRow is ALWAYS valid
+                          --- haven't proven this, FIXME
+                          !candidateRow= {-linearSearchUp-}
+                               basicHybridSearchUp nonZeroRow 0 (y_row_range-1 )
+
+
+                          {- FIXME, to get the right complexity
+                          to linear search on first log #rows + 1 slots, then fall
+                          back to binary search
+                          punting for now because this probably wont matter than often
+
+                          the solution will be to replace linearSearchUp
+                          with a hybridSearchUp
+                           -}
+                          nonZeroRow =
+                              \ !row_ix ->
+                                   -- the first row to satisfy this property
+                                  (rowStartIndex V.! (row_ix+1) >  rowStartIndex V.! row_ix)
+                                  -- if the start index is >0, already past the min address row!
+                                    ||  (rowStartIndex V.! row_ix) - addrShift > 0
+
+                                  --else  maxIxP1 >  rowStartIndex V.! row_ix
+                        in Just $! SparseAddress  candidateRow $! 0
+
+      maxAddress  =
+        \(FormatContiguousCompressedSparseRow
+            (FormatContiguousCompressedSparseInternal   y_row_range x_col_range
+                                                        columnIndex rowStartIndex)) ->
+                if  y_row_range < 1  || x_col_range < 1|| (V.length columnIndex  < 1)
+                  then Nothing
+                  else
+                  -- the value buffer has the invariant the the end points
+                  -- of the buffer MUST be valid  in bounds values if length buffer > 0
+                --SparseAddress $! 0 $! 0
+
+                -- hoisted where into if branch as let so lets could be strict
+                    let
+                      !addrShift = columnIndex V.! 0
+                      !maxIxP1 = V.length columnIndex
+
+                      -- for now assuming candidateRow is ALWAYS valid
+                      --- haven't proven this, FIXME
+                      !candidateRow= {-linearSearchDown-}
+                          basicHybridSearchDown nonZeroRow 0 (y_row_range-1 )
+
+
+                      {- FIXME, to get the right complexity
+                      to linear search on last log #rows + 1 slots, then fall
+                      back to binary search
+                      punting for now because this probably wont matter than often
+
+                      the solution will be to replace linearSearchDown
+                      with a hybridSearchDown
+                       -}
+                      nonZeroRow =
+                          \ !row_ix ->
+                       -- the first row to satisfy this property (going down from last row)
+                              (rowStartIndex V.! (row_ix+1) >  rowStartIndex V.! row_ix)
+                      -- if the start index is >= maxIxP1, havent gone down to max addres yet
+                      -- if < maxIxp1, we're at or below the max address
+                                ||  (rowStartIndex V.! row_ix) - addrShift < maxIxP1
+
+                              --else  maxIxP1 >  rowStartIndex V.! row_ix
+                    in
+                        Just $!
+                         SparseAddress  candidateRow $! (V.length columnIndex - 1 )
+
+       -- \ (FormatContiguousCompressedSparseRow
+       -- (FormatContiguousCompressedSparseInternal _ y_range
+       --          columnIndex _)) ->
+       --       SparseAddress (y_range - 1) (V.length columnIndex - 1 )
+
+  {-#  INLINE basicAddressAsInt #-}
+  basicAddressAsInt = \ _ (SparseAddress _ addr)-> addr
+
+  {-# INLINE basicToIndex #-}
+  basicToIndex =
+        \ (FormatContiguousCompressedSparseRow
+            (FormatContiguousCompressedSparseInternal  _ _ columnIndex _))
+            (SparseAddress outer inner) ->
+              (columnIndex V.! inner ) :* outer :*  Nil
+          -- outer is the row (y index) and inner is the lookup position for the x index
+
+
+{-
+theres 3 cases for contiguous next address:
+in the middle of a run on a fixed outer dimension,
+need to bump the outer dimension, or we're at the end of the entire array
+
+we make the VERY strong assumption that no illegal addresses are ever made!
+
+note that for very very small sparse matrices, the branching will have some
+overhead, but in general branch prediction should work out ok.
+-}
+  {-# INLINE basicNextAddress #-}
+  basicNextAddress =
+         \ (FormatContiguousCompressedSparseRow
+            (FormatContiguousCompressedSparseInternal  _ _
+              columnIndex rowStartIndex))
+            (SparseAddress outer inner) ->
+              if  inner < (V.length columnIndex -1)
+               -- can advance further
+                 -- && ( outer == (y_row_range-1)
+                  --- either last row
+                  || ((inner +1) < (rowStartIndex V.! (outer + 1)  - (rowStartIndex V.! 0 )))
+                     -- or our address is before the next row starts
+                     -- 3 vector CSR has a +1 slot at the end of the rowStartIndex
+
+                then
+                  Just (SparseAddress outer (inner+1))
+                else
+                  if inner == (V.length columnIndex -1)
+                    then Nothing
+                    else Just (SparseAddress (outer + 1) (inner + 1 ) )
+
+
+  -- {-# INLINE basicToAddress #-}
+  basicToAddress =
+        \ (FormatContiguousCompressedSparseRow
+            (FormatContiguousCompressedSparseInternal  y_row_range x_col_range
+              columnIndex rowStartIndex))
+          (ix_x:*ix_y :* _ ) ->
+            if  not (ix_x >= x_col_range ||  ix_y >=y_row_range )
+              then
+              -- slightly different logic when ix_y < range_y-1 vs == range_y-1
+              -- because contiguous, don't need the index space shift though!
+                let
+                  shift = (rowStartIndex V.! 0)
+                  checkIndex i =
+                      if  (columnIndex V.!i) == ix_x
+                        then Just i
+                        else Nothing
+                in
+                 (SparseAddress ix_y  <$>) $!
+                    checkIndex =<<
+                 --- FIXME  : need to check
+                      lookupExactRange columnIndex ix_x
+                          ((rowStartIndex V.! ix_y) - shift)
+                          ((rowStartIndex V.! (ix_y+1) ) - shift)
+
+              else   (Nothing :: Maybe SparseAddress )
+
+
+  -- {-# INLINE basicNextIndex #-}
+  {-  because nextIndex acts like a range query
+      it doesn't make sense for inner loops
+  -}
+  basicNextIndex =
+     \_form@(FormatContiguousCompressedSparseRow
+              (FormatContiguousCompressedSparseInternal
+                y_row_range x_col_range _columnIndex _rowStartIndex))
+      _ix@(innerX :* outerY :*Nil) mebeSparseAddress ->
+        if  not $ (innerX >=0 && innerX  < x_col_range ) && (outerY >= 0 && outerY < y_row_range)
+          -- checking if index is inbounds for logical shape
+          -- return Nothing if its out of bounds
+          -- QUESTION: should it throw an error instead of returning nothing?
+        then Nothing
+        else
+          case mebeSparseAddress of
+            Nothing -> error "finish me "
+              where
+              {- Okay here we check if the proposed current index is manifest, or not
+                Is it the right Row To search for the next index,
+                Or if We need to search further along. This is the way that
+                enables Usage of operations That give a complexity that is O(1)
+                in the average/best case and O(log N )in the worst case
+
+                The logical we do is roughly first check If there is an element
+                strictly Greater than ix in next we are doing the successor
+                That is within that Row And if so we can directly
+                  do a binary search therein
+                -}
+                _resRow = error "finish me "
+
+            (Just (SparseAddress _innerix _outerix) )
+                -> error "really finish me"
+
+
+        --case mebeAddress of
+        --  Nothing ->
+        --    let
+        --    resAddr = Address $! bsearchUp  (\lix-> ix < ((lut V.! lix)-shift) )
+        --                0 (V.length lut )
+        --  in
+        --   resAddr `seq` (Just (basicToIndex form resAddr ,  resAddr))
+        --        -- Q: do i want the Index part of the tuple to be strict or not?
+        --        -- leaving it lazy for now
+        --        -- TODO / FIX / AUDIT ME / NOT SURE
+        --      -- this is the fall back binary search based lookup
+        --  Just (Address adr)->
+        --  -- make sure the address hint is in bounds and
+        --  -- is <= the current position
+        --      if adr >0 && adr < (V.length lut -1) && ix >=((lut V.! adr )-shift)
+        --      then
+        --        -- by construction we know theres at least one applicable index
+        --        -- thats
+        --        let !nextAddr = Address $!
+        --                        basicHybridSearchUp
+        --                          (\lix-> ix <  ((lut V.! lix)-shift ) )
+        --                          adr (V.length lut -1)
+        --          in  Just (basicToIndex form nextAddr ,  nextAddr)
+        --      else
+        --        resAddr `seq` (Just (basicToIndex form resAddr ,  resAddr))
+
+
+
+
+--type instance Transposed (Format CompressedSparseRow InnerContiguous (S (S Z)) rep )=
+--    (Format CompressedSparseColumn InnerContiguous (S (S Z)) rep )
+
+--type instance Transposed (Format CompressedSparseColumn InnerContiguous (S (S Z)) rep )=
+--    (Format CompressedSparseRow InnerContiguous (S (S Z)) rep )
+
+
+--instance Layout (Format CompressedSparseRow InnerContiguous (S (S Z)) rep ) (S (S Z)) where
+--  transposedLayout  = \(FormatInnerContiguousCompressedSparseRow a b c d e f) ->
+--    (FormatInnerContiguousCompressedSparseColumn a b c d e f)
+--  {-# INLINE transposedLayout #-}
+--  basicFormShape = \ form -> logicalRowShapeInnerContiguousCSR form  :*
+--         logicalColumnShapeInnerContiguousCSR form :* Nil
+--  {-# INLINE basicFormShape #-}
+--  basicCompareIndex = \ _ as  bs ->shapeCompareRightToLeft as bs
+--  {-# INLINE basicCompareIndex#-}
+
+
+
+--instance  (V.Vector (BufferPure rep) Int )
+--  => SparseLayout (Format CompressedSparseRow InnerContiguous (S (S Z)) rep ) (S (S Z)) where
+
+--      type LayoutAddress (Format CompressedSparseRow
+--          InnerContiguous (S (S Z)) rep ) = SparseAddress
+
+--      {-# INLINE minSparseAddress #-}
+--      minSparseAddress = \_ -> SparseAddress 0 0
+
+--      {-# INLINE maxSparseAddress#-}
+--      maxSparseAddress  =
+--       \ (FormatInnerContiguousCompressedSparseInternal _ outer_dim_range _
+--          innerDimIndex _) ->
+--              SparseAddress (outer_dim_range - 1) (V.length innerDimIndex - 1 )
+
+
+--      {-#INLINE basicToIndex #-}
+--      basicToIndex =
+--       \ (FormatInnerContiguousCompressedSparseInternal _ _  _ innerDimIndex _)
+--          (SparseAddress outer inner) -> (innerDimIndex V.! inner ) :* outer :*  Nil
+--          -- outer is the row (y index) and inner is the lookup position for the x index
+
+
+
+--theres 3 cases for contiguous next address:
+--in the middle of a run on a fixed outer dimension,
+--need to bump the outer dimension, or we're at the end of the entire array
+
+--we make the VERY strong assumption that no illegal addresses are ever made!
+
+--note that for very very small sparse matrices, the branching will have some
+--overhead, but in general branch prediction should work out ok.
+
+--      {-# INLINE basicNextAddress #-}
+--      basicNextAddress =
+--         \ (FormatInnerContiguousCompressedSparseRow
+--                (FormatInnerContiguousCompressedSparseInternal _ _ _
+--                                                         columnIndex rowstartIndex))
+--            (SparseAddress outer inner) ->
+--              if not  (inner == (V.length columnIndex -1)
+--                                          {- && outer == (y_range-1) -}
+--                     || (inner +1) == (rowstartIndex V.! (outer + 1)))
+--                then
+--                  Just (SparseAddress outer (inner+1))
+--                else
+--                  if inner == (V.length columnIndex -1)
+--                    then Nothing
+--                    else Just (SparseAddress (outer + 1) (inner + 1 ) )
+
+--        --  error "finish me damn it"
+--      {-# INLINE basicToSparseAddress #-}
+--      basicToSparseAddress =
+--        \ (FormatInnerContiguousCompressedSparseRow
+--            (FormatInnerContiguousCompressedSparseInternal x_range y_range addrShift
+--                      columnIndex rowstartIndex))
+--          (ix_x:*ix_y :* _ ) ->
+--            if  not (ix_x >= x_range ||  ix_y >=y_range )
+--              then
+--              -- slightly different logic when ix_y < range_y-1 vs == range_y-1
+--              -- because contiguous, don't need the index space shift though!
+--                       SparseAddress ix_y   <$>
+--                          lookupExactRange columnIndex ix_x
+--                              -- ((rowstartIndex V.! ix_y) - addrShift)
+--                            (if  ix_y < (y_range-1)
+--                              -- addr shift is for correcting from a major axis slice
+--                              then  (rowstartIndex V.! (ix_y+1) ) - addrShift
+--                              else V.length columnIndex  - 1 )
+--              else   (Nothing :: Maybe SparseAddress )
diff --git a/src/Numerical/Array/Locality.hs b/src/Numerical/Array/Locality.hs
new file mode 100644
--- /dev/null
+++ b/src/Numerical/Array/Locality.hs
@@ -0,0 +1,81 @@
+
+
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE DeriveDataTypeable #-}
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE StandaloneDeriving #-}
+
+
+module Numerical.Array.Locality(Locality(..),LocalityMax,LocalityMin) where
+
+import Data.Data
+
+data  Locality = Contiguous | Strided  | InnerContiguous
+  deriving (Eq,Show,Read,Typeable,Data)
+
+#if defined(__GLASGOW_HASKELL__) && ( __GLASGOW_HASKELL__ >= 707) && ( __GLASGOW_HASKELL__ < 709)
+deriving instance Typeable 'Strided
+deriving instance Typeable 'InnerContiguous
+deriving instance  Typeable 'Contiguous
+#endif
+
+#if defined(__GLASGOW_HASKELL__) && ( __GLASGOW_HASKELL__ >= 707)
+type family LocalityMax (a :: Locality) (b :: Locality)  :: Locality where
+  LocalityMax 'Contiguous 'Contiguous = 'Contiguous
+  LocalityMax 'Contiguous  'InnerContiguous = 'Contiguous
+  LocalityMax 'Contiguous  'Strided = 'Contiguous
+  LocalityMax 'InnerContiguous  'Contiguous  = 'Contiguous
+  LocalityMax 'Strided  'Contiguous  = 'Contiguous
+  LocalityMax 'InnerContiguous  'InnerContiguous  = 'InnerContiguous
+  LocalityMax 'InnerContiguous  'Strided  = 'InnerContiguous
+  LocalityMax 'Strided 'InnerContiguous  = 'InnerContiguous
+  LocalityMax 'Strided 'Strided = 'Strided
+type family LocalityMin (a::Locality) (b ::Locality) :: Locality where
+  LocalityMin 'Contiguous 'Contiguous = 'Contiguous
+  LocalityMin 'Contiguous  'InnerContiguous = 'InnerContiguous
+  LocalityMin 'Contiguous  'Strided = 'Strided
+  LocalityMin 'InnerContiguous  'Contiguous  = 'InnerContiguous
+  LocalityMin 'Strided  'Contiguous  = 'Strided
+  LocalityMin 'InnerContiguous  'InnerContiguous  = 'InnerContiguous
+  LocalityMin 'InnerContiguous  'Strided  = 'Strided
+  LocalityMin 'Strided 'InnerContiguous  = 'Strided
+  LocalityMin 'Strided 'Strided = 'Strided
+
+#else
+type family LocalityMax (a :: Locality) (b :: Locality)  :: Locality
+type instance  LocalityMax  a b = LocalityMaxPrivate a b
+
+type family LocalityMaxPrivate (a :: Locality) (b :: Locality)  :: Locality
+type instance  LocalityMaxPrivate  'Contiguous 'Contiguous = 'Contiguous
+type instance  LocalityMaxPrivate 'Contiguous  'InnerContiguous = 'Contiguous
+type instance  LocalityMaxPrivate 'Contiguous  'Strided = 'Contiguous
+type instance  LocalityMaxPrivate 'InnerContiguous  'Contiguous  = 'Contiguous
+type instance  LocalityMaxPrivate 'Strided  'Contiguous  = 'Contiguous
+type instance  LocalityMaxPrivate 'InnerContiguous  'InnerContiguous  = 'InnerContiguous
+type instance  LocalityMaxPrivate 'InnerContiguous  'Strided  = 'InnerContiguous
+type instance  LocalityMaxPrivate 'Strided 'InnerContiguous  = 'InnerContiguous
+type instance  LocalityMaxPrivate 'Strided 'Strided = 'Strided
+
+type family LocalityMin (a::Locality) (b ::Locality) :: Locality
+type instance  LocalityMin a b = LocalityMinPrivate a b
+
+
+type family LocalityMinPrivate (a::Locality) (b ::Locality) :: Locality
+type instance  LocalityMinPrivate  'Contiguous 'Contiguous = 'Contiguous
+type instance  LocalityMinPrivate 'Contiguous  'InnerContiguous = 'InnerContiguous
+type instance  LocalityMinPrivate 'Contiguous  'Strided = 'Strided
+type instance  LocalityMinPrivate 'InnerContiguous  'Contiguous  = 'InnerContiguous
+type instance  LocalityMinPrivate 'Strided  'Contiguous  = 'Strided
+type instance  LocalityMinPrivate 'InnerContiguous  'InnerContiguous  = 'InnerContiguous
+type instance  LocalityMinPrivate 'InnerContiguous  'Strided  = 'Strided
+type instance  LocalityMinPrivate 'Strided 'InnerContiguous  = 'Strided
+type instance  LocalityMinPrivate 'Strided 'Strided = 'Strided
+
+
+#endif
+
+
+
+
+
diff --git a/src/Numerical/Array/Mutable.hs b/src/Numerical/Array/Mutable.hs
new file mode 100644
--- /dev/null
+++ b/src/Numerical/Array/Mutable.hs
@@ -0,0 +1,411 @@
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE ScopedTypeVariables#-}
+{-# LANGUAGE DeriveDataTypeable #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE ConstraintKinds #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE FunctionalDependencies #-}
+
+
+module Numerical.Array.Mutable(
+    MArray(..)
+    ,Array(..)
+    ,RectilinearArray(..)
+    ,DenseArrayBuilder(..)
+    ,DenseArray(..)
+    ,Boxed
+    ,Unboxed
+    ,Stored
+    --,module Numerical.Array.Layout
+    ,module Numerical.Array.Shape
+    ) where
+
+import Control.Monad.Primitive ( PrimMonad, PrimState )
+--import qualified Numerical.Array.DenseLayout as L
+import Numerical.Array.Address
+import qualified Numerical.Array.Layout as L
+import Numerical.Array.Layout (Layout,Locality(..),LayoutAddress,Format(..),Range(..),AffineRange(..))
+import Numerical.Array.Shape
+--import Numerical.Nat
+--import GHC.Prim(Constraint)
+import Numerical.World
+--import Numerical.Array.Range
+--import Numerical.Array.Storage(Boxed,Unboxed,Stored)
+--import Numerical.Array.Locality
+
+import qualified Numerical.Array.Pure as P
+import qualified Numerical.Array.Storage as S
+import Numerical.Array.Storage (Buffer,Boxed,Unboxed,Stored)
+import  qualified Data.Vector.Generic as VG
+import qualified  Data.Vector.Generic.Mutable as VGM
+
+import Control.Monad (liftM)
+--import qualified Data.Vector.Storable.Mutable as SM
+--import qualified Data.Vector.Unboxed.Mutable as UM
+--import qualified Data.Vector.Mutable as BM
+
+{-
+For now we're going to just crib the vector style api and Lift it
+up into a multi dimensional setting.
+
+the tentative design is to have something like
+
+
+
+you'd think that the following array type is ``right''
+but then you'll hit problems supporting
+-}
+
+-- data MArray world rep lay (view:: Locality) rank elm where
+--      MArray
+--          {_marrBuffer :: {-# UNPACK #!(MBuffer  world rep elm)
+--          ,_marrForm :: {-# UNPACK #- } !(Form lay loc rank)
+--          --,_marrShift :: {-# UNPACK #- } !Address
+--          }
+
+
+-- shift will be zero for most reps i'll ever care about, but in certain cases,
+-- might not be. So for now not including it, but might be needed later,
+-- though likely in regards to some sparse format of some sort.
+--Omitting it for now, but may need to revisit later!
+--
+--For now any 'Address' shift will need to be via the buffer
+--
+-- One ssue in the formats is ``logical'' vs ``manifest'' Address.
+--
+--
+--we eedto have 'RepConstraint' be decoupled from the type class instances
+-- because we to sometimes have things that are world parametric
+--
+-- indexing should be oblivious to locality,
+
+
+
+--NB: one important assumption we'll have for now, is that every
+
+
+-- dsfdf
+--type family RepConstraint world  rep el :: Constraint
+--type instance MArrayElem
+
+{- | 'MArray' is the generic data family that
+-}
+data family MArray world rep lay (view::Locality) (rank :: Nat ) st  el
+
+data instance  MArray Native rep lay locality rank st el =
+  MutableNativeArray {
+          nativeBuffer  :: ! (S.BufferMut rep st el  )
+          ,nativeFormat :: ! (Format lay locality rank rep)
+    }
+
+
+-- | Every 'MutableArray'  instance has a contiguous version
+-- of itself, This contiguous version will ALWAYS have a Builder instance.
+type family MutableArrayContiguous (marr :: * -> * -> *) :: * ->  * -> *
+type instance  MutableArrayContiguous (MArray world rep layout locality rank)= MArray world rep layout 'Contiguous rank
+
+-- | Sadly 'ArrMutable'  will have to have instances written by hand for now
+-- May later migrate the freeze / thaw machinery to Array.Phased, but lets
+type  family  ArrMutable ( arr :: * -> * )  :: * -> * -> *
+
+class P.PureArray (ArrPure marr)  rank a => Array marr (rank:: Nat)  a | marr -> rank  where
+
+    type   ArrPure (marr :: * -> * -> * ) :: * -> *
+
+    -- the type of the underlying storage buffer
+    --type MutableArrayBuffer marr :: * -> * -> *
+
+    -- really shouldnt appear in end user code, will only
+    -- come up in writing new combinators
+    -- the abstraction here is a reflection of the need for
+    type MArrayAddress (marr :: * -> * -> * ) ::  *
+
+    -- | 'basicUnsafeAffineAddressShift' is needed to handle abstracting access in popcount space
+    basicUnsafeAffineAddressShift :: (address ~ MArrayAddress marr) => marr st a -> Int -> address -> address
+  -- question, should the type be  -> address or  -> Maybe address
+
+    -- | Unsafely convert a mutable Array to its immutable version without copying.
+    -- The mutable Array may not be used after this operation. Assumed O(1) complexity
+    basicUnsafeFreeze :: (PrimMonad m, arr ~ ArrPure marr, marr ~ ArrMutable arr)
+        => marr (PrimState m) a -> m (arr a)
+
+    -- | Unsafely convert a pure Array to its mutable version without copying.
+    -- the pure array may not be used after this operation. Assumed O(1) complexity
+    basicUnsafeThaw :: (PrimMonad m, marr ~ ArrMutable arr, arr ~ ArrPure marr )
+        => arr a -> m (marr (PrimState m) a)
+
+    -- | gives the shape, a 'rank' length list of the dimensions
+    basicShape :: marr st    a -> Index rank
+
+    -- | 'basicCardinality' reports the number of manifest addresses/entries are
+    -- in the array in a given address sub range.
+    -- This is useful for determining when to switch from a recursive algorithm
+    -- to a direct algorithm.
+    -- Should this be renamed to something like basicPopCount/
+    basicCardinality ::(address ~ MArrayAddress marr) => marr st a -> Range address  -> Int
+
+    --basicUnsafeRead  :: PrimMonad m => marr  (PrimState m)   a -> Shape rank Int -> m (Maybe a)
+
+    --  | basicMutableSparseIndexToAddres checks if a index is present or not
+    -- helpful primitive for authoring codes for (un)structured sparse array format
+    basicSparseIndexToAddress :: (address ~ MArrayAddress marr)
+      => marr s   a -> Index rank  ->  Maybe address
+
+    -- | 'basicMutableAddressToIndex' assumes you only give it legal manifest addresses
+    basicAddressToIndex :: (address ~ MArrayAddress marr) =>marr s   a -> address ->    Index rank
+
+    -- |  return the smallest and largest valid logical address
+    basicAddressRange :: (address ~ MArrayAddress marr)=> marr st   a ->  Maybe (Range address)
+
+
+    -- | gives the next valid logical address
+    -- undefined on invalid addresses and the greatest valid address.
+    -- Note that for invalid addresses in between minAddress and maxAddress,
+    -- will return the next valid address.
+
+    basicSparseNextAddress :: (address ~ MArrayAddress marr)=> marr st  a -> address -> Maybe address
+
+
+    -- I think the case could be made for a basicPreviousAddress opeeration
+
+    -- | gives the next valid array index, the least valid index that is
+    -- or
+    basicSparseNextIndex ::(address ~ MArrayAddress marr)=>
+         marr st  a ->  Index rank -> Maybe address  -> Maybe ( Index rank, address)
+
+
+    -- | for a given valid address, @'basicAddressRegion' addr @ will return an AddressInterval
+    -- that contains @addr@. This will be a singleton when the "maximal uniform stride interval"
+    -- containing @addr@ has strictly less than 3 elements. Otherwise will return an Address range
+    -- covering the maximal interval that will have cardinality at least 3.
+    basicLocalAffineAddressRegion ::(address ~ MArrayAddress marr)
+          => marr st a ->address ->  AffineRange address
+
+    -- | this doesn't quite fit in this class, but thats ok, will deal with that later
+    basicOverlaps :: marr st   a -> marr st   a -> Bool
+
+    -- | Reset all elements of the vector to some undefined value, clearing all
+    -- references to external objects. This is usually a noop for unboxed
+    -- vectors. This method should not be called directly, use 'clear' instead.
+    basicClear :: PrimMonad m => marr (PrimState m)   a -> m ()
+
+    ---- | Yield the element at the given position. This method should not be
+    ---- called directly, use 'unsafeRead' instead.
+    basicUnsafeAddressRead  :: (PrimMonad m ,address ~ MArrayAddress marr) =>
+        marr  (PrimState m)   a -> address-> m a
+
+    ---- | Replace the element at the given position. This method should not be
+    ---- called directly, use 'unsafeAddressWrite' instead.
+    basicUnsafeAddressWrite :: (PrimMonad m ,address ~ MArrayAddress marr) =>
+         marr  (PrimState m)   a -> address  -> a -> m ()
+
+
+    --note  the sparsewrite and sparse read are "fused" versions of basicManifestAddress
+    -- and address read and write. probably needs to be benchmarked! TODO
+
+    -- | Yield the element at the given position. This method should not be
+    -- called directly, use 'unsafeSparseRead' instead.
+    basicUnsafeSparseRead :: PrimMonad m => marr  (PrimState m)   a ->
+       Index rank -> m (Maybe a)
+
+    --  Replace the element at the given position. This method should not be
+    -- called directly, use 'unsafeWrite' instead.
+    -- the following is the type that normal Array indexing,
+    -- as folks are used to, lookslike
+    -- its wrong
+    --basicUnsafeSparseWrite :: PrimMonad m => marr (PrimState m) a ->
+    --  Index rank -> m( Maybe (a -> m ()))
+-- this might get axed
+
+
+instance (Buffer rep el, Layout (Format  lay locality  rank rep) rank )
+  =>Array (MArray Native rep lay locality rank) rank el  where
+
+    type ArrPure (MArray Native rep lay locality rank)= P.ImmArray Native rep lay locality rank
+
+    type MArrayAddress (MArray Native rep lay locality rank)= LayoutAddress (Format  lay locality  rank rep)
+
+    {-# INLINE basicShape #-}
+    basicShape =  L.basicLogicalShape . nativeFormat
+
+    {-# NOINLINE basicUnsafeFreeze #-}
+    basicUnsafeFreeze = \marr -> do
+        pureBuffer <- VG.unsafeFreeze $ nativeBuffer marr
+        return $ P.ImMutableNativeArray pureBuffer $ nativeFormat marr
+
+    {-#  NOINLINE basicUnsafeThaw #-}
+    basicUnsafeThaw = \parr -> do
+        mutBuffer <- VG.unsafeThaw $ P.nativeBufferPure parr
+        return $ MutableNativeArray mutBuffer $ P.nativeFormatPure parr
+
+    {-# INLINE basicSparseIndexToAddress #-}
+    basicSparseIndexToAddress = \ marr  -> L.basicToAddress (nativeFormat marr)
+
+    {-# INLINE basicAddressToIndex #-}
+    basicAddressToIndex = \ marr  -> L.basicToIndex (nativeFormat marr)
+
+    {-# INLINE basicSparseNextAddress #-}
+    basicSparseNextAddress = \marr -> L.basicNextAddress (nativeFormat marr)
+
+    {-# INLINE basicSparseNextIndex #-}
+    basicSparseNextIndex = \marr -> L.basicNextIndex (nativeFormat marr)
+
+    basicOverlaps = \marr1 marr2 -> VGM.overlaps (nativeBuffer marr1) (nativeBuffer marr2)
+
+    basicClear = \marr -> VGM.clear (nativeBuffer marr)
+
+    {-# INLINE basicUnsafeAddressRead #-}
+    basicUnsafeAddressRead = \marr addr ->
+      VGM.unsafeRead (nativeBuffer marr) (L.basicAddressAsInt (nativeFormat marr) addr)
+
+    {-# INLINE basicUnsafeAddressWrite #-}
+    basicUnsafeAddressWrite = \marr addr v->
+      VGM.unsafeWrite (nativeBuffer marr) (L.basicAddressAsInt (nativeFormat marr) addr) v
+
+    {-# INLINE basicUnsafeSparseRead #-}
+    basicUnsafeSparseRead = \marr ix  ->  do
+      maddr <- return $ basicSparseIndexToAddress marr ix
+      maybe (return Nothing) (\addr -> liftM Just $  basicUnsafeAddressRead marr addr ) maddr
+
+    {-# INLINE basicAddressRange #-}
+    basicAddressRange = \marr -> L.basicAddressRange (nativeFormat marr)
+
+    basicCardinality = \marr -> L.basicAddressPopCount (nativeFormat marr)
+
+
+    basicUnsafeAffineAddressShift = error "carter needs to add this"
+    basicLocalAffineAddressRegion = error "crter needs to add this"
+{-
+
+
+type ArrPure marr :: * -> *
+
+type MArrayAddress marr :: *
+
+basicUnsafeAffineAddressShift :: (address ~ MArrayAddress marr) => marr st a -> Int -> address -> address
+
+basicLocalAffineAddressRegion :: (address ~ MArrayAddress marr) => marr st a -> address -> AffineRange address
+
+-}
+
+
+
+
+
+
+
+class ( Array marr rank a, P.PureDenseArray (ArrPure marr) rank a  )=>
+            DenseArray marr rank a | marr -> rank   where
+    -- | for Dense arrays, it is always easy to check if a given index is valid.
+    -- this operation better have  O(1) complexity or else!
+    basicIndexInBounds :: marr st a -> Index rank  -> Bool
+
+
+    --basicUnsafeAddressDenseRead  :: PrimMonad m => marr  (PrimState m)   a -> Address-> m a
+
+    -- i already have dense address indexing ?
+    --basicUnsafeAddressDenseWrite :: PrimMonad m => marr  (PrimState m)   a -> Address -> a -> m ()
+
+    -- | Yield the element at the given position. This method should not be
+    -- called directly, use 'unsafeRead' instead.
+    basicUnsafeDenseRead  :: PrimMonad m => marr  (PrimState m)   a -> Index rank -> m a
+
+    -- | Replace the element at the given position. This method should not be
+    -- called directly, use 'unsafeWrite' instead.
+    basicUnsafeDenseWrite :: PrimMonad m => marr (PrimState m)   a -> Index rank   -> a -> m ()
+
+
+    -- | gives the next valid logical address
+    -- undefined on invalid addresses and the greatest valid address.
+    -- Note that for invalid addresses in between minAddress and maxAddress,
+    -- will return the next valid address.
+
+    basicNextAddress ::  marr st  a -> Address ->  Address
+
+
+    -- I think the case could be made for a basicPreviousAddress opeeration
+
+    -- | gives the next valid array index
+    -- undefined on invalid indices and the greatest valid index
+    basicNextIndex :: marr st  a -> Index rank  -> Index rank
+
+
+
+
+
+
+
+{-
+
+Mutable (Dense) Array Builder will only have contiguous instances
+and only makes sense for dense arrays afaik
+
+BE VERY THOUGHTFUL about what instances you write, or i'll be mad
+
+
+not including the general sparse building in the first release,
+will include subsequently
+-}
+
+--class MutableArray marr (rank:: Nat) a => MutableArrayBuilder marr rank a where
+    --basicBuildArray:: Index rank -> b
+
+class DenseArray marr rank a => DenseArrayBuilder marr rank a where
+    basicUnsafeNew :: PrimMonad m => Index rank -> m (marr (PrimState m)   a)
+    basicUnsafeReplicate :: PrimMonad m => Index rank  -> a -> m (marr (PrimState m)  a)
+
+
+
+
+class RectilinearArray marr rank a | marr -> rank   where
+
+    -- | @'MutableRectilinearOrientation' marr@ should equal Row or Column for any sane choice
+    -- of instance, because every MutableRectilinear instance will have a notion of
+    -- what the nominal major axix will be.
+    -- The intended use case is side condition constraints like
+    -- @'MutableRectilinearOrientation' marr~Row)=> marr -> b @
+    -- for operations where majorAxix projections are correct only for Row
+    -- major formats. Such  as Row based forward/backward substitution (triangular solvers)
+    type MutableRectilinearOrientation marr :: *
+
+    type MutableArrayDownRank  marr ( st:: * ) a
+
+
+    -- | MutableInnerContigArray is the "meet" (minimum) of the locality level of marr and InnerContiguous.
+    -- Thus both Contiguous and InnerContiguous are made InnerContiguous, and Strided stays Strided
+    -- for now this makes sense to have in the MutableRectilinear class, though that may change.
+    -- This could also be thought of as being the GLB (greatest lower bound) on locality
+    type MutableInnerContigArray (marr :: * ->  * -> *)  st  a
+
+
+
+    --type MutableArrayBuffer
+    --not implementing this .. for now
+
+    -- | @'basicSliceMajorAxis' arr (x,y)@ returns the sub array of the same rank,
+    -- with the outermost (ie major axis) dimension of arr restricted to the
+    -- (x,y) is an inclusive interval, MUST satisfy x<y , and be a valid
+    -- subinterval of the major axis of arr.
+    basicMutableSliceMajorAxis :: PrimMonad m => marr (PrimState m)  a ->
+      (Int,Int)-> m (marr (PrimState m)  a)
+    --but  should it be primmonadic? nah, tis pure!
+
+    --  |  semantically, 'basicProjectMajorAxis' arr ix, is the rank reducing version of what
+    -- basicSliceMajorAxis arr (ix,ix) would mean _if_ the (ix,ix) tuple was a legal major axis slice
+    basicMutableProjectMajorAxis :: PrimMonad m =>marr (PrimState m)  a
+        -> Int -> m (MutableArrayDownRank marr (PrimState m)  a )
+
+    -- | @'basicMutableSlice' arr ix1 ix2@  picks out the (hyper) rectangle in dimension @rank@
+    -- where ix1 is the minimal corner and ix2
+    basicMutableSlice :: PrimMonad m => marr (PrimState m)  a -> Index rank -> Index rank
+        -> m (MutableInnerContigArray marr (PrimState m)  a )
+
+
+
+
diff --git a/src/Numerical/Array/Pure.hs b/src/Numerical/Array/Pure.hs
new file mode 100644
--- /dev/null
+++ b/src/Numerical/Array/Pure.hs
@@ -0,0 +1,170 @@
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE ScopedTypeVariables#-}
+{-# LANGUAGE DeriveDataTypeable #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE ConstraintKinds #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE FunctionalDependencies #-}
+
+module Numerical.Array.Pure      where
+
+
+--import Numerical.Array.Address
+import qualified Numerical.Array.Layout as L
+
+import Numerical.Array.Locality
+import Numerical.Array.Shape
+import Numerical.Array.Range
+import Numerical.Array.Storage as S
+import Numerical.World
+
+import qualified Data.Vector.Generic as VG
+
+
+{-
+a general question that you might ask is "what primops need have a monad constraint"
+
+ie rather than having type a -> b, why are they type Monad m => a -> m b ?
+
+the answer boils down to the following: most array types have
+a *PURE* header data structure that can't be mutated,
+that contains the Shape, extent, some handle/pointer to the associated underlying
+buffer/datastructure.  Any (even nominally pure) access to  that potentially
+mutable buffer should be mediated by a monad.
+
+I further assume that the *structure* and *extent* of this underlying buffer cannot change.
+
+That is, A valid address will always stay valid, even if after some mutation it may
+    correspond to a *different* index than it did before.
+-}
+
+{-
+Fix ME, these names are lame
+
+ImmArray == immutable array
+
+-}
+data family ImmArray world rep lay (view::Locality) (rank :: Nat )   el
+
+data instance  ImmArray Native rep lay locality rank  el =
+  ImMutableNativeArray {
+          nativeBufferPure  :: ! (S.BufferPure rep  el  )
+          ,nativeFormatPure :: ! (L.Format lay locality rank rep)
+    }
+
+
+class  PureArray arr   (rank:: Nat)   a |  arr -> rank   where
+    type PureArrayAddress (arr :: *  -> * ) ::  *
+
+    -- | gives the shape, a 'rank' length list of the dimensions
+    basicShape :: arr   a -> Index rank
+
+    --basicUnsafeRead  :: PrimMonad m => marr  (PrimState m)   a -> Shape rank Int -> m (Maybe a)
+
+    --  | basicMutableSparseIndexToAddres checks if a index is present or not
+    -- helpful primitive for authoring codes for (un)structured sparse array format
+    -- FIXME : THIS IS A TERRIBLE NAME
+    basicSparseIndexToAddress :: ( address ~PureArrayAddress  arr) => arr a -> Index rank  -> (Maybe address)
+
+    -- |
+    basicAddressToIndex :: (address ~PureArrayAddress  arr) => arr a -> address ->  (Index rank  )
+
+    -- |  return the Range of valid logical addresses
+    basicAddressRange :: (address ~PureArrayAddress  arr)=>  arr a -> Maybe (Range address)
+
+
+
+    -- | gives the next valid logical address
+    -- undefined on invalid addresses and the greatest valid address.
+    -- Note that for invalid addresses in between minAddress and maxAddress,
+    -- will return the next valid address
+    basicNextAddress :: (address ~PureArrayAddress  arr)=>  arr a -> address -> Maybe address
+
+    -- I think the case could be made for a basicPreviousAddress opeeration
+
+    -- | gives the next valid array index
+    -- undefined on invalid indices and the greatest valid index
+    basicNextIndex :: (address ~PureArrayAddress  arr)=>
+      arr a ->  Index rank -> Maybe address  -> Maybe ( Index rank, address)
+
+
+    -- | for a given valid address, @'basicAddressRegion' addr @ will return an AddressInterval
+    -- that contains @addr@. This will be a singleton when the "maximal uniform stride interval"
+    -- containing @addr@ has strictly less than 3 elements. Otherwise will return an Address range
+    -- covering the maximal interval that will have cardinality at least 3.
+
+
+    --basicAddressRegion :: (address ~PureArrayAddress  arr)=>  arr   a -> address ->  UniformAddressInterval address
+
+    ---- | Yield the element at the given position. This method should not be
+    ---- called directly, use 'unsafeRead' instead.
+    basicUnsafeAddressRead  :: (Monad m , address ~PureArrayAddress  arr)=>  arr   a -> address-> m  a
+
+
+
+    -- | Yield the element at the given position. This method should not be
+    -- called directly, use 'unsafeSparseRead' instead.
+    basicUnsafeSparseRead  :: Monad m => arr   a -> Index rank  -> m (Maybe a)
+
+-- the catch all layout instance
+
+instance (Buffer rep el , L.Layout (L.Format  lay locality  rank rep) rank)
+  =>PureArray (ImmArray Native rep lay locality rank )   rank el   where
+      type PureArrayAddress (ImmArray Native rep lay locality rank )
+        =L.LayoutAddress (L.Format lay locality rank rep)
+
+      {-# INLINE basicShape #-}
+      basicShape = L.basicLogicalShape . nativeFormatPure
+
+      {-# INLINE basicSparseIndexToAddress #-}
+      basicSparseIndexToAddress= L.basicToAddress . nativeFormatPure
+
+      {-#  INLINE basicAddressToIndex #-}
+      basicAddressToIndex = L.basicToIndex . nativeFormatPure
+
+      {-# INLINE basicAddressRange #-}
+      basicAddressRange = L.basicAddressRange . nativeFormatPure
+
+      {-# INLINE basicNextAddress #-}
+      basicNextAddress= L.basicNextAddress . nativeFormatPure
+
+      {-# INLINE basicNextIndex #-}
+      basicNextIndex = L.basicNextIndex . nativeFormatPure
+
+      {-# INLINE basicUnsafeSparseRead #-}
+      basicUnsafeSparseRead =
+          \ arr ix -> case basicSparseIndexToAddress arr ix of
+                          Nothing -> return Nothing
+                          (Just addr) ->  basicUnsafeAddressRead arr addr >>= ( return . Just)
+
+      {-# INLINE basicUnsafeAddressRead #-}
+      basicUnsafeAddressRead =
+          \ arr  addr ->
+            VG.basicUnsafeIndexM (nativeBufferPure arr)
+                (L.basicAddressAsInt (nativeFormatPure arr) $ addr)
+
+class PureArray arr rank a => PureDenseArray arr rank a where
+
+    -- | 'basicIndexInBounds' is an O(1) bounds check.
+    basicIndexInBounds :: arr a -> Index rank -> Bool
+
+    -- |
+    basicUnsafeAddressDenseRead  :: (address ~PureArrayAddress  arr,Monad m) => arr  a -> address-> m a
+
+
+
+    -- | Yield the element at the given position. This method should not be
+    -- called directly, use 'unsafeRead' instead.
+    basicUnsafeDenseReadM  :: Monad m =>  arr     a -> Index rank  -> m a
+
+
+
+
+
+
diff --git a/src/Numerical/Array/Range.hs b/src/Numerical/Array/Range.hs
new file mode 100644
--- /dev/null
+++ b/src/Numerical/Array/Range.hs
@@ -0,0 +1,74 @@
+
+{-# LANGUAGE DeriveDataTypeable #-}
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE FunctionalDependencies #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE DeriveFunctor, DeriveGeneric, DeriveFunctor #-}
+{-# LANGUAGE DeriveFoldable,DeriveTraversable #-}
+module Numerical.Array.Range (
+    Range(..)
+    ,AffineRange(..)
+    ,HasRange(..)
+    ,affineRangeStride) where
+
+import Data.Data
+import GHC.Generics
+#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ < 709
+import Data.Foldable
+import Data.Traversable
+#endif
+
+
+{-
+not quite the right module for this notion of range, but lets
+fix that later
+-}
+-- | whenever you are  tempted to do a (lo,hi) tuple, use this instead
+--  This should perhaps be made lazy, but strict for now.
+data Range a =Range {_RangeMin :: !a
+                      ,_RangeMax :: !a}
+        deriving (Eq,Show,Data,Typeable,Generic ,Foldable,Traversable,Functor)
+
+class HasRange r a | r -> a where
+  rangeMin:: Functor f => (a -> f a )-> r -> f r
+  rangeMax:: Functor f => (a -> f a )-> r -> f r
+
+instance HasRange (Range a) a where
+  rangeMax = _rangeMax
+  {-# INLINE rangeMax#-}
+
+  rangeMin = _rangeMin
+  {-# INLINE rangeMin #-}
+
+instance HasRange (AffineRange a) a where
+  rangeMin = _affineRangeMin
+  {-# INLINE rangeMin #-}
+
+  rangeMax = _affineRangeMax
+  {-# INLINE rangeMax #-}
+
+_rangeMin :: Functor f => (a -> f a)-> Range a -> f (Range a)
+_rangeMin = \ fun rec  -> fmap (\mup -> rec{_RangeMin= mup}) $ fun (_RangeMin rec )
+{-# INLINE _rangeMin#-}
+
+_rangeMax :: Functor f => (a -> f a) -> Range a -> f (Range a)
+_rangeMax =  \ fun rec -> fmap (\mup -> rec{_RangeMax= mup}) $ fun (_RangeMax rec )
+{-# INLINE _rangeMax #-}
+
+-- | this is uniform address interval by any other name
+data AffineRange a = AffineRange{_AffineRangeMin :: !a
+                                ,_AffineRangeStride :: ! Int
+                                ,_AffineRangeMax :: !a}
+        deriving (Eq,Show,Data,Generic,Typeable,Functor,Foldable,Traversable )
+
+_affineRangeMin :: Functor f => (a-> f a) -> AffineRange a -> f (AffineRange a)
+_affineRangeMin= \ fun rec -> fmap (\mup -> rec{_AffineRangeMin=mup}) $ fun (_AffineRangeMin rec)
+{-# INLINE _affineRangeMin#-}
+
+_affineRangeMax :: Functor f => (a -> f a) -> AffineRange a -> f (AffineRange a)
+_affineRangeMax= \ fun rec -> fmap (\mup -> rec{_AffineRangeMax=mup}) $ fun (_AffineRangeMax rec)
+{-# INLINE _affineRangeMax #-}
+
+affineRangeStride :: Functor f => (Int -> f Int) -> AffineRange a -> f (AffineRange a)
+affineRangeStride = \fun rec  -> fmap (\mup -> rec{_AffineRangeStride=mup}) $ fun (_AffineRangeStride rec)
+{-# INLINE affineRangeStride #-}
diff --git a/src/Numerical/Array/Shape.hs b/src/Numerical/Array/Shape.hs
new file mode 100644
--- /dev/null
+++ b/src/Numerical/Array/Shape.hs
@@ -0,0 +1,671 @@
+{-# LANGUAGE DataKinds, GADTs, TypeFamilies #-}
+{-# LANGUAGE ScopedTypeVariables  #-}
+{-# LANGUAGE ExplicitForAll  #-}
+{-# LANGUAGE DeriveDataTypeable #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE UnboxedTuples #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FunctionalDependencies #-}
+{-# LANGUAGE UndecidableInstances #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE StandaloneDeriving #-}
+{-# LANGUAGE DeriveFunctor #-}
+{-# LANGUAGE NoImplicitPrelude #-}
+
+
+module Numerical.Array.Shape(
+  -- * Shape
+    Shape(..)
+    -- * Shape Utilities
+    ,foldl
+    ,foldr
+    ,foldl'
+    ,foldl1
+    ,foldr1
+    ,map
+    ,map2
+    ,reverseShape
+    ,Nat(..)
+    ,shapeSize
+    ,SNat(..)
+    ,weaklyDominates
+    ,strictlyDominates
+    ,shapeToList
+    ,Index
+    ,backwards
+    -- * Unboxed Vector Morphism
+    ,UnBoxedShapeMorphism(..)
+    ,unShapeVector
+    ,reShapeVector
+    ,T.traverse
+    --,T.Traversable(..)
+
+    )
+    where
+
+--import Data.Data
+import Data.Typeable
+import Data.Data
+
+import qualified Data.Functor as Fun
+import qualified  Data.Foldable as F
+import qualified Control.Applicative as A
+import Control.Monad (liftM)
+import Control.Monad.ST (runST)
+import qualified Data.Traversable as T
+
+--import Control.NumericalMonad.State.Strict
+import Control.NumericalApplicative.Backwards
+
+
+import Numerical.Nat
+import qualified Data.Monoid as Monoid
+
+import Prelude hiding  (map,foldl,foldr,init,scanl,scanr,scanl1,scanr1,foldl1,foldr1)
+
+import qualified Foreign.Storable  as Store
+import qualified Foreign.Ptr as Ptr
+
+import qualified Data.Vector.Unboxed as UV
+import  qualified  Data.Vector.Unboxed.Mutable as UVM
+import qualified Data.Vector.Generic as GV
+import qualified Data.Vector.Generic.Mutable as GMV
+
+{-
+Shape may get renamed to Index in the near future!
+
+PSA: do not take the INLINE pragmas as a style suggestion.
+The only reason for the INLINEs, SPECIALIZE and the
+nonrecursive type class definitions of operations
+in this module are because shape will be used in the inner loops of
+array indexing heavy computations,
+
+-}
+
+
+ {-
+the concern basically boils down to "will it specialize / inline well"
+
+ -}
+
+
+{-
+should explore using the Reverse and Backwards transformers in the
+Transformers package, but not right now
+
+note also the *Tup operations could be done with a more general State monad
+for the tupled accumulation parameter. If theres no perf regression, should
+move to using that instead.
+
+-}
+
+
+infixr 3 :*
+
+type Index rank = Shape rank Int
+
+data Shape (rank :: Nat) a where
+    Nil  :: Shape 'Z a
+    (:*) ::  !(a) -> !(Shape r a ) -> Shape  ('S r) a
+
+
+deriving instance Typeable Shape
+
+
+nilShapeConstrRep :: Constr
+nilShapeConstrRep    = mkConstr shapeDataTypeRep "Nil" [] Prefix
+consShapeConstrRep :: Constr
+consShapeConstrRep   = mkConstr shapeDataTypeRep ":*" [] Infix
+
+shapeDataTypeRep :: DataType
+shapeDataTypeRep = mkDataType "Numerical.Array.Shape.Shape" [nilShapeConstrRep,consShapeConstrRep]
+
+
+--deriving instance (Data a, Typeable n ) => Data (Shape n a)
+--  gfoldl f z xs = gfoldl f z (shapeToList xs)
+
+--  --gfoldl _ z Nil = z Nil
+--  --gfoldl f z (x :* xs) = z (:*)  `f`  x `f` xs
+
+-- I would like to have (Data (Shape n a)) but that seems tricky
+
+instance (Data a,Typeable 'Z) =>  Data (Shape 'Z a) where
+    gfoldl _ z Nil = z Nil
+    gunfold _ z _  = z Nil -- not sure if _ z _ is the right one, but typechecks
+    dataTypeOf _  = shapeDataTypeRep
+    toConstr _ = nilShapeConstrRep
+
+instance (Data a, Data (Shape n a), Typeable ('S n))=> Data (Shape ('S n) a ) where
+    gfoldl k z (a :* b) = (z (:*) `k` a) `k` b
+    gunfold k z _ = k (k (z (:*)))
+    dataTypeOf _ = shapeDataTypeRep
+    toConstr _   = consShapeConstrRep
+
+
+-- figure this out!
+--look at  http://hackage.haskell.org/package/HList-0.3.4.1/docs/src/Data-HList-Data.html
+--and https://hackage.haskell.org/package/base-4.3.1.0/docs/Data-Data.html#nilConstr
+-- for examples
+--instance Data a => Data (Shape Z a) where
+  --gfoldl
+
+--deriving instance Data (Shape Z a)
+--deriving instance (Data a,Data (Shape n a))=> (Data (Shape (S n) a))
+
+{-
+too much work to do data instance with pre 7.8 typeable
+-}
+
+
+
+instance  Eq (Shape 'Z a) where
+    (==) _ _ = True
+    {-#INLINE (==)#-}
+instance (Eq a,Eq (Shape s a))=> Eq (Shape ('S s) a )  where
+    (==)  (a:* as) (b:* bs) =  (a == b) && (as == bs )
+    {-#INLINE (==)#-}
+instance  Show (Shape 'Z a) where
+    show _ = "Nil"
+
+instance (Show a, Show (Shape s a))=> Show (Shape ('S s) a) where
+    show (a:* as) = show a  ++ " :* " ++ show as
+
+-- at some point also try data model that
+-- has layout be dynamically reified, but for now
+-- keep it phantom typed for sanity / forcing static dispatch.
+-- NB: may need to make it more general at some future point
+--data Strided r a lay = Strided {   getStrides :: Shape r a   }
+
+-- may want to typeclassify this?
+
+
+shapeToList :: Shape n a -> [a]
+shapeToList Nil = []
+shapeToList (a:* as) = a : (shapeToList as )
+
+
+{-
+the Traversable instance for shape needs both Z and S Z base
+cases to interact nicely with the instances defined for
+foldable
+-}
+
+instance T.Traversable (Shape 'Z) where
+  traverse = \ _ Nil -> A.pure Nil
+  {-# INLINE traverse #-}
+  sequenceA = T.traverse id
+  mapM f = A.unwrapMonad . T.traverse (A.WrapMonad . f)
+  sequence = T.mapM id
+  {-#INLINE sequenceA #-}
+  {-#INLINE mapM #-}
+  {-#INLINE sequence #-}
+
+
+instance  T.Traversable (Shape ('S 'Z)) where
+  traverse = \ f (a:* as) ->  (:*) A.<$> f a A.<*> T.traverse f as
+  {-# INLINE traverse #-}
+  sequenceA = T.traverse id
+  mapM f = A.unwrapMonad . T.traverse (A.WrapMonad . f)
+  sequence = T.mapM id
+  {-#INLINE sequenceA #-}
+  {-#INLINE mapM #-}
+  {-#INLINE sequence #-}
+
+instance T.Traversable (Shape ('S n)) => T.Traversable (Shape ('S ('S n))) where
+  traverse = \ f (a:* as) ->  (:*) A.<$> f a A.<*> T.traverse f as
+  {-#INLINE traverse #-}
+  sequenceA = T.traverse id
+  mapM f = A.unwrapMonad . T.traverse (A.WrapMonad . f)
+  sequence = T.mapM id
+  {-#INLINE sequenceA #-}
+  {-#INLINE mapM #-}
+  {-#INLINE sequence #-}
+
+backwards :: (T.Traversable t, A.Applicative f) =>
+      ((a -> Backwards f b) -> t a -> Backwards f (t b))
+        -> ((a -> f b) -> t a -> f (t b))
+backwards= \ traver f container ->
+    forwards $ traver  (\x -> Backwards $ f x) container
+{-#INLINE backwards #-}
+
+
+
+
+--instance Fun.Functor (Shape r) where
+--    fmap = mapShape
+--    {-#INLINE fmap #-}
+
+instance Fun.Functor (Shape 'Z) where
+    fmap  = \ _ Nil -> Nil
+    {-# INLINE  fmap #-}
+
+instance  (Fun.Functor (Shape r)) => Fun.Functor (Shape ('S r)) where
+    fmap  = \ f (a :* rest) -> f a :* ( Fun.fmap f rest )
+    {-# INLINE  fmap  #-}
+
+instance  A.Applicative (Shape 'Z) where
+    pure = \ _ -> Nil
+    {-# INLINE  pure  #-}
+    (<*>) = \ _  _ -> Nil
+    {-# INLINE  (<*>) #-}
+
+instance  A.Applicative (Shape r)=> A.Applicative (Shape ('S r)) where
+    pure = \ a -> a :* (A.pure a)
+    {-# INLINE pure #-}
+    (<*>) = \ (f:* fs) (a :* as) ->  f a :* ((A.<*>)) fs as
+    {-# INLINE  (<*>) #-}
+
+{-
+only doing Foldable for ranks >= 1 does mean that
+we dont get the cute "rank zero arrays are references"
+property. But want foldr1 and foldl1 to always succeed
+
+lets try having rank 0 anyways, i'll be happier if i can support it
+
+-}
+
+instance    F.Foldable (Shape  'Z) where
+    foldl' = \ _  !init _->  init
+    foldr'  = \ _ !init _ ->  init
+    foldl  = \ _ init _->  init
+    foldr  = \ _ init _->   init
+    foldMap = \ _f _col -> mempty
+    {-# INLINE foldMap  #-}
+    {-#  INLINE foldl #-}
+    {-#  INLINE foldr  #-}
+    {-# INLINE foldl' #-}
+    {-#  INLINE foldr'  #-}
+    foldr1 = \ _ _ -> error "you can't call foldr1 on a rank Z(ero) Shape"
+    foldl1 =  \_ _  ->  error "you can't call foldl1 on a rank Z(ero) Shape"
+
+
+instance    F.Foldable (Shape  ('S 'Z)) where
+    foldl' = \ f !init (a:*Nil)->  f init a
+    foldr'  = \ f !init (a:*Nil)->  f a init
+    foldl  = \ f init (a:*Nil)->  f init a
+    foldr  = \ f init (a:*Nil)->  f a init
+    foldMap = \ f (a :* Nil ) -> f a
+    {-# INLINE foldMap  #-}
+    {-#  INLINE foldl #-}
+    {-#  INLINE foldr  #-}
+    {-# INLINE foldl' #-}
+    {-#  INLINE foldr'  #-}
+    foldr1 = \ _ (a:* Nil) -> a
+    foldl1 =  \ _ (a:* Nil) -> a
+    {-#  INLINE foldl1 #-}
+    {-#  INLINE foldr1 #-}
+instance ( F.Foldable (Shape ('S r)) )=> F.Foldable (Shape ('S ('S r))) where
+    foldl' = \ f  init (a:* as) -> F.foldl' f (f init a) as
+    foldr' = \f !init (a :* as ) -> f a $!  F.foldr' f init as
+    foldl  = \ f  init (a:* as) -> F.foldl' f (f init a) as
+    foldr  = \ f  init (a:* as) ->   f a $!  F.foldr f init as
+    foldl1 = \ f (a:* as) -> F.foldl' f a as
+    foldr1 = \ f (a :* as) -> F.foldr' f a as
+    foldMap = \ f (a :* as ) -> f a Monoid.<> F.foldMap f as
+    {-# INLINE foldMap  #-}
+    {-# INLINE foldl #-}
+    {-# INLINE foldr  #-}
+    {-# INLINE foldl' #-}
+    {-# INLINE foldr'  #-}
+    {-# INLINE foldl1 #-}
+    {-# INLINE foldr1 #-}
+
+instance (Semigroup a, A.Applicative (Shape n))=> (Semigroup (Shape n a)) where
+  (<>) = \ a b -> A.pure (<>) A.<*> a A.<*> b
+
+instance (Monoid.Monoid a, A.Applicative (Shape n))=> Monoid.Monoid (Shape n a) where
+  mempty = A.pure Monoid.mempty
+  mappend = \ a b -> A.pure Monoid.mappend A.<*> a A.<*> b
+
+
+
+
+{- when you lift a toral order onto vectors, you get
+interesting partial order -}
+
+-- | 'weaklyDominates' is the '<=' operator lifted onto  a sized vector to
+-- induce a partial order relation
+weaklyDominates :: (Ord a, A.Applicative  (Shape n), F.Foldable (Shape n) )=>
+                        Shape n a -> Shape n a -> Bool
+weaklyDominates = \major minor -> foldl (&&) True $! map2 (>=)  major minor
+{-# INLINE weaklyDominates #-}
+
+-- | 'strictlyDominates' is the '<' operator lifted onto  a sized vector to
+-- induce a partial order relation
+strictlyDominates :: (Ord a, A.Applicative  (Shape n), F.Foldable (Shape n) )=>
+                        Shape n a -> Shape n a -> Bool
+
+strictlyDominates  = \major minor -> foldl (&&) True $! map2 (>)  major minor
+{-# INLINE strictlyDominates #-}
+
+{-# INLINE reverseShape #-}
+reverseShape :: Shape n a -> Shape n a
+reverseShape Nil = Nil
+reverseShape r@(_ :* Nil)= r
+reverseShape (a:* b :* Nil) = b:* a :* Nil
+reverseShape (a:* b :* c:* Nil )=  c :* b :* a :* Nil
+reverseShape (a:* b :* c :* d :* Nil)= d :* c :* b :* a :* Nil
+reverseShape list = go SZero Nil list
+  where
+    go :: SNat n1 -> Shape n1  a-> Shape n2 a -> Shape (n1 + n2) a
+    go snat acc Nil = gcastWith (plus_id_r snat) acc
+    go snat acc (h :* (t :: Shape n3 a)) =
+      gcastWith (plus_succ_r snat (Proxy :: Proxy n3))
+              (go (SSucc snat) (h :* acc) t)
+
+
+{-
+TODO: abstract out all the different unrolled cases i have
+
+
+-}
+
+
+
+
+{-# INLINE map2 #-}
+map2 :: forall a b c r . (A.Applicative (Shape r))=>   (a->b ->c) -> (Shape r a) -> (Shape r b) -> (Shape r c )
+map2  = \ f shpa shpb -> f A.<$> shpa  A.<*> shpb
+
+
+{-# INLINE map #-}
+map:: forall a b r . (A.Applicative (Shape r))=> (a->b) -> (Shape r a )->( Shape r b)
+map  =  \ f shp -> f A.<$> shp
+
+
+
+{-# INLINE  foldr #-}
+foldr :: forall a b r . (F.Foldable (Shape r))=>  (a->b-> b) -> b -> Shape r a -> b
+foldr  = \ f init shp -> F.foldr  f init shp
+
+
+
+
+{-# INLINE  foldl #-}
+foldl :: forall a b r. (F.Foldable (Shape r))=> (b-> a -> b) -> b -> Shape r a -> b
+foldl  = \ f init shp -> F.foldl f init shp
+
+
+{-# INLINE foldl' #-}
+foldl' :: forall a b r . (F.Foldable (Shape r))=> (b-> a -> b) -> b -> Shape r a -> b
+foldl' = \ f init shp -> F.foldl' f init shp
+
+{-# INLINE  foldr1 #-}
+foldr1 :: forall b r . (F.Foldable (Shape ('S r)))=>  (b->b-> b)  -> Shape ('S r) b -> b
+foldr1  = \ f  shp -> F.foldr1  f  shp
+
+
+
+
+{-# INLINE  foldl1 #-}
+foldl1 :: forall  b r. (F.Foldable (Shape ('S r)))=> (b-> b -> b)  -> Shape ('S r) b -> b
+foldl1  = \ f  shp -> F.foldl1 f  shp
+
+
+
+
+
+instance Store.Storable a =>Store.Storable (Shape ('S 'Z) a) where
+    {-#INLINE sizeOf#-}
+    sizeOf = \ _ ->  (Store.sizeOf (undefined :: a))
+    -- might want to boost the alignment, but ignore for now
+    {-# INLINE alignment #-}
+    alignment = \ _ -> Store.alignment (undefined :: a )
+    {-# INLINE peek #-}
+    peek = \ptr -> do  res <- Store.peek (Ptr.castPtr ptr) ; return (res :* Nil)
+    {-# INLINE poke #-}
+    poke = \ptr (a:*_) -> Store.poke (Ptr.castPtr ptr) a
+    {-# INLINE pokeElemOff #-}
+    {-# INLINE peekElemOff #-}
+    peekElemOff = \ ptr off -> Store.peekByteOff ptr (off * Store.sizeOf (undefined ::  a ))
+    pokeElemOff ptr off val = Store.pokeByteOff ptr (off * Store.sizeOf val) val
+
+    peekByteOff ptr off = Store.peek (ptr `Ptr.plusPtr` off)
+    pokeByteOff ptr off = Store.poke (ptr `Ptr.plusPtr` off)
+    {-# INLINE peekByteOff #-}
+    {-# INLINE pokeByteOff #-}
+
+
+instance (Store.Storable a,Store.Storable (Shape ('S n) a)) =>Store.Storable (Shape ('S ('S n)) a) where
+    {-#INLINE sizeOf#-}
+    sizeOf = \ _ ->  Store.sizeOf (undefined :: a)  + Store.sizeOf (undefined :: (Shape ('S n) a ))
+    -- might want to boost the alignment, but ignore for now
+    {-# INLINE alignment #-}
+    alignment = \ _ -> Store.alignment (undefined :: a )
+    {-# INLINE peek #-}
+    peek = \ptr -> do
+                a <- Store.peek (Ptr.castPtr ptr) ;
+                as <- Store.peek (ptr `Ptr.plusPtr` Store.sizeOf (undefined :: a ))
+                return (a:* as)
+    {-# INLINE poke #-}
+    poke = \ptr (a:*as ) -> do
+                        Store.poke (Ptr.castPtr ptr) a
+                        Store.poke (ptr `Ptr.plusPtr` Store.sizeOf (undefined :: a )) as
+    {-# INLINE pokeElemOff #-}
+    {-# INLINE peekElemOff #-}
+    peekElemOff = \ ptr off -> Store.peekByteOff ptr (off * Store.sizeOf (undefined :: (Shape ('S ('S n)) a) ))
+    pokeElemOff ptr off val = Store.pokeByteOff ptr (off * Store.sizeOf val) val
+
+    peekByteOff ptr off = Store.peek (ptr `Ptr.plusPtr` off)
+    pokeByteOff ptr off = Store.poke (ptr `Ptr.plusPtr` off)
+    {-# INLINE peekByteOff #-}
+    {-# INLINE pokeByteOff #-}
+
+-- this instance is a bit weird and should never be used
+-- but probably legal
+instance Store.Storable a =>Store.Storable (Shape 'Z a) where
+    {-#INLINE sizeOf#-}
+    sizeOf = \ _ ->  Store.sizeOf (undefined :: a )
+    -- might want to boost the alignment, but ignore for now
+    {-# INLINE alignment #-}
+    alignment = \ _ -> Store.alignment (undefined :: a )
+    {-# INLINE peek #-}
+    peek = \ _  -> return Nil
+    {-# INLINE poke #-}
+    poke = \ _  _-> return ()
+    {-# INLINE pokeElemOff #-}
+    {-# INLINE peekElemOff #-}
+    peekElemOff = \ _ _  -> return Nil
+    pokeElemOff = \ _ _ _  -> return ()
+
+    peekByteOff  = \ _ _ -> return Nil
+    pokeByteOff  = \ _ _ _ -> return ()
+    {-# INLINE peekByteOff #-}
+    {-# INLINE pokeByteOff #-}
+
+{-# INLINE shapeSize #-}
+shapeSize :: F.Foldable (Shape n)=>Shape n a -> Int
+shapeSize  = \ as -> ( F.foldl (\ct _ -> ct +1) 0 as )
+
+unShapeVector ::(UnBoxedShapeMorphism n a, T.Traversable (Shape n), UV.Unbox a) => UV.Vector  (Shape n a) -> (Int, Shape n (UV.Vector  a))
+unShapeVector vs = runST  $
+            do  (l,mvs) <- fmap unShapeMVector $  UV.unsafeThaw vs
+                shpvs <- T.traverse UV.unsafeFreeze mvs
+                return (l,shpvs)
+
+
+reShapeVector::(UnBoxedShapeMorphism n a, T.Traversable (Shape n), UV.Unbox a)=>
+      (Int, Shape n (UV.Vector   a)) -> UV.Vector (Shape n a)
+reShapeVector (l,vs) = runST $
+          do  mShapeV <-  T.traverse UV.unsafeThaw  vs
+              mvShp <- return $ reShapeMVector (l,mShapeV)
+              UV.unsafeFreeze mvShp
+
+
+{- THis is a convenience type class so i dont have to export the constructors -}
+class (UV.Unbox (Shape n a)) => UnBoxedShapeMorphism n a  where
+   --unShapeVector :: UV.Vector (Shape n a) -> (Int, Shape n (UV.Vector a))
+   --reShapeVector :: (Int, Shape n (UV.Vector a)) -> UV.Vector (Shape n a)
+
+   unShapeMVector :: UVM.MVector s (Shape n a) -> (Int, Shape n (UV.MVector s a))
+   reShapeMVector :: (Int, Shape n (UVM.MVector s  a)) -> UVM.MVector s (Shape n a)
+
+instance (UV.Unbox a)=>  UnBoxedShapeMorphism 'Z a where
+  --unShapeVector (V_ShapeZ l)= (l,Nil)
+  unShapeMVector (MV_ShapeZ l) = (l,Nil )
+
+  --reShapeVector (l,Nil)  = (V_ShapeZ l)
+  reShapeMVector (l,Nil ) = (MV_ShapeZ l)
+
+
+instance (UV.Unbox a)=>  UnBoxedShapeMorphism ('S 'Z) a  where
+  --unShapeVector (V_ShapeSZ v)= (UV.length v, v :* Nil)
+
+  unShapeMVector (MV_ShapeSZ v) = (UVM.length v,v:* Nil )
+
+  --reShapeVector (l,v :* Nil)  = (V_ShapeSZ v)
+  reShapeMVector (_,v :* _ ) = (MV_ShapeSZ v)
+
+--UV.V_2
+--UVM.MV_2
+instance ((UV.Unbox a),UnBoxedShapeMorphism ('S n) a )=> UnBoxedShapeMorphism ('S ('S n)) a where
+  --unShapeVector (V_ShapeSSN (UV.V_2 l vhead vtail))= (l, vhead :* snd (unShapeVector vtail)  )
+  unShapeMVector (MV_ShapeSSN (UVM.MV_2 l vhead vtail)) = (l,vhead:*  snd (unShapeMVector vtail ))
+
+  --reShapeVector (l,vh :* vt)  = (V_ShapeSSN (UV.V_2 l vh (reShapeVector (l,vt) )))
+  reShapeMVector (l,vh :* vt ) = (MV_ShapeSSN (UVM.MV_2 l vh (reShapeMVector (l,vt) )))
+
+newtype instance UV.MVector s (Shape 'Z a)  = MV_ShapeZ  Int
+newtype instance UV.Vector    (Shape 'Z a) = V_ShapeZ  Int
+
+newtype instance UV.MVector s (Shape ('S 'Z) a)  = MV_ShapeSZ (UV.MVector s a)
+newtype instance UV.Vector    (Shape ('S 'Z) a) = V_ShapeSZ  (UV.Vector    a)
+
+newtype instance UV.MVector s (Shape ('S ('S n)) a)  = MV_ShapeSSN (UV.MVector s (a, Shape ('S n) a) )
+newtype instance UV.Vector    (Shape ('S ('S n)) a) = V_ShapeSSN  (UV.Vector   (a, Shape ('S n) a) )
+
+
+instance UV.Unbox a => UV.Unbox (Shape 'Z a)
+instance UV.Unbox a =>  UV.Unbox (Shape ('S 'Z) a)
+instance (UV.Unbox a,UV.Unbox (Shape ('S n) a) )=> UV.Unbox (Shape ('S ('S n)) a)
+
+
+
+instance UV.Unbox a => GMV.MVector UV.MVector  (Shape 'Z a) where
+  {-# INLINE basicLength #-}
+  {-# INLINE basicUnsafeSlice #-}
+  {-# INLINE basicOverlaps #-}
+  {-# INLINE basicUnsafeNew #-}
+  {-# INLINE basicUnsafeRead #-}
+  {-# INLINE basicUnsafeWrite #-}
+  {-# INLINE basicClear #-}
+  {-# INLINE basicSet #-}
+  {-# INLINE basicUnsafeCopy #-}
+  {-# INLINE basicUnsafeGrow #-}
+  {-# INLINE basicInitialize #-}
+  basicInitialize = \ (MV_ShapeZ _n) -> return ()
+  basicLength  = \ (MV_ShapeZ n) ->  n
+  basicUnsafeSlice  = \ _ m (MV_ShapeZ _) -> MV_ShapeZ m
+  basicOverlaps = \ _ _  ->  False
+  basicUnsafeNew  = \ n -> return (MV_ShapeZ n)
+  basicUnsafeRead  = \ (MV_ShapeZ _) _ ->  return Nil
+  basicUnsafeWrite  = \ (MV_ShapeZ _) _ Nil ->  return ()
+  basicClear = \ _ -> return ()
+  basicSet  =  \ (MV_ShapeZ _) Nil -> return ()
+  basicUnsafeCopy  = \ (MV_ShapeZ _) (MV_ShapeZ _) ->  return ()
+  basicUnsafeGrow  = \ (MV_ShapeZ n) m -> return $ MV_ShapeZ (n+m)
+
+instance UV.Unbox a => GV.Vector UV.Vector  (Shape 'Z a) where
+  {-# INLINE basicUnsafeFreeze #-}
+  basicUnsafeFreeze  = \ (MV_ShapeZ n) ->  return $ V_ShapeZ n
+  {-# INLINE basicUnsafeThaw #-}
+  basicUnsafeThaw = \ (V_ShapeZ n)  -> return $ MV_ShapeZ n
+  {-# INLINE basicLength #-}
+  basicLength  = \(V_ShapeZ n) ->  n
+  {-# INLINE basicUnsafeSlice #-}
+  basicUnsafeSlice  = \ _ m (V_ShapeZ _) ->  V_ShapeZ m
+  {-# INLINE basicUnsafeIndexM #-}
+  basicUnsafeIndexM = \ (V_ShapeZ _) _  ->  return Nil
+  {-# INLINE basicUnsafeCopy #-}
+  basicUnsafeCopy  = \ (MV_ShapeZ _) (V_ShapeZ _) ->  return ()
+  {-# INLINE elemseq #-}
+  elemseq  =  \ _ -> seq
+
+instance (UV.Unbox a) => GMV.MVector UV.MVector (Shape ('S 'Z) a) where
+  {-# INLINE basicLength #-}
+  {-# INLINE basicUnsafeSlice #-}
+  {-# INLINE basicOverlaps #-}
+  {-# INLINE basicUnsafeNew #-}
+  {-# INLINE basicUnsafeReplicate #-}
+  {-# INLINE basicUnsafeRead #-}
+  {-# INLINE basicUnsafeWrite #-}
+  {-# INLINE basicClear #-}
+  {-# INLINE basicSet #-}
+  {-# INLINE basicUnsafeCopy #-}
+  {-# INLINE basicUnsafeGrow #-}
+  {-# INLINE basicInitialize #-}
+  basicInitialize = \ (MV_ShapeSZ v) ->  GMV.basicInitialize v
+  basicLength  = \(MV_ShapeSZ v)-> GMV.basicLength v
+  basicUnsafeSlice  = \ i n (MV_ShapeSZ v) ->  MV_ShapeSZ $ GMV.basicUnsafeSlice i n v
+  basicOverlaps = \ (MV_ShapeSZ v1) (MV_ShapeSZ v2)  ->  GMV.basicOverlaps v1 v2
+  basicUnsafeNew  = \ n ->  MV_ShapeSZ `liftM` GMV.basicUnsafeNew n
+  basicUnsafeReplicate  = \ n (a:*_) ->  MV_ShapeSZ `liftM` GMV.basicUnsafeReplicate n a
+  basicUnsafeRead  = \ (MV_ShapeSZ v) i ->  ( :* Nil ) `liftM` GMV.basicUnsafeRead v i
+  basicUnsafeWrite  = \ (MV_ShapeSZ v) i (a:* _) ->  GMV.basicUnsafeWrite v i a
+  basicClear = \ (MV_ShapeSZ v)  ->  GMV.basicClear v
+  basicSet =  \ (MV_ShapeSZ v) (a:*_)  ->  GMV.basicSet v a
+  basicUnsafeCopy  = \ (MV_ShapeSZ v1) (MV_ShapeSZ v2) ->  GMV.basicUnsafeCopy v1 v2
+  basicUnsafeMove  = \ (MV_ShapeSZ v1) (MV_ShapeSZ v2) -> GMV.basicUnsafeMove v1 v2
+  basicUnsafeGrow  = \ (MV_ShapeSZ v) n ->  MV_ShapeSZ `liftM` GMV.basicUnsafeGrow v n
+
+instance ( UV.Unbox a) => GV.Vector UV.Vector (Shape ('S 'Z) a ) where
+  {-# INLINE basicUnsafeFreeze #-}
+  {-# INLINE basicUnsafeThaw #-}
+  {-# INLINE basicLength #-}
+  {-# INLINE basicUnsafeSlice #-}
+  {-# INLINE basicUnsafeIndexM #-}
+  {-# INLINE elemseq #-}
+  basicUnsafeFreeze = \ (MV_ShapeSZ v)  ->  V_ShapeSZ `liftM` GV.basicUnsafeFreeze v
+  basicUnsafeThaw = \ (V_ShapeSZ v)  ->  MV_ShapeSZ`liftM` GV.basicUnsafeThaw v
+  basicLength  = \ (V_ShapeSZ v)-> GV.basicLength v
+  basicUnsafeSlice  = \ i n (V_ShapeSZ v) ->  V_ShapeSZ $ GV.basicUnsafeSlice i n v
+  basicUnsafeIndexM  = \ (V_ShapeSZ v) i -> ( :* Nil ) `liftM` GV.basicUnsafeIndexM v i
+  basicUnsafeCopy   = \ (MV_ShapeSZ mv) (V_ShapeSZ v) -> GV.basicUnsafeCopy mv v
+  elemseq  = \ _ (a:*_) z ->    GV.elemseq (undefined :: UV.Vector a) a z
+
+
+instance (UV.Unbox a,UV.Unbox (Shape ('S n) a)) => GMV.MVector UV.MVector (Shape ('S ('S n)) a) where
+  {-# INLINE basicLength #-}
+  {-# INLINE basicUnsafeSlice #-}
+  {-# INLINE basicOverlaps #-}
+  {-# INLINE basicUnsafeNew #-}
+  {-# INLINE basicUnsafeReplicate #-}
+  {-# INLINE basicUnsafeRead #-}
+  {-# INLINE basicUnsafeWrite #-}
+  {-# INLINE basicClear #-}
+  {-# INLINE basicSet #-}
+  {-# INLINE basicUnsafeCopy #-}
+  {-# INLINE basicUnsafeGrow #-}
+  {-# INLINE basicInitialize #-}
+  basicInitialize = \ (MV_ShapeSSN v) ->  GMV.basicInitialize v
+  basicLength  = \ (MV_ShapeSSN v) -> GMV.basicLength v
+  basicUnsafeSlice  = \ i n (MV_ShapeSSN v) -> MV_ShapeSSN $ GMV.basicUnsafeSlice i n v
+  basicOverlaps  = \ (MV_ShapeSSN v1) (MV_ShapeSSN v2) -> GMV.basicOverlaps v1 v2
+  basicUnsafeNew = \ n -> MV_ShapeSSN `liftM` GMV.basicUnsafeNew n
+  basicUnsafeReplicate  = \ n (a :* as) ->  MV_ShapeSSN `liftM` GMV.basicUnsafeReplicate n (a,as)
+  basicUnsafeRead = \ (MV_ShapeSSN v) i  ->  uncurry (:*) `liftM` GMV.basicUnsafeRead v i
+  basicUnsafeWrite = \(MV_ShapeSSN v) i (a :* as )  -> GMV.basicUnsafeWrite v i (a,as)
+  basicClear = \ (MV_ShapeSSN v)  ->  GMV.basicClear v
+  basicSet  = \ (MV_ShapeSSN v) (a :* as) ->  GMV.basicSet v (a,as)
+  basicUnsafeCopy  = \ (MV_ShapeSSN v1) (MV_ShapeSSN v2) -> GMV.basicUnsafeCopy v1 v2
+  basicUnsafeMove  = \ (MV_ShapeSSN v1) (MV_ShapeSSN v2) ->  GMV.basicUnsafeMove v1 v2
+  basicUnsafeGrow = \ (MV_ShapeSSN v) n  -> MV_ShapeSSN `liftM` GMV.basicUnsafeGrow v n
+
+
+instance (UV.Unbox a,UV.Unbox (Shape ('S n) a)) =>  GV.Vector UV.Vector (Shape ('S ('S n)) a) where
+  {-# INLINE basicUnsafeFreeze #-}
+  {-# INLINE basicUnsafeThaw #-}
+  {-# INLINE basicLength #-}
+  {-# INLINE basicUnsafeSlice #-}
+  {-# INLINE basicUnsafeIndexM #-}
+  {-# INLINE elemseq #-}
+  basicUnsafeFreeze  = \ (MV_ShapeSSN v) ->  V_ShapeSSN `liftM` GV.basicUnsafeFreeze v
+  basicUnsafeThaw = \ (V_ShapeSSN v)  ->  MV_ShapeSSN `liftM` GV.basicUnsafeThaw v
+  basicLength = \ (V_ShapeSSN v)  -> GV.basicLength v
+  basicUnsafeSlice = \ i n (V_ShapeSSN v)  -> V_ShapeSSN $ GV.basicUnsafeSlice i n v
+  basicUnsafeIndexM  = \ (V_ShapeSSN v)  i -> uncurry (:*) `liftM` GV.basicUnsafeIndexM v i
+  basicUnsafeCopy    =  \ (MV_ShapeSSN mv) (V_ShapeSSN v) -> GV.basicUnsafeCopy mv v
+  elemseq = \  _ (a :* as) z ->  GV.elemseq (undefined :: UV.Vector a) a
+                       $ GV.elemseq (undefined :: UV.Vector (Shape ('S n) a)) as z
diff --git a/src/Numerical/Array/Storage.hs b/src/Numerical/Array/Storage.hs
new file mode 100644
--- /dev/null
+++ b/src/Numerical/Array/Storage.hs
@@ -0,0 +1,267 @@
+
+{-# LANGUAGE TypeFamilies,FlexibleInstances,MultiParamTypeClasses,FlexibleContexts #-}
+{-# LANGUAGE UndecidableInstances,StandaloneDeriving,  DeriveDataTypeable #-}
+{-# LANGUAGE DeriveFunctor, DeriveFoldable, DeriveTraversable, DeriveGeneric #-}
+{-# LANGUAGE CPP #-}
+module Numerical.Array.Storage(
+  Boxed
+  ,Unboxed
+  ,Stored
+  ,BufferPure(..)
+  ,BufferMut(..)
+  ,Buffer
+  ,MBuffer
+  ,unsafeBufferThaw
+  ,unsafeBufferFreeze) where
+
+
+import Control.Monad.Primitive ( PrimMonad, PrimState )
+
+import qualified Data.Vector.Generic as VG
+import qualified Data.Vector.Generic.Mutable as VGM
+import qualified Data.Vector as BV
+import qualified Data.Vector.Storable as SV
+import qualified Data.Vector.Unboxed as UV
+
+--import qualified Data.Functor as F hiding (Functor)
+--import qualified Data.Foldable as F hiding (Foldable)
+--import qualified Data.Traversable  as T hiding (Traversable)
+#if  defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ < 709
+--import Data.Functor (Functor)
+import Data.Foldable (Foldable)
+import Data.Traversable (Traversable)
+#endif
+
+import Data.Typeable
+import Data.Data
+import GHC.Generics
+
+
+{-
+FIXME : should i require that the element type and
+mode are both instance of Typeable for Buffers?
+
+-}
+
+
+{-
+FIX MEEEEE REMINDERS
+make the allocators for   Storable Buffers  do AVX sized alignment
+-}
+
+-- | The class instance @'Buffer' mode a@ is a shorthand for saying that a given buffer representation @mode@
+-- has a 'VG.Vector' instance for both 'BufferPure'  and  'BufferMut'.
+class (VG.Vector (BufferPure mode) a, VGM.MVector (BufferMut mode) a)=> Buffer mode a
+
+instance (VG.Vector (BufferPure mode) a, VGM.MVector (BufferMut mode) a)=> Buffer mode a
+
+-- not sure if MBuffer class should exist, fixme. if/when removed, this
+class VGM.MVector (BufferMut mode) a=> MBuffer mode a
+
+-- not sure if MBuffer should exist, FIXME
+instance VGM.MVector (BufferMut mode) a=> MBuffer mode a
+
+-- | 'Boxed' is the type index for `Buffer`s that use the  boxed data structure `Data.Vector.Vector`
+-- as the underlying storage representation.
+data Boxed
+  deriving Typeable
+
+deriving instance Data Boxed
+
+-- | 'Unboxed' is the type index for 'Buffer's that use the unboxed data structure
+-- 'Data.Vector.Unboxed.Vector' as the underlying storage representation.
+data Unboxed
+  deriving Typeable
+
+deriving instance Data Unboxed
+
+-- | 'Stored' is the type index for 'Buffer's that use the 'Foreign.Storable'
+-- for values, in pinned byte array  buffers, provided by 'Data.Vector.Storable'
+data Stored
+  deriving Typeable
+
+deriving instance Data Stored
+
+type instance VG.Mutable (BufferPure sort) = BufferMut sort
+
+
+data family   BufferPure sort  elem
+
+deriving instance Typeable BufferPure
+
+newtype instance BufferPure Boxed elem = BoxedBuffer (BV.Vector elem)
+  deriving (Show,Data,Generic,Functor,Foldable,Traversable)
+
+
+
+newtype instance BufferPure Unboxed elem = UnboxedBuffer (UV.Vector elem)
+  deriving (Show,Data,Generic)
+--deriving instance Typeable a => Typeable (BufferPure Unboxed a)
+
+newtype instance BufferPure Stored elem = StorableBuffer (SV.Vector elem)
+  deriving (Show,Data,Generic)
+
+data family   BufferMut sort st elem
+deriving instance Typeable BufferMut
+
+
+newtype instance BufferMut Boxed st   elem = BoxedBufferMut (BV.MVector st elem)
+  --deriving (Show,Data,Generic)
+newtype instance BufferMut Unboxed st elem = UnboxedBufferMut (UV.MVector st elem)
+  --deriving (Show,Data,Generic)
+newtype instance BufferMut Stored st  elem = StorableBufferMut (SV.MVector st elem)
+
+-- | 'unsafeBufferFreeze'
+unsafeBufferFreeze :: (Buffer rep a,PrimMonad m) => BufferMut rep (PrimState m )  a -> m (BufferPure rep a)
+unsafeBufferFreeze =  VG.basicUnsafeFreeze
+
+unsafeBufferThaw :: (Buffer rep a,PrimMonad m) => (BufferPure rep a) -> m (BufferMut rep (PrimState m )  a)
+unsafeBufferThaw = VG.basicUnsafeThaw
+
+instance (VGM.MVector BV.MVector elem) => VGM.MVector (BufferMut Boxed)  elem where
+  basicInitialize = \(BoxedBufferMut v) -> VGM.basicInitialize v
+  basicLength = \(BoxedBufferMut v) -> VGM.basicLength v
+  basicUnsafeSlice =
+    \ ix1 ix2 (BoxedBufferMut bv) ->
+      BoxedBufferMut $ VGM.basicUnsafeSlice ix1 ix2 bv
+  basicOverlaps =
+    \ (BoxedBufferMut bv1) (BoxedBufferMut bv2) -> VGM.basicOverlaps bv1 bv2
+  basicUnsafeNew = \ size ->
+      do
+        res<- VGM.basicUnsafeNew size
+        return  (BoxedBufferMut res)
+  basicUnsafeRead= \(BoxedBufferMut bv) ix -> VGM.basicUnsafeRead bv ix
+  basicUnsafeWrite = \(BoxedBufferMut bv ) ix val -> VGM.basicUnsafeWrite bv ix val
+
+  {-Q/todo/check fixme, do these other operations need be provided in a pass through way too?
+  or will there be no difference in the derived code perf ? -}
+--  basicUnsafeClear
+--  basicUnsafeSet
+--  basicUnsafeCopy
+--  basicUnsafeMove
+--  basicUnsafeGrow
+--  basicUnsafeReplicate
+  {-# INLINE basicInitialize #-}
+  {-# INLINE basicLength #-}
+  {-# INLINE basicUnsafeSlice #-}
+  {-# INLINE basicOverlaps #-}
+  {-# INLINE basicUnsafeNew #-}
+  {-# INLINE basicUnsafeRead #-}
+  {-# INLINE basicUnsafeWrite #-}
+
+--  {-# INLINE basicUnsafeClear#-}
+--  {-# INLINE basicUnsafeSet#-}
+--  {-# INLINE basicUnsafeCopy#-}
+--  {-# INLINE basicUnsafeMove#-}
+--  {-# INLINE basicUnsafeGrow#-}
+--  {-# INLINE basicUnsafeReplicate#-}
+
+instance (SV.Storable elem) => VGM.MVector (BufferMut Stored)  elem where
+  basicInitialize = \(StorableBufferMut v) -> VGM.basicInitialize v
+  basicLength = \(StorableBufferMut v) -> VGM.basicLength v
+  basicUnsafeSlice =
+    \ ix1 ix2 (StorableBufferMut bv) ->
+      StorableBufferMut $ VGM.basicUnsafeSlice ix1 ix2 bv
+  basicOverlaps =
+    \ (StorableBufferMut bv1) (StorableBufferMut bv2) -> VGM.basicOverlaps bv1 bv2
+  basicUnsafeNew = \ size ->
+      do
+        res<- VGM.basicUnsafeNew size
+        return  (StorableBufferMut res)
+  basicUnsafeRead= \(StorableBufferMut bv) ix -> VGM.basicUnsafeRead bv ix
+  basicUnsafeWrite = \(StorableBufferMut bv ) ix val -> VGM.basicUnsafeWrite bv ix val
+  {-# INLINE basicInitialize #-}
+  {-# INLINE basicLength #-}
+  {-# INLINE basicUnsafeSlice #-}
+  {-# INLINE basicOverlaps #-}
+  {-# INLINE basicUnsafeNew #-}
+  {-# INLINE basicUnsafeRead #-}
+  {-# INLINE basicUnsafeWrite #-}
+
+
+instance (VGM.MVector UV.MVector elem) => VGM.MVector (BufferMut Unboxed)  elem where
+  {-# INLINE basicInitialize #-}
+  basicInitialize = \(UnboxedBufferMut v) -> VGM.basicInitialize v
+  basicLength = \(UnboxedBufferMut v) -> VGM.basicLength v
+  basicUnsafeSlice =
+    \ ix1 ix2 (UnboxedBufferMut bv) ->
+      UnboxedBufferMut $ VGM.basicUnsafeSlice ix1 ix2 bv
+  basicOverlaps =
+    \ (UnboxedBufferMut bv1) (UnboxedBufferMut bv2) -> VGM.basicOverlaps bv1 bv2
+  basicUnsafeNew = \ size ->
+      do
+        res<- VGM.basicUnsafeNew size
+        return  (UnboxedBufferMut res)
+  basicUnsafeRead= \(UnboxedBufferMut bv) ix -> VGM.basicUnsafeRead bv ix
+  basicUnsafeWrite = \(UnboxedBufferMut bv ) ix val -> VGM.basicUnsafeWrite bv ix val
+
+  {-# INLINE basicLength #-}
+  {-# INLINE basicUnsafeSlice #-}
+  {-# INLINE basicOverlaps #-}
+  {-# INLINE basicUnsafeNew #-}
+  {-# INLINE basicUnsafeRead #-}
+  {-# INLINE basicUnsafeWrite #-}
+
+----
+----
+instance VG.Vector BV.Vector  a  => VG.Vector (BufferPure Boxed) a   where
+
+  basicUnsafeFreeze =
+     \(BoxedBufferMut mv) ->(\ x->return $ BoxedBuffer x) =<<  VG.basicUnsafeFreeze mv
+  basicUnsafeThaw= \(BoxedBuffer v) ->(\x -> return $ BoxedBufferMut x ) =<< VG.basicUnsafeThaw v
+  basicLength = \(BoxedBuffer v) -> VG.basicLength v
+  basicUnsafeSlice =
+    \ start len (BoxedBuffer v) ->  BoxedBuffer $! VG.basicUnsafeSlice start len v
+  basicUnsafeIndexM =
+    \ (BoxedBuffer v) ix  -> VG.basicUnsafeIndexM v ix
+  elemseq = \ (BoxedBuffer v) a b -> VG.elemseq v a b
+
+
+  {-# INLINE basicUnsafeFreeze #-}
+  {-# INLINE basicUnsafeThaw #-}
+  {-# INLINE basicLength #-}
+  {-# INLINE basicUnsafeSlice #-}
+  {-# INLINE basicUnsafeIndexM #-}
+  {-# INLINE elemseq  #-}
+
+
+instance (SV.Storable a)  => VG.Vector (BufferPure Stored) a   where
+
+  basicUnsafeFreeze =
+     \(StorableBufferMut mv) -> (\x ->return $StorableBuffer x) =<<  VG.basicUnsafeFreeze mv
+  basicUnsafeThaw=
+    \(StorableBuffer v) -> (\x -> return $ StorableBufferMut x) =<< VG.basicUnsafeThaw v
+  basicLength = \(StorableBuffer v) -> VG.basicLength v
+  basicUnsafeSlice =
+    \ start len (StorableBuffer v) ->  StorableBuffer $! VG.basicUnsafeSlice start len v
+  basicUnsafeIndexM =
+    \ (StorableBuffer v) ix  -> VG.basicUnsafeIndexM v ix
+  elemseq = \ (StorableBuffer v) a b -> VG.elemseq v a b
+
+
+  {-# INLINE basicUnsafeFreeze #-}
+  {-# INLINE basicUnsafeThaw #-}
+  {-# INLINE basicLength #-}
+  {-# INLINE basicUnsafeSlice #-}
+  {-# INLINE basicUnsafeIndexM #-}
+  {-# INLINE elemseq  #-}
+
+
+instance VG.Vector UV.Vector  a  => VG.Vector (BufferPure Unboxed) a   where
+
+  basicUnsafeFreeze = \(UnboxedBufferMut mv) -> (\x -> return $ UnboxedBuffer x) =<<  VG.basicUnsafeFreeze mv
+  basicUnsafeThaw= \(UnboxedBuffer v) ->(\x -> return $  UnboxedBufferMut x) =<< VG.basicUnsafeThaw v
+  basicLength = \(UnboxedBuffer v) -> VG.basicLength v
+  basicUnsafeSlice =
+    \ start len (UnboxedBuffer v) ->  UnboxedBuffer $! VG.basicUnsafeSlice start len v
+  basicUnsafeIndexM =
+    \ (UnboxedBuffer v) ix  -> VG.basicUnsafeIndexM v ix
+  elemseq = \ (UnboxedBuffer v) a b -> VG.elemseq v a b
+
+
+  {-# INLINE basicUnsafeFreeze #-}
+  {-# INLINE basicUnsafeThaw #-}
+  {-# INLINE basicLength #-}
+  {-# INLINE basicUnsafeSlice #-}
+  {-# INLINE basicUnsafeIndexM #-}
+  {-# INLINE elemseq  #-}
diff --git a/src/Numerical/Data/Vector/HPair.hs b/src/Numerical/Data/Vector/HPair.hs
new file mode 100644
--- /dev/null
+++ b/src/Numerical/Data/Vector/HPair.hs
@@ -0,0 +1,235 @@
+{- | This  module is pretty cool because it gives you a way to talk about
+heterogeneous representations for different columns!
+
+might be replaced with an HList of Vectors approach
+-}
+
+
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeFamilyDependencies#-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE PolyKinds #-}
+{-# LANGUAGE MultiParamTypeClasses ,FlexibleInstances  , FlexibleContexts,UndecidableInstances #-}
+
+module  Numerical.Data.Vector.HPair(
+    VHProd(..)
+    ,vHPair
+    ,vUnHPair
+    ,MVHProd(..)
+    ,HProd(..)
+    ,MutableHProdTree
+    ,TransformHProdTree
+    --,mvUnPair
+    --,mvPair
+      ) where
+
+import qualified Data.Vector.Generic as V
+import qualified Data.Vector.Generic.Mutable as MV
+
+import Control.Monad.Primitive (PrimMonad)
+
+
+
+--type instance V.Mutable (VPair v) = MVPair (V.Mutable v)
+
+
+{-
+currently primmonad doesn't get its free applicative/functor powers :*(
+
+-}
+
+(<$$$>) :: PrimMonad m => (a->b) -> m a -> m b
+(<$$$>) f mv = do v <- mv ; return (f v )
+{-# INLINE (<$$$>) #-}
+
+(<***>) :: PrimMonad m => m (a->b) -> m a -> m b
+(<***>) mf mv =  do f <- mf ; v <- mv ; return (f v)
+{-# INLINE (<***>) #-}
+
+{-
+probably should just
+
+-}
+
+{-
+the names are terrible, fix them later!
+HProd , HPair, HUnit, VHPro
+-}
+
+
+data HProd a  where
+    HPair :: HProd a-> HProd a  -> HProd a
+    HUnit :: a -> HProd a
+
+data  VHProd  (prd:: HProd ( * -> * )) val where
+    VHLeaf ::  !(v a) -> VHProd   ('HUnit v) a
+    VHNode  :: !(VHProd  pra a) -> !(VHProd  prb b ) ->VHProd  ('HPair  pra prb) (a,b)
+
+data  MVHProd   (prd:: HProd (* -> * -> *) ) (st :: * ) val where
+    MVHLeaf :: !(mv  st a) -> MVHProd   ('HUnit mv) st  a
+    MVHNode  :: !(MVHProd pra st a) -> !(MVHProd   prb   st b ) -> MVHProd  ('HPair pra prb) st (a,b)
+
+
+vHPair :: (va a,vb b)->VHProd ('HPair ('HUnit va) ('HUnit vb)) (a,b)
+vHPair  = \ (va,vb) ->  VHNode (VHLeaf va) (VHLeaf vb)
+{-# INLINE vHPair #-}
+
+vUnHPair  :: VHProd  ('HPair ('HUnit va) ('HUnit vb)) (a,b) -> (va a, vb b)
+vUnHPair = \ (VHNode (VHLeaf va) (VHLeaf vb))-> (va,vb)
+{-# INLINE vUnHPair #-}
+
+type instance  V.Mutable (VHProd  prod)= MVHProd  (MutableHProdTree prod)
+
+type family MutableHProdTree (a :: HProd (* -> *))  = r | r -> a where
+  MutableHProdTree ('HUnit v ) = 'HUnit (V.Mutable v)
+  MutableHProdTree ('HPair left right) = 'HPair (MutableHProdTree left) (MutableHProdTree right )
+
+type family TransformHProdTree (f :: k-> m) (a :: HProd k) :: HProd m where
+  TransformHProdTree f ('HUnit v)= 'HUnit (f v)
+  TransformHProdTree f ('HPair left right) = 'HPair (TransformHProdTree f left) (TransformHProdTree f right)
+
+
+
+
+--mvPair :: (mv st a,mv st b)->MVPair mv st (a,b)
+--mvPair  = \ (mva, mvb) ->  TheMVPair mva mvb
+--{-# INLINE mvPair #-}
+
+--mvUnPair  :: MVPair mv st  (a,b) -> (mv st a,mv st b)
+--mvUnPair = \ (TheMVPair mva mvb)-> (mva,mvb)
+--{-# INLINE mvUnPair #-}
+
+instance  (MV.MVector (MVHProd  (MutableHProdTree ('HPair pa pb )) ) (a,b) ,
+  V.Vector (VHProd  pa) a, V.Vector (VHProd  pb) b)
+  => V.Vector (VHProd  ('HPair pa pb )) (a,b)  where
+    {-# INLINE  basicUnsafeFreeze #-}
+    {-# INLINE basicUnsafeThaw #-}
+    {-# INLINE basicLength #-}
+    {-# INLINE basicUnsafeSlice #-}
+    {-# INLINE basicUnsafeIndexM #-}
+
+
+    basicUnsafeFreeze = \(MVHNode  mva mvb) ->
+      VHNode <$$$> V.basicUnsafeFreeze mva <***> V.basicUnsafeFreeze mvb
+
+
+    basicUnsafeThaw = \(VHNode va vb) ->
+      MVHNode <$$$> V.basicUnsafeThaw va <***> V.basicUnsafeThaw vb
+
+
+    basicLength = \(VHNode va _) -> V.basicLength va
+
+
+    basicUnsafeSlice = \start len (VHNode va vb) ->
+      VHNode (V.basicUnsafeSlice start len va) (V.basicUnsafeSlice start len vb)
+
+
+    basicUnsafeIndexM = \(VHNode va vb) ix ->
+      do
+          a <- V.basicUnsafeIndexM va ix
+          b <- V.basicUnsafeIndexM vb ix
+          return (a,b)
+
+instance  (MV.MVector (MVHProd  ('HUnit (V.Mutable v))  ) a ,V.Vector v a)
+  => V.Vector (VHProd  ('HUnit v)) a  where
+
+    {-# INLINE  basicUnsafeFreeze #-}
+    {-# INLINE basicUnsafeThaw #-}
+    {-# INLINE basicLength #-}
+    {-# INLINE basicUnsafeSlice #-}
+    {-# INLINE basicUnsafeIndexM #-}
+
+    basicUnsafeFreeze = \(MVHLeaf mva) ->
+      VHLeaf <$$$> V.basicUnsafeFreeze mva
+    basicUnsafeThaw = \(VHLeaf va ) ->
+      MVHLeaf <$$$> V.basicUnsafeThaw va
+    basicLength = \(VHLeaf va ) -> V.basicLength va
+    basicUnsafeSlice = \start len (VHLeaf va ) ->
+      VHLeaf(V.basicUnsafeSlice start len va)
+    basicUnsafeIndexM = \(VHLeaf va) ix ->  V.basicUnsafeIndexM va ix
+
+
+instance (MV.MVector mv a) => MV.MVector (MVHProd  ('HUnit mv )) a where
+  basicLength = \ (MVHLeaf mva) -> MV.basicLength mva
+  {-# INLINE basicLength #-}
+
+  basicUnsafeSlice = \ start len (MVHLeaf mva  )->
+    MVHLeaf (MV.basicUnsafeSlice start len mva)
+  {-# INLINE basicUnsafeSlice #-}
+
+  basicOverlaps = \ (MVHLeaf mva ) (MVHLeaf mva2 )-> (MV.basicOverlaps mva mva2)
+  {-# INLINE basicOverlaps #-}
+
+  basicUnsafeNew =
+      \ size ->
+          MVHLeaf <$$$> MV.basicUnsafeNew size
+  {-# INLINE basicUnsafeNew #-}
+
+  basicUnsafeReplicate =
+      \ size a ->
+         MVHLeaf <$$$>
+            MV.basicUnsafeReplicate size a
+  {-# INLINE basicUnsafeReplicate #-}
+
+  basicUnsafeRead = \(MVHLeaf mva ) ix ->   MV.basicUnsafeRead mva ix
+  {-#INLINE basicUnsafeRead #-}
+
+  basicUnsafeWrite = \ (MVHLeaf mva ) ix a  ->
+    do
+      MV.basicUnsafeWrite mva ix a
+      return ()
+  {-#INLINE basicUnsafeWrite #-}
+
+  {-#INLINE basicUnsafeGrow #-}
+  basicUnsafeGrow = \ (MVHLeaf mva ) growth ->
+      MVHLeaf <$$$> MV.basicUnsafeGrow mva growth
+
+
+
+instance (MV.MVector (MVHProd pra) a,MV.MVector (MVHProd  prb) b)
+  => MV.MVector (MVHProd  ('HPair pra prb)) (a,b) where
+
+  basicLength = \ (MVHNode mva _) -> MV.basicLength mva
+  {-# INLINE basicLength #-}
+
+  basicUnsafeSlice = \ start len (MVHNode mva mvb )->
+    MVHNode (MV.basicUnsafeSlice start len mva) (MV.basicUnsafeSlice start len mvb)
+  {-# INLINE basicUnsafeSlice #-}
+
+  basicOverlaps = \ (MVHNode mva mvb) (MVHNode mva2 mvb2)-> (MV.basicOverlaps mva mva2) || (MV.basicOverlaps mvb mvb2)
+  {-# INLINE basicOverlaps #-}
+
+  basicUnsafeNew =
+      \ size ->
+          MVHNode <$$$> MV.basicUnsafeNew size <***> MV.basicUnsafeNew size
+  {-# INLINE basicUnsafeNew #-}
+
+  basicUnsafeReplicate =
+      \ size (a,b) ->
+         MVHNode <$$$>
+            MV.basicUnsafeReplicate size a <***>
+            MV.basicUnsafeReplicate size b
+  {-# INLINE basicUnsafeReplicate #-}
+
+  basicUnsafeRead = \(MVHNode mva mvb) ix ->
+    (,) <$$$>  MV.basicUnsafeRead mva ix <***> MV.basicUnsafeRead mvb ix
+
+  {-#INLINE basicUnsafeRead #-}
+
+  basicUnsafeWrite = \ (MVHNode mva mvb) ix (a,b) ->
+    do
+      MV.basicUnsafeWrite mva ix a
+      MV.basicUnsafeWrite mvb ix b
+      return ()
+  {-#INLINE basicUnsafeWrite #-}
+
+  {-#INLINE basicUnsafeGrow #-}
+  basicUnsafeGrow = \ (MVHNode mva mvb) growth ->
+      MVHNode <$$$> MV.basicUnsafeGrow mva growth <***>
+          MV.basicUnsafeGrow mvb growth
+
+
+
+
+
diff --git a/src/Numerical/Data/Vector/Pair.hs b/src/Numerical/Data/Vector/Pair.hs
new file mode 100644
--- /dev/null
+++ b/src/Numerical/Data/Vector/Pair.hs
@@ -0,0 +1,205 @@
+
+{- | This  module is pretty cool because it gives you a way to talk about
+open struct of arrays style vectors
+
+might be replaced with an HList of Vectors approach
+
+
+-}
+
+
+
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE MultiParamTypeClasses ,FlexibleInstances  , FlexibleContexts,UndecidableInstances #-}
+
+module  Numerical.Data.Vector.Pair(
+    VProd(..)
+    ,vPair
+    ,vUnPair
+    ,MVProd(..)
+    --,mvUnPair
+    ,Prod(..)
+    --,mvPair
+      ) where
+
+import qualified Data.Vector.Generic as V
+import qualified Data.Vector.Generic.Mutable as MV
+
+
+--type instance V.Mutable (VPair v) = MVPair (V.Mutable v)
+
+
+
+
+data Prod = Pair Prod Prod | Unit
+
+
+data family   VProd  (vect :: * -> * ) (prd:: Prod ) val  -- where
+data instance VProd v 'Unit a where
+    VLeaf ::  !(v a) -> VProd v   'Unit a
+
+data instance VProd v ('Pair pra prb )  (a,b) where
+    VPair  :: !(VProd v pra a) -> !(VProd v prb b ) ->VProd v ('Pair  pra prb) (a,b)
+
+data family   MVProd  (vect :: * -> * -> * )  (prd:: Prod ) (st :: * ) val  -- where
+data instance   MVProd mv 'Unit  st a where
+  MVLeaf :: !(mv  st a) -> MVProd mv  'Unit st  a
+data instance   MVProd mv ('Pair pra prb)  st (a,b) where
+    MVPair  :: !(MVProd mv pra st a) -> !(MVProd mv  prb   st b ) -> MVProd mv  ('Pair pra prb) st (a,b)
+
+
+vPair :: (v a,v b)->VProd v ('Pair 'Unit 'Unit) (a,b)
+vPair  = \ (va,vb) ->  VPair (VLeaf va) (VLeaf vb)
+{-# INLINE vPair #-}
+
+vUnPair  :: VProd v ('Pair 'Unit 'Unit) (a,b) -> (v a, v b)
+vUnPair = \ (VPair (VLeaf va) (VLeaf vb))-> (va,vb)
+{-# INLINE vUnPair #-}
+
+type instance  V.Mutable (VProd vec prod)= MVProd (V.Mutable vec) prod
+
+
+--mvPair :: (mv st a,mv st b)->MVPair mv st (a,b)
+--mvPair  = \ (mva, mvb) ->  TheMVPair mva mvb
+--{-# INLINE mvPair #-}
+
+--mvUnPair  :: MVPair mv st  (a,b) -> (mv st a,mv st b)
+--mvUnPair = \ (TheMVPair mva mvb)-> (mva,mvb)
+--{-# INLINE mvUnPair #-}
+
+instance  (MV.MVector (MVProd (V.Mutable v) ('Pair pa pb )  ) (a,b) ,V.Vector (VProd v pa) a,V.Vector (VProd v pb) b)
+  => V.Vector (VProd v ('Pair pa pb )) (a,b)  where
+    {-# INLINE  basicUnsafeFreeze #-}
+    basicUnsafeFreeze = \(MVPair mva mvb) ->
+      VPair <$> V.basicUnsafeFreeze mva <*> V.basicUnsafeFreeze mvb
+
+    {-# INLINE basicUnsafeThaw #-}
+    basicUnsafeThaw = \(VPair va vb) ->
+      MVPair <$> V.basicUnsafeThaw va <*> V.basicUnsafeThaw vb
+
+    {-# INLINE basicLength #-}
+    basicLength = \(VPair va _) -> V.basicLength va
+
+    {-# INLINE basicUnsafeSlice #-}
+    basicUnsafeSlice = \start len (VPair va vb) ->
+      VPair (V.basicUnsafeSlice start len va) (V.basicUnsafeSlice start len vb)
+
+    {-# INLINE basicUnsafeIndexM #-}
+    basicUnsafeIndexM = \(VPair va vb) ix ->
+      do
+          a <- V.basicUnsafeIndexM va ix
+          b <- V.basicUnsafeIndexM vb ix
+          return (a,b)
+
+instance  (MV.MVector (MVProd (V.Mutable v) 'Unit  ) a ,V.Vector v a)
+  => V.Vector (VProd v 'Unit) a  where
+
+    {-# INLINE  basicUnsafeFreeze #-}
+    {-# INLINE basicUnsafeThaw #-}
+    {-# INLINE basicLength #-}
+    {-# INLINE basicUnsafeSlice #-}
+    {-# INLINE basicUnsafeIndexM #-}
+
+    basicUnsafeFreeze = \(MVLeaf mva) ->
+      VLeaf <$> V.basicUnsafeFreeze mva
+    basicUnsafeThaw = \(VLeaf va ) ->
+      MVLeaf <$> V.basicUnsafeThaw va
+    basicLength = \(VLeaf va ) -> V.basicLength va
+    basicUnsafeSlice = \start len (VLeaf va ) ->
+      VLeaf(V.basicUnsafeSlice start len va)
+    basicUnsafeIndexM = \(VLeaf va) ix ->  V.basicUnsafeIndexM va ix
+
+
+instance (MV.MVector mv a) => MV.MVector (MVProd mv 'Unit) a where
+  basicLength = \ (MVLeaf mva) -> MV.basicLength mva
+  {-# INLINE basicLength #-}
+
+
+  basicInitialize = \ (MVLeaf mva) -> MV.basicInitialize mva
+  {-# INLINE basicInitialize #-}
+
+  basicUnsafeSlice = \ start len (MVLeaf mva  )->
+    MVLeaf (MV.basicUnsafeSlice start len mva)
+  {-# INLINE basicUnsafeSlice #-}
+
+  basicOverlaps = \ (MVLeaf mva ) (MVLeaf mva2 )-> (MV.basicOverlaps mva mva2)
+  {-# INLINE basicOverlaps #-}
+
+  basicUnsafeNew =
+      \ size ->
+          MVLeaf <$> MV.basicUnsafeNew size
+  {-# INLINE basicUnsafeNew #-}
+
+  basicUnsafeReplicate =
+      \ size a ->
+         MVLeaf <$>
+            MV.basicUnsafeReplicate size a
+  {-# INLINE basicUnsafeReplicate #-}
+
+  basicUnsafeRead = \(MVLeaf mva ) ix ->   MV.basicUnsafeRead mva ix
+  {-#INLINE basicUnsafeRead #-}
+
+  basicUnsafeWrite = \ (MVLeaf mva ) ix a  ->
+    do
+      MV.basicUnsafeWrite mva ix a
+      return ()
+  {-#INLINE basicUnsafeWrite #-}
+
+  {-#INLINE basicUnsafeGrow #-}
+  basicUnsafeGrow = \ (MVLeaf mva ) growth ->
+      MVLeaf <$> MV.basicUnsafeGrow mva growth
+
+
+
+instance (MV.MVector (MVProd mv pra) a,MV.MVector (MVProd mv prb) b) => MV.MVector (MVProd mv ('Pair pra prb)) (a,b) where
+  basicLength = \ (MVPair mva _) -> MV.basicLength mva
+  {-# INLINE basicLength #-}
+
+
+  basicInitialize = \ (MVPair mva mvb) ->
+                        do  MV.basicInitialize mva ;
+                            MV.basicInitialize mvb
+  {-# INLINE basicInitialize #-}
+
+  basicUnsafeSlice = \ start len (MVPair mva mvb )->
+    MVPair (MV.basicUnsafeSlice start len mva) (MV.basicUnsafeSlice start len mvb)
+  {-# INLINE basicUnsafeSlice #-}
+
+  basicOverlaps = \ (MVPair mva mvb) (MVPair mva2 mvb2)-> (MV.basicOverlaps mva mva2) || (MV.basicOverlaps mvb mvb2)
+  {-# INLINE basicOverlaps #-}
+
+  basicUnsafeNew =
+      \ size ->
+          MVPair <$> MV.basicUnsafeNew size <*> MV.basicUnsafeNew size
+  {-# INLINE basicUnsafeNew #-}
+
+  basicUnsafeReplicate =
+      \ size (a,b) ->
+         MVPair <$>
+            MV.basicUnsafeReplicate size a <*>
+            MV.basicUnsafeReplicate size b
+  {-# INLINE basicUnsafeReplicate #-}
+
+  basicUnsafeRead = \(MVPair mva mvb) ix ->
+    (,) <$>  MV.basicUnsafeRead mva ix <*> MV.basicUnsafeRead mvb ix
+
+  {-# INLINE basicUnsafeRead #-}
+
+  basicUnsafeWrite = \ (MVPair mva mvb) ix (a,b) ->
+    do
+      MV.basicUnsafeWrite mva ix a
+      MV.basicUnsafeWrite mvb ix b
+      return ()
+  {-#INLINE basicUnsafeWrite #-}
+
+  {-#INLINE basicUnsafeGrow #-}
+  basicUnsafeGrow = \ (MVPair mva mvb) growth ->
+      MVPair <$> MV.basicUnsafeGrow mva growth <*>
+          MV.basicUnsafeGrow mvb growth
+
+
+
+
+
diff --git a/src/Numerical/InternalUtils.hs b/src/Numerical/InternalUtils.hs
new file mode 100644
--- /dev/null
+++ b/src/Numerical/InternalUtils.hs
@@ -0,0 +1,14 @@
+{-# LANGUAGE NoImplicitPrelude#-}
+module Numerical.InternalUtils(
+    error
+) where
+
+--import GHC.Stack (errorWithStackTrace,currentCallStack,whoCreated)
+import Prelude (error)
+
+
+
+{-
+note well: the stack traces only exist
+when doing a profiling build in GHC < 7.9/7.10
+-}
diff --git a/src/Numerical/Matrix/Basic.hs b/src/Numerical/Matrix/Basic.hs
new file mode 100644
--- /dev/null
+++ b/src/Numerical/Matrix/Basic.hs
@@ -0,0 +1,32 @@
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE ScopedTypeVariables#-}
+{-# LANGUAGE DeriveDataTypeable #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE ConstraintKinds #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FunctionalDependencies #-}
+
+
+module Numerical.Matrix.Basic where
+{-
+note, the contents of this module will probably be completely relocated elsewhere
+at some point
+
+-}
+--import Numerical.Array.Mutable as Mut
+
+
+--cleverDotProduct  :: (Address ~ MArrayAddress mvecL
+--                    , Address ~ MArrayAddress mvecR
+--                    ,Array mvecL (S Z) a
+--                    ,Array mvecR (S Z) a
+--                    ,Num a)=>
+
+--naiveDotProduct
+
+
diff --git a/src/Numerical/Nat.hs b/src/Numerical/Nat.hs
new file mode 100644
--- /dev/null
+++ b/src/Numerical/Nat.hs
@@ -0,0 +1,107 @@
+{-# LANGUAGE DataKinds, GADTs, TypeFamilies, TypeOperators,
+             ConstraintKinds, ScopedTypeVariables, RankNTypes #-}
+{-# LANGUAGE UndecidableInstances #-}
+{-# LANGUAGE StandaloneDeriving #-}
+{-# LANGUAGE DeriveDataTypeable#-}
+{-# LANGUAGE CPP #-}
+
+module Numerical.Nat(Nat(..),N0,N1,N2,N3,N4,N5,N6,N7,N8,N9,N10
+    ,SNat(..), type (+),plus_id_r,plus_succ_r,gcastWith,Proxy(..),LitNat,U)  where
+import Data.Typeable
+import Data.Data
+import qualified GHC.TypeLits as TL
+
+
+#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 707
+import Data.Type.Equality(gcastWith)
+#else
+import Data.Proxy
+#endif
+
+type LitNat = TL.Nat
+
+data Nat = S !Nat  | Z
+    deriving (Eq,Show,Read,Typeable,Data)
+
+#if defined(__GLASGOW_HASKELL__) && ( __GLASGOW_HASKELL__ >= 707) && ( __GLASGOW_HASKELL__ < 709)
+deriving instance Typeable 'Z
+deriving instance Typeable 'S
+#endif
+
+{-
+use closed type families when available,
+need to test that the
+-}
+
+
+type family U (n:: TL.Nat) :: Nat  where
+  U 0 = 'Z
+  U n = 'S (U (((TL.-)) n  1))
+
+
+type family n1 + n2 where
+  'Z + n2 = n2
+  ('S n1') + n2 = 'S (n1' + n2)
+
+
+--  ghc 7.6 instances
+
+--type family U (n:: (TL.Nat)) :: Nat
+
+---- can't induct, hence crippled
+--type instance U n = Z
+
+--type family (n1::Nat) + (n2::Nat) :: Nat
+--type instance Z + n2 = n2
+--type instance  (S n1) + n2 = S (n1 + n2)
+--gcastWith :: (a :~: b) -> ((a ~ b) => r) -> r
+--gcastWith Refl x = x
+--data a :~: b where
+--  Refl :: a :~: a
+
+
+
+
+
+-- singleton for Nat
+
+
+data SNat :: Nat -> * where
+  SZero :: SNat 'Z
+  SSucc :: SNat n -> SNat ('S n)
+
+
+
+-- inductive proof of right-identity of +
+plus_id_r :: SNat n -> ((n + 'Z) :~: n)
+plus_id_r SZero = Refl
+plus_id_r (SSucc n) = gcastWith (plus_id_r n) Refl
+
+-- inductive proof of simplification on the rhs of +
+plus_succ_r :: SNat n1 -> Proxy n2 -> ((n1 + ('S n2)) :~: ('S (n1 + n2)))
+plus_succ_r SZero _ = Refl
+plus_succ_r (SSucc n1) proxy_n2 = gcastWith (plus_succ_r n1 proxy_n2) Refl
+
+
+
+type N0 = 'Z
+
+type N1 = 'S N0
+
+type N2 = 'S N1
+
+type N3 = 'S N2
+
+type N4 = 'S N3
+
+type N5 = 'S N4
+
+type N6 = 'S N5
+
+type N7 = 'S N6
+
+type N8 = 'S N7
+
+type N9 = 'S N8
+
+type N10 = 'S N9
diff --git a/src/Numerical/World.hs b/src/Numerical/World.hs
new file mode 100644
--- /dev/null
+++ b/src/Numerical/World.hs
@@ -0,0 +1,18 @@
+
+
+
+module Numerical.World where 
+
+
+
+{-| 
+Every numerical algorithm runs somewhere.
+
+This could be on a CPU, a GPU,  
+
+-}
+-- Native is Just Haskell and Cbits, no external Deps
+data Native
+
+-- ForeignNative can have foreign lib deps, 
+data ForeignNative
diff --git a/tests/Main.hs b/tests/Main.hs
new file mode 100644
--- /dev/null
+++ b/tests/Main.hs
@@ -0,0 +1,32 @@
+module Main where
+
+
+import   NumericalUnit.Layout
+import   NumericalUnit.Shape
+
+
+
+import Data.List
+import Data.Ord
+
+import Test.Hspec
+import Control.Exception (evaluate)
+
+main :: IO ()
+main = hspec  $ do
+  describe "Shape Unit Tests" $  unitTestShape
+
+--main = defaultMain tests
+
+--tests :: Spec
+--tests = testGroup "Unit Tests" [unitTestShape] -- , unitTestLayout ]
+
+
+--unitTests = testGroup "Unit tests"
+--  [ testCase "List comparison (different length)" $
+--      [1, 2, 3] `compare` [1,2] @?= GT
+
+--  -- the following test does not hold
+--  , testCase "List comparison (same length)" $
+--      [1, 2, 3] `compare` [1,2,2] @?= LT
+--  ]
diff --git a/tests/NumericalUnit/Layout.hs b/tests/NumericalUnit/Layout.hs
new file mode 100644
--- /dev/null
+++ b/tests/NumericalUnit/Layout.hs
@@ -0,0 +1,7 @@
+module NumericalUnit.Layout(unitTestLayout) where 
+
+import Test.HUnit
+
+unitTestLayout =  [
+        
+        ]
diff --git a/tests/NumericalUnit/Shape.hs b/tests/NumericalUnit/Shape.hs
new file mode 100644
--- /dev/null
+++ b/tests/NumericalUnit/Shape.hs
@@ -0,0 +1,53 @@
+{-# LANGUAGE DataKinds, GADTs, TypeFamilies #-}
+{-# LANGUAGE ScopedTypeVariables  #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FunctionalDependencies #-}
+
+module NumericalUnit.Shape(unitTestShape) where
+
+
+import Numerical.Array.Shape as S
+import qualified Data.Vector.Storable as SV
+import qualified Data.Vector.Unboxed as UV
+import Prelude as P
+import Test.Hspec.Expectations
+import Test.Hspec
+
+unitTestShape :: Spec
+unitTestShape = describe "unit tests for Shape" $ do
+        specify "foldl on shape" $  S.foldl (+) 0 (1:* 2:* 3 :* Nil )  `shouldBe`   P.foldl   (+) 0  [1,2,3]
+        specify "foldr on shape" $  S.foldr (+) 0 (1:* 2:* 3 :* Nil )  `shouldBe`   P.foldr  (+) 0  [1,2,3]
+        specify "foldl1 on shape" $  S.foldl1 (+) (1:* 2:* 3 :* Nil )  `shouldBe`   P.foldl1   (+)  [1,2,3]
+        specify "foldr1 on shape" $  S.foldr1 (+) (1:* 2:* 3 :* Nil )  `shouldBe`   P.foldr1  (+)  [1,2,3]
+
+        specify "shapeToList on shape" $  S.shapeToList (1:* 2 :* 3 :* Nil) `shouldBe` [1,2,3]
+
+        specify "Show on Nil shape" $ show Nil `shouldBe` "Nil"
+        specify "Show on 1:* Nil"  $ show (1:* Nil) `shouldBe` "1 :* Nil"
+
+        specify "storable on size 0 shape" $
+          do a <- return (svFromList [Nil,Nil :: Shape Z Int]) ; SV.toList a `shouldBe` [Nil,Nil]
+        specify "storable on size 1 shape" $
+          do a <- return (svFromList [1:*Nil,2:*Nil :: Shape (S Z) Int]) ; SV.toList a `shouldBe` [1:*Nil,2:*Nil]
+        specify "storable on size 2 shape" $
+          do  a <- return (svFromList [3:* 4:* Nil,1:*2:*Nil :: Shape (S (S Z)) Int]) ;
+              SV.toList a `shouldBe` [3:* 4:* Nil,1:*2:*Nil]
+
+        specify "unboxed on size 0 shape" $
+          do a <- return (uvFromList [Nil,Nil :: Shape Z Int]) ; UV.toList a `shouldBe` [Nil,Nil]
+        specify "unboxed on size 1 shape" $
+          do a <- return (uvFromList [1:*Nil,2:*Nil :: Shape (S Z) Int]) ; UV.toList a `shouldBe` [1:*Nil,2:*Nil]
+        specify "unboxed on size 2 shape" $
+          do  a <- return (uvFromList [3:* 4:* Nil,1:*2:*Nil :: Shape (S (S Z)) Int]) ;
+                UV.toList a `shouldBe` [3:* 4:* Nil,1:*2:*Nil]
+
+    where
+        {- The NOINLINE is need to properly check storable /unboxed instances, otherwise fusion removes the allocation! -}
+        svFromList  :: SV.Storable a => [a] -> SV.Vector a
+        svFromList = SV.fromList
+        {-# NOINLINE svFromList #-}
+
+        uvFromList :: UV.Unbox a => [a] -> UV.Vector a
+        uvFromList = UV.fromList
+        {-# NOINLINE uvFromList#-}