diff --git a/LICENSE b/LICENSE
--- a/LICENSE
+++ b/LICENSE
@@ -2,7 +2,7 @@
 BSD3 Full Text:
 --------------------------------------------------------------------------------
 
-Copyright (c) 2011-2013, Joel Svensson 
+Copyright (c) 2011-2014, Joel Svensson 
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without
diff --git a/Obsidian.cabal b/Obsidian.cabal
--- a/Obsidian.cabal
+++ b/Obsidian.cabal
@@ -1,13 +1,14 @@
 Name:           Obsidian
-Version:        0.0.0.5
+Version:        0.1.0.0
 
+
 License:                BSD3
 License-file:           LICENSE
 Stability:              Beta
 Maintainer:		Joel Svensson<svenssonjoel@yahoo.se>
 Author:			Joel Svensson<svenssonjoel@yahoo.se>
 
-Copyright:              Copyright (c) 2011-2013 Joel Svensson 		  			
+Copyright:              Copyright (c) 2011-2014 Joel Svensson 		  			
   
 Synopsis:               Embedded language for GPU Programming 
 HomePage:               https://github.com/svenssonjoel/Obsidian
@@ -21,23 +22,30 @@
 
 build-type: Simple
 
---extra-source-files:
---     examples/tests/Test_DotProd.hs
-
-
 source-repository head
   type:     git
   location: git://github.com/svenssonjoel/Obsidian.git
 
 ----------------------------------------------------------------------------------------------------
 Library
+  -- Stable packages, no upper bounds
   build-depends: base >= 4 && < 5
-               , mtl >= 2.0 
-               , value-supply >= 0.6
-               , containers >= 0.4.2.1
+               , vector           >= 0.10.9.1 
+               , mtl              >= 2.0 
+               , containers       >= 0.4.2.1
+               , text             >= 0.11.3.1
+               , process          >= 1.1.0.2
+               , rdtsc            == 1.3.0.0  
+  -- Less-stable packages, upper bounds on next major version:
+  build-depends: value-supply     >= 0.6      && < 0.7
+               , language-c-quote >= 0.7.2    && < 0.8
+               , mainland-pretty  >= 0.2.6    && < 0.3
+               , mwc-random       >= 0.13.1.1 && < 0.14
+               , cuda             >= 0.5.1.1  && < 0.6
              
   exposed-modules: Obsidian
-                 
+                 ,  Obsidian.CodeGen.CUDA
+                 ,  Obsidian.Run.CUDA.Exec                    
             
   other-modules:  Obsidian.Array
                ,  Obsidian.Atomic
@@ -46,25 +54,20 @@
                ,  Obsidian.Force
                ,  Obsidian.Globs
                ,  Obsidian.Library
-               ,  Obsidian.LibraryG 
---               ,  Obsidian.Lift
                ,  Obsidian.Memory
+               ,  Obsidian.Mutable
                ,  Obsidian.Names
                ,  Obsidian.Program 
                ,  Obsidian.SeqLoop 
                ,  Obsidian.Types 
---	       ,  Obsidian.CodeGen.C
-               ,  Obsidian.CodeGen.CUDA
-               ,  Obsidian.CodeGen.Common
-               ,  Obsidian.CodeGen.InOut
+               ,  Obsidian.CodeGen.CompileIM
                ,  Obsidian.CodeGen.Liveness
                ,  Obsidian.CodeGen.Memory
---               ,  Obsidian.CodeGen.OpenCL
-               ,  Obsidian.CodeGen.PP 
                ,  Obsidian.CodeGen.Program
-               ,  Obsidian.CodeGen.SPMDC
+               ,  Obsidian.CodeGen.Reify 
+              
 
- 	       
+        
   GHC-Options: 
 -- -O2 
 
diff --git a/Obsidian.hs b/Obsidian.hs
--- a/Obsidian.hs
+++ b/Obsidian.hs
@@ -4,13 +4,14 @@
                  module Obsidian.Types, 
                  module Obsidian.Force, 
                  module Obsidian.Library,
-                 module Obsidian.LibraryG,
-                 module Obsidian.CodeGen.InOut,
+                 module Obsidian.CodeGen.Reify,
+                 module Obsidian.CodeGen.CompileIM,
                  module Obsidian.CodeGen.CUDA,
                  module Obsidian.Atomic, 
                  module Obsidian.SeqLoop, 
                  module Obsidian.Memory, 
-                 module Obsidian.Names ) where
+                 module Obsidian.Names,
+                 module Obsidian.Mutable) where
 
 
 
@@ -19,11 +20,13 @@
 import Obsidian.Types
 import Obsidian.Array
 import Obsidian.Library
-import Obsidian.LibraryG
 import Obsidian.Force
-import Obsidian.CodeGen.InOut
-import Obsidian.CodeGen.CUDA 
+import Obsidian.CodeGen.Reify
+import Obsidian.CodeGen.CompileIM
+import Obsidian.CodeGen.CUDA
 import Obsidian.Atomic
 import Obsidian.SeqLoop
 import Obsidian.Memory
 import Obsidian.Names
+import Obsidian.Mutable
+
diff --git a/Obsidian/Array.hs b/Obsidian/Array.hs
--- a/Obsidian/Array.hs
+++ b/Obsidian/Array.hs
@@ -1,25 +1,41 @@
-{-# LANGUAGE MultiParamTypeClasses,  
-             FlexibleInstances, FlexibleContexts,
-             GADTs, 
-             TypeFamilies,
-             RankNTypes #-} 
+{-# LANGUAGE MultiParamTypeClasses,
+             FlexibleInstances,
+             GADTs  #-} 
 
 {- Joel Svensson 2012
 
    Notes:
+    2014-04-08: Experimenting with API 
+    ---- OUTDATED ----
+    2013-08-26: Experimenting with warp programs.
+                These do not fit that well in established Idioms!
+                TODO: Improve this situation. 
+    ---- OUTDATED ----
     2013-01-08: Removed number-of-blocks field from Distribs
     2012-12-10: Drastically shortened. 
 -}
 
-module Obsidian.Array  where
+module Obsidian.Array (Pull, Push, SPull, DPull, SPush, DPush,
+                       Pushable, 
+                       mkPull,
+                       mkPush,
+                       push,
+                       setSize,
+                       (!),
+                       (<:),
+                       Array(..),
+                       ArrayLength(..),
+                       ASize(..),
+                       namedGlobal,
+                       undefinedGlobal) where
 
 import Obsidian.Exp 
 import Obsidian.Types
 import Obsidian.Globs
 import Obsidian.Program
-import Obsidian.Names
 
-import Data.List
+import Prelude hiding (replicate) 
+import Data.List hiding (replicate) 
 import Data.Word
 
 ---------------------------------------------------------------------------
@@ -31,15 +47,18 @@
 type SPush t a = Push t Word32 a
 type DPush t a = Push t EWord32 a 
 ---------------------------------------------------------------------------
--- Create arrats
+-- Create arrays
 ---------------------------------------------------------------------------
+-- | An undefined array. Use as placeholder when generating code
 undefinedGlobal n = Pull n $ \gix -> undefined
+-- | A named global array. 
 namedGlobal name n = Pull n $ \gix -> index name gix
-namedPull name n = Pull n $ \gix -> index name gix
+-- namedPull name n = Pull n $ \gix -> index name gix
 
 ---------------------------------------------------------------------------
 -- Class ArraySize
---------------------------------------------------------------------------- 
+---------------------------------------------------------------------------
+-- | ASize provides conversion to Exp Word32 for array sizes
 class (Integral a, Num a) => ASize a where
   sizeConv :: a ->  Exp Word32
 
@@ -52,128 +71,171 @@
 ---------------------------------------------------------------------------
 -- Push and Pull arrays
 ---------------------------------------------------------------------------
+-- | Push array. Parameterised over Program type and size type.
 data Push p s a =
-  Push s ((a -> Exp Word32 -> TProgram NameInfo) -> Program p NameInfo)
+  Push s ((a -> EWord32 -> TProgram ()) -> Program p ())
 
+-- | Pull array.
 data Pull s a = Pull {pullLen :: s, 
                       pullFun :: EWord32 -> a}
 
-mkPushArray :: s -> ((a -> EWord32 -> TProgram NameInfo)
-                             -> Program t NameInfo) -> Push t s a
-mkPushArray n p = Push n p 
-mkPullArray n p = Pull n p  
+-- | Create a push array. 
+mkPush :: s
+       -> ((a -> EWord32 -> TProgram ()) -> Program t ())
+       -> Push t s a
+mkPush n p = Push n p 
 
+-- | Create a pull array. 
+mkPull n p = Pull n p 
+
+-- Fix this.
+--   * you cannot safely resize either push or pull arrays
+--   * you can shorten pull arrays safely.  
+setSize :: l -> Pull l a -> Pull l a
+setSize n (Pull _ ixf) = mkPull n ixf
+
+---------------------------------------------------------------------------
+-- Array Class 
+---------------------------------------------------------------------------
+class ArrayLength a where
+  -- | Get the length of an array.
+  len :: a s e -> s
+
+instance ArrayLength Pull where
+  len    (Pull n ixf) = n
+
+instance ArrayLength (Push t) where
+  len  (Push s p) = s
+
 class Array a where
-  resize :: r -> a s e -> a r e
-  len    :: ASize s => a s e -> s
-  aMap   :: (e -> e') -> a s e -> a s e'
-  ixMap  :: (Exp Word32 -> Exp Word32)
-            -> a s e -> a s e
+  -- | Array of consecutive integers
+  iota      :: ASize s => s -> a s EWord32
+  -- | Create an array by replicating an element. 
+  replicate :: ASize s => s -> e -> a s e 
+
+  -- | Map a function over an array. 
+  aMap      :: (e -> e') -> a s e -> a s e'
+  -- | Perform arbitrary permutations (dangerous). 
+  ixMap     :: (EWord32 -> EWord32)
+               -> a s e -> a s e
+  -- -- | Reduce an array using a provided operator. 
+  -- fold1     :: (e -> e -> e) -> a Word32 e -> a Word32 e  
+
+  -- would require Choice !
+  -- | Append two arrays. 
+  append    :: (ASize s, Choice e) => a s e -> a s e -> a s e 
   
-instance Array Pull where 
-  resize m (Pull _ ixf) = Pull m ixf
-  len      (Pull s _)   = s
+  -- technicalities
+  -- | Statically sized array to dynamically sized array.
+  toDyn     :: a Word32 e -> a EW32 e
+  -- | Dynamically sized array to statically sized array. 
+  fromDyn   :: Word32 -> a EW32 e -> a Word32 e 
+  
+instance Array Pull where
+  iota   s = Pull s $ \ix -> ix 
+  replicate s e = Pull s $ \_ -> e
+
   aMap   f (Pull n ixf) = Pull n (f . ixf)
   ixMap  f (Pull n ixf) = Pull n (ixf . f) 
+
+  append a1 a2 = Pull (n1+n2)
+               $ \ix -> ifThenElse (ix <* (sizeConv n1)) 
+                       (a1 ! ix) 
+                       (a2 ! (ix - (sizeConv n1)))
+    where 
+      n1 = len a1
+      n2 = len a2 
+
+  -- technicalities
+  toDyn (Pull n ixf) = Pull (fromIntegral n) ixf
+  fromDyn n (Pull _ ixf) = Pull n ixf 
+   
   
-instance Array (Push t) where 
-  resize m (Push _ p) = Push m p
-  len      (Push s _) = s
+instance Array (Push t) where
+  iota s = Push s $ \wf ->
+    do
+      forAll (sizeConv s) $ \ix -> wf ix ix 
+  replicate s e = Push s $ \wf ->
+    do
+      forAll (sizeConv s) $ \ix -> wf e ix 
   aMap   f (Push s p) = Push s $ \wf -> p (\e ix -> wf (f e) ix)
   ixMap  f (Push s p) = Push s $ \wf -> p (\e ix -> wf e (f ix))
-  
 
-class Indexible a where 
-  access :: a s e -> Exp Word32 -> e 
-  
-instance Indexible Pull where
-  access p ix = pullFun p ix
+  -- unfortunately a Choice constraint. 
+  append p1 p2  =
+    Push (n1 + n2) $ \wf ->
+      do p1 <: wf
+         p2 <: \a i -> wf a (sizeConv n1 + i) 
+           where 
+             n1 = len p1
+             n2 = len p2 
 
+   -- technicalities
+  toDyn (Push n p) = Push (fromIntegral n) p 
+  fromDyn n (Push _ p) = Push n p 
+ 
+
 ---------------------------------------------------------------------------
 -- Functor instance Pull/Push arrays
 ---------------------------------------------------------------------------
 instance Array arr => Functor (arr w) where 
   fmap = aMap
 
-
 ---------------------------------------------------------------------------
 -- Pushable
 ---------------------------------------------------------------------------
-class Pushable a where 
-  push  :: ASize s => PT t -> a s e -> Push t s e
-  -- Push using m threads
-  --  m must be a divisor of nm  (TODO: error otherwise) 
-  --pushN :: Word32 -> a e -> Push e
+class Pushable t where
+  push :: ASize s => Pull s e -> Push t s e 
 
-  -- push grouped elements to adjacent indices using
-  -- one thread per group. 
-  --pushF :: a [e] -> Push e 
- 
-instance Pushable (Push Thread) where 
-  push Thread = id
-  push Block  = error "not implemented: Program transformation!"
-  push Grid   = error "not implemented: Program transformation!" 
-instance Pushable (Push Block) where 
-  push Block = id 
-  push Thread = error "not implemented: Program transformations!"
-  push Grid = error "not implemented: Program transformations!" 
-instance Pushable (Push Grid) where 
-  push Grid = id
-  push Thread = error "not implemented: Program transformations!"
-  push Block = error "not implemented: Program transformations!" 
-  
-instance Pushable Pull where
-  push Thread (Pull n ixf) =
+instance Pushable Thread where
+  push (Pull n ixf) =
     Push n $ \wf -> seqFor (sizeConv n) $ \i -> wf (ixf i) i
-  push Block (Pull n ixf) =
-    Push n $ \wf -> ForAll (sizeConv n) $ \i -> wf (ixf i) i 
-  push Grid (Pull n ixf) =
-    Push n $ \wf -> ForAllThreads (sizeConv n) $ \i -> wf (ixf i) i 
 
-{-  pushN m (Pull nm ixf) =
-    Push nm -- There are still nm elements (info for alloc) 
-    $ \wf ->
-    ForAll (Just m) 
-    $ \i ->
-    sequence_ [wf (ixf (i + fromIntegral (j * n))) (i + (fromIntegral (j * n)))
-              | j <- [0..n]] 
-    -- Force can still Allocate n elements for this Push array.
-    where
-      n = fromIntegral (nm `div` m)
-
-  pushF (Pull n ixf) =
-    Push (n * m) $ \wf ->
-    ForAll (Just n) $ \i ->
-    sequence_ [wf ((ixf i) !! fromIntegral j)  (i * fromIntegral m + fromIntegral j)
-              | j <- [0..m]]
-    where 
-      m = fromIntegral$ length $ ixf 0
--} 
-
-{- ------------------------------------------------------------------------
-
-     m     m     m     m     m 
-  |-----|-----|-----|-----|-----| (n * m)
-
-   01234 01234 01234 01234 01234  k 
+instance Pushable Warp where
+  push (Pull n ixf) =
+    Push n $ \wf ->
+      forAll (sizeConv n) $ \i -> wf (ixf i) i
 
-  0     1     2     3     4       j
+instance Pushable Block where
+  push (Pull n ixf) =
+    Push n $ \wf ->
+      forAll (sizeConv n) $ \i -> wf (ixf i) i
 
-  m threads, each writing n elements:
-  [(tid + (j*m) | j <- [0..n]] 
+instance Pushable Grid where
+  push (Pull n ixf) =
+    Push n $ \wf ->
+      forAll (sizeConv n) $ \i -> wf (ixf i) i 
+  
+-- class PushableN t where
+--   pushN :: ASize s => Word32 -> Pull s e -> Push t s e
 
-  n threads, each writing m elements:
-  [(tid * m + k | k <- [0..m]] 
+-- instance PushableN Block where
+--   pushN n (Pull m ixf) =
+--     Push m $ \ wf -> forAll (sizeConv (m `div` fromIntegral n)) $ \tix ->
+--     warpForAll 1 $ \_ -> 
+--     seqFor (fromIntegral n) $ \ix -> wf (ixf (tix * fromIntegral n + ix))
+--                                              (tix * fromIntegral n + ix) 
+ 
+-- instance PushableN Grid where
+--   pushN n (Pull m ixf) =
+--     Push m $ \ wf -> forAll (sizeConv (m `div` fromIntegral n)) $ \bix ->
+--     forAll (fromIntegral n) $ \tix -> wf (ixf (bix * fromIntegral n + tix))
+--                                               (bix * fromIntegral n + tix) 
 
------------------------------------------------------------------------- -}   
----------------------------------------------------------------------------
+--------------------------------------------------------------------------
 -- Indexing, array creation.
 ---------------------------------------------------------------------------
-namedArray name n = mkPullArray n (\ix -> index name ix)
-indexArray n      = mkPullArray n (\ix -> ix)
 
-pushApp (Push _ p) a = p a 
+pushApp (Push _ p) a = p a
 
+infixl 9 <:
+(<:) :: Push t s a
+        -> (a -> EWord32 -> Program Thread ())
+        -> Program t ()
+(<:) = pushApp 
+
 infixl 9 ! 
-(!) :: Indexible a => a s e -> Exp Word32 -> e 
-(!) = access
+(!) :: Pull s e -> Exp Word32 -> e 
+(!) arr = pullFun arr 
+
+
diff --git a/Obsidian/Atomic.hs b/Obsidian/Atomic.hs
--- a/Obsidian/Atomic.hs
+++ b/Obsidian/Atomic.hs
@@ -7,8 +7,31 @@
 
 import Obsidian.Exp
 import Data.Word
-       
+import Data.Int
 
+
+-- Anyone can extend these with new instances.
+-- Not good. (I need to think about how to separate
+-- low level CUDA concerns, from programmer level concerns) 
+class Scalar a => AtomicInc a
+instance AtomicInc Word32
+
+class Scalar a => AtomicAdd a
+instance AtomicAdd Word32
+instance AtomicAdd Int32
+instance AtomicAdd Word64
+
+class Scalar a => AtomicSub a
+instance AtomicSub Word32
+instance AtomicSub Int32
+
+class Scalar a => AtomicExch a
+instance AtomicExch Word32
+instance AtomicExch Word64
+instance AtomicExch Int32
+
+
+
 ---------------------------------------------------------------------------
 -- Atomic operations 
 ---------------------------------------------------------------------------
@@ -16,6 +39,12 @@
 
   -- Cuda only allows AtomicInc on the Int type
   --  (todo: figure out if CUDA int is 32 or 64 bit) 
-  AtomicInc :: Atomic Word32 
+  AtomicInc :: AtomicInc a => Atomic a
 
+  AtomicAdd :: AtomicAdd a => Exp a -> Atomic a
+
+  AtomicSub :: AtomicSub a => Exp a -> Atomic a 
+
+  AtomicExch :: AtomicExch a => Exp a -> Atomic a
+  
 printAtomic AtomicInc = "atomicInc"
diff --git a/Obsidian/CodeGen/CUDA.hs b/Obsidian/CodeGen/CUDA.hs
--- a/Obsidian/CodeGen/CUDA.hs
+++ b/Obsidian/CodeGen/CUDA.hs
@@ -1,209 +1,30 @@
-{- Joel Svensson 2012, 2013 -}
 
-{-# LANGUAGE GADTs #-} 
-module Obsidian.CodeGen.CUDA 
-       (genKernel) where  
 
-import Data.List
-import Data.Word 
-import Data.Monoid
-import qualified Data.Map as Map
-import Control.Monad.State
+module Obsidian.CodeGen.CUDA (genKernel) where 
 
-import Obsidian.Array
-import Obsidian.Exp 
 
-import Obsidian.Types
-import Obsidian.Globs
-import Obsidian.Atomic 
-
-import Obsidian.CodeGen.PP
-import Obsidian.CodeGen.Common
-import Obsidian.CodeGen.InOut 
-import Obsidian.CodeGen.Memory
+import Obsidian.CodeGen.Reify
+import Obsidian.CodeGen.CompileIM
 import Obsidian.CodeGen.Liveness
-
--- New imports
-import Obsidian.CodeGen.Program 
-import qualified Obsidian.Program as P 
-
-import Obsidian.CodeGen.SPMDC
-
----------------------------------------------------------------------------
--- a gc
---------------------------------------------------------------------------- 
-gc = genConfig "" ""
-
+import Obsidian.CodeGen.Memory
+import Text.PrettyPrint.Mainland
 
+import qualified Data.Map as M
+import Data.Word
 ---------------------------------------------------------------------------
--- C style function "header"
+-- Generate CUDA kernels
 ---------------------------------------------------------------------------
 
-kernelHead :: Name -> 
-              [(String,Type)] -> 
-              [(String,Type)] -> 
-              PP () 
-kernelHead name ins outs = 
-  do 
-    line ("__global__ void " ++ name ++ "(" ++ types ++ ")" )   
-  where 
-    types = concat (intersperse "," (typeList (ins ++ outs)))
-    typeList :: [(String,Type)] -> [String] 
-    typeList []              = [] 
-    typeList ((a,t):xs)      = (genType gc t ++ a) : typeList xs
- 
----------------------------------------------------------------------------
--- genKernel 
----------------------------------------------------------------------------
-    
-genKernel :: ToProgram a b => String -> (a -> b) -> Ips a b -> String 
-genKernel name kernel a = proto ++ ts ++ cuda 
-  where
-    (ins,im) = toProgram 0 kernel a
-
-    outs = getOutputs im
-    
-    lc  = computeLiveness im 
-    
-    -- Creates (name -> memory address) map      
-    (m,mm) = mmIM lc sharedMem Map.empty
-             
-    -- What if its Right ??? (I DONT KNOW!) 
-    (Left threadBudget) = numThreads im
-    ts = "/* number of threads needed " ++ show threadBudget ++ "*/\n"
-
-    spmd = imToSPMDC threadBudget im
-    
-    
-    body' = (if size m > 0 then (shared :) else id)  $ mmSPMDC mm spmd
-
-    em = snd $ execState (collectExps body') ( 0, Map.empty)
-    (decls,body'') = replacePass em body'
-    spdecls = declsToSPMDC decls 
-
-    body = spdecls ++ body''
-              
-    swap (x,y) = (y,x)
-    inputs = map ((\(t,n) -> (typeToCType t,n)) . swap) ins
-    outputs = map ((\(t,n) -> (typeToCType t,n)) . swap) outs 
-    
-    ckernel = CKernel CQualifyerKernel CVoid name (inputs++outputs) body
-    shared = CDecl (CQualified CQualifyerExtern (CQualified CQualifyerShared ((CQualified (CQualifyerAttrib (CAttribAligned 16)) (CArray []  (CWord8)))))) "sbase"
-
-    proto = getProto name ins outs 
-    cuda = printCKernel (PPConfig "__global__" "" "" "__syncthreads()") ckernel 
-
-
----------------------------------------------------------------------------
--- Generate a function prototype
---------------------------------------------------------------------------- 
-getProto :: Name -> [(String,Type)] -> [(String,Type)] -> String
-getProto name ins outs =
-  runPP (
-    do 
-      line "extern \"C\" "
-      kernelHead name ins outs
-      line ";"
-      newline) 0 
-
----------------------------------------------------------------------------
--- generate a sbase CExpr
----------------------------------------------------------------------------
-sbaseCExpr 0    = cVar "sbase" (CPointer CWord8) 
-sbaseCExpr addr = cBinOp CAdd (cVar "sbase" (CPointer CWord8)) 
-                              (cLiteral (Word32Val addr) CWord32) 
-                              (CPointer CWord8) 
----------------------------------------------------------------------------
--- Memory map the arrays in an SPMDC
----------------------------------------------------------------------------
-mmSPMDC :: MemMap -> [SPMDC] -> [SPMDC] 
-mmSPMDC mm [] = [] 
-mmSPMDC mm (x:xs) = mmSPMDC' mm x : mmSPMDC mm xs
-
-mmSPMDC' :: MemMap -> SPMDC -> SPMDC
-mmSPMDC' mm (CAssign e1 es e2) = 
-  cAssign (mmCExpr mm e1) 
-          (map (mmCExpr mm) es)    
-          (mmCExpr mm e2)
-mmSPMDC' mm (CAtomic op e1 e2 e3) = cAtomic op (mmCExpr mm e1)
-                                               (mmCExpr mm e2)
-                                               (mmCExpr mm e3) 
-mmSPMDC' mm (CFunc name es) = cFunc name (map (mmCExpr mm) es) 
-mmSPMDC' mm CSync           = CSync
-mmSPMDC' mm (CIf   e s1 s2) = cIf (mmCExpr mm e) (mmSPMDC mm s1) (mmSPMDC mm s2)
-mmSPMDC' mm (CFor name e s) = cFor name (mmCExpr mm e) (mmSPMDC mm s)
-mmSPMDC' mm (CDeclAssign t nom e) = cDeclAssign t nom (mmCExpr mm e)
-mmSPMDC' mm a@(CDecl t nom) = a
-mmSPMDC' mm a = error $ "mmSPMDC': " ++ show a
----------------------------------------------------------------------------
--- Memory map the arrays in an CExpr
----------------------------------------------------------------------------
-mmCExpr mm (CExpr (CVar nom t)) =  
-  case Map.lookup nom mm of 
-    Just (addr,t) -> 
-      let core = sbaseCExpr addr 
-          cast c = cCast  c (typeToCType t)
-      in cast core
-    
-    Nothing -> cVar nom t
-mmCExpr mm (CExpr (CIndex (e1,es) t)) = cIndex (mmCExpr mm e1, map (mmCExpr mm) es) t
-mmCExpr mm (CExpr (CBinOp op e1 e2 t)) = cBinOp op (mmCExpr mm e1) (mmCExpr mm e2) t
-mmCExpr mm (CExpr (CUnOp op e t)) = cUnOp op (mmCExpr mm e) t 
-mmCExpr mm (CExpr (CFuncExpr nom exprs t)) = cFuncExpr nom (map (mmCExpr mm) exprs) t
-mmCExpr mm (CExpr (CCast e t)) = cCast (mmCExpr mm e) t
-mmCExpr mm (CExpr (CCond e1 e2 e3 t)) = cCond (mmCExpr mm e1)
-                                              (mmCExpr mm e2)
-                                              (mmCExpr mm e3)
-                                              t
-mmCExpr mm a = a 
-          
-  
----------------------------------------------------------------------------
--- New IM to SPCMD
----------------------------------------------------------------------------
-atomicOpToCAtomicOp AtomicInc = CAtomicInc
-
-imToSPMDC :: Word32 -> IMList a -> [SPMDC]
-imToSPMDC nt im = concatMap (process nt) im
+genKernel :: ToProgram prg
+             => Word32
+             -> String
+             -> prg
+             -> String
+genKernel nt kn prg = prgStr
   where
-    process nt (SAssign name [] e,_) =
-      [cAssign (cVar name (typeToCType (typeOf e))) [] (expToCExp e)]
-
-    process nt (SAssign name [ix] e,_) = 
-      [cAssign (cVar name (typeToCType (Pointer (typeOf e)))) [expToCExp ix] (expToCExp e)]
-
-    process nt (SAtomicOp res arr e op,_) = 
-      [cAtomic (atomicOpToCAtomicOp op)
-               (cVar res (typeToCType (typeOf e)))
-               (cVar arr (typeToCType (Pointer (typeOf e))))
-               (expToCExp e)]
-
-    process nt (SCond bexp im,_) =
-      [cIf (expToCExp bexp) (imToSPMDC nt im) []]
-
-    process nt (SSeqFor name e im,_) =
-      [cFor name (expToCExp e) (imToSPMDC nt im)]
-
-    process nt (SForAll (Literal n) im,_) =
-      if (n < nt) 
-      then 
-        [cIf (cBinOp CLt (cThreadIdx X)  (cLiteral (Word32Val n) CWord32) CInt)
-         code []]
-      else 
-        code 
-      where 
-        code = imToSPMDC nt im
-
-    -- This one is tricky (since no corresponding CUDA construct exists) 
-    process nt (SForAllBlocks n im,_) =
-      -- TODO: there should be "number of blocks"-related conditionals here (possibly) 
-      imToSPMDC nt im
-    -- This one is even more tricky
-    process nt (SForAllThreads n im,_) =
-      imToSPMDC nt im 
-    process nt (SAllocate name size t,_) = []
-    process nt (SDeclare name t,_) =
-      [cDecl (typeToCType t) name]
-    process nt (SOutput name t,_) = [] -- RIGHT!
-    process nt (SSynchronize,_)   = [CSync]
-    
+    prgStr = pretty 75 $ ppr $ compile PlatformCUDA (Config nt bytesShared) kn (a,rim) 
+    (a,im) = toProgram_ 0 prg
+    iml = computeLiveness im
+    (m,mm) = mmIM iml sharedMem (M.empty)
+    bytesShared = size m 
+    rim = renameIM mm iml
diff --git a/Obsidian/CodeGen/Common.hs b/Obsidian/CodeGen/Common.hs
deleted file mode 100644
--- a/Obsidian/CodeGen/Common.hs
+++ /dev/null
@@ -1,176 +0,0 @@
-{-# LANGUAGE RankNTypes, GADTs  #-}
-
-{- Joel Svensson 2012 -} 
-module Obsidian.CodeGen.Common where 
-
-import Data.List
-import Data.Word
-import qualified Data.Map as Map 
-
-
-import Obsidian.Exp 
-import Obsidian.Types
-import Obsidian.Globs
-
-import Obsidian.CodeGen.PP
-import Obsidian.CodeGen.Memory
-
----------------------------------------------------------------------------
-data GenConfig = GenConfig { global :: String,
-                             local  :: String };
-  
-genConfig = GenConfig
-
-
----------------------------------------------------------------------------
--- Helpers
-
-mappedName :: Name -> Bool 
-mappedName name = isPrefixOf "arr" name
-
-tid :: Exp Word32
-tid = ThreadIdx X
-
-genType _ Int = "int "
-genType _ Int8 = "int8_t "
-genType _ Int16 = "int16_t "
-genType _ Int32 = "int32_t "
-genType _ Int64 = "int64_t "
-genType _ Float = "float "
-genType _ Double = "double "
-genType _ Bool = "int " 
-genType _ Word8 = "uint8_t "
-genType _ Word16 = "uint16_t "
-genType _ Word32 = "uint32_t "
-genType _ Word64 = "uint64_t "
-
-genType gc (Pointer t) = genType gc t ++ "*"
-genType gc (Global t) = global gc ++" "++ genType gc t  -- "__global " ++ genType t
-genType gc (Local t)  = local gc  ++" "++ genType gc t 
-
-genCast gc t = "(" ++ genType gc t ++ ")"
-
-parens s = '(' : s ++ ")"
- 
----------------------------------------------------------------------------
--- genExp C-style 
-genExp :: Scalar a => GenConfig -> MemMap -> Exp a -> [String]
-
--- Cheat and do CUDA printing here as well
-genExp gc _ (BlockDim X) = ["blockDim.x"]
-genExp gc _ (BlockIdx X) = ["blockIdx.x"]
-genExp gc _ (BlockIdx Y) = ["blockIdx.y"]
-genExp gc _ (BlockIdx Z) = ["blockIdx.z"]
-genExp gc _ (ThreadIdx X) = ["threadIdx.x"]
-genExp gc _ (ThreadIdx Y) = ["threadIdx.y"]
-genExp gc _ (ThreadIdx Z) = ["threadIdx.z"]
-
-
-genExp gc _ (Literal a) = [show a] 
-genExp gc _ (Index (name,[])) = [name]
-genExp gc mm exp@(Index (name,es)) = 
-  [name' ++ genIndices gc mm es]
-  where 
-    (offs,t)  = 
-      case Map.lookup name mm of  
-        Nothing -> error "array does not excist in map" 
-        (Just x) -> x
-    name' = if mappedName name 
-            then parens$ genCast gc t ++ 
-                 if offs > 0 
-                 then "(sbase+" ++ show offs ++ ")"             
-                 else "sbase"
-            else name
-
-   
-genExp gc mm (BinOp op e1 e2) = 
-  [genOp op (genExp gc mm e1 ++ genExp gc mm e2)]
-
-genExp gc mm (UnOp op e) = 
-  [genOp op (genExp gc mm e)] 
-  
-genExp gc mm (If b e1 e2) =   
-  [genIf (genExp gc mm b ++ 
-          genExp gc mm e1 ++ 
-          genExp gc mm e2 )] 
-
----------------------------------------------------------------------------
---
-genIndices gc mm es = concatMap (pIndex mm) es  
-  where 
-    pIndex mm e = "[" ++ concat (genExp gc mm e) ++ "]"
-
-
-genIf         [b,e1,e2] = "(" ++ b ++ " ? " ++ e1 ++ " : " ++ e2 ++ ")"
-
----------------------------------------------------------------------------
--- genOp
-genOp :: Op a -> [String] -> String
-genOp Add     [a,b] = oper "+" a b 
-genOp Sub     [a,b] = oper "-" a b 
-genOp Mul     [a,b] = oper "*" a b 
-genOp Div     [a,b] = oper "/" a b 
-
-genOp Mod     [a,b] = oper "%" a b 
-
-genOp Sin     [a]   = func "sin" a 
-genOp Cos     [a]   = func "cos" a 
--- Bool ops
-genOp Eq      [a,b] = oper "==" a b 
-genOp Lt      [a,b] = oper "<" a b 
-genOp LEq     [a,b] = oper "<=" a b 
-genOp Gt      [a,b] = oper ">" a b
-genOp GEq     [a,b] = oper ">=" a b
-
--- Bitwise ops
-genOp BitwiseAnd [a,b] = oper "&" a b 
-genOp BitwiseOr  [a,b] = oper "|" a b 
-genOp BitwiseXor [a,b] = oper "^" a b 
-genOp BitwiseNeg [a]   = unOp "~" a 
-genOp ShiftL     [a,b] = oper "<<" a b 
-genOp ShiftR     [a,b] = oper ">>" a b 
-
-
--- built-ins 
-genOp Min      [a,b] = func "min" (a ++ "," ++ b) 
-genOp Max      [a,b] = func "max" (a ++ "," ++ b) 
-
-genOp Int32ToWord32 [a]  = func "(uint32_t)" a
-genOp Word32ToInt32 [a]  = func "(int32_t)" a 
-
-func  f a = f ++ "(" ++ a ++ ")" 
-oper  f a b = "(" ++ a ++ f ++ b ++ ")" 
-unOp  f a   = "(" ++ f ++ a ++ ")"
-
----------------------------------------------------------------------------
--- Configurations, threads,memorymap 
-
-data Config = Config {configThreads  :: NumThreads, 
-                      configMM       :: MemMap,
-                      configLocalMem :: Word32} 
-config = Config
-
-
-assign :: Scalar a => GenConfig -> MemMap -> Exp a -> Exp a -> PP () 
-assign gc mm name val = line ((concat (genExp gc mm name)) ++ 
-                           " = " ++  concat (genExp gc mm val) ++ 
-                           ";") 
-                                                    
-cond :: GenConfig -> MemMap -> Exp Bool -> PP ()  
-cond gc mm e = line ("if " ++ concat (genExp gc mm e))  
-
-
-
--- used in both OpenCL and CUDA generation
-potentialCond gc mm n nt pp
-  | n < nt = 
-    do
-      cond gc mm (tid <* (fromIntegral n))
-      begin
-      pp       
-      end 
-  | n == nt = pp
-              
-  | otherwise = error "potentialCond: should not happen"
-
-
diff --git a/Obsidian/CodeGen/CompileIM.hs b/Obsidian/CodeGen/CompileIM.hs
new file mode 100644
--- /dev/null
+++ b/Obsidian/CodeGen/CompileIM.hs
@@ -0,0 +1,565 @@
+
+{-# LANGUAGE QuasiQuotes #-}
+{-# LANGUAGE PackageImports #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-} 
+
+{-
+
+   Joel Svensson 2013
+
+
+-} 
+
+module Obsidian.CodeGen.CompileIM where 
+
+import Language.C.Quote hiding (Block)
+import Language.C.Quote.CUDA as CU
+
+import qualified Language.C.Quote.OpenCL as CL 
+
+import qualified "language-c-quote" Language.C.Syntax as C 
+
+import Obsidian.Exp (IExp(..),IBinOp(..),IUnOp(..))
+import Obsidian.Types
+import Obsidian.DimSpec 
+import Obsidian.CodeGen.Program
+
+import Control.Monad.State
+
+import Data.Word
+import Data.Int
+
+{- TODOs:
+    
+   * Pass a target "platform" to code generator.
+      - CUDA
+      - OpenCL
+      - Sequential C
+      - C with OpenMP ? 
+   * rewrite some functions here to use  a reader monad. 
+
+
+
+   * TODO: Make sure tid always has correct Value 
+-} 
+
+---------------------------------------------------------------------------
+-- Platform
+---------------------------------------------------------------------------
+data Platform = PlatformCUDA
+              | PlatformOpenCL
+              | PlatformC
+
+data Config = Config { configThreadsPerBlock :: Word32,
+                       configSharedMem :: Word32}
+
+
+
+
+---------------------------------------------------------------------------
+-- compileExp (maybe a bad name)
+---------------------------------------------------------------------------
+compileExp :: IExp -> Exp 
+compileExp (IVar name t) = [cexp| $id:name |]
+
+
+-- TODO: Fix all this! 
+compileExp (IBlockIdx X) = [cexp| $id:("bid")|] -- [cexp| $id:("blockIdx.x") |]
+compileExp (IBlockIdx Y) = [cexp| $id:("blockIdx.y") |]
+compileExp (IBlockIdx Z) = [cexp| $id:("blockIdx.z") |]
+
+compileExp (IThreadIdx X) = [cexp| $id:("threadIdx.x") |]
+compileExp (IThreadIdx Y) = [cexp| $id:("threadIdx.y") |]
+compileExp (IThreadIdx Z) = [cexp| $id:("threadIdx.z") |]
+
+compileExp (IBlockDim X) = [cexp| $id:("blockDim.x") |]
+compileExp (IBlockDim Y) = [cexp| $id:("blockDim.y") |]
+compileExp (IBlockDim Z) = [cexp| $id:("blockDim.z") |]
+
+compileExp (IGridDim X) = [cexp| $id:("GridDim.x") |]
+compileExp (IGridDim Y) = [cexp| $id:("GridDim.y") |]
+compileExp (IGridDim Z) = [cexp| $id:("GridDim.z") |]
+
+compileExp (IBool True) = [cexp|1|]
+compileExp (IBool False) = [cexp|0|]
+compileExp (IInt8 n) = [cexp| $int:(toInteger n) |]
+compileExp (IInt16 n) = [cexp| $int:(toInteger n) |]
+compileExp (IInt32 n) = [cexp| $int:(toInteger n) |]
+compileExp (IInt64 n) = [cexp| $lint:(toInteger n) |]
+
+compileExp (IWord8 n) = [cexp| $uint:(toInteger n) |]
+compileExp (IWord16 n) = [cexp| $uint:(toInteger n) |]
+compileExp (IWord32 n) = [cexp| $uint:(toInteger n) |]
+compileExp (IWord64 n) = [cexp| $ulint:(toInteger n) |]
+
+compileExp (IFloat n) = [cexp| $float:(toRational n) |]
+compileExp (IDouble n) = [cexp| $double:(toRational n) |]
+
+compileExp (IIndex (i1,[e]) t) = [cexp| $(compileExp i1)[$(compileExp e)] |] 
+
+compileExp (ICond e1 e2 e3 t) = [cexp| $(compileExp e1) ? $(compileExp e2) : $(compileExp e3) |]
+
+compileExp (IBinOp op e1 e2 t) = go op 
+  where
+    x = compileExp e1
+    y = compileExp e2
+    go IAdd = [cexp| $x + $y |]
+    go ISub = [cexp| $x - $y |]
+    go IMul = [cexp| $x * $y |]
+    go IDiv = [cexp| $x / $y |]
+    go IMod = [cexp| $x % $y |]
+    go IEq = [cexp| $x == $y |]
+    go INotEq = [cexp| $x !=  $y |]
+    go ILt = [cexp| $x < $y |]
+    go IGt = [cexp| $x > $y |]
+    go IGEq = [cexp| $x >= $y |]
+    go ILEq = [cexp| $x <=  $y |]
+    go IAnd = [cexp| $x && $y |]
+    go IOr = [cexp| $x  || $y |]
+    go IPow = case t of
+                Float ->  [cexp|powf($x,$y) |]
+                Double -> [cexp|pow($x,$y)  |] 
+    go IBitwiseAnd = [cexp| $x & $y |]
+    go IBitwiseOr = [cexp| $x | $y |]
+    go IBitwiseXor = [cexp| $x ^ $y |]
+    go IShiftL = [cexp| $x << $y |]
+    go IShiftR = [cexp| $x >> $y |]
+compileExp (IUnOp op e t) = go op
+  where
+    x = compileExp e
+    go IBitwiseNeg = [cexp| ~$x|]
+    go INot        = [cexp| !$x|]
+    
+compileExp (IFunCall name es t) = [cexp| $fc |]
+  where
+    es' = map compileExp es
+    fc  = [cexp| $id:(name)($args:(es')) |]
+
+compileExp (ICast e t) = [cexp| ($ty:(compileType t)) $e' |]
+  where
+    e' = compileExp e
+   
+compileType (Int8) = [cty| typename int8_t |]
+compileType (Int16) = [cty| typename int16_t |]
+compileType (Int32) = [cty| typename int32_t |]
+compileType (Int64) = [cty| typename int64_t |]
+compileType (Word8) = [cty| typename uint8_t |]
+compileType (Word16) = [cty| typename uint16_t |]
+compileType (Word32) = [cty| typename uint32_t |]
+compileType (Word64) = [cty| typename uint64_t |]
+compileType (Float) = [cty| float |]
+compileType (Double) = [cty| double |]
+compileType (Pointer t) = [cty| $ty:(compileType t)* |]
+compileType (Volatile t) =  [cty| volatile $ty:(compileType t)|] 
+
+--compileType' (Volatile t) = [cty| volatile $ty:(compileType t)|] 
+--compileType' t = compileType t 
+
+{-
+
+  Solve the volatile issue.
+  When operating in a warp without syncs Add Volatile to pointers
+  and only there!
+  Need a way to convey this information all the way from the force functions
+  down to the code generator.
+  
+  Also, Change code generation to declare names of pointers into
+  shared memory at different types (that are used in the program)
+    __shared__ uint8_t sbase[X]
+    uint32_t *sbaseU32 = sbase;
+    float *sbaseF = sbase;
+  This changes how offsets are computed within code. May need names for each
+  intermediate array used in the program
+    __shared__ uint8_t sbase[X]
+    uint32_t *arr1 = sbase;
+    uint32_t *arr2 = (sbase + 16U); 
+    float *arr3 = sbase;
+    float *arr4 = (sbase + 4U); 
+   
+
+-} 
+
+---------------------------------------------------------------------------
+-- **     
+-- Compile IM
+-- ** 
+--------------------------------------------------------------------------- 
+
+-- newtype CInfo a = CInfo (State CInfoState a)
+--                 deriving (Monad, MonadState CInfoState) 
+
+-- data CInfoState = CInfoState { cInfoTid  :: (Bool,Exp),
+--                                cInfoWarpID     :: (Bool,Exp),
+--                                cInfoWarpIx     :: (Bool,Exp) } 
+-- evalCInfo (CInfo s) = evalState s                     
+
+---------------------------------------------------------------------------
+-- Statement t to Stm
+---------------------------------------------------------------------------
+
+
+compileStm :: Platform -> Config -> Statement t -> [Stm]
+compileStm p c (SAssign name [] e) =
+   [[cstm| $(compileExp name) = $(compileExp e);|]]
+compileStm p c (SAssign name [ix] e) = 
+   [[cstm| $(compileExp name)[$(compileExp ix)] = $(compileExp e); |]]
+compileStm p c (SAtomicOp name ix atop) = 
+  case atop of
+    AtInc -> [[cstm| atomicInc(&$(compileExp name)[$(compileExp ix)],0xFFFFFFFF); |]]
+    AtAdd e -> [[cstm| atomicAdd(&$(compileExp name)[$(compileExp ix)],$(compileExp e));|]]
+    AtSub e -> [[cstm| atomicSub(&$(compileExp name)[$(compileExp ix)],$(compileExp e));|]]
+    AtExch e -> [[cstm| atomicExch(&$(compileExp name)[$(compileExp ix)],$(compileExp e));|]]
+
+compileStm p c (SCond be im) = [[cstm| if ($(compileExp be)) { $stms:body } |]]
+  where 
+    body = compileIM p c im  -- (compileIM p c im)
+compileStm p c (SSeqFor loopVar n im) = 
+    [[cstm| for (int $id:loopVar = 0; $id:loopVar < $(compileExp n); ++$id:loopVar) 
+              { $stms:body } |]]
+  where
+    body = compileIM p c im -- (compileIM p c im)
+
+
+-- Just relay to specific compileFunction
+compileStm p c a@(SForAll lvl n im) = compileForAll p c a
+
+compileStm p c a@(SDistrPar lvl n im) = compileDistr p c a 
+
+compileStm p c (SSeqWhile b im) =
+  [[cstm| while ($(compileExp b)) { $stms:body}|]]
+  where
+    body = compileIM p c im 
+
+compileStm p c SSynchronize 
+  = case p of
+      PlatformCUDA -> [[cstm| __syncthreads(); |]]
+      PlatformOpenCL -> [[cstm| barrier(CLK_LOCAL_MEM_FENCE); |]]
+
+compileStm _ _ (SAllocate _ _ _) = []
+compileStm _ _ (SDeclare name t) = []
+
+compileStm _ _ a = error  $ "compileStm: missing case "
+
+---------------------------------------------------------------------------
+-- DistrPar 
+---------------------------------------------------------------------------
+compileDistr :: Platform -> Config -> Statement t -> [Stm] 
+compileDistr PlatformCUDA c (SDistrPar Block n im) =  codeQ ++ codeR
+  -- New here is BLOCK virtualisation
+  where
+    cim = compileIM PlatformCUDA c im
+    
+    numBlocks = [cexp| $id:("gridDim.x") |]
+    
+    blocksQ = [cexp| $exp:(compileExp n) / $exp:numBlocks|]
+    blocksR = [cexp| $exp:(compileExp n) % $exp:numBlocks|] 
+    
+    codeQ = [[cstm| for (int b = 0; b < $exp:blocksQ; ++b) { $stms:bodyQ }|]]
+                
+    bodyQ = [cstm| $id:("bid") = blockIdx.x * $exp:blocksQ + b;|] : cim  ++  
+            [[cstm| bid = blockIdx.x;|],
+             [cstm| __syncthreads();|]] -- yes no ? 
+         
+    codeR = [[cstm| bid = ($exp:numBlocks * $exp:blocksQ) + blockIdx.x;|], 
+             [cstm| if (blockIdx.x < $exp:blocksR) { $stms:cim }|],
+             [cstm| bid = blockIdx.x;|], 
+             [cstm| __syncthreads();|]] -- yes no ? 
+                    
+-- Can I be absolutely sure that 'n' here is statically known ? 
+-- I must look over the functions that can potentially create this IM. 
+-- Can make a separate case for unknown 'n' but generate worse code.
+-- (That is true for all levels)  
+compileDistr PlatformCUDA c (SDistrPar Warp (IWord32 n) im) = codeQ  ++ codeR 
+  -- Here the 'im' should be distributed over 'n'warps.
+  -- 'im' uses a warpID variable to identify what warp it is.
+  -- 'n' may be higher than the actual number of warps we have!
+  -- So GPU warp virtualisation is needed. 
+  where
+    cim = compileIM PlatformCUDA c im
+
+    nWarps   = fromIntegral $ configThreadsPerBlock c `div` 32
+    numWarps = [cexp| $int:nWarps|] 
+
+    warpsQ   = [cexp| $int:(n `div` nWarps)|]
+    warpsR   = [cexp| $int:(n `mod` nWarps)|]
+    
+    codeQ = [[cstm| for (int w = 0; w < $exp:warpsQ; ++w) { $stms:bodyQ } |]]
+    
+    bodyQ = [cstm| warpID = (threadIdx.x / 32) * $exp:warpsQ + w;|] : cim ++
+            --[cstm| warpID = w * $exp:warpsQ + (threadIdx.x / 32);|] : cim ++ 
+            [[cstm| warpID = threadIdx.x / 32;|]] 
+
+    codeR = case (n `mod` nWarps)  of 
+             0 -> [] 
+             n -> [[cstm| warpID = ($exp:numWarps * $exp:warpsQ)+ (threadIdx.x / 32);|],
+                   [cstm| if (threadIdx.x / 32 < $exp:warpsR) { $stms:cim } |], 
+                   [cstm| warpID = threadIdx.x / 32; |], 
+                   [cstm| __syncthreads();|]]
+
+---------------------------------------------------------------------------
+-- ForAll is compiled differently for different platforms
+---------------------------------------------------------------------------
+compileForAll :: Platform -> Config -> Statement t -> [Stm]
+compileForAll PlatformCUDA c (SForAll Warp  (IWord32 n) im) = codeQ ++ codeR
+  where
+    nt = 32
+
+    q = n `div` nt
+    r = n `mod` nt
+
+    cim = compileIM PlatformCUDA c im 
+    
+    codeQ =
+      case q of
+        0 -> []
+        1 -> cim
+        n -> [[cstm| for ( int vw = 0; vw < $int:q; ++vw) { $stms:body } |], 
+              [cstm| $id:("warpIx") = threadIdx.x % 32; |]]
+              -- [cstm| __syncthreads();|]]
+             where 
+               body = [cstm|$id:("warpIx") = vw*$int:nt + (threadIdx.x % 32); |] : cim
+               --body = [cstm|$id:("warpIx") = (threadIdx.x % 32) * q + vw; |] : cim
+
+    codeR = 
+      case r of 
+        0 -> [] 
+        n -> [[cstm| if ((threadIdx.x % 32) < $int:r) { 
+                            $id:("warpIx") = $int:(q*32) + (threadIdx.x % 32);  
+                            $stms:cim } |],
+                  -- [cstm| __syncthreads();|],
+                  [cstm| $id:("warpIx") = threadIdx.x % 32; |]]
+
+compileForAll PlatformCUDA c (SForAll Block (IWord32 n) im) = goQ ++ goR 
+  where
+    cim = compileIM PlatformCUDA c im
+   
+    nt = configThreadsPerBlock c 
+
+    q  = n `quot` nt
+    r  = n `rem`  nt 
+
+    -- q is the number full "passes" needed to cover the iteration
+    -- space given we have nt threads. 
+    goQ =
+      case q of
+        0 -> []
+        1 -> cim -- [cstm|$id:loopVar = threadIdx.x; |]:cim
+            --do
+            --  stm <- updateTid [cexp| threadIdx.x |]
+            --  return $ [cstm| $id:loopVar = threadIdx.x; |] : cim 
+        n -> [[cstm| for ( int i = 0; i < $int:q; ++i) { $stms:body } |], 
+              [cstm| $id:("tid") = threadIdx.x; |]]
+             where 
+               body = [cstm|$id:("tid") =  i*$int:nt + threadIdx.x; |] : cim
+   
+    -- r is the number of elements left. 
+    -- This generates code for when fewer threads are 
+    -- needed than available. (some threads shut down due to the conditional). 
+    goR = 
+      case (r,q) of 
+        (0,_) -> [] 
+        --(n,0) -> [[cstm| if (threadIdx.x < $int:n) { 
+        --                    $stms:cim } |]] 
+        (n,m) -> [[cstm| if (threadIdx.x < $int:n) { 
+                            $id:("tid") = $int:(q*nt) + threadIdx.x;  
+                            $stms:cim } |], 
+                  [cstm| $id:("tid") = threadIdx.x; |]]
+
+compileForAll PlatformCUDA c (SForAll Grid n im) = cim
+  -- The grid case is special. May need more thought
+  --error "compileForAll: Grid"
+  where
+    cim = compileIM PlatformCUDA c im
+
+compileForAll PlatformC c (SForAll lvl (IWord32 n) im) = go
+  where
+    body = compileIM PlatformC c im 
+    go  = [ [cstm| for (int i = 0; i <$int:n; ++i) { $stms:body } |] ] 
+      
+
+---------------------------------------------------------------------------
+-- compileWarp (needs fixing so that warpIx and warpID are always correct) 
+---------------------------------------------------------------------------
+-- compileWarp :: Platform -> Config -> IExp -> IMList t -> [Stm]
+-- compileWarp PlatformCUDA c (IWord32 warps) im = 
+--     concatMap (go . fst)  im
+--   where
+--     go (SAllocate nom n t) = []  -- Skip allocations at this point. Been handled earlier. 
+--     go (SForAll Warp (IWord32 n) im) =
+--         {-
+--            warps : number of warps that should execute this code
+--            n     : number of threads per each of these warps
+
+--            warps may be a larger number than there are actual warps
+--            n map be a larger number than 32 (the actual number of threads per warp)
+
+--            c : is the configuration. It knows the actual number of threads
+--                we are generating code for. 
+--         -} 
+--         if (wholeRealWarps <= 0)
+--         then error "compileWarp: Atleast one full warp of real threads needed!"
+--         else 
+--           case (wholeRealWarps `compare` warps) of
+--             {-
+--                3 potential outcomes.
+--                #1 we have more real warps than we need
+--                #2 we have exactly the amount of warps that we need
+--                #3 we have fewer real warps than requested.
+--             -} 
+--             GT -> [[cstm| if (threadIdx.x < $int:(warps*32)) {
+--                            $stms:(goQ ++ goR) } |]]   
+--             EQ -> goQR
+--             LT -> wQ ++ wR 
+             
+--        where
+--         cim = compileIM PlatformCUDA c im
+ 
+--         nt = configThreadsPerBlock c
+--         wholeRealWarps = nt `quot` 32
+--         threadsPartialWarp = nt `rem` 32 -- Do not use partial warps! 
+       
+--         -- Maybe something else than s for goR
+--         -- It will depend on how 
+--         goQR = goQ  ++ goR 
+  
+--         threadQ = n `quot` 32   -- Set up virtual threads within warp
+--         threadR = n `rem`  32   -- 
+
+
+-- ###################################################################### 
+-- CLEAN THIS MESS UP ! !!! ! ! ! !! ! !!!! ! !! ! ! ! ! !  ! !!! ! ! ! !
+-- ######################################################################
+
+
+        -- -- Compile for warps, potentially with virtual threads  
+        -- goQ :: [Stm]
+        -- goQ = case threadQ of 
+        --         0 -> [] 
+        --         1 -> cim 
+        --         n -> [[cstm| for (int i = 0; i < $int:threadQ; ++i) {
+        --                                 $stms:body } |],
+        --                          [cstm| $id:("warpIx") = threadIdx.x % 32; |]]
+        --             where 
+        --               body = [cstm| $id:("warpIx") = i*32 + threadIdx.x % 32; |] : cim
+        -- goR :: [Stm] 
+        -- goR = case threadR of 
+        --         0 -> [] 
+        --         _ -> [[cstm| if ( threadIdx.x % 32 < $int:threadR) { 
+        --                        $id:("warpIx") = $int:(threadQ*32) + (threadIdx.x % 32);
+        --                        $stms:cim } |],
+        --                 [cstm| $id:("warpIx") = threadIdx.x % 32;|]]
+
+        -- -- Compile for virtual warps 
+        -- warpQ = warps `quot` wholeRealWarps -- wholeRealWarps may be zero! 
+        -- warpR = warps `rem`  wholeRealWarps
+        -- -- each warp will pretend to be warpQ number of warps.        
+         
+        -- -- warp jump size 
+        -- -- warpQ is the number of virtual warps that EACH real warp 
+        -- -- will pretend to be. 
+        
+        -- -- warpR is a number of leftover warps (less than wholeRealWarps) 
+        -- -- that need to be processed afterwards. 
+        
+        -- wQ = case warpQ of 
+        --        0 -> [] 
+        --        1 -> goQR  
+        --        n -> [[cstm| for (int vw = 0; vw < $int:warpQ; ++vw) { $stms:body } |],
+        --              [cstm| $id:("warpID") = threadIdx.x / 32; |]]
+        --            where 
+        --              --- vw * wholeRealWarps is incorrect. 
+        --              --- 
+        --              body = [cstm| $id:("warpID") = vw + ((threadIdx.x / 32)*$int:warpQ); |] : goQR
+
+        -- wR = case warpR of 
+        --        0 -> [] 
+        --        -- Skip one division, 
+        --        -- Premature optimisation? CUDA compiler will share the threadIdx.x / 32 result. 
+        --        n -> [[cstm| if (threadIdx.x < $int:(n * 32)) {
+        --                       $id:("warpID") = $int:(warpQ*wholeRealWarps) + threadIdx.x / 32; 
+        --                       $stms:(goQR) } |], -- cim
+        --                     [cstm| $id:("warpID") = threadIdx.x / 32;|]]
+--------------------------------------------------------------------------- 
+-- CompileIM to list of Stm 
+--------------------------------------------------------------------------- 
+compileIM :: Platform -> Config -> IMList a -> [Stm]
+compileIM pform conf im = concatMap ((compileStm pform conf) . fst) im
+
+
+---------------------------------------------------------------------------
+-- Generate entire Kernel 
+---------------------------------------------------------------------------
+type Parameters = [(String,Obsidian.Types.Type)]
+
+compile :: Platform -> Config -> String -> (Parameters,IMList a) -> Definition
+compile pform config kname (params,im)
+  = go pform 
+  where
+    stms = compileIM pform config im
+    
+    ps = compileParams pform params
+    go PlatformCUDA
+      = [cedecl| extern "C" __global__ void $id:kname($params:ps) {$items:cudabody} |]
+    go PlatformOpenCL
+      = [CL.cedecl| __kernel void $id:kname($params:ps) {$stms:stms} |]
+    go PlatformC
+      = [cedecl| extern "C" void $id:kname($params:ps) {$items:cbody} |] 
+
+    cudabody = (if (configSharedMem config > 0)
+                -- then [BlockDecl [cdecl| extern volatile __shared__  typename uint8_t sbase[]; |]] 
+                then [BlockDecl [cdecl| __shared__  typename uint8_t sbase[$uint:(configSharedMem config)]; |]] 
+                else []) ++
+                --[BlockDecl [cdecl| typename uint32_t tid = threadIdx.x; |]] ++
+                --[BlockDecl [cdecl| typename uint32_t warpID = threadIdx.x / 32; |],
+                --       BlockDecl [cdecl| typename uint32_t warpIx = threadIdx.x % 32; |]] ++
+--                [BlockDecl [cdecl| typename uint32_t bid = blockIdx.x; |]] ++
+               (if (usesGid im) 
+                then [BlockDecl [cdecl| typename uint32_t gid = blockIdx.x * blockDim.x + threadIdx.x; |]]
+                else []) ++ 
+               (if (usesBid im) 
+                then [BlockDecl [cdecl| typename uint32_t bid = blockIdx.x; |]]
+                else []) ++ 
+               (if (usesTid im) 
+                then [BlockDecl [cdecl| typename uint32_t tid = threadIdx.x; |]]
+                else []) ++
+               (if (usesWarps im) 
+                then  [BlockDecl [cdecl| typename uint32_t warpID = threadIdx.x / 32; |],
+                       BlockDecl [cdecl| typename uint32_t warpIx = threadIdx.x % 32; |]] 
+                else []) ++
+                -- All variables used will be unique and can be declared 
+                -- at the top level 
+                concatMap declares im ++ 
+                -- Not sure if I am using language.C correctly. 
+                -- Maybe compileSTM should create BlockStms ?
+                -- TODO: look how Nikola does it. 
+                map BlockStm stms
+
+    cbody = -- add memory allocation 
+            map BlockStm stms
+
+-- Declare variables. 
+declares (SDeclare name t,_) = [BlockDecl [cdecl| $ty:(compileType t)  $id:name;|]]
+declares (SCond _ im,_) = concatMap declares im 
+declares (SSeqWhile _ im,_) = concatMap declares im
+declares (SForAll _ _ im,_) = concatMap declares im
+declares (SDistrPar _ _ im,_) = concatMap declares im 
+-- declares (SForAllBlocks _ im,_) = concatMap declares im
+-- declares (SNWarps _ im,_) = concatMap declares im
+-- declares (SWarpForAll _ im,_) = concatMap declares im 
+declares _ = []
+
+
+--------------------------------------------------------------------------- 
+-- Parameter lists for functions  (kernel head) 
+---------------------------------------------------------------------------
+compileParams :: Platform -> Parameters -> [Param]
+compileParams PlatformOpenCL = map go
+  where
+    go (name,Pointer t) = [CL.cparam| global  $ty:(compileType t) $id:name |]
+    go (name, t)        = [CL.cparam| $ty:(compileType t) $id:name |]
+
+-- C or CUDA 
+compileParams _ = map go
+  where
+    go (name,t) = [cparam| $ty:(compileType t) $id:name |]
+ 
diff --git a/Obsidian/CodeGen/InOut.hs b/Obsidian/CodeGen/InOut.hs
deleted file mode 100644
--- a/Obsidian/CodeGen/InOut.hs
+++ /dev/null
@@ -1,155 +0,0 @@
-{-# LANGUAGE FlexibleInstances,
-             OverlappingInstances,
-             UndecidableInstances,
-             FlexibleContexts,
-             MultiParamTypeClasses,
-             TypeOperators,
-             TypeFamilies ,
-             ScopedTypeVariables
-             #-}
-
-{- Joel Svensson 2012, 2013
-   Niklas Ulvinge 2013
-
-  Notes:
-
-  2013-01-24: Changes with the new Array types in mind
-  2013-01-08: Edited
-  2012-12-10: Edited
-
--} 
-
-
-module Obsidian.CodeGen.InOut where 
-
-import Obsidian.Exp 
-import Obsidian.Array
-
-import Obsidian.Types
-import Obsidian.Globs 
-import Obsidian.Program
-import Obsidian.Force
-import Obsidian.Memory
-
-import qualified Obsidian.CodeGen.Program as CG 
-
-import Data.Word
-import Data.Int
-      
----------------------------------------------------------------------------
--- New approach (hopefully)
----------------------------------------------------------------------------
--- "reify" Haskell functions into CG.Programs
-
-{-
-   Blocks needs to be of specific sizes (a design choice we've made).
-   Because of this a prototypical input array needs to be provided
-   that has a static block size (the number of blocks is dynamic).
-
-   To make things somewhat general a heterogeneous list of input arrays
-   that has same shape as the actual parameter list of the function
-   is passed into toProgram (the reifyer). 
-
--} 
-  
-type Inputs = [(Name,Type)] 
-
-class ToProgram a b where
-  toProgram :: Int -> (a -> b) -> Ips a b -> (Inputs,CG.IM)
-  
-typeOf_ a = typeOf (Literal a)
-
-
----------------------------------------------------------------------------
--- Base cases
---------------------------------------------------------------------------- 
-instance (Scalar t) => ToProgram (Exp t) (GProgram b) where
-  toProgram i f a = ([(nom,t)],CG.compileStep1 (f input))
-    where nom = "s" ++ show i
-          input = variable nom
-          t = typeOf_ (undefined :: t)
-
-instance (Scalar t) => ToProgram (Pull (Exp Word32) (Exp t)) (GProgram a) where
-  toProgram i f (Pull n ixf) = ([(nom,Pointer t),(n,Word32)],CG.compileStep1 (f input)) 
-      where nom = "input" ++ show i
-            n   = "n" ++ show i 
-            lengthVar = variable n
-            input = namedGlobal nom lengthVar
-            t = typeOf_ (undefined :: t)
-
-instance (Scalar t) => ToProgram (Pull Word32 (Exp t)) (GProgram a) where
-  toProgram i f (Pull n ixf) = ([(nom,Pointer t){-,(n,Word32)-}],CG.compileStep1 (f input)) 
-      where nom = "input" ++ show i
-            --n   = "n" ++ show i 
-            --lengthVar = variable n
-            input = namedGlobal nom n -- lengthVar
-            t = typeOf_ (undefined :: t)
-
----------------------------------------------------------------------------
--- More natural to work with these in some cases
----------------------------------------------------------------------------
-instance (ToProgram b (GProgram ()),
-          GlobalMemoryOps a)
-          => ToProgram b (Push Grid Word32 a) where
-  toProgram i f arr = toProgram i (forceG . f)  arr
-
-instance (ToProgram b (GProgram ()),
-          GlobalMemoryOps a)
-          => ToProgram b (Push Grid EWord32 a) where
-  toProgram i f arr = toProgram i (forceG . f)  arr
-
-
----------------------------------------------------------------------------
--- Recursive cases
---------------------------------------------------------------------------- 
-instance (Scalar t, ToProgram b c) => ToProgram (Exp t) (b -> c) where
-  toProgram i f (a :-> rest) = ((nom,t):ins,prg)
-    where
-      (ins,prg) = toProgram (i+1) (f input) rest
-      nom = "s" ++ show i
-      input = variable nom
-      t = typeOf_ (undefined :: t)
-
-instance (Scalar t, ToProgram b c) => ToProgram (Pull (Exp Word32) (Exp t)) (b -> c) where
-  toProgram i f ((Pull n ixf) :-> rest) = ((nom,Pointer t):(n,Word32):ins,prg)
-    where
-      (ins,prg) = toProgram (i+1) (f input) rest
-      nom = "input" ++ show i
-      n   = "n" ++ show i
-      lengthVar = variable n
-      input = namedGlobal nom lengthVar
-      t = typeOf_ (undefined :: t)
-
-
-instance (Scalar t, ToProgram b c) => ToProgram (Pull Word32 (Exp t)) (b -> c) where
-  toProgram i f ((Pull n ixf) :-> rest) = ((nom,Pointer t){-:(n,Word32)-}:ins,prg)
-    where
-      (ins,prg) = toProgram (i+1) (f input) rest
-      nom = "input" ++ show i
-      --n   = "n" ++ show i
-      --lengthVar = variable n
-      input = namedGlobal nom n --lengthVar
-      t = typeOf_ (undefined :: t)
-
-
----------------------------------------------------------------------------
--- heterogeneous lists of inputs 
----------------------------------------------------------------------------
-data head :-> tail = head :-> tail
-
-infixr 5 :->
-
-
----------------------------------------------------------------------------
--- Function types to input list types. 
---------------------------------------------------------------------------- 
-type family Ips a b
- 
--- type instance Ips a (GlobArray b) = Ips' a -- added Now 26
--- type instance Ips a (Final (GProgram b)) = a
-type instance Ips a (Push Grid l b) = a
-type instance Ips a (Pull l b) = a
-type instance Ips a (GProgram b) = a
-type instance Ips a (b -> c) =  a :-> Ips b c
-
-
diff --git a/Obsidian/CodeGen/Liveness.hs b/Obsidian/CodeGen/Liveness.hs
--- a/Obsidian/CodeGen/Liveness.hs
+++ b/Obsidian/CodeGen/Liveness.hs
@@ -62,16 +62,30 @@
     process (SAssign nom ixs e,_) =
       do
         s <- get
-        let arrays = collectArrays e
-            living = Set.fromList (nom:arrays) `Set.union` s
+        let arrays = collectArraysI "arr" e
+            arrays1 = collectArraysI "arr" nom
+            living = Set.fromList (arrays1++arrays) `Set.union` s
         
         put living  -- update state   
         return (SAssign nom ixs e,living)
-    
-    process (SAtomicOp n1 n2 ixs op,_) =
+    process (SAtomicOp nom ix atop,_) =
       do
         s <- get
-        return (SAtomicOp n1 n2 ixs op,s)
+        let arrays =
+              case atop of
+                AtInc -> []
+                AtAdd e -> collectArraysI "arr" e
+                AtSub e -> collectArraysI "arr" e
+                AtExch e -> collectArraysI "arr" e
+            arrays1 = collectArraysI "arr" nom
+            living = Set.fromList (arrays1++arrays) `Set.union` s
+        put living
+        return (SAtomicOp nom ix atop, living) 
+    
+--    process (SAtomicOp n1 n2 ixs op,_) =
+--      do
+--        s <- get
+--        return (SAtomicOp n1 n2 ixs op,s)
         
     process (SAllocate name size t,_) =
       do
@@ -84,11 +98,6 @@
         s <- get 
         return (SDeclare name t,s)
 
-    process (SOutput name t,_) = 
-      do 
-        s <- get 
-        return (SOutput name t,s)
-
     process (SSynchronize,_) = 
       do 
         s <- get
@@ -113,32 +122,68 @@
             ns  = s `Set.union` l
         put ns
         return (SSeqFor nom n iml,ns) 
-
-
-    process (SForAll n im,_) =  
+    process (SSeqWhile b im,_) = 
       do 
         s <- get 
         let iml = computeLiveness1 s im 
             l   = safeHead iml 
             ns  = s `Set.union` l
-        put ns
-        return (SForAll n iml,ns) 
-    
-    process (SForAllBlocks n im,_) = 
+        put ns 
+        return (SSeqWhile b iml,ns)
+    process (SBreak,_) = 
       do 
         s <- get 
+        return (SBreak,s)
+
+
+    process (SForAll lvl n im,_) =  
+      do 
+        s <- get 
         let iml = computeLiveness1 s im 
             l   = safeHead iml 
             ns  = s `Set.union` l
         put ns
-        return (SForAllBlocks n iml,ns)
-
-    process (SForAllThreads n im,_) = 
+        return (SForAll lvl n iml,ns) 
+    process (SDistrPar lvl n im,_) = 
       do 
         s <- get 
         let iml = computeLiveness1 s im 
             l   = safeHead iml 
-            ns  = s `Set.union` l
+            ns  = s `Set.union` l 
         put ns 
-        return (SForAllThreads n iml,ns)
+        return (SDistrPar lvl n iml,ns) 
+    -- process (SForAllBlocks n im,_) = 
+    --   do 
+    --     s <- get 
+    --     let iml = computeLiveness1 s im 
+    --         l   = safeHead iml 
+    --         ns  = s `Set.union` l
+    --     put ns
+    --     return (SForAllBlocks n iml,ns)
+    -- process (SNWarps n im,_) = 
+    --   do 
+    --     s <- get
+    --     let iml = computeLiveness1 s im 
+    --         l   = safeHead iml 
+    --         ns  = s `Set.union` l
+    --     put ns
+    --     return (SNWarps n iml,ns)
+    -- process (SWarpForAll n im,_) =  
+    --   do 
+    --     s <- get 
+    --     let iml = computeLiveness1 s im 
+    --         l   = safeHead iml 
+    --         ns  = s `Set.union` l
+    --     put ns
+    --     return (SWarpForAll n iml,ns) 
+
+
+    -- process (SForAllThreads n im,_) = 
+    --   do 
+    --     s <- get 
+    --     let iml = computeLiveness1 s im 
+    --         l   = safeHead iml 
+    --         ns  = s `Set.union` l
+    --     put ns 
+    --     return (SForAllThreads n iml,ns)
 
diff --git a/Obsidian/CodeGen/Memory.hs b/Obsidian/CodeGen/Memory.hs
--- a/Obsidian/CodeGen/Memory.hs
+++ b/Obsidian/CodeGen/Memory.hs
@@ -15,9 +15,8 @@
         sharedMem,  
         Address,
         Bytes,
-        --mapMemory,
-        -- NEW
-        mmIM) 
+        mmIM,
+        renameIM ) 
        where 
 
 import qualified Data.List as List
@@ -33,8 +32,9 @@
 
 import qualified Data.Map as Map 
 
-------------------------------------------------------------------------------
+---------------------------------------------------------------------------
 -- Memory layout
+---------------------------------------------------------------------------
 
 type MemMap = Map.Map Name (Word32,Type)
 
@@ -73,6 +73,7 @@
                                 else fl
                       in  (updateMax (m {freeList = fl', 
                                          allocated = (a,b):allocated m}) ,a)
+    [] -> error "out of shared memory"
    
                           
 free :: Memory -> Address -> Memory
@@ -127,26 +128,122 @@
               (Just as) -> freeAll m' (map fst as)
               Nothing   -> m'
       in r xs (mNew,mm')
-         
-    process (SAllocate name size t,_) m mm = (m',mm') 
-      where (m',addr) = allocate m size
-            mm' = case Map.lookup name mm of
-                      Nothing -> Map.insert name (addr,t) mm
-                      (Just (a, t)) -> error $ "mmIm: " ++ name ++ " is already mapped to " ++ show a
+mmIM' :: IML -> Memory -> MemMap -> (Memory, MemMap)
+mmIM' im memory memmap = r im (memory,memmap)
+  where
+    r [] m = m
+    r (x:xs) (m,mm) =
+      let
+          (m',mm') = process x m mm
+           
+          freeable = getFreeableSet x xs
+          freeableAddrs = mapM (flip Map.lookup mm') (filter dontMap (Set.toList freeable))
+          dontMap name = not ((List.isPrefixOf "input" name) || 
+                              (List.isPrefixOf "output" name))
+          mNew =
+            case freeableAddrs of
+              (Just as) -> m' -- freeAll m' (map fst as)
+              Nothing   -> m'
+      in r xs (mNew,mm')
 
+    
+process (SAllocate name size t,_) m mm = (m',mm') 
+  where (m',addr) = allocate m size
+        mm' = case Map.lookup name mm of
+          Nothing -> Map.insert name (addr,t) mm
+          (Just (a, t)) -> error $ "mmIm: " ++ name ++ " is already mapped to " ++ show a
+
     -- A tricky case.                      
-    process (SForAllBlocks n im,_) m mm = mmIM im m mm
+--    process (SForAllBlocks n im,_) m mm = mmIM im m mm
     -- Another tricky case. 
-    process (SSeqFor _ n im,_) m mm = mmIM im m mm
+process (SSeqFor _ n im,_) m mm = mmIM im m mm
+process (SSeqWhile b im,_) m mm = mmIM im m mm 
     -- Yet another tricky case.
-    process (SForAll n im,_) m mm = mmIM im m mm 
+process (SForAll _ n im,_) m mm = mmIM im m mm
+process (SDistrPar Warp n im,_) m mm = mmIM' im m mm -- mmIM im m mm
+process (SDistrPar Block n im,_) m mm = mmIM im m mm 
     -- The worst of them all.
-    process (SForAllThreads n im,_) m mm = mmIM im m mm
+--    process (SForAllThreads n im,_) m mm = mmIM im m mm
+--    process (SNWarps _ im,_) m mm = mmIM im m mm
+--    process (SWarpForAll _ im,_) m mm = mmIM im m mm 
 
-    process (_,_) m mm = (m,mm) 
+  --  process im m mm = error $ printStm im -- "process: WHat!"
+process (_,_) m mm = (m,mm) 
 
 -- Friday (2013 Mars 29, discovered bug) 
 getFreeableSet :: (Statement Liveness,Liveness) -> IML -> Liveness 
 getFreeableSet (_,l) [] = Set.empty -- not l ! 
 getFreeableSet (_,l) ((_,l1):_) = l Set.\\ l1
+
+---------------------------------------------------------------------------
+-- Rename arrays in IM
+--------------------------------------------------------------------------- 
+
+renameIM :: MemMap -> IML -> IMList ()
+renameIM mm im = zip (map (go . fst) im) (repeat ())
+  where
+    go (SAssign name ix e) = SAssign (renameIVar mm name)
+                                     (map (renameIExp mm) ix)
+                                     (renameIExp mm e)
+    go (SAtomicOp name ix atop) = SAtomicOp (renameIVar mm name)
+                                            (renameIExp mm ix)
+                                            (renameAtOp mm atop) 
+    go (SCond be im) = SCond (renameIExp mm be)
+                             (renameIM mm im)
+    go (SSeqFor str n im) = SSeqFor str (renameIExp mm n)
+                                        (renameIM mm im)
+    go SBreak = SBreak
+    go (SSeqWhile n im) = SSeqWhile (renameIExp mm n)
+                                    (renameIM mm im)
+    go (SForAll lvl n im)   = SForAll lvl (renameIExp mm n)
+                                          (renameIM mm im)
+    go (SDistrPar lvl n im) = SDistrPar lvl (renameIExp mm n)
+                                            (renameIM mm im) 
+
+--    go (SForAllBlocks n im) = SForAllBlocks (renameIExp mm n)
+--                                            (renameIM mm im)
+--    go (SNWarps n im) = SNWarps (renameIExp mm n)
+--                                (renameIM mm im)
+--    go (SWarpForAll n im) = SWarpForAll (renameIExp mm n)
+--                                        (renameIM mm im) 
+    -- Strip this out earlier. 
+    go (SAllocate name n t)  = SAllocate name n t 
+    go (SDeclare name t) = SDeclare name t
+    go SSynchronize      = SSynchronize 
+
+---------------------------------------------------------------------------
+-- Memory map the arrays in an CExpr
+---------------------------------------------------------------------------
+renameIExp mm e@(IVar nom t) =  renameIVar mm e 
+renameIExp mm (IIndex (e1,es) t) = IIndex (renameIExp mm e1, map (renameIExp mm) es) t
+renameIExp mm (IBinOp op e1 e2 t) = IBinOp op (renameIExp mm e1) (renameIExp mm e2) t
+renameIExp mm (IUnOp op e t) = IUnOp op (renameIExp mm e) t 
+renameIExp mm (IFunCall nom exprs t) = IFunCall nom (map (renameIExp mm) exprs) t
+renameIExp mm (ICast e t) = ICast (renameIExp mm e) t
+renameIExp mm (ICond e1 e2 e3 t) = ICond (renameIExp mm e1)
+                                         (renameIExp mm e2)
+                                         (renameIExp mm e3)
+                                         t
+renameIExp mm a = a
+
+
+renameIVar mm (IVar name t) =
+    case Map.lookup name mm of 
+    Just (addr,t) -> 
+      let core = sbaseIExp addr 
+          cast c = ICast c t
+      in cast core
+    
+    Nothing -> IVar name t
+    where
+      sbaseIExp 0    = IVar "sbase" (Pointer Word8) 
+      sbaseIExp addr = IBinOp IAdd (IVar "sbase" (Pointer Word8)) 
+                                   (IWord32 addr) 
+                                   (Pointer Word8) 
+
+renameAtOp mm AtInc = AtInc
+renameAtOp mm (AtAdd e) = AtAdd (renameIExp mm e)
+renameAtOp mm (AtSub e) = AtSub (renameIExp mm e)
+renameAtOp mm (AtExch e) = AtExch (renameIExp mm e) 
+
 
diff --git a/Obsidian/CodeGen/PP.hs b/Obsidian/CodeGen/PP.hs
deleted file mode 100644
--- a/Obsidian/CodeGen/PP.hs
+++ /dev/null
@@ -1,55 +0,0 @@
-
-{- Joel Svensson 2012 -}
-module Obsidian.CodeGen.PP where 
-
-
-import Control.Monad.State
-------------------------------------------------------------------------------
--- print and indent and stuff... 
---  This is probably very ugly 
-
--- TODO: There is a chapter about this pretty printing in "implementing functional lang..." 
---       Look at that and learn 
-
-
-type PP a = State (Int,String) a  
-
-indent :: PP ()
-indent = 
-  do 
-    (i,s) <- get 
-    put (i+1,s) 
-    
-unindent :: PP () 
-unindent = 
-  do 
-    (i,s) <- get 
-    if i <= 0 then error "Whats going on" else put (i-1,s) 
-
-line :: String -> PP () 
-line str = 
-  do 
-    (i,s) <- get 
-    put (i,s ++ str) 
-
-  
-newline :: PP () 
-newline = 
-  do 
-    (i,s) <- get 
-    let ind = replicate (i*2) ' '
-    put (i,s ++ "\n" ++ ind)
-    
-runPP :: PP a -> Int -> String
-runPP pp i = snd$ execState pp (i,"")
-
-begin :: PP () 
-begin = line "{" >> indent >> newline
-
-end :: PP () 
-end =  unindent >> newline >> line "}" >> newline
-
-space   = line " " 
-cTermLn = line ";" >> newline
-
-wrap s e p = line s >> p >> line e 
diff --git a/Obsidian/CodeGen/Program.hs b/Obsidian/CodeGen/Program.hs
--- a/Obsidian/CodeGen/Program.hs
+++ b/Obsidian/CodeGen/Program.hs
@@ -1,6 +1,8 @@
 {-# LANGUAGE GADTs,
              ExistentialQuantification,
-             FlexibleInstances #-}
+             FlexibleInstances  #-}
+{-# LANGUAGE TypeSynonymInstances #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
 
 {- CodeGen.Program.
 
@@ -8,10 +10,6 @@
 
    Notes:
      2013-03-17: Codegeneration is changing
-
-
-
-
 -} 
 
 
@@ -22,7 +20,8 @@
 import Obsidian.Types
 import Obsidian.Atomic
 
-import qualified Obsidian.Program as P 
+import qualified Obsidian.Program as P
+import Obsidian.Program (Step,Zero)
 
 import Data.Word
 import Data.Supply
@@ -30,6 +29,9 @@
 
 import System.IO.Unsafe
 
+import Control.Monad.State
+
+
 ---------------------------------------------------------------------------
 -- New Intermediate representation
 ---------------------------------------------------------------------------
@@ -41,129 +43,273 @@
 -- out :: 
 out a = [(a,())]
 
+-- Atomic operations
+data AtOp = AtInc
+          | AtAdd  IExp 
+          | AtSub  IExp 
+          | AtExch IExp 
 
-data Statement t = forall a. (Show a, Scalar a) => SAssign Name [Exp Word32] (Exp a)
-               | forall a. (Show a, Scalar a) => SAtomicOp Name Name (Exp Word32) (Atomic a)
-               | SCond (Exp Bool) (IMList t) 
-               | SSeqFor String (Exp Word32) (IMList t)
-                 -- See if it is possible to get away
-                 -- with only one kind of ForAll (plus maybe some flag) 
-               | SForAll (Exp Word32) (IMList t) 
-               | SForAllBlocks (Exp Word32) (IMList t)
-                 -- a special loop over all threads..
-               | SForAllThreads (Exp Word32) (IMList t)
 
-                 -- Memory Allocation..
-               | SAllocate Name Word32 Type
-               | SDeclare  Name Type
-               | SOutput   Name Type
+data HLevel = Thread
+            | Warp
+            | Block
+            | Grid
+            
 
-                 -- Synchronisation
-               | SSynchronize
+-- Statements 
+data Statement t = SAssign IExp [IExp] IExp
+                 | SAtomicOp IExp IExp AtOp
+                 | SCond IExp (IMList t) 
+                 | SSeqFor String IExp (IMList t)
+                 | SBreak
+                 | SSeqWhile IExp (IMList t)
 
-                 -- ProgramPar and ProgramSeq does not exist
-                -- at this level (the par or seq info is lost!)
-               
+                 --               Iters   Body
+                 | SForAll HLevel IExp    (IMList t)
+                 | SDistrPar HLevel IExp  (IMList t)
+                   
+--                 | SForAllBlocks IExp (IMList t)
+--                 | SNWarps IExp (IMList t)
+--                 | SWarpForAll  IExp (IMList t)                    
+--                 | SWarpForAll String   String IExp (IMList t) 
 
-compileStep1 :: P.Program t a -> IM
-compileStep1 p = snd $ cs1 ns p
-  where
-    ns = unsafePerformIO$ newEnumSupply
+    -- Memory Allocation..
+                 | SAllocate Name Word32 Type
+                 | SDeclare  Name Type
 
-cs1 :: Supply Int -> P.Program t a -> (a,IM) 
-cs1 i P.Identifier = (supplyValue i, [])
+    -- Synchronisation
+                 | SSynchronize
 
-cs1 i (P.Assign name ix e) =
-  ((),out (SAssign name ix e))
- 
-cs1 i (P.AtomicOp name ix at) = (v,out im)
-  where 
-    nom = "a" ++ show (supplyValue i)
-    v = variable nom
-    im = SAtomicOp nom name ix at
-      
-cs1 i (P.Cond bexp p) = ((),out (SCond bexp im)) 
-  where ((),im) = cs1 i p
 
-cs1 i (P.SeqFor n f) = (a,out (SSeqFor nom n im))
-  where
-    (i1,i2) = split2 i
-    nom = "i" ++ show (supplyValue i1)
-    v = variable nom
-    p = f v
-    (a,im) = cs1 i2 p 
-    
-cs1 i (P.ForAll n f) = (a,out (SForAll n im))
-  where
-    p = f (ThreadIdx X)  
-    (a,im) = cs1 i p 
+--------------------------------------------------------------------------- 
+-- Collect and pass around data during first step compilation
+data Context = Context { ctxNWarps :: Maybe Word32,
+                         ctxGLBUsesTid :: Bool,
+                         ctxGLBUsesWid :: Bool}
 
-cs1 i (P.ForAllBlocks n f) = (a,out (SForAllBlocks n im)) 
-  where
-    p = f (BlockIdx X)
-    (a,im) = cs1 i p
+newtype CM a = CM (State Context a)
+   deriving (Monad, MonadState Context)
 
+runCM :: CM a -> Context -> a 
+--runCM (CM cm) ctx = evalState cm ctx
+runCM (CM cm) = evalState  cm 
 
--- Warning: Every thread will ALWAYS need to perform a conditional
---     (Only in special case is the conditional not needed) 
--- TRY To express all library functions using ForAllBlocks + ForAll
--- For more flexibility and probably in the end performance. 
-cs1 i (P.ForAllThreads n f) = (a,out (SForAllThreads n im)) 
-  where
-    p = f (BlockIdx X * BlockDim X + ThreadIdx X)
-    (a,im) = cs1 i p
+evalCM :: CM a -> Context -> (a, Context)
+evalCM (CM cm) = runState cm 
 
+setUsesTid :: CM ()
+setUsesTid = modify $ \ctx -> ctx { ctxGLBUsesTid = True } 
 
-cs1 i (P.Allocate id n t) = ((),out (SAllocate id n t))
-cs1 i (P.Declare  id t)   = ((),out (SDeclare id t))
--- Output works in a different way! (FIX THIS!)
---  Uniformity! (Allocate Declare Output) 
-cs1 i (P.Output   t)      = (nom,out (SOutput nom t))
-  where nom = "output" ++ show (supplyValue i) 
-cs1 i (P.Sync)            = ((),out (SSynchronize))
+setUsesWid :: CM ()
+setUsesWid = modify $ \ctx -> ctx { ctxGLBUsesWid = True } 
 
+enterWarp :: Word32 -> CM ()
+enterWarp  n = modify $ \ctx -> ctx { ctxNWarps = Just n }
 
-cs1 i (P.Bind p f) = (b,im1 ++ im2) 
+clearWarp :: CM ()
+clearWarp = modify $ \ctx -> ctx {ctxNWarps = Nothing}
+
+getNWarps :: CM (Maybe Word32)
+getNWarps = do
+  ctx <- get
+  return $ ctxNWarps ctx 
+
+emptyCtx = Context Nothing False False
+---------------------------------------------------------------------------
+
+
+-- Sort these out and improve! 
+usesWarps :: IMList t -> Bool
+usesWarps = any (go . fst)
   where
-    (s1,s2) = split2 i
-    (a,im1) = cs1 s1 p
-    (b,im2) = cs1 s2 (f a)
+    go (SDistrPar _ _ im) = usesWarps im 
+    go (SForAll Warp _ _) = True
+    go _ = False
 
-cs1 i (P.Return a) = (a,[])
+usesTid :: IMList t -> Bool
+usesTid = any (go . fst)
+  where
+    go (SDistrPar _ _ im) = usesTid im 
+    go (SForAll Block _ _) = True
+    go _ = False
+usesBid :: IMList t -> Bool
+usesBid = any (go . fst)
+  where
+    go (SDistrPar Block _ im) = True -- usesBid im
+    -- go (SForAll Block _ _) = True
+    go _ = False
+usesGid :: IMList t -> Bool
+usesGid = any (go . fst)
+  where
+    go (SForAll Grid _ _) = True
+    go _ = False
+ 
 
 
 ---------------------------------------------------------------------------
--- Analysis
+-- COmpilation of Program to IM
 --------------------------------------------------------------------------- 
-numThreads :: IMList a -> Either Word32 (EWord32)
-numThreads im = foldl maxCheck (Left 0) $ map process im
+
+compileStep1 :: Compile t => P.Program t a -> IM
+compileStep1 p = snd $ runCM (compile ns p) emptyCtx
   where
-    process (SCond bexp im,_) = numThreads im
-    process (SSeqFor _ _ _,_) = Left 1
-    process (SForAll (Literal n) _,_) = Left n
-    process (SForAll n _,_) = Right n
-    process (SForAllBlocks _ im,_) = numThreads im
-    process (SForAllThreads n im,_) = Right (variable "UNKNOWN") --fix this!
-    process a = Left 0 -- ok ? 
+    ns = unsafePerformIO$ newEnumSupply
 
-    maxCheck (Left a) (Right b)  = Right $ max (fromIntegral a) b
-    maxCheck (Right a) (Left b)  = Right $ max a (fromIntegral b)
-    maxCheck (Left a) (Left  b)  = Left  $ max a b
-    maxCheck (Right a) (Right b) = Right $ max a b
 
+class Compile t where
+  compile :: Supply Int -> P.Program t a -> CM (a,IM)
 
-getOutputs :: IMList a -> [(Name,Type)]
-getOutputs im = concatMap process im
-  where
-    process (SOutput name t,_)      = [(name,t)]
-    process (SSeqFor _ _ im,_)      = getOutputs im
-    process (SForAll _ im,_)        = getOutputs im
-    process (SForAllBlocks _ im,_)  = getOutputs im
-    process (SForAllThreads _ im,_) = getOutputs im
-    process a = []
+-- Compile Thread program 
+instance Compile P.Thread  where
+  -- Can add cases for P.ForAll here.
+  -- Turn into sequential loop. Could be important to make push
+  -- operate uniformly across entire hierarchy.
+  compile s (P.ForAll n f) = 
+    do
+      let (i1,i2) = split2 s
+          nom = "i" ++ show (supplyValue i1)
+          v = variable nom
+          p = f v
+      (a,im) <-  compile i2 p
+
+      return ((),out $ SSeqFor nom (expToIExp n) im)
+                             
+  compile s p = cs s p 
+
+-- Compile Warp program
+instance Compile P.Warp where
+  compile s (P.DistrPar n f) =
+    error "Currently not supported to distribute over the threads, use ForAll instead!"
+  compile s (P.ForAll n f) = do
+    let p = f (variable "warpIx") 
+    (a,im) <- compile s p 
+    return (a, out $ SForAll Warp (expToIExp n) im)
+    --undefined -- compile a warp program that iterates over a space n large
+  compile s (P.Allocate nom n t) = do
+    (Just nw) <- getNWarps -- Must be a Just here, or something is wrong!
+    return ((),out $ SAllocate nom (nw*n) t)
+  compile s (P.Bind p f) = do
+    let (s1,s2) = split2 s
+    (a,im1) <- compile s1 p
+    (b,im2) <- compile s2 (f a)
+    return (b,(im1 ++ im2))
+  compile s (P.Return a) = return (a,[])
+  compile s (P.Identifier) = return (supplyValue s, [])
+
+-- Compile Block program 
+instance Compile P.Block where
+  compile s (P.ForAll n f) = do
+    let nom = "tid"
+        v   = variable nom
+        p   = f v
+    setUsesTid 
+    (a,im) <-  compile s p
+    -- in this case a could be () (since it is guaranteed to be anyway). a
+    return (a,out (SForAll Block (expToIExp n) im))
     
+  compile s (P.DistrPar n'@(Literal n) f) = do
+    
+    {- Distribute work over warps! -} 
+    enterWarp n
+    -- Need to generate some IM here that the backend can read the
+    -- Number of desired warps from. 
+    (a,im) <- compile s (f (variable "warpID"))
+    return (a, out (SDistrPar Warp (expToIExp n') im))
+  compile s p = cs s p
 
+-- Compile a Grid Program 
+instance Compile P.Grid where
+  compile s (P.ForAll n f) = do
+    
+    -- Incorrect, need to compute global thread ids and apply 
+    let p = f gid -- (BlockIdx X)
+        gid = variable "gid" 
+    (a,im) <- compile s p 
+    return (a, out (SForAll Grid (expToIExp n) im))
+
+  {- Distribute over blocks -}
+  compile s (P.DistrPar n f) = do
+    -- Need to generate IM here that the backend can read desired number of blocks from
+    let p = f (BlockIdx X) 
+    
+    (a, im) <- compile s p -- (f (BlockIdx X)) 
+    return (a, out (SDistrPar Block (expToIExp n) im))
+  compile s p = cs s p
+
+
+
+
 ---------------------------------------------------------------------------
+-- General compilation
+---------------------------------------------------------------------------
+cs :: forall t a . Compile t => Supply Int -> P.Program t a -> CM (a,IM) 
+cs i P.Identifier = return $ (supplyValue i, [])
+cs i (P.Assign name ix e) =
+  return $ ((),out (SAssign (IVar name (typeOf e)) (map expToIExp ix) (expToIExp e)))
+
+cs i (P.AtomicOp name ix atom) =
+  case atom of
+    AtomicInc -> return $ ((),out (SAtomicOp (IVar name Word32) (expToIExp ix) AtInc))
+    AtomicAdd e -> undefined
+    AtomicSub e -> undefined
+    AtomicExch e -> undefined 
+  -- (vres, out im)
+  -- where
+  --   res = "a" ++ show (supplyValue i)
+  --   vres = IVar res (typeOf (undefined :: a))
+  --   vname = IVar name (typeOf (undefined :: a))
+  --   im = SAtomicOp vres vname ix atom
+
+--cs i (P.AtomicOp name ix at) = (v,out im)
+--  where 
+--    nom = "a" ++ show (supplyValue i)
+--    v = variable nom
+--    im = SAtomicOp nom name ix at
+      
+cs i (P.Cond bexp p) = do
+  ((),im) <-  compile i p
+  return ((),out (SCond (expToIExp bexp) im)) 
+ 
+
+cs i (P.SeqFor n f) = do
+  let (i1,i2) = split2 i
+      nom = "i" ++ show (supplyValue i1)
+      v = variable nom
+      p = f v
+  (a,im) <-  compile i2 p
+  
+  return (a,out (SSeqFor nom (expToIExp n) im))
+
+    
+cs i (P.SeqWhile b p) = do
+  (a,im) <-  compile i p
+  return (a, out (SSeqWhile (expToIExp b) im))
+
+    
+
+cs i (P.Break) = return ((), out SBreak)
+
+cs i (P.Allocate id n t) = return ((),out (SAllocate id n t))
+cs i (P.Declare  id t)   = return ((),out (SDeclare id t))
+
+cs i (P.Sync)            = return ((),out (SSynchronize))
+
+
+cs i (P.Bind p f) = do 
+  let (s1,s2) = split2 i
+  (a,im1) <- compile s1 p
+  (b,im2) <- compile s2 (f a)
+
+  return (b,im1 ++ im2) 
+ 
+ 
+cs i (P.Return a) = return (a,[])
+
+
+-- Unhandled cases 
+cs i p = error $ "#Program.hs# unhandled in cs: " ++ P.printPrg p -- compile i p 
+
+---------------------------------------------------------------------------
 -- Turning IM to strings
 ---------------------------------------------------------------------------
 
@@ -173,19 +319,17 @@
 -- Print a Statement with metadata 
 printStm :: Show a => (Statement a,a) -> String
 printStm (SAssign name [] e,m) =
-  name ++ " = " ++ printExp e ++ ";" ++ meta m
+  show name ++ " = " ++ show e ++ ";" ++ meta m
 printStm (SAssign name ix e,m) =
-  name ++ "[" ++ concat (intersperse "," (map printExp ix)) ++ "]" ++
-  " = " ++ printExp e ++ ";" ++ meta m
-printStm (SAtomicOp res arr ix op,m) =
-  res ++ " = " ++
-  printAtomic op ++ "(" ++ arr ++ "[" ++ printExp ix ++ "]);" ++ meta m
+  show name ++ "[" ++ concat (intersperse "," (map show ix)) ++ "]" ++
+  " = " ++ show e ++ ";" ++ meta m
+--printStm (SAtomicOp res arr ix op,m) =
+--  res ++ " = " ++
+--  printAtomic op ++ "(" ++ arr ++ "[" ++ show ix ++ "]);" ++ meta m
 printStm (SAllocate name n t,m) =
   name ++ " = malloc(" ++ show n ++ ");" ++ meta m
 printStm (SDeclare name t,m) =
   show t ++ " " ++ name ++ ";" ++ meta m
-printStm (SOutput name t,m) =
-  show t ++ " " ++ name ++ ";" ++ meta m
 printStm (SCond bexp im,m) =
   "if " ++ show bexp ++ "{\n" ++ 
   concatMap printStm im ++ "\n};" ++ meta m
@@ -197,20 +341,41 @@
   "for " ++ name  ++ " in [0.." ++ show n ++"] do" ++ meta m ++ 
   concatMap printStm im ++ "\ndone;\n"
 
-printStm (SForAll n im,m) =
-  "forAll i in [0.." ++ show n ++"] do" ++ meta m ++
+
+printStm (SForAll Warp n im,m) =
+  "forAll wid" ++ "  in [0.." ++ show n ++"] do" ++ meta m ++
   concatMap printStm im ++ "\ndone;\n"
 
-printStm (SForAllBlocks n im,m) =
-  "forAllBlocks i in [0.." ++ show n ++"] do" ++ meta m ++
+printStm (SForAll Block n im,m) =
+  "forAll tid" ++ "  in [0.." ++ show n ++"] do" ++ meta m ++
   concatMap printStm im ++ "\ndone;\n"
-printStm (SForAllThreads n im,m) =
-  "forAllThreads i in [0.." ++ show n ++"] do" ++ meta m ++ 
+
+printStm (SForAll Grid n im,m) =
+  "forAll gid in [0.." ++ show n ++"] do" ++ meta m ++
   concatMap printStm im ++ "\ndone;\n"
 
+printStm (SDistrPar lvl n im,m) = 
+  "forAll gid in [0.." ++ show n ++"] do" ++ meta m ++
+  concatMap printStm im ++ "\ndone;\n"
 
+-- printStm (SWarpForAll n im ,m) =
+--   "forAll(InWarp) tid" ++ "  in [0.." ++ show n ++"] do" ++ meta m ++
+--   concatMap printStm im ++ "\ndone;\n"
+
+-- printStm (SNWarps n im,m) = "Run " ++ show n++ " Warps {\n" ++
+--                             printIM im ++ "\n }"
+--printStm (SForAllThreads n im,m) =
+--  "forAllThreads i in [0.." ++ show n ++"] do" ++ meta m ++ 
+--  concatMap printStm im ++ "\ndone;\n"
+
+
   
 -- printStm (a,m) = error $ show m 
 
 meta :: Show a => a -> String
 meta m = "\t//" ++ show m ++ "\n" 
+
+
+
+      
+
diff --git a/Obsidian/CodeGen/Reify.hs b/Obsidian/CodeGen/Reify.hs
new file mode 100644
--- /dev/null
+++ b/Obsidian/CodeGen/Reify.hs
@@ -0,0 +1,215 @@
+{-# LANGUAGE FlexibleInstances,
+             OverlappingInstances,
+             UndecidableInstances,
+             FlexibleContexts,
+             MultiParamTypeClasses,
+             TypeOperators,
+             TypeFamilies ,
+             ScopedTypeVariables
+             #-}
+
+{- Joel Svensson 2012, 2013
+   Niklas Ulvinge 2013
+
+ 
+-} 
+
+
+module Obsidian.CodeGen.Reify (ToProgram(..)) where 
+
+import Obsidian.Exp 
+import Obsidian.Array
+import Obsidian.Mutable 
+
+import Obsidian.Types
+import Obsidian.Globs 
+import Obsidian.Program
+import Obsidian.Force
+import Obsidian.Memory
+
+import Obsidian.Names
+import Obsidian.Library
+
+import qualified Obsidian.CodeGen.Program as CG
+import Obsidian.CodeGen.CompileIM
+
+import Data.Word
+import Data.Int
+import qualified Data.Map as M 
+      
+---------------------------------------------------------------------------
+-- New approach (hopefully)
+---------------------------------------------------------------------------
+-- "reify" Haskell functions into CG.Programs
+
+
+---------------------------------------------------------------------------
+--
+--------------------------------------------------------------------------- 
+class ToProgram a where
+  toProgram :: Int -> a -> InputList a -> (Parameters,CG.IM)
+  toProgram_ :: Int -> a -> (Parameters, CG.IM)
+
+
+typeOf_ a = typeOf (Literal a)
+
+
+---------------------------------------------------------------------------
+-- Base cases
+---------------------------------------------------------------------------
+
+-- This instance is incorrect
+instance ToProgram (GProgram ()) where
+  -- toProgram i prg () = toProgram $ pJoin prg
+  toProgram i prg () = ([],CG.compileStep1 prg) 
+  -- Needs to deal with GProgram () and GProgram (Push a), GProgram (Pull a)
+  -- in different ways.
+
+  toProgram_ i prg = ([],CG.compileStep1 prg) 
+  
+-- This instance might fix the problem with empty kernels being generated
+instance (ToProgram (Push Grid l a)) => ToProgram (GProgram (Push Grid l a)) where
+  toProgram i p a = toProgram i (runPush p) a
+
+  toProgram_ i p = toProgram_ i (runPush p) 
+
+-- No ToProgram (GProgram (Pull a)) instance is needed. These programs
+-- cannot currently be created using the API. The reason is that GProgram (Pull a)
+-- implies a capability that GPUs do not have. The pulling from an array computed globally.
+-- That kind of computation can not be synced and its result would be undefined. 
+
+
+instance Scalar a => ToProgram (Push Grid l (Exp a)) where
+  toProgram i p a =
+    let outT = Pointer $ typeOf_ (undefined :: a)
+        outN = "output" ++ show i
+        
+        prg = p <: assignOut outN
+        
+        (inputs,im) = toProgram (i+1) prg a
+        
+    in (inputs++[(outN,outT)],im) 
+     
+    where
+      assignOut out a ix = Assign out [ix] a
+  toProgram_ i p = toProgram i p () 
+
+instance (Scalar a, Scalar b) => ToProgram (Push Grid l (Exp a,Exp b)) where
+  toProgram i p a =
+    let   outT1 = Pointer $ typeOf_ (undefined :: a)
+          outT2 = Pointer $ typeOf_ (undefined :: b)
+          outN1 = "output" ++ show i
+          outN2 = "output" ++ show (i+1)
+          
+
+          prg = p <: assignOut (outN1,outN2) 
+            
+          (inputs,im) = toProgram (i+2) prg a
+          
+    in (inputs++[(outN1,outT1),(outN2,outT2)],im)
+    where
+      assignOut (o1,o2) (a,b) ix =
+        do
+          Assign o1 [ix] a
+          Assign o2 [ix] b
+  toProgram_ i p = toProgram i p () 
+
+---------------------------------------------------------------------------
+-- Recursive
+---------------------------------------------------------------------------
+          
+instance (ToProgram b, Scalar t) => ToProgram (Pull EWord32 (Exp t) -> b) where
+  toProgram i f (a :- rest) = ((nom,Pointer t):(n,Word32):ins,prg)
+    where
+      (ins,prg) = toProgram (i+1) (f input) rest
+      nom  = "input" ++ show i
+      n    = "n" ++ show i
+      lengthVar = variable n
+      input = namedGlobal nom lengthVar
+      t     = typeOf_ (undefined :: t)
+  toProgram_ i f = ((nom,Pointer t):(n,Word32):ins,prg)
+    where
+      (ins,prg) = toProgram_ (i+1) (f input)
+      nom  = "input" ++ show i
+      n    = "n" ++ show i
+      lengthVar = variable n
+      input = namedGlobal nom lengthVar
+      t     = typeOf_ (undefined :: t)
+
+instance (ToProgram b, Scalar t) => ToProgram (Pull Word32 (Exp t) -> b) where
+  toProgram i f (a :- rest) = ((nom,Pointer t):ins,prg)
+    where
+      (ins,prg) = toProgram (i+1) (f input) rest
+      nom  = "input" ++ show i
+      input = namedGlobal nom  (len a) 
+      t     = typeOf_ (undefined :: t)
+  toProgram_ _ _ = error "toProgram_: static length" 
+
+
+instance (ToProgram b, Scalar t) => ToProgram (Mutable Global EWord32 (Exp t) -> b) where
+  toProgram i f (a :- rest) = ((nom,Pointer t):(n,Word32):ins,prg)
+    where
+      (ins,prg) = toProgram (i+1) (f input) rest
+      nom  = "input" ++ show i
+      n    = "n" ++ show i
+      lengthVar = variable n
+      input = namedMutable nom lengthVar
+      t     = typeOf_ (undefined :: t)
+  toProgram_ i f = ((nom,Pointer t):(n,Word32):ins,prg)
+    where
+      (ins,prg) = toProgram_ (i+1) (f input)
+      nom  = "input" ++ show i
+      n    = "n" ++ show i
+      lengthVar = variable n
+      input = namedMutable nom lengthVar
+      t     = typeOf_ (undefined :: t)
+
+
+instance (ToProgram b, Scalar t) => ToProgram ((Exp t) -> b) where
+  toProgram i f (a :- rest) = ((nom,t):ins,prg)
+    where
+      (ins,prg) = toProgram (i+1) (f input) rest
+      nom  = "input" ++ show i
+      input = variable nom -- namedGlobal nom (len a) 
+      t     = typeOf_ (undefined :: t)
+  toProgram_ i f = ((nom,t):ins,prg)
+    where
+      (ins,prg) = toProgram_ (i+1) (f input) 
+      nom  = "input" ++ show i
+      input = variable nom -- namedGlobal nom (len a) 
+      t     = typeOf_ (undefined :: t)
+
+
+
+    
+---------------------------------------------------------------------------
+-- heterogeneous lists of inputs 
+---------------------------------------------------------------------------
+data head :- tail = head :- tail
+
+infixr 5 :-
+
+
+---------------------------------------------------------------------------
+-- Function types to input list types. 
+---------------------------------------------------------------------------
+
+type family InputList a
+
+type instance InputList (a -> b)        = a :- (InputList b)
+type instance InputList (Push Grid l b) = ()
+type instance InputList (GProgram b)    = () 
+
+-- genKernelSM :: ToProgram prg => Word32 -> String -> prg -> (String, Word32)
+-- genKernelSM = genKernelSpecsNL
+
+-- genKernelSpecsNL :: ToProgram prg => Word32 -> String -> prg -> (String, Word32)
+-- genKernelSpecsNL nt kn prg = (prgStr,bytesShared) 
+--   where
+--     prgStr = pretty 75 $ ppr $ compile PlatformCUDA (Config nt bytesShared) kn (a,rim) 
+--     (a,im) = toProgram_ 0 prg
+--     iml = computeLiveness im
+--     (m,mm) = mmIM iml sharedMem (M.empty)
+--     bytesShared = size m 
+--     rim = renameIM mm iml
+
diff --git a/Obsidian/CodeGen/SPMDC.hs b/Obsidian/CodeGen/SPMDC.hs
deleted file mode 100644
--- a/Obsidian/CodeGen/SPMDC.hs
+++ /dev/null
@@ -1,744 +0,0 @@
-
-{- Joel Svensson 2012,2013 -} 
-module Obsidian.CodeGen.SPMDC where
-
-import Obsidian.Globs
-import Obsidian.DimSpec
-
-import Obsidian.CodeGen.PP
-
-import Data.Word
-import Data.Int
-
-import qualified Data.List as L
-import qualified Data.Map as M
-import qualified Data.Set as S
-import Data.Tuple
-
-import Control.Monad.State
-
-import Data.Maybe
-
--- TODO: Add Atomic ops 
-
----------------------------------------------------------------------------
--- A C LIKE AST (SPMDC - Single Program Multiple Data C)  
---------------------------------------------------------------------------- 
-data Value = IntVal Int         -- allow ? 
-           | Int8Val Int8
-           | Int16Val Int16
-           | Int32Val Int32
-           | Int64Val Int64
-           | FloatVal Float 
-           | DoubleVal Double
-           | WordVal   Word     -- allow ? 
-           | Word8Val  Word8
-           | Word16Val Word16
-           | Word32Val Word32
-           | Word64Val Word64
-           deriving (Eq,Ord,Show)
-             
-data CType = CVoid | CInt | CFloat | CDouble
-           | CInt8 | CInt16 | CInt32 | CInt64 
-           | CWord | CWord8 | CWord16 | CWord32 | CWord64
-           | CPointer CType -- *type
-           | CArray [CExpr] CType -- type[e1][e2][e3]..[en] or type[] 
-           | CQualified CQualifyer CType 
-           deriving (Eq,Ord,Show)
-             
-data CQualifyer = CQualifyerGlobal  -- CUDA: ""           OpenCL: "__global" 
-                | CQualifyerLocal   -- CUDA: ""           OpenCL: "__local"
-                | CQualifyerKernel  -- CUDA: "__global__" OpenCL: "__kernel"  
-                | CQualifyerShared  -- CUDA: "__shared__" OpenCL: "__local" 
-                | CQualifyerExtern  -- extern   
-                | CQualifyerAttrib CQAttribute
-                deriving (Eq,Ord,Show)
-
-data CQAttribute = CAttribAligned Word32
-                   deriving (Eq,Ord,Show)
-
-
-data CExprP e  = CVar Name CType 
-               -- Threads, Blocks, Grids (All of type Word32) 
-               | CBlockIdx  DimSpec 
-               | CThreadIdx DimSpec
-               | CBlockDim  DimSpec
-               | CGridDim   DimSpec
-                 
-               | CLiteral Value CType
-               | CIndex (e,[e]) CType
-               | CCond e e e CType
-               | CBinOp CBinOp e e  CType
-               | CUnOp  CUnOp  e    CType
-               | CFuncExpr Name [e] CType  -- min, max, sin, cos 
-               | CCast e CType             -- cast expr to type 
-               deriving (Eq,Ord,Show)
-cTypeOfP (CVar _ t) = t
-cTypeOfP (CBlockIdx d) = CWord32
-cTypeOfP (CThreadIdx d) = CWord32
-cTypeOfP (CBlockDim d) = CWord32
-cTypeOfP (CGridDim d) = CWord32
-cTypeOfP (CLiteral _ t) = t
-cTypeOfP (CIndex _ t) = t
-cTypeOfP (CCond  _ _ _ t) = t
-cTypeOfP (CBinOp _ _ _ t) = t
-cTypeOfP (CUnOp _ _ t) = t
-cTypeOfP (CFuncExpr _ _ t) = t
-cTypeOfP (CCast _ t) = t
-
-cSizeOf (CExpr (CIndex (e,es) _))  = 1 + max (cSizeOf e) (maximum (map cSizeOf es))
-cSizeOf (CExpr (CCond e1 e2 e3 _)) = 1 + maximum [cSizeOf e1, cSizeOf e2, cSizeOf e3] 
-cSizeOf (CExpr (CFuncExpr _ es _)) = 1 + maximum (map cSizeOf es) 
-cSizeOf (CExpr (CUnOp  _ e _)) = 1 + cSizeOf e 
-cSizeOf (CExpr (CBinOp _ e1 e2 _ )) = 1 + cSizeOf e1 + cSizeOf e2 
-cSizeOf e = 0
-
-
-data CBinOp = CAdd | CSub | CMul | CDiv | CMod  
-            | CEq | CNotEq | CLt | CLEq | CGt | CGEq 
-            | CAnd | COr
-            | CPow
-            | CBitwiseAnd | CBitwiseOr | CBitwiseXor 
-            | CShiftL | CShiftR 
-            deriving (Eq,Ord,Show) 
-                     
-data CUnOp = CBitwiseNeg
-           deriving (Eq,Ord,Show)
-
-data CAtomicOp = CAtomicAdd | CAtomicInc
-               deriving (Eq, Ord, Show) 
-
----------------------------------------------------------------------------
--- SPMDC
----------------------------------------------------------------------------
-data SPMDC = CAssign CExpr [CExpr] CExpr  -- array or scalar assign
-           | CAtomic CAtomicOp CExpr CExpr CExpr 
-           | CDecl CType Name             -- Declare but no assign
-           | CDeclAssign CType Name CExpr -- declare variable and assign a value 
-           | CFunc   Name  [CExpr]                    
-           | CSync                  -- CUDA: "__syncthreads()" OpenCL: "barrier(CLK_LOCAL_MEM_FENCE)"
-           | CThreadFence
-           | CThreadFenceBlock      -- these could be taken care of with a simple
-                                    -- application of the CFunc constructor
-                                    -- but since sync,threadfence etc are special
-                                    -- and might need attention during code gen
-                                    -- I give them specific constructors. 
-           | CFor    Name CExpr [SPMDC]  -- very simple loop for now.
-           | CIf     CExpr [SPMDC] [SPMDC]
-           deriving (Eq,Ord,Show)
-                    
---                                ret_t       param list     body
-data CKernel = CKernel CQualifyer CType Name [(CType,Name)] [SPMDC] 
-             deriving (Eq,Show)
-           
-----------------------------------------------------------------------------
--- CExpr 
-newtype CExpr = CExpr (CExprP CExpr)
-             deriving (Eq,Ord,Show)
-
-cTypeOf (CExpr e) = cTypeOfP e 
-                      
-----------------------------------------------------------------------------                      
--- DAGs
-type NodeID = Int               
-newtype CENode = CENode (CExprP NodeID) 
-               deriving (Show, Ord, Eq)
-                        
-----------------------------------------------------------------------------
--- Helpers 
-
-cexpr1 exp a       = CExpr $ exp a 
-cexpr2 exp a b     = CExpr $ exp a b 
-cexpr3 exp a b c   = CExpr $ exp a b c 
-cexpr4 exp a b c d = CExpr $ exp a b c d  
-
-cWarpSize  = CExpr $ CVar "warpSize" CWord32 
-cBlockIdx  = cexpr1 CBlockIdx
-cThreadIdx = cexpr1 CThreadIdx
-cBlockDim  = cexpr1 CBlockDim
-cGridDim   = cexpr1 CGridDim 
-cVar       = cexpr2 CVar 
-cLiteral   = cexpr2 CLiteral 
-cIndex     = cexpr2 CIndex 
-cCond      = cexpr4 CCond   
-cFuncExpr  = cexpr3 CFuncExpr 
-cBinOp     = cexpr4 CBinOp 
-cUnOp      = cexpr3 CUnOp 
-cCast      = cexpr2 CCast 
-
-cAssign     = CAssign
-cAtomic     = CAtomic 
-cFunc       = CFunc  
-cDecl       = CDecl
-cSync       = CSync
-cThreadFence = CThreadFence
-cThreadFenceBlock = CThreadFenceBlock
-cDeclAssign = CDeclAssign 
-cIf         = CIf 
-cFor        = CFor 
---------------------------------------------------------------------------
--- Printing 
-data PPConfig = PPConfig {ppKernelQ :: String, 
-                          ppGlobalQ :: String, 
-                          ppLocalQ  :: String,
-                          ppSyncLine :: String} 
-
-printCKernel :: PPConfig -> CKernel -> String 
-printCKernel ppc kern = runPP (ppCKernel ppc  kern ) 0 
-
-ppCKernel :: PPConfig -> CKernel -> PP () 
-ppCKernel ppc (CKernel q t nom ins body) = 
-  ppCQual ppc q >> space >> ppCType ppc t >> space >> line nom >> ppCommaSepList ppIns "(" ")" ins >> 
-  begin >> indent >> newline >> 
-  ppSPMDCList ppc body >>  unindent >> newline >>
-  end 
-  where 
-    ppIns (t,nom) = ppCType ppc t >> space >> line nom
-  
-----------------------------------------------------------------------------
-ppCQual ppc CQualifyerGlobal = line$ ppGlobalQ ppc 
-ppCQual ppc CQualifyerLocal  = line$ ppLocalQ ppc 
-ppCQual ppc CQualifyerKernel = line$ ppKernelQ ppc 
-ppCQual ppc CQualifyerExtern = line$ "extern" 
-ppCQual ppc CQualifyerShared = line$ "__shared__" -- should this be same as local ?
-ppCQual ppc (CQualifyerAttrib a) = ppCAttrib ppc a
-
-ppCAttrib ppc (CAttribAligned x) = line$ "__attribute__ ((aligned(" ++ show x ++ ")))" 
-----------------------------------------------------------------------------
-ppCType ppc CVoid    = line "void"
-ppCType ppc CInt     = line "int"
-ppCType ppc CInt8    = line "int8_t"
-ppCType ppc CInt16   = line "int16_t"
-ppCType ppc CInt32   = line "int32_t"
-ppCType ppc CInt64   = line "int64_t"
-ppCType ppc CFloat   = line "float"
-ppCType ppc CDouble  = line "double"            
-ppCType ppc CWord8   = line "uint8_t"
-ppCType ppc CWord16  = line "uint16_t"
-ppCType ppc CWord32  = line "uint32_t"
-ppCType ppc CWord64  = line "uint64_t" 
-ppCType ppc (CPointer t) = ppCType ppc t >> line "*"
-ppCType ppc (CQualified q t) = ppCQual ppc q >> space >> ppCType ppc t
-
--- a hack (whats the correct way to handle C's t[] ?)
--- Breaks down already for a[][], i think.
-ppCTypedName ppc CVoid   nom = line "void" >> space >> line nom
-ppCTypedName ppc CInt    nom = line "int" >> space >> line nom
-ppCTypedName ppc CFloat  nom = line "float" >> space >> line nom
-ppCTypedName ppc CDouble nom = line "double" >> space >> line nom     
-ppCTypedName ppc CWord8  nom = line "uint8_t" >> space >> line nom
-ppCTypedName ppc CWord16 nom = line "uint16_t" >> space >> line nom
-ppCTypedName ppc CWord32 nom = line "uint32_t" >> space >> line nom
-ppCTypedName ppc CWord64 nom = line "uint64_t" >> space >> line nom
-ppCTypedName ppc (CPointer t) nom = ppCType ppc t >> line "*" >> line nom
-ppCTypedName ppc (CArray [] t) nom = ppCType ppc t >> space >> line nom >> line "[]"
-ppCTypedName ppc (CQualified q t) nom = ppCQual ppc q >> space >> ppCTypedName ppc t nom 
-
-----------------------------------------------------------------------------
-ppValue (IntVal i)    = line$ show i
-ppValue (Int8Val i)   = line$ show i
-ppValue (Int16Val i)  = line$ show i
-ppValue (Int32Val i)  = line$ show i
-ppValue (Int64Val i)  = line$ show i
-ppValue (FloatVal f)  = line$ show f 
-ppValue (DoubleVal d) = line$ show d
-ppValue (Word8Val  w) = line$ show w 
-ppValue (Word16Val w) = line$ show w
-ppValue (Word32Val w) = line$ show w
-ppValue (Word64Val w) = line$ show w 
-
-----------------------------------------------------------------------------
-ppBinOp CAdd = line$ "+"
-ppBinOp CSub = line$ "-"
-ppBinOp CMul = line$ "*"
-ppBinOp CDiv = line$ "/"
-ppBinOp CMod = line$ "%" 
-ppBinOp CEq  = line$ "=="
-ppBinOp CLt  = line$ "<" 
-ppBinOp CLEq = line$ "<="
-ppBinOp CGt  = line$ ">" 
-ppBinOp CGEq = line$ ">="
-ppBinOp CNotEq = line$ "/=" 
-ppBinOp CAnd   = line$ "&&"
-ppBinOp COr    = line$ "||" 
-ppBinOp CBitwiseAnd = line$ "&"  
-ppBinOp CBitwiseOr  = line$ "|" 
-ppBinOp CBitwiseXor = line$ "^" 
-ppBinOp CShiftL     = line$ "<<" 
-ppBinOp CShiftR     = line$ ">>"
-                     
-ppUnOp CBitwiseNeg = line$ "~"       
--- May be incorrect.
---ppUnOp CInt32ToWord32 = line$ "(uint32_t)"
---ppUnOp CWord32ToInt32 = line$ "(int32_t)" 
-
----------------------------------------------------------------------------
---
----------------------------------------------------------------------------
-ppCommaSepList ppElt s e xs = 
-  line s >>  
-  sequence_ (L.intersperse (line ",") (commaSepList' xs)) >> line e
-  where 
-    commaSepList' [] = [] 
-    commaSepList' (x:xs) = ppElt x : commaSepList' xs
-  
----------------------------------------------------------------------------
---
----------------------------------------------------------------------------
-ppSPMDCList ppc xs = sequence_ (map (ppSPMDC ppc) xs) 
-
-
-ppSPMDC :: PPConfig -> SPMDC -> PP () 
-ppSPMDC ppc (CAssign e [] expr) =
-  ppCExpr ppc e >> 
-  line " = " >> 
-  ppCExpr ppc expr >> 
-  cTermLn
-ppSPMDC ppc (CAssign e exprs expr) =
-  ppCExpr ppc e >> 
-  ppCommaSepList (ppCExpr ppc) "[" "]" exprs >> 
-  line " = " >> 
-  ppCExpr ppc expr >> 
-  cTermLn
-ppSPMDC ppc (CAtomic op res arr e) =
-  --ppCExpr ppc res >>
-  --line " = " >>
-  ppAtomicOp ppc op >>
-  wrap "(" ")" (ppCExpr ppc arr >> line ", " >> ppCExpr ppc e ) >>
-  cTermLn 
-
-ppSPMDC ppc (CDecl t n) = ppCTypedName ppc t n  >> cTermLn
-ppSPMDC ppc (CDeclAssign t n e) =
-  ppCTypedName ppc t n >>
-  line " = " >>
-  ppCExpr ppc e >> cTermLn
-ppSPMDC ppc (CFunc nom args) =
-  line nom >>
-  ppCommaSepList (ppCExpr ppc) "(" ")" args >> cTermLn
-ppSPMDC ppc  CSync = line (ppSyncLine ppc) >> cTermLn 
-ppSPMDC ppc (CIf e [] []) = return ()
-ppSPMDC ppc (CIf e xs []) =
-  line "if " >> 
-  wrap "(" ")" (ppCExpr ppc e) >> 
-  begin >> indent >> newline  >> 
-  ppSPMDCList ppc xs >>  unindent >> end
-ppSPMDC ppc (CIf e xs ys) =
-  line "if " >> 
-  wrap "(" ")" (ppCExpr ppc e) >> 
-  begin >> indent >> newline >> 
-  ppSPMDCList ppc xs >>  unindent >> end >> 
-  line "else " >> begin >> indent >> newline >> 
-  ppSPMDCList ppc ys >>  unindent >> end
--- TODO: Clean up here
-ppSPMDC ppc (CFor name e s) =
-  line "for " >>
-  wrap "(" ")" (line ("int " ++ name ++ " = 0;") >>
-                line (name ++ " < ") >> (ppCExpr ppc e) >>
-                line (";") >> line (name ++ "++")) >>
-  begin >> indent >> newline >> 
-  ppSPMDCList ppc s >> unindent >> end
-
-
-ppAtomicOp :: PPConfig -> CAtomicOp -> PP ()
-ppAtomicOp ppc CAtomicInc = line "atomicInc" 
-
-----------------------------------------------------------------------------
---
-ppCExpr :: PPConfig -> CExpr -> PP ()  
--- Cheat and do CUDA print for now!
-  -- should do lookup in PPConfig and figure out how to 
-  -- print these for CUDA/OpenCL
-ppCExpr ppc (CExpr (CBlockIdx X)) = line "blockIdx.x" 
-ppCExpr ppc (CExpr (CBlockIdx Y)) = line "blockIdx.y" 
-ppCExpr ppc (CExpr (CBlockIdx Z)) = line "blockIdx.z" 
-ppCExpr ppc (CExpr (CThreadIdx X)) = line "threadIdx.x" 
-ppCExpr ppc (CExpr (CThreadIdx Y)) = line "threadIdx.y" 
-ppCExpr ppc (CExpr (CThreadIdx Z)) = line "threadIdx.z" 
-ppCExpr ppc (CExpr (CBlockDim X)) = line "blockDim.x" 
-ppCExpr ppc (CExpr (CBlockDim Y)) = line "blockDim.y" 
-ppCExpr ppc (CExpr (CBlockDim Z)) = line "blockDim.z" 
-ppCExpr ppc (CExpr (CGridDim X)) = line "gridDim.x" 
-ppCExpr ppc (CExpr (CGridDim Y)) = line "gridDim.y" 
-ppCExpr ppc (CExpr (CGridDim Z)) = line "gridDim.z" 
-
-ppCExpr ppc (CExpr (CVar nom _)) = line nom
-ppCExpr ppc (CExpr (CLiteral v _)) = ppValue v 
-ppCExpr ppc (CExpr (CIndex (e,[]) _)) = ppCExpr ppc e 
-ppCExpr ppc (CExpr (CIndex (e,xs) _)) =
-  ppCExpr ppc e  >>  
-  ppCommaSepList (ppCExpr ppc) "[" "]" xs
-ppCExpr ppc (CExpr (CCond e1 e2 e3 _)) =
-  wrap "(" ")" 
-  (ppCExpr ppc e1 >> 
-   line " ? " >> 
-   ppCExpr ppc e2 >> 
-   line " : " >>  
-   ppCExpr ppc e3
-  )
-ppCExpr ppc (CExpr (CBinOp bop e1 e2 _)) =
-  wrap "(" ")"
-  (
-   ppCExpr ppc e1 >> 
-   ppBinOp bop >> 
-   ppCExpr ppc e2 
-  ) 
-ppCExpr ppc (CExpr (CUnOp  uop  e _)) =
-  wrap "(" ")" 
-  (
-   ppUnOp uop >> 
-   ppCExpr ppc e 
-  )
-ppCExpr ppc (CExpr (CFuncExpr nom args _)) =
-  line nom >> 
-  ppCommaSepList (ppCExpr ppc) "(" ")" args
-ppCExpr ppc (CExpr (CCast e t)) =
-  line "((" >> 
-  ppCType ppc t >> 
-  line ")" >> 
-  ppCExpr ppc e >> 
-  line ")"
-
----------------------------------------------------------------------------
--- Optimize for complicated indexing expressions
----------------------------------------------------------------------------
-
--- TODO: #1: Discover all expressions that represent an index into an array
---       #2: Count usages of them (is more complicated than expected)
---       #3: For "Complicated" expressions used more than once
---           declare a new name for the index and compute it once. (if not data dependent) 
---
---       Possible approach is two passes over the SPMDC structure.
---       The first discovers expressions
---         The in-between create small SPMDC code that declares variables. 
---       The second replaces some of them by a variable
---
-
--- Assign with all expressions an integer 
-type ExpMap = M.Map CExpr (Int,Int) 
-
-type Decl = (Int,CExpr) 
-
-
--- Insert, but only if size is right! 
-insert :: CExpr -> State (Int,ExpMap) () 
-insert e | cSizeOf e >= 2 =
-  do
-    (i,m) <- get
-    case M.lookup e m of
-      (Just (id,count)) ->
-        do
-          let m' = M.insert e (id,count+1) m
-          put (i,m')
-      Nothing           ->
-        do
-          let m' = M.insert e (i,1) m
-          put (i+1,m')
-insert e = return () 
-
--- Decide if an expression is safe or not to move to
--- function prelude.
--- Simply put, it checks for any data dependency.
--- (This code is unused! ) 
-safeExp :: S.Set Name -> CExpr -> Bool
-safeExp s (CExpr (CVar name _)) = S.member name s
-safeExp s (CExpr (CIndex (e,es) _)) = safeExp s e && all (safeExp s) es
-safeExp s (CExpr (CCond e1 e2 e3 _)) = safeExp s e1 && safeExp s e2 && safeExp s e3
-safeExp s (CExpr (CBinOp _ e1 e2 _)) = safeExp s e2 && safeExp s e2
-safeExp s (CExpr (CUnOp _ e _)) = safeExp s e
-safeExp s (CExpr (CFuncExpr _ es _)) = all (safeExp s) es
-safeExp s (CExpr (CCast e _)) = safeExp s e 
-safeExp _ _ = True 
-          
-collectExps :: [SPMDC] -> State (Int,ExpMap) () 
-collectExps sp =  mapM_ process sp
-  where
-    process (CAssign _ ixs e) =
-      do
-        mapM_ processE ixs 
-        processE e
-    process (CDeclAssign _ _ e) = processE e
-    process (CFunc _ es) = mapM_ processE es
-    process (CFor  _ e sp) =
-      do 
-        processE e
-        collectExps sp
-    process (CIf bexp sp1 sp2) =
-      do
-        processE bexp
-        collectExps sp1
-        collectExps sp2 
-    process a = return () 
-
-
-    processE (CExpr (CVar _ _))     = return () -- too simple
-    processE (CExpr (CBlockIdx d))  = return () 
-    processE (CExpr (CThreadIdx d)) = return ()
-    processE (CExpr (CBlockDim d))  = return ()
-    processE (CExpr (CGridDim d))   = return ()
-    processE (CExpr (CLiteral _ _)) = return ()
-    processE e@(CExpr (CIndex (e1,es) _)) =
-      do 
-        -- insert e
-        processE e1
-        mapM_ processE es
-    processE e@(CExpr (CCond e1 e2 e3 _)) =
-      do
-        insert e
-        mapM_ processE [e1,e2,e3]
-    processE e@(CExpr (CBinOp _ e1 e2 _)) =
-      do
-        insert e
-        processE e1
-        processE e2
-    processE e@(CExpr (CUnOp _ e1 _)) =
-      do
-        insert e
-        processE e1
-    processE e@(CExpr (CFuncExpr _ es _)) =
-      do
-        insert e
-        mapM_ processE es
-    processE e@(CExpr (CCast e1 _)) =
-      do
-        -- refine this step. Only insert if e1 is nonsimple
-        insert e
-        processE e1
-    
-
-
--- REMEMBER TO KEEP IT SIMPLE.
-replacePass :: ExpMap -> [SPMDC] -> ([Decl],[SPMDC])
-replacePass _ []     = ([],[])
-replacePass m (x:xs) = let (decls,x') = process m x
-                           (rest, xs') = replacePass m xs
-                         
-                       in  (L.nubBy fstEq (decls ++ rest), x':xs')
-  where
-    fstEq :: (Int,a) -> (Int,a) -> Bool
-    fstEq a b = fst a == fst b
-
-    process m (CFor name e sp) = (decls,CFor name e' sp')
-      where
-        (decls1, e') = processE m e
-        (decls2, sp') = replacePass m sp
-        decls = L.nubBy fstEq (decls1++decls2) 
-    process m (CAssign name es e) = (decls,CAssign name es' e')  
-      where
-        (decls1,es') = processEList m es
-        (decls2,e')  = processE m e
-        decls = L.nubBy fstEq (decls1 ++ decls2)
-    process m s = ([],s)    
-
-    processEList m [] = ([],[])
-    processEList m (e:es) =
-      let (decls1,e') = processE m e
-          (decls2,es') = processEList m es
-      in  (L.nubBy fstEq (decls1 ++ decls2),e':es')
-
-
-    processE :: ExpMap ->  CExpr -> ([Decl],CExpr)
-    processE m e@(CExpr (CIndex (e1,es) t)) =
-      case M.lookup e m of
-        Nothing ->
-          let (d1,es') = processEList m es
-          in (L.nubBy fstEq d1, CExpr (CIndex (e1,es') t))
-           
-        (Just _) -> error "Just in CIndex case"
-
-    processE m e@(CExpr (CCond e1 e2 e3 t)) =
-      case M.lookup e m of
-        Nothing ->
-          let 
-            (d1,e1') = processE m e1
-            (d2,e2') = processE m e2
-            (d3,e3') = processE m e3
-          in (L.nubBy fstEq (d1++d2++d3), CExpr (CCond e1' e2' e3' t))
-        Just (id,1) ->
-          let 
-            (d1,e1') = processE m e1
-            (d2,e2') = processE m e2
-            (d3,e3') = processE m e3
-          in (L.nubBy fstEq (d1++d2++d3), CExpr (CCond e1' e2' e3' t))
-        Just (id,n) -> 
-          ([(id,e)],CExpr (CVar ("t" ++ show id) (cTypeOf e)))
-        
-    processE m e@(CExpr (CBinOp op e1 e2 t))  =
-      case M.lookup e m of
-        Nothing -> 
-           let (d1,e1') = processE m e1
-               (d2,e2') = processE m e2
-           in (L.nubBy fstEq (d1++d2), CExpr (CBinOp op e1' e2' t))
-             
-        (Just (id,1)) -> 
-           let (d1,e1') = processE m e1
-               (d2,e2') = processE m e2
-           in (L.nubBy fstEq (d1++d2), CExpr (CBinOp op e1' e2' t))
-          
-        (Just (id,n)) ->
-          --let (decls1,e1') = processE m e1
-          --    (decls2,e2') = processE m e2
-          --    e' = CExpr (CBinOp op e1' e2' t)
-          --in (L.nubBy fstEq (decls1++decls2++[(id,e')]),CExpr (CVar ("t" ++ show id) (cTypeOf e)))
-          ([(id,e)],CExpr (CVar ("t" ++ show id) (cTypeOf e)))
-          --let (decls,e1') = performCSE e
-          --in  (decls,e1')
-    processE m e@(CExpr (CUnOp op e1 t)) =
-      case M.lookup e m of
-        Nothing -> -- e is not a candidate for hoisting
-          let (d1,e1') = processE m e1
-          in (L.nubBy fstEq d1,CExpr (CUnOp op e1' t))
-        Just (id,1) -> -- e occurs only once, do not replace
-          let (d1,e1') = processE m e1
-          in (L.nubBy fstEq d1,CExpr (CUnOp op e1' t))
-        Just (id,n) -> -- e occurs n times, replace it! 
-          ([(id,e)],CExpr (CVar ("t" ++ show id) (cTypeOf e)))
-    processE m e@(CExpr (CCast e1 t)) = (id,CExpr (CCast e1' t))
-      where
-        (id,e1') = processE m e1 
-
-    processE m e =
-      case M.lookup e m of
-        Nothing -> ([],e)
-        (Just (id,1)) -> ([],e)
-        (Just (id,n)) -> ([(id,e)],CExpr (CVar ("t" ++ show id) (cTypeOf e)))
-    
-
-
-
-declsToSPMDC :: [Decl] -> [SPMDC]
-declsToSPMDC decls = map process decls
-  where
-    process (i,e) = CDeclAssign (cTypeOf e) ("t" ++ show i) e 
-
-
-
----------------------------------------------------------------------------
--- Custom form of CSE 
--- 
----------------------------------------------------------------------------
-
-
-performCSE :: CExpr -> ([Decl],CExpr)
-performCSE exp = dagToExp dag a 
-  where (a,(i,dag,cpd)) = runState (buildDAG exp) (0,M.empty,M.empty)
-        
-
-type DAG = M.Map NodeID CENode
-type CPD = M.Map CENode NodeID
-
---isComputed :: [Decl] -> CExpr -> Maybe NodeID
---isComputed c e = L.lookup e (map swap c) 
-
--- ensure that node ids stay separate from already generated names.
--- By using the integer from the collectPass as initial state
-buildDAG :: CExpr -> State (Int,DAG,CPD) NodeID
-buildDAG (CExpr (CVar n t)) =
-  addNode (CENode (CVar n t))
-buildDAG (CExpr (CBlockIdx d)) =
-  addNode (CENode (CBlockIdx d))
-buildDAG (CExpr (CThreadIdx d)) =
-  addNode (CENode (CThreadIdx d))
-buildDAG (CExpr (CBlockDim d)) =
-  addNode (CENode (CBlockDim d))
-buildDAG (CExpr (CGridDim d)) =
-  addNode (CENode (CGridDim d))
-buildDAG (CExpr (CLiteral v t)) =
-  addNode (CENode (CLiteral v t))
-buildDAG (CExpr (CIndex (e,es) t)) =
-  do
-    e' <- buildDAG e
-    es' <- mapM buildDAG es
-    addNode (CENode (CIndex (e',es') t))
-buildDAG (CExpr (CCond e1 e2 e3 t)) =
-  do
-    [e1',e2',e3'] <- mapM buildDAG [e1,e2,e3]
-    addNode (CENode (CCond e1' e2' e3' t))
-buildDAG (CExpr (CBinOp op e1 e2 t)) =
-  do
-    [e1',e2'] <- mapM buildDAG [e1,e2]
-    addNode (CENode (CBinOp op e1' e2' t))
-buildDAG (CExpr (CUnOp op e t)) =
-  do 
-    e' <- buildDAG e
-    addNode (CENode (CUnOp op e' t))
-buildDAG (CExpr (CFuncExpr n es t)) =
-  do
-    es' <- mapM buildDAG es
-    addNode (CENode (CFuncExpr n es' t))
-buildDAG (CExpr (CCast e t)) =
-  do
-    e' <- buildDAG e
-    addNode (CENode (CCast e' t)) 
-  
-
-addNode :: CENode -> State (Int, DAG, CPD) NodeID
-addNode node = do
-  (i,d,c) <- get
-  case M.lookup node c of
-    Nothing ->
-      do 
-        put (i+1,M.insert i node d,M.insert node i c)
-        return i
-    (Just n) -> return n
-    
----------------------------------------------------------------------------
--- Reconstruct a CExpr + [decls]
----------------------------------------------------------------------------
-
-dagToExp :: DAG -> NodeID -> ([Decl], CExpr)
-dagToExp dag nid =
-  case M.lookup nid dag of
-    Nothing -> error "dagToExp: Broken DAG"
-    (Just (CENode (CVar nom t)))   -> ([],CExpr (CVar nom t))
-    (Just (CENode (CBlockIdx d)))  -> ([],CExpr (CBlockIdx d))
-    (Just (CENode (CThreadIdx d))) -> ([],CExpr (CThreadIdx d)) 
-    (Just (CENode (CGridDim d)))   -> ([],CExpr (CGridDim d))
-    (Just (CENode (CLiteral v t))) -> ([],CExpr (CLiteral v t))
-
-    -- Bit of a special case 
-    (Just (CENode (CIndex (e,es) t))) ->
-      let
-        (d,e') = dagToExp dag e
-        r       = map (dagToExp dag) es
-        ds      = concatMap fst r
-        es'     = map snd r
-      in (d ++ ds, CExpr (CIndex (e',es') t))
-
-    -- Normal cases      
-    (Just (CENode (CCond e1 e2 e3 t))) ->
-      let (d1,e1') = dagToExp dag e1
-          (d2,e2') = dagToExp dag e2
-          (d3,e3') = dagToExp dag e3
-          exp = CExpr (CCond e1' e2' e3' t)
-          this = (nid,tmp nid t)
-      in (d1++d2++d3++[this], exp)
-    (Just (CENode (CBinOp op e1 e2 t))) ->
-      let (d1,e1') = dagToExp dag e1
-          (d2,e2') = dagToExp dag e2
-          exp = CExpr (CBinOp op e1' e2' t)
-          this = (nid,exp)
-      in (d1++d2++[this],tmp nid t)
-    (Just (CENode (CUnOp op e t))) ->
-      let (d,e') = dagToExp dag e
-          exp = CExpr (CUnOp op e' t)
-          this = (nid, exp)
-      in (d++[this],tmp nid t)
-    (Just (CENode (CFuncExpr nom es t))) ->
-      let r = map (dagToExp dag) es
-          ds = concatMap fst r
-          es' = map snd r
-          exp = CExpr (CFuncExpr nom es' t)
-          this = (nid,exp)
-      in (ds++[this],tmp nid t)
-    -- never pull out just for a cast
-    (Just (CENode (CCast e t))) ->
-      let (d,e') = dagToExp dag e
-          exp = CExpr (CCast e' t)
-          --this = (nid,exp)
-      in (d,exp) 
-
-         
-  
-
-
-tmp nid t = CExpr $ CVar ("t" ++ show nid) t 
diff --git a/Obsidian/Exp.hs b/Obsidian/Exp.hs
--- a/Obsidian/Exp.hs
+++ b/Obsidian/Exp.hs
@@ -3,6 +3,7 @@
              FlexibleContexts,
              FlexibleInstances, 
              UndecidableInstances,
+             OverlappingInstances,
              RankNTypes #-} 
 
 {- Joel Svensson 2012 -} 
@@ -27,11 +28,9 @@
 import Obsidian.Types
 import Obsidian.Globs
 
-import Obsidian.CodeGen.SPMDC
-
 ---------------------------------------------------------------------------
 -- some synonyms
-type Data a = Exp a 
+--type Data a = Exp a 
 
 
 type EInt    = Exp Int      
@@ -47,6 +46,17 @@
 type EWord32  = Exp Word32 
 type EWord64  = Exp Word64 
 
+type EI8   = Exp Int8
+type EI16  = Exp Int16
+type EI32  = Exp Int32
+type EI64  = Exp Int64
+
+type EW8   = Exp Word8   
+type EW16  = Exp Word16 
+type EW32  = Exp Word32 
+type EW64  = Exp Word64 
+
+
 type EFloat  = Exp Float  
 type EDouble = Exp Double 
 type EBool   = Exp Bool    
@@ -56,7 +66,7 @@
 ---------------------------------------------------------------------------
 -- Class Scalar. All the things we can handle code generation for 
 
-class (Eq a, ExpToCExp a, Show a) => Scalar a where 
+class (Eq a, ExpToIExp a, Show a) => Scalar a where 
   sizeOf :: Exp a -> Int   --  
   typeOf :: Exp a -> Type  --   Good enough for me ... 
 
@@ -191,6 +201,8 @@
   -- Boolean 
   And :: Op ((Bool,Bool) -> Bool) 
   Or  :: Op ((Bool,Bool) -> Bool)
+
+  Not :: Op (Bool -> Bool) 
   
   -- Bitwise 
   BitwiseAnd :: Bits a => Op ((a,a) -> a) 
@@ -200,7 +212,7 @@
 
   -- I DO NOT EVEN KNOW WHAT THIS MEANS: work around it! 
   ShiftL     :: forall a b. (Num b, Bits a) => Op ((a, b) -> a)  
-  ShiftR     :: forall a b .(Num b, Bits a) => Op ((a, b) -> a)  
+  ShiftR     :: forall a b. (Num b, Bits a) => Op ((a, b) -> a)  
   
   -- built-ins
   Min        :: Ord a => Op ((a,a) -> a) 
@@ -226,10 +238,12 @@
   ATanH :: Floating a => Op (a -> a) -- "atanhf"
   ACosH :: Floating a => Op (a -> a) -- "acoshf"
   -- There is no "div" in "Num" but it's already defined above. 
-  FDiv :: Floating a => Op ((a, a) -> a) -- "acoshf"
+  FDiv :: Floating a => Op ((a, a) -> a) 
 
   Int32ToWord32 :: Op (Int32 -> Word32)
-  Word32ToInt32 :: Op (Word32 -> Int32) 
+  Word32ToInt32 :: Op (Word32 -> Int32)
+  Word32ToFloat :: Op (Word32 -> Float)
+  Word32ToWord8 :: Op (Word32 -> Word8) 
   
 
 ---------------------------------------------------------------------------
@@ -242,48 +256,28 @@
 warpSize = WarpSize
 
 ---------------------------------------------------------------------------
--- Collect array names
+-- Typecasts
+---------------------------------------------------------------------------
+i32ToW32 = UnOp Int32ToWord32
+w32ToI32 = UnOp Word32ToInt32
 
-collectArrays :: Scalar a => Exp a -> [Name]
-collectArrays (Literal _) = []
-collectArrays (ThreadIdx _) = []
-collectArrays (BlockIdx _) = [] 
-collectArrays (Index (name,[])) = []
-collectArrays (Index (name,_)) = [name]
-collectArrays (BinOp _ e1 e2) = collectArrays e1 ++ collectArrays e2
-collectArrays (UnOp  _ e) = collectArrays e
-collectArrays (If b e1 e2) = collectArrays b ++ 
-                             collectArrays e1 ++ 
-                             collectArrays e2
--- collectArrays a = error $ show a
+w32ToF = UnOp Word32ToFloat
 
-collectArrayIndexPairs :: Scalar a => Exp a -> [(Name,Exp Word32)]
-collectArrayIndexPairs (Literal _) = []
-collectArrayIndexPairs (Index (name,[])) = []
-collectArrayIndexPairs (Index (name,[ix])) = [(name,ix)]
-collectArrayIndexPairs (BinOp _ e1 e2) = collectArrayIndexPairs e1 ++ collectArrayIndexPairs e2
-collectArrayIndexPairs (UnOp  _ e) = collectArrayIndexPairs e
-collectArrayIndexPairs (If b e1 e2) = collectArrayIndexPairs b ++ 
-                                      collectArrayIndexPairs e1 ++ 
-                                      collectArrayIndexPairs e2
+w32ToW8 = UnOp Word32ToWord8
 
 
----------------------------------------------------------------------------
--- Typecasts
----------------------------------------------------------------------------
-int32ToWord32 = UnOp Int32ToWord32
-word32ToInt32 = UnOp Word32ToInt32
 
+
 ---------------------------------------------------------------------------
 -- 
 instance Scalar a => Show (Exp a) where 
   show = printExp 
 
--- Look this over. Do I really need a types expression data type ?
+-- Look this over. Do I really need a typed expression data type ?
 --  (No real need for a Exp GADT I think. Go back to keeping it simple!) 
 instance (Eq a, Scalar a) => Eq (Exp a) where
-  (==) a b = -- error $ "equality test between exps: " ++ show a ++ " " ++ show b --
-    expToCExp a == expToCExp b
+  (==) a b = 
+    expToIExp a == expToIExp b
     -- Maybe not efficient! But simple.
 
   
@@ -292,15 +286,24 @@
     max a b = BinOp Max a b
 
 ---------------------------------------------------------------------------
--- INT Instances
+-- Num instance Exp a?
 ---------------------------------------------------------------------------
-instance Num (Exp Int) where 
+instance (Scalar a ,Num a) => Num (Exp a) where 
   (+) a (Literal 0) = a
   (+) (Literal 0) a = a
   (+) (Literal a) (Literal b) = Literal (a+b)
-  -- Added 2 Oct 2012
+
+  -- Added 15 Jan 2013
+  (+) (BinOp Mul (BinOp Div x (Literal a)) (Literal b))
+       (BinOp Mod y (Literal c))
+        | x == y && a == b && b == c = x 
+      -- This spots the kind of indexing that occurs from 
+      --  converting a bix tix view to and from gix view
+        
+  -- Added 2 oct 2012
   (+) (BinOp Sub b (Literal a)) (Literal c) | a == c  = b 
   (+) (Literal b) (BinOp Sub a (Literal c)) | b == c  = a 
+ 
   (+) a b = BinOp Add a b  
   
   (-) a (Literal 0) = a 
@@ -314,10 +317,53 @@
   (*) (Literal a) (Literal b) = Literal (a*b) 
   (*) a b = BinOp Mul a b 
   
-  signum = error "signum: not implemented for Exp Int" 
-  abs = error "abs: not implemented for Exp Int" 
+  signum = error "signum: not implemented for Exp a"
+  abs = error "abs: not implemented for Exp a" 
   fromInteger a = Literal (fromInteger a) 
+   
+instance (Scalar a, Real a) => Real (Exp a) where 
+  toRational = error "toRational: not implemented for Exp a" 
   
+
+instance (Scalar a, Enum a) => Enum (Exp a) where
+  toEnum = error "toEnum: not implemented for Exp a" 
+  fromEnum = error "fromEnum: not implemented for Exp a" 
+
+instance (Scalar a, Integral a) => Integral (Exp a) where
+  mod (Literal a) (Literal b) = Literal (a `mod` b) 
+  mod a b = BinOp Mod a b
+  div _ (Literal 0) = error "Division by zero in expression" 
+  div a b = BinOp Div a b
+  quotRem   = error "quotRem: not implemented for Exp a" 
+  toInteger = error "toInteger: not implemented for Exp a"
+ 
+---------------------------------------------------------------------------
+-- INT Instances
+---------------------------------------------------------------------------
+-- instance Num (Exp Int) where 
+--   (+) a (Literal 0) = a
+--   (+) (Literal 0) a = a
+--   (+) (Literal a) (Literal b) = Literal (a+b)
+--   -- Added 2 Oct 2012
+--   (+) (BinOp Sub b (Literal a)) (Literal c) | a == c  = b 
+--   (+) (Literal b) (BinOp Sub a (Literal c)) | b == c  = a 
+--   (+) a b = BinOp Add a b  
+  
+--   (-) a (Literal 0) = a 
+--   (-) (Literal a) (Literal b) = Literal (a - b) 
+--   (-) a b = BinOp Sub a b 
+  
+--   (*) a (Literal 1) = a 
+--   (*) (Literal 1) a = a
+--   (*) _ (Literal 0) = Literal 0
+--   (*) (Literal 0) _ = Literal 0
+--   (*) (Literal a) (Literal b) = Literal (a*b) 
+--   (*) a b = BinOp Mul a b 
+  
+--   signum = error "signum: not implemented for Exp Int" 
+--   abs = error "abs: not implemented for Exp Int" 
+--   fromInteger a = Literal (fromInteger a) 
+  
 -- Added new cases for literal 0 (2012/09/25)
 instance Bits (Exp Int) where  
   (.&.) x (Literal 0) = Literal 0
@@ -344,47 +390,48 @@
 
 
 -- TODO: change undefined to some specific error.
-instance Real (Exp Int) where
-  toRational = error "toRational: not implemented for Exp Int)"  
+-- instance Real (Exp Int) where
+--   toRational = error "toRational: not implemented for Exp Int)"  
 
-instance Enum (Exp Int) where
-  toEnum = error "toEnum: not implemented for Exp Int" 
-  fromEnum = error "fromEnum: not implemented for Exp Int"
+-- instance Enum (Exp Int) where
+--   toEnum = error "toEnum: not implemented for Exp Int" 
+--   fromEnum = error "fromEnum: not implemented for Exp Int"
          
-instance Integral (Exp Int) where
-  mod (Literal a) (Literal b) = Literal (a `mod` b) 
-  mod a b = BinOp Mod a b
-  div _ (Literal 0) = error "Division by zero in expression" 
-  div a b = BinOp Div a b
-  quotRem = error "quotRem: not implemented for Exp Int" 
-  toInteger = error "toInteger: not implemented for Exp Int" 
+-- instance Integral (Exp Int) where
+--   mod (Literal a) (Literal b) = Literal (a `mod` b) 
+--   mod a b = BinOp Mod a b
+--   div _ (Literal 0) = error "Division by zero in expression" 
+--   div a b = BinOp Div a b
+--   quotRem = error "quotRem: not implemented for Exp Int" 
+--   toInteger = error "toInteger: not implemented for Exp Int" 
 
+
 ---------------------------------------------------------------------------
 -- Int32
 ---------------------------------------------------------------------------
-instance Num (Exp Int32) where 
-  (+) a (Literal 0) = a
-  (+) (Literal 0) a = a
-  (+) (Literal a) (Literal b) = Literal (a+b)
-  -- Added 2 Oct 2012
-  (+) (BinOp Sub b (Literal a)) (Literal c) | a == c  = b 
-  (+) (Literal b) (BinOp Sub a (Literal c)) | b == c  = a 
-  (+) a b = BinOp Add a b  
+-- instance Num (Exp Int32) where 
+--   (+) a (Literal 0) = a
+--   (+) (Literal 0) a = a
+--   (+) (Literal a) (Literal b) = Literal (a+b)
+--   -- Added 2 Oct 2012
+--   (+) (BinOp Sub b (Literal a)) (Literal c) | a == c  = b 
+--   (+) (Literal b) (BinOp Sub a (Literal c)) | b == c  = a 
+--   (+) a b = BinOp Add a b  
   
-  (-) a (Literal 0) = a 
-  (-) (Literal a) (Literal b) = Literal (a - b) 
-  (-) a b = BinOp Sub a b 
+--   (-) a (Literal 0) = a 
+--   (-) (Literal a) (Literal b) = Literal (a - b) 
+--   (-) a b = BinOp Sub a b 
   
-  (*) a (Literal 1) = a 
-  (*) (Literal 1) a = a
-  (*) _ (Literal 0) = 0
-  (*) (Literal 0) _ = 0 
-  (*) (Literal a) (Literal b) = Literal (a*b) 
-  (*) a b = BinOp Mul a b 
+--   (*) a (Literal 1) = a 
+--   (*) (Literal 1) a = a
+--   (*) _ (Literal 0) = 0
+--   (*) (Literal 0) _ = 0 
+--   (*) (Literal a) (Literal b) = Literal (a*b) 
+--   (*) a b = BinOp Mul a b 
   
-  signum = error "signum: not implemented for Exp Int32"
-  abs = error "abs: not implemented for Exp Int32" 
-  fromInteger a = Literal (fromInteger a) 
+--   signum = error "signum: not implemented for Exp Int32"
+--   abs = error "abs: not implemented for Exp Int32" 
+--   fromInteger a = Literal (fromInteger a) 
   
 -- Added new cases for literal 0 (2012/09/25)
 instance Bits (Exp Int32) where  
@@ -412,57 +459,57 @@
 
 
 -- TODO: change undefined to some specific error.
-instance Real (Exp Int32) where
-  toRational = error "toRational: not implemented for Exp Int32"
+-- instance Real (Exp Int32) where
+--   toRational = error "toRational: not implemented for Exp Int32"
 
-instance Enum (Exp Int32) where
-  toEnum = error "toEnum: not implemented for Exp Int32" 
-  fromEnum = error "fromEnum: not implemented for Exp Int32" 
+-- instance Enum (Exp Int32) where
+--   toEnum = error "toEnum: not implemented for Exp Int32" 
+--   fromEnum = error "fromEnum: not implemented for Exp Int32" 
          
-instance Integral (Exp Int32) where
-  mod (Literal a) (Literal b) = Literal (a `mod` b) 
-  mod a b = BinOp Mod a b
-  div _ (Literal 0) = error "Division by zero in expression" 
-  div a b = BinOp Div a b
-  quotRem = error "quotRem: not implemented for Exp Int32" 
-  toInteger = error "toInteger: not implemented for Exp Int32" 
+-- instance Integral (Exp Int32) where
+--   mod (Literal a) (Literal b) = Literal (a `mod` b) 
+--   mod a b = BinOp Mod a b
+--   div _ (Literal 0) = error "Division by zero in expression" 
+--   div a b = BinOp Div a b
+--   quotRem = error "quotRem: not implemented for Exp Int32" 
+--   toInteger = error "toInteger: not implemented for Exp Int32" 
 
 
 ---------------------------------------------------------------------------
 -- Word32 Instances
 ---------------------------------------------------------------------------
-instance Num (Exp Word32) where 
-  (+) a (Literal 0) = a
-  (+) (Literal 0) a = a
-  (+) (Literal a) (Literal b) = Literal (a+b)
+-- instance Num (Exp Word32) where 
+--   (+) a (Literal 0) = a
+--   (+) (Literal 0) a = a
+--   (+) (Literal a) (Literal b) = Literal (a+b)
 
-  -- Added 15 Jan 2013
-  (+) (BinOp Mul (BinOp Div x (Literal a)) (Literal b))
-       (BinOp Mod y (Literal c))
-        | x == y && a == b && b == c = x 
-      -- This spots the kind of indexing that occurs from 
-      --  converting a bix tix view to and from gix view
+--   -- Added 15 Jan 2013
+--   (+) (BinOp Mul (BinOp Div x (Literal a)) (Literal b))
+--        (BinOp Mod y (Literal c))
+--         | x == y && a == b && b == c = x 
+--       -- This spots the kind of indexing that occurs from 
+--       --  converting a bix tix view to and from gix view
         
-  -- Added 2 oct 2012
-  (+) (BinOp Sub b (Literal a)) (Literal c) | a == c  = b 
-  (+) (Literal b) (BinOp Sub a (Literal c)) | b == c  = a 
+--   -- Added 2 oct 2012
+--   (+) (BinOp Sub b (Literal a)) (Literal c) | a == c  = b 
+--   (+) (Literal b) (BinOp Sub a (Literal c)) | b == c  = a 
  
-  (+) a b = BinOp Add a b  
+--   (+) a b = BinOp Add a b  
   
-  (-) a (Literal 0) = a 
-  (-) (Literal a) (Literal b) = Literal (a - b) 
-  (-) a b = BinOp Sub a b 
+--   (-) a (Literal 0) = a 
+--   (-) (Literal a) (Literal b) = Literal (a - b) 
+--   (-) a b = BinOp Sub a b 
   
-  (*) a (Literal 1) = a 
-  (*) (Literal 1) a = a
-  (*) _ (Literal 0) = Literal 0
-  (*) (Literal 0) _ = Literal 0
-  (*) (Literal a) (Literal b) = Literal (a*b) 
-  (*) a b = BinOp Mul a b 
+--   (*) a (Literal 1) = a 
+--   (*) (Literal 1) a = a
+--   (*) _ (Literal 0) = Literal 0
+--   (*) (Literal 0) _ = Literal 0
+--   (*) (Literal a) (Literal b) = Literal (a*b) 
+--   (*) a b = BinOp Mul a b 
   
-  signum = error "signum: not implemented for Exp Word32"
-  abs = error "abs: not implemented for Exp Word32" 
-  fromInteger a = Literal (fromInteger a) 
+--   signum = error "signum: not implemented for Exp Word32"
+--   abs = error "abs: not implemented for Exp Word32" 
+--   fromInteger a = Literal (fromInteger a) 
   
 
 -- adding special shift operators for when both inputs are 
@@ -499,44 +546,46 @@
   testBit = error "testBit: is undefined for Exp Word32"
   popCount = error "popCoint: is undefined for Exp Word32"
 
-instance Real (Exp Word32) where 
-  toRational = error "toRational: not implemented for Exp Word32" 
+-- instance Real (Exp Word32) where 
+--   toRational = error "toRational: not implemented for Exp Word32" 
   
 
-instance Enum (Exp Word32) where
-  toEnum = error "toEnum: not implemented for Exp Word32" 
-  fromEnum = error "fromEnum: not implemented for Exp Word32" 
+-- instance Enum (Exp Word32) where
+--   toEnum = error "toEnum: not implemented for Exp Word32" 
+--   fromEnum = error "fromEnum: not implemented for Exp Word32" 
 
-instance Integral (Exp Word32) where
-  mod (Literal a) (Literal b) = Literal (a `mod` b) 
-  mod a b = BinOp Mod a b
-  div _ (Literal 0) = error "Division by zero in expression" 
-  div a b = BinOp Div a b
-  quotRem = error "quotRem: not implemented for Exp Word32" 
-  toInteger = error "toInteger: not implemented for Exp Word32"
+-- instance Integral (Exp Word32) where
+--   mod (Literal a) (Literal b) = Literal (a `mod` b) 
+--   mod a b = BinOp Mod a b
+--   div _ (Literal 0) = error "Division by zero in expression"
+--   div (Literal a) (Literal b) = Literal (a `div` b) 
+--   div a b = BinOp Div a b
+--   quotRem = error "quotRem: not implemented for Exp Word32" 
+--   toInteger = error "toInteger: not implemented for Exp Word32"
   
-instance Num (Exp Float) where
-  (+) a (Literal 0) = a
-  (+) (Literal 0) a = a
-  (+) (Literal a) (Literal b) = Literal (a + b)
-  (+) a b = BinOp Add a b
+-- instance Num (Exp Float) where
+--   (+) a (Literal 0) = a
+--   (+) (Literal 0) a = a
+--   (+) (Literal a) (Literal b) = Literal (a + b)
+--   (+) a b = BinOp Add a b
   
-  (-) a (Literal 0) = a
-  (-) (Literal a) (Literal b) = Literal (a - b)
-  (-) a b = BinOp Sub a b
+--   (-) a (Literal 0) = a
+--   (-) (Literal a) (Literal b) = Literal (a - b)
+--   (-) a b = BinOp Sub a b
   
-  (*) a (Literal 1) = a
-  (*) (Literal 1) a = a
-  (*) _ (Literal 0) = Literal 0
-  (*) (Literal 0) _ = Literal 0
-  (*) (Literal a) (Literal b) = Literal (a * b)
-  (*) a b = BinOp Mul a b
+--   (*) a (Literal 1) = a
+--   (*) (Literal 1) a = a
+--   (*) _ (Literal 0) = Literal 0
+--   (*) (Literal 0) _ = Literal 0
+--   (*) (Literal a) (Literal b) = Literal (a * b)
+--   (*) a b = BinOp Mul a b
   
-  signum = undefined
-  abs = undefined
-  fromInteger a = Literal (fromInteger a)
+--   signum = undefined
+--   abs = undefined
+--   fromInteger a = Literal (fromInteger a)
 
 instance Fractional (Exp Float) where
+  (/) (Literal a) (Literal b) = Literal (a/b)
   (/) a b = BinOp FDiv a b
   recip a = (Literal 1) / a
   fromRational a = Literal (fromRational a)
@@ -583,9 +632,6 @@
   -- other compilers have this).
   --(/) (Literal 1) (UnOp Sqrt b) = UnOp RSqrt b -- Optimisation.
 
-  
-  
-  
 ---------------------------------------------------------------------------
   
 infix 4 ==*, /=*, <*, >*, >=*, <=* 
@@ -605,6 +651,7 @@
 (&&*) a b = BinOp And a b 
 (||*) a b = BinOp Or a b 
 
+notE = UnOp Not
 ---------------------------------------------------------------------------
 -- Choice class
 ---------------------------------------------------------------------------
@@ -618,7 +665,13 @@
   
 instance (Choice a, Choice b) => Choice (a,b) where
   ifThenElse b (e1,e1') (e2,e2') = (ifThenElse b e1 e2,
-                                    ifThenElse b e1' e2') 
+                                    ifThenElse b e1' e2')
+
+instance (Choice a, Choice b, Choice c) => Choice (a,b,c) where
+  ifThenElse b (e1,e1',e1'') (e2,e2',e2'') = (ifThenElse b e1 e2,
+                                              ifThenElse b e1' e2',
+                                              ifThenElse b e1'' e2'')
+
   
 ---------------------------------------------------------------------------
 -- Print Expressions
@@ -654,7 +707,7 @@
 printOp GEq = " >= " 
 
 printOp And = " && "
-printOp Or  = " || " 
+printOp Or  = " || "
 
 printOp Min = " Min "
 printOp Max = " Max " 
@@ -667,149 +720,208 @@
 printOp BitwiseXor = " ^ " 
 printOp BitwiseNeg = " ~ "  
 
+
 ---------------------------------------------------------------------------
--- Turn expressions into backend-expressions
+-- Internal exp (not a GADT) 
 ---------------------------------------------------------------------------
-class ExpToCExp a where 
-  expToCExp :: Exp a -> CExpr 
 
+data IExp = IVar Name Type
+          | IBlockIdx  DimSpec
+          | IThreadIdx DimSpec
+          | IBlockDim  DimSpec
+          | IGridDim   DimSpec
 
-instance  ExpToCExp Bool where 
-  expToCExp (Literal True) = cLiteral (IntVal 1) CInt 
-  expToCExp (Literal False) = cLiteral (IntVal 0) CInt
-  expToCExp a = expToCExpGeneral a 
+          | IBool Bool 
+          | IInt8 Int8 | IInt16 Int16 | IInt32 Int32 | IInt64 Int64
+          | IWord8 Word8 | IWord16 Word16 | IWord32 Word32 | IWord64 Word64
+          | IFloat Float | IDouble Double
+                           
+          | IIndex (IExp,[IExp]) Type
+          | ICond IExp IExp IExp Type
+          | IBinOp IBinOp IExp IExp Type
+          | IUnOp  IUnOp  IExp Type
+          | IFunCall Name [IExp] Type
+          | ICast IExp Type
+          deriving (Eq, Ord, Show) 
+ 
 
-instance ExpToCExp Int where 
-  expToCExp (Literal a) = cLiteral (IntVal a) CInt
-  expToCExp a = expToCExpGeneral a  
+data IBinOp = IAdd | ISub | IMul | IDiv | IMod
+            | IEq | INotEq | ILt | IGt | IGEq | ILEq
+            | IAnd | IOr | IPow
+            | IBitwiseAnd | IBitwiseOr | IBitwiseXor
+            | IShiftL | IShiftR
+            deriving (Eq, Ord, Show) 
 
-instance ExpToCExp Int8 where 
-  expToCExp (Literal a) = cLiteral (Int8Val a) CInt8
-  expToCExp a = expToCExpGeneral a  
+data IUnOp = IBitwiseNeg | INot
+           deriving (Eq, Ord, Show)
 
-instance ExpToCExp Int16 where 
-  expToCExp (Literal a) = cLiteral (Int16Val a) CInt16
-  expToCExp a = expToCExpGeneral a  
 
-instance ExpToCExp Int32 where 
-  expToCExp (Literal a) = cLiteral (Int32Val a) CInt32
-  expToCExp a = expToCExpGeneral a  
 
-instance ExpToCExp Int64 where 
-  expToCExp (Literal a) = cLiteral (Int64Val a) CInt64
-  expToCExp a = expToCExpGeneral a  
+---------------------------------------------------------------------------
+-- Remove type info from operations
+---------------------------------------------------------------------------
 
-instance ExpToCExp Float where 
-  expToCExp (Literal a) = cLiteral (FloatVal a) CFloat
-  expToCExp a = expToCExpGeneral a 
+binOpToIBinOp :: Op t -> IBinOp
+binOpToIBinOp Add = IAdd
+binOpToIBinOp Sub = ISub
+binOpToIBinOp Mul = IMul
+binOpToIBinOp Div = IDiv 
+binOpToIBinOp FDiv = IDiv -- (???)  
+binOpToIBinOp Mod = IMod
 
-instance ExpToCExp Double where 
-  expToCExp (Literal a) = cLiteral (DoubleVal a) CDouble
-  expToCExp a = expToCExpGeneral a 
+binOpToIBinOp Eq  = IEq 
+binOpToIBinOp NotEq = INotEq
+binOpToIBinOp Lt  = ILt 
+binOpToIBinOp LEq = ILEq
+binOpToIBinOp Gt  = IGt 
+binOpToIBinOp GEq = IGEq 
 
-instance ExpToCExp Word where 
-  expToCExp (Literal a) = cLiteral (WordVal a) CWord
-  expToCExp a = expToCExpGeneral a 
+binOpToIBinOp And = IAnd 
+binOpToIBinOp Or  = IOr 
 
-instance ExpToCExp Word8 where 
-  expToCExp (Literal a) = cLiteral (Word8Val a) CWord8
-  expToCExp a = expToCExpGeneral a 
+binOpToIBinOp Pow = IPow
 
-instance ExpToCExp Word16 where 
-  expToCExp (Literal a) = cLiteral (Word16Val a) CWord16
-  expToCExp a = expToCExpGeneral a 
+binOpToIBinOp BitwiseAnd = IBitwiseAnd
+binOpToIBinOp BitwiseOr  = IBitwiseOr
+binOpToIBinOp BitwiseXor = IBitwiseXor
+binOpToIBinOp ShiftL     = IShiftL 
+binOpToIBinOp ShiftR     = IShiftR
 
-instance ExpToCExp Word32 where 
-  expToCExp (Literal a) = cLiteral (Word32Val a) CWord32
-  expToCExp a = expToCExpGeneral a 
+unOpToIUnOp :: Op t -> IUnOp
+unOpToIUnOp   BitwiseNeg = IBitwiseNeg
+unOpToIUnOp   Not = INot 
 
-instance ExpToCExp Word64 where 
-  expToCExp (Literal a) = cLiteral (Word64Val a) CWord64
-  expToCExp a = expToCExpGeneral a 
 
-  
-expToCExpGeneral :: ExpToCExp a  => Exp a -> CExpr
-expToCExpGeneral WarpSize      = cWarpSize 
-expToCExpGeneral (BlockIdx d)  = cBlockIdx d
-expToCExpGeneral (BlockDim d)  = cBlockDim d 
-expToCExpGeneral (ThreadIdx d) = cThreadIdx d
+---------------------------------------------------------------------------
+-- Turn Exp a to IExp with type information. 
+---------------------------------------------------------------------------
 
-expToCExpGeneral e@(Index (name,[])) = cVar name (typeToCType (typeOf e))
-expToCExpGeneral e@(Index (name,xs)) = cIndex (cVar name (CPointer (typeToCType (typeOf e))),map expToCExp xs) (typeToCType (typeOf e)) 
-expToCExpGeneral e@(If b e1 e2)      = cCond  (expToCExp b) (expToCExp e1) (expToCExp e2) (typeToCType (typeOf e))
 
-expToCExpGeneral (UnOp Word32ToInt32 e) = cCast (expToCExp e) CInt32
-expToCExpGeneral (UnOp Int32ToWord32 e) = cCast (expToCExp e) CWord32
+class ExpToIExp a where 
+  expToIExp :: Exp a -> IExp 
 
-expToCExpGeneral e@(BinOp Min e1 e2) = cFuncExpr "min" [expToCExp e1, expToCExp e2] (typeToCType (typeOf e)) 
-expToCExpGeneral e@(BinOp Max e1 e2) = cFuncExpr "max" [expToCExp e1, expToCExp e2] (typeToCType (typeOf e)) 
-expToCExpGeneral e@(BinOp op e1 e2)  = cBinOp (binOpToCBinOp op) (expToCExp e1) (expToCExp e2) (typeToCType (typeOf e)) 
 
-expToCExpGeneral (UnOp Exp e)        = cFuncExpr "exp" [expToCExp e] (typeToCType (typeOf e))
-expToCExpGeneral (UnOp Sqrt e)       = cFuncExpr "sqrt" [expToCExp e] (typeToCType (typeOf e))
-expToCExpGeneral (UnOp Log e)        = cFuncExpr "log" [expToCExp e] (typeToCType (typeOf e))
-expToCExpGeneral (UnOp Log2 e)       = cFuncExpr "log2" [expToCExp e] (typeToCType (typeOf e))
-expToCExpGeneral (UnOp Log10 e)      = cFuncExpr "log10" [expToCExp e] (typeToCType (typeOf e))
-  
--- Floating trig
-expToCExpGeneral (UnOp Sin e)        = cFuncExpr "sin" [expToCExp e] (typeToCType (typeOf e))
-expToCExpGeneral (UnOp Cos e)        = cFuncExpr "cos" [expToCExp e] (typeToCType (typeOf e))
-expToCExpGeneral (UnOp Tan e)        = cFuncExpr "tan" [expToCExp e] (typeToCType (typeOf e))
-expToCExpGeneral (UnOp ASin e)       = cFuncExpr "asin" [expToCExp e] (typeToCType (typeOf e))
-expToCExpGeneral (UnOp ACos e)       = cFuncExpr "acos" [expToCExp e] (typeToCType (typeOf e))
-expToCExpGeneral (UnOp ATan e)       = cFuncExpr "atan" [expToCExp e] (typeToCType (typeOf e))
-expToCExpGeneral (UnOp SinH e)       = cFuncExpr "sinh" [expToCExp e] (typeToCType (typeOf e))
-expToCExpGeneral (UnOp CosH e)       = cFuncExpr "cosh" [expToCExp e] (typeToCType (typeOf e))
-expToCExpGeneral (UnOp TanH e)       = cFuncExpr "tanh" [expToCExp e] (typeToCType (typeOf e))
-expToCExpGeneral (UnOp ASinH e)      = cFuncExpr "asinh" [expToCExp e] (typeToCType (typeOf e))
-expToCExpGeneral (UnOp ACosH e)      = cFuncExpr "acosh" [expToCExp e] (typeToCType (typeOf e))
-expToCExpGeneral (UnOp ATanH e)      = cFuncExpr "atanh" [expToCExp e] (typeToCType (typeOf e))
- 
-expToCExpGeneral e@(UnOp op e1)     = cUnOp  (unOpToCUnOp op) (expToCExp e1) (typeToCType (typeOf e))
- 
+instance  ExpToIExp Bool where 
+  expToIExp (Literal True) = IBool True 
+  expToIExp (Literal False) = IBool False 
+  expToIExp a = expToIExpGeneral a 
 
-typeToCType Bool = CInt 
-typeToCType Int  = CInt
-typeToCType Int8  = CInt8
-typeToCType Int16  = CInt16
-typeToCType Int32  = CInt32
-typeToCType Int64  = CInt64
-typeToCType Float = CFloat
-typeToCType Double = CDouble
-typeToCType Word8 = CWord8
-typeToCType Word16 = CWord16
-typeToCType Word32 = CWord32
-typeToCType Word64 = CWord64
-typeToCType (Pointer t) = CPointer (typeToCType t)
-typeToCType (Global t)  = CQualified CQualifyerGlobal (typeToCType t) 
-typeToCType (Local t)  = CQualified CQualifyerLocal (typeToCType t) 
+-- This is strange. 
+instance ExpToIExp Int where 
+  expToIExp (Literal a) = IInt32 (fromIntegral a)  
+  expToIExp a = expToIExpGeneral a  
 
--- maybe unnecessary
-binOpToCBinOp Add = CAdd
-binOpToCBinOp Sub = CSub
-binOpToCBinOp Mul = CMul
-binOpToCBinOp Div = CDiv 
-binOpToCBinOp FDiv = CDiv -- (???)  
-binOpToCBinOp Mod = CMod
+instance ExpToIExp Int8 where 
+  expToIExp (Literal a) = IInt8 a 
+  expToIExp a = expToIExpGeneral a  
 
-binOpToCBinOp Eq  = CEq 
-binOpToCBinOp NotEq = CNotEq
-binOpToCBinOp Lt  = CLt 
-binOpToCBinOp LEq = CLEq
-binOpToCBinOp Gt  = CGt 
-binOpToCBinOp GEq = CGEq 
+instance ExpToIExp Int16 where 
+  expToIExp (Literal a) = IInt16 a 
+  expToIExp a = expToIExpGeneral a  
 
-binOpToCBinOp And = CAnd 
-binOpToCBinOp Or  = COr 
+instance ExpToIExp Int32 where 
+  expToIExp (Literal a) = IInt32 a 
+  expToIExp a = expToIExpGeneral a  
 
-binOpToCBinOp Pow = CPow
+instance ExpToIExp Int64 where 
+  expToIExp (Literal a) = IInt64 a 
+  expToIExp a = expToIExpGeneral a  
 
-binOpToCBinOp BitwiseAnd = CBitwiseAnd
-binOpToCBinOp BitwiseOr  = CBitwiseOr
-binOpToCBinOp BitwiseXor = CBitwiseXor
-binOpToCBinOp ShiftL     = CShiftL 
-binOpToCBinOp ShiftR     = CShiftR
--- notice min and max is not here ! 
+instance ExpToIExp Float where 
+  expToIExp (Literal a) = IFloat a 
+  expToIExp a = expToIExpGeneral a 
 
-unOpToCUnOp   BitwiseNeg = CBitwiseNeg
+instance ExpToIExp Double where 
+  expToIExp (Literal a) = IDouble a 
+  expToIExp a = expToIExpGeneral a 
+
+-- This is strange. 
+instance ExpToIExp Word where 
+  expToIExp (Literal a) = IWord32 (fromIntegral a)
+  expToIExp a = expToIExpGeneral a 
+
+instance ExpToIExp Word8 where 
+  expToIExp (Literal a) = IWord8 a 
+  expToIExp a = expToIExpGeneral a 
+
+instance ExpToIExp Word16 where 
+  expToIExp (Literal a) = IWord16 a 
+  expToIExp a = expToIExpGeneral a 
+
+instance ExpToIExp Word32 where 
+  expToIExp (Literal a) = IWord32 a 
+  expToIExp a = expToIExpGeneral a 
+
+instance ExpToIExp Word64 where 
+  expToIExp (Literal a) = IWord64 a 
+  expToIExp a = expToIExpGeneral a
+
+
   
+expToIExpGeneral :: ExpToIExp a  => Exp a -> IExp
+expToIExpGeneral WarpSize      = IVar "warpsize" Word32 
+expToIExpGeneral (BlockIdx d)  = IBlockIdx d
+expToIExpGeneral (BlockDim d)  = IBlockDim d 
+expToIExpGeneral (ThreadIdx d) = IThreadIdx d
+
+expToIExpGeneral e@(Index (name,[])) = IVar name  (typeOf e)
+expToIExpGeneral e@(Index (name,xs))
+  = IIndex (IVar name (Pointer (typeOf e)),map expToIExp xs) (typeOf e) 
+expToIExpGeneral e@(If b e1 e2)
+  = ICond  (expToIExp b) (expToIExp e1) (expToIExp e2) (typeOf e)
+
+
+expToIExpGeneral (UnOp Word32ToInt32 e) = ICast (expToIExp e) Int32
+expToIExpGeneral (UnOp Int32ToWord32 e) = ICast (expToIExp e) Word32
+expToIExpGeneral (UnOp Word32ToFloat e) = ICast (expToIExp e) Float
+expToIExpGeneral (UnOp Word32ToWord8 e) = ICast (expToIExp e) Word8 
+
+expToIExpGeneral e@(BinOp Min e1 e2)
+  = IFunCall "min" [expToIExp e1, expToIExp e2] (typeOf e)
+    
+expToIExpGeneral e@(BinOp Max e1 e2)
+  = IFunCall "max" [expToIExp e1, expToIExp e2] (typeOf e)
+    
+expToIExpGeneral e@(BinOp op e1 e2)
+  = IBinOp (binOpToIBinOp op) (expToIExp e1) (expToIExp e2) (typeOf e)
+
+
+expToIExpGeneral (UnOp Exp e)        = IFunCall "exp" [expToIExp e] (typeOf e)
+expToIExpGeneral (UnOp Sqrt e)       = IFunCall "sqrt" [expToIExp e] (typeOf e)
+expToIExpGeneral (UnOp Log e)        = IFunCall "log" [expToIExp e]  (typeOf e)
+expToIExpGeneral (UnOp Log2 e)       = IFunCall "log2" [expToIExp e] (typeOf e)
+expToIExpGeneral (UnOp Log10 e)      = IFunCall "log10" [expToIExp e] (typeOf e)
+  
+-- Floating trig
+expToIExpGeneral (UnOp Sin e)        = IFunCall "sin" [expToIExp e] (typeOf e)
+expToIExpGeneral (UnOp Cos e)        = IFunCall "cos" [expToIExp e] (typeOf e)
+expToIExpGeneral (UnOp Tan e)        = IFunCall "tan" [expToIExp e] (typeOf e)
+expToIExpGeneral (UnOp ASin e)       = IFunCall "asin" [expToIExp e] (typeOf e)
+expToIExpGeneral (UnOp ACos e)       = IFunCall "acos" [expToIExp e] (typeOf e)
+expToIExpGeneral (UnOp ATan e)       = IFunCall "atan" [expToIExp e] (typeOf e)
+expToIExpGeneral (UnOp SinH e)       = IFunCall "sinh" [expToIExp e] (typeOf e)
+expToIExpGeneral (UnOp CosH e)       = IFunCall "cosh" [expToIExp e] (typeOf e)
+expToIExpGeneral (UnOp TanH e)       = IFunCall "tanh" [expToIExp e]  (typeOf e)
+expToIExpGeneral (UnOp ASinH e)      = IFunCall "asinh" [expToIExp e] (typeOf e)
+expToIExpGeneral (UnOp ACosH e)      = IFunCall "acosh" [expToIExp e] (typeOf e)
+expToIExpGeneral (UnOp ATanH e)      = IFunCall "atanh" [expToIExp e] (typeOf e)
+
+expToIExpGeneral e@(UnOp op e1) = IUnOp  (unOpToIUnOp op) (expToIExp e1) (typeOf e)
+ 
+
+---------------------------------------------------------------------------
+-- Collect arrays from an IExp 
+---------------------------------------------------------------------------
+
+collectArraysI :: String -> IExp -> [Name]
+collectArraysI pre e = go e
+  where
+    go (IVar name _) = if isPrefixOf pre name then [name] else []
+    go (IIndex (ne,es) _) = go ne ++ concatMap go es 
+    go (IBinOp _ e1 e2 _) = go e1 ++ go e2
+    go (IUnOp  _ e _) = go e
+    go (ICond b e1 e2 _) = go b ++ go e1 ++ go e2
+    go (IFunCall _ es _) = concatMap go es
+    go (ICast e _) = go e 
+    go _ = [] 
+
diff --git a/Obsidian/Force.hs b/Obsidian/Force.hs
--- a/Obsidian/Force.hs
+++ b/Obsidian/Force.hs
@@ -1,13 +1,21 @@
+{-# LANGUAGE ScopedTypeVariables #-} 
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE UndecidableInstances #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
 
-{-# LANGUAGE MultiParamTypeClasses,
-             FlexibleInstances,
-             ScopedTypeVariables,
-             TypeFamilies,
-             GADTs  #-}
 
+
 {- Joel Svensson 2012, 2013 
 
    Notes:
+   2014-03-28: Changed API.
+               Not using Obsidian.Mutable currently, it needs more work. 
+   2013-06-24: Changed code. uses Obsidian.Mutable now
+   2013-05-02: Removing things to do with forceG
+               Removed the extensions (no longer needed) 
+   2013-04-27: Something is broken. 
    2013-04-10: Looking at force and threads
    2013-01-27: globArrays nolonger exist
    2013-01-02: Added simple forceG for globArrays
@@ -15,11 +23,8 @@
 
 -}
 
---  write_ should be internal use only
-module Obsidian.Force (write,
-                       force,
-                       forceG
-                      ) where
+module Obsidian.Force (Forceable, force, forcePull, unsafeForce) where
+-- Write, force, forcePull, unsafeForce, unsafeWritePush) where 
 
 
 import Obsidian.Program
@@ -27,42 +32,97 @@
 import Obsidian.Array
 import Obsidian.Types
 import Obsidian.Globs
-import Obsidian.Memory
+import Obsidian.Memory 
+
 import Obsidian.Names
 
 import Data.Word
+
+import Control.Monad
 ---------------------------------------------------------------------------
--- Force local (requires static lengths!) 
+-- Force local (requires static lengths!)
 ---------------------------------------------------------------------------
 
-write :: forall p a. (Array p, Pushable p, MemoryOps a) => p Word32 a -> BProgram (Pull Word32 a)
-write arr = do 
-  -- snames <- names "arr" (undefined :: a)
-
-  -- Here I know that this pattern match will succeed
-  let n = len arr
+class Write t where
+  unsafeWritePush :: Storable a => Bool -> Push t Word32 a -> Program t (Pull Word32 a)
+  -- unsafeWritePull :: MemoryOps a => Pull Word32 a -> Program t (Pull Word32 a) 
   
-  -- allocateArray snames (undefined :: a) n
+instance Write Warp where
+  unsafeWritePush volatile p  =
+    do
+      let n = len p
+      names <- names "arr"
+      allocateVolatileArray names n
+     
+      p <: warpAssignArray names (variable "warpID") n 
+      return $ warpPullFrom names (variable "warpID") n
 
-  let (Push m p) = push Block arr
+instance Write Block where
+  unsafeWritePush volatile p =
+    do
+      let  n = len p
+      names <- names "arr"
+      if (volatile)
+        then allocateVolatileArray names n
+        else allocateArray names n
+             
+      p <: assignArray names 
+      return $ pullFrom names n
 
-  snames <- p (assignArrayN n) 
+instance Write Thread where
+  unsafeWritePush volatile p =
+    do
+      (snames :: Names a)  <- names "arr" 
+
+      -- Here I know that this pattern match will succeed
+      let n = len p
+    
+      allocateArray snames  n
+      p <: assignArray snames 
       
-  return $ pullFromS snames n
+      return $ pullFrom snames n
 
-  
-force :: forall p a. (Array p, Pushable p, MemoryOps a) =>  p Word32 a -> BProgram (Pull Word32 a)
+---------------------------------------------------------------------------
+-- Force functions 
+---------------------------------------------------------------------------
+
+-- | It is possible to force at level T if we can Write and Sync at that level. 
+class (Sync t, Write t) => Forceable t
+
+instance Forceable Thread
+instance Forceable Warp
+instance Forceable Block
+
+
+-- | force turns a @Push@ array to a @Program@ generating a @Pull@ array.
+--   The returned array represents reading from an array manifest in memory.
+force :: (Storable a, Forceable t)
+         => Push t Word32 a -> Program t (Pull Word32 a)      
 force arr = do
-  rval <- write arr
-  Sync
+  rval <- unsafeWritePush False arr
+  sync
   return rval
 
+-- | Make a @Pull@ array manifest in memory. 
+forcePull :: (Storable a, Forceable t, Pushable t) 
+             => Pull Word32 a -> Program t (Pull Word32 a)  
+forcePull = unsafeForce . push 
 
-forceG :: forall l a. GlobalMemoryOps a
-        => Push Grid l a
-        -> GProgram () 
-forceG (Push _ p)  =
-  do
-    output <- outputs (undefined :: a) 
-    p (\a e -> do {assignOut output a e; return (Single (Var,""))}) 
-    return ()
+-- | unsafeForce is dangerous on @Push@ arrays as it does not
+--   insert synchronization primitives. 
+unsafeForce :: (Storable a, Forceable t) =>
+         Push t  Word32 a -> Program t (Pull Word32 a)      
+unsafeForce arr = 
+  if (len arr <= 32)
+  then do
+    rval <- unsafeWritePush True arr
+    return rval
+  else do
+    rval <- unsafeWritePush False arr 
+    sync
+    return rval 
+
+
+
+
+
diff --git a/Obsidian/Library.hs b/Obsidian/Library.hs
--- a/Obsidian/Library.hs
+++ b/Obsidian/Library.hs
@@ -1,7 +1,8 @@
-{- Joel Svensson 2012, 2013 
+{- Joel Svensson 2012, 2013, 2014 
    Mary Sheeran  2012
 
    Notes:
+   2014-03-31: Merged Library and LibraryG 
    2013-01-24: GlobPull nolonger exists
                GlobPush is Push Grid
               
@@ -27,159 +28,226 @@
 import Obsidian.Force
 
 -- needed for threadsPerBlock analysis 
-import qualified Obsidian.CodeGen.Program as P 
+-- import qualified Obsidian.CodeGen.Program as P 
 
+import Control.Monad
+
 import Data.Bits 
 import Data.Word
 
-import Prelude hiding (splitAt,zipWith,replicate)
-
-
-
+import Prelude hiding (splitAt,zipWith,replicate,reverse)
 
 ---------------------------------------------------------------------------
 -- Reverse an array by indexing in it backwards
 ---------------------------------------------------------------------------
-  
-rev :: ASize l => Pull l a -> Pull l a 
-rev arr = mkPullArray n (\ix -> arr ! ((sizeConv m) - ix))  
+
+-- | Reverses a Pull array.
+reverse :: ASize l => Pull l a -> Pull l a 
+reverse arr = mkPull n (\ix -> arr ! ((sizeConv m) - ix))  
    where m = n-1
          n = len arr
          
 ---------------------------------------------------------------------------
 -- splitAt (name clashes with Prelude.splitAt)
 ---------------------------------------------------------------------------
+
+-- | Splits a Pull array at a given point. Performs no bounds checks. 
 splitAt :: (Integral i, ASize l) => i -> Pull l a -> (Pull l a, Pull l a) 
-splitAt n arr = (mkPullArray m (\ix -> arr ! ix), 
-                 mkPullArray  (len arr - m) (\ix -> arr ! (ix + pos)))
+splitAt n arr = (mkPull m (\ix -> arr ! ix), 
+                 mkPull  (len arr - m) (\ix -> arr ! (ix + pos)))
   where pos = fromIntegral n
         m   = fromIntegral n
 
-
+-- | Splits a Pull array in the middle. 
+halve :: ASize l => Pull l a -> (Pull l a, Pull l a) 
 halve arr = splitAt n2 arr
   where 
     n = len arr
     n2 = n `div` 2
 
+-- | Splits a Pull array into chunks of size n. Result is a Pull of Pull arrays.
+splitUp :: ASize l => Word32 -> Pull l a -> Pull l (SPull a)
+splitUp n arr {-(Pull m ixf)-} =
+  mkPull (len arr `div` fromIntegral n) $ \i ->
+    mkPull n $ \j -> arr ! (i * (sizeConv n) + j)                                               
+
+-- | Same as @splitUp@ but also performs a permutation of the elements. 
+coalesce :: ASize l
+         => Word32 -> Pull l a -> Pull l (Pull Word32 a)
+coalesce n arr =
+  mkPull s $ \i ->
+    mkPull n $ \j -> arr ! (i + (sizeConv s) * j)
+  where s = len arr `div` fromIntegral n
+        
 ---------------------------------------------------------------------------
--- replicate 
+-- elements at even indices to fst output, odd to snd.
 ---------------------------------------------------------------------------
-replicate n a = mkPullArray n (\ix -> a)
+-- | Split a Pull array into its even and odd indexed elements.
+evenOdds :: ASize l => Pull l a -> (Pull l a, Pull l a)
+evenOdds arr = (mkPull (n-n2) (\ix -> arr ! (2*ix)) ,
+                mkPull n2     (\ix -> arr ! (2*ix + 1)))
+  where
+    n  = len arr
+    n2 = div n 2
+-- | Extract the elements at even indices from a Pull array
+evens :: ASize l => Pull l a -> Pull l a
+evens = fst . evenOdds
 
-singleton a = replicate 1 a 
+-- | Extract the elements at odd indices from a Pull array 
+odds :: ASize l => Pull l a -> Pull l a
+odds  = snd . evenOdds
 
 ---------------------------------------------------------------------------
--- last
+-- everyNth 
 ---------------------------------------------------------------------------
+-- | Extract every nth element from a Pull array.
+everyNth :: ASize l => Word32 -> Word32 -> Pull l a -> Pull l a
+everyNth n m arr = mkPull n' $ \ix -> arr ! (ix * (fromIntegral n) + fromIntegral m)
+  where
+    n' = len arr `div` (fromIntegral n) 
+  
 
+---------------------------------------------------------------------------
+-- replicate 
+---------------------------------------------------------------------------
+-- | Generates a Pull array of length one, containing @a@. 
+singleton :: (Array a, ASize l) => e -> a l e
+singleton a = replicate 1 a
+
+-- | Generate a pull or push array usign a function from Index to element.
+generate :: (Functor (a s), Array a, ASize s)
+         => s -> (EWord32 -> b) -> a s b
+generate n f = fmap f (iota n) 
+---------------------------------------------------------------------------
+-- last and first 
+---------------------------------------------------------------------------
+-- | Extract last element from a Pull array.
+last :: ASize l => Pull l a -> a 
 last arr = arr ! fromIntegral ( len arr - 1)
 
+-- | Extract the first element from a Pull array.
+first :: ASize l => Pull l a -> a
+first arr = arr ! 0 
 
 ---------------------------------------------------------------------------
 -- Take and Drop (what about strange sizes ? fix) 
 ---------------------------------------------------------------------------
+-- | Take the first @n@ elements from a Pull array
 take :: ASize l => l -> Pull l a -> Pull l a
-take n arr = resize n arr
+take n arr = setSize n arr
 
+-- | Drop the first @n@ elements from a Pull array
 drop :: ASize l => l -> Pull l a -> Pull l a
-drop n arr = resize (len arr - n) $ ixMap (\ix -> ix + sizeConv n) arr
+drop n arr = setSize (len arr - n) $ ixMap (\ix -> ix + sizeConv n) arr
 
 ---------------------------------------------------------------------------
--- Shift arrays
+-- fold (sequential , unrolled)  
 ---------------------------------------------------------------------------
-shiftRight :: (ASize l, Choice a) => Word32 -> a -> Pull l a -> Pull l a
-shiftRight dist elt arr = resize (len arr)
-                          $ replicate (fromIntegral dist) elt `conc` arr
+-- | Fold a nonempty pull array using a given operator. The result a singleton array (push or pull). 
+fold1 :: Array a => (e -> e -> e) -> Pull Word32 e -> a Word32 e
+fold1  f arr = replicate 1
+               $ foldl1 f [arr ! (fromIntegral i) | i <- [0..(n-1)]]   
+  where n = len arr
 
--- TODO: incorrect! 
-shiftLeft :: (ASize l, Choice a) => Word32 -> a -> Pull l a -> Pull l a
-shiftLeft dist elt arr = mkPullArray (len arr)
-                         $ \ix -> (arr `conc`  replicate (fromIntegral dist) elt)
-                                  ! (ix + fromIntegral dist) 
-                         
+
 ---------------------------------------------------------------------------
--- elements at even indices to fst output, odd to snd.
+-- Shift arrays
 ---------------------------------------------------------------------------
-evenOdds :: ASize l => Pull l a -> (Pull l a, Pull l a)
-evenOdds arr = (mkPullArray (n-n2) (\ix -> arr ! (2*ix)) ,
-                mkPullArray n2     (\ix -> arr ! (2*ix + 1)))
-  where
-    n  = len arr
-    n2 = div n 2
-    
-evens, odds :: ASize l => Pull l a -> Pull l a
-evens = fst . evenOdds
-odds  = snd . evenOdds
+-- shiftRight :: (ASize l, Choice a) => Word32 -> a -> Pull l a -> Pull l a
+-- shiftRight dist elt arr = setSize (len arr)
+--                           $ replicate (fromIntegral dist) elt `conc` arr
 
+-- -- TODO: incorrect! 
+-- shiftLeft :: (ASize l, Choice a) => Word32 -> a -> Pull l a -> Pull l a
+-- shiftLeft dist elt arr = mkPull (len arr)
+--                          $ \ix -> (arr `conc`  replicate (fromIntegral dist) elt)
+--                                   ! (ix + fromIntegral dist) 
+                         
+
 -- opposite of evenOdds 
-shuffle :: ASize l => PT t -> Pull l a -> Pull l a -> Push t l a
-shuffle Block a1 a2 =
-  Push (len a1 + len a2) $
-    \ wf -> ForAll (sizeConv (len a1)) $
-            \ tid -> do
-              wf (a1 ! tid) (tid * 2) 
-              wf (a2 ! tid) (tid * 2 + 1) 
+--shuffle :: ASize l => PT t -> Pull l a -> Pull l a -> Push t l a
+--shuffle Block a1 a2 =
+--  Push (len a1 + len a2) $
+--    \ wf -> ForAll (sizeConv (len a1)) $
+--            \ tid -> do
+--              wf (a1 ! tid) (tid * 2) 
+--              wf (a2 ! tid) (tid * 2 + 1) 
 
 
 ---------------------------------------------------------------------------
 -- Concatenate the arrays
 ---------------------------------------------------------------------------
-conc :: (ASize l, Choice a) => Pull l a -> Pull l a -> Pull l a 
-conc a1 a2 = mkPullArray (n1+n2)
-               $ \ix -> ifThenElse (ix <* (fromIntegral n1)) 
-                       (a1 ! ix) 
-                       (a2 ! (ix - (fromIntegral n1)))
-  where 
-    n1 = len a1
-    n2 = len a2 
+-- -- | Concatenate two Pull arrays (Potentially inefficient).
+-- conc :: (ASize l, Choice a) => Pull l a -> Pull l a -> Pull l a 
+-- conc a1 a2 = mkPull (n1+n2)
+--                $ \ix -> ifThenElse (ix <* (sizeConv n1)) 
+--                        (a1 ! ix) 
+--                        (a2 ! (ix - (sizeConv n1)))
+--   where 
+--     n1 = len a1
+--     n2 = len a2 
 
     
 ---------------------------------------------------------------------------
 -- zipp unzipp
 ---------------------------------------------------------------------------
-unzipp :: ASize l =>  Pull l (a,b) -> (Pull l a, Pull l b)       
-unzipp arr = (mkPullArray (len arr) (\ix -> fst (arr ! ix)) ,
-              mkPullArray (len arr) (\ix -> snd (arr ! ix)) )
-              
-zipp :: ASize l => (Pull l a, Pull l b) -> Pull l (a, b)             
-zipp (arr1,arr2) =  Pull (min (len arr1) (len arr2))
-                      $ \ix -> (arr1 ! ix, arr2 ! ix) 
+-- | Unzip implemented on Pull arrays    
+unzip :: ASize l =>  Pull l (a,b) -> (Pull l a, Pull l b)       
+unzip arr = (mkPull (len arr) (\ix -> fst (arr ! ix)) ,
+             mkPull (len arr) (\ix -> snd (arr ! ix)) )
 
-unzipp3 :: ASize l => Pull l (a,b,c) 
+-- | Zip implemented on Pull arrays
+zip :: ASize l => Pull l a -> Pull l b -> Pull l (a, b)             
+zip arr1 arr2 = mkPull (min (len arr1) (len arr2))
+                     $ \ix -> (arr1 ! ix, arr2 ! ix) 
+
+-- | Unzip tripples. 
+unzip3 :: ASize l => Pull l (a,b,c) 
            -> (Pull l a, Pull l b, Pull l c)       
-unzipp3 arr = (fmap (\(x,_,_) -> x) arr,
+unzip3 arr = (fmap (\(x,_,_) -> x) arr,
                fmap (\(_,y,_) -> y) arr,
-               fmap (\(_,_,z) -> z)  arr) 
-
-
-zipp3 :: ASize l =>  (Pull l a, Pull l b, Pull l c) 
-         -> Pull l (a,b,c)             
-zipp3 (arr1,arr2,arr3) = 
-  mkPullArray (minimum [len arr1, len arr2, len arr3])
+               fmap (\(_,_,z) -> z)  arr)
+             
+-- | Zip three arrays
+zip3 :: ASize l
+        => Pull l a
+        -> Pull l b
+        -> Pull l c 
+        -> Pull l (a,b,c)             
+zip3 arr1 arr2 arr3 = 
+  mkPull (minimum [len arr1, len arr2, len arr3])
   (\ix -> (arr1 ! ix, arr2 ! ix, arr3 ! ix))
     
 
-
+-- | Perform elementwise operation.
 zipWith :: ASize l => (a -> b -> c) -> Pull l a -> Pull l b -> Pull l c
 zipWith op a1 a2 =  
-  mkPullArray (min (len a1) (len a2))
+  mkPull (min (len a1) (len a2))
   (\ix -> (a1 ! ix) `op` (a2 ! ix))
-                                      
+
+-- | Perform elementwise operation. 
+zipWith3 :: ASize l => (a -> b -> c-> d) -> Pull l a -> Pull l b -> Pull l c -> Pull l d
+zipWith3 f a1 a2 a3 =  
+  mkPull (minimum [len a1,len a2,len a3]) $ 
+    \ix -> f (a1 ! ix)  (a2 ! ix) (a3 ! ix)
+
+  
 ---------------------------------------------------------------------------
 -- pair 
 ---------------------------------------------------------------------------
+-- | Pair up consecutive elements in a Pull array.
 pair :: ASize l => Pull l a -> Pull l (a,a)
-pair (Pull n ixf) = 
-  mkPullArray n' (\ix -> (ixf (ix*2),ixf (ix*2+1))) 
+pair arr = 
+  mkPull n' (\ix -> (arr ! (ix*2),arr ! (ix*2+1))) 
   where 
-    n' = n `div` 2 
-
+    n' = len arr `div` 2 
 
 
+-- | Flatten a Pull array of pairs. 
 unpair :: ASize l => Choice a => Pull l (a,a) -> Pull l a
 unpair arr = 
     let n = len arr
-    in  mkPullArray (2*n) (\ix -> ifThenElse ((mod ix 2) ==* 0) 
+    in  mkPull (2*n) (\ix -> ifThenElse ((mod ix 2) ==* 0) 
                                   (fst (arr ! (ix `shiftR` 1)))
                                   (snd (arr ! (ix `shiftR` 1)))) 
 
@@ -188,21 +256,26 @@
 -- twoK (untested for proper functionality) 
 ---------------------------------------------------------------------------
 
+-- | Recursively split an array in the middle. Apply an array to array computation
+--   on each part. @binSplit 3@ divides the array into 8 pieces. 
+binSplit ::Int
+           -> (Pull Word32 a -> Pull Word32 b)
+           -> Pull Word32 a
+           -> Pull Word32 b 
 binSplit = twoK
 
 -- See if this should be specifically for Static size pull arrays
-twoK ::Int -> (Pull Word32 a -> Pull Word32 b) -> Pull Word32 a -> Pull Word32 b 
+twoK :: Int -> (Pull Word32 a -> Pull Word32 b) -> Pull Word32 a -> Pull Word32 b 
 twoK 0 f = f  -- divide 0 times and apply f
-twoK n f =  (\arr -> 
-              let arr' = mkPullArray lt (\i -> (f (mkPullArray  m (\j -> (arr ! (g i j)))) ! (h i))) 
+twoK n f = \arr -> 
+              let arr' = mkPull lt (\i -> (f (mkPull  m (\j -> (arr ! (g i j)))) ! (h i))) 
                   m    = (len arr `shiftR` n)   --pow of two           
                   g i j = i .&. (fromIntegral (complement (m-1))) .|. j  
                   h i   = i .&. (fromIntegral (nl2-1))   -- optimize 
 
-                  nl2   = (len (f (mkPullArray  m (\j -> arr ! variable "X"))))
+                  nl2   = len (f (mkPull  m (\j -> arr ! variable "X")))
                   lt    = nl2 `shiftL` n 
-              in arr')
- 
+              in arr'
 
 ---------------------------------------------------------------------------
 -- ***                          PUSHY LIBRARY                       *** ---
@@ -211,396 +284,204 @@
 ---------------------------------------------------------------------------
 -- Concatenate on Push arrays 
 ---------------------------------------------------------------------------
-concP :: (Array arr1, Array arr2, ASize l,
-          Pushable arr1, Pushable arr2)
-         => PT t -> arr1 l a -> arr2 l a -> Push t l a     
-concP pt arr1 arr2 = 
-  mkPushArray  (n1 + n2)
-  $ \wf ->
-  do
-    parr1 wf
-    parr2 $ \a i -> wf a (sizeConv n1 + i)
-  
-  where
-    n1 = len arr1 
-    n2 = len arr2 
-    (Push _ parr1) = push pt arr1
-    (Push _ parr2) = push pt arr2
 
-
--- More general versions of this can be imagined 
-mergeL :: (EWord32 -> a -> a -> a) -> Pull Word32 a -> Pull Word32 a -> Push Block Word32 a
-mergeL _ arr1 arr2 | len arr1 /= len arr2 = error "incorrect lengths" 
-mergeL f arr1 arr2 =
-  Push (len arr1) $ \wf ->
+-- | Concatenate two push arrays.
+concP :: ASize l
+         => Push t l a -> Push t l a -> Push t l a 
+concP p1 p2  =
+  mkPush (n1 + n2) $ \wf ->
   do
-    ForAll (sizeConv (len arr1)) $
-      \ tid -> wf (f tid (arr1 ! tid) (arr2 ! tid)) tid 
+    p1 <: wf
+    p2 <: \a i -> wf a (sizeConv n1 + i) 
+ where 
+   n1 = len p1
+   n2 = len p2 
 
 
-----------------------------------------------------------------------------
---
---unpairP :: Pushable arr => arr (a,a) -> Push a
---unpairP arr =
---  Push n $ \k -> pushf (everyOther k)
---  where
---    parr@(Push n pushf) = push arr
+-- Implement unpair on pusharrays,
+-- Impement for more tuples than pairs. 
+--unpair :: ASize l => Choice a => Pull l (a,a) -> Pull l a
+--unpair arr = 
+--    let n = len arr
+--    in  mkPull (2*n) (\ix -> ifThenElse ((mod ix 2) ==* 0) 
+--                                  (fst (arr ! (ix `shiftR` 1)))
+--                                  (snd (arr ! (ix `shiftR` 1)))) 
 
---everyOther :: (a -> Exp Word32 -> TProgram ()) 
---              -> (a,a) -> Exp Word32 -> TProgram ()
---everyOther wf (a,b) ix = wf a (ix * 2) >> wf b (ix * 2 + 1)  
 
----------------------------------------------------------------------------
--- zipP
----------------------------------------------------------------------------
---zipP :: Pushable arr  => arr a -> arr a -> Push a  
---zipP arr1 arr2 =
---  Push (n1+n2)
---  $ \func -> p1 (\a ix -> func a (2*ix)) >>
---             p2 (\a ix -> func a (2*ix + 1))
---         
---  where 
---    (Push n1 p1) = push arr1
---    (Push n2 p2) = push arr2
 
 
 ---------------------------------------------------------------------------
--- ***                          GLOBAL ARRAYS                        *** --
+-- load / Store 
 ---------------------------------------------------------------------------
-{- 
--- A very experimenental instance of mapG 
-mapG' :: (Pull a -> BProgram (Pull b))
-        -> Word32
-        -> GlobPull a
-        -> GlobPush b
-mapG' f n (GlobPull ixf)  =
-  GlobPush 
-        $ \wf ->
-          ForAllBlocks 
-           $ \bix ->
-             do -- BProgram do block 
-               let pully = Pull n (\ix -> ixf (bix * fromIntegral n + ix))
+load :: Word32 -> Pull Word32 a -> Push Block Word32 a 
+load n arr =
+  mkPush m (\wf ->
+  forAll (fromIntegral n') (\tid ->
+  do
+    seqFor (fromIntegral n) (\ix -> 
+      wf (arr ! (tid + (ix*fromIntegral n'))) (tid + (ix*fromIntegral n')))))
 
-               let res' = f pully 
-               res <- res' 
-               let numThreads = P.threadsPerBlock $ P.convPrg res'
-                   elemsPerThread = len res `div` numThreads
+  where
+    m = len arr
+    n' = m `div` n
 
-               if (numThreads == 0 || len res `mod` numThreads /= 0)
-                 then 
-                 ForAll (Just n) $ \ix -> wf (res ! ix) (bix * fromIntegral n + ix)
-                 else
-                 ForAll (Just n) $ \ix ->
-                 sequence_ [wf (res ! (ix + fromIntegral (numThreads * j)))
-                               (bix * fromIntegral n + (ix + fromIntegral (numThreads * j )))
-                            | j <- [0..elemsPerThread]]
+store :: Word32 -> SPull a -> SPush Block a 
+store = load 
 
 
--- Old fasioned mapG
-mapG :: (Pull a -> BProgram (Pull b))
-        -> Word32 -- BlockSize ! 
-        -> GlobPull a
-        -> GlobPush b
-mapG f n (GlobPull ixf)  =
-  GlobPush 
-        $ \wf ->
-          ForAllBlocks 
-           $ \bix ->
-             do -- BProgram do block 
-               let pully = Pull n (\ix -> ixf (bix * fromIntegral n + ix))
-               res <- f pully
-               ForAll (Just (len res)) $ \ix -> wf (res ! ix) (bix * fromIntegral n + ix)
+-- ########################################################################
+--
+-- Programming the Hierarchy
+--
+-- ########################################################################
 
 
--- I Think this one has more potential for generalisations. 
-mapD :: (Pull a -> BProgram b)
-        -> Word32
-        -> GlobPull a
-        -> DistPull b
-mapD f n (GlobPull ixf) =
-  DistPull $ \bix ->
-    do
-      let pully = Pull n (\ix -> ixf (bix * fromIntegral n + ix))
-      f pully 
 
-mapDist :: (Pull a -> BProgram b)
-           -> Word32
-           -> GlobPull a
-           -> (Exp Word32 -> BProgram b)
-mapDist f n (GlobPull ixf) bix =
-  let pully = Pull n (\ix -> ixf (bix * fromIntegral n + ix))
-  in  f pully 
-
--} 
-
 ---------------------------------------------------------------------------
--- See if having an Assignable/Allocable class is enough to generalise
--- the code below to pairs etc of Exp
+-- Parallel concatMap  
 ---------------------------------------------------------------------------
--- Experimental (Improve upon this if it works)
--- I think this code looks quite horrific right now.
--- The potentially unnecessary assignments are pretty bad.
 
-{- 
-    
-class ToGProgram a where
-  type Global a
-  toGProgram :: (Exp Word32 -> BProgram a) -> GProgram (Global a)
-
-instance Scalar a => ToGProgram (Pull (Exp a)) where
-  type Global (Pull (Exp a)) = GlobPush (Exp a)
-  toGProgram f =
-    do      
-      let (pulla,_) = runPrg 0 $ f (BlockIdx X)
-      let n = len pulla 
-
-      shared <- uniqueSM
-
-      ForAllBlocks $ \bix ->
-        do
-          res <- f bix -- compute res.
+pConcatMap f = pConcat . fmap f
+pUnCoalesceMap f = pUnCoalesce . fmap f
+pConcatMapJoin f = pConcat . fmap (runPush.f)
+pUnCoalesceMapJoin f = pUnCoalesce . fmap (runPush.f)
+pCoalesceMap n f = pUnCoalesce . fmap f . coalesce n
+pSplitMap n f = pConcat . fmap f . splitUp n
 
-          -- Sync
-  
-          Allocate shared (n * fromIntegral (sizeOf (undefined :: Exp a)))
-                          (Pointer (typeOf (undefined :: Exp a)))
+---------------------------------------------------------------------------
+-- Step into the Hierarchy by distributing a
+-- Thread program parameterized on a threadId over the threads
+-- at a specific level in the Hierarchy. 
+---------------------------------------------------------------------------
 
-          ForAll (Just n) $ \tix ->
-            -- potentially unnessecary assignment...
-            -- if the last thing the local computation does is force. 
-            Assign shared tix (res ! tix)
+-- | A way to enter into the hierarchy
+-- A bunch of Thread computations, spread across the threads of either
+-- a Warp, block or grid. (or performed sequentially in a single thread) 
+tConcat :: ASize l
+           => Pull l (SPush Thread b)
+           -> Push t l b  
+tConcat arr =
+  mkPush (n * fromIntegral s) $ \wf -> do
+    forAll (sizeConv n) $ \tid -> 
+       let wf' a ix = wf a (tid * sizeConv s + ix)
+           p = arr ! tid -- f tid 
+       in p <: wf'
+  where
+    n = len arr
+    s  = len (arr ! 0) --(f (variable "tid")) -- arr
 
+-- | Variant of @tConcat@. 
+tDistribute :: ASize l
+               => l
+               -> (EWord32 -> SPush Thread b)
+               -> Push t l b
+tDistribute n f = tConcat (mkPull n f) 
 
-          Sync
-          
-      return $
-        GlobPush $ \wf ->
-        do
-          ForAllBlocks $ \bix-> 
-            ForAll (Just n) $ \tix ->
-              wf (index shared tix)
-                 (bix * fromIntegral n + tix)
-        
-                
-                        
       
-instance (Scalar a, Scalar b) => ToGProgram (Pull (Exp a), Pull (Exp b)) where
-  type Global (Pull (Exp a),Pull (Exp b)) = (GlobPush (Exp a), GlobPush (Exp b))
-  toGProgram f = 
-    do      
-      let ((pulla,pullb),_) = runPrg 0 $ f (BlockIdx X)
-      let n1 = len pulla 
-      let n2 = len pullb  
-
-
-      shared1 <- uniqueSM
-      shared2 <- uniqueSM 
-
-      ForAllBlocks $ \bix ->
-        do
-          -- This is the heart of the matter. I want the f bix Program 
-          --  to only occur once in the generated complete Program.
-          (res1,res2) <- f bix -- compute results.
-
-          --  Sync
-      
-          ------------------------------------------------------------------
-          Allocate shared1 (n1 * fromIntegral (sizeOf (undefined :: Exp a)))
-                           (Pointer (typeOf (undefined :: Exp a)))
-                         
-          ForAll (Just n1) $ \tix ->
-            -- potentially unnessecary assignment...
-            -- if the last thing the local computation does is force. 
-            Assign shared1 tix (res1 ! tix)
-
-          ------------------------------------------------------------------
-          Allocate shared2 (n2 * fromIntegral (sizeOf (undefined :: Exp b)))
-                   (Pointer (typeOf (undefined :: Exp b)))
-  
-          ForAll (Just n2) $ \tix ->
-            -- potentially unnessecary assignment...
-            -- if the last thing the local computation does is force. 
-            Assign shared2 tix (res2 ! tix)
-
-          Sync
-            
-      let gp1 = 
-            GlobPush $ \wf ->
-              do
-                ForAllBlocks $ \bix-> 
-                  ForAll (Just n1) $ \tix ->
-                  wf (index shared1 tix)
-                  (bix * fromIntegral n1 + tix) 
-
-      let gp2 =
-            GlobPush $ \wf ->
-            do
-              ForAllBlocks $ \bix-> 
-                ForAll (Just n2) $ \tix ->
-                wf (index shared2 tix)
-                   (bix * fromIntegral n2 + tix)
-
+-- | Distribute work across the parallel resources at a given level of the GPU hiearchy
+pConcat :: ASize l => Pull l (SPush t a) -> Push (Step t) l a
+pConcat arr =
+  mkPush (n * fromIntegral rn) $ \wf ->
+    distrPar (sizeConv n) $ \bix ->
+      let p = arr ! bix 
+          wf' a ix = wf a (bix * sizeConv rn + ix) 
           
-      return (gp1,gp2)
--} 
---------------------------------------------------------------------------- 
-        
-
--- The Problem is that I cannot share computations that
--- take place on the gridlevel (I think). 
--- Experiment
-{- 
-mapE :: forall a b . (Scalar a, Scalar b)
-        => (Pull (Exp a) -> BProgram (Pull (Exp b)))
-        -> Word32
-        -> GlobPull (Exp a)
-        -> GProgram (Pull (Exp b)) 
-mapE f n (GlobPull ixf) =
-  do
-
-
-    shared <- uniqueSM 
-    -- Allocate bytes in all shared memories and obtain a name
-    
-    --return undefined 
-    ForAllBlocks $ \bix ->
-      do 
-        let pully = Pull n (\ix -> ixf (bix * fromIntegral n + ix))
-        res <- f pully
-
-        Allocate shared (n * fromIntegral (sizeOf (undefined :: Exp b)))
-                       (Pointer (typeOf (undefined :: Exp b)) )
-        ForAll (Just n) $ \ tid ->        
-          do 
-            Assign shared tid (res ! tid)
-        
-    return $ Pull 32 $ \ix -> (index shared (ix `mod` (fromIntegral n)))
--} 
-  
- -- Assume Pull a here is one the special distributed pulls from above
-{- 
-pushBlocks :: Pull a -> GlobPush a
-pushBlocks (Pull n ixf) =
-  GlobPush $ \wf ->
-     ForAllBlocks $ \bix ->
-       ForAll (Just n) $ \tix -> wf (ixf tix) (bix * fromIntegral n + tix) 
-
+      in p <: wf'
+  where
+    n  = len arr
+    rn = len (arr ! 0) -- All arrays are same length
 
-experiment :: GlobPull (Exp Int) -> GProgram (Pull (Exp Int), Pull (Exp Int))
-experiment input =
-  do 
-    arr <- mapE force 32 input
-    return (arr, singleton (arr ! 31))
+-- | Distribute work across the parallel resources at a given level of the GPU hierarchy
+pDistribute :: ASize l => l -> (EWord32 -> SPush t a) -> Push (Step t) l a
+pDistribute n f = pConcat (mkPull n f) 
 
-experiment2 :: GlobPull (Exp Int) -> GProgram (GlobPush (Exp Int), GlobPush (Exp Int))
-experiment2 input =
+-- | Sequential concatenation of a Pull of Push.
+sConcat :: ASize l => Pull l (SPush t a) -> Push t l a
+sConcat arr =
+  mkPush (n * fromIntegral rn) $ \wf ->
   do
-    (arr1,arr2) <- experiment input
-    return (pushBlocks arr1, pushBlocks arr2)
--} 
-                                 
-{-
-  I think something like mapE above is needed to expressed shared computations.
-  Bad things about mapE is its very specific type and
-  that it allocates a new shared memory array and potentially performs a
-  completely unnecessary copy from the old shared memory array.
+    seqFor (sizeConv n) $ \bix ->
+      let p = arr ! bix -- (Push _ p) = arr ! bix
+          wf' a ix = wf a (bix * sizeConv rn + ix)              
+      in p <: wf'
+  where 
+    n  = len arr
+    rn = len $ arr ! 0
 
-  To generalise mapE typeclasses are probably needed.
+-- | Variant of sConcat.
+sDistribute :: ASize l => l -> (EWord32 -> SPush t a) -> Push t l a
+sDistribute n f = sConcat (mkPull n f) 
 
-  A way to skip the unnecessary copy would be if it was possible
-  to hand the name over to the local computation. "compute this and
-  store the result here".
-  This sounds like some notion of a mutable array...
-  What would it mean if we could express such local computations ?
+-- pUnCoalesce adapted from Niklas branch.
+-- | Combines work that was distributed in a Coalesced way.
+-- | Applies a permutation on stores.
+pUnCoalesce :: ASize l => Pull l (SPush t a) -> Push (Step t) l a
+pUnCoalesce arr =
+  mkPush (n * fromIntegral rn) $ \wf ->
+  distrPar (sizeConv n) $ \bix ->
+    let p = arr ! bix
+        wf' a ix = wf a (bix * sizeConv rn + ix)
+    in p <: (g wf')
+  where
+    n  = len arr
+    rn = len $ arr ! 0
+    s  = sizeConv rn 
+    g wf a i = wf a (i `div` s + (i`mod`s)*(sizeConv n))
 
-  sklansky :: Int -> (Exp a -> Exp a -> Exp a) -> Pull (Exp a) -> Name -> BProgram (Pull (Exp a) 
-  sklansky 0 op arr res =
-    forAllN (len arr) $ \ix ->
-      Assign res ix (arr ! ix)
-  sklansky n op arr res =
-    let arr1 = twoK (n-1) (fan op) arr
-    arr2 <- force arr1
-    sklansky (n-1) op arr2 
+---------------------------------------------------------------------------
+-- RunPush 
+---------------------------------------------------------------------------
 
-  Not pure... But then, who says our dsl must be ? (isn't that part of the beauty
-  of Haskell and embedded languages?). We are already in a wierd Program Monad
-  and allow really dangerous push arrays... 
+-- | Fuses the program that computes a Push array into the Push array. 
+runPush :: Program t (Push t s a) -> Push t s a
+runPush prg =
+  mkPush n $ \wf -> do
+    parr <- prg
+    parr <: wf
+    -- It is a bit scary that I need to "evaluate" programs here. 
+  where n = len $ fst $ runPrg 0 prg
 
-  Going even further then maybe force should take a mutable array and a
-  push/pull array and they push the elements into the mutable array..
-  This means that creation of mutable arrays must also be in the hands of
-  the programmer.
-  If array creation and force targets are in the hand of the programmer then
-  in-placeness is also in the hands of the programmer + about a trillion
-  new ways to shoot ones foot off. 
+-- | Lifts @runPush@ to one input functions.
+runPush1 :: (a -> Program t (Push t s b)) -> a -> Push t s b
+runPush1 f a = runPush (f a)
 
-  I imagine that a mutable array could be simply a:
-  data MArray a = MArray Word32 Name -- length and identifier.
--} 
+-- | Lifts @runPush@ to two input functions.
+runPush2 :: (a -> b -> Program t (Push t s c)) -> a -> b -> Push t s c
+runPush2 f a b = runPush (f a b)  
 
+-- | Converts a program computing a pull Array to a Push array
+runPull :: (Pushable t, ASize s) => Program t (Pull s a) -> Push t s a
+runPull = runPush . liftM push 
 
+-- | Lifts @runPull@ to one input functions.
+runPull1 :: (Pushable t, ASize s) => (a -> Program t (Pull s b)) -> a -> Push t s b
+runPull1 f a = runPull (f a)
 
--- mapG2 is really hard to get right. 
-{- 
-mapG2 :: (Pull a -> BProgram (Pull b, Pull c))
-         -> Word32
-         -> GlobPull a
-         -> (GlobPush b, GlobPush c) -- hard to get this right!
-                                     -- without repeating computations. 
-mapG2 f n (GlobPull ixf) =
-  (GlobPush
-   $ \wf -> ForAllBlock
-            $ \bix ->
-            do
-              let pully = Pull n (\ixf (bix * fromIntegral n + ix))
-                  res = f pully 
-   , GlobPush 
--} 
+-- | Lifts @runPull@ to two input functions.
+runPull2 :: (Pushable t, ASize s) => (a -> b -> Program t (Pull s c)) -> a -> b -> Push t s c
+runPull2 f a b = runPull (f a b)
 
 ---------------------------------------------------------------------------
--- From Distributed array of blocks to a global push array
+-- 
 ---------------------------------------------------------------------------
---toGlobPush :: Distrib (BProgram (Pull a))
---               -> GlobPush a               
---toGlobPush inp@(Distrib bixf) =
---  GlobPush bs $
---    \wf -> ForAllBlocks 
---           $ \bix ->
---           do -- BProgram do block 
---             arr <- bixf bix 
---             ForAll bs $ \ix -> wf (arr ! ix) bix ix 
---  where
---    -- Is this Ok?! 
---    bs = len $ fst $ runPrg 0 $ bixf 0
+pushPrg :: Program t a -> SPush t a
+pushPrg = singletonPush
 
----------------------------------------------------------------------------
--- Create a global array that pushes to global
--- memory N elements per thread. 
---------------------------------------------------------------------------- 
---toGlobPushN :: Word32
---                -> Distrib (BProgram (Pull a))
---                -> GlobPush a
---toGlobPushN n dist =
---  GlobPush bs $ 
---  \wf -> ForAllBlocks 
---         $ \bix ->
---         do -- BProgram do block
---             arr <- getBlock dist bix
---             ForAll (bs `div` n) $
---               \ix ->
---                    sequence_ 
---                    -- correct indexing ? 
---                    [wf (arr ! (ix * n' + i')) bix (ix * n' + i')
---                   | i <- [0..n-1]
---                    , let n' = fromIntegral n
---                    , let i' = fromIntegral i]
---           
---                  
--- where
---    bs = len $ fst $ runPrg 0 $ getBlock dist 0 
---    -- nb = numBlocks dist
 
-    
 ---------------------------------------------------------------------------
---
+-- Singleton push arrays 
 ---------------------------------------------------------------------------
+
+-- Danger! use only with Scalar a's 
+-- -- | Create a singleton Push array.
+--singletonPush :: a -> SPush t a
+--singletonPush = singletonPushP . return 
+
+-- | Monadic version of @singleton@.
+singletonPush :: Program t a -> SPush t a
+singletonPush prg =
+  mkPush 1 $ \wf -> do
+    a <- prg
+    forAll 1 $ \ix -> 
+      wf a 0
+
diff --git a/Obsidian/LibraryG.hs b/Obsidian/LibraryG.hs
deleted file mode 100644
--- a/Obsidian/LibraryG.hs
+++ /dev/null
@@ -1,131 +0,0 @@
-
-{-# LANGUAGE ScopedTypeVariables #-}
-module Obsidian.LibraryG where
-
-import Obsidian.Array
-import Obsidian.Program
-import Obsidian.Exp
-import Obsidian.Memory
-
-import Data.Word
-
-
----------------------------------------------------------------------------
--- 
---------------------------------------------------------------------------- 
-mapG :: ASize l => (SPull a -> BProgram (SPull b))
-        -> Pull l (SPull a)
-        -> Push Grid l b
-mapG kern as =
-  Push (blocks * fromIntegral rn) $
-  \wf ->
-    do
-      forAllBlocks (sizeConv blocks) $ \bix -> do
-        res <- kern (as ! bix)
-        let (Push _ p) = push Block res
-            wf' a ix = wf a (bix * sizeConv rn + ix)
-        p wf'
-  where
-    blocks = len as
-    -- TODO: ensure this is not insane (runPrg functionality) 
-    rn = len $ fst $ runPrg 0 (kern (as ! 0))
-    n = len (as ! 0)
-
-
-mapT :: (SPull a -> TProgram (SPull b))
-        -> SPull (SPull a)
-        -> SPush Block b
-mapT threadf as =
-  Push (n * m) $
-  \wf ->
-    do
-      forAll (sizeConv n) $ \tix -> do
-        res <- threadf (as ! tix)
-        let (Push _ p) = push Thread res
-            wf' a ix = wf a (tix * sizeConv m + ix)
-        p wf'      
-
-  where
-    n = len as
-    m = len (as ! 0) 
----------------------------------------------------------------------------
--- 
---------------------------------------------------------------------------- 
-zipWithG :: ASize l => (SPull a -> SPull b -> BProgram (SPull c))
-           -> Pull l (SPull a)
-           -> Pull l (SPull b)
-           -> Push Grid l c
-zipWithG kern as bs =
-  Push (blocks * fromIntegral rn) $
-  \wf ->
-    do
-      forAllBlocks (sizeConv blocks) $ \bix -> do
-        res <- kern (as ! bix) (bs ! bix)
-        let (Push _ p) = push Block res
-            wf' a ix = wf a (bix * sizeConv rn + ix)
-        p wf'
-  where
-    -- Is this ok?! (does it break?) 
-    rn = len $ fst $ runPrg 0 (kern (as ! 0) (bs ! 0))
-    n = min m k
-    -- assume uniformity
-    m = len (as ! 0)
-    k = len (bs ! 0)
-    blocks = min (len as) (len bs) 
-     
-
-zipWithT :: (SPull a -> SPull b -> TProgram (SPull c))
-        -> SPull (SPull a)
-        -> SPull (SPull b) 
-        -> SPush Block c
-zipWithT threadf as bs =
-  Push (threads * rn) $
-  \wf ->
-    do
-      forAll (sizeConv threads) $ \tix -> do
-        res <- threadf (as ! tix) (bs ! tix) 
-        let (Push _ p) = push Thread res
-            wf' a ix = wf a (tix * sizeConv n + ix)
-        p wf'      
-
-  where
-
-    -- Is this ok?! (does it break?) 
-    rn = len $ fst $ runPrg 0 (threadf (as ! 0) (bs ! 0))
-    n = min m k 
-
-    m  = len (as ! 0)
-    k  = len (bs ! 0)
-    threads = min (len as) (len bs) 
-
----------------------------------------------------------------------------
--- Experimental
----------------------------------------------------------------------------
-zipWithG' :: forall a b c l. (ASize l, MemoryOps c)
-             => (SPull a -> SPull b -> BProgram (SPull c))
-             -> Pull l (SPull a)
-             -> Pull l (SPull b)
-             -> GProgram (Push Grid l c)
-zipWithG' kern as bs =
-  do
-    snames <- forAllBlocks (sizeConv n) $ \bix ->
-      do
-        res <- kern (as ! bix) (bs ! bix)
-        let (Push _ p) = push Block res
-        p (assignArrayN n)
-    let pully = Pull blocks $ \bix -> (pullFromS snames n :: Pull Word32 c)
-        
-    return $ Push (blocks * fromIntegral n) $
-      \wf ->
-      do
-        forAllBlocks (sizeConv blocks) $ \bix ->
-          forAll (sizeConv n) $ \tix ->
-            do
-              wf ((pully ! bix) ! tix) (bix * (sizeConv n) + tix)
-      
-  where
-    n = min m k 
-    -- Assume uniformity
-    m = len (as ! 0)
-    k = len (bs ! 0)
-    blocks = (min (len as) (len bs))
diff --git a/Obsidian/Memory.hs b/Obsidian/Memory.hs
--- a/Obsidian/Memory.hs
+++ b/Obsidian/Memory.hs
@@ -6,9 +6,11 @@
    This Module became quite messy.
    TODO: CLEAN IT UP! 
 
+   notes: 2013-05-02: Cleaned out inspect. 
+
 -} 
 
-module Obsidian.Memory (MemoryOps(..),GlobalMemoryOps(..),assignArrayN)  where
+module Obsidian.Memory (Storable(..))  where
 
 
 import Obsidian.Program
@@ -20,165 +22,182 @@
 
 import Data.Word
 
-class Inspect a where
-  inspect :: a -> Tree (Either (Kind,Name) Type)
-
-instance Scalar a => Inspect (Exp a) where
-  inspect (Index (name,[])) = Single (Left (Var,name))
-  inspect (Index (name,[ix])) | isid ix =  Single (Left (Arr,name))
-  inspect _ = Single $ Right (typeOf (undefined :: Exp a)) 
-
-instance (Inspect a, Inspect b) => Inspect (a,b) where
-  inspect (a,b) = Tuple [inspect a, inspect b]
-
-instance (Inspect a, Inspect b, Inspect c) => Inspect (a,b,c) where
-  inspect (a,b,c) = Tuple [inspect a, inspect b, inspect c]
-
-isid (ThreadIdx X) = True
-isid _ = False
+-- class MemoryOps a => Storable a 
 
 ---------------------------------------------------------------------------
 -- Local Memory
 ---------------------------------------------------------------------------
-class Inspect a => MemoryOps a where
-  names          :: String -> a -> Program t Names
-  allocateArray  :: Names -> a -> Word32 -> Program t ()
-  allocateScalar :: Names -> a -> Program t () 
-  assignArray    :: Names -> a -> Exp Word32 -> TProgram ()
-  assignArrayS   :: Tree (Either (Kind,Name) (Kind,Name))
-                    -> a
-                    -> Exp Word32
-                    -> TProgram (Tree (Kind,Name)) 
-  assignScalar   :: Names -> a -> TProgram () 
-  pullFrom       :: Names -> Word32 -> Pull Word32 a
-  readFrom       :: Names -> a
+class Storable a where
+  -- | Obtain new names for variables / arrays 
+  names          :: String -> Program t (Names a) 
 
-  pullFromS      :: Tree (Kind,Name) -> Word32 -> Pull Word32 a
+  -- Array operations 
+  assignArray    :: Names a -> a -> Exp Word32 -> Program Thread ()
+  allocateArray  :: Names a -> Word32 -> Program t ()
+  pullFrom       :: Names a -> Word32 -> Pull Word32 a
 
----------------------------------------------------------------------------
--- Derived
----------------------------------------------------------------------------
-assignArrayN :: MemoryOps a =>
-                Word32 -> a -> Exp Word32 -> TProgram (Tree (Kind,Name))
-assignArrayN n a ix  =
-  do
-    names <- allocateNeeded n insp
-    assignArrayS names a ix
-    -- return names
-  where
-    insp = inspect a
+  
 
+  -- Scalar operations 
+  assignScalar   :: Names a -> a -> Program Thread ()
+  allocateScalar :: Names a ->  Program t () 
+  readFrom       :: Names a -> a
+  
+  
+  -- Warp level operations   
+  warpAssignArray   :: Names a
+                      -> EWord32
+                      -> Word32
+                      -> a
+                      -> EWord32
+                      -> Program Thread ()
+  warpPullFrom      :: Names a -> EWord32 -> Word32 -> Pull Word32 a
+  
+  -- Extra
+  allocateVolatileArray :: Names a -> Word32 -> Program t ()
 
-allocateNeeded :: Word32 -> Tree (Either (Kind,Name) Type) -> TProgram (Tree (Either (Kind,Name) (Kind,Name)))
-allocateNeeded n None = return None
-allocateNeeded n a@(Single (Left (k,nom))) = return $ Single (Left (k,nom))
-allocateNeeded n (Single (Right t)) =
-  do
-    name <- uniqueNamed "arr"
-    Allocate name (n * typeSize t) (Pointer t)
-    return $ Single (Right (Arr,name))
-allocateNeeded n (Tuple xs) =
-  do
-    ns <- mapM (allocateNeeded n) xs
-    return $ Tuple ns 
-    
 
+
+
 ---------------------------------------------------------------------------
 -- Instances
 ---------------------------------------------------------------------------
-instance Scalar a => MemoryOps (Exp a) where
-  names pre a = do {i <- uniqueNamed pre; return (Single i)}
-  allocateArray (Single name) a n = 
-    Allocate name (n * fromIntegral (sizeOf a))
-                  (Pointer (typeOf a))
-  allocateScalar (Single name) a =
-    Declare name (typeOf a) 
+instance Scalar a => Storable (Exp a) where
+
+  -- Names 
+  names pre = do {i <- uniqueNamed pre; return (Single i)}
+
+  --Array ops 
+  allocateArray (Single name) n = 
+    Allocate name (n * fromIntegral (sizeOf (undefined :: Exp a)))
+                  (Pointer (typeOf (undefined :: Exp a)))
   assignArray  (Single name) a ix = Assign name [ix] a
-  assignArrayS (Single (Left (k,name))) a ix = return $ Single (k,name)
-  assignArrayS (Single (Right (k,name))) a ix =
-    do 
-      Assign name [ix] a -- ? 
-      return $ Single (k,name)
-    
+  pullFrom (Single name) n = mkPull n (\i -> index name i)
   
+  -- Scalar ops 
+  allocateScalar (Single name) =
+    Declare name (typeOf (undefined :: Exp a)) 
+  assignScalar (Single name) a    = Assign name [] a
+  readFrom  (Single name) = variable name
 
-  
+  -- Warp ops   
+  warpAssignArray (Single name) warpID step a ix =
+    Assign name [warpID * fromIntegral step + ix] a 
 
-  assignScalar (Single name) a    = Assign name [] a  
-  pullFrom (Single name) n = Pull n (\i -> index name i) 
-  readFrom (Single name) = variable name
+  warpPullFrom (Single name) warpID n
+    = mkPull n (\i -> index name (warpID * fromIntegral n + i))
 
-  pullFromS (Single (Var,name)) n = Pull n $ \_ -> variable name
-  pullFromS (Single (Arr,name)) n = Pull n (\i -> index name i) 
+  -- Extra 
+  allocateVolatileArray (Single name) n = 
+    Allocate name (n * fromIntegral (sizeOf (undefined :: Exp a)))
+                  (Volatile (Pointer (typeOf (undefined :: Exp a))))
+  
 
-instance (MemoryOps a, MemoryOps b) => MemoryOps (a, b) where
-  names pre (a,b) =
-    do
-      a' <- names pre a
-      b' <- names pre b
-      return $ Tuple [a', b']
-  allocateArray (Tuple [ns1,ns2]) (a,b) n =
-    do 
-      allocateArray ns1 a n
-      allocateArray ns2 b n
-  allocateScalar (Tuple [ns1,ns2]) (a,b) =
-    do
-      allocateScalar ns1 a
-      allocateScalar ns2 b 
-  assignArray (Tuple [ns1,ns2]) (a,b) ix =
-    do
-      assignArray ns1 a ix 
-      assignArray ns2 b ix
 
-  assignArrayS (Tuple [ns1,ns2]) (a,b) ix =
-    do
-      nas1 <- assignArrayS ns1 a ix 
-      nas2 <- assignArrayS ns2 b ix
-      return $ Tuple [nas1,nas2]
-  assignScalar (Tuple [ns1,ns2]) (a,b) =
+instance (Storable a, Storable b) => Storable (a, b) where
+  names pre =
     do
-      assignScalar ns1 a 
+      (a' :: Names a) <- names pre
+      (b' :: Names b) <- names pre 
+      return $ Tuple a' b'
+  allocateArray (Tuple ns1 ns2)  n =
+      allocateArray ns1 n >> 
+      allocateArray ns2 n
+      
+  allocateVolatileArray (Tuple ns1 ns2)  n =
+      allocateVolatileArray ns1 n >> 
+      allocateVolatileArray ns2 n
+
+      
+  allocateScalar (Tuple ns1 ns2) =
+      allocateScalar ns1 >> 
+      allocateScalar ns2
+      
+  assignArray (Tuple ns1 ns2) (a,b) ix =
+      assignArray ns1 a ix >>
+      assignArray ns2 b ix
+      
+  warpAssignArray (Tuple ns1 ns2) warpID step (a,b) ix =
+      warpAssignArray ns1 warpID step a ix >>
+      warpAssignArray ns2 warpID step b ix
+      
+  
+  assignScalar (Tuple ns1 ns2) (a,b) =
+      assignScalar ns1 a >>
       assignScalar ns2 b  
-  pullFrom (Tuple [ns1,ns2]) n =
+
+  pullFrom (Tuple ns1 ns2) n =
     let p1 = pullFrom ns1 n
         p2 = pullFrom ns2 n
-    in Pull n (\ix -> (p1 ! ix, p2 ! ix))
-  readFrom (Tuple [ns1,ns2])  =
+    in mkPull n (\ix -> (p1 ! ix, p2 ! ix))
+       
+  warpPullFrom (Tuple ns1 ns2) warpID n
+    = let p1 = warpPullFrom ns1 warpID n
+          p2 = warpPullFrom ns2 warpID n
+      in mkPull n (\ix -> (p1 ! ix, p2 ! ix)) 
+
+  readFrom (Tuple ns1 ns2)  =
     let p1 = readFrom ns1
         p2 = readFrom ns2
     in (p1,p2)
 
-  pullFromS (Tuple [ns1,ns2]) n =
-    let p1 = pullFromS ns1 n
-        p2 = pullFromS ns2 n
-    in Pull n (\ix -> (p1 ! ix, p2 ! ix))
-
-
----------------------------------------------------------------------------
--- Global Memory
----------------------------------------------------------------------------
-
-class GlobalMemoryOps a where
-  outputs   :: a -> GProgram Names
-  assignOut :: Names -> a -> Exp Word32 -> Program Thread ()
+  
 
 
-instance Scalar a => GlobalMemoryOps (Exp a) where
-  outputs a =
+instance (Storable a, Storable b, Storable c) => Storable (a, b, c) where
+  names pre =
     do
-      name <- Output $ Pointer $ typeOf a
-      return (Single name) 
-  assignOut (Single name) a ix = Assign name [ix] a
+      (a :: Names a) <- names pre 
+      (b :: Names b) <- names pre 
+      (c :: Names c) <- names pre 
+      return $ Triple a b c
+      
+  allocateArray (Triple ns1 ns2 ns3) n =
+      allocateArray ns1 n >>
+      allocateArray ns2 n >> 
+      allocateArray ns3 n
+      
+  allocateVolatileArray (Triple ns1 ns2 ns3) n =
+      allocateVolatileArray ns1 n >> 
+      allocateVolatileArray ns2 n >> 
+      allocateVolatileArray ns3 n 
+      
+  allocateScalar (Triple ns1 ns2 ns3) =
+      allocateScalar ns1 >>
+      allocateScalar ns2 >>
+      allocateScalar ns3
+      
+  assignArray (Triple ns1 ns2 ns3) (a,b,c) ix =
+      assignArray ns1 a ix >>
+      assignArray ns2 b ix >>
+      assignArray ns3 c ix
+      
+  warpAssignArray (Triple ns1 ns2 ns3) warpID step (a,b,c) ix =
+      warpAssignArray ns1 warpID step a ix >>
+      warpAssignArray ns2 warpID step b ix >>
+      warpAssignArray ns3 warpID step c ix 
+  
 
+  assignScalar (Triple ns1 ns2 ns3) (a,b,c) =
+      assignScalar ns1 a >>
+      assignScalar ns2 b >>
+      assignScalar ns3 c
+      
+  pullFrom (Triple ns1 ns2 ns3) n =
+    let p1 = pullFrom ns1 n
+        p2 = pullFrom ns2 n
+        p3 = pullFrom ns3 n
+    in mkPull n (\ix -> (p1 ! ix, p2 ! ix,p3 ! ix))
+       
+  warpPullFrom (Triple ns1 ns2 ns3) warpID n
+    = let p1 = warpPullFrom ns1 warpID n
+          p2 = warpPullFrom ns2 warpID n
+          p3 = warpPullFrom ns3 warpID n
+      in mkPull n (\ix -> (p1 ! ix, p2 ! ix, p3 ! ix)) 
 
-instance (GlobalMemoryOps a, GlobalMemoryOps b)
-         => GlobalMemoryOps (a,b) where
-  outputs (a,b) =
-    do
-      na <- outputs a
-      nb <- outputs b
-      return (Tuple [na,nb]) 
-  assignOut (Tuple [n1,n2]) (a,b) ix =
-    do
-      assignOut n1 a ix 
-      assignOut n2 b ix
+  readFrom (Triple ns1 ns2 ns3)  =
+    let p1 = readFrom ns1
+        p2 = readFrom ns2
+        p3 = readFrom ns3
+    in (p1,p2,p3)
+ 
diff --git a/Obsidian/Mutable.hs b/Obsidian/Mutable.hs
new file mode 100644
--- /dev/null
+++ b/Obsidian/Mutable.hs
@@ -0,0 +1,246 @@
+{-# LANGUAGE ScopedTypeVariables,
+             TypeFamilies,
+             EmptyDataDecls,
+             FlexibleInstances #-}
+
+{- Joel Svensson 2013 -}
+
+module Obsidian.Mutable ( Mutable(Mutable) 
+                        , Shared
+                        , Global 
+                        , newS
+                        , forceTo
+                        , writeTo
+                        , pullFrom
+                        , atomicInc
+--                         , mutlen -- hack
+                        , namedMutable
+                        , undefinedMutable
+                        )  where 
+
+
+
+import Obsidian.Exp
+import Obsidian.Types
+import Obsidian.Globs
+import Obsidian.Program
+import Obsidian.Memory
+import Obsidian.Names
+import Obsidian.Array
+import Obsidian.Atomic 
+
+import Data.Word
+
+{-
+  Todo: Think about Global vs Shared.
+  Todo: Add creation of mutable global arrays. 
+
+  Todo: Make mutable interface (atomic ops) very low-level
+
+
+  TODO: Rethink. Have two sepparate types of mutable arrays.
+        Also Skip the Type family magic if possible.
+        Make both kinds of Mutable arrays an instance of Array 
+-} 
+
+---------------------------------------------------------------------------
+-- Mutable arrays 
+---------------------------------------------------------------------------
+--
+-- Global mutable arrays can only be passed as inputs to a function. 
+-- Shared mutable arrays may be created using newS 
+--
+
+data Shared
+data Global
+  
+-- A mutable array has an attached location.
+-- Either it recides in Global or in Shared memory. 
+data Mutable mloc s a = Mutable s (Names a)
+
+type MShared a = Mutable Shared Word32 a
+type MGlobal a = Mutable Global EWord32 a
+
+instance ArrayLength (Mutable Shared) where
+  len (Mutable n _) = n
+
+  
+
+
+namedMutable s v = Mutable v (Single s)
+undefinedMutable v = Mutable v undefined 
+---------------------------------------------------------------------------
+-- Create Mutable Shared memory arrays
+--   # allocates shared memory
+---------------------------------------------------------------------------
+
+newS :: Storable a => SPush Block a -> Program Block (Mutable Shared Word32 a)
+newS arr = do
+  (snames :: Names a) <- names "arr"
+  allocateArray snames n
+  let mut = Mutable n snames
+  writeTo mut arr
+  return $ mut -- Mutable n snames 
+  where
+    n = len arr
+
+---------------------------------------------------------------------------
+-- forceTo & writeTo
+---------------------------------------------------------------------------
+-- Much Hacking here 
+writeTo :: Storable a
+           => Mutable Shared Word32 a
+           -> Push Block Word32 a
+           -> Program Block ()
+writeTo (Mutable n snames) p 
+  | n <= m =  p <: assignArray snames
+  | otherwise = error "WriteTo: Incompatible sizes" 
+  where
+    m = len p
+   
+    
+-- Add forceTo with offsets (why? just thought it might be useful)
+forceTo :: Storable a
+           => Mutable Shared Word32 a
+           -> Push Block Word32 a
+           -> Program Block ()
+forceTo m arr =
+  do
+    writeTo m arr
+    Sync 
+---------------------------------------------------------------------------
+-- pullFrom 
+---------------------------------------------------------------------------
+
+toPull :: Storable a => Mutable Shared Word32 a -> SPull  a
+toPull (Mutable n snames) = pullFrom snames n  
+
+
+---------------------------------------------------------------------------
+-- Atomics
+---------------------------------------------------------------------------
+-- | Increment atomically 
+atomicInc :: forall mloc a s t . AtomicInc a
+             => EWord32  
+             -> Mutable mloc s (Exp a)
+             -> TProgram ()
+atomicInc ix (Mutable n noms) = mapNamesM_ f noms
+  where
+    f nom = atomicOp nom ix (AtomicInc  :: Atomic a) >> return ()
+  
+
+-- | Add atomically 
+atomicAdd :: forall mloc a s. AtomicAdd a
+             => EWord32
+             -> Exp a 
+             -> Mutable mloc s  (Exp a)
+             -> TProgram ()
+atomicAdd ix v (Mutable n noms) = mapNamesM_ f noms
+  where
+    f nom = atomicOp nom ix (AtomicAdd v) >> return ()
+  
+
+-- | Subtract atomically     
+atomicSub :: forall mloc a s. AtomicSub a
+             => EWord32
+             -> Exp a
+             -> Mutable mloc s (Exp a)
+             -> TProgram ()
+atomicSub ix v (Mutable n noms) = mapNamesM_ f noms
+                               
+  where
+    f nom = atomicOp nom ix (AtomicSub v) >> return ()
+
+
+-- Special case ? No. 
+atomicExch :: forall mloc a s. AtomicExch a
+              => EWord32
+              -> Exp a
+              -> Mutable mloc s (Exp a)
+              -> TProgram ()
+atomicExch ix v  (Mutable n (Single nom)) = f nom
+  where
+    f nom = atomicOp nom ix (AtomicExch v) 
+                         
+  
+{-
+
+---------------------------------------------------------------------------
+atomicExch()
+
+int atomicExch(int* address, int val);
+unsigned int atomicExch(unsigned int* address,
+                        unsigned int val);
+unsigned long long int atomicExch(unsigned long long int* address,
+                                  unsigned long long int val);
+float atomicExch(float* address, float val);
+
+---------------------------------------------------------------------------
+atomicMin()
+
+int atomicMin(int* address, int val);
+unsigned int atomicMin(unsigned int* address,
+                       unsigned int val);
+unsigned long long int atomicMin(unsigned long long int* address,
+                                 unsigned long long int val);
+
+---------------------------------------------------------------------------
+atomicMax()
+
+int atomicMax(int* address, int val);
+unsigned int atomicMax(unsigned int* address,
+                       unsigned int val);
+unsigned long long int atomicMax(unsigned long long int* address,
+                                 unsigned long long int val);
+
+
+---------------------------------------------------------------------------
+atomicInc()
+
+unsigned int atomicInc(unsigned int* address,
+                       unsigned int val);
+
+---------------------------------------------------------------------------
+atomicDec()
+
+unsigned int atomicDec(unsigned int* address,
+                       unsigned int val);
+
+---------------------------------------------------------------------------
+atomicCAS()
+
+int atomicCAS(int* address, int compare, int val);
+unsigned int atomicCAS(unsigned int* address,
+                       unsigned int compare,
+                       unsigned int val);
+unsigned long long int atomicCAS(unsigned long long int* address,
+                                 unsigned long long int compare,
+                                 unsigned long long int val);
+
+---------------------------------------------------------------------------
+atomicAnd()
+
+int atomicAnd(int* address, int val);
+unsigned int atomicAnd(unsigned int* address,
+                       unsigned int val);
+unsigned long long int atomicAnd(unsigned long long int* address,
+                                 unsigned long long int val);
+
+---------------------------------------------------------------------------
+atomicOr()
+
+int atomicOr(int* address, int val);
+unsigned int atomicOr(unsigned int* address,
+                      unsigned int val);
+unsigned long long int atomicOr(unsigned long long int* address,
+                                unsigned long long int val);
+---------------------------------------------------------------------------
+atomicXor()
+
+int atomicXor(int* address, int val);
+unsigned int atomicXor(unsigned int* address,
+                       unsigned int val);
+unsigned long long int atomicXor(unsigned long long int* address,
+                                 unsigned long long int val);
+
+-} 
diff --git a/Obsidian/Names.hs b/Obsidian/Names.hs
--- a/Obsidian/Names.hs
+++ b/Obsidian/Names.hs
@@ -1,20 +1,23 @@
-
+{-# LANGUAGE GADTs #-}
 
-module Obsidian.Names where
+module Obsidian.Names (Names(..), mapNamesM_) where
 
 import Obsidian.Globs
 
-data Tree a = None
-            | Single a
-            | Tuple [Tree a]
 
-type Names = Tree Name
-
-type NameInfo = Tree (Kind,Name)
+data Names a where
+  Single :: Name -> Names a
+  Tuple  :: Names a -> Names b -> Names (a,b)
+  Triple :: Names a -> Names b -> Names c -> Names (a,b,c) 
 
-data Kind = Var | Arr
+---------------------------------------------------------------------------
+-- helpers
+---------------------------------------------------------------------------
+mapNamesM_ :: Monad m => (Name -> m ()) -> Names a -> m ()
+mapNamesM_ f (Single nom)  = f nom
+mapNamesM_ f (Tuple n1 n2) = mapNamesM_ f n1 >>
+                             mapNamesM_ f n2
+mapNamesM_ f (Triple n1 n2 n3) = mapNamesM_ f n1 >>
+                                 mapNamesM_ f n2 >>
+                                 mapNamesM_ f n3
 
-instance Functor Tree where
-  fmap f None = None
-  fmap f (Single a) = Single $ f a
-  fmap f (Tuple ts) = Tuple $ map (fmap f) ts
diff --git a/Obsidian/Program.hs b/Obsidian/Program.hs
--- a/Obsidian/Program.hs
+++ b/Obsidian/Program.hs
@@ -1,16 +1,43 @@
 {- Joel Svensson 2012,2013
 
    Notes:
+   2013-04-02: Added a Break statement to the language.
+               Use it to break out of sequential loops.
    2013-01-08: removed number-of-blocks field from ForAllBlocks
 
 -}
 
-{-# LANGUAGE GADTs,
-             FlexibleInstances,
-             EmptyDataDecls #-} 
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE EmptyDataDecls #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE ScopedTypeVariables #-}
 
 
-module Obsidian.Program  where 
+
+
+module Obsidian.Program  (
+  -- Hierarchy 
+  Thread, Block, Grid, Step, Zero, Warp, 
+  -- Program type
+  -- CoreProgram(..),
+  Program(..), -- all exported.. for now
+  TProgram, BProgram, GProgram, WProgram(..), 
+
+  -- Class
+  Sync, 
+  
+  -- helpers
+  printPrg,
+  runPrg,
+  uniqueNamed, uniqueNamed_,  
+
+  assign, allocate, declare,
+  atomicOp, 
+  -- Programming interface
+  seqFor, forAll, forAll2, seqWhile, sync, distrPar,
+  ) where 
  
 import Data.Word
 import Data.Monoid
@@ -25,26 +52,33 @@
 import Data.Supply
 import System.IO.Unsafe
 
+import Control.Monad
+import Control.Applicative
+
+
 ---------------------------------------------------------------------------
 -- Thread/Block/Grid 
 ---------------------------------------------------------------------------
-data Thread
-data Block
-data Grid 
 
-data PT a where
-  Thread :: PT Thread
-  Block  :: PT Block
-  Grid   :: PT Grid 
 
-type Identifier = Int 
 
----------------------------------------------------------------------------
--- Obsidian
+-- A hierarchy! 
+data Step a -- A step in the hierarchy
+data Zero 
+  
+type Thread = Zero
+type Warp   = Step Thread
+type Block  = Step Warp
+type Grid   = Step Block
+
+
+
+
+
 ---------------------------------------------------------------------------
-data Obsidian a = Obsidian (Program Grid a)
-      
 
+type Identifier = Int
+                  
 ---------------------------------------------------------------------------
 -- Program datatype
 --------------------------------------------------------------------------
@@ -58,77 +92,95 @@
             -> (Exp a)
             -> Program Thread ()
            
-            
+  -- 4 March 2014, Changed so that AtOp does not return a result. 
+  -- Change this back later if an application requires. 
   AtomicOp :: Scalar a
-              => Name 
-              -> Exp Word32
-              -> Atomic a
-              -> Program Thread (Exp a)
+              => Name        -- Array name 
+              -> Exp Word32  -- Index to operate on 
+              -> Atomic a    -- Atomic operation to perform 
+              -> Program Thread ()
 
   Cond :: Exp Bool
-          -> Program t ()
-          -> Program t ()
-  
-  -- DONE: Code generation for this.
-  -- TODO: Generalize this loop! (Replace Thread with t) 
-  SeqFor :: Exp Word32 -> (Exp Word32 -> Program t a)
-            -> Program t a
- 
-  ForAll :: (Exp Word32) 
-            -> (Exp Word32 -> Program Thread a)
-            -> Program Block a 
+          -> Program Thread ()
+          -> Program Thread ()
+                
+  SeqWhile :: Exp Bool ->
+              Program Thread () ->
+              Program Thread () 
+              
+  Break  :: Program Thread () 
 
-  {-
-     I'm not sure about this constructor.
-     As I see it programs from which we generate a kernel
-     must be wrapped in one of these ForAllBlocks.
-     Programs with sequences of 'ForAllBlocks' are problematic.
 
-     Maybe a (ForAllBlocks n f *>* ForAllBlocks m g) Program
-     should be split into two kernels. 
-  -} 
-  ForAllBlocks :: (Exp Word32) -> (Exp Word32 -> Program Block a) 
-                  -> Program Grid a
+  -- use threads along one level
+  -- Warp, Block, Grid. 
+  ForAll ::  EWord32 
+            -> (EWord32 -> Program Thread ())
+            -> Program t () -- (really atleast Step t) ! 
 
-  ForAllThreads :: (Exp Word32) -> (Exp Word32 -> Program Thread a)
-                   -> Program Grid a 
+  
+  -- Distribute over Warps yielding a Block
+  -- Distribute over Blocks yielding a Grid 
+  DistrPar :: EWord32
+           -> (EWord32 -> Program t ())
+           -> Program (Step t) ()
+           -- (really Step t -> Step (Step t) ) 
+                           
+  SeqFor :: EWord32 -> (EWord32 -> Program t ())
+            -> Program t ()
 
-  -- Allocate shared memory in each MP
+                           
+ 
+  --        #w          warpId     
+  --NWarps :: EWord32 -> (EWord32 -> Program Warp ()) -> Program Block () 
 
-  -- TODO: Change the Liveness analysis to a two-pass algo
-  --       and remove the Allocate constructor. 
+  --WarpForAll :: EWord32 
+  --              -> (EWord32 -> Program Thread ()) 
+  --              -> Program Warp ()
+  -- WarpAllocate :: Name -> Word32 -> Type -> Program Warp ()  -- For now. 
+
+  -- Allocate shared memory in each MP
   Allocate :: Name -> Word32 -> Type -> Program t () 
 
   -- Automatic Variables
   Declare :: Name -> Type -> Program t () 
-              
-  {- About Output (Creates a named output array). 
-     This is similar to Allocate but concerning global arrays.
-
-     Since we cannot synchronize writes to a global array inside of an
-     kernel, global arrays will only be written as outputs of the kernel
-  -} 
-  Output   :: Type -> Program Grid Name
-  -- (Output may be replaced by AllocateG) 
-  
+                
   Sync     :: Program Block ()
-  -- Two very experimental threadfence constructs.
-  -- should correspond to cuda __threadfence();
-  -- and __threadfenceBlock(); 
-  ThreadFence :: Program Grid ()
-  ThreadFenceBlock :: Program Block () 
 
   -- Parallel composition of Programs
   -- TODO: Will I use this ? 
-  Par :: Program p () ->
-         Program p () ->
-         Program p () 
+  --Par :: Program p () ->
+  --       Program p () ->
+  --       Program p () 
 
   -- Monad
   Return :: a -> Program t a
   Bind   :: Program t a -> (a -> Program t b) -> Program t b
 
 ---------------------------------------------------------------------------
+-- Aliases 
+---------------------------------------------------------------------------
+type TProgram = Program Thread
+type WProgram = Program Warp 
+type BProgram = Program Block
+type GProgram = Program Grid 
+
+-- -- Programs are a reader monad 
+-- newtype Program t a = Program (EWord32 -> CoreProgram t a)
+
+-- instance Monad (Program t )where
+--     return x = Program $ \ _ -> return x
+--     -- :: WProgram a -> (a -> WProgram b) -> WProgram b
+--     (Program h) >>= f = Program
+--                          $ \w ->
+--                          do
+--                            a <- h w
+--                            let (Program g) = f a
+--                            g w
+
+
+-- core :: Program t a -> EWord32 ->  CoreProgram t a
+-- core (Program f) id = f id 
+---------------------------------------------------------------------------
 -- Helpers 
 --------------------------------------------------------------------------- 
 uniqueSM = do
@@ -139,46 +191,76 @@
   id <- Identifier
   return $ pre ++ show id 
 
+uniqueNamed_ pre = do id <- Identifier
+                      return $ pre ++ show id 
 ---------------------------------------------------------------------------
--- forAll and forAllN
+-- Memory 
 ---------------------------------------------------------------------------
---forAll :: (Exp Word32 -> Program Thread ()) -> Program Block () 
---forAll f = ForAll Nothing f
+assign :: Scalar a => Name -> [Exp Word32] -> (Exp a) -> Program Thread ()
+assign nom ix e = Assign nom ix e 
 
-forAll :: Exp Word32 -> (Exp Word32 -> Program Thread a) -> Program Block a
-forAll n f = ForAll n f
+allocate :: Name -> Word32 -> Type -> Program t () 
+allocate nom l t = Allocate nom l t 
 
-(*||*) = Par
+declare :: Name -> Type -> Program t ()
+declare nom t = Declare nom t 
 
 ---------------------------------------------------------------------------
--- SeqFor
+-- atomicOp 
 ---------------------------------------------------------------------------
-seqFor :: Exp Word32 -> (Exp Word32 -> Program t a)
-            -> Program t a
-seqFor (Literal 1) f = f 0
-seqFor n f = SeqFor n f
-
+atomicOp :: Scalar a
+            => Name        -- Array name 
+            -> Exp Word32  -- Index to operate on 
+            -> Atomic a    -- Atomic operation to perform 
+            -> Program Thread ()
+atomicOp nom ix atop = AtomicOp nom ix atop 
 
 ---------------------------------------------------------------------------
--- forAllT
---------------------------------------------------------------------------- 
--- When we know that all threads are independent and
--- independent of "blocksize".
--- Also any allocation of local storage is impossible.
--- Composition of something using forAllT with something
--- that performs local computations is impossible.
--- Using the hardcoded BlockDim may turn out to be a problem when
--- we want to compute more than one thing per thread (may be fine though). 
-forAllT :: (Exp Word32) -> (Exp Word32 -> Program Thread a)
-           -> Program Grid a
-forAllT n f = ForAllThreads n 
-            $ \gtid -> f gtid 
+-- forAll 
+---------------------------------------------------------------------------
+forAll :: EWord32 -> (EWord32 -> Program Thread ()) -> Program t ()
+forAll n f = ForAll n $ \ix -> f ix
+  
+-- forAll :: EWord32 -> (EWord32 -> Program t ()) -> Program (Step t) ()
+-- forAll n f = Program $ \id -> ForAll n $ \ix -> core (f ix) id
 
+forAll2 :: EWord32
+           -> EWord32
+           -> (EWord32 -> EWord32 -> Program Thread ())
+           -> Program (Step (Step t)) ()
+forAll2 b n f =  DistrPar b $ \bs -> ForAll n $ \ix -> f bs ix 
 
+-- forAll2 :: EWord32
+--            -> EWord32
+--            -> (EWord32 -> EWord32 -> Program t ())
+--            -> Program (Step (Step t)) ()
+-- forAll2 b n f =  forAll b $ \bs -> forAll n (f bs)
 
-forAllBlocks = ForAllBlocks
+distrPar :: EWord32
+           -> (EWord32 -> Program t ())
+           -> Program (Step t) ()
+distrPar b f = DistrPar b $ \bs -> f bs
 
 ---------------------------------------------------------------------------
+-- warpForAll 
+---------------------------------------------------------------------------
+--warpForAll :: EWord32 -> (EWord32 -> Program Thread ()) -> Program Warp ()
+-- warpForAll n f = Program $ \id -> WarpForAll n $ \ix -> core (f ix) id 
+ 
+---------------------------------------------------------------------------
+-- seqFor
+---------------------------------------------------------------------------
+seqFor :: EWord32 -> (EWord32 -> Program t ()) -> Program t ()
+seqFor (Literal 1) f = f 0
+seqFor n f = SeqFor n $ \ix -> f ix
+
+---------------------------------------------------------------------------
+-- seqWhile
+---------------------------------------------------------------------------
+seqWhile :: Exp Bool -> Program Thread () -> Program Thread ()
+seqWhile b prg = SeqWhile b prg 
+
+---------------------------------------------------------------------------
 -- Monad
 --------------------------------------------------------------------------
 instance Monad (Program t) where
@@ -186,53 +268,98 @@
   (>>=) = Bind
 
 ---------------------------------------------------------------------------
--- Aliases 
+-- Functor
 ---------------------------------------------------------------------------
-type TProgram = Program Thread
-type BProgram = Program Block
-type GProgram = Program Grid 
+instance Functor (Program t) where
+  fmap g fa = do {a <- fa; return $ g a}
 
 ---------------------------------------------------------------------------
--- runPrg (fix types here, Integer!)
+-- Applicative 
 ---------------------------------------------------------------------------
+instance Applicative (Program t) where
+  pure = return
+  ff <*> fa = 
+    do
+      f <- ff
+      fmap f fa
+
+---------------------------------------------------------------------------
+-- Class Sync
+---------------------------------------------------------------------------
+class Sync t where
+  sync :: Program t () 
+
+instance Sync Warp where
+  sync = return ()
+
+instance Sync Thread where
+  sync = return ()
+
+instance Sync Block where
+  sync = Sync
+
+instance Sync Grid where
+  sync = error "sync: not implemented on grid computations" 
+  -- (implement this using counters and locks)
+
+---------------------------------------------------------------------------
+-- runPrg (RETHINK!) (Works for Block programs, but all?)
+---------------------------------------------------------------------------
 runPrg :: Int -> Program t a -> (a,Int)
 runPrg i Identifier = (i,i+1)
+
+-- Maybe these two are the most interesting cases!
+-- Return may for example give an array. 
 runPrg i (Return a) = (a,i)
 runPrg i (Bind m f) =
   let (a,i') = runPrg i m
-  in runPrg i' (f a) 
+  in runPrg i' (f a)
+     
 runPrg i (Sync) = ((),i)
 runPrg i (ForAll n ixf) =
   let (p,i') = runPrg i (ixf (variable "tid")) 
-  in  (p,i') 
-runPrg i (Allocate id _ _ ) = ((),i)
+  in  (p,i')
+runPrg i (DistrPar n f) =
+  let (p,i') = runPrg i (f (variable "DUMMY"))
+  in (p,i')
+-- What can this boolean depend upon ? its quite general!
+--  (we know p returns a ()... ) 
+runPrg i (Cond b p) = ((),i) 
+runPrg i (Declare _ _) = ((),i)
+runPrg i (Allocate _ _ _ ) = ((),i)
 runPrg i (Assign _ _ a) = ((),i) -- Probaby wrong.. 
-runPrg i (AtomicOp _ _ _) = (variable ("new"++show i),i+1)
-     
+runPrg i (AtomicOp _ _ _) = ((),i) -- variable ("new"++show i),i+1)
+
+{- What do I want from runPrg ?
+
+   # I want to it to "work" for all block programs (no exceptions)
+   # I want a BProgram (Pull a) to return a Pull array of "correct length)
+-}
+
+                            
 ---------------------------------------------------------------------------
--- printPrg
+-- printPrg (REIMPLEMENT) xs
 ---------------------------------------------------------------------------
+printPrg :: Program t a -> String
 printPrg prg = (\(_,x,_) -> x) $ printPrg' 0 prg
 
 printPrg' :: Int -> Program t a -> (a,String,Int)
 printPrg' i Identifier = (i,"getId;\n",i+1) 
--- printPrg' i Skip = ((),";\n", i)
 printPrg' i (Assign n ix e) =
   ((),n ++ "[" ++ show ix ++ "] = " ++ show e ++ ";\n", i) 
 printPrg' i (AtomicOp n ix e) =
   let newname = "r" ++ show i
-  in (variable newname,
-      newname ++ " = " ++ printAtomic e ++
-      "( " ++ n ++ "[" ++ show ix ++ "])\n",i+1)
+  --in (variable newname,
+  --    newname ++ " = " ++ printAtomic e ++
+  --    "( " ++ n ++ "[" ++ show ix ++ "])\n",i+1)
+  in ((), printAtomic e ++
+          "( " ++ n ++ "[" ++ show ix ++ "])\n",i+1)
 printPrg' i (Allocate id n t) =
   let newname = id -- "arr" ++ show id
   in ((),newname ++ " = malloc(" ++ show n ++ ");\n",i+1)
 printPrg' i (Declare id t) =
   let newname = id -- "arr" ++ show id
   in ((),show t ++ " " ++ newname ++ "\n",i+1)
-printPrg' i (Output t) =
-  let newname = "globalOut" ++ show i
-  in (newname,newname ++ " = new Global output;\n",i+1)
 printPrg' i (SeqFor n f) =
   let (a,prg2,i') = printPrg' i (f (variable "i"))
       
@@ -248,12 +375,12 @@
        "par (i in 0.." ++ show n ++ ")" ++
        "{\n" ++ prg2 ++ "\n}",
        i')
-printPrg' i (ForAllBlocks n f) =
-  let (d,prg2,i') = printPrg' i (f (variable "BIX"))
-  in (d, 
-      "blocks (i)" ++
-      "{\n" ++ prg2 ++ "\n}",
-      i')
+--printPrg' i (ForAllBlocks n f) =
+--  let (d,prg2,i') = printPrg' i (f (variable "BIX"))
+--  in (d, 
+--      "blocks (i)" ++
+--      "{\n" ++ prg2 ++ "\n}",
+--      i')
 printPrg' i (Return a) = (a,"MonadReturn;\n",i)
 printPrg' i (Bind m f) =
   let (a1, str1,i1) = printPrg' i m
diff --git a/Obsidian/Run/CUDA/Exec.hs b/Obsidian/Run/CUDA/Exec.hs
new file mode 100644
--- /dev/null
+++ b/Obsidian/Run/CUDA/Exec.hs
@@ -0,0 +1,387 @@
+{-# LANGUAGE TypeOperators,
+             ScopedTypeVariables,
+             TypeFamilies,
+             TypeSynonymInstances,
+             FlexibleInstances #-} 
+
+module Obsidian.Run.CUDA.Exec where
+
+---------------------------------------------------------------------------
+--
+-- Low level interface to CUDA functionality from Obsidian
+--
+---------------------------------------------------------------------------
+
+
+import qualified Foreign.CUDA.Driver as CUDA
+import qualified Foreign.CUDA.Driver.Device as CUDA
+import qualified Foreign.CUDA.Analysis.Device as CUDA
+import qualified Foreign.CUDA.Driver.Stream as CUDAStream
+
+import Obsidian.CodeGen.Reify
+import Obsidian.CodeGen.CUDA
+
+import Obsidian.Types -- experimental
+import Obsidian.Exp
+import Obsidian.Array
+import Obsidian.Program (Grid, GProgram)
+import Obsidian.Mutable
+
+import Foreign.Marshal.Array
+import Foreign.ForeignPtr.Unsafe -- (req GHC 7.6 ?)
+import Foreign.ForeignPtr hiding (unsafeForeignPtrToPtr)
+
+import qualified Data.Vector.Storable as V 
+
+import Data.Word
+import Data.Int
+import Data.Supply
+import Data.List
+import qualified Data.Map as M
+import Data.Maybe
+
+import System.IO.Unsafe
+import System.Process
+import System.Random.MWC
+import System.CPUTime.Rdtsc
+
+import Control.Monad.State
+
+
+
+debug = False
+
+---------------------------------------------------------------------------
+-- Tools 
+---------------------------------------------------------------------------
+mkRandomVec :: forall a.(V.Storable a, Variate a) => Int -> IO (V.Vector a)
+mkRandomVec k = withSystemRandom $ \g -> uniformVector g k :: IO (V.Vector a) 
+
+---------------------------------------------------------------------------
+-- An array located in GPU memory
+---------------------------------------------------------------------------
+data CUDAVector a = CUDAVector {cvPtr :: CUDA.DevicePtr a,
+                                cvLen :: Word32} 
+
+---------------------------------------------------------------------------
+-- Get a list of devices from the CUDA driver
+---------------------------------------------------------------------------
+getDevices :: IO [(CUDA.Device,CUDA.DeviceProperties)]
+getDevices = do
+  num <- CUDA.count 
+  devs <- mapM CUDA.device [0..num-1]
+  props <- mapM CUDA.props devs
+  return $ zip devs props
+
+---------------------------------------------------------------------------
+-- Print a Summary of a device's properties. 
+---------------------------------------------------------------------------
+propsSummary :: CUDA.DeviceProperties -> String
+propsSummary props = unlines
+  ["Device Name: " ++ CUDA.deviceName props,
+   "Compute Cap: " ++ show (CUDA.computeCapability props),
+   "Global Mem:  " ++ show (CUDA.totalGlobalMem props),
+   "Shared Mem/Block: " ++ show (CUDA.sharedMemPerBlock props),
+   "Registers/Block: "  ++ show (CUDA.regsPerBlock props),
+   "Warp Size: " ++ show (CUDA.warpSize props),
+   "Max threads/Block: " ++ show (CUDA.maxThreadsPerBlock props),
+   "Max threads/MP: " ++ show (CUDA.maxThreadsPerMultiProcessor props),
+   "Clock rate: " ++ show (CUDA.clockRate props),
+   "Num MP: " ++ show (CUDA.multiProcessorCount props),
+   "Mem bus width: " ++ show (CUDA.memBusWidth props)] 
+
+--------------------------------------------------------------------------
+-- Environment to run CUDA computations in.
+--  # Needs to keep track of generated and loaded functions etc. 
+---------------------------------------------------------------------------
+
+data CUDAState = CUDAState { csIdent :: Int,
+                             csCtx   :: CUDA.Context,
+                             csProps :: CUDA.DeviceProperties}
+
+type CUDA a =  StateT CUDAState IO a
+
+data Kernel = Kernel {kFun :: CUDA.Fun,
+                      kThreadsPerBlock :: Word32,
+                      kSharedBytes :: Word32}
+
+-- Change so that the type parameter to KernelT
+-- represents the "captured" type, with CUDAVectors instead of Pull, Push vectors. 
+data KernelT a = KernelT {ktFun :: CUDA.Fun,
+                          ktThreadsPerBlock :: Word32,
+                          ktSharedBytes :: Word32,
+                          ktInputs :: [CUDA.FunParam],
+                          ktOutput :: [CUDA.FunParam] }
+
+---------------------------------------------------------------------------
+-- Kernel Input and Output classes
+---------------------------------------------------------------------------
+class KernelI a where
+  type KInput a 
+  addInParam :: KernelT (KInput a -> b) -> a -> KernelT b
+
+class KernelM a where
+  type KMutable a
+  addMutable :: KernelT (KMutable a -> b) -> a -> KernelT b 
+  
+class KernelO a where
+  type KOutput a 
+  addOutParam :: KernelT (KOutput a) -> a -> KernelT ()
+
+instance KernelI Int32 where
+  type KInput Int32 = Exp Int32
+  addInParam (KernelT f t s i o) a =
+    KernelT f t s (i ++ [CUDA.IArg $ fromIntegral a]) o 
+
+instance KernelI Word32 where
+  type KInput Word32 = Exp Word32
+  addInParam (KernelT f t s i o) a =
+    KernelT f t s (i ++ [CUDA.IArg $ fromIntegral a]) o 
+
+instance Scalar a => KernelI (CUDAVector a) where
+  type KInput (CUDAVector a) = DPull (Exp a) 
+  addInParam (KernelT f t s i o) b =
+    KernelT f t s (i ++ [CUDA.VArg (cvPtr b),
+                         CUDA.IArg $ fromIntegral (cvLen b)]) o
+
+instance Scalar a => KernelM (CUDAVector a) where
+  type KMutable (CUDAVector a) = Mutable Global EW32 (Exp a) 
+  addMutable (KernelT f t s i o) b =
+    KernelT f t s (i ++ [CUDA.VArg (cvPtr b),
+                         CUDA.IArg $ fromIntegral (cvLen b)]) o
+
+instance Scalar a => KernelO (CUDAVector a) where
+  type KOutput (CUDAVector a) = DPush Grid (Exp a) 
+  addOutParam (KernelT f t s i o) b =
+    KernelT f t s i (o ++ [CUDA.VArg (cvPtr b)])
+
+---------------------------------------------------------------------------
+-- (<>) apply a kernel to an input
+---------------------------------------------------------------------------
+(<>) :: KernelI a
+        => (Word32,KernelT (KInput a -> b)) -> a -> (Word32,KernelT b)
+(<>) (blocks,kern) a = (blocks,addInParam kern a)
+
+---------------------------------------------------------------------------
+-- Assign a mutable input/output to a kernel
+--------------------------------------------------------------------------- 
+(<:>) :: KernelM a
+         => (Word32, KernelT (KMutable a -> b)) -> a -> (Word32, KernelT b)
+(<:>) (blocks,kern) a = (blocks, addMutable kern a)
+
+
+---------------------------------------------------------------------------
+-- Execute a kernel and store output to an array
+---------------------------------------------------------------------------
+(<==) :: KernelO b => b -> (Word32, KernelT (KOutput b)) -> CUDA ()
+(<==) o (nb,kern) =
+  do
+    let k = addOutParam kern o
+    lift $ CUDA.launchKernel
+      (ktFun k)
+      (fromIntegral nb,1,1)
+      (fromIntegral (ktThreadsPerBlock k), 1, 1)
+      (fromIntegral (ktSharedBytes k))
+      Nothing -- stream
+      (ktInputs k ++ ktOutput k) -- params
+    
+
+-- | A variant that returns kernel timing information as well.
+-- TODO: We probably can remove this:
+-- (<==!) :: KernelO b => b -> (Word32, KernelT (KOutput b)) -> CUDA Word64
+-- (<==!) o (nb,kern) =
+--   do
+--     let k = addOutParam kern o
+--     t1 <- lift rdtsc 
+--     lift $ CUDA.launchKernel
+--       (ktFun k)
+--       (fromIntegral nb,1,1)
+--       (fromIntegral (ktThreadsPerBlock k), 1, 1)
+--       (fromIntegral (ktSharedBytes k))
+--       Nothing -- stream
+--       (ktInputs k ++ ktOutput k) -- params
+--     lift $ CUDA.sync
+--     t2 <- lift rdtsc
+--     return (t2 - t1) 
+syncAll :: CUDA () 
+syncAll = lift $ CUDA.sync
+
+-- Tweak these 
+infixl 4 <>
+infixl 3 <==
+-- infixl 3 <==!
+
+---------------------------------------------------------------------------
+-- Execute a kernel that has no Output ( a -> GProgram ()) KernelT (GProgram ()) 
+---------------------------------------------------------------------------
+exec :: (Word32, KernelT (GProgram ())) -> CUDA ()
+exec (nb, k) =
+ lift $ CUDA.launchKernel
+     (ktFun k)
+     (fromIntegral nb,1,1)
+     (fromIntegral (ktThreadsPerBlock k), 1, 1)
+     0 {- (fromIntegral (ktSharedBytes k)) -}
+     Nothing -- stream
+     (ktInputs k)
+---------------------------------------------------------------------------
+-- Get a fresh identifier
+---------------------------------------------------------------------------
+newIdent :: CUDA Int
+newIdent =
+  do
+    i <- gets csIdent
+    modify (\s -> s {csIdent = i+1 }) 
+    return i
+
+
+---------------------------------------------------------------------------
+-- Run a CUDA computation
+---------------------------------------------------------------------------
+withCUDA p =
+  do
+    CUDA.initialise []
+    devs <- getDevices
+    case devs of
+      [] -> error "No CUDA device found!" 
+      (x:xs) ->
+        do 
+          ctx <- CUDA.create (fst x) [CUDA.SchedAuto] 
+          (a,_) <- runStateT p (CUDAState 0 ctx (snd x)) 
+          CUDA.destroy ctx
+          return a
+
+---------------------------------------------------------------------------
+-- Capture without an inputlist! 
+---------------------------------------------------------------------------    
+
+-- | Compile a program to a CUDA kernel and then load it into memory.
+capture :: ToProgram prg
+        => Word32  -- ^ Threads per block
+        -> prg     -- ^ Program to capture
+        -> CUDA (KernelT prg) 
+capture threadsPerBlock f =
+  do
+    i <- newIdent
+
+    props <- return . csProps =<< get
+    
+    let kn     = "gen" ++ show i
+        fn     = kn ++ ".cu"
+        cub    = fn ++ ".cubin"
+
+        prgstr = genKernel threadsPerBlock kn f
+        header = "#include <stdint.h>\n" -- more includes ? 
+
+    when debug $ 
+      do 
+        lift $ putStrLn $ prgstr
+
+    let arch = archStr props
+        
+    lift $ storeAndCompile arch (fn) (header ++ prgstr)
+    
+    mod <- liftIO $ CUDA.loadFile cub
+    fun <- liftIO $ CUDA.getFun mod kn 
+
+    {- After loading the binary into the running process
+       can I delete the .cu and the .cu.cubin ? -} 
+           
+    return $! KernelT fun threadsPerBlock 0 {-bytesShared-} [] []
+
+---------------------------------------------------------------------------
+-- useVector: Copies a Data.Vector from "Haskell" onto the GPU Global mem 
+--------------------------------------------------------------------------- 
+useVector :: V.Storable a =>
+             V.Vector a -> (CUDAVector a -> CUDA b) -> CUDA b
+useVector v f =
+  do
+    let (hfptr,n) = V.unsafeToForeignPtr0 v
+    
+    dptr <- lift $ CUDA.mallocArray n
+    let hptr = unsafeForeignPtrToPtr hfptr
+    lift $ CUDA.pokeArray n hptr dptr
+    let cvector = CUDAVector dptr (fromIntegral (V.length v)) 
+    b <- f cvector -- dptr
+    lift $ CUDA.free dptr
+    return b
+
+---------------------------------------------------------------------------
+-- allocaVector: allocates room for a vector in the GPU Global mem
+---------------------------------------------------------------------------
+allocaVector :: V.Storable a => 
+                Int -> (CUDAVector a -> CUDA b) -> CUDA b                
+allocaVector n f =
+  do
+    dptr <- lift $ CUDA.mallocArray n
+    let cvector = CUDAVector dptr (fromIntegral n)
+    b <- f cvector -- dptr
+    lift $ CUDA.free dptr
+    return b 
+
+---------------------------------------------------------------------------
+-- Allocate and fill with default value
+---------------------------------------------------------------------------
+allocaFillVector :: V.Storable a => 
+                Int -> a -> (CUDAVector a -> CUDA b) -> CUDA b                
+allocaFillVector n a f =
+  do
+    dptr <- lift $ CUDA.mallocArray n
+    lift $ CUDA.memset dptr n a 
+    let cvector = CUDAVector dptr (fromIntegral n)
+    b <- f cvector -- dptr
+    lift $ CUDA.free dptr
+    return b 
+
+---------------------------------------------------------------------------
+-- Fill a Vector
+---------------------------------------------------------------------------
+fill :: V.Storable a =>
+        CUDAVector a -> 
+        a -> CUDA ()
+fill (CUDAVector dptr n) a =
+  lift $ CUDA.memset dptr (fromIntegral n) a 
+
+---------------------------------------------------------------------------
+-- Peek in a CUDAVector (Simple "copy back")
+---------------------------------------------------------------------------
+peekCUDAVector :: V.Storable a => CUDAVector a -> CUDA [a]
+peekCUDAVector (CUDAVector dptr n) = 
+    lift $ CUDA.peekListArray (fromIntegral n) dptr
+    
+copyOut :: V.Storable a => CUDAVector a -> CUDA (V.Vector a)
+copyOut (CUDAVector dptr n) =
+  do
+    (fptr :: ForeignPtr a) <- lift $ mallocForeignPtrArray (fromIntegral n)
+    let ptr = unsafeForeignPtrToPtr fptr
+    lift $ CUDA.peekArray (fromIntegral n) dptr ptr
+    return $ V.unsafeFromForeignPtr fptr 0 (fromIntegral n)
+
+---------------------------------------------------------------------------
+-- Get the "architecture" of the present CUDA device
+---------------------------------------------------------------------------
+  
+archStr :: CUDA.DeviceProperties -> String
+archStr props = "-arch=sm_" ++ archStr' (CUDA.computeCapability props)
+  where
+    -- Updated for Cuda bindings version 0.5.0.0
+    archStr' (CUDA.Compute h l) = show h ++ show l
+    --archStr' (CUDA.Compute 1 0) = "10"
+    --archStr' (CUDA.Compute 1 2) = "12"
+    --archStr' (CUDA.Compute 2 0) = "20" 
+    --archStr' (CUDA.Compute 3 0) = "30"
+    --archStr' x = error $ "Unknown architecture: " ++ show x 
+    
+---------------------------------------------------------------------------
+-- Compile to Cubin (interface with nvcc)
+---------------------------------------------------------------------------
+storeAndCompile :: String -> FilePath -> String -> IO FilePath
+storeAndCompile arch fp code =
+  do
+    writeFile fp code
+    
+    let nfp = fp ++  ".cubin"
+
+    (_,_,_,pid) <-
+      createProcess (shell ("nvcc " ++ arch ++ " -cubin -o " ++ nfp ++ " " ++ fp))
+    exitCode <- waitForProcess pid
+    return nfp
diff --git a/Obsidian/SeqLoop.hs b/Obsidian/SeqLoop.hs
--- a/Obsidian/SeqLoop.hs
+++ b/Obsidian/SeqLoop.hs
@@ -15,88 +15,116 @@
 import Obsidian.Memory
 import Obsidian.Names
 
--- TODO: Add suitable allocs
+import Data.Word
+
 -- TODO: Rename module to something better
 
 ---------------------------------------------------------------------------
--- seqFold (actually reduce) 
+-- seqReduce 
 ---------------------------------------------------------------------------
-
-seqFold :: forall l a. (ASize l, MemoryOps a)
+-- | Sequential reduction of Pull array. Results in a for loop in generated code. 
+seqReduce :: Storable a
            => (a -> a -> a)
-           -> a
-           -> Pull l a
+           -> SPull a
            -> Program Thread a
-seqFold op init arr = do
-  ns  <- names "v" (undefined :: a) 
-  allocateScalar ns (undefined :: a)
-
-  assignScalar ns init  
-  -- Assign nom [] init  
-  SeqFor n $ \ ix ->
-    do
-      assignScalar ns (readFrom ns `op`  (arr ! ix))
-      return None
+seqReduce op arr =
+  do
+    (ns :: Names a)  <- names "v" 
+    allocateScalar ns 
 
+    assignScalar ns init  
+ 
+    SeqFor (n-1) $ \ ix ->
+      do
+        assignScalar ns (readFrom ns `op`  (arr ! (ix + 1)))
     
-  return $ readFrom ns
-  where 
+    return $ readFrom ns
+  where
     n = sizeConv$ len arr
+    init = arr ! 0 
 
 
 ---------------------------------------------------------------------------
--- Sequential scan
+-- Iterate
 ---------------------------------------------------------------------------
+-- | iterate a function. Results in a for loop in generated code. 
+seqIterate :: Storable a
+              => EWord32
+              -> (EWord32 -> a -> a)
+              -> a
+              -> Program Thread a
+seqIterate n f init =
+  do
+    (ns :: Names a)  <- names "v" 
+    allocateScalar ns 
 
-seqScan :: forall l a. (ASize l, MemoryOps a)
-           => (a -> a -> a)
-           -> Pull l a
-           -> Push Thread l a
-seqScan op (Pull n ixf)  =
-  Push n $ \wf -> do
-    ns <- names "v" (undefined :: a) 
-    allocateScalar ns (undefined :: a)
-    assignScalar ns (ixf 0)
-    wf (readFrom ns) 0 
-    SeqFor (sizeConv (n-1)) $ \ix -> do
-      wf (readFrom ns) ix                  
-      assignScalar ns  $ readFrom ns `op` (ixf (ix + 1))
-      return None
-                 
----------------------------------------------------------------------------
--- Sequential Map (here for uniformity) 
----------------------------------------------------------------------------
+    assignScalar ns init
+    SeqFor n $ \ix ->
+      do
+        assignScalar ns $ f ix (readFrom ns)
 
-seqMap :: forall l a b. ASize l
-          => (a -> b)
-          -> Pull l a
-          -> Push Thread l b
-seqMap f arr =
-  Push (len arr) $ \wf -> do
-    SeqFor (sizeConv (len arr)) $ \ix ->
-      wf (f (arr ! ix)) ix 
+    return $ readFrom ns
 
+---------------------------------------------------------------------------
+-- 
+---------------------------------------------------------------------------
+-- | iterate a function until a condition holds. Results in a while loop
+-- with a break in the generated code. 
+seqUntil :: Storable a
+            => (a -> a)
+            -> (a -> EBool)
+            -> a
+            -> Program Thread a
+seqUntil f p init =
+  do 
+    (ns :: Names a) <- names "v" 
+    allocateScalar ns 
 
+    assignScalar ns init
+    SeqWhile (p (readFrom ns)) $ 
+      do
+        (tmp :: Names a) <- names "t"
+        allocateScalar tmp
+        assignScalar tmp (readFrom ns) 
+        assignScalar ns $ f (readFrom tmp)
+    return $ readFrom ns
+  
 ---------------------------------------------------------------------------
--- Sequential Map and scan (generalisation of map + accum) 
+-- Sequential scan
 ---------------------------------------------------------------------------
+-- | Sequential scan over the elements in a pull array. Results in a for loop
+-- in the generated code. 
+seqScan :: Storable a
+           => (a -> a -> a)
+           -> SPull a
+           -> SPush Thread a
+seqScan op arr {-(Pull n ixf)-}  =
+  mkPush n $ \wf ->  do
+    (ns :: Names a) <- names "v" -- (ixf 0) 
+    allocateScalar ns -- (ixf 0)
+    assignScalar ns (arr ! 0)
+    wf (readFrom ns) 0
+    SeqFor (sizeConv (n-1)) $ \ix -> do
+      assignScalar ns  $ readFrom ns `op` (arr ! (ix + 1))
+      wf (readFrom ns) (ix+1)
+    where
+      n = len arr
 
-seqMapScan :: forall l a b acc. (ASize l, MemoryOps acc, MemoryOps b)
-              => (acc -> a -> (acc,b))
-              -> acc 
-              -> Pull l a
-              -> Push Thread l (acc,b)
-seqMapScan op acc (Pull n ixf)  =
-  Push n $ \wf -> do
-    ns <- names "v" (undefined :: b) 
-    allocateScalar ns (undefined :: b)
-    nacc <- names "v" (undefined :: acc)
-    allocateScalar nacc (undefined :: acc)
+-- | Sequential scan that takes a carry-in. 
+seqScanCin :: Storable a
+           => (a -> a -> a)
+           -> a -- cin  
+           -> SPull a
+           -> SPush Thread a
+seqScanCin op a arr {-(Pull n ixf)-} =
+  mkPush n $ \wf ->  do
+    (ns :: Names a) <- names "v" -- (ixf 0) 
+    allocateScalar ns -- (ixf 0)
+    assignScalar ns a -- (ixf 0)
+    -- wf (readFrom ns) 0 
+    SeqFor (sizeConv  n) $ \ix -> do
+      assignScalar ns  $ readFrom ns `op` (arr ! ix)
+      wf (readFrom ns) ix
+  where
+    n = len arr
     
-    assignScalar nacc acc
-
-    SeqFor (sizeConv n) $ \ix -> do
-      let (a,b) = op (readFrom nacc) (ixf ix)
-      wf (a,b) ix                  
-      assignScalar nacc a
-      return None
diff --git a/Obsidian/Types.hs b/Obsidian/Types.hs
--- a/Obsidian/Types.hs
+++ b/Obsidian/Types.hs
@@ -1,6 +1,15 @@
-module Obsidian.Types where 
+{-|
+Module      : Types
+Description : Type information, used internally by Obsídian.
+Copyright   : (c) Joel Svensson, 2014
+License     : BSD
+Maintainer  : bo.joel.svensson@gmail.com
+Stability   : experimental
 
+-}
 
+module Obsidian.Types where 
+
 ---------------------------------------------------------------------------
 -- Types
 ---------------------------------------------------------------------------
@@ -14,11 +23,12 @@
   | Word8 | Word16 | Word32 | Word64 
   | Float | Double                     
             
--- Used by CUDA, C And OpenCL generators        
-  | Pointer Type   -- C thing 
+-- Used by CUDA, C And OpenCL generators
+  | Volatile Type  -- For warp local computations. 
+  | Pointer Type   -- Pointer to a @type@ 
   | Global Type    -- OpenCL thing
   | Local Type     -- OpenCL thing
-  deriving (Eq, Show)
+  deriving (Eq, Ord, Show)
            
 typeSize Int8 = 1
 typeSize Int16 = 2