diff --git a/BLAS/Internal.hs b/BLAS/Internal.hs
--- a/BLAS/Internal.hs
+++ b/BLAS/Internal.hs
@@ -24,6 +24,15 @@
     checkMatMatOp,
     checkMatVecMult,
     checkMatMatMult,
+    checkMatVecMultAdd,
+    checkMatMatMultAdd,
+    checkMatVecSolv,
+    checkMatMatSolv,
+    checkMatVecSolvTo,
+    checkMatMatSolvTo,
+    checkSquare,
+    checkFat,
+    checkTall,
     diagStart,
     diagLen,    
     ) where
@@ -104,14 +113,14 @@
     | otherwise =
         sub o n'
 
-checkVecVecOp :: String -> Int -> Int -> IO ()
+checkVecVecOp :: String -> Int -> Int -> a -> a
 checkVecVecOp name n1 n2
     | n1 /= n2 =
-        ioError $ userError $ printf
+        error $ printf
             ("%s: x and y have different dimensions.  x has dimension `%d',"
              ++ " and y has dimension `%d'") name n1 n2
-    | otherwise =
-        return ()
+    | otherwise = id
+{-# INLINE checkVecVecOp #-}
 
 checkedRow ::  (Int,Int) -> (Int -> v) -> Int -> v
 checkedRow (m,n) row i 
@@ -168,31 +177,114 @@
         sub (i,j) (m',n')
 
 
-checkMatMatOp :: String -> (Int,Int) -> (Int,Int) -> IO ()
+checkMatMatOp :: String -> (Int,Int) -> (Int,Int) -> a -> a
 checkMatMatOp name mn1 mn2
     | mn1 /= mn2 =
-        ioError $ userError $ printf
+        error $ printf
             ("%s: x and y have different shapes.  x has shape `%s',"
              ++ " and y has shape `%s'") name (show mn1) (show mn2)
-    | otherwise =
-        return ()
+    | otherwise = id
         
-checkMatVecMult :: (Int,Int) -> Int -> IO ()
+checkMatVecMult :: (Int,Int) -> Int -> a -> a
 checkMatVecMult mn n
     | snd mn /= n =
-        ioError $ userError $ printf
+        error $ printf
             ("Tried to multiply a matrix with shape `%s' by a vector of dimension `%d'")
             (show mn) n
-    | otherwise =
-        return ()
+    | otherwise = id
         
-checkMatMatMult :: (Int,Int) -> (Int,Int) -> IO ()
+checkMatMatMult :: (Int,Int) -> (Int,Int) -> a -> a
 checkMatMatMult mk kn
     | snd mk /= fst kn =
-        ioError $ userError $ printf
+        error $ printf
             ("Tried to multiply a matrix with shape `%s' by a matrix with shape `%s'")
             (show mk) (show kn)
-    | otherwise =
-        return ()
+    | otherwise = id
 
-        
+checkMatVecMultAdd :: (Int,Int) -> Int -> Int -> a -> a
+checkMatVecMultAdd mn n m
+    | snd mn /= n =
+        error $ printf
+            ("Tried to multiply a matrix with shape `%s' by a vector of dimension `%d'")
+            (show mn) n
+    | fst mn /= m =
+        error $ printf
+            ("Tried to add a vector of dimension `%d' to a vector of dimension `%d'")
+            (fst mn) m
+    | otherwise = id
+
+checkMatMatMultAdd :: (Int,Int) -> (Int,Int) -> (Int,Int) -> a -> a
+checkMatMatMultAdd mk kn mn
+    | snd mk /= fst kn =
+        error $ printf
+            ("Tried to multiply a matrix with shape `%s' by a matrix with shape `%s'")
+            (show mk) (show kn)
+    | (fst mk, snd kn) /= mn =
+        error $ printf
+            ("Tried to add a matrix with shape `%s' to a matrix with shape `%s'")
+            (show (fst mk, snd kn)) (show mn)
+    | otherwise = id
+
+checkMatVecSolv :: (Int,Int) -> Int -> a -> a
+checkMatVecSolv mn m
+    | fst mn /= m =
+        error $ printf
+            ("Tried to solve a matrix with shape `%s' for a vector of dimension `%d'")
+            (show mn) m
+    | otherwise = id
+
+checkMatVecSolvTo :: (Int,Int) -> Int -> Int -> a -> a
+checkMatVecSolvTo mn m n
+    | fst mn /= m =
+        error $ printf
+            ("Tried to solve a matrix with shape `%s' for a vector of dimension `%d'")
+            (show mn) m
+    | snd mn /= n =
+        error $ printf
+            ("Tried to store a vector of dimension `%d' in a vector of dimension `%d'")
+            (show $ snd mn) n
+    | otherwise = id
+
+checkMatMatSolv :: (Int,Int) -> (Int,Int) -> a -> a
+checkMatMatSolv mn mk
+    | fst mn /= fst mk =
+        error $ printf
+            ("Tried to solve a matrix with shape `%s' for a matrix with shape `%s'")
+            (show mn) (show mk)
+    | otherwise = id
+
+checkMatMatSolvTo :: (Int,Int) -> (Int,Int) -> (Int,Int) -> a -> a
+checkMatMatSolvTo mk mn kn
+    | fst mn /= fst mk =
+        error $ printf
+            ("Tried to solve a matrix with shape `%s' for a matrix with shape `%s'")
+            (show mk) (show mn)
+    | kn /= (snd mk, snd mn) =
+        error $ printf
+            ("Tried to store a matrix with shape `%s' in a matrix with shape `%s'")
+            (show (snd mk, snd mn)) (show kn)
+    | otherwise = id
+
+checkSquare :: (Int,Int) -> a -> a
+checkSquare (m,n)
+    | m /= n =
+        error $ printf
+            ("Expected a square matrix but got one with shape `%s'")
+            (show (m,n))
+    | otherwise = id
+
+checkFat :: (Int,Int) -> a -> a
+checkFat (m,n)
+    | m > n =
+        error $ printf
+            ("Expected a fat matrix but got one with shape `%s'")
+            (show (m,n))
+    | otherwise = id
+
+checkTall :: (Int,Int) -> a -> a
+checkTall (m,n)
+    | m < n =
+        error $ printf
+            ("Expected a tall matrix but got one with shape `%s'")
+            (show (m,n))
+    | otherwise = id
diff --git a/BLAS/Matrix/Base.hs b/BLAS/Matrix/Base.hs
--- a/BLAS/Matrix/Base.hs
+++ b/BLAS/Matrix/Base.hs
@@ -1,3 +1,4 @@
+{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module     : BLAS.Matrix.Base
@@ -8,10 +9,14 @@
 --
 
 module BLAS.Matrix.Base (
-    Matrix(..)
+    Matrix(..),
+    isSquare,
+    isFat,
+    isTall,
     ) where
 
 import BLAS.Elem.Base ( Elem )
+import BLAS.Tensor.Base
 
 -- | A base class for matrices.
 class Matrix a where
@@ -24,3 +29,20 @@
     -- | Creates a new matrix view that conjugates and transposes the 
     -- given matrix.
     herm :: Elem e => a (m,n) e -> a (n,m) e
+
+
+isSquare :: (Matrix a) => a (m,n) e -> Bool
+isSquare a = numRows a == numCols a
+
+isFat :: (Matrix a) => a (m,n) e -> Bool
+isFat a = numRows a <= numCols a
+
+isTall :: (Matrix a) => a (m,n) e -> Bool
+isTall a = numRows a >= numCols a
+
+instance (Matrix a) => Tensor (a (m,n)) (Int,Int) e where
+    shape a = (numRows a, numCols a)
+    
+    bounds a = ((0,0), (m-1,n-1))
+      where (m,n) = shape a
+      
diff --git a/BLAS/Matrix/Immutable.hs b/BLAS/Matrix/Immutable.hs
--- a/BLAS/Matrix/Immutable.hs
+++ b/BLAS/Matrix/Immutable.hs
@@ -1,4 +1,5 @@
 {-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}
+{-# OPTIONS_GHC -fglasgow-exts #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module     : BLAS.Matrix.Immutable
@@ -14,21 +15,55 @@
 
 import BLAS.Access
 import BLAS.Elem ( BLAS3 )
-import qualified BLAS.Matrix.Base as Base
+import BLAS.Internal ( checkMatVecMult, checkMatMatMult )
+import BLAS.Matrix.ReadOnly
 import Data.Vector.Dense
 import Data.Matrix.Dense.Internal
-import Data.Matrix.Dense.Operations ( apply, applyMat )
 
-infixl 7 <*>, <**>
+import System.IO.Unsafe ( unsafePerformIO )
 
-class Base.Matrix a => IMatrix a e where
+infixr 7 <*>, <**>
+
+class (RMatrix a e) => IMatrix a e where
     -- | Apply to a vector
     (<*>) :: a (m,n) e -> Vector n e -> Vector m e
-    
+    (<*>) a x = 
+        checkMatVecMult (shape a) (dim x) $ unsafeApply a x
+        
     -- | Apply to a matrix
     (<**>) :: a (m,k) e -> Matrix (k,n) e -> Matrix (m,n) e
+    (<**>) a b = 
+        checkMatMatMult (numRows a, numCols a) 
+                        (numRows b, numCols b) $ unsafeApplyMat a b
 
-instance (BLAS3 e) => IMatrix (DMatrix Imm) e where
-    (<*>) = apply
-    (<**>)  = applyMat
+    sapply :: e -> a (m,n) e -> Vector n e -> Vector m e
+    sapply k a x =
+        checkMatVecMult (numRows a, numCols a) (dim x) $ unsafeSApply k a x
+        
+    sapplyMat :: e -> a (m,k) e -> Matrix (k,n) e -> Matrix (m,n) e    
+    sapplyMat k a b = 
+        checkMatMatMult (numRows a, numCols a) 
+                        (numRows b, numCols b) $ unsafeSApplyMat k a b
+
+
+    unsafeApply :: a (m,n) e -> Vector n e -> Vector m e
+    unsafeApply = unsafeSApply 1
     
+    unsafeApplyMat :: a (m,k) e -> Matrix (k,n) e -> Matrix (m,n) e
+    unsafeApplyMat = unsafeSApplyMat 1
+    
+    unsafeSApply :: e -> a (m,n) e -> Vector n e -> Vector m e
+    unsafeSApply alpha a x = unsafePerformIO $ unsafeGetSApply alpha a x
+    {-# NOINLINE unsafeSApply #-}
+
+    unsafeSApplyMat :: e -> a (m,k) e -> Matrix (k,n) e -> Matrix (m,n) e
+    unsafeSApplyMat alpha a b = unsafePerformIO $ unsafeGetSApplyMat alpha a b
+    {-# NOINLINE unsafeSApplyMat #-}
+
+
+instance (BLAS3 e) => IMatrix (DMatrix Imm) e where
+
+{-# RULES
+"scale.apply/sapply"       forall k a x. (<*>) (k *> a) x = sapply k a x
+"scale.applyMat/sapplyMat" forall k a b. (<**>) (k *> a) b = sapplyMat k a b
+  #-}
diff --git a/BLAS/Matrix/ReadOnly.hs b/BLAS/Matrix/ReadOnly.hs
--- a/BLAS/Matrix/ReadOnly.hs
+++ b/BLAS/Matrix/ReadOnly.hs
@@ -9,22 +9,251 @@
 --
 
 module BLAS.Matrix.ReadOnly (
-    RMatrix(..)
+    RMatrix(..),
+    
+    -- * Matrix and vector multiplication
+    getApply,
+    getApplyMat,
+    getSApply,
+    getSApplyMat,
+    
+    -- * In-place multiplication
+    doApply,
+    doApplyAdd,
+    doApplyMat,
+    doApplyAddMat,
+    doSApply,
+    doSApplyAdd,
+    doSApplyMat,
+    doSApplyAddMat,
+    doApply_,
+    doApplyMat_,
+    doSApply_,
+    doSApplyMat_,
+    
+    -- * Unsafe operations
+    unsafeGetApply,
+    unsafeGetApplyMat,
+    unsafeGetSApply,
+    unsafeGetSApplyMat,
+    
     ) where
 
-import BLAS.Elem ( BLAS3 )
+import BLAS.Elem ( BLAS1, BLAS3 )
+import BLAS.Internal ( checkMatVecMult, checkMatMatMult, checkMatVecMultAdd,
+    checkMatMatMultAdd, checkSquare )
 import qualified BLAS.Matrix.Base as Base
 import Data.Vector.Dense.Internal
+import Data.Vector.Dense.Operations
+import qualified Data.Vector.Dense.Operations as V
 import Data.Matrix.Dense.Internal
+import Data.Matrix.Dense.Operations
 import qualified Data.Matrix.Dense.Operations as M
 
-class Base.Matrix a => RMatrix a e where
-    -- | Apply to a vector
-    getApply :: a (m,n) e -> DVector t n e -> IO (DVector r m e)
+import Unsafe.Coerce
+
+-- | Minimal complete definition: (unsafeDoSApplyAdd, unsafeDoSApplyAddMat)
+-- or (unsafeDoSApply, unsafeDoSApplyMat)
+class (Base.Matrix a, BLAS1 e) => RMatrix a e where
+    unsafeDoApply :: a (m,n) e -> DVector t n e -> IOVector m e -> IO ()
+    unsafeDoApply =
+        unsafeDoSApply 1
+        
+    unsafeDoApplyAdd :: a (m,n) e -> DVector t n e -> e -> IOVector m e -> IO ()
+    unsafeDoApplyAdd a x beta y =
+        unsafeDoSApplyAdd 1 a x beta y
+
+    unsafeDoSApply :: e -> a (m,n) e -> DVector t n e -> IOVector m e -> IO ()
+    unsafeDoSApply alpha a x y =
+        unsafeDoSApplyAdd alpha a x 0 y
+
+    unsafeDoApplyMat :: a (m,k) e -> DMatrix t (k,n) e -> IOMatrix (m,n) e -> IO ()
+    unsafeDoApplyMat a b c =
+        unsafeDoSApplyMat 1 a b c
+        
+    unsafeDoApplyAddMat :: a (m,k) e -> DMatrix t (k,n) e -> e -> IOMatrix (m,n) e -> IO ()
+    unsafeDoApplyAddMat a b beta c =
+        unsafeDoSApplyAddMat 1 a b beta c
+
+    unsafeDoSApplyMat :: e -> a (m,k) e -> DMatrix t (k,n) e -> IOMatrix (m,n) e -> IO ()
+    unsafeDoSApplyMat alpha a b c =
+        unsafeDoSApplyAddMat alpha a b 0 c
     
-    -- | Apply to a matrix
-    getApplyMat :: a (m,k) e -> DMatrix t (k,n) e -> IO (DMatrix r (m,n) e)
+    unsafeDoApply_ :: a (n,n) e -> IOVector n e -> IO ()
+    unsafeDoApply_ a x =
+        unsafeDoSApply_ 1 a x
+        
+    unsafeDoApplyMat_ :: a (m,m) e -> IOMatrix (m,n) e -> IO ()
+    unsafeDoApplyMat_ a b =
+        unsafeDoSApplyMat_ 1 a b
 
+    unsafeDoSApply_ :: e -> a (n,n) e -> IOVector n e -> IO ()
+    unsafeDoSApply_ alpha a x = do
+        y <- newVector_ (dim x)
+        unsafeDoSApply alpha a x y
+        unsafeCopyVector x y
+
+    unsafeDoSApplyMat_ :: e -> a (m,m) e -> IOMatrix (m,n) e -> IO ()
+    unsafeDoSApplyMat_ alpha a b = do
+        c <- newMatrix_ (shape b)
+        unsafeDoSApplyMat alpha a b c
+        unsafeCopyMatrix b c
+        
+    unsafeDoSApplyAdd :: e -> a (m,n) e -> DVector t n e -> e -> IOVector m e -> IO ()
+    unsafeDoSApplyAdd alpha a x beta y = do
+        y' <- unsafeGetApply a x
+        V.scaleBy beta y
+        V.unsafeAxpy alpha y' y
+        
+    unsafeDoSApplyAddMat :: e -> a (m,k) e -> DMatrix t (k,n) e -> e -> IOMatrix (m,n) e -> IO ()
+    unsafeDoSApplyAddMat alpha a b beta c = do
+        c' <- unsafeGetApplyMat a b
+        M.scaleBy beta c
+        M.axpy alpha c' c
+
+
+unsafeGetApply :: (RMatrix a e) =>
+    a (m,n) e -> DVector t n e -> IO (DVector r m e)
+unsafeGetApply = unsafeGetSApply 1
+
+unsafeGetApplyMat :: (RMatrix a e) =>
+    a (m,k) e -> DMatrix t (k,n) e -> IO (DMatrix r (m,n) e)
+unsafeGetApplyMat = unsafeGetSApplyMat 1
+
+unsafeGetSApply :: (RMatrix a e) =>
+    e -> a (m,n) e -> DVector t n e -> IO (DVector r m e)
+unsafeGetSApply alpha a x = do
+    y <- newZero (numRows a)
+    unsafeDoSApply alpha a x y
+    return (unsafeCoerce y)
+
+unsafeGetSApplyMat :: (RMatrix a e) =>
+    e -> a (m,k) e -> DMatrix t (k,n) e -> IO (DMatrix r (m,n) e)
+unsafeGetSApplyMat alpha a b = do
+    c <- newZero (numRows a, numCols b)
+    unsafeDoSApplyMat alpha a b c
+    return (unsafeCoerce c)
+
+
+
+
+-- | Apply to a vector
+getApply :: (RMatrix a e) =>
+    a (m,n) e -> DVector t n e -> IO (DVector r m e)
+getApply a x =
+    checkMatVecMult (shape a) (dim x) $ do
+        unsafeGetApply a x
+
+-- | Apply to a matrix
+getApplyMat :: (RMatrix a e) =>
+    a (m,k) e -> DMatrix t (k,n) e -> IO (DMatrix r (m,n) e)
+getApplyMat a b =
+    checkMatMatMult (shape a) (shape b) $
+        unsafeGetApplyMat a b
+
+-- | Scale and apply to a vector
+getSApply :: (RMatrix a e) =>
+    e -> a (m,n) e -> DVector t n e -> IO (DVector r m e)
+getSApply k a x =
+    checkMatVecMult (shape a) (dim x) $ 
+        unsafeGetSApply k a x
+
+-- | Scale and apply to a matrix
+getSApplyMat :: (RMatrix a e) =>
+    e -> a (m,k) e -> DMatrix t (k,n) e -> IO (DMatrix r (m,n) e)
+getSApplyMat k a b =
+    checkMatMatMult (shape a) (shape b) $
+        unsafeGetSApplyMat k a b
+    
+-- | Apply to a vector and store the result in another vector
+doApply :: (RMatrix a e) =>
+    a (m,n) e -> DVector t n e -> IOVector m e -> IO ()
+doApply a x y =
+    checkMatVecMultAdd (numRows a, numCols a) (dim x) (dim y) $
+        unsafeDoApply a x y
+    
+-- | @y := a x + beta y@
+doApplyAdd :: (RMatrix a e) =>
+    a (m,n) e -> DVector t n e -> e -> IOVector m e -> IO ()
+doApplyAdd a x beta y =
+    checkMatVecMultAdd (shape a) (dim x) (dim y) $
+        unsafeDoApplyAdd a x beta y
+
+-- | @y := alpha a x@
+doSApply :: (RMatrix a e) =>
+    e -> a (m,n) e -> DVector t n e -> IOVector m e -> IO ()
+doSApply alpha a x y =
+    checkMatVecMultAdd (shape a) (dim x) (dim y) $
+        unsafeDoSApply alpha a x y
+
+-- | @y := alpha a x + beta y@    
+doSApplyAdd :: (RMatrix a e) =>
+    e -> a (m,n) e -> DVector t n e -> e -> IOVector m e -> IO ()
+doSApplyAdd alpha a x beta y =
+    checkMatVecMultAdd (shape a) (dim x) (dim y) $
+        unsafeDoSApplyAdd alpha a x beta y
+
+-- | Apply to a matrix and store the result in another matrix
+doApplyMat :: (RMatrix a e) =>
+    a (m,k) e -> DMatrix t (k,n) e -> IOMatrix (m,n) e -> IO ()
+doApplyMat a b c =
+    checkMatMatMultAdd (shape a) (shape b) (shape c) $
+        unsafeDoApplyMat a b c
+    
+-- | @c := a b + beta c@
+doApplyAddMat :: (RMatrix a e) =>
+    a (m,k) e -> DMatrix t (k,n) e -> e -> IOMatrix (m,n) e -> IO ()
+doApplyAddMat a b beta c =
+    checkMatMatMultAdd (shape a) (shape b) (shape c)
+        unsafeDoApplyAddMat a b beta c
+
+-- | Scale and apply to a matrix and store the result in another matrix
+doSApplyMat :: (RMatrix a e) => 
+    e -> a (m,k) e -> DMatrix t (k,n) e -> IOMatrix (m,n) e -> IO ()
+doSApplyMat alpha a b c =
+    checkMatMatMultAdd (shape a) (shape b) (shape c) $
+        unsafeDoSApplyMat alpha a b c
+        
+-- | @c := alpha a b + beta c@
+doSApplyAddMat :: (RMatrix a e) =>
+    e -> a (m,k) e -> DMatrix t (k,n) e -> e -> IOMatrix (m,n) e -> IO ()
+doSApplyAddMat alpha a b beta c =
+    checkMatMatMultAdd (shape a) (shape b) (shape c)
+        unsafeDoSApplyAddMat alpha a b beta c
+
+-- | @x := a x@    
+doApply_ :: (RMatrix a e) =>
+    a (n,n) e -> IOVector n e -> IO ()
+doApply_ a x =
+    checkSquare (shape a) $
+        checkMatVecMult (shape a) (dim x) $
+            unsafeDoApply_ a x
+
+-- | @ b := a b@
+doApplyMat_ :: (RMatrix a e) =>
+    a (m,m) e -> IOMatrix (m,n) e -> IO ()
+doApplyMat_ a b =
+    checkSquare (shape a) $
+        checkMatMatMult (shape a) (shape b) $
+            unsafeDoApplyMat_ a b
+
+-- | @ x := alpha a x@        
+doSApply_ :: (RMatrix a e) =>
+    e -> a (n,n) e -> IOVector n e -> IO ()
+doSApply_ alpha a x =
+    checkSquare (shape a) $
+        checkMatVecMult (shape a) (dim x) $
+            unsafeDoSApply_ alpha a x
+
+-- | @ b := alpha a b@
+doSApplyMat_ :: (RMatrix a e) =>
+    e -> a (m,m) e -> IOMatrix (m,n) e -> IO ()
+doSApplyMat_ alpha a b =
+    checkSquare (shape a) $
+        checkMatMatMult (shape a) (shape b) $
+            unsafeDoSApplyMat_ alpha a b
+
+
 instance (BLAS3 e) => RMatrix (DMatrix t) e where
-    getApply    = M.getApply
-    getApplyMat = M.getApplyMat
+    unsafeDoSApplyAdd    = M.gemv
+    unsafeDoSApplyAddMat = M.gemm
diff --git a/BLAS/Matrix/Solve/Immutable.hs b/BLAS/Matrix/Solve/Immutable.hs
--- a/BLAS/Matrix/Solve/Immutable.hs
+++ b/BLAS/Matrix/Solve/Immutable.hs
@@ -9,20 +9,67 @@
 --
 
 module BLAS.Matrix.Solve.Immutable (
-    ISolve(..)
+    ISolve(..),
+    (<\>),
+    (<\\>),
+    ssolve,
+    ssolveMat,
+    
     ) where
 
-import BLAS.Matrix.Immutable
+import BLAS.Internal ( checkMatVecSolv, checkMatMatSolv )
+import BLAS.Matrix.Solve.ReadOnly
 
-import Data.Vector.Dense ( Vector )
-import Data.Matrix.Dense ( Matrix )
+import Data.Vector.Dense ( Vector, dim )
+import Data.Matrix.Dense ( Matrix, shape )
 
-infixl 7 <\>, <\\>
+import System.IO.Unsafe ( unsafePerformIO )
 
-class IMatrix a e => ISolve a e where
-    -- | Solve for a vector
-    (<\>) :: a (m,n) e -> Vector m e -> Vector n e
+infixr 7 <\>, <\\>
+
+class RSolve a e => ISolve a e where
+    unsafeSolve :: a (m,n) e -> Vector m e -> Vector n e
+    unsafeSolve a y =
+        unsafePerformIO $ unsafeGetSolve a y
+    {-# NOINLINE unsafeSolve #-}
     
-    -- | Solve for a matrix
-    (<\\>) :: a (m,n) e -> Matrix (m,k) e -> Matrix (n,k) e
+    unsafeSolveMat :: a (m,n) e -> Matrix (m,k) e -> Matrix (n,k) e
+    unsafeSolveMat a c =
+        unsafePerformIO $ unsafeGetSolveMat a c
+    {-# NOINLINE unsafeSolveMat #-}
+
+    unsafeSSolve :: e -> a (m,n) e -> Vector m e -> Vector n e
+    unsafeSSolve alpha a y =
+        unsafePerformIO $ unsafeGetSSolve alpha a y
+    {-# NOINLINE unsafeSSolve #-}
+    
+    unsafeSSolveMat :: e -> a (m,n) e -> Matrix (m,k) e -> Matrix (n,k) e
+    unsafeSSolveMat alpha a c =
+        unsafePerformIO $ unsafeGetSSolveMat alpha a c
+    {-# NOINLINE unsafeSSolveMat #-}
+
+
+-- | Solve for a vector
+(<\>) :: (ISolve a e) => a (m,n) e -> Vector m e -> Vector n e
+(<\>) a y =
+    checkMatVecSolv (shape a) (dim y) $
+        unsafeSolve a y
+
+-- | Solve for a matrix
+(<\\>) :: (ISolve a e) => a (m,n) e -> Matrix (m,k) e -> Matrix (n,k) e
+(<\\>) a b =
+    checkMatMatSolv (shape a) (shape b) $
+        unsafeSolveMat a b
+
+-- | Solve for a vector and scale
+ssolve :: (ISolve a e) => e -> a (m,n) e -> Vector m e -> Vector n e
+ssolve alpha a y =
+    checkMatVecSolv (shape a) (dim y) $
+        unsafeSSolve alpha a y
+
+-- | Solve for a matrix and scale
+ssolveMat :: (ISolve a e) => e -> a (m,n) e -> Matrix (m,k) e -> Matrix (n,k) e
+ssolveMat alpha a b =
+    checkMatMatSolv (shape a) (shape b) $
+        unsafeSSolveMat alpha a b
 
diff --git a/BLAS/Matrix/Solve/ReadOnly.hs b/BLAS/Matrix/Solve/ReadOnly.hs
--- a/BLAS/Matrix/Solve/ReadOnly.hs
+++ b/BLAS/Matrix/Solve/ReadOnly.hs
@@ -9,17 +9,183 @@
 --
 
 module BLAS.Matrix.Solve.ReadOnly (
-    RSolve(..)
+    RSolve(..),
+    
+    -- * Matrix and vector solving
+    getSolve,
+    getSolveMat,
+    getSSolve,
+    getSSolveMat,
+    
+    -- * In-place solving
+    doSolve,
+    doSolveMat,
+    doSSolve,
+    doSSolveMat,
+    doSolve_,
+    doSolveMat_,
+    doSSolve_,
+    doSSolveMat_,
+    
+    -- * unsafe operations
+    unsafeGetSolve,
+    unsafeGetSolveMat,
+    unsafeGetSSolve,
+    unsafeGetSSolveMat,
+
     ) where
 
-import BLAS.Matrix.ReadOnly
+import BLAS.Internal ( checkMatVecSolv, checkMatMatSolv, checkMatVecSolvTo,
+    checkMatMatSolvTo, checkSquare )
 
-import Data.Vector.Dense.IO ( DVector )
-import Data.Matrix.Dense.IO ( DMatrix )
+import Data.Vector.Dense.IO ( DVector, IOVector, dim, newVector_ )
+import qualified Data.Vector.Dense.IO as V
+import Data.Matrix.Dense.IO ( DMatrix, IOMatrix, RMatrix, numCols, shape, newMatrix_ )
+import qualified Data.Matrix.Dense.IO as M
 
+import Unsafe.Coerce
+
+
 class RMatrix a e => RSolve a e where
-    -- | Solve for a vector
-    getSolve :: a (m,n) e -> DVector t m e -> IO (DVector r n e)
+    unsafeDoSolve :: a (m,n) e -> DVector t m e -> IOVector n e -> IO ()
+    unsafeDoSolve = unsafeDoSSolve 1
     
-    -- | Solve for a matrix
-    getSolveMat :: a (m,n) e -> DMatrix t (m,k) e -> IO (DMatrix r (n,k) e)
+    unsafeDoSolveMat :: a (m,n) e -> DMatrix t (m,k) e -> IOMatrix (n,k) e -> IO ()
+    unsafeDoSolveMat = unsafeDoSSolveMat 1
+    
+    unsafeDoSSolve :: e -> a (m,n) e -> DVector t m e -> IOVector n e -> IO ()
+    unsafeDoSSolve alpha a y x = do
+        unsafeDoSolve a y x
+        V.scaleBy alpha x
+    
+    unsafeDoSSolveMat :: e -> a (m,n) e -> DMatrix t (m,k) e -> IOMatrix (n,k) e -> IO ()
+    unsafeDoSSolveMat alpha a c b = do
+        unsafeDoSolveMat a c b
+        M.scaleBy alpha b
+
+    unsafeDoSolve_ :: a (n,n) e -> IOVector n e -> IO ()
+    unsafeDoSolve_ = unsafeDoSSolve_ 1
+
+    unsafeDoSSolve_ :: e -> a (n,n) e -> IOVector n e -> IO ()
+    unsafeDoSSolve_ alpha a x = do
+        V.scaleBy alpha x
+        unsafeDoSolve_ a x
+        
+    unsafeDoSolveMat_ :: a (m,m) e -> IOMatrix (m,n) e -> IO ()
+    unsafeDoSolveMat_ = unsafeDoSSolveMat_ 1
+        
+    unsafeDoSSolveMat_ :: e -> a (m,m) e -> IOMatrix (m,n) e -> IO ()
+    unsafeDoSSolveMat_ alpha a b = do
+        M.scaleBy alpha b
+        unsafeDoSolveMat_ a b
+
+
+unsafeGetSolve :: (RSolve a e) => 
+    a (m,n) e -> DVector t m e -> IO (DVector r n e)
+unsafeGetSolve a y = do
+    x <- newVector_ (numCols a)
+    unsafeDoSolve a y x
+    return (unsafeCoerce x)
+    
+unsafeGetSSolve :: (RSolve a e) => 
+    e -> a (m,n) e -> DVector t m e -> IO (DVector r n e)
+unsafeGetSSolve alpha a y = do
+    x <- newVector_ (numCols a)
+    unsafeDoSSolve alpha a y x
+    return (unsafeCoerce x)
+    
+unsafeGetSolveMat :: (RSolve a e) => 
+    a (m,n) e -> DMatrix t (m,k) e -> IO (DMatrix r (n,k) e)
+unsafeGetSolveMat a c = do
+    b <- newMatrix_ (numCols a, numCols c)
+    unsafeDoSolveMat a c b
+    return (unsafeCoerce b)
+
+unsafeGetSSolveMat :: (RSolve a e) => 
+    e -> a (m,n) e -> DMatrix t (m,k) e -> IO (DMatrix r (n,k) e)
+unsafeGetSSolveMat alpha a c = do
+    b <- newMatrix_ (numCols a, numCols c)
+    unsafeDoSSolveMat alpha a c b
+    return (unsafeCoerce b)
+
+-- | Solve for a vector
+getSolve :: (RSolve a e) => 
+    a (m,n) e -> DVector t m e -> IO (DVector r n e)
+getSolve a y = 
+    checkMatVecSolv (shape a) (dim y) $
+        unsafeGetSolve a y
+
+-- | Solve for a vector and scale
+getSSolve :: (RSolve a e) => 
+    e -> a (m,n) e -> DVector t m e -> IO (DVector r n e)
+getSSolve alpha a y = 
+    checkMatVecSolv (shape a) (dim y) $
+        unsafeGetSSolve alpha a y
+
+-- | Solve for a matrix
+getSolveMat :: (RSolve a e) => 
+    a (m,n) e -> DMatrix t (m,k) e -> IO (DMatrix r (n,k) e)
+getSolveMat a b =
+    checkMatMatSolv (shape a) (shape b) $
+            unsafeGetSolveMat a b
+            
+-- | Solve for a matrix and scale
+getSSolveMat :: (RSolve a e) => 
+    e -> a (m,n) e -> DMatrix t (m,k) e -> IO (DMatrix r (n,k) e)
+getSSolveMat alpha a b =
+    checkMatMatSolv (shape a) (shape b) $
+            unsafeGetSSolveMat alpha a b
+
+
+doSolve :: (RSolve a e) => 
+    a (m,n) e -> DVector t m e -> IOVector n e -> IO ()
+doSolve a y x =
+    checkMatVecSolvTo (shape a) (dim y) (dim x) $
+        unsafeDoSolve a y x
+        
+doSolveMat :: (RSolve a e) => 
+    a (m,n) e -> DMatrix t (m,k) e -> IOMatrix (n,k) e -> IO ()
+doSolveMat a c b =
+    checkMatMatSolvTo (shape a) (shape c) (shape b) $
+        unsafeDoSolveMat a c b
+    
+doSSolve :: (RSolve a e) => 
+    e -> a (m,n) e -> DVector t m e -> IOVector n e -> IO ()
+doSSolve alpha a y x =
+    checkMatVecSolvTo (shape a) (dim y) (dim x) $
+        unsafeDoSSolve alpha a y x
+
+doSSolveMat :: (RSolve a e) => 
+    e -> a (m,n) e -> DMatrix t (m,k) e -> IOMatrix (n,k) e -> IO ()
+doSSolveMat alpha a c b =
+    checkMatMatSolvTo (shape a) (shape c) (shape b) $
+        unsafeDoSSolveMat alpha a c b
+
+doSolve_ :: (RSolve a e) => 
+    a (n,n) e -> IOVector n e -> IO ()
+doSolve_ a x =
+    checkSquare (shape a) $
+        checkMatVecSolv (shape a) (dim x) $
+            unsafeDoSolve_ a x
+
+doSSolve_ :: (RSolve a e) => 
+    e -> a (n,n) e -> IOVector n e -> IO ()
+doSSolve_ alpha a x =
+    checkSquare (shape a) $
+        checkMatVecSolv (shape a) (dim x) $
+            unsafeDoSSolve_ alpha a x
+
+doSolveMat_ :: (RSolve a e) => 
+    a (m,m) e -> IOMatrix (m,n) e -> IO ()
+doSolveMat_ a b =
+    checkSquare (shape a) $
+        checkMatMatSolv (shape a) (shape b) $
+            unsafeDoSolveMat_ a b
+
+doSSolveMat_ :: (RSolve a e) => 
+    e -> a (m,m) e -> IOMatrix (m,n) e -> IO ()
+doSSolveMat_ alpha a b =
+    checkSquare (shape a) $
+        checkMatMatSolv (shape a) (shape b) $
+            unsafeDoSSolveMat_ alpha a b
+
diff --git a/BLAS/Tensor/Dense/Immutable.hs b/BLAS/Tensor/Dense/Immutable.hs
--- a/BLAS/Tensor/Dense/Immutable.hs
+++ b/BLAS/Tensor/Dense/Immutable.hs
@@ -21,3 +21,8 @@
     
     -- | Get a new constant tensor of the given shape.
     constant :: i -> e -> x e
+
+    -- | Apply a function to pairs of elements of tensors that are the 
+    -- same shape.
+    azipWith :: (ITensor x i f, ITensor x i g) => (e -> f -> g) -> x e -> x f -> x g
+    
diff --git a/BLAS/Tensor/Immutable.hs b/BLAS/Tensor/Immutable.hs
--- a/BLAS/Tensor/Immutable.hs
+++ b/BLAS/Tensor/Immutable.hs
@@ -49,10 +49,6 @@
     -- | Apply a function elementwise to a tensor.
     amap :: (ITensor x i e') => (e -> e') -> x e -> x e'
     
-    -- | Apply a function to pairs of elements of tensors that are the 
-    -- same shape.
-    azipWith :: (ITensor x i f, ITensor x i g) => (e -> f -> g) -> x e -> x f -> x g
-    
     -- ixmap :: i -> (i -> i) -> x e -> x e
     -- unsafeIxMap
 
diff --git a/Data/Matrix/Banded.hs b/Data/Matrix/Banded.hs
new file mode 100644
--- /dev/null
+++ b/Data/Matrix/Banded.hs
@@ -0,0 +1,69 @@
+{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module     : Data.Matrix.Banded
+-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+
+module Data.Matrix.Banded (
+    -- * Banded matrix type
+    Banded,
+    
+    module BLAS.Matrix.Base,
+    module BLAS.Tensor.Base,
+    module BLAS.Tensor.Dense.Immutable,
+    module BLAS.Tensor.Immutable,
+    module BLAS.Tensor.Scalable,
+    module Data.Matrix.Banded.Operations,
+
+    -- * Creating banded matrices
+    banded, 
+    listsBanded,
+    
+    -- * Properties
+    numLower,
+    numUpper,
+    bandwidth,
+    
+    -- * Rows and columns
+    row,
+    col,
+    rows,
+    cols,
+
+    -- * Diagonals
+    diag,
+    toLists,
+
+    -- * Casting matrices
+    coerceBanded,
+    
+    -- * Converting between vectors and matrices
+    
+    -- * Unsafe operations
+    unsafeBanded,
+    unsafeRow,
+    unsafeCol,
+    unsafeDiag,
+    
+    ) where
+
+import BLAS.Access
+import BLAS.Elem ( BLAS1 )
+import BLAS.Matrix.Base hiding ( Matrix )
+import BLAS.Tensor.Base
+import BLAS.Tensor.Dense.Immutable
+import BLAS.Tensor.Immutable
+import BLAS.Tensor.Scalable
+
+import Data.Matrix.Banded.Internal
+import qualified Data.Matrix.Banded.Internal as B
+import Data.Matrix.Banded.Operations hiding ( RMatrix(..), getScaled, 
+    getInvScaled, doConj, scaleBy, invScaleBy )
+import Data.Vector.Dense hiding ( scale, invScale )
+
+instance (BLAS1 e) => Scalable (BMatrix Imm (m,n)) e where
+    (*>) = scale
diff --git a/Data/Matrix/Banded/IO.hs b/Data/Matrix/Banded/IO.hs
new file mode 100644
--- /dev/null
+++ b/Data/Matrix/Banded/IO.hs
@@ -0,0 +1,72 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module     : Data.Matrix.Banded.IO
+-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+
+module Data.Matrix.Banded.IO (
+    -- * The mutable banded matrix data type
+    BMatrix(..),
+    IOBanded,
+    
+    module BLAS.Matrix.Base,
+    module BLAS.Tensor.Base,
+    module BLAS.Tensor.Dense.ReadOnly,
+    module BLAS.Tensor.ReadOnly,
+    module BLAS.Tensor.Mutable,
+    
+    -- * Creating new matrices
+    newBanded,
+    newBanded_,
+    newListsBanded,
+
+    -- * Views
+    -- ** Rows and columns
+    rowView,
+    colView,
+    getRow,
+    getCol,
+    
+    -- ** Diagonals
+    diag,
+    
+    -- * Operations
+    module Data.Matrix.Banded.Operations,
+    
+    -- * Converting to and from banded matrices
+    -- ** @ForeignPtr@s
+    toForeignPtr,
+    fromForeignPtr,
+
+    -- * Bandwith properties
+    bandwidth,
+    numLower,
+    numUpper,
+    
+    -- * Coercing
+    coerceBanded,
+    
+    -- * Unsafe operations
+    unsafeNewBanded,
+    unsafeWithElemPtr,
+    unsafeRowView,
+    unsafeColView,
+    unsafeGetRow,
+    unsafeGetCol,
+    unsafeDiag,
+    unsafeFreeze,
+    unsafeThaw,
+    
+    ) where
+
+import Data.Matrix.Banded.Internal
+import Data.Matrix.Banded.Operations hiding ( IMatrix(..), scale, invScale )
+    
+import BLAS.Matrix.Base hiding ( Matrix )
+import BLAS.Tensor.Base
+import BLAS.Tensor.Dense.ReadOnly
+import BLAS.Tensor.ReadOnly
+import BLAS.Tensor.Mutable
diff --git a/Data/Matrix/Banded/Internal.hs b/Data/Matrix/Banded/Internal.hs
--- a/Data/Matrix/Banded/Internal.hs
+++ b/Data/Matrix/Banded/Internal.hs
@@ -1,4 +1,5 @@
 {-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}
+{-# OPTIONS_GHC -fglasgow-exts #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module     : Data.Matrix.Banded.Internal
@@ -20,12 +21,10 @@
     -- * Converting to and from foreign pointers
     toForeignPtr,
     fromForeignPtr,
-    ldaOf,
-    isHerm,
     
     -- * To and from the underlying storage matrix
-    toMatrix,
-    fromMatrix,
+    toRawMatrix,
+    fromRawMatrix,
     
     -- * Bandwith properties
     bandwidth,
@@ -44,8 +43,11 @@
     -- * Getting rows and columns
     row,
     col,
+    rows,
+    cols,
     getRow,
     getCol,
+    toLists,
     
     -- * Vector views
     diag,
@@ -53,7 +55,7 @@
     colView,
     
     -- * Casting matrices
-    coerceMatrix,
+    coerceBanded,
     
     -- * Unsafe operations
     unsafeBanded,
@@ -61,6 +63,7 @@
     unsafeFreeze,
     unsafeThaw,
     unsafeWithElemPtr,
+    unsafeWithBasePtr,
     unsafeDiag,
     unsafeGetRow,
     unsafeGetCol,
@@ -74,6 +77,7 @@
 import Control.Monad ( zipWithM_ )
 import Data.Ix ( inRange, range )
 import Data.List ( foldl' )
+import Data.Maybe ( fromJust )
 import Foreign
 import System.IO.Unsafe
 import Unsafe.Coerce          
@@ -88,6 +92,8 @@
 import qualified BLAS.Matrix.Base as C
 import BLAS.Tensor
 
+import Data.AEq
+
 import Data.Matrix.Dense.Internal ( DMatrix )
 import qualified Data.Matrix.Dense.Internal as M
                                 
@@ -96,34 +102,25 @@
 
         
 data BMatrix t mn e 
-    = BM { fptr   :: !(ForeignPtr e)
-         , offset :: !Int
-         , size1  :: !Int
-         , size2  :: !Int
-         , lowBW  :: !Int
-         , upBW   :: !Int
-         , lda    :: !Int
+    = BM { fptrOf   :: {-# UNPACK #-} !(ForeignPtr e)
+         , offsetOf :: {-# UNPACK #-} !Int
+         , size1    :: {-# UNPACK #-} !Int
+         , size2    :: {-# UNPACK #-} !Int
+         , lowBW    :: {-# UNPACK #-} !Int
+         , upBW     :: {-# UNPACK #-} !Int
+         , ldaOf    :: {-# UNPACK #-} !Int
+         , isHerm   :: {-# UNPACK #-} !Bool
          }
-    | H !(BMatrix t mn e) 
 
 type Banded = BMatrix Imm
 type IOBanded = BMatrix Mut
 
-fromForeignPtr :: ForeignPtr e -> Int -> (Int,Int) -> (Int,Int) -> Int 
+fromForeignPtr :: ForeignPtr e -> Int -> (Int,Int) -> (Int,Int) -> Int -> Bool
     -> BMatrix t (m,n) e
-fromForeignPtr f o (m,n) (kl,ku) l = BM f o m n kl ku l
-
-toForeignPtr :: BMatrix t (m,n) e -> (ForeignPtr e, Int, (Int,Int), (Int,Int), Int)
-toForeignPtr (H a) = toForeignPtr a
-toForeignPtr (BM f o m n kl ku l) = (f, o, (m,n), (kl,ku), l)
-
-ldaOf :: BMatrix t (m,n) e -> Int
-ldaOf (H a) = ldaOf a
-ldaOf a     = lda a
+fromForeignPtr f o (m,n) (kl,ku) l h = BM f o m n kl ku l h
 
-isHerm :: BMatrix t (m,n) e -> Bool
-isHerm (H a) = not (isHerm a)
-isHerm _     = False
+toForeignPtr :: BMatrix t (m,n) e -> (ForeignPtr e, Int, (Int,Int), (Int,Int), Int, Bool)
+toForeignPtr (BM f o m n kl ku l h) = (f, o, (m,n), (kl,ku), l, h)
 
 unsafeFreeze :: BMatrix t mn e -> Banded mn e
 unsafeFreeze = unsafeCoerce
@@ -132,27 +129,38 @@
 unsafeThaw = unsafeCoerce
 
 -- | Coerce the phantom shape type from one type to another.
-coerceMatrix :: BMatrix t mn e -> BMatrix t kl e
-coerceMatrix = unsafeCoerce
+coerceBanded :: BMatrix t mn e -> BMatrix t kl e
+coerceBanded = unsafeCoerce
 
-toMatrix :: (Elem e) => BMatrix t (m,n) e -> (DMatrix t (m',n') e, (Int,Int), (Int,Int))
-toMatrix (H a) = 
-    case toMatrix (herm a) of (b, (m,n), (kl,ku)) -> (herm b, (n,m), (ku,kl))
-toMatrix (BM f o m n kl ku ld) = 
-    (M.fromForeignPtr f o (kl+1+ku,n) ld, (m,n), (kl,ku))
+toLists :: (BLAS1 e) => Banded (m,n) e -> ((Int,Int), (Int,Int),[[e]])
+toLists a = ( (m,n)
+            , (kl,ku)
+            , map paddedDiag [(-kl)..ku]
+            )
+  where
+    (m,n)   = shape a
+    (kl,ku) = (numLower a, numUpper a)
+    
+    padBegin i = replicate (max (-i) 0)    0
+    padEnd   i = replicate (max (m-n+i) 0) 0
+    paddedDiag i = (padBegin i) ++ (elems $ diag a i) ++ (padEnd i)
 
-fromMatrix :: (Elem e) => DMatrix t (m,n) e -> (Int,Int) -> (Int,Int) -> BMatrix t (m',n') e
-fromMatrix a (m,n) (kl,ku) = case a of
-    (M.H a') -> herm (fromMatrix a' (n,m) (ku,kl))
-    _        -> 
-        let (f,o,(m',n'),ld) = M.toForeignPtr a
-        in case undefined of
-            _ | m' /= kl+1+ku -> 
-                error $ "fromMatrix: number of rows must be equal to number of diagonals"
-            _ | n' /= n ->
-                error $ "fromMatrix: numbers of columns must be equal"
-            _ ->
-                BM f o m n kl ku ld
+toRawMatrix :: (Elem e) => BMatrix t (m,n) e -> ((Int,Int), (Int,Int), DMatrix t (m',n') e, Bool)
+toRawMatrix (BM f o m n kl ku ld h) = 
+    ((m,n), (kl,ku), M.fromForeignPtr f o (kl+1+ku,n) ld False, h)
+
+fromRawMatrix :: (Elem e) => (Int,Int) -> (Int,Int) -> DMatrix t (m,n) e -> Bool -> Maybe (BMatrix t (m',n') e)
+fromRawMatrix (m,n) (kl,ku) a h = 
+    if M.isHerm a 
+        then Nothing
+        else let (f,o,(m',n'),ld,_) = M.toForeignPtr a
+             in case undefined of
+                 _ | m' /= kl+1+ku -> 
+                     error $ "fromMatrix: number of rows must be equal to number of diagonals"
+                 _ | n' /= n ->
+                     error $ "fromMatrix: numbers of columns must be equal"
+                 _ ->
+                     Just $ BM f o m n kl ku ld h
                 
 
 bandwidth :: BMatrix t (m,n) e -> (Int,Int)
@@ -161,12 +169,12 @@
     in (negate kl, ku)
 
 numLower :: BMatrix t (m,n) e -> Int
-numLower (H a) = numUpper a
-numLower a     = lowBW a
+numLower a | isHerm a  = upBW a
+           | otherwise = lowBW a
 
 numUpper :: BMatrix t (m,n) e -> Int
-numUpper (H a) = numLower a
-numUpper a     = upBW a
+numUpper a | isHerm a  = lowBW a
+           | otherwise = upBW a
 
 
 newBanded_ :: (Elem e) => (Int,Int) -> (Int,Int) -> IO (BMatrix t (m,n) e)
@@ -185,9 +193,10 @@
         let off = 0
             m'  = kl + 1 + ku
             l   = m'
+            h   = False
         in do    
             ptr <- mallocForeignPtrArray (m' * n)
-            return $ fromForeignPtr ptr off (m,n) (kl,ku) l
+            return $ fromForeignPtr ptr off (m,n) (kl,ku) l h
     where
       err s = ioError $ userError $ 
                   "newBanded_ " ++ show (m,n) ++ " " ++ show (kl,ku) ++ ": " ++ s
@@ -211,7 +220,7 @@
     -> (Int,Int) -> (Int,Int) -> [((Int,Int),e)] -> IO (BMatrix t (m,n) e)
 newBandedHelp set (m,n) (kl,ku) ijes = do
     x <- newBanded_ (m,n) (kl,ku)
-    withForeignPtr (fptr x) $ flip clearArray ((kl+1+ku)*n)
+    withForeignPtr (fptrOf x) $ flip clearArray ((kl+1+ku)*n)
     mapM_ (uncurry $ set $ unsafeThaw x) ijes
     return x
 
@@ -232,33 +241,47 @@
         in zipWithM_ (unsafeWriteElem d) [0..(dim d - 1)] es'
 
 unsafeDiag :: (Elem e) => BMatrix t (m,n) e -> Int -> DVector t k e
-unsafeDiag (H a) d = 
-    conj $ unsafeDiag a (negate d)
-unsafeDiag a d =
-    let f      = fptr a
-        off    = indexOf a (diagStart d)
-        len    = diagLen (shape a) d
-        stride = lda a
-    in V.fromForeignPtr f off len stride
+unsafeDiag a d
+    | isHerm a = conj $ unsafeDiag (herm a) (negate d)
+    | otherwise =
+        let f      = fptrOf a
+            off    = indexOf a (diagStart d)
+            len    = diagLen (shape a) d
+            stride = ldaOf a
+            c      = False
+        in V.fromForeignPtr f off len stride c
         
 diag :: (Elem e) => BMatrix t (m,n) e -> Int -> DVector t k e
 diag a = checkedDiag (shape a) (unsafeDiag a) 
 
 
-indexOf :: BMatrix t mn e -> (Int,Int) -> Int
-indexOf (H a) (i,j) = indexOf a (j,i)
-indexOf (BM _ off _ _ _ ku ld) (i,j) =
-    off + ku + (i - j) + j * ld
-    --off + i * tda + (j - i + kl)
+indexOf :: BMatrix t (m,n) e -> (Int,Int) -> Int
+indexOf (BM _ off _ _ _ ku ld h) (i,j) =
+    let (i',j') = if h then (j,i) else (i,j)
+    in off + ku + (i' - j') + j' * ld
+    --off + i' * tda + (j' - i' + kl)
            
 
+hasStorage :: BMatrix t (m,n) e -> (Int,Int) -> Bool
+hasStorage (BM _ _ m n kl ku _ h) (i,j) =
+    let (i',j') = if h then (j,i) else (i,j)
+    in (  inRange (0,m-1) i'
+       && inRange (0,n-1) j'
+       && inRange (max 0 (j'-ku), min (m-1) (j'+kl)) i'
+       )
+              
             
 unsafeWithElemPtr :: (Elem e) => BMatrix t (m,n) e -> (Int,Int) -> (Ptr e -> IO a) -> IO a
-unsafeWithElemPtr a (i,j) f = case a of
-    (H a') -> unsafeWithElemPtr a' (j,i) f
-    _      -> withForeignPtr (fptr a) $ \ptr ->
-                  f $ ptr `advancePtr` (indexOf a (i,j))
+unsafeWithElemPtr a (i,j) f
+    | isHerm a  = unsafeWithElemPtr (herm a) (j,i) f
+    | otherwise = withForeignPtr (fptrOf a) $ \ptr ->
+                      f $ ptr `advancePtr` (indexOf a (i,j))
 
+unsafeWithBasePtr :: (Elem e) => BMatrix t (m,n) e -> (Ptr e -> IO a) -> IO a
+unsafeWithBasePtr a f =
+    withForeignPtr (fptrOf a) $ \ptr ->
+        f $ ptr `advancePtr` (offsetOf a)
+
 row :: (BLAS1 e) => Banded (m,n) e -> Int -> Vector n e
 row a = checkedRow (shape a) (unsafeRow a)
 
@@ -290,38 +313,44 @@
 unsafeGetCol :: (BLAS1 e) => BMatrix t (m,n) e -> Int -> IO (DVector r m e)
 unsafeGetCol a j = unsafeGetRow (herm a) j >>= return . conj
 
+rows :: (BLAS1 e) => Banded (m,n) e -> [Vector n e]
+rows a = [ unsafeRow a i | i <- [0..(numRows a - 1)] ]
 
+cols :: (BLAS1 e) => Banded (m,n) e -> [Vector m e]
+cols a = [ unsafeCol a i | i <- [0..(numCols a - 1)] ]
+
 unsafeColView :: (Elem e) => BMatrix t (m,n) e -> Int -> (Int, DVector t k e, Int)
-unsafeColView (BM f off m _ kl ku ld) j =
-    let nb     = max (j - ku)         0
-        na     = max (m - 1 - j - kl) 0
-        r      = max (ku - j) 0 
-        c      = j 
-        off'   = off + r + c * ld
-        stride = 1
-        len    = m - (nb + na)
-    in if len >= 0
-        then (nb, V.fromForeignPtr f off' len stride, na)
-        else (m , V.fromForeignPtr f off' 0   stride,  0)
-      
-unsafeColView a j = 
-    case unsafeRowView (herm a) j of (nb, v, na) -> (nb, conj v, na)
+unsafeColView a@(BM f off m _ kl ku ld _) j
+    | isHerm a = 
+        case unsafeRowView (herm a) j of (nb, v, na) -> (nb, conj v, na)
+    | otherwise =
+        let nb     = max (j - ku)         0
+            na     = max (m - 1 - j - kl) 0
+            r      = max (ku - j) 0 
+            c      = j 
+            off'   = off + r + c * ld
+            stride = 1
+            len    = m - (nb + na)
+        in if len >= 0
+            then (nb, V.fromForeignPtr f off' len stride False, na)
+            else (m , V.fromForeignPtr f off' 0   stride False,  0)
 
-unsafeRowView :: (Elem e) => BMatrix t (m,n) e -> Int -> (Int, DVector t k e, Int)
-unsafeRowView (BM f off _ n kl ku ld) i =
-    let nb     = max (i - kl)         0
-        na     = max (n - 1 - i - ku) 0
-        r      = min (ku + i)         (kl + ku)
-        c      = max (i - kl)         0 
-        off'   = off + r + c * ld
-        stride = ld - 1
-        len    = n - (nb + na)
-    in if len >= 0 
-        then (nb, V.fromForeignPtr f off' len stride, na)
-        else (n , V.fromForeignPtr f off' 0   stride,  0)
 
-unsafeRowView a i =
-    case unsafeColView (herm a) i of (nb, v, na) -> (nb, conj v, na)
+unsafeRowView :: (Elem e) => BMatrix t (m,n) e -> Int -> (Int, DVector t k e, Int)
+unsafeRowView a@(BM f off _ n kl ku ld _) i
+    | isHerm a =
+        case unsafeColView (herm a) i of (nb, v, na) -> (nb, conj v, na)        
+    | otherwise =
+        let nb     = max (i - kl)         0
+            na     = max (n - 1 - i - ku) 0
+            r      = min (ku + i)         (kl + ku)
+            c      = max (i - kl)         0 
+            off'   = off + r + c * ld
+            stride = ld - 1
+            len    = n - (nb + na)
+        in if len >= 0 
+            then (nb, V.fromForeignPtr f off' len stride False, na)
+            else (n , V.fromForeignPtr f off' 0   stride False,  0)
 
 
 rowView :: (Elem e) => BMatrix t (m,n) e -> Int -> (Int, DVector t k e, Int)
@@ -332,19 +361,15 @@
 
 
 instance C.Matrix (BMatrix t) where
-    numRows = fst . shape
-    numCols = snd . shape
+    numRows a | isHerm a  = size2 a
+              | otherwise = size1 a
+              
+    numCols a | isHerm a  = size1 a
+              | otherwise = size2 a
 
-    herm a = case a of
-        (H a')   -> coerceMatrix a'
-        _        -> H (coerceMatrix a)
-        
-instance Tensor (BMatrix t (m,n)) (Int,Int) e where
-    shape a = case a of
-        (H a')   -> case shape a' of (m,n) -> (n,m)
-        _        -> (size1 a, size2 a)
+    herm a = let h' = (not . isHerm) a
+             in coerceBanded $ a{ isHerm=h' }
 
-    bounds a = let (m,n) = shape a in ((0,0), (m-1,n-1))
     
 instance (BLAS1 e) => ITensor (BMatrix Imm (m,n)) (Int,Int) e where
     size = inlinePerformIO . getSize
@@ -358,7 +383,7 @@
     (//)          = replaceHelp writeElem
     unsafeReplace = replaceHelp unsafeWriteElem
 
-    amap f a = banded (shape a) (bandwidth a) ies
+    amap f a = banded (shape a) (numLower a, numUpper a) ies
       where
         ies = map (second f) (assocs a)
 
@@ -376,20 +401,23 @@
     
 instance (BLAS1 e) => RTensor (BMatrix t (m,n)) (Int,Int) e IO where
     newCopy b = 
-        let (a,mn,kl) = toMatrix b
+        let (mn,kl,a,h) = toRawMatrix b
         in do
             a' <- newCopy a
-            return $ fromMatrix a' mn kl
+            return $ fromJust $ fromRawMatrix mn kl a' h
     
-    getSize a = case a of
-        (H a')               -> getSize a'
-        (BM _ _ m n kl ku _) -> return $ foldl' (+) 0 $ 
-                                         map (diagLen (m,n)) [(-kl)..ku]
+    getSize (BM _ _ m n kl ku _ _) = 
+        return $ foldl' (+) 0 $ 
+             map (diagLen (m,n)) [(-kl)..ku]
     
-    unsafeReadElem a (i,j) = case a of
-        (H a') -> unsafeReadElem a' (j,i) >>= return . E.conj
-        _      -> withForeignPtr (fptr a) $ \ptr ->
-                      peekElemOff ptr (indexOf a (i,j))
+    unsafeReadElem a (i,j)
+        | isHerm a = 
+            unsafeReadElem (herm a) (j,i) >>= return . E.conj
+        | hasStorage a (i,j) =
+            withForeignPtr (fptrOf a) $ \ptr ->
+                peekElemOff ptr (indexOf a (i,j))
+        | otherwise =
+            return 0
 
     getIndices a =
         return $ filter (\ij -> inlinePerformIO $ canModifyElem (unsafeThaw a) ij)
@@ -402,18 +430,48 @@
         mapM (\i -> unsafeReadElem a i >>= \e -> return (i,e)) is
 
 instance (BLAS1 e) => MTensor (BMatrix Mut (m,n)) (Int,Int) e IO where
-    setZero a = case toMatrix a of (a',_,_) -> setZero a'
-    setConstant e a = case toMatrix a of (a',_,_) -> setConstant e a'
-        
-    canModifyElem a (i,j) = case a of
-        (H a')               -> canModifyElem a' (j,i)
-        (BM _ _ m n kl ku _) -> return $ inRange (0,m-1) i && 
-                                         inRange (0,n-1) j &&
-                                         inRange (max 0 (j-ku), min (m-1) (j+kl)) i
+    setZero a = case toRawMatrix a of (_,_,a',_) -> setZero a'
     
-    unsafeWriteElem a (i,j) e = case a of
-        (H a') -> unsafeWriteElem a' (j,i) (E.conj e)
-        _      -> withForeignPtr (fptr a) $ \ptr ->
-                      pokeElemOff ptr (indexOf a (i,j)) e
+    setConstant e a = 
+        let (_,_,a',h) = toRawMatrix a
+            e' = if h then E.conj e else e
+        in setConstant e' a'
+        
+    canModifyElem a ij = return $ hasStorage a ij
 
-    modifyWith f a = case toMatrix a of (a',_,_) -> modifyWith f a'
+    unsafeWriteElem a (i,j) e
+        | isHerm a = 
+            unsafeWriteElem (herm a) (j,i) (E.conj e)
+        | otherwise =
+            withForeignPtr (fptrOf a) $ \ptr ->
+                pokeElemOff ptr (indexOf a (i,j)) e
+
+    modifyWith f a = 
+        let (_,_,a',h) = toRawMatrix a
+        in if h then modifyWith f a'
+                else modifyWith f (herm a)
+
+instance (BLAS1 e) => Show (BMatrix Imm (m,n) e) where
+    show a 
+        | isHerm a = 
+           "herm (" ++ show (herm a) ++ ")"
+        | otherwise = 
+             let (mn,kl,es) = toLists a 
+             in "listsBanded " ++ show mn ++ " " ++ show kl ++ " " ++ show es
+       
+       
+compareHelp :: (BLAS1 e) => 
+    (e -> e -> Bool) -> Banded (m,n) e -> Banded (m,n) e -> Bool
+compareHelp cmp x y
+    | isHerm x && isHerm y =
+        compareHelp cmp (herm x) (herm y)
+compareHelp cmp x y =
+    (shape x == shape y) && (and $ zipWith cmp (elems x) (elems y))
+
+instance (BLAS1 e, Eq e) => Eq (BMatrix Imm (m,n) e) where
+    (==) = compareHelp (==)
+
+instance (BLAS1 e, AEq e) => AEq (BMatrix Imm (m,n) e) where
+    (===) = compareHelp (===)
+    (~==) = compareHelp (~==)       
+             
diff --git a/Data/Matrix/Banded/Operations.hs b/Data/Matrix/Banded/Operations.hs
new file mode 100644
--- /dev/null
+++ b/Data/Matrix/Banded/Operations.hs
@@ -0,0 +1,148 @@
+{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module     : Data.Matrix.Banded.Operations
+-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+
+module Data.Matrix.Banded.Operations (
+    module BLAS.Matrix.Immutable,
+    module BLAS.Matrix.ReadOnly,
+    
+    -- * Matrix Arithmetic
+    -- ** Pure
+    scale,
+    invScale,
+    
+    -- ** Impure
+    getScaled,
+    getInvScaled,
+    
+    -- * In-place operations
+    doConj,
+    scaleBy,
+    invScaleBy,
+
+    ) where
+
+import System.IO.Unsafe
+import Unsafe.Coerce
+
+import Data.Matrix.Banded.Internal
+import Data.Matrix.Dense.Internal ( DMatrix, IOMatrix)
+import Data.Vector.Dense.Internal hiding ( unsafeWithElemPtr, unsafeThaw, 
+    unsafeFreeze )
+import qualified Data.Vector.Dense.Operations as V
+import qualified Data.Vector.Dense.Internal as V
+import qualified Data.Matrix.Dense.Internal as M
+import qualified Data.Matrix.Dense.Operations as M
+
+import BLAS.Access
+import BLAS.C ( CBLASTrans, colMajor, noTrans, conjTrans )
+import qualified BLAS.C as BLAS
+import BLAS.Elem ( BLAS1, BLAS2  )
+import qualified BLAS.Elem as E
+import BLAS.Matrix.Immutable
+import BLAS.Matrix.ReadOnly
+
+infixl 7 `scale`, `invScale`
+
+
+-- | Form a new matrix by multiplying every element by a value.
+getScaled :: (BLAS1 e) => e -> BMatrix t (m,n) e -> IO (BMatrix r (m,n) e)
+getScaled k = unaryOp (scaleBy k)
+
+-- | Form a new matrix by dividing every element by a value.
+getInvScaled :: (BLAS1 e) => e -> BMatrix t (m,n) e -> IO (BMatrix r (m,n) e)
+getInvScaled k = unaryOp (invScaleBy k)
+
+-- | Conjugate every element in a matrix.
+doConj  :: (BLAS1 e) => IOBanded (m,n) e -> IO ()
+doConj a = let (_,_,a',_) = toRawMatrix a
+           in M.doConj a'
+
+-- | Scale every element in a matrix by the given value.
+scaleBy :: (BLAS1 e) => e -> IOBanded (m,n) e -> IO ()
+scaleBy k a = 
+    let (_,_,a',h) = toRawMatrix a
+        k' = if h then E.conj k else k
+    in M.scaleBy k' a'
+    
+-- | Divide every element by the given value.
+invScaleBy :: (BLAS1 e) => e -> IOBanded (m,n) e -> IO ()
+invScaleBy k a = 
+    let (_,_,a',h) = toRawMatrix a
+        k' = if h then E.conj k else k
+    in M.invScaleBy k' a'
+
+blasTransOf :: BMatrix t (m,n) e -> CBLASTrans
+blasTransOf a = 
+    case (isHerm a) of
+          False -> noTrans
+          True  -> conjTrans
+
+flipShape :: (Int,Int) -> (Int,Int)
+flipShape (m,n) = (n,m)
+
+
+-- | @gbmv alpha a x beta y@ replaces @y := alpha a * x + beta y@
+gbmv :: (BLAS2 e) => e -> BMatrix s (m,n) e -> DVector t n e -> e -> IOVector m e -> IO ()
+gbmv alpha a x beta y
+    | numRows a == 0 || numCols a == 0 =
+        return ()
+    | isConj x = do
+        x' <- V.getConj (conj x)
+        gbmv alpha a x' beta y
+    | isConj y = do
+        V.doConj y
+        gbmv alpha a x beta (conj y)
+        V.doConj y
+    | otherwise =
+        let order  = colMajor
+            transA = blasTransOf a
+            (m,n)  = case (isHerm a) of
+                         False -> shape a
+                         True  -> (flipShape . shape) a
+            (kl,ku) = case (isHerm a) of
+                          False -> (numLower a, numUpper a)
+                          True  -> (numUpper a, numLower a)
+            ldA    = ldaOf a
+            incX   = V.strideOf x
+            incY   = V.strideOf y
+        in unsafeWithBasePtr a $ \pA ->
+               V.unsafeWithElemPtr x 0 $ \pX ->
+                    V.unsafeWithElemPtr y 0 $ \pY -> do
+                        BLAS.gbmv order transA m n kl ku alpha pA ldA pX incX beta pY incY
+
+-- | @gbmm alpha a b beta c@ replaces @c := alpha a * b + beta c@.
+gbmm :: (BLAS2 e) => e -> BMatrix s (m,k) e -> DMatrix t (k,n) e -> e -> IOMatrix (m,n) e -> IO ()
+gbmm alpha a b beta c =
+    sequence_ $
+        zipWith (\x y -> gbmv alpha a x beta y) (M.cols b) (M.cols c)
+
+unaryOp :: (BLAS1 e) => (IOBanded (m,n) e -> IO ()) 
+    -> BMatrix t (m,n) e -> IO (BMatrix r (m,n) e)
+unaryOp f a = do
+    a' <- newCopy a
+    f (unsafeThaw a')
+    return (unsafeCoerce a')
+
+-- | Create a new matrix by scaling another matrix by the given value.
+scale :: (BLAS1 e) => e -> Banded (m,n) e -> Banded (m,n) e
+scale k a = unsafePerformIO $ getScaled k a
+{-# NOINLINE scale #-}
+
+-- | Form a new matrix by dividing every element by a value.
+invScale :: (BLAS1 e) => e -> Banded (m,n) e -> Banded (m,n) e
+invScale k a = unsafePerformIO $ getInvScaled k a
+{-# NOINLINE invScale #-}
+
+instance (BLAS2 e) => RMatrix (BMatrix s) e where
+    unsafeDoSApplyAdd    = gbmv
+    unsafeDoSApplyAddMat = gbmm
+
+instance (BLAS2 e) => IMatrix (BMatrix Imm) e
+
diff --git a/Data/Matrix/Dense.hs b/Data/Matrix/Dense.hs
--- a/Data/Matrix/Dense.hs
+++ b/Data/Matrix/Dense.hs
@@ -40,12 +40,6 @@
     -- * Augmenting matrices
     submatrix,
 
-    -- * Matrix multiplication
-    apply,
-    applyMat,
-    sapply,
-    sapplyMat,
-
     -- * Matrix arithmetic
     shift,
     scale,
@@ -81,8 +75,8 @@
 
 import Data.Matrix.Dense.Internal
 import qualified Data.Matrix.Dense.Internal as M
-import Data.Matrix.Dense.Operations ( apply, applyMat, sapply, sapplyMat,
-    shift, scale, invScale, plus, minus, times, divide )
+import Data.Matrix.Dense.Operations ( shift, scale, invScale, plus, minus, 
+    times, divide )
 import Data.Vector.Dense hiding ( scale, invScale, shift )
 
 
diff --git a/Data/Matrix/Dense/IO.hs b/Data/Matrix/Dense/IO.hs
--- a/Data/Matrix/Dense/IO.hs
+++ b/Data/Matrix/Dense/IO.hs
@@ -59,8 +59,6 @@
     -- ** @ForeignPtr@s
     toForeignPtr,
     fromForeignPtr,
-    ldaOf,
-    isHerm,
     
     -- ** Coercing
     coerceMatrix,
@@ -78,8 +76,8 @@
     ) where
 
 import Data.Matrix.Dense.Internal
-import Data.Matrix.Dense.Operations hiding ( gemv, gemm, apply, applyMat,
-    sapply, sapplyMat, add, plus, minus, times, divide, getApply, getApplyMat )
+import Data.Matrix.Dense.Operations hiding ( add, plus, minus, times, divide, 
+    gemv, gemm )
     
 import BLAS.Matrix.Base hiding ( Matrix )
 import BLAS.Matrix.ReadOnly
diff --git a/Data/Matrix/Dense/Internal.hs b/Data/Matrix/Dense/Internal.hs
--- a/Data/Matrix/Dense/Internal.hs
+++ b/Data/Matrix/Dense/Internal.hs
@@ -1,4 +1,5 @@
 {-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}
+{-# OPTIONS_GHC -fglasgow-exts #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module     : Data.Matrix.Dense.Internal
@@ -20,8 +21,6 @@
     -- * Converting to and from foreign pointers
     toForeignPtr,
     fromForeignPtr,
-    ldaOf,
-    isHerm,
 
     -- * Creating new matrices
     newMatrix,
@@ -80,7 +79,7 @@
 import Data.AEq
 
 import Data.Vector.Dense.Internal hiding ( toForeignPtr, fromForeignPtr,
-    unsafeFreeze, unsafeThaw, fptr, offset, unsafeWithElemPtr )
+    unsafeFreeze, unsafeThaw, storageOf, offsetOf, unsafeWithElemPtr )
 import qualified Data.Vector.Dense.Internal as V
 import qualified Data.Vector.Dense.Operations as V
 
@@ -99,13 +98,13 @@
 -- column-major order, or provide a view into another matrix.  The view 
 -- transposes and conjugates the underlying matrix.
 data DMatrix t mn e =
-      DM { fptr   :: !(ForeignPtr e) -- ^ a pointer to the storage region
-         , offset :: !Int            -- ^ an offset (in elements, not bytes) to the first element in the matrix. 
-         , size1  :: !Int            -- ^ the number of rows in the matrix
-         , size2  :: !Int            -- ^ the number of columns in the matrix
-         , lda    :: !Int            -- ^ the leading dimension size of the matrix
+      DM { storageOf :: {-# UNPACK #-} !(ForeignPtr e) -- ^ a pointer to the storage region
+         , offsetOf  :: {-# UNPACK #-} !Int            -- ^ an offset (in elements, not bytes) to the first element in the matrix. 
+         , size1     :: {-# UNPACK #-} !Int            -- ^ the number of rows in the matrix
+         , size2     :: {-# UNPACK #-} !Int            -- ^ the number of columns in the matrix
+         , ldaOf     :: {-# UNPACK #-} !Int            -- ^ the leading dimension size of the matrix
+         , isHerm    :: {-# UNPACK #-} !Bool           -- ^ indicates whether or not the matrix is transposed and conjugated
          }
-    | H !(DMatrix t mn e)           -- ^ a transposed and conjugated matrix
 
 type Matrix = DMatrix Imm
 type IOMatrix = DMatrix Mut
@@ -120,42 +119,26 @@
 coerceMatrix :: DMatrix t mn e -> DMatrix t kl e
 coerceMatrix = unsafeCoerce
 
--- | @fromForeignPtr f o mn l@ creates a matrix view of the data pointed to
--- by @f@ starting at offset @o@ and having shape @mn@ and lda @l@.
-fromForeignPtr :: ForeignPtr e -> Int -> (Int,Int) -> Int -> DMatrix t (m,n) e
-fromForeignPtr f o (m,n) l = DM f o m n l
-
--- | Convert a dense matrix to a pointer, offset, size, and lda.  Note that this
--- does not give the conjugacy/transpose information.  For that, use 'isHerm'.
-toForeignPtr :: DMatrix t (m,n) e -> (ForeignPtr e, Int, (Int,Int), Int)
-toForeignPtr   (H a)            = toForeignPtr a
-toForeignPtr a@(DM _ _ _ _ _)   = (fptr a, offset a, (size1 a, size2 a), lda a)
-
--- | Get the lda of a matrix, defined as the number of elements in the underlying
--- array that separate two consecutive elements in the same row of the matrix.
-ldaOf :: DMatrix t (m,n) e -> Int
-ldaOf   (H a)          = ldaOf a
-ldaOf a@(DM _ _ _ _ _) = lda a
+-- | @fromForeignPtr f o mn l h@ creates a matrix view of the data pointed to
+-- by @f@ starting at offset @o@ and having shape @mn@ and lda @l@.  If @h@
+-- is @True@ the matrix is interpreted as transposed and conjugated.
+fromForeignPtr :: ForeignPtr e -> Int -> (Int,Int) -> Int -> Bool -> DMatrix t (m,n) e
+fromForeignPtr f o (m,n) l h = DM f o m n l h
 
+-- | Convert a dense matrix to a pointer, offset, size, lda, herm.
+toForeignPtr :: DMatrix t (m,n) e -> (ForeignPtr e, Int, (Int,Int), Int, Bool)
+toForeignPtr a = (storageOf a, offsetOf a, (size1 a, size2 a), ldaOf a, isHerm a)
 
 indexOf :: DMatrix t (m,n) e -> (Int,Int) -> Int
-indexOf   (H a)          (i,j) = indexOf a (j,i)
-indexOf a@(DM _ _ _ _ _) (i,j) = 
-    let o = offset a
-        l = lda a
-    in o + i + j*l
+indexOf a (i,j) = 
+    let (i',j') = case isHerm a of
+                        True  -> (j,i)
+                        False -> (i,j)
+        o = offsetOf a
+        l = ldaOf a
+    in o + i' + j'*l
+{-# INLINE indexOf #-}
     
--- | Get the storage order of the matrix.  If 'isTrans' is true, this 
--- will be 'RowMajor'.  Otherwise, it will be 'ColMajor'.
-orderOf :: DMatrix t (m,n) e -> Order
-orderOf (H a)          = flipOrder (orderOf a)
-orderOf (DM _ _ _ _ _) = ColMajor
-
--- | Get whether or not the matrix is transposed and conjugated.
-isHerm :: DMatrix t (m,n) e -> Bool
-isHerm (H a)          = not (isHerm a)
-isHerm (DM _ _ _ _ _) = False
-
 -- | Create a new matrix of the given size and initialize the given elements to
 -- the given values.  All other elements get initialized to zero.
 newMatrix :: (BLAS1 e) => (Int,Int) -> [((Int,Int), e)] -> IO (DMatrix t (m,n) e)
@@ -181,7 +164,7 @@
             "Tried to create a matrix with shape `" ++ show (m,n) ++ "'"
     | otherwise = do
         f <- mallocForeignPtrArray (m*n)
-        return $ fromForeignPtr f 0 (m,n) (max 1 m)
+        return $ fromForeignPtr f 0 (m,n) (max 1 m) False
 
 -- | Create a new matrix with the given elements in column-major order.
 newListMatrix :: (Elem e) => (Int,Int) -> [e] -> IO (DMatrix t (m,n) e)
@@ -232,9 +215,8 @@
 -- at the given index.  It may be necessary to conjugate or scale values before
 -- reading or writing to or from the location.
 unsafeWithElemPtr :: (Elem e) => DMatrix t (m,n) e -> (Int,Int) -> (Ptr e -> IO a) -> IO a
-unsafeWithElemPtr   (H a) (i,j) f = unsafeWithElemPtr a (j,i) f
-unsafeWithElemPtr a@(DM _ _ _ _ _) ij f =
-    withForeignPtr (fptr a) $ \ptr ->
+unsafeWithElemPtr a ij f =
+    withForeignPtr (storageOf a) $ \ptr ->
         let ptr' = ptr `advancePtr` (indexOf a ij)
         in f ptr'
         
@@ -256,23 +238,29 @@
 
 -- | Same as 'row', but does not do any bounds checking.
 unsafeRow ::  (Elem e) => DMatrix t (m,n) e -> Int -> DVector t n e
-unsafeRow a@(H _)          i = conj $ unsafeCol (herm a) i
-unsafeRow a@(DM _ _ _ _ _) i =
-    let f = fptr a
-        o = indexOf a (i,0)
-        n = size2 a
-        s = lda a
-    in V.fromForeignPtr f o n s
+unsafeRow a i
+    | isHerm a =
+        conj $ unsafeCol (herm a) i
+    | otherwise =
+        let f = storageOf a
+            o = indexOf a (i,0)
+            n = numCols a
+            s = ldaOf a
+            c = False
+        in V.fromForeignPtr f o n s c
 
 -- | Same as 'col', but does not do any bounds checking.
 unsafeCol :: (Elem e) => DMatrix t (m,n) e -> Int -> DVector t m e
-unsafeCol a@(H _)          j = conj $ unsafeRow (herm a) j
-unsafeCol a@(DM _ _ _ _ _) j =
-    let f = fptr a
-        o = indexOf a (0,j)
-        m = size1 a
-        s = 1
-    in V.fromForeignPtr f o m s
+unsafeCol a j 
+    | isHerm a =
+        conj $ unsafeRow (herm a) j
+    | otherwise =
+        let f = storageOf a
+            o = indexOf a (0,j)
+            m = numRows a
+            s = 1
+            c = False
+        in V.fromForeignPtr f o m s c
 
 -- | @diag a 0@ gets a vector view of the main diagonal of @a@.  @diag a k@ for 
 -- @k@ positive gets a view of the @k@th superdiagonal.  For @k@ negative, it
@@ -282,13 +270,16 @@
 
 -- | Same as 'diag', but does not do any bounds checking.
 unsafeDiag :: (Elem e) => DMatrix t (m,n) e -> Int -> DVector t k e
-unsafeDiag   (H a)          i = conj $ unsafeDiag a (negate i)
-unsafeDiag a@(DM _ _ _ _ _) i =
-    let f = fptr a
-        o = indexOf a (diagStart i)
-        n = diagLen (shape a) i
-        s = lda a + 1
-    in V.fromForeignPtr f o n s
+unsafeDiag a i 
+    | isHerm a = 
+        conj $ unsafeDiag (herm a) (negate i)
+    | otherwise =            
+        let f = storageOf a
+            o = indexOf a (diagStart i)
+            n = diagLen (shape a) i
+            s = ldaOf a + 1
+            c = False
+        in V.fromForeignPtr f o n s c
 
 -- | @submatrix a ij mn@ returns a view of the submatrix of @a@ with element @(0,0)@
 -- being element @ij@ in @a@, and having shape @mn@.
@@ -297,73 +288,74 @@
 
 -- | Same as 'submatrix' but does not do any bounds checking.
 unsafeSubmatrix :: (Elem e) => DMatrix t (m,n) e -> (Int,Int) -> (Int,Int) -> DMatrix t (k,l) e
-unsafeSubmatrix a@(H _)          (i,j) (m',n') = herm $ unsafeSubmatrix (herm a) (j,i) (n',m')
-unsafeSubmatrix a@(DM _ _ _ _ _) (i,j) mn' =
-    let f = fptr a
-        o = indexOf a (i,j)
-        l = lda a
-    in fromForeignPtr f o mn' l
+unsafeSubmatrix a (i,j) (m,n)
+    | isHerm a  = 
+        herm $ unsafeSubmatrix (herm a) (j,i) (n,m)
+    | otherwise =
+        let f = storageOf a
+            o = indexOf a (i,j)
+            l = ldaOf a
+        in fromForeignPtr f o (m,n) l False
     
 
 -- | Create a matrix view of a row vector.  This will fail if the
--- stride is not @1@ and the vector is conjugated.
+-- vector is conjugated and the stride is not @1@.
 maybeFromRow :: (Elem e) => DVector t m e -> Maybe (DMatrix t (one,m) e)
-maybeFromRow (V.C   (V.C x))   = maybeFromRow x
-maybeFromRow (V.C x@(V.DV _ _ _ _))
-    | V.stride x == 1 =
-        let f = V.fptr x
-            o = V.offset x
-            n = V.dim x
+maybeFromRow x
+    | isConj x && strideOf x == 1 =
+        let f = V.storageOf x
+            o = V.offsetOf x
+            n = dim x
             l = max 1 n
-        in Just $ herm $ fromForeignPtr f o (n,1) l
+            h = True
+        in Just $ fromForeignPtr f o (n,1) l h
+    | (not . isConj) x =
+        let f = V.storageOf x
+            o = V.offsetOf x
+            n = dim x
+            s = strideOf x
+            l = max 1 s
+            h = False
+        in Just $ fromForeignPtr f o (1,n) l h
     | otherwise =
         Nothing
-maybeFromRow x@(V.DV _ _ _ _) =
-    let f = V.fptr x
-        o = V.offset x
-        n = V.dim x
-        s = V.stride x
-        l = max 1 s
-    in Just $ fromForeignPtr f o (1,n) l
 
 
 -- | Possibly create a matrix view of a column vector.  This will fail
 -- if the stride of the vector is not @1@ and the vector is not conjugated.
 maybeFromCol :: (Elem e) => DVector t n e -> Maybe (DMatrix t (n,one) e)
-maybeFromCol   (V.C x)   = maybeFromRow x >>= return . herm
-maybeFromCol x@(V.DV _ _ _ _)
-    | V.stride x == 1 =
-        let f = V.fptr x
-            o = V.offset x
+maybeFromCol x
+    | isConj x = maybeFromRow (conj x) >>= return . herm
+    | strideOf x == 1 =
+        let f = V.storageOf x
+            o = V.offsetOf x
             m = dim x
             l = max 1 m
-        in Just $ fromForeignPtr f o (m,1) l
+        in Just $ fromForeignPtr f o (m,1) l False
     | otherwise =
         Nothing
 
 maybeToVector :: (Elem e) => DMatrix t (m,n) e -> Maybe (Order, DVector t k e)
-maybeToVector (H a)   = maybeToVector a >>= (\(o,x) -> return (flipOrder o, conj x))
-maybeToVector (DM f o m n ld)
+maybeToVector a@(DM f o m n ld h) 
+    | h = 
+        maybeToVector (herm a) >>= (\(ord,x) -> return (flipOrder ord, conj x))
     | ld == m =
-        Just $ (ColMajor, V.fromForeignPtr f o (m*n) 1)
+        Just $ (ColMajor, V.fromForeignPtr f o (m*n) 1  False)
     | m == 1 =
-        Just $ (ColMajor, V.fromForeignPtr f o n    ld)
+        Just $ (ColMajor, V.fromForeignPtr f o n     ld False)
     | otherwise =
         Nothing
 
--- | Modify each element in-place by applying a function to it.
--- modifyWith :: (Elem e) => (e -> e) -> IOMatrix (m,n) e -> IO ()
 
-
 -- | Take a unary elementwise vector operation and apply it to the elements of a matrix.
 liftV :: (Elem e) => (DVector t k e -> IO ()) -> DMatrix t (m,n) e -> IO ()
 liftV f a =
     case maybeToVector a of
         Just (_,x) -> f x
         _ -> 
-            let xs  = case orderOf a of
-                          RowMajor -> rows (coerceMatrix a)
-                          ColMajor -> cols (coerceMatrix a)
+            let xs  = case isHerm a of
+                          True  -> rows (coerceMatrix a)
+                          False -> cols (coerceMatrix a)
             in mapM_ f xs
 
 -- | Take a binary elementwise vector operation and apply it to the elements of a pair
@@ -375,27 +367,22 @@
         (Just (RowMajor,x), Just (RowMajor,y)) -> f x y
         (Just (ColMajor,x), Just (ColMajor,y)) -> f x y
         _ -> 
-            let (xs,ys) = case orderOf a of
-                               RowMajor -> (rows (coerceMatrix a), rows (coerceMatrix b))
-                               ColMajor -> (cols (coerceMatrix a), cols (coerceMatrix b))
+            let (xs,ys) = case isHerm a of
+                               True  -> (rows (coerceMatrix a), rows (coerceMatrix b))
+                               False -> (cols (coerceMatrix a), cols (coerceMatrix b))
             in zipWithM_ f xs ys
 
 
 instance C.Matrix (DMatrix t) where
-    numRows = fst . shape
-    numCols = snd . shape
+    numRows a | isHerm a  = size2 a
+              | otherwise = size1 a
 
-    herm a = case a of
-        (H a')   -> coerceMatrix a'
-        _        -> H (coerceMatrix a)
-    
-instance Tensor (DMatrix t (m,n)) (Int,Int) e where
-    shape a = case a of
-        (H a')   -> case shape a' of (m,n) -> (n,m)
-        _        -> (size1 a, size2 a)
-    
-    bounds a = let (m,n) = shape a in ((0,0), (m-1,n-1))
+    numCols a | isHerm a  = size1 a
+              | otherwise = size2 a
 
+    herm a = a{ isHerm=(not . isHerm) a }
+
+
 instance (BLAS1 e) => ITensor (DMatrix Imm (m,n)) (Int,Int) e where
     size a = (numRows a * numCols a)
     
@@ -417,17 +404,7 @@
 
     amap f a = listMatrix (shape a) (map f $ elems a)
 
-    azipWith f a b
-        | shape b /= mn =
-            error ("azipWith: matrix shapes differ; first has shape `"
-                   ++ show mn ++ "' and second has shape `" 
-                   ++ show (shape b) ++ "'")
-        | otherwise =
-            listMatrix mn (zipWith f (elems a) (elems b))
-        where
-            mn = shape a
 
-
 replaceHelp :: (BLAS1 e) => 
        (IOMatrix (m,n) e -> (Int,Int) -> e -> IO ())
     -> Matrix (m,n) e -> [((Int,Int), e)] -> Matrix (m,n) e
@@ -446,21 +423,31 @@
     constant mn = unsafePerformIO . newConstant mn
     {-# NOINLINE constant #-}
      
+    azipWith f a b
+        | shape b /= mn =
+            error ("azipWith: matrix shapes differ; first has shape `"
+                   ++ show mn ++ "' and second has shape `" 
+                   ++ show (shape b) ++ "'")
+        | otherwise =
+            listMatrix mn (zipWith f (elems a) (elems b))
+        where
+            mn = shape a
      
+     
 instance (BLAS1 e) => RTensor (DMatrix t (m,n)) (Int,Int) e IO where
     getSize a = return (numRows a * numCols a)
     
-    newCopy a = case a of
-        (H a')   -> newCopy a' >>= return . H
-        _        -> do
-            a' <- newMatrix_ (shape a)
-            liftV2 V.copyVector (unsafeThaw a') a
-            return a'
-    
-    unsafeReadElem a (i,j) = case a of
-        (H a')   -> unsafeReadElem a' (j,i) >>= return . E.conj
-        _        -> withForeignPtr (fptr a) $ \ptr ->
-                        peekElemOff ptr (indexOf a (i,j))
+    newCopy a | isHerm a  = 
+                  newCopy (herm a) >>= return . herm
+              | otherwise = do
+                  a' <- newMatrix_ (shape a)
+                  liftV2 V.copyVector (unsafeThaw a') a
+                  return a'
+           
+    unsafeReadElem a (i,j)
+        | isHerm a  = unsafeReadElem (herm a) (j,i) >>= return . E.conj
+        | otherwise = withForeignPtr (storageOf a) $ \ptr ->
+                          peekElemOff ptr (indexOf a (i,j))
     {-# INLINE unsafeReadElem #-}
     
     getIndices = return . indices . unsafeFreeze
@@ -506,18 +493,19 @@
         return $ inRange (bounds a) ij
     {-# INLINE canModifyElem #-}
 
-    unsafeWriteElem a (i,j) e = case a of
-        (H a')   -> unsafeWriteElem a' (j,i) $ E.conj e
-        _        -> withForeignPtr (fptr a) $ \ptr ->
-                        pokeElemOff ptr (indexOf a (i,j)) e
+    unsafeWriteElem a (i,j) e
+        | isHerm a  = unsafeWriteElem (herm a) (j,i) $ E.conj e
+        | otherwise = withForeignPtr (storageOf a) $ \ptr ->
+                          pokeElemOff ptr (indexOf a (i,j)) e
     
     modifyWith f = liftV (modifyWith f)       
 
 
 instance (BLAS1 e, Show e) => Show (DMatrix Imm (m,n) e) where
-    show a = case a of
-        (H a')   -> "herm (" ++ show a' ++ ")"
-        _        -> "listMatrix " ++ show (shape a) ++ " " ++ show (elems a)
+    show a | isHerm a = 
+                "herm (" ++ show (herm a) ++ ")"
+           | otherwise =
+                "listMatrix " ++ show (shape a) ++ " " ++ show (elems a)
         
 compareHelp :: (BLAS1 e) => 
     (e -> e -> Bool) -> Matrix (m,n) e -> Matrix (m,n) e -> Bool
diff --git a/Data/Matrix/Dense/Operations.hs b/Data/Matrix/Dense/Operations.hs
--- a/Data/Matrix/Dense/Operations.hs
+++ b/Data/Matrix/Dense/Operations.hs
@@ -13,19 +13,6 @@
     copyMatrix,
     swapMatrices,
 
-    -- * Matrix multiplication
-    -- ** Pure
-    apply,
-    applyMat,
-    sapply,
-    sapplyMat,
-    
-    -- ** Impure
-    getApply,
-    getApplyMat,
-    getSApply,
-    getSApplyMat,
-
     -- * Matrix Arithmetic
     -- ** Pure
     shift,
@@ -61,12 +48,17 @@
     gemv,
     gemm,
     
+    -- * Unsafe operations
+    unsafeCopyMatrix,
+    unsafeSwapMatrices,
+    
     ) where
 
+import Data.Maybe ( fromJust )
 import System.IO.Unsafe
 import Unsafe.Coerce
 
-import BLAS.Internal ( checkMatMatOp, checkMatVecMult, checkMatMatMult )
+import BLAS.Internal ( checkMatMatOp )
 import Data.Matrix.Dense.Internal
 import Data.Vector.Dense.Internal hiding ( unsafeWithElemPtr, unsafeThaw, 
     unsafeFreeze )
@@ -78,52 +70,26 @@
 import BLAS.Elem ( BLAS1, BLAS2, BLAS3 )
 import qualified BLAS.Elem as E
 
-infixl 7 `apply`, `applyMat`, `scale`, `invScale`
+infixl 7 `scale`, `invScale`
 infixl 6 `shift`
 infixl 1 +=, -=, *=, //=
 
 -- | @copy dst src@ copies the elements from the second argument to the first.
 copyMatrix :: (BLAS1 e) => IOMatrix (m,n) e -> DMatrix t (m,n) e -> IO ()
-copyMatrix a b = checkMatMatOp "copyMatrix" (shape a) (shape b) >> unsafeCopyMatrix a b
+copyMatrix a b = 
+    checkMatMatOp "copyMatrix" (shape a) (shape b) $ unsafeCopyMatrix a b
 
 unsafeCopyMatrix :: (BLAS1 e) => IOMatrix (m,n) e -> DMatrix t (m,n) e -> IO ()
 unsafeCopyMatrix = liftV2 (V.unsafeCopyVector)
 
 -- | @swap a b@ exchanges the elements stored in two matrices.
 swapMatrices :: (BLAS1 e) => IOMatrix (m,n) e -> IOMatrix (m,n) e -> IO ()
-swapMatrices a b = checkMatMatOp "swapMatrices" (shape a) (shape b) >> unsafeSwapMatrices a b
+swapMatrices a b = 
+    checkMatMatOp "swapMatrices" (shape a) (shape b) $ unsafeSwapMatrices a b
 
 unsafeSwapMatrices :: (BLAS1 e) => IOMatrix (m,n) e -> IOMatrix (m,n) e -> IO ()
 unsafeSwapMatrices = liftV2 (V.unsafeSwapVectors)
 
--- | Multiply a matrix by a vector.
-getApply :: (BLAS2 e) => DMatrix s (m,n) e -> DVector t n e -> IO (DVector r m e)
-getApply = getSApply 1
-
--- | Multiply a scaled matrix by a vector.
-getSApply :: (BLAS2 e) => e -> DMatrix s (m,n) e -> DVector t n e -> IO (DVector r m e)
-getSApply alpha a x = checkMatVecMult (shape a) (V.dim x) >> unsafeGetSApply alpha a x
-
-unsafeGetSApply :: (BLAS2 e) => e -> DMatrix s (m,n) e -> DVector t n e -> IO (DVector r m e)
-unsafeGetSApply alpha a x = do
-    y <- V.newZero (numRows a)
-    gemv alpha a x 0 y
-    return (unsafeCoerce y)
-
--- | Multiply a matrix by a matrix.
-getApplyMat :: (BLAS3 e) => DMatrix s (m,k) e -> DMatrix t (k,n) e -> IO (DMatrix r (m,n) e)
-getApplyMat = getSApplyMat 1
-
--- | Multiply a scaled matrix by a matrix.
-getSApplyMat :: (BLAS3 e) => e -> DMatrix s (m,k) e -> DMatrix t (k,n) e -> IO (DMatrix r (m,n) e)
-getSApplyMat alpha a b = checkMatMatMult (shape a) (shape b) >> unsafeGetSApplyMat alpha a b
-
-unsafeGetSApplyMat :: (BLAS3 e) => e -> DMatrix s (m,k) e -> DMatrix t (k,n) e -> IO (DMatrix r (m,n) e)
-unsafeGetSApplyMat alpha a b = do
-    c <- newZero (numRows a, numCols b)
-    gemm alpha a b 0 c
-    return (unsafeCoerce c)
-
 -- | Form a new matrix by adding a value to every element in a matrix.
 getShifted :: (BLAS1 e) => e -> DMatrix t (m,n) e -> IO (DMatrix r (m,n) e)
 getShifted k = unaryOp (shiftBy k)
@@ -138,21 +104,22 @@
 
 -- | Create a new matrix by taking the elementwise sum of two matrices.
 getSum :: (BLAS1 e) => e -> DMatrix s (m,n) e -> e -> DMatrix t (m,n) e -> IO (DMatrix r (m,n) e)
-getSum alpha a beta b = checkMatMatOp "getSum" (shape a) (shape b) 
-    >> unsafeGetSum alpha a beta b
+getSum alpha a beta b = 
+    checkMatMatOp "getSum" (shape a) (shape b) $ unsafeGetSum alpha a beta b
 
 unsafeGetSum :: (BLAS1 e) => e -> DMatrix s (m,n) e -> e -> DMatrix t (m,n) e -> IO (DMatrix r (m,n) e)
-unsafeGetSum alpha a@(H _) beta b = do
-    s <- unsafeGetSum (E.conj alpha) (herm a) (E.conj beta) (herm b)
-    return (herm s)
-unsafeGetSum alpha a@(DM _ _ _ _ _) beta b = do
-    s <- getScaled alpha a
-    axpy beta b (unsafeThaw s)
-    return (unsafeCoerce s)
+unsafeGetSum alpha a beta b 
+    | isHerm a = do
+        s <- unsafeGetSum (E.conj alpha) (herm a) (E.conj beta) (herm b)
+        return (herm s)
+    | otherwise = do
+        s <- getScaled alpha a
+        axpy beta b (unsafeThaw s)
+        return (unsafeCoerce s)
 
 -- | Create a new matrix by taking the elementwise difference of two matrices.
 getDiff :: (BLAS1 e) => DMatrix s (m,n) e -> DMatrix t (m,n) e -> IO (DMatrix r (m,n) e)
-getDiff a b = checkMatMatOp "minus" (shape a) (shape b) >> unsafeGetSum 1 a (-1) b
+getDiff a b = checkMatMatOp "minus" (shape a) (shape b) $ unsafeGetSum 1 a (-1) b
 
 -- | Create a new matrix by taking the elementwise product of two matrices.
 getProduct :: (BLAS2 e) => DMatrix s (m,n) e -> DMatrix t (m,n) e -> IO (DMatrix r (m,n) e)
@@ -163,13 +130,9 @@
 getRatio = binaryOp "getRatio" (//=)
 
 -- | Conjugate every element in a matrix.
-doConj  :: (BLAS1 e) => IOMatrix (m,n) e -> IO (IOMatrix (m,n) e)
-doConj (H a) = do
-    a' <- doConj a
-    return (H a')
-doConj a@(DM _ _ _ _ _) = do
-    liftV (\x -> V.doConj x >> return ()) a
-    return a
+doConj  :: (BLAS1 e) => IOMatrix (m,n) e -> IO ()
+doConj a | isHerm a  = doConj (herm a)
+         | otherwise = liftV (\x -> V.doConj x >> return ()) a
 
 -- | Scale every element in a matrix by the given value.
 scaleBy :: (BLAS1 e) => e -> IOMatrix (m,n) e -> IO ()
@@ -211,19 +174,27 @@
 flipShape :: (Int,Int) -> (Int,Int)
 flipShape (m,n) = (n,m)
 
--- | @gemv alpha a x beta y@ replaces @y := alpha a * x + beta y@
-gemv :: (BLAS2 e) => e -> DMatrix s (m,n) e -> DVector t n e -> e -> IOVector m e -> IO ()
+-- | @gemv alpha a x beta y@ replaces @y := alpha a * x + beta y@.
+gemv :: (BLAS3 e) => e -> DMatrix s (m,n) e -> DVector t n e -> e -> IOVector m e -> IO ()
 gemv alpha a x beta y
     | numRows a == 0 || numCols a == 0 =
         return ()
-    | V.isConj y = do
-        V.doConj y
-        gemv alpha a x beta (V.conj y)
-        V.doConj y
-    | V.isConj x = do
-        x' <- V.newCopy (V.unsafeThaw x)
-        V.doConj x'
-        gemv alpha a (V.conj x') beta y
+    | isConj x =
+        let b = fromJust $ maybeFromCol x
+        in case maybeFromCol y of
+               Just c  -> gemm alpha a b beta c
+               Nothing -> do
+                   V.doConj y
+                   gemm alpha a b beta (fromJust $ maybeFromCol $ V.conj y)
+                   V.doConj y
+    | isConj y = 
+        let c = fromJust $ maybeFromCol y
+        in case maybeFromCol x of
+               Just b  -> gemm alpha a b beta c
+               Nothing -> do
+                   x' <- V.newCopy (V.unsafeThaw x)
+                   V.doConj x'
+                   gemm alpha a (fromJust $ maybeFromCol $ V.conj x') beta c
     | otherwise =
         let order  = colMajor
             transA = blasTransOf a
@@ -269,30 +240,11 @@
 binaryOp :: (BLAS1 e) => String -> (IOMatrix (m,n) e -> DMatrix t (m,n) e -> IO ()) 
     -> DMatrix s (m,n) e -> DMatrix t (m,n) e -> IO (DMatrix r (m,n) e)
 binaryOp name f a b = 
-    checkMatMatOp name (shape a) (shape b) >> do
+    checkMatMatOp name (shape a) (shape b) $ do
         a' <- newCopy a
         f (unsafeThaw a') b
         return (unsafeCoerce a')
-        
-        
--- | Multiply a matrix by a vector.
-apply :: (BLAS2 e) => Matrix (m,n) e -> Vector n e -> Vector m e
-apply = sapply 1
-
--- | Multiply a scaled matrix by a vector.
-sapply :: (BLAS2 e) => e -> Matrix (m,n) e -> Vector n e -> Vector m e
-sapply alpha a x = unsafePerformIO $ getSApply alpha a x
-{-# NOINLINE sapply #-}
-
--- | Multiply a scaled matrix by a matrix.
-sapplyMat :: (BLAS3 e) => e -> Matrix (m,k) e -> Matrix (k,n) e -> Matrix (m,n) e
-sapplyMat alpha a b = unsafePerformIO $ getSApplyMat alpha a b
-{-# NOINLINE sapplyMat #-}
-    
--- | Multiply a matrix by a matrix.
-applyMat :: (BLAS3 e) => Matrix (m,k) e -> Matrix (k,n) e -> Matrix (m,n) e
-applyMat = sapplyMat 1
-
+            
 -- | Create a new matrix by scaling another matrix by the given value.
 scale :: (BLAS1 e) => e -> Matrix (m,n) e -> Matrix (m,n) e
 scale k a = unsafePerformIO $ getScaled k a
@@ -340,8 +292,5 @@
 "scale.minus/add"  forall k l x y. minus (scale k x) (scale l y) = add k x (-l) y
 "scale1.minus/add" forall k x y.   minus (scale k x) y = add k x (-1) y
 "scale2.minus/add" forall k x y.   minus x (scale k y) = add 1 x (-k) y
-
-"scale.apply/sapply"       forall k a x. apply (scale k a) x = sapply k a x
-"scale.applyMat/sapplyMat" forall k a b. applyMat (scale k a) b = sapplyMat k a b
   #-}
   
diff --git a/Data/Matrix/Diag.hs b/Data/Matrix/Diag.hs
new file mode 100644
--- /dev/null
+++ b/Data/Matrix/Diag.hs
@@ -0,0 +1,180 @@
+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module     : Data.Matrix.Diag
+-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+
+module Data.Matrix.Diag (
+    module BLAS.Matrix,
+    module BLAS.Tensor,
+
+    -- * The diagonal matrix types
+    DiagMatrix,
+    Diag,
+    IODiag,
+    
+    -- * Converting to and from @Vector@s
+    fromVector,
+    toVector,
+
+    ) where
+
+import BLAS.Access
+import BLAS.Elem ( BLAS1, BLAS2 )
+import BLAS.Matrix hiding ( Matrix )
+import qualified BLAS.Matrix as Base
+import BLAS.Tensor
+
+import Control.Monad ( zipWithM_ )
+
+import Data.AEq
+import Data.Matrix.Dense.Internal ( rows )
+import Data.Matrix.Dense.Operations ( unsafeCopyMatrix )
+import qualified Data.Matrix.Dense.Operations as M
+import qualified Data.Matrix.Dense.Internal as M
+import Data.Vector.Dense.Internal
+import qualified Data.Vector.Dense.Internal as V
+import qualified Data.Vector.Dense.Operations as V
+import Data.Vector.Dense.Operations ( (//=), scaleBy, invScaleBy,
+    unsafeCopyVector )
+
+import System.IO.Unsafe
+import Unsafe.Coerce
+
+newtype DiagMatrix t mn e = Diag (DVector t mn e)
+type Diag   = DiagMatrix Imm
+type IODiag = DiagMatrix Mut
+
+coerceDiag :: DiagMatrix t mn e -> DiagMatrix t mn' e
+coerceDiag = unsafeCoerce
+
+fromVector :: DVector t n e -> DiagMatrix t (n,n) e
+fromVector = Diag . unsafeCoerce
+
+toVector :: DiagMatrix t (n,n) e -> DVector t n e
+toVector (Diag x) = unsafeCoerce x
+
+instance (BLAS1 e) => Scalable (DiagMatrix Imm (n,n)) e where
+    (*>) k (Diag x) = Diag $ k *> x
+
+
+replaceHelp :: (BLAS1 e) => 
+       (Vector n e -> [(Int,e)] -> Vector n e) 
+    -> Diag (n,n) e -> [((Int,Int),e)] -> Diag (n,n) e
+replaceHelp f a ijes =
+    let iies = filter (\((i,j),_) -> i == j) ijes
+        ies  = map (\((i,_),e) -> (i,e)) iies
+        x'   = f (toVector a) ies
+    in fromVector x'
+    
+    
+instance (BLAS1 e) => ITensor (DiagMatrix Imm (n,n)) (Int,Int) e where
+    size = size . toVector
+    
+    assocs a =
+        let ies = assocs $ toVector a
+        in map (\(i,e) -> ((i,i),e)) ies
+    
+    (//) = replaceHelp (//)
+    
+    amap f a = fromVector (amap f $ toVector a)
+    
+    unsafeAt a (i,j) | i /= j = 0
+                     | otherwise = unsafeAt (toVector a) i
+                     
+    unsafeReplace = replaceHelp unsafeReplace
+    
+    
+instance (BLAS1 e) => RTensor (DiagMatrix t (n,n)) (Int,Int) e IO where
+    getSize = getSize . toVector
+    
+    newCopy a = do
+        x' <- newCopy $ toVector a
+        return $ fromVector x'
+    
+    unsafeReadElem a (i,j)
+        | i /= j    = return 0
+        | otherwise = unsafeReadElem (toVector a) i
+        
+        
+instance (BLAS1 e) => MTensor (DiagMatrix Mut (n,n)) (Int,Int) e IO where
+    setZero = setZero . toVector
+    
+    setConstant k = setConstant k . toVector
+    
+    canModifyElem a (i,j) = return (i == j && i >= 0 && i < numRows a)
+    
+    unsafeWriteElem a (i,_) = unsafeWriteElem (toVector a) i
+    
+    modifyWith f a = modifyWith f (toVector a)
+    
+        
+instance Base.Matrix (DiagMatrix t) where
+    numRows (Diag x) = dim x
+    numCols = numRows
+
+    herm (Diag x) = unsafeCoerce $ Diag (conj x)
+
+
+instance (BLAS2 e) => IMatrix (DiagMatrix Imm) e
+
+instance (BLAS2 e) => RMatrix (DiagMatrix t) e where
+    unsafeDoSApplyAdd alpha a x beta y = do
+        x' <- newCopy x
+        unsafeDoSApply_ 1 (coerceDiag a) (V.unsafeThaw x')
+        scaleBy beta y
+        V.unsafeAxpy alpha x' (V.coerceVector y)
+
+    unsafeDoSApplyAddMat alpha a b beta c = do
+        M.scaleBy beta c
+        ks <- unsafeInterleaveIO $ getElems (toVector $ coerceDiag a)
+        let (kxs) = zip ks (rows b)
+            ys    = rows c
+        zipWithM_ (\(k,x) y -> V.unsafeAxpy (alpha*k) x y) kxs ys
+
+    unsafeDoSApply_ alpha a x = do
+        V.unsafeTimesEquals x (toVector a)
+        V.scaleBy alpha x
+
+    unsafeDoSApplyMat_ alpha a b = do
+        ks <- unsafeInterleaveIO $ getElems (toVector a)
+        zipWithM_ (\k r -> scaleBy (alpha*k) r) ks (rows b)
+
+
+instance (BLAS2 e) => ISolve (DiagMatrix Imm) e
+
+
+instance (BLAS2 e) => RSolve (DiagMatrix Imm) e where
+    unsafeDoSSolve alpha a y x = do
+        unsafeCopyVector x (unsafeCoerce y)
+        unsafeDoSSolve_ alpha (coerceDiag a) x
+        
+    unsafeDoSSolveMat alpha a c b = do
+        unsafeCopyMatrix b (unsafeCoerce c)
+        unsafeDoSSolveMat_ alpha (coerceDiag a) b
+    
+    unsafeDoSSolve_ alpha a x = do
+        V.scaleBy alpha x
+        x //= (toVector a)
+
+    unsafeDoSSolveMat_ alpha a b = do
+        M.scaleBy alpha b
+        ks <- unsafeInterleaveIO $ getElems (toVector a)
+        zipWithM_ (\k r -> invScaleBy k r) ks (rows b)
+
+
+instance (Show e, BLAS1 e) => Show (DiagMatrix Imm (n,n) e) where
+    show (Diag x) = "fromVector (" ++ show x ++ ")"
+
+
+instance (Eq e, BLAS1 e) => Eq (DiagMatrix Imm (n,n) e) where
+    (==) (Diag x) (Diag y) = (==) x y
+
+
+instance (AEq e, BLAS1 e) => AEq (DiagMatrix Imm (n,n) e) where
+    (===) (Diag x) (Diag y) = (===) x y
+    (~==) (Diag x) (Diag y) = (~==) x y
diff --git a/Data/Matrix/Herm.hs b/Data/Matrix/Herm.hs
--- a/Data/Matrix/Herm.hs
+++ b/Data/Matrix/Herm.hs
@@ -19,46 +19,43 @@
     hermL,
     hermU,
 
+    coerceHerm,
+
     ) where
 
 import Unsafe.Coerce
 
-import qualified BLAS.Elem as E
 import BLAS.Matrix
-import BLAS.Tensor
 import BLAS.Types ( UpLo(..) )
 
-data Herm a nn e = Herm UpLo e (a nn e)
+data Herm a nn e = Herm UpLo (a nn e)
 
+coerceHerm :: Herm a mn e -> Herm a mn' e
+coerceHerm = unsafeCoerce
+
 mapHerm :: (a (n,n) e -> b (n,n) e) -> Herm a (n,n) e -> Herm b (n,n) e
-mapHerm f (Herm u e a) = Herm u e $ f a
+mapHerm f (Herm u a) = Herm u $ f a
 
-fromBase :: UpLo -> e -> a (n,n) e -> Herm a (n,n) e
+fromBase :: UpLo -> a (n,n) e -> Herm a (n,n) e
 fromBase = Herm
         
-toBase :: Herm a (n,n) e -> (UpLo, e, a (n,n) e)
-toBase (Herm u e a) = (u,e,a)
+toBase :: Herm a (n,n) e -> (UpLo, a (n,n) e)
+toBase (Herm u a) = (u,a)
 
-hermL :: (Num e) => a (n,n) e -> Herm a (n,n) e
-hermL = Herm Lower 1
+hermL :: a (n,n) e -> Herm a (n,n) e
+hermL = Herm Lower
 
-hermU :: (Num e) => a (n,n) e -> Herm a (n,n) e
-hermU = Herm Upper 1
+hermU :: a (n,n) e -> Herm a (n,n) e
+hermU = Herm Upper
       
 instance Matrix a => Matrix (Herm a) where
-    numRows (Herm _ _ a) = numRows a
-    numCols (Herm _ _ a) = numCols a
-    herm    (Herm u e a) = Herm u (E.conj e) (unsafeCoerce a)
+    numRows (Herm _ a) = numRows a
+    numCols (Herm _ a) = numCols a
+    herm = coerceHerm
     
-instance (Num e) => Scalable (Herm a nn) e where
-    (*>) k (Herm u e a) = Herm u (k*e) a
-
-instance (Show (a mn e), Show e, Num e) => Show (Herm a mn e) where
-    show (Herm u k a) 
-        | k /= 1 = "(" ++ show k ++ ") *> " ++ show (Herm u 1 a)
-        | otherwise =
-            constructor ++ " (" ++ show a ++ ")"
-        where
-          constructor = case u of
-              Lower -> "hermL"
-              Upper -> "hermU"
+instance Show (a mn e) => Show (Herm a mn e) where
+    show (Herm u a) = constructor ++ " (" ++ show a ++ ")"
+      where
+        constructor = case u of
+            Lower -> "hermL"
+            Upper -> "hermU"
diff --git a/Data/Matrix/Herm/Banded.hs b/Data/Matrix/Herm/Banded.hs
new file mode 100644
--- /dev/null
+++ b/Data/Matrix/Herm/Banded.hs
@@ -0,0 +1,84 @@
+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module     : Data.Matrix.Herm.Banded
+-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+
+module Data.Matrix.Herm.Banded (
+    module Data.Matrix.Herm,
+    module BLAS.Matrix.Immutable,
+    module BLAS.Matrix.ReadOnly,
+    ) where
+
+import Control.Monad ( zipWithM_ )
+
+import BLAS.Access
+import BLAS.C ( BLAS2, colMajor, cblasUpLo )
+import BLAS.Types ( flipUpLo )
+import qualified BLAS.C as BLAS
+
+import Data.Matrix.Banded.Internal
+import Data.Matrix.Dense.Internal ( DMatrix, IOMatrix, coerceMatrix )
+import Data.Vector.Dense.Internal
+import qualified Data.Matrix.Banded.Internal as B
+import qualified Data.Matrix.Dense.Internal as M
+import qualified Data.Vector.Dense.Internal as V
+import qualified Data.Vector.Dense.Operations as V
+
+import Data.Matrix.Herm
+import BLAS.Matrix.Immutable
+import BLAS.Matrix.ReadOnly
+
+instance (BLAS2 e) => IMatrix (Herm (BMatrix Imm)) e where
+
+instance (BLAS2 e) => RMatrix (Herm (BMatrix s)) e where
+    unsafeDoSApplyAdd alpha a x beta y = 
+        hbmv alpha (coerceHerm a) x beta (coerceVector y)
+        
+    unsafeDoSApplyAddMat alpha a b beta c = 
+        hbmm alpha (coerceHerm a) b beta (coerceMatrix c)
+
+
+hbmv :: (BLAS2 e) => e -> Herm (BMatrix t) (n,n) e -> DVector s n e -> e -> IOVector n e -> IO ()
+hbmv alpha h x beta y
+    | numRows h == 0 =
+        return ()
+    | isConj y = do
+        V.doConj y
+        hbmv alpha h x beta (V.conj y)
+        V.doConj y
+    | isConj x = do
+        x' <- newCopy x
+        V.doConj (V.unsafeThaw x')
+        hbmv alpha h (conj x') beta y
+    | otherwise =
+        let order = colMajor
+            (u,a) = toBase h
+            n     = numCols a
+            k     = case u of 
+                        Upper -> numUpper a
+                        Lower -> numLower a      
+            u'    = case (isHerm a) of
+                        True  -> flipUpLo u
+                        False -> u
+            uploA = cblasUpLo u'
+            ldA   = ldaOf a
+            incX  = strideOf x
+            incY  = strideOf y
+            withPtrA 
+                  = case u' of Upper -> B.unsafeWithBasePtr a
+                               Lower -> B.unsafeWithElemPtr a (0,0)
+                    
+        in withPtrA $ \pA ->
+               V.unsafeWithElemPtr x 0 $ \pX ->
+                    V.unsafeWithElemPtr y 0 $ \pY -> do
+                        BLAS.hbmv order uploA n k alpha pA ldA pX incX beta pY incY
+
+
+hbmm :: (BLAS2 e) => e -> Herm (BMatrix t) (m,m) e -> DMatrix s (m,n) e -> e -> IOMatrix (m,n) e -> IO ()
+hbmm alpha h b beta c =
+    zipWithM_ (\x y -> hbmv alpha h x beta y) (M.cols b) (M.cols c)
diff --git a/Data/Matrix/Herm/Dense.hs b/Data/Matrix/Herm/Dense.hs
--- a/Data/Matrix/Herm/Dense.hs
+++ b/Data/Matrix/Herm/Dense.hs
@@ -12,20 +12,13 @@
     module Data.Matrix.Herm,
     module BLAS.Matrix.Immutable,
     module BLAS.Matrix.ReadOnly,
-
-    hemv,
-    hemm,
     ) where
 
 import Control.Monad ( zipWithM_ )
-import System.IO.Unsafe
-import Unsafe.Coerce
 
 import BLAS.Access
-import BLAS.Elem ( BLAS2, BLAS3 )
-import BLAS.C ( colMajor, rightSide, leftSide, cblasUpLo )
+import BLAS.C ( BLAS2, BLAS3, colMajor, rightSide, leftSide, cblasUpLo )
 import BLAS.Types ( flipUpLo )
-import qualified BLAS.Elem as E
 import qualified BLAS.C as BLAS
 
 import Data.Matrix.Dense.Internal
@@ -40,23 +33,13 @@
 
 
 instance (BLAS3 e) => IMatrix (Herm (DMatrix Imm)) e where
-    (<*>) h x = unsafePerformIO $ getApply h x
-    {-# NOINLINE (<*>) #-}
 
-    (<**>) h a = unsafePerformIO $ getApplyMat h a
-    {-# NOINLINE (<**>) #-}
-
-
 instance (BLAS3 e) => RMatrix (Herm (DMatrix s)) e where
-    getApply h x = do
-        y <- newZero (dim x)
-        hemv 1 (unsafeCoerce h) x 1 y
-        return (unsafeCoerce y)
-    
-    getApplyMat h a = do
-        b <- newZero (shape a)
-        hemm 1 (unsafeCoerce h) a 1 b
-        return (unsafeCoerce b)
+    unsafeDoSApplyAdd alpha a x beta y = 
+        hemv alpha (coerceHerm a) x beta (coerceVector y)
+        
+    unsafeDoSApplyAddMat alpha a b beta c = 
+        hemm alpha (coerceHerm a) b beta (coerceMatrix c)
 
 
 hemv :: (BLAS2 e) => e -> Herm (DMatrix t) (n,n) e -> DVector s n e -> e -> IOVector n e -> IO ()
@@ -72,22 +55,20 @@
         V.doConj (V.unsafeThaw x')
         hemv alpha h (conj x') beta y
     | otherwise =
-        let order   = colMajor
-            (u,e,a) = toBase h
-            alpha'  = alpha * e
-            n       = numCols a
-            (u',alpha'') 
-                    = case (isHerm a) of
-                          True  -> (flipUpLo u, E.conj alpha')
-                          False -> (u, alpha')
-            uploA   = cblasUpLo u'
-            ldA     = ldaOf a
-            incX    = strideOf x
-            incY    = strideOf y
+        let order = colMajor
+            (u,a) = toBase h
+            n     = numCols a
+            u'    = case isHerm a of
+                        True  -> flipUpLo u
+                        False -> u
+            uploA = cblasUpLo u'
+            ldA   = ldaOf a
+            incX  = strideOf x
+            incY  = strideOf y
         in M.unsafeWithElemPtr a (0,0) $ \pA ->
                V.unsafeWithElemPtr x 0 $ \pX ->
                     V.unsafeWithElemPtr y 0 $ \pY -> do
-                        BLAS.hemv order uploA n alpha'' pA ldA pX incX beta pY incY
+                        BLAS.hemv order uploA n alpha pA ldA pX incX beta pY incY
 
 hemm :: (BLAS3 e) => e -> Herm (DMatrix t) (m,m) e -> DMatrix s (m,n) e -> e -> IOMatrix (m,n) e -> IO ()
 hemm alpha h b beta c
@@ -97,11 +78,10 @@
     | otherwise =
         let order   = colMajor
             (m,n)   = shape c
-            alpha'  = alpha * e
-            (side,u',m',n', alpha'')
+            (side,u',m',n')
                     = if isHerm a
-                          then (rightSide, flipUpLo u, n, m, E.conj alpha')
-                          else (leftSide,  u, m, n, alpha')
+                          then (rightSide, flipUpLo u, n, m)
+                          else (leftSide,  u,          m, n)
             uploA   = cblasUpLo u'
             ldA     = ldaOf a
             ldB     = ldaOf b
@@ -109,6 +89,6 @@
         in M.unsafeWithElemPtr a (0,0) $ \pA ->
                M.unsafeWithElemPtr b (0,0) $ \pB ->
                    M.unsafeWithElemPtr c (0,0) $ \pC ->
-                       BLAS.hemm order side uploA m' n' alpha'' pA ldA pB ldB beta pC ldC
+                       BLAS.hemm order side uploA m' n' alpha pA ldA pB ldB beta pC ldC
     where
-      (u,e,a) = toBase h
+      (u,a) = toBase h
diff --git a/Data/Matrix/Perm.hs b/Data/Matrix/Perm.hs
new file mode 100644
--- /dev/null
+++ b/Data/Matrix/Perm.hs
@@ -0,0 +1,187 @@
+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module     : Data.Matrix.Perm
+-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+
+module Data.Matrix.Perm (
+    module BLAS.Matrix.Base,
+    module BLAS.Matrix.Immutable,
+    module BLAS.Matrix.ReadOnly,
+    module BLAS.Matrix.Solve.Immutable,
+    module BLAS.Matrix.Solve.ReadOnly,
+
+    Perm(..),
+    
+    -- * The identity permutation
+    identity,
+    
+    -- * Converting to/from @Permutation@s
+    fromPermutation,
+    toPermutation,
+    
+    -- * Coercing
+    coercePerm
+
+    ) where
+
+import Control.Monad ( forM_ )
+import Foreign ( peek, poke )
+
+import BLAS.Elem ( Elem, BLAS1 )
+import qualified BLAS.Elem as E
+
+import BLAS.Matrix.Base hiding ( Matrix )
+import qualified BLAS.Matrix.Base as Base
+import BLAS.Matrix.Immutable
+import BLAS.Matrix.ReadOnly
+import BLAS.Matrix.Solve.Immutable
+import BLAS.Matrix.Solve.ReadOnly
+
+import Data.AEq
+
+import Data.Vector.Dense.IO ( dim, isConj, conj, unsafeSwapVectors, 
+    unsafeCopyVector, unsafeWithElemPtr, unsafeReadElem, unsafeWriteElem, 
+    coerceVector )
+import qualified Data.Vector.Dense.IO as V
+    
+import Data.Matrix.Dense.IO ( unsafeRow, unsafeCopyMatrix, coerceMatrix )
+import qualified Data.Matrix.Dense.IO as M
+
+import Data.Permutation ( Permutation )
+import qualified Data.Permutation as P
+
+import Unsafe.Coerce
+
+
+data Perm mn e = 
+      P { baseOf :: !Permutation
+        , isHerm :: !Bool
+        }
+    | I !Int
+
+identity :: Int -> Perm (n,n) e
+identity = I
+
+fromPermutation :: Permutation -> Perm (n,n) e
+fromPermutation = flip P False
+
+toPermutation :: Perm (n,n) e -> Permutation
+toPermutation (I n)       = P.identity n
+toPermutation (P sigma h) = if h then P.inverse sigma else sigma
+
+coercePerm :: Perm mn e -> Perm mn' e
+coercePerm = unsafeCoerce
+
+          
+instance Base.Matrix Perm where
+    numRows (P sigma _) = P.size sigma
+    numRows (I n)       = n
+    
+    numCols a = numRows a
+    
+    herm a@(P _ _) = a{ isHerm=(not . isHerm) a }
+    herm a@(I _)   = (unsafeCoerce a)
+
+
+instance (BLAS1 e) => IMatrix Perm e where
+          
+         
+instance (BLAS1 e) => RMatrix Perm e where
+    unsafeDoApply_   (I _)   _ = return ()
+    unsafeDoApply_ p@(P _ _) x
+        | isHerm p  = P.invertWith swap sigma
+        | otherwise = P.applyWith swap sigma
+      where
+        sigma = baseOf p
+        swap i j = 
+            unsafeWithElemPtr x i $ \pI ->
+                unsafeWithElemPtr x j $ \pJ -> do
+                    xi <- peek pI
+                    xj <- peek pJ
+                    poke pI xj
+                    poke pJ xi
+
+
+    unsafeDoApplyMat_   (I _)   _ = return ()
+    unsafeDoApplyMat_ p@(P _ _) a
+        | isHerm p  = P.invertWith swap sigma
+        | otherwise = P.applyWith swap sigma
+      where
+        sigma = baseOf p
+        swap i j = unsafeSwapVectors (unsafeRow a i) (unsafeRow a j)
+
+
+    unsafeDoSApply k (I _) x y = do
+        unsafeCopyVector y (coerceVector x)
+        V.scaleBy k y
+    unsafeDoSApply k p x y
+        | isConj x =
+            unsafeDoSApply (E.conj k) p (conj x) (conj y)
+        | otherwise =
+            let n     = dim x
+                sigma = baseOf p
+            in do
+                forM_ [0..(n-1)] $ \i ->
+                    let i' = P.unsafeApply sigma i
+                    in case (isHerm p) of
+                           False -> unsafeReadElem x i  
+                                    >>= unsafeWriteElem y i'
+                           True  -> unsafeReadElem x i' 
+                                    >>= unsafeWriteElem y i
+                V.scaleBy k y
+
+
+    unsafeDoSApplyMat alpha (I _) b c = do
+        unsafeCopyMatrix c (coerceMatrix b)
+        M.scaleBy alpha c
+    unsafeDoSApplyMat alpha p b c =
+        let m     = numCols p
+            sigma = baseOf p
+        in do
+            forM_ [0..(m-1)] $ \i ->
+                let i' = P.unsafeApply sigma i
+                in case (isHerm p) of
+                       False -> unsafeCopyVector (unsafeRow c i') 
+                                                 (unsafeRow b i)
+                       True  -> unsafeCopyVector (unsafeRow c i)
+                                                 (unsafeRow b i')
+            M.scaleBy alpha c
+
+
+instance (BLAS1 e) => ISolve Perm e where
+    
+instance (BLAS1 e) => RSolve Perm e where    
+    unsafeDoSolve_ p          = unsafeDoApply_ (herm p)
+    unsafeDoSolveMat_ p       = unsafeDoApplyMat_ (herm p)
+    unsafeDoSSolve alpha p    = unsafeDoSApply alpha (coercePerm $ herm p)
+    unsafeDoSSolveMat alpha p = unsafeDoSApplyMat alpha (coercePerm $ herm p)
+
+
+instance (Elem e) => Show (Perm (n,n) e) where
+    show (I n)         = "identity " ++ show n
+    show p | isHerm p  = "herm (" ++ show (herm p) ++ ")"
+           | otherwise = "fromPermutation (" ++ show (baseOf p) ++ ")"
+    
+    
+instance (Elem e) => Eq (Perm (n,n) e) where
+    (==) (I n) (I n') = n == n'
+    (==) (I n) p
+        | isHerm p   = (==) (I n) (herm p)
+        | otherwise  = (==) (fromPermutation $ P.identity n) p
+    (==) p     (I n) = (==) (I n) p
+
+    (==) (P sigma h) (P sigma' h') 
+        | h == h'   = sigma == sigma'
+        | otherwise = P.size sigma == P.size sigma'
+                      && sigma == (P.inverse sigma')
+
+
+instance (Elem e) => AEq (Perm (n,n) e) where
+    (===) = (==)
+    (~==) = (==)
+    
diff --git a/Data/Matrix/Tri.hs b/Data/Matrix/Tri.hs
--- a/Data/Matrix/Tri.hs
+++ b/Data/Matrix/Tri.hs
@@ -17,57 +17,109 @@
     mapTri,
 
     lower,
+    lowerFat,
+    lowerTall,
+    
     lowerU,
+    lowerUFat,
+    lowerUTall,
+    
     upper,
+    upperFat,
+    upperTall,
+    
     upperU,
+    upperUFat,
+    upperUTall,
 
+    coerceTri,
+
     ) where
 
-import qualified BLAS.Elem as E
+import BLAS.Internal ( checkSquare, checkFat, checkTall )
 import BLAS.Matrix
 import BLAS.Tensor
 import BLAS.Types ( UpLo(..), Diag(..), flipUpLo )
 
-data Tri a nn e = Tri UpLo Diag e (a nn e)
+import Unsafe.Coerce
 
-mapTri :: (a (n,n) e -> b (n,n) e) -> Tri a (n,n) e -> Tri b (n,n) e
-mapTri f (Tri u d n a) = Tri u d n $ f a
+data Tri a mn e = Tri UpLo Diag (a mn e)
 
-fromBase :: UpLo -> Diag -> e -> a (n,n) e -> Tri a (n,n) e
+-- | Coerce the shape type.
+coerceTri :: Tri a mn e -> Tri a mn' e
+coerceTri = unsafeCoerce
+
+mapTri :: (a (m,n) e -> b (m,n) e) -> Tri a (m,n) e -> Tri b (m,n) e
+mapTri f (Tri u d a) = Tri u d $ f a
+
+fromBase :: UpLo -> Diag -> a (m,n) e -> Tri a (m,n) e
 fromBase = Tri
         
-toBase :: Tri a (n,n) e -> (UpLo, Diag, e, a (n,n) e)
-toBase (Tri u d e a) = (u,d,e,a)
+toBase :: Tri a (m,n) e -> (UpLo, Diag, a (m,n) e)
+toBase (Tri u d a) = (u,d,a)
 
-lower :: (Num e) => a (n,n) e -> Tri a (n,n) e
-lower = Tri Lower NonUnit 1
 
-lowerU :: (Num e) => a (n,n) e -> Tri a (n,n) e
-lowerU = Tri Lower Unit 1
+lower :: (Matrix a) => a (n,n) e -> Tri a (n,n) e
+lower a = checkSquare (shape a) $ Tri Lower NonUnit a
 
-upper :: (Num e) => a (n,n) e -> Tri a (n,n) e
-upper = Tri Upper NonUnit 1
+lowerFat :: (Matrix a) => a (m,n) e -> Tri a (m,m) e
+lowerFat a = checkFat (shape a) $ Tri Lower NonUnit (unsafeCoerce a)
 
-upperU :: (Num e) => a (n,n) e -> Tri a (n,n) e
-upperU = Tri Upper Unit 1
+lowerTall :: (Matrix a) => a (m,n) e -> Tri a (m,n) e
+lowerTall a = checkTall (shape a) $ Tri Lower NonUnit a
+
+
+lowerU :: (Matrix a) => a (n,n) e -> Tri a (n,n) e
+lowerU a = checkSquare (shape a) $ Tri Lower Unit a
+
+lowerUFat :: (Matrix a) => a (m,n) e -> Tri a (m,m) e
+lowerUFat a = checkFat (shape a) $ Tri Lower Unit (unsafeCoerce a)
+
+lowerUTall :: (Matrix a) => a (m,n) e -> Tri a (m,n) e
+lowerUTall a = checkTall (shape a) $ Tri Lower Unit a
+
+
+upper :: (Matrix a) => a (n,n) e -> Tri a (n,n) e
+upper a = checkSquare (shape a) $ Tri Upper NonUnit a
+
+upperFat :: (Matrix a) => a (m,n) e -> Tri a (m,n) e
+upperFat a = checkFat (shape a) $ Tri Upper NonUnit a
+
+upperTall :: (Matrix a) => a (m,n) e -> Tri a (n,n) e
+upperTall a = checkTall (shape a) $ Tri Upper NonUnit (unsafeCoerce a)
+
+
+upperU :: (Matrix a) => a (n,n) e -> Tri a (n,n) e
+upperU a = checkSquare (shape a) $ Tri Upper Unit a
+
+upperUFat :: (Matrix a) => a (m,n) e -> Tri a (m,n) e
+upperUFat a = checkFat (shape a) $ Tri Upper Unit a
+
+upperUTall :: (Matrix a) => a (m,n) e -> Tri a (n,n) e
+upperUTall a = checkTall (shape a) $ Tri Upper Unit (unsafeCoerce a)
+
       
 instance Matrix a => Matrix (Tri a) where
-    numRows (Tri _ _ _ a) = numRows a
-    numCols (Tri _ _ _ a) = numCols a
-    herm    (Tri u d e a) = Tri (flipUpLo u) d (E.conj e) (herm a)
+    numRows (Tri Lower _ a) = numRows a
+    numRows (Tri Upper _ a) = min (numRows a) (numCols a)
     
-instance (Num e) => Scalable (Tri a nn) e where
-    (*>) k (Tri u d e a) = Tri u d (k*e) a
+    numCols (Tri Lower _ a) = min (numRows a) (numCols a)
+    numCols (Tri Upper _ a) = numCols a
+    
+    herm (Tri u d a) = Tri (flipUpLo u) d (herm a)
 
-instance (Show (a mn e), Show e, Num e) => Show (Tri a mn e) where
-    show (Tri u d k a) 
-        | k /= 1 = "(" ++ show k ++ ") *> " ++ show (Tri u d 1 a)
-        | otherwise =
-            constructor ++ " (" ++ show a ++ ")"
+
+instance (Show (a (m,n) e), Matrix a) => Show (Tri a (m,n) e) where
+    show (Tri u d a) =
+        constructor ++ suffix ++ " (" ++ show a ++ ")"
         where
           constructor = case (u,d) of
               (Lower, NonUnit) -> "lower"
               (Lower, Unit   ) -> "lowerU"
               (Upper, NonUnit) -> "upper"
               (Upper, Unit   ) -> "upperU"
-        
+
+          suffix = case undefined of
+                       _ | isSquare a -> ""
+                       _ | isFat a    -> "Fat"
+                       _              -> "Tall"
diff --git a/Data/Matrix/Tri/Banded.hs b/Data/Matrix/Tri/Banded.hs
new file mode 100644
--- /dev/null
+++ b/Data/Matrix/Tri/Banded.hs
@@ -0,0 +1,128 @@
+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module     : Data.Matrix.Tri.Banded
+-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+
+module Data.Matrix.Tri.Banded (
+    module Data.Matrix.Tri,
+    module BLAS.Matrix.Immutable,
+    module BLAS.Matrix.ReadOnly,
+    module BLAS.Matrix.Solve,
+    ) where
+
+
+import Data.Matrix.Banded.Internal
+import qualified Data.Matrix.Banded.Internal as B
+import Data.Matrix.Dense.IO ( IOMatrix, coerceMatrix, unsafeCopyMatrix )
+import qualified Data.Matrix.Dense.IO as M
+import Data.Vector.Dense.IO ( IOVector, isConj, strideOf, conj, coerceVector, 
+    unsafeCopyVector )
+import qualified Data.Vector.Dense.IO as V
+
+import BLAS.Access
+import BLAS.Types ( flipUpLo )
+
+import BLAS.C ( BLAS2, cblasDiag, cblasUpLo, colMajor, noTrans, conjTrans )                                   
+import qualified BLAS.C as BLAS
+
+import BLAS.Matrix.Immutable
+import BLAS.Matrix.ReadOnly
+import BLAS.Matrix.Solve
+
+import Data.Matrix.Tri
+
+instance (BLAS2 e) => IMatrix (Tri (BMatrix Imm)) e where
+
+instance (BLAS2 e) => RMatrix (Tri (BMatrix s)) e where
+
+    unsafeDoSApply alpha a x y = do
+        unsafeCopyVector (coerceVector y) x
+        unsafeDoSApply_ alpha (coerceTri a) y
+    
+    unsafeDoSApplyMat alpha a b c = do
+        unsafeCopyMatrix (coerceMatrix c) b
+        unsafeDoSApplyMat_ alpha (coerceTri a) c
+        
+    unsafeDoSApply_    = tbmv
+    unsafeDoSApplyMat_ = tbmm
+
+instance (BLAS2 e) => ISolve (Tri (BMatrix Imm)) e where
+
+instance (BLAS2 e) => RSolve (Tri (BMatrix s)) e where
+    unsafeDoSSolve alpha a y x = do
+        unsafeCopyVector (coerceVector x) y
+        unsafeDoSSolve_ alpha (coerceTri a) x
+    
+    unsafeDoSSolveMat alpha a c b = do
+        unsafeCopyMatrix (coerceMatrix b) c
+        unsafeDoSSolveMat_ alpha (coerceTri a) b
+
+    unsafeDoSSolve_    = tbsv
+    unsafeDoSSolveMat_ = tbsm
+
+
+tbmv :: (BLAS2 e) => e -> Tri (BMatrix t) (n,n) e -> IOVector n e -> IO ()
+tbmv alpha t x | isConj x = do
+    V.doConj x
+    tbmv alpha t (conj x)
+    V.doConj x
+
+tbmv alpha t x =
+    let (u,d,a) = toBase t
+        order     = colMajor
+        (transA,u') = if isHerm a then (conjTrans, flipUpLo u) else (noTrans, u)
+        uploA     = cblasUpLo u'
+        diagA     = cblasDiag d
+        n         = numCols a
+        k         = case u of Upper -> numUpper a 
+                              Lower -> numLower a
+        ldA       = ldaOf a
+        incX      = strideOf x
+        withPtrA  = case u' of 
+                        Upper -> B.unsafeWithBasePtr a
+                        Lower -> B.unsafeWithElemPtr a (0,0)
+    in withPtrA $ \pA ->
+           V.unsafeWithElemPtr x 0 $ \pX -> do
+               V.scaleBy alpha x
+               BLAS.tbmv order uploA transA diagA n k pA ldA pX incX
+              
+               
+tbmm :: (BLAS2 e) => e -> Tri (BMatrix t) (m,m) e -> IOMatrix (m,n) e -> IO ()
+tbmm 1     t b = mapM_ (\x -> tbmv 1 t x) (M.cols b)
+tbmm alpha t b = M.scaleBy alpha b >> tbmm 1 t b
+
+
+tbsv :: (BLAS2 e) => e -> Tri (BMatrix t) (n,n) e -> IOVector n e -> IO ()
+tbsv alpha t x | isConj x = do
+    V.doConj x
+    tbsv alpha t (conj x)
+    V.doConj x
+    
+tbsv alpha t x = 
+    let (u,d,a) = toBase t
+        order     = colMajor
+        (transA,u') = if isHerm a then (conjTrans, flipUpLo u) else (noTrans, u)
+        uploA     = cblasUpLo u'
+        diagA     = cblasDiag d
+        n         = numCols a
+        k         = case u of Upper -> numUpper a 
+                              Lower -> numLower a        
+        ldA       = ldaOf a
+        incX      = strideOf x
+        withPtrA  = case u' of 
+                        Upper -> B.unsafeWithBasePtr a
+                        Lower -> B.unsafeWithElemPtr a (0,0)
+    in withPtrA $ \pA ->
+           V.unsafeWithElemPtr x 0 $ \pX -> do
+               V.scaleBy alpha x
+               BLAS.tbsv order uploA transA diagA n k pA ldA pX incX
+
+
+tbsm :: (BLAS2 e) => e -> Tri (BMatrix t) (m,m) e -> IOMatrix (m,n) e -> IO ()
+tbsm 1     t b = mapM_ (\x -> tbsv 1 t x) (M.cols b)
+tbsm alpha t b = M.scaleBy alpha b >> tbsm 1 t b
diff --git a/Data/Matrix/Tri/Dense.hs b/Data/Matrix/Tri/Dense.hs
--- a/Data/Matrix/Tri/Dense.hs
+++ b/Data/Matrix/Tri/Dense.hs
@@ -17,22 +17,24 @@
     trmv,
     trsv,
     trmm,    
-    trsm
+    trsm,
+    
     ) where
 
 import Control.Monad ( when )
 import Data.Maybe ( fromJust )
-import System.IO.Unsafe
-import Unsafe.Coerce
 
-import Data.Matrix.Dense.Internal
+import Data.Matrix.Dense.Internal hiding ( diag )
 import qualified Data.Matrix.Dense.Internal as M
+import qualified Data.Matrix.Dense.Operations as M
 import Data.Vector.Dense.Internal
+import Data.Vector.Dense.Operations( unsafeCopyVector )
+import Data.Matrix.Dense.Operations( unsafeCopyMatrix )
 import qualified Data.Vector.Dense.Internal as V
 import qualified Data.Vector.Dense.Operations as V
 
 import BLAS.Access
-import BLAS.Elem ( BLAS3 )
+import BLAS.Elem ( Elem, BLAS3 )
 import qualified BLAS.Elem as E
 import BLAS.C.Types ( cblasDiag, cblasUpLo, cblasTrans, colMajor, 
     noTrans, conjTrans, leftSide, rightSide )
@@ -48,69 +50,179 @@
 import Data.Matrix.Tri
 
 instance (BLAS3 e) => IMatrix (Tri (DMatrix Imm)) e where
-    (<*>) t x = unsafePerformIO $ getApply t x
-    {-# NOINLINE (<*>) #-}
+instance (BLAS3 e) => ISolve (Tri (DMatrix Imm)) e where
 
-    (<**>) t a = unsafePerformIO $ getApplyMat t a
-    {-# NOINLINE (<**>) #-}
 
-instance (BLAS3 e) => ISolve (Tri (DMatrix Imm)) e where
-    (<\>) t x = unsafePerformIO $ getSolve t x
-    {-# NOINLINE (<\>) #-}
+toLower :: (Elem e) => Diag -> DMatrix s (m,n) e 
+        -> Either (Tri (DMatrix s) (m,m) e) 
+                  (Tri (DMatrix s) (n,n) e, DMatrix s (d,n) e)
+toLower diag a =
+    if m <= n
+        then Left $  fromBase Lower diag (unsafeSubmatrix a (0,0) (m,m))
+        else let t = fromBase Lower diag (unsafeSubmatrix a (0,0) (n,n))
+                 r = unsafeSubmatrix a (n,0) (d,n)
+             in Right (t,r)
+  where
+    (m,n) = shape a
+    d     = m - n
+    
+toUpper :: (Elem e) => Diag -> DMatrix s (m,n) e
+        -> Either (Tri (DMatrix s) (n,n) e)
+                  (Tri (DMatrix s) (m,m) e, DMatrix s (m,d) e)
+toUpper diag a =
+    if n <= m
+        then Left $  fromBase Upper diag (unsafeSubmatrix a (0,0) (n,n))
+        else let t = fromBase Upper diag (unsafeSubmatrix a (0,0) (m,m))
+                 r = unsafeSubmatrix a (0,m) (m,d)
+             in Right (t,r)
+  where
+    (m,n) = shape a
+    d     = n - m
 
-    (<\\>) t a = unsafePerformIO $ getSolveMat t a
-    {-# NOINLINE (<\\>) #-}
 
 instance (BLAS3 e) => RMatrix (Tri (DMatrix s)) e where
-    getApply t x = do
-        x' <- newCopy x
-        trmv (unsafeCoerce t) (V.unsafeThaw x')
-        return (unsafeCoerce x')
+    unsafeDoSApply_    = trmv
+    unsafeDoSApplyMat_ = trmm
     
-    getApplyMat t a = do
-        a' <- newCopy a
-        trmm (unsafeCoerce t) (M.unsafeThaw a')
-        return (unsafeCoerce a')
+    unsafeDoSApply alpha t x y =
+        case (u, toLower d a, toUpper d a) of
+            (Lower,Left t',_) -> do
+                unsafeCopyVector y (coerceVector x)
+                trmv alpha t' y
+                
+            (Lower,Right (t',r),_) -> do
+                let y1 = unsafeSubvector y 0            (numRows t')
+                    y2 = unsafeSubvector y (numRows t') (numRows r)
+                unsafeCopyVector y1 x
+                trmv alpha t' y1
+                unsafeDoSApply alpha r x y2
+                
+            (Upper,_,Left t') -> do
+                unsafeCopyVector (coerceVector y) x
+                trmv alpha t' (coerceVector y)
 
-instance (BLAS3 e) => RSolve (Tri (DMatrix s)) e where
-    getSolve t x = do
-        x' <- newCopy x
-        trsv (unsafeCoerce t) (V.unsafeThaw x')
-        return (unsafeCoerce x')
+            (Upper,_,Right (t',r)) ->
+                let x1 = unsafeSubvector x 0            (numCols t')
+                    x2 = unsafeSubvector x (numCols t') (numCols r)
+                in do
+                    unsafeCopyVector y x1
+                    trmv alpha t' y
+                    unsafeDoSApplyAdd alpha r x2 1 y
+      where
+        (u,d,a) = toBase t
+
+
+    unsafeDoSApplyMat alpha t b c =
+        case (u, toLower d a, toUpper d a) of
+            (Lower,Left t',_) -> do
+                unsafeCopyMatrix c (coerceMatrix b)
+                trmm alpha t' c
+                
+            (Lower,Right (t',r),_) -> do
+                let c1 = unsafeSubmatrix c (0,0)          (numRows t',numCols c)
+                    c2 = unsafeSubmatrix c (numRows t',0) (numRows r ,numCols c)
+                unsafeCopyMatrix c1 b
+                trmm alpha t' c1
+                unsafeDoSApplyMat alpha r b c2
+                
+            (Upper,_,Left t') -> do
+                unsafeCopyMatrix (coerceMatrix c) b
+                trmm alpha t' (coerceMatrix c)
+
+            (Upper,_,Right (t',r)) ->
+                let b1 = unsafeSubmatrix b (0,0)          (numCols t',numCols b)
+                    b2 = unsafeSubmatrix b (numCols t',0) (numCols r ,numCols b)
+                in do
+                    unsafeCopyMatrix c b1
+                    trmm alpha t' c
+                    unsafeDoSApplyAddMat alpha r b2 1 c
+      where
+        (u,d,a) = toBase t
+        
     
-    getSolveMat t a = do
-        a' <- newCopy a
-        trsm (unsafeCoerce t) (M.unsafeThaw a')
-        return (unsafeCoerce a')
+instance (BLAS3 e) => RSolve (Tri (DMatrix s)) e where
+    unsafeDoSSolve_    = trsv
+    unsafeDoSSolveMat_ = trsm
 
+    unsafeDoSSolve alpha t y x =
+        case (u, toLower d a, toUpper d a) of
+            (Lower,Left t',_) -> do
+                unsafeCopyVector x (coerceVector y)
+                trsv alpha t' (coerceVector x)
+                
+            (Lower,Right (t',_),_) -> do
+                let y1 = unsafeSubvector y 0            (numRows t')
+                unsafeCopyVector x y1
+                trsv alpha t' x
+                
+            (Upper,_,Left t') -> do
+                unsafeCopyVector x (coerceVector y)
+                trsv alpha t' x
 
-trmv :: (BLAS3 e) => Tri (DMatrix t) (n,n) e -> IOVector n e -> IO ()
-trmv _ x
-    | dim x == 0 = return ()
-trmv t x
+            (Upper,_,Right (t',r)) ->
+                let x1 = unsafeSubvector x 0            (numCols t')
+                    x2 = unsafeSubvector x (numCols t') (numCols r)
+                in do
+                    unsafeCopyVector x1 y
+                    trsv alpha t' x1
+                    setZero x2
+      where
+        (u,d,a) = toBase t
+
+
+    unsafeDoSSolveMat alpha t c b =
+        case (u, toLower d a, toUpper d a) of
+            (Lower,Left t',_) -> do
+                unsafeCopyMatrix b (coerceMatrix c)
+                trsm alpha t' (coerceMatrix b)
+                
+            (Lower,Right (t',_),_) -> do
+                let c1 = unsafeSubmatrix c (0,0)          (numRows t',numCols c)
+                unsafeCopyMatrix b c1
+                trsm alpha t' b
+                
+            (Upper,_,Left t') -> do
+                unsafeCopyMatrix (coerceMatrix b) c
+                trsm alpha t' (coerceMatrix b)
+
+            (Upper,_,Right (t',r)) ->
+                let b1 = unsafeSubmatrix b (0,0)          (numCols t',numCols b)
+                    b2 = unsafeSubmatrix b (numCols t',0) (numCols r ,numCols b)
+                in do
+                    unsafeCopyMatrix b1 c
+                    trsm alpha t' b1
+                    setZero b2
+      where
+        (u,d,a) = toBase t
+
+
+trmv :: (BLAS3 e) => e -> Tri (DMatrix t) (n,n) e -> IOVector n e -> IO ()
+trmv alpha t x 
+    | dim x == 0 = 
+        return ()
     | isConj x =
         let b = fromJust $ maybeFromCol x
-        in trmm t b
-trmv t x =
-    let (u,d,alpha,a) = toBase t
-        order     = colMajor
-        (transA,u') = if isHerm a then (conjTrans, flipUpLo u) else (noTrans, u)
-        uploA     = cblasUpLo u'
-        diagA     = cblasDiag d
-        n         = numCols a
-        ldA       = ldaOf a
-        incX      = strideOf x
-    in M.unsafeWithElemPtr a (0,0) $ \pA ->
-           V.unsafeWithElemPtr x 0 $ \pX -> do
-               BLAS.trmv order uploA transA diagA n pA ldA pX incX
-               when (alpha /= 1) $ V.scaleBy alpha x
-               
+        in trmm alpha t b
+    | otherwise =
+        let (u,d,a)   = toBase t
+            order     = colMajor
+            (transA,u') = if isHerm a then (conjTrans, flipUpLo u) else (noTrans, u)
+            uploA     = cblasUpLo u'
+            diagA     = cblasDiag d
+            n         = dim x
+            ldA       = ldaOf a
+            incX      = strideOf x
+        in M.unsafeWithElemPtr a (0,0) $ \pA ->
+               V.unsafeWithElemPtr x 0 $ \pX -> do
+                   when (alpha /= 1) $ V.scaleBy alpha x
+                   BLAS.trmv order uploA transA diagA n pA ldA pX incX
+
                
-trmm :: (BLAS3 e) => Tri (DMatrix t) (m,m) e -> IOMatrix (m,n) e -> IO ()
-trmm _ b
+trmm :: (BLAS3 e) => e -> Tri (DMatrix t) (m,m) e -> IOMatrix (m,n) e -> IO ()
+trmm _ _ b
     | M.numRows b == 0 || M.numCols b == 0 = return ()
-trmm t b =
-    let (u,d,alpha,a) = toBase t
+trmm alpha t b =
+    let (u,d,a)   = toBase t
         order     = colMajor
         (h,u')    = if isHerm a then (ConjTrans, flipUpLo u) else (NoTrans, u)
         (m,n)     = shape b
@@ -127,32 +239,34 @@
            M.unsafeWithElemPtr b (0,0) $ \pB ->
                BLAS.trmm order side uploA transA diagA m' n' alpha' pA ldA pB ldB
 
-trsv :: (BLAS3 e) =>Tri (DMatrix t) (n,n) e -> IOVector n e -> IO ()
-trsv _ x
+
+trsv :: (BLAS3 e) => e -> Tri (DMatrix t) (n,n) e -> IOVector n e -> IO ()
+trsv _ _ x
     | dim x == 0 = return ()
-trsv t x
+trsv alpha t x
     | isConj x =
         let b = fromJust $ maybeFromCol x
-        in trsm t b
-trsv t x =
-    let (u,d,alpha,a) = toBase t
+        in trsm alpha t b
+trsv alpha t x =
+    let (u,d,a) = toBase t
         order     = colMajor
         (transA,u') = if isHerm a then (conjTrans, flipUpLo u) else (noTrans, u)
         uploA     = cblasUpLo u'
         diagA     = cblasDiag d
-        n         = numCols a
+        n         = dim x
         ldA       = ldaOf a
         incX      = strideOf x
     in M.unsafeWithElemPtr a (0,0) $ \pA ->
            V.unsafeWithElemPtr x 0 $ \pX -> do
+               when (alpha /= 1) $ V.scaleBy alpha x
                BLAS.trsv order uploA transA diagA n pA ldA pX incX
-               when (alpha /= 1) $ V.invScaleBy alpha x
 
-trsm :: (BLAS3 e) => Tri (DMatrix t) (m,m) e -> IOMatrix (m,n) e -> IO ()
-trsm _ b
+
+trsm :: (BLAS3 e) => e -> Tri (DMatrix t) (m,m) e -> IOMatrix (m,n) e -> IO ()
+trsm _ _ b
     | M.numRows b == 0 || M.numCols b == 0 = return ()
-trsm t b =
-    let (u,d,alpha,a) = toBase t
+trsm alpha t b =
+    let (u,d,a)   = toBase t
         order     = colMajor
         (h,u')    = if isHerm a then (ConjTrans, flipUpLo u) else (NoTrans, u)
         (m,n)     = shape b
@@ -167,5 +281,6 @@
         ldB       = ldaOf b
     in M.unsafeWithElemPtr a (0,0) $ \pA ->
            M.unsafeWithElemPtr b (0,0) $ \pB -> do
-               BLAS.trsm order side uploA transA diagA m' n' (1/alpha') pA ldA pB ldB
+               BLAS.trsm order side uploA transA diagA m' n' alpha' pA ldA pB ldB
+
                
diff --git a/Data/Vector/Dense/IO.hs b/Data/Vector/Dense/IO.hs
--- a/Data/Vector/Dense/IO.hs
+++ b/Data/Vector/Dense/IO.hs
@@ -36,8 +36,6 @@
     -- * Conversion to and from @ForeignPtr@s
     fromForeignPtr,
     toForeignPtr,
-    isConj,
-    strideOf,
     
     -- * Converting between mutable and immutable vectors
     unsafeFreeze,
diff --git a/Data/Vector/Dense/Internal.hs b/Data/Vector/Dense/Internal.hs
--- a/Data/Vector/Dense/Internal.hs
+++ b/Data/Vector/Dense/Internal.hs
@@ -1,4 +1,5 @@
 {-# LANGUAGE BangPatterns, FlexibleInstances, MultiParamTypeClasses #-}
+{-# OPTIONS_GHC -fglasgow-exts #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module     : Data.Vector.Dense.Internal
@@ -20,8 +21,6 @@
     -- * Conversion to and from @ForeignPtr@s.
     fromForeignPtr,
     toForeignPtr,
-    isConj,
-    strideOf,
     
     -- * Creating new vectors
     newVector, 
@@ -76,12 +75,12 @@
 -- element type.  A @DVector@ @x@ stores @dim x@ elements.  Indices into
 -- the vector are @0@-based.
 data DVector t n e = 
-      DV { fptr   :: !(ForeignPtr e) -- ^ a pointer to the storage region
-         , offset :: !Int            -- ^ an offset (in elements, not bytes) to the first element in the vector. 
-         , len    :: !Int            -- ^ the length of the vector
-         , stride :: !Int            -- ^ the stride (in elements, not bytes) between elements.
+      DV { storageOf :: {-# UNPACK #-} !(ForeignPtr e) -- ^ a pointer to the storage region
+         , offsetOf  :: {-# UNPACK #-} !Int            -- ^ an offset (in elements, not bytes) to the first element in the vector. 
+         , lengthOf  :: {-# UNPACK #-} !Int            -- ^ the length of the vector
+         , strideOf  :: {-# UNPACK #-} !Int            -- ^ the stride (in elements, not bytes) between elements.
+         , isConj    :: {-# UNPACK #-} !Bool           -- ^ indicates whether or not the vector is conjugated
          }
-    | C !(DVector t n e)            -- ^ a conjugated vector
 
 type Vector n e = DVector Imm n e
 type IOVector n e = DVector Mut n e
@@ -89,38 +88,20 @@
 -- | Cast the phantom length type.
 coerceVector :: DVector t n e -> DVector t m e
 coerceVector = unsafeCoerce
-
--- | Gets the pointer to the storage block
-storageOf :: DVector t n e -> ForeignPtr e
-storageOf (C x)          = storageOf x
-storageOf x@(DV _ _ _ _) = fptr x
-
--- | Gets the stride of the vector.
-strideOf :: DVector t n e -> Int
-strideOf (C x)          = strideOf x
-strideOf x@(DV _ _ _ _) = stride x
-{-# INLINE strideOf #-}
-
--- | Indicates whether or not the vector has been conjugated.  For 
--- newly-created vectors, this will be @False@.
-isConj :: (Elem e) => DVector t n e -> Bool
-isConj (C x)        = not (isConj x)
-isConj (DV _ _ _ _) = False
-{-# INLINE isConj #-}
+{-# INLINE coerceVector #-}
 
--- | @fromForeignPtr fptr offset n inc@ creates a vector view of a
+-- | @fromForeignPtr fptr offset n inc c@ creates a vector view of a
 -- region in memory starting at the given offset and having dimension @n@,
--- with a stride of @inc@.
-fromForeignPtr :: ForeignPtr e -> Int -> Int -> Int -> DVector t n e
+-- with a stride of @inc@, and with @isConj@ set to @c@.
+fromForeignPtr :: ForeignPtr e -> Int -> Int -> Int -> Bool -> DVector t n e
 fromForeignPtr = DV
 {-# INLINE fromForeignPtr #-}
 
--- | Gets the tuple @(fptr,offset,n,inc)@, where @n@ is the dimension and 
--- @inc@ is the stride of the vector.  Note that this does not return the
--- conjugacy information of the vector.  For that information, use @isConj@.
-toForeignPtr :: DVector t n e -> (ForeignPtr e, Int, Int, Int)
-toForeignPtr (C x)        = toForeignPtr x
-toForeignPtr (DV f o n s) = (f, o, n, s)
+-- | Gets the tuple @(fptr,offset,n,inc,c)@, where @n@ is the dimension and 
+-- @inc@ is the stride of the vector, and @c@ indicates whether or not the
+-- vector is conjugated.
+toForeignPtr :: DVector t n e -> (ForeignPtr e, Int, Int, Int, Bool)
+toForeignPtr (DV f o n s c) = (f, o, n, s, c)
 {-# INLINE toForeignPtr #-}
 
 -- | @subvector x o n@ creates a subvector view of @x@ starting at index @o@ 
@@ -140,14 +121,14 @@
     
 -- | Same as 'subvectorWithStride' but arguments are not range-checked.
 unsafeSubvectorWithStride :: Int -> DVector t n e -> Int -> Int -> DVector t m e
-unsafeSubvectorWithStride s (C x)   o n = C   $ unsafeSubvectorWithStride s x o n
-unsafeSubvectorWithStride s x@(DV _ _ _ _) o n =
-    let f  = fptr x
+unsafeSubvectorWithStride s x o n =
+    let f  = storageOf x
         o' = indexOf x o
         n' = n
-        s' = s * (stride x)
+        s' = s * (strideOf x)
+        c  = isConj x
     in 
-        fromForeignPtr f o' n' s'
+        fromForeignPtr f o' n' s' c
 
 -- | Creates a new vector of the given length.  The elements will be 
 -- uninitialized.
@@ -158,7 +139,7 @@
             "Tried to create a vector with `" ++ show n ++ "' elements."
     | otherwise = do
         arr <- mallocForeignPtrArray n
-        return $ fromForeignPtr arr 0 n 1
+        return $ fromForeignPtr arr 0 n 1 False
 
 -- | Creates a new vector of the given dimension with the given elements.
 -- If the list has length less than the passed-in dimenson, the tail of
@@ -166,7 +147,7 @@
 newListVector :: (Elem e) => Int -> [e] -> IO (DVector t n e)
 newListVector n es = do
     x <- newVector_ n
-    withForeignPtr (fptr x) $ flip pokeArray $ take n es
+    withForeignPtr (storageOf x) $ flip pokeArray $ take n es
     return x
 
 -- | @listVector n es@ is equivalent to @vector n (zip [0..(n-1)] es)@, except
@@ -217,8 +198,7 @@
         unsafeWriteElem x i 1 
 
 indexOf :: DVector t n e -> Int -> Int
-indexOf (C x)   i         = indexOf x i
-indexOf x@(DV _ _ _ _) i  = offset x + i * stride x
+indexOf x i = offsetOf x + i * strideOf x
 {-# INLINE indexOf #-}
 
 -- | Evaluate a function with a pointer to the value stored at the given
@@ -240,28 +220,14 @@
 unsafeThaw = unsafeCoerce
 
 instance C.Vector (DVector t) where
-    dim x = case x of
-        (C x')   -> dim x'
-        _        -> len x
+    dim = lengthOf
     {-# INLINE dim #-}
 
-    conj x = case x of
-        (C x')   -> x'
-        _        -> C x
+    conj x = let c' = (not . isConj) x 
+             in x { isConj=c' }
     {-# INLINE conj #-}
 
-{-# RULES 
-    "conj/Float"  conj = conjFloat 
-    "conj/Double" conj = conjDouble
-  #-}
-conjFloat :: DVector t n Float -> DVector t n Float
-conjFloat = id
 
-conjDouble :: DVector t n Double -> DVector t n Double
-conjDouble = id
-
-
-
 instance Tensor (DVector t n) Int e where
     shape = dim
     {-# INLINE shape #-}
@@ -283,16 +249,6 @@
     
     amap f x = listVector (dim x) (map f $ elems x)
     
-    azipWith f x y
-        | dim y /= n =
-            error ("amap2: vector lengths differ; first has length `" ++
-                    show n ++ "' and second has length `" ++
-                    show (dim y) ++ "'")
-        | otherwise =
-            listVector n (zipWith f (elems x) (elems y))
-      where
-        n = dim x
-    
     (//) = replaceHelp writeElem
     unsafeReplace = replaceHelp unsafeWriteElem
 
@@ -314,13 +270,24 @@
     constant n e = unsafeFreeze $ unsafePerformIO $ newConstant n e
     {-# NOINLINE constant #-}
 
+    azipWith f x y
+        | dim y /= n =
+            error ("amap2: vector lengths differ; first has length `" ++
+                    show n ++ "' and second has length `" ++
+                    show (dim y) ++ "'")
+        | otherwise =
+            listVector n (zipWith f (elems x) (elems y))
+      where
+        n = dim x
+    
 
 instance (BLAS1 e) => RTensor (DVector t n) Int e IO where
     getSize = return . dim
     
-    newCopy x = case x of
-        (C x')   -> newCopy x' >>= return . C
-        _        -> do
+    newCopy x
+        | isConj x = 
+            newCopy (conj x) >>= return . conj
+        | otherwise = do
             y <- newVector_ (dim x)
             unsafeWithElemPtr x 0 $ \pX ->
                 unsafeWithElemPtr y 0 $ \pY ->
@@ -333,17 +300,21 @@
     getIndices = return . indices . unsafeFreeze
     {-# INLINE getIndices #-}
     
-    unsafeReadElem x i = case x of
-        (C x')   -> unsafeReadElem x' i >>= return . E.conj
-        _        -> withForeignPtr (fptr x) $ \ptr ->
-                        peekElemOff ptr (indexOf x i) 
+    unsafeReadElem x i
+        | isConj x = 
+            unsafeReadElem (conj x) i >>= return . E.conj
+        | otherwise =
+            withForeignPtr (storageOf x) $ \ptr ->
+                peekElemOff ptr (indexOf x i) 
 
-    getAssocs x = case x of
-        (C x')   -> getAssocs x' >>= return . map (\(i,e) -> (i,E.conj e))
-        _        -> let (f,o,n,incX) = toForeignPtr x
-                        ptr = (unsafeForeignPtrToPtr f) `advancePtr` o
-                    in return $ go n f incX ptr 0
-          where
+    getAssocs x
+        | isConj x =
+            getAssocs (conj x) >>= return . map (\(i,e) -> (i,E.conj e))
+        | otherwise =
+            let (f,o,n,incX,_) = toForeignPtr x
+                ptr = (unsafeForeignPtrToPtr f) `advancePtr` o
+            in return $ go n f incX ptr 0
+      where
             go !n !f !incX !ptr !i 
                 | i >= n = 
                      -- This is very important since we are doing unsafe IO.
@@ -372,9 +343,9 @@
         | strideOf x == 1 = unsafeWithElemPtr x 0 $ flip clearArray (dim x)
         | otherwise       = setConstant 0 x
 
-    setConstant e x = case x of
-        (C x')   -> setConstant (E.conj e) x'
-        _        -> unsafeWithElemPtr x 0 $ go (dim x)
+    setConstant e x 
+        | isConj x  = setConstant (E.conj e) (conj x)
+        | otherwise = unsafeWithElemPtr x 0 $ go (dim x)
       where
         go !n !ptr | n <= 0 = return ()
                    | otherwise = let ptr' = ptr `advancePtr` (strideOf x)
@@ -382,17 +353,18 @@
                                  in poke ptr e >> 
                                     go n' ptr'
     
-    unsafeWriteElem x i e = case x of
-        (C x')   -> unsafeWriteElem x' i $ E.conj e
-        _        -> withForeignPtr (fptr x) $ \ptr -> 
-                        pokeElemOff ptr (indexOf x i) e
+    unsafeWriteElem x i e =
+        let e' = if isConj x then E.conj e else e
+        in withForeignPtr (storageOf x) $ \ptr -> 
+               pokeElemOff ptr (indexOf x i) e'
                         
     canModifyElem x i = return $ inRange (bounds x) i
     {-# INLINE canModifyElem #-}
     
-    modifyWith f x = case x of
-        (C x')   -> modifyWith (E.conj . f . E.conj) x'
-        _        -> withForeignPtr (fptr x) $ go (dim x)
+    modifyWith f x
+        | isConj x  = modifyWith (E.conj . f . E.conj) (conj x)
+        | otherwise = withForeignPtr (storageOf x) $ 
+                          \ptr -> go (dim x) (ptr `advancePtr` offsetOf x)
       where
         go !n !ptr | n <= 0 = return ()
                    | otherwise = do
@@ -403,10 +375,11 @@
     
 compareHelp :: (BLAS1 e) => 
     (e -> e -> Bool) -> Vector n e -> Vector n e -> Bool
-compareHelp cmp (C x) (C y) =
-    compareHelp cmp x y
-compareHelp cmp x y =
-    (dim x == dim y) && (and $ zipWith cmp (elems x) (elems y))
+compareHelp cmp x y
+    | isConj x && isConj y =
+        compareHelp cmp (conj x) (conj y)
+    | otherwise =
+        (dim x == dim y) && (and $ zipWith cmp (elems x) (elems y))
 
 instance (BLAS1 e, Eq e) => Eq (DVector Imm n e) where
     (==) = compareHelp (==)
@@ -416,7 +389,7 @@
     (~==) = compareHelp (~==)
 
 instance (BLAS1 e, Show e) => Show (DVector Imm n e) where
-    show x = case x of
-        (C x')   -> "conj (" ++ show x' ++ ")"
-        _        -> "listVector " ++ show (dim x) ++ " " ++ show (elems x)
+    show x
+        | isConj x  = "conj (" ++ show (conj x) ++ ")"
+        | otherwise = "listVector " ++ show (dim x) ++ " " ++ show (elems x)
     
diff --git a/Data/Vector/Dense/Operations.hs b/Data/Vector/Dense/Operations.hs
--- a/Data/Vector/Dense/Operations.hs
+++ b/Data/Vector/Dense/Operations.hs
@@ -37,6 +37,7 @@
     divide,
     
     -- ** Impure
+    getConj,
     getShifted,
     getScaled,
     getInvScaled,
@@ -54,13 +55,19 @@
     (-=),
     (*=),
     (//=),
+
+    -- * BLAS calls
+    axpy,
     
     -- * Unsafe operations
     unsafeCopyVector,
     unsafeSwapVectors,
-    
-    -- * BLAS calls
-    axpy,
+    unsafeGetDot,
+    unsafeAxpy,
+    unsafePlusEquals,
+    unsafeMinusEquals,
+    unsafeTimesEquals,
+    unsafeDivideEquals,
     
     ) where
   
@@ -87,36 +94,39 @@
 -- | @copyVector dst src@ replaces the elements in @dst@ with the values from @src@.
 -- This may result in a loss of precision.
 copyVector :: (BLAS1 e) => IOVector n e -> DVector t n e -> IO ()
-copyVector x y = checkVecVecOp "copyVector" (dim x) (dim y) >> unsafeCopyVector x y
+copyVector x y = 
+    checkVecVecOp "copyVector" (dim x) (dim y) $ unsafeCopyVector x y
 
 -- | Same as 'copyVector' but does not check the dimensions of the arguments.
 unsafeCopyVector :: (BLAS1 e) => IOVector n e -> DVector t n e -> IO ()
-unsafeCopyVector (C x) (C y) =
-    unsafeCopyVector x y
-unsafeCopyVector x@(DV _ _ _ _) y@(DV _ _ _ _) =
-    call2 BLAS.copy y x
-unsafeCopyVector x y = do
-    forM_ [0..(dim x - 1)] $ \i -> do
-        unsafeReadElem y i >>= unsafeWriteElem x i
-
+unsafeCopyVector x y 
+    | isConj x && isConj y =
+        unsafeCopyVector (conj x) (conj y)
+    | isConj x || isConj y =
+        forM_ [0..(dim x - 1)] $ \i -> do
+            unsafeReadElem y i >>= unsafeWriteElem x i
+    | otherwise =
+        call2 BLAS.copy y x
 
 -- | @swapVectors x y@ replaces the elements in @x@ with the values from @y@, and 
 -- replaces the elements in @y@ with the values from @x@.  This may result in 
 -- a loss of precision.
 swapVectors :: (BLAS1 e) => IOVector n e -> IOVector n e -> IO ()
-swapVectors x y = checkVecVecOp "swapVectors" (dim x) (dim y) >> unsafeSwapVectors x y
+swapVectors x y = 
+    checkVecVecOp "swapVectors" (dim x) (dim y) $ unsafeSwapVectors x y
 
 -- | Same as 'swap' but does not check the dimensions of the arguments.
 unsafeSwapVectors :: (BLAS1 e) => IOVector n e -> IOVector n e -> IO ()
-unsafeSwapVectors (C x) (C y) =
-    unsafeSwapVectors x y
-unsafeSwapVectors x@(DV _ _ _ _) y@(DV _ _ _ _) =
-    call2 BLAS.swap x y
-unsafeSwapVectors x y = do
-    forM_ [0..(dim x - 1)] $ \i -> do
-        tmp <- unsafeReadElem x i
-        unsafeReadElem y i >>= unsafeWriteElem x i
-        unsafeWriteElem y i tmp
+unsafeSwapVectors x y
+    | isConj x && isConj y =
+        unsafeSwapVectors (conj x) (conj y)
+    | isConj x || isConj y =
+        forM_ [0..(dim x - 1)] $ \i -> do
+            tmp <- unsafeReadElem x i
+            unsafeReadElem y i >>= unsafeWriteElem x i
+            unsafeWriteElem y i tmp
+    | otherwise =
+        call2 BLAS.swap x y
 
 -- | Gets the sum of the absolute values of the vector entries.
 getSumAbs :: (BLAS1 e) => DVector t n e -> IO Double
@@ -140,20 +150,37 @@
 
 -- | Computes the dot product of two vectors.
 getDot :: (BLAS1 e) => DVector s n e -> DVector t n e -> IO e
-getDot x y = checkVecVecOp "dot" (dim x) (dim y) >> unsafeGetDot x y
+getDot x y = checkVecVecOp "dot" (dim x) (dim y) $ unsafeGetDot x y
+{-# INLINE getDot #-}
 
 unsafeGetDot :: (BLAS1 e) => DVector s n e -> DVector t n e -> IO e
-unsafeGetDot x@(DV _ _ _ _) y@(DV _ _ _ _) =
-    call2 dotc x y
-unsafeGetDot (C x@(DV _ _ _ _)) (y@(DV _ _ _ _)) =
-    call2 dotu x y
-unsafeGetDot (x@(DV _ _ _ _)) (C y@(DV _ _ _ _)) =
-    call2 dotu x y >>= return . E.conj
-unsafeGetDot x@(DV _ _ _ _) (C (C y)) = 
-    unsafeGetDot x y
-unsafeGetDot (C x) y = 
-    unsafeGetDot x (conj y) >>= return . E.conj
+unsafeGetDot x y =
+    case (isConj x, isConj y) of
+        (False, False) -> call2 dotc x y
+        (True , False) -> call2 dotu x y
+        (False, True ) -> call2 dotu x y >>= return . E.conj
+        (True , True)  -> call2 dotc x y >>= return . E.conj
+{-# INLINE unsafeGetDot #-}
 
+unsafeGetDotDouble :: DVector s n Double -> DVector t n Double -> IO Double
+unsafeGetDotDouble x y = call2 dotc x y
+{-# INLINE unsafeGetDotDouble #-}
+
+{-# RULES "unsafeGetDot/Double" unsafeGetDot = unsafeGetDotDouble #-}
+
+
+-- | Create a new vector by taking the conjugate of a vector.  The
+-- resulting vector will have 'isConj' equal to 'False'.
+getConj :: (BLAS1 e) => DVector t n e -> IO (DVector r n e)
+getConj x 
+    | isConj x = do
+        y <- newCopy (conj x)
+        return (unsafeCoerce y)
+    | otherwise = do
+        y <- newCopy x
+        doConj (unsafeThaw y)
+        return (unsafeCoerce y)
+
 -- | Create a new vector by adding a value to every element in a vector.
 getShifted :: (BLAS1 e) => e -> DVector t n e -> IO (DVector r n e)
 getShifted k x = do
@@ -178,23 +205,24 @@
 
 -- | Computes the sum of two vectors.
 getSum :: (BLAS1 e) => e -> DVector s n e -> e -> DVector t n e -> IO (DVector r n e)
-getSum alpha x beta y = checkVecVecOp "getSum" (dim x) (dim y) >> unsafeGetSum alpha x beta y
+getSum alpha x beta y = checkVecVecOp "getSum" (dim x) (dim y) $ unsafeGetSum alpha x beta y
 
 unsafeGetSum :: (BLAS1 e) => e -> DVector s n e -> e -> DVector t n e -> IO (DVector r n e)
 unsafeGetSum 1 x beta y
     | beta /= 1 = unsafeGetSum beta y 1 x
-unsafeGetSum alpha (C x) beta y = do
-    s <- unsafeGetSum (E.conj alpha) x (E.conj beta) (conj y)
-    return (conj s)
-unsafeGetSum alpha x@(DV _ _ _ _) beta y = do
-    s <- newCopy y
-    scaleBy beta (unsafeThaw s)
-    axpy alpha x (unsafeThaw s)
-    return (unsafeCoerce s)
+unsafeGetSum alpha x beta y
+    | isConj x = do
+        s <- unsafeGetSum (E.conj alpha) (conj x) (E.conj beta) (conj y)
+        return (conj s)
+    | otherwise = do
+        s <- newCopy y
+        scaleBy beta (unsafeThaw s)
+        axpy alpha x (unsafeThaw s)
+        return (unsafeCoerce s)
             
 -- | Computes the difference of two vectors.
 getDiff :: (BLAS1 e) => DVector s n e -> DVector t n e -> IO (DVector r n e)
-getDiff x y = checkVecVecOp "getDiff" (dim x) (dim y) >> unsafeGetSum 1 x (-1) y
+getDiff x y = checkVecVecOp "getDiff" (dim x) (dim y) $ unsafeGetSum 1 x (-1) y
 
 -- | Computes the elementwise product of two vectors.
 getProduct :: (BLAS2 e) => DVector s n e -> DVector t n e -> IO (DVector r n e)
@@ -210,74 +238,79 @@
             
 -- | Add a value to every element in a vector.
 shiftBy :: (BLAS1 e) => e -> IOVector n e -> IO ()
-shiftBy alpha (C x) = shiftBy (E.conj alpha) x
-shiftBy alpha x     = modifyWith (alpha+) x
+shiftBy alpha x | isConj x  = shiftBy (E.conj alpha) (conj x)
+                | otherwise = modifyWith (alpha+) x
 
 -- | Scale every element in the vector by the given value, and return a view
 -- to the scaled vector.  See also 'scale'.
 scaleBy :: (BLAS1 e) => e -> IOVector n e -> IO ()
 scaleBy 1 _ = return ()
-scaleBy k (C x) = do
-    scaleBy (E.conj k) x
-scaleBy k x@(DV _ _ _ _) =
-    call (flip scal k) x
+scaleBy k x | isConj x  = scaleBy (E.conj k) (conj x)
+            | otherwise = call (flip scal k) x
 
 -- | Divide every element by a value.
 invScaleBy :: (BLAS1 e) => e -> IOVector n e -> IO ()
-invScaleBy k (C x) = invScaleBy (E.conj k) x
-invScaleBy k x     = modifyWith (/k) x
+invScaleBy k x | isConj x  = invScaleBy (E.conj k) (conj x)
+               | otherwise = modifyWith (/k) x
 
+axpy :: (BLAS1 e) => e -> DVector t n e -> IOVector n e -> IO ()
+axpy alpha x y = checkVecVecOp "axpy" (dim x) (dim y) $ unsafeAxpy alpha x y
+
+unsafeAxpy :: (BLAS1 e) => e -> DVector t n e -> IOVector n e -> IO ()
+unsafeAxpy alpha x y
+    | isConj y =
+        axpy (E.conj alpha) (conj x) (conj y)
+    | isConj x =
+        call2 (flip BLAS.acxpy alpha) x y
+    | otherwise =
+        call2 (flip BLAS.axpy alpha) x y
+
 -- | @y += x@ replaces @y@ by @y + x@.
 (+=) :: (BLAS1 e) => IOVector n e -> DVector t n e -> IO ()
-(+=) y x = checkVecVecOp "(+=)" (dim y) (dim x) >> axpy 1 x y
+(+=) y x = checkVecVecOp "(+=)" (dim y) (dim x) $ unsafePlusEquals y x
 
-axpy :: (BLAS1 e) => e -> DVector t n e -> IOVector n e -> IO ()
-axpy alpha x@(DV _ _ _ _) y@(DV _ _ _ _) =
-    call2 (flip BLAS.axpy alpha) x y
-axpy alpha (C x@(DV _ _ _ _))  y@(DV _ _ _ _) =
-    call2 (flip BLAS.acxpy alpha) x y
-axpy alpha (C (C x)) y = 
-    axpy alpha x y
-axpy alpha x (C y) =
-    axpy (E.conj alpha) (conj x) y
+unsafePlusEquals :: (BLAS1 e) => IOVector n e -> DVector t n e -> IO ()
+unsafePlusEquals y x = axpy 1 x y
 
 -- | @y -= x@ replaces @y@ by @y - x@.
 (-=) :: (BLAS1 e) => IOVector n e -> DVector t n e -> IO ()
-(-=) y x = checkVecVecOp "(-=)" (dim y) (dim x) >> axpy (-1) x y
+(-=) y x = checkVecVecOp "(-=)" (dim y) (dim x) $ unsafeMinusEquals y x
 
+unsafeMinusEquals :: (BLAS1 e) => IOVector n e -> DVector t n e -> IO ()
+unsafeMinusEquals y x = unsafeAxpy (-1) x y
+
 -- | @y *= x@ replaces @y@ by @x * y@, the elementwise product.
 (*=) :: (BLAS2 e) => IOVector n e -> DVector t n e -> IO ()
-(*=) y x = checkVecVecOp "(*=)" (dim y) (dim x) >> timesEquals y x
+(*=) y x = checkVecVecOp "(*=)" (dim y) (dim x) $ unsafeTimesEquals y x
 
-timesEquals :: (BLAS2 e) => IOVector n e -> DVector t n e -> IO ()
-timesEquals y@(DV _ _ _ _) x@(DV _ _ _ _) =
-    call2 (flip (tbmv T.colMajor T.upper T.noTrans T.nonUnit) 0) x y
-timesEquals y@(DV _ _ _ _) (C x@(DV _ _ _ _)) =
-    call2 (flip (tbmv T.colMajor T.upper T.conjTrans T.nonUnit) 0) x y    
-timesEquals y@(DV _ _ _ _) (C (C x)) = 
-    timesEquals y x
-timesEquals (C y) x =
-    timesEquals y (conj x)
+unsafeTimesEquals :: (BLAS2 e) => IOVector n e -> DVector t n e -> IO ()
+unsafeTimesEquals y x
+    | isConj y =
+        unsafeTimesEquals (conj y) (conj x)
+    | isConj x =
+        call2 (flip (tbmv T.colMajor T.upper T.conjTrans T.nonUnit) 0) x y    
+    | otherwise =
+        call2 (flip (tbmv T.colMajor T.upper T.noTrans T.nonUnit) 0) x y
 
 -- | @y //= x@ replaces @y@ by @y / x@, the elementwise ratio.
 (//=) :: (BLAS2 e) => IOVector n e -> DVector t n e -> IO ()
-(//=) y x = checkVecVecOp "(//=)" (dim y) (dim x) >> divideEquals y x
+(//=) y x = checkVecVecOp "(//=)" (dim y) (dim x) $ unsafeDivideEquals y x
 
-divideEquals :: (BLAS2 e) => IOVector n e -> DVector t n e -> IO ()
-divideEquals y@(DV _ _ _ _) x@(DV _ _ _ _) =
-    call2 (flip (tbsv T.colMajor T.upper T.noTrans T.nonUnit) 0) x y
-divideEquals y@(DV _ _ _ _) (C x@(DV _ _ _ _)) =
-    call2 (flip (tbsv T.colMajor T.upper T.conjTrans T.nonUnit) 0) x y
-divideEquals y@(DV _ _ _ _) (C (C x)) = 
-    divideEquals y x
-divideEquals (C y) x =
-    divideEquals y (conj x)
+unsafeDivideEquals :: (BLAS2 e) => IOVector n e -> DVector t n e -> IO ()
+unsafeDivideEquals y x
+    | isConj y =
+        unsafeDivideEquals (conj y) (conj x)
+    | isConj x =
+        call2 (flip (tbsv T.colMajor T.upper T.conjTrans T.nonUnit) 0) x y
+    | otherwise =
+        call2 (flip (tbsv T.colMajor T.upper T.noTrans T.nonUnit) 0) x y
                
 call :: (Elem e) => (Int -> Ptr e -> Int -> IO a) -> DVector t n e -> IO a
 call f x =
     let n    = dim x
         incX = strideOf x
     in unsafeWithElemPtr x 0 $ \pX -> f n pX incX
+{-# INLINE call #-}
 
 call2 :: (Elem e) => 
        (Int -> Ptr e -> Int -> Ptr e -> Int -> IO a) 
@@ -289,11 +322,12 @@
     in unsafeWithElemPtr x 0 $ \pX ->
            unsafeWithElemPtr y 0 $ \pY ->
                f n pX incX pY incY
+{-# INLINE call2 #-}
 
 binaryOp :: (BLAS1 e) => String -> (IOVector n e -> DVector t n e -> IO ()) 
     -> DVector s n e -> DVector t n e -> IO (DVector r n e)
 binaryOp name f x y =
-    checkVecVecOp name (dim x) (dim y) >> do
+    checkVecVecOp name (dim x) (dim y) $ do
         x' <- newCopy x >>= return . unsafeThaw
         f x' y
         return $! (unsafeCoerce x')
diff --git a/INSTALL b/INSTALL
--- a/INSTALL
+++ b/INSTALL
@@ -17,12 +17,12 @@
 
 Installing blas requires a working CBLAS to be installed on your system.  If you
 do not have a CBLAS installed and you are not too concerned about performance,
-probably the easiest one to install is the one that comes with the GNU 
-Scientific Library (GSL), available at http://www.gnu.org/software/gsl/.  If
-you care about performance, a better option is to use ATLAS, available at
+probably the easiest to install is the one that comes with the GNU Scientific
+Library (GSL), available at http://www.gnu.org/software/gsl/.  If you care about
+performance, a better option is to use ATLAS, available at 
 http://math-atlas.sourceforge.net/.  Other options include the Goto BLAS and
 Intel's MKL library.  If you are running Mac OS X, you can use vecLib, which
-is part of the Accelerate framework.
+is part of the Accelerate framework and is installed by default.
 
 
 II. CONFIGURING
@@ -75,7 +75,8 @@
 If you have verified that blas is installed correctly, you can try building and
 running the tests.  These are in the "tests" directory, and are in the form of
 QuickCheck properties.  To build and run them, just run "make" in the tests
-directory.
+directory.  You must have the "pqc" package installed (available on hackage)
+to build and run the tests.
 
 
 VI. CUSTOM CBLAS CONFIGURATIONS
diff --git a/NEWS b/NEWS
new file mode 100644
--- /dev/null
+++ b/NEWS
@@ -0,0 +1,24 @@
+
+Changes in 0.5:
+
+* Add Banded matrix data type, as well as Tri Banded and Herm Banded.
+
+* Add support for trapezoidal dense matrices (Tri Matrix (m,n) e, where
+  m is not the same as n).  Note that trapezoidal banded matrices are
+  *NOT* supported.
+
+* Add Diag matrix data type for diagonal matrices.
+
+* Add Perm matrix data type, for permutation matrices.
+
+* Enhance the RMatrix and RSolve type classes with an API that allows 
+  specifying where to store the result of a computation.
+  
+* Enhance the IMatrix, RMatrix, ISolve, and RSolve type classes to add
+  "scale and multiply" operations.
+  
+* Remove the scale parameter for Tri and Herm matrix data types.
+
+* Flatten the data types for DVector and DMatrix.
+
+* Some inlining and unpacking performance improvements.
diff --git a/Setup.lhs b/Setup.lhs
--- a/Setup.lhs
+++ b/Setup.lhs
@@ -11,6 +11,16 @@
 >     system "runhaskell -lcblas -DCOMPLEX tests/HermMatrix.hs"
 >     system "runhaskell -lcblas -DREAL    tests/TriMatrix.hs"
 >     system "runhaskell -lcblas -DCOMPLEX tests/TriMatrix.hs"
+>     system "runhaskell -lcblas -DREAL    tests/Banded.hs"
+>     system "runhaskell -lcblas -DCOMPLEX tests/Banded.hs"
+>     system "runhaskell -lcblas -DREAL    tests/HermBanded.hs"
+>     system "runhaskell -lcblas -DCOMPLEX tests/HermBanded.hs"
+>     system "runhaskell -lcblas -DREAL    tests/TriBanded.hs"
+>     system "runhaskell -lcblas -DCOMPLEX tests/TriBanded.hs"
+>     system "runhaskell -lcblas -DREAL    tests/Perm.hs"
+>     system "runhaskell -lcblas -DCOMPLEX tests/Perm.hs"
+>     system "runhaskell -lcblas -DREAL    tests/Diag.hs"
+>     system "runhaskell -lcblas -DCOMPLEX tests/Diag.hs"
 >     return ()
 >
 > main = defaultMainWithHooks defaultUserHooks{ runTests=testing }
diff --git a/Test/QuickCheck/Matrix.hs b/Test/QuickCheck/Matrix.hs
--- a/Test/QuickCheck/Matrix.hs
+++ b/Test/QuickCheck/Matrix.hs
@@ -39,3 +39,10 @@
             return $ Assocs (m,n) (zip ijs' es)
         
     coarbitrary (Assocs mn ijes) = coarbitrary (mn,ijes)
+
+matrixSized :: (Int -> Gen a) -> Gen a
+matrixSized f = 
+    sized $ \s ->
+        let s' = ceiling (sqrt $ fromInteger $ toInteger s :: Double)
+        in f s'
+        
diff --git a/Test/QuickCheck/Matrix/Banded.hs b/Test/QuickCheck/Matrix/Banded.hs
new file mode 100644
--- /dev/null
+++ b/Test/QuickCheck/Matrix/Banded.hs
@@ -0,0 +1,177 @@
+{-# LANGUAGE FlexibleInstances #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module     : Test.QuickCheck.Matrix.Banded
+-- Copyright  : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+
+module Test.QuickCheck.Matrix.Banded
+    where
+
+import Control.Monad( forM )
+
+import Test.QuickCheck hiding ( vector )
+import qualified Test.QuickCheck as QC
+import Test.QuickCheck.Vector.Dense ( TestVector(..), dvector )
+import Test.QuickCheck.Matrix.Dense ( dmatrix )
+
+import Data.Vector.Dense ( Vector, dim )
+import Data.Matrix.Dense ( Matrix )
+import Data.Matrix.Banded
+import BLAS.Elem ( Elem, BLAS1 )
+
+newtype TestBanded m n e = TestBanded (Banded (m,n) e) deriving (Eq, Show)
+
+bmatrix :: (BLAS1 e, Arbitrary e) => 
+    (Int,Int) -> (Int,Int) -> Gen (Banded (m,n) e)
+bmatrix mn kl = frequency [ (3, rawBanded mn kl)  
+                          , (2, hermedBanded mn kl)
+                          ]
+
+rawBanded :: (BLAS1 e, Arbitrary e) => 
+    (Int,Int) -> (Int,Int) -> Gen (Banded (m,n) e)
+rawBanded (m,n) (k,l) = do
+    xs <- QC.vector ((k+1+l)*len)
+    return $ listsBanded (m,n) (k,l) (splitDiags xs)
+  where
+    len = min m n
+    
+    splitDiags [] = [[]]
+    splitDiags xs = let (d,xs') = splitAt len xs
+                    in d:(splitDiags xs')
+
+
+hermedBanded :: (BLAS1 e, Arbitrary e) => 
+    (Int,Int) -> (Int,Int) -> Gen (Banded (m,n) e)
+hermedBanded (m,n) (l,u) = do
+    x <- rawBanded (n,m) (u,l)
+    return $ (herm x)
+
+instance (Arbitrary e, BLAS1 e) => Arbitrary (TestBanded m n e) where
+    arbitrary = sized $ \k -> 
+        let k' = ceiling (sqrt $ fromInteger $ toInteger k :: Double)
+        in do
+            m <- choose (0,k')
+            n <- choose (0,k')
+            kl <- if m == 0 then return 0 else choose (0,m-1)
+            ku <- if n == 0 then return 0 else choose (0,n-1)
+            bmatrix (m,n) (kl,ku) >>= return . TestBanded
+        
+    coarbitrary (TestBanded x) =
+        coarbitrary (elems x)
+
+data BandedAt m n e = BandedAt (Banded (m,n) e) (Int,Int) deriving (Eq, Show)
+instance (Arbitrary e, BLAS1 e) => Arbitrary (BandedAt m n e) where
+    arbitrary = sized $ \k ->
+        let k' = ceiling (sqrt $ fromInteger $ toInteger k :: Double)
+        in do
+            m  <- choose (1,k'+1)
+            n  <- choose (1,k'+1)
+            kl <- choose (0,m-1)
+            ku <- choose (0,n-1)
+            i  <- choose (0,m-1)
+            j  <- choose (0,n-1)
+            a  <- bmatrix (m,n) (kl,ku)
+            
+            return $ BandedAt a (i,j)
+
+    coarbitrary = undefined
+
+
+data ListsBanded e = ListsBanded (Int,Int) (Int,Int) [[e]] deriving (Eq,Show)
+instance (Arbitrary e, Elem e) => Arbitrary (ListsBanded e) where
+    arbitrary = sized $ \k ->
+       let k' = ceiling (sqrt $ fromInteger $ toInteger k :: Double)
+       in do
+           m <- choose (1,k'+1)
+           n <- choose (1,k'+1)
+           kl <- choose (0,m-1)
+           ku <- choose (0,n-1)
+          
+           ds <- forM [(-kl)..ku] $ \i ->
+                     let beginPad = max (-i)    0
+                         endPad   = max (m-n+i) 0
+                         len      = m - (beginPad+endPad)
+                     in do
+                         xs <- QC.vector len
+                         return $ replicate beginPad 0 ++ xs ++ replicate endPad 0
+         
+           return $ ListsBanded (m,n) (kl,ku) ds
+          
+    coarbitrary (ListsBanded mn lu ds) = coarbitrary (mn,lu,ds)
+   
+{-
+                    
+data MatrixPair m n e = Pair (Matrix (m,n) e) (Matrix (m,n) e) deriving (Eq, Show)
+
+instance (Arbitrary e, Elem e) => Arbitrary (MatrixPair m n e) where
+    arbitrary = sized $ \k -> 
+        let k' = ceiling (sqrt $ fromInteger $ toInteger k :: Double)
+        in do
+            m <- choose (0,k')
+            n <- choose (0,k')
+            a <- dmatrix (m,n)
+            b <- dmatrix (m,n)
+            return $ Pair a b
+        
+    coarbitrary = undefined
+-}  
+  
+data BandedMV m n e = BandedMV (Banded (m,n) e) (Vector n e) deriving (Eq, Show)
+
+instance (Arbitrary e, BLAS1 e) => Arbitrary (BandedMV m n e) where
+    arbitrary = sized $ \k -> 
+        let k' = ceiling (sqrt $ fromInteger $ toInteger k :: Double)
+        in do
+            m <- choose (0,k')
+            n <- choose (0,k')
+            kl <- if m == 0 then return 0 else choose (0,m-1)
+            ku <- if n == 0 then return 0 else choose (0,n-1)
+            a <- bmatrix (m,n) (kl,ku)             
+            x <- dvector n
+            return $ BandedMV a x
+            
+    coarbitrary = undefined
+
+data BandedMVPair m n e = BandedMVPair (Banded (m,n) e) (Vector n e) (Vector n e) 
+    deriving (Eq, Show)
+    
+instance (Arbitrary e, BLAS1 e) => Arbitrary (BandedMVPair m n e) where
+    arbitrary = do
+        (BandedMV a x) <- arbitrary
+        y <- dvector (dim x)
+        return $ BandedMVPair a x y
+        
+    coarbitrary (BandedMVPair a x y) =
+        coarbitrary (BandedMV a x, TestVector y)
+        
+data BandedMM m n k e = BandedMM (Banded (m,k) e) (Matrix (k,n) e) deriving (Eq, Show)
+
+instance (Arbitrary e, BLAS1 e) => Arbitrary (BandedMM m n k e) where
+    arbitrary = sized $ \s ->
+        let s' = ceiling (sqrt $ fromInteger $ toInteger s :: Double)
+        in do
+            m <- choose (0,s')
+            k <- choose (0,s')
+            n <- choose (0,s')
+            kl <- if m == 0 then return 0 else choose (0,m-1)
+            ku <- if k == 0 then return 0 else choose (0,k-1)
+            a <- bmatrix (m,k) (kl,ku)             
+            b <- dmatrix (k,n)
+            return $ BandedMM a b
+            
+    coarbitrary = undefined
+        
+data BandedMMPair m n k e = BandedMMPair (Banded (m,k) e) (Matrix (k,n) e) (Matrix (k,n) e)
+    deriving (Eq, Show)
+    
+instance (Arbitrary e, BLAS1 e) => Arbitrary (BandedMMPair m n k e) where
+    arbitrary = do
+        (BandedMM a b) <- arbitrary
+        c <- dmatrix (shape b)
+        return $ BandedMMPair a b c
+        
+    coarbitrary = undefined
diff --git a/Test/QuickCheck/Matrix/Diag.hs b/Test/QuickCheck/Matrix/Diag.hs
new file mode 100644
--- /dev/null
+++ b/Test/QuickCheck/Matrix/Diag.hs
@@ -0,0 +1,100 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module     : Test.QuickCheck.Matrix.Diag
+-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+
+module Test.QuickCheck.Matrix.Diag
+    where
+
+import Test.QuickCheck hiding ( vector )
+import qualified Test.QuickCheck as QC
+
+import Test.QuickCheck.Vector.Dense ( dvector )
+import Test.QuickCheck.Matrix.Dense ( dmatrix )
+import Test.QuickCheck.Matrix ( matrixSized )
+
+import BLAS.Elem ( BLAS1, BLAS2 )
+
+import Data.Vector.Dense ( Vector )
+import Data.Matrix.Dense ( Matrix, matrix, diag )
+import Data.Matrix.Diag 
+
+diagMatrix :: (Arbitrary e, BLAS1 e) => Int -> Gen (Matrix (n,n) e)
+diagMatrix n = do
+    es <- QC.vector n
+    let a = matrix (n,n) [ ((i,i),e) | (i,e) <- zip [0..(n-1)] es ]
+    elements [ a, herm a ]
+
+
+data TestDiag n e = 
+    TestDiag (Diag (n,n) e) 
+             (Matrix (n,n) e) 
+    deriving Show
+
+instance (Arbitrary e, BLAS1 e) => Arbitrary (TestDiag n e) where
+    arbitrary = matrixSized $ \s -> do
+        n <- choose (0,s)
+        a <- diagMatrix n
+        return $ TestDiag (fromVector $ diag a 0) a
+        
+    coarbitrary = undefined
+
+
+data DiagMV n e = 
+    DiagMV (Diag (n,n) e)
+           (Matrix (n,n) e)
+           (Vector n e)
+    deriving Show
+    
+instance (Arbitrary e, BLAS1 e) => Arbitrary (DiagMV n e) where
+    arbitrary = do
+        (TestDiag d a) <- arbitrary
+        x <- dvector (numCols a)
+        return $ DiagMV d a x
+
+    coarbitrary = undefined
+
+
+data DiagMM m n e = 
+    DiagMM (Diag (m,m) e)
+           (Matrix (m,m) e)
+           (Matrix (m,n) e)
+    deriving Show
+
+instance (Arbitrary e, BLAS1 e) => Arbitrary (DiagMM m n e) where
+    arbitrary = matrixSized $ \s -> do
+        (TestDiag d a) <- arbitrary
+        n <- choose (0,s)
+        b <- dmatrix (numCols a, n)
+        return $ DiagMM d a b
+
+    coarbitrary = undefined
+    
+    
+data DiagSV n e =
+    DiagSV (Diag (n,n) e)
+           (Vector n e)
+    deriving Show
+
+instance (Arbitrary e, BLAS2 e) => Arbitrary (DiagSV n e) where
+    arbitrary = do
+        (DiagMV d _ x) <- arbitrary
+        return $ DiagSV d (d <*> x)
+        
+    coarbitrary = undefined
+
+data DiagSM m n e =
+    DiagSM (Diag (m,m) e)
+           (Matrix (m,n) e)
+    deriving Show
+    
+instance (Arbitrary e, BLAS2 e) => Arbitrary (DiagSM m n e) where
+    arbitrary = do
+        (DiagMM d _ b) <- arbitrary
+        return $ DiagSM d (d <**> b)
+
+    coarbitrary = undefined
diff --git a/Test/QuickCheck/Matrix/Herm/Banded.hs b/Test/QuickCheck/Matrix/Herm/Banded.hs
new file mode 100644
--- /dev/null
+++ b/Test/QuickCheck/Matrix/Herm/Banded.hs
@@ -0,0 +1,116 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module     : Test.QuickCheck.Matrix.Herm.Banded
+-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+
+module Test.QuickCheck.Matrix.Herm.Banded 
+    where
+
+import Control.Monad ( replicateM )
+
+import Test.QuickCheck hiding ( vector )
+import qualified Test.QuickCheck as QC
+import Test.QuickCheck.Vector.Dense ( dvector )
+import Test.QuickCheck.Matrix ( matrixSized )
+import Test.QuickCheck.Matrix.Dense ( dmatrix )
+
+import BLAS.Elem ( Elem, BLAS2, toReal, fromReal, conj )
+import BLAS.Types ( flipUpLo )
+
+import Data.Vector.Dense ( Vector )
+import Data.Matrix.Banded
+import Data.Matrix.Dense ( Matrix )
+import Data.Matrix.Herm
+
+
+
+hermBanded :: (BLAS2 e, Arbitrary e) => Int -> Int -> Gen (Banded (n,n) e)
+hermBanded n k 
+    | n < 0 = 
+        error $ "hermBanded: n must be non-negative"
+    | n == 0 =
+        return $ listsBanded (0,0) (0,0) []
+    | k >= n =
+        error $ "hermBanded: k must be less than n"
+    | k < 0 = 
+        error $ "hermBanded: k must be non-negative"
+    | k == 0 = do
+        d <- QC.vector n >>= return . realPart
+        return $ listsBanded (n,n) (0,0) [d]
+    | otherwise = do
+        a <- hermBanded n (k-1)
+        let (_,_,ds) = toLists a
+        
+        d <- QC.vector (n-k)
+        let d'  = map conj d
+            pad = replicate k 0
+            ds' = [pad ++ d] ++ ds ++ [d' ++ pad]
+
+        return $ listsBanded (n,n) (k,k) ds'
+        
+  where
+    realPart :: Elem e => [e] -> [e]
+    realPart = map (fromReal . toReal)
+
+data HermBanded n e =
+    HermBanded (Herm Banded (n,n) e)
+               (Banded (n,n) e)
+    deriving Show
+    
+instance (Arbitrary e, BLAS2 e) => Arbitrary (HermBanded n e) where
+    arbitrary = matrixSized $ \s -> do
+        n <- choose (0,s)
+        k <- if n == 0 then return 0 else choose (0,n-1)
+        l <- if n == 0 then return 0 else choose (0,n-1)
+            
+        a <- hermBanded n k
+            
+        junk <- replicateM l $ QC.vector n
+        let (_,_,ds) = toLists a
+            (u ,b ) = (Upper, listsBanded (n,n) (l,k) $ junk ++ (drop k ds))
+            (u',b') = (Lower, listsBanded (n,n) (k,l) $ (take (k+1) ds) ++ junk)
+        
+        h <- elements [ fromBase u             b
+                      , fromBase (flipUpLo u)  (herm b)
+                      , fromBase u'            b'
+                      , fromBase (flipUpLo u') (herm b')
+                      ]
+            
+        return $ HermBanded h a
+
+    coarbitrary = undefined
+
+data HermBandedMV n e = 
+    HermBandedMV (Herm Banded (n,n) e) 
+                 (Banded (n,n) e) 
+                 (Vector n e) 
+    deriving Show
+
+instance (Arbitrary e, BLAS2 e) => Arbitrary (HermBandedMV n e) where
+    arbitrary = do
+        (HermBanded h a) <- arbitrary
+        x <- dvector (numCols a)
+        return $ HermBandedMV h a x
+
+    coarbitrary = undefined
+    
+    
+data HermBandedMM m n e = 
+    HermBandedMM (Herm Banded (m,m) e) 
+                 (Banded (m,m) e) 
+                 (Matrix (m,n) e) 
+    deriving Show
+    
+instance (Arbitrary e, BLAS2 e) => Arbitrary (HermBandedMM m n e) where
+    arbitrary = matrixSized $ \s -> do
+        (HermBanded a h) <- arbitrary
+        n <- choose (0,s)
+        b <- dmatrix (numCols h,n)
+
+        return $ HermBandedMM a h b
+            
+    coarbitrary = undefined
diff --git a/Test/QuickCheck/Matrix/Herm/Dense.hs b/Test/QuickCheck/Matrix/Herm/Dense.hs
--- a/Test/QuickCheck/Matrix/Herm/Dense.hs
+++ b/Test/QuickCheck/Matrix/Herm/Dense.hs
@@ -13,9 +13,11 @@
 import Test.QuickCheck hiding ( vector )
 import qualified Test.QuickCheck as QC
 import Test.QuickCheck.Vector.Dense ( dvector )
+import Test.QuickCheck.Matrix ( matrixSized )
 import Test.QuickCheck.Matrix.Dense ( dmatrix, rawMatrix )
 
 import BLAS.Elem ( BLAS2 )
+import BLAS.Types ( flipUpLo )
 
 import Data.Vector.Dense
 import Data.Matrix.Dense
@@ -26,47 +28,61 @@
 hermMatrix :: (BLAS2 e, Arbitrary e) => Int -> Gen (Matrix (n,n) e)
 hermMatrix n  = do
     a <- rawMatrix (n,n)
-    return $ (a + herm a)
+    let h = (a + herm a)
+    elements [ h, herm h ]
 
+
+data HermMatrix n e = 
+    HermMatrix (Herm Matrix (n,n) e)
+               (Matrix (n,n) e)
+    deriving Show
+
+instance (Arbitrary e, BLAS2 e) => Arbitrary (HermMatrix n e) where
+    arbitrary = matrixSized $ \k -> do
+        n <- choose (0,k)
+        a <- hermMatrix n
+        
+        junk <- QC.vector (n*n)
+        let (u ,b ) = (Upper, a // zip (filter (uncurry (>)) $ indices a) junk)
+            (u',b') = (Lower, a // zip (filter (uncurry (<)) $ indices a) junk)
+
+        h <- elements [ fromBase u             b
+                      , fromBase (flipUpLo u)  (herm b)
+                      , fromBase u'            b'
+                      , fromBase (flipUpLo u') (herm b')
+                      ]
+        return $ HermMatrix h a
+        
+    coarbitrary = undefined
+
+
 data HermMatrixMV n e = 
-    HermMatrixMV UpLo (Matrix (n,n) e) (Matrix (n,n) e) (Vector n e) deriving (Show)
+    HermMatrixMV (Herm Matrix (n,n) e) 
+                 (Matrix (n,n) e) 
+                 (Vector n e) 
+    deriving Show
 
 instance (Arbitrary e, BLAS2 e) => Arbitrary (HermMatrixMV n e) where
-    arbitrary = sized $ \k ->
-        let k' = ceiling (sqrt $ fromInteger $ toInteger k :: Double)
-        in do
-            u <- elements [ Upper, Lower ]
-            n <- choose (0,k')
-            h <- hermMatrix n
-            let f  = case u of
-                         Upper -> \(i,j) -> i > j
-                         Lower -> \(i,j) -> i < j
-                zs = zip (filter f (indices h)) (repeat 0)
-                a  = h // zs
-                
-            x <- dvector n
-            return $ HermMatrixMV u a h x
+    arbitrary = do
+        (HermMatrix h a) <- arbitrary
+        x <- dvector (numCols a)
+        return $ HermMatrixMV h a x
+        
     coarbitrary = undefined
+
     
 data HermMatrixMM m n e = 
-    HermMatrixMM UpLo (Matrix (m,m) e) (Matrix (m,m) e) (Matrix (m,n) e) deriving (Show)
+    HermMatrixMM (Herm Matrix (m,m) e) 
+                 (Matrix (m,m) e) 
+                 (Matrix (m,n) e) 
+    deriving Show
     
 instance (Arbitrary e, BLAS2 e) => Arbitrary (HermMatrixMM m n e) where
-    arbitrary = sized $ \k ->
-        let k' = ceiling (sqrt $ fromInteger $ toInteger k :: Double)
-        in do
-            u <- elements [ Upper, Lower ]
-            m <- choose (0,k')
-            n <- choose (0,k')
-            h <- hermMatrix m
-            let f  = case u of
-                         Upper -> \(i,j) -> i > j
-                         Lower -> \(i,j) -> i < j
-                zs = zip (filter f (indices h)) (repeat 0)
-                a  = h // zs
-
-            b <- dmatrix (m,n)
-            return $ HermMatrixMM u a h b
+    arbitrary = matrixSized $ \k -> do
+        (HermMatrix h a) <- arbitrary
+        n <- choose (0,k)
+        b <- dmatrix (numCols a,n)
+        return $ HermMatrixMM h a b
             
     coarbitrary = undefined
         
diff --git a/Test/QuickCheck/Matrix/Perm.hs b/Test/QuickCheck/Matrix/Perm.hs
new file mode 100644
--- /dev/null
+++ b/Test/QuickCheck/Matrix/Perm.hs
@@ -0,0 +1,104 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module     : Test.QuickCheck.Matrix.Perm
+-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+
+module Test.QuickCheck.Matrix.Perm
+    where
+
+import Test.QuickCheck hiding ( vector )
+import qualified Test.QuickCheck as QC
+
+import Test.QuickCheck.Permutation
+import Test.QuickCheck.Vector.Dense ( dvector )
+import Test.QuickCheck.Matrix.Dense ( dmatrix )
+import Test.QuickCheck.Matrix ( matrixSized )
+
+import BLAS.Elem ( Elem, BLAS1 )
+
+import Data.Vector.Dense ( Vector )
+import Data.Matrix.Dense ( Matrix, shape )
+import Data.Matrix.Perm
+
+pmatrix :: (Elem e) => Int -> Gen (Perm (n,n) e)
+pmatrix n = frequency [ (2, rawPerm n)
+                      , (2, hermedPerm n)
+                      , (1, idPerm n)
+                      ]
+
+rawPerm :: Int -> Gen (Perm (n,n) e)
+rawPerm n = permutation n >>= return . fromPermutation
+
+hermedPerm :: (Elem e) => Int -> Gen (Perm (n,n) e)
+hermedPerm n = rawPerm n >>= return . herm
+
+idPerm :: Int -> Gen (Perm (n,n) e)
+idPerm = return . identity
+
+
+newtype TestPerm n e = TestPerm (Perm (n,n) e) deriving Show
+
+instance (Elem e) => Arbitrary (TestPerm n e) where
+    arbitrary = matrixSized $ \s -> do
+        n <- choose (0,s)
+        p <- pmatrix n
+        return $ TestPerm p
+        
+    coarbitrary = undefined
+
+data PermMBasis n e = PermMBasis (Perm (n,n) e) Int deriving Show
+
+instance (Elem e) => Arbitrary (PermMBasis n e) where
+    arbitrary = do
+        (TestPerm p) <- arbitrary
+        i <- choose (0, numCols p - 1)
+        return $ PermMBasis p i
+        
+    coarbitrary = undefined
+    
+data PermMV n e = PermMV (Perm (n,n) e) (Vector n e) deriving Show
+
+instance (BLAS1 e, Arbitrary e) => Arbitrary (PermMV n e) where
+    arbitrary = do
+        (TestPerm p) <- arbitrary
+        x <- dvector (numCols p)
+        return $ PermMV p x
+    
+    coarbitrary = undefined
+    
+data PermMVPair n e = 
+    PermMVPair (Perm (n,n) e) (Vector n e) (Vector n e) deriving Show
+    
+instance (BLAS1 e, Arbitrary e) => Arbitrary (PermMVPair n e) where
+    arbitrary = do
+        (PermMV p x) <- arbitrary
+        y <- dvector (numCols p)
+        return $ PermMVPair p x y
+        
+    coarbitrary = undefined
+    
+data PermMM m n e = PermMM (Perm (m,m) e) (Matrix (m,n) e) deriving Show
+
+instance (BLAS1 e, Arbitrary e) => Arbitrary (PermMM m n e) where
+    arbitrary = matrixSized $ \s -> do
+        (TestPerm p) <- arbitrary
+        n <- choose (0,s)
+        a <- dmatrix (numCols p, n)
+        return $ PermMM p a
+        
+    coarbitrary = undefined
+    
+data PermMMPair m n e = 
+    PermMMPair (Perm (m,m) e) (Matrix (m,n) e) (Matrix (m,n) e) deriving Show
+    
+instance (BLAS1 e, Arbitrary e) => Arbitrary (PermMMPair m n e) where
+    arbitrary = do
+        (PermMM p a)<- arbitrary
+        b <- dmatrix (shape a)
+        return $ PermMMPair p a b
+        
+    coarbitrary = undefined
diff --git a/Test/QuickCheck/Matrix/Tri/Banded.hs b/Test/QuickCheck/Matrix/Tri/Banded.hs
new file mode 100644
--- /dev/null
+++ b/Test/QuickCheck/Matrix/Tri/Banded.hs
@@ -0,0 +1,151 @@
+{-# LANGUAGE FlexibleInstances #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module     : Test.QuickCheck.Matrix.Tri.Banded
+-- Copyright  : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+
+module Test.QuickCheck.Matrix.Tri.Banded
+    where
+
+import Control.Monad ( replicateM )
+
+import Test.QuickCheck hiding ( vector )
+import qualified Test.QuickCheck as QC
+import Test.QuickCheck.Vector.Dense ( TestVector(..), dvector )
+import Test.QuickCheck.Matrix.Dense ( dmatrix )
+import Test.QuickCheck.Matrix ( matrixSized )
+
+import Data.Vector.Dense ( Vector )
+import Data.Matrix.Dense ( Matrix )
+import Data.Matrix.Banded
+import BLAS.Elem ( BLAS1, BLAS2 )
+
+import Data.Matrix.Tri.Banded ( Tri, UpLo(..), Diag(..), fromBase )
+
+triBanded :: (BLAS1 e, Arbitrary e) => UpLo -> Diag -> Int -> Int -> Gen (Banded (n,n) e)
+triBanded Upper NonUnit n k = do
+    a <- triBanded Upper Unit n k
+    d <- QC.vector n
+    let (_,_,(_:ds)) = toLists a
+    return $ listsBanded (n,n) (0,k) (d:ds)
+
+triBanded Lower NonUnit n k = do
+    a <- triBanded Lower Unit n k
+    d <- QC.vector n
+    let (_,_,ds) = toLists a
+        ds' = (init ds) ++ [d]
+    return $ listsBanded (n,n) (k,0) ds'
+    
+triBanded _ Unit n 0 = do
+    return $ listsBanded (n,n) (0,0) [replicate n 1]
+    
+triBanded Upper Unit n k = do
+    a <- triBanded Upper Unit n (k-1)
+    let (_,_,ds) = toLists a
+    
+    d <- QC.vector (n-k) >>= \xs -> return $ xs ++ replicate k 0
+    
+    return $ listsBanded (n,n) (0,k) $ ds ++ [d]
+    
+triBanded Lower Unit n k = do
+    a <- triBanded Lower Unit n (k-1)
+    let (_,_,ds) = toLists a
+
+    d <- QC.vector (n-k) >>= \xs -> return $ replicate k 0 ++ xs
+    
+    return $ listsBanded (n,n) (k,0) $ [d] ++ ds
+    
+    
+data TriBanded n e = 
+    TriBanded (Tri Banded (n,n) e) (Banded (n,n) e) deriving Show
+
+instance (Arbitrary e, BLAS1 e) => Arbitrary (TriBanded n e) where
+    arbitrary = matrixSized $ \s -> do
+        u <- elements [ Upper, Lower  ]
+        d <- elements [ Unit, NonUnit ]
+        n <- choose (0,s)
+        k <- if n == 0 then return 0 else choose (0,n-1)
+        a <- triBanded u d n k
+
+        l    <- if n == 0 then return 0 else choose (0,n-1)
+        junk <- replicateM l $ QC.vector n
+        diagJunk <- QC.vector n
+        let (_,_,ds) = toLists a
+            t = case (u,d) of 
+                    (Upper,NonUnit) -> 
+                        listsBanded (n,n) (l,k) $ junk ++ ds
+                    (Upper,Unit) ->
+                        listsBanded (n,n) (l,k) $ junk ++ [diagJunk] ++ tail ds
+                    (Lower,NonUnit) -> 
+                        listsBanded (n,n) (k,l) $ ds ++ junk
+                    (Lower,Unit) -> 
+                        listsBanded (n,n) (k,l) $ init ds ++ [diagJunk] ++ junk
+            t' = fromBase u d t 
+
+        return $ TriBanded t' a
+            
+    coarbitrary = undefined
+
+data TriBandedMV n e = 
+    TriBandedMV (Tri Banded (n,n) e) (Banded (n,n) e) (Vector n e) deriving Show
+
+instance (Arbitrary e, BLAS1 e) => Arbitrary (TriBandedMV n e) where
+    arbitrary = do
+        (TriBanded t a) <- arbitrary
+        x <- dvector (numCols t)
+        return $ TriBandedMV t a x
+        
+    coarbitrary (TriBandedMV t a x) =
+        coarbitrary (TriBanded t a, TestVector x)
+        
+data TriBandedMM m n e = 
+    TriBandedMM (Tri Banded (m,m) e) (Banded (m,m) e) (Matrix (m,n) e) deriving Show
+
+instance (Arbitrary e, BLAS1 e) => Arbitrary (TriBandedMM m n e) where
+    arbitrary = matrixSized $ \s -> do
+        (TriBanded t a) <- arbitrary
+        n <- choose (0,s)
+        b <- dmatrix (numCols t, n)
+        return $ TriBandedMM t a b
+            
+    coarbitrary = undefined
+        
+data TriBandedSV n e = 
+    TriBandedSV (Tri Banded (n,n) e) (Vector n e) deriving (Show)
+    
+instance (Arbitrary e, BLAS2 e) => Arbitrary (TriBandedSV n e) where
+    arbitrary = do
+        (TriBanded t _) <- arbitrary
+        t' <- elements [ t
+                       , herm t
+                       ]
+        x <- dvector (numCols t')
+        let y = t' <*> x
+        return (TriBandedSV t' y)
+        
+    coarbitrary = undefined
+
+
+data TriBandedSM m n e = 
+    TriBandedSM (Tri Banded (m,m) e) (Matrix (m,n) e) 
+    deriving (Show)
+    
+instance (Arbitrary e, BLAS2 e) => Arbitrary (TriBandedSM m n e) where
+    arbitrary = matrixSized $ \s -> do
+        (TriBandedSV t _) <- arbitrary
+        t' <- elements [ t
+                       , herm t
+                       ]
+
+        n <- choose (0, s)
+        a <- dmatrix (numCols t, n)
+        
+        let b = t' <**> a
+        return (TriBandedSM t' b)
+        
+    coarbitrary = undefined
+    
diff --git a/Test/QuickCheck/Matrix/Tri/Dense.hs b/Test/QuickCheck/Matrix/Tri/Dense.hs
--- a/Test/QuickCheck/Matrix/Tri/Dense.hs
+++ b/Test/QuickCheck/Matrix/Tri/Dense.hs
@@ -15,118 +15,125 @@
 
 import Test.QuickCheck hiding ( vector )
 import qualified Test.QuickCheck as QC
-import Test.QuickCheck.Vector.Dense ( TestVector(..), dvector )
+import Test.QuickCheck.Vector.Dense ( dvector )
 import Test.QuickCheck.Matrix.Dense ( dmatrix )
+import Test.QuickCheck.Matrix ( matrixSized )
 
 import Data.Vector.Dense
 import Data.Matrix.Dense
 import BLAS.Elem ( BLAS1, BLAS3 )
 
-import Data.Matrix.Tri.Dense ( Tri, UpLo(..), Diag(..), upper, lower, upperU, lowerU )
+import Data.Matrix.Tri.Dense ( Tri, UpLo(..), Diag(..), fromBase )
 
-triMatrix :: (BLAS1 e, Arbitrary e) => UpLo -> Diag -> Int -> Gen (Matrix (n,n) e)
-triMatrix u d n = 
-    let nz = case d of
-                 NonUnit -> n * (n + 1) `div` 2
-                 Unit    -> n * (n - 1) `div` 2
-    in do
-        h <- arbitrary
-        let f = case (h,u,d) of
-                 (False, Upper, NonUnit) -> \(i,j) -> i <= j
-                 (False, Upper,    Unit) -> \(i,j) -> i <  j
-                 (False, Lower, NonUnit) -> \(i,j) -> i >= j
-                 (False, Lower,    Unit) -> \(i,j) -> i >  j
-                 (True,  Upper, NonUnit) -> \(i,j) -> i >= j
-                 (True,  Upper,    Unit) -> \(i,j) -> i >  j
-                 (True,  Lower, NonUnit) -> \(i,j) -> i <= j
-                 (True,  Lower,    Unit) -> \(i,j) -> i <  j
-            ijs = filter f $ range ((0,0), (n-1,n-1))
-        es <- QC.vector nz
-        let a = matrix (n,n) $ zip ijs es
+import Unsafe.Coerce
 
-        a' <- case h of
-                  False -> return a
-                  True  -> return (herm a)
-        return a'
+triMatrix :: (BLAS1 e, Arbitrary e) => UpLo -> Diag -> (Int,Int) -> Gen (Matrix (m,n) e)
+triMatrix u d (m,n) =
+    let ijs = filter (isTriIndex u d) $ range ((0,0), (m-1,n-1))
+    in do
+        es <- QC.vector (m*n)
+        let a = matrix (m,n) $ zip ijs es
+            a' = case d of
+                    NonUnit -> a
+                    Unit    -> a // [ ((i,i),1) | i <- [0..(mn-1)] ]
+        return $ a'
+  where
+    mn = min m n
 
+isTriIndex :: UpLo -> Diag -> (Int,Int) -> Bool
+isTriIndex Upper NonUnit (i,j) = i <= j
+isTriIndex Upper Unit    (i,j) = i <  j
+isTriIndex Lower NonUnit (i,j) = i >= j
+isTriIndex Lower Unit    (i,j) = i >  j
 
-data TriMatrix n e = TriMatrix UpLo Diag (Matrix (n,n) e) deriving (Eq, Show)
+-- | A triangular matrix and an equivalent dense matrix
+data TriMatrix m n e = 
+    TriMatrix (Tri Matrix (m,n) e) 
+              (Matrix (m,n) e) 
+    deriving Show
 
-instance (Arbitrary e, BLAS1 e) => Arbitrary (TriMatrix n e) where
-    arbitrary = sized $ \k ->
-        let k' = ceiling (sqrt $ fromInteger $ toInteger k :: Double)
-        in do
-            u <- elements [ Upper, Lower  ]
-            d <- elements [ Unit, NonUnit ]
-            n <- choose (0,k')
-            a <- triMatrix u d n
-            return $ TriMatrix u d a
+instance (Arbitrary e, BLAS1 e) => Arbitrary (TriMatrix m n e) where
+    arbitrary = matrixSized $ \k -> do
+        u <- elements [ Upper, Lower  ]
+        d <- elements [ Unit, NonUnit ]
+        m <- choose (0,k)
+        n <- choose (0,k)
+        a <- triMatrix u d (m,n)
+        
+        junk <- QC.vector (m*n)
+        let ijs = [ (i,j) | i <- [0..(m-1)]
+                          , j <- [0..(n-1)]
+                          , (not . (isTriIndex u d)) (i,j) ]
+            t   = fromBase u d $ a // zip ijs junk
             
+        (t',a') <- elements [ (t,a), unsafeCoerce (herm t, herm a) ]
+            
+        return $ TriMatrix t' $ submatrix a' (0,0) (shape t')
+            
     coarbitrary = undefined
 
-data TriMatrixMV n e = 
-    TriMatrixMV UpLo Diag (Matrix (n,n) e) (Vector n e) deriving (Eq, Show)
 
-instance (Arbitrary e, BLAS1 e) => Arbitrary (TriMatrixMV n e) where
+-- | A triangular matrix, and equivalent dense matrix, and a vector in
+-- their domain.        
+data TriMatrixMV m n e = 
+    TriMatrixMV (Tri Matrix (m,n) e) 
+                (Matrix (m,n) e) 
+                (Vector n e) 
+    deriving Show
+
+instance (Arbitrary e, BLAS1 e) => Arbitrary (TriMatrixMV m n e) where
     arbitrary = do
-        (TriMatrix u d a) <- arbitrary
+        (TriMatrix t a) <- arbitrary
         x <- dvector (numCols a)
-        return $ TriMatrixMV u d a x
-        
-    coarbitrary (TriMatrixMV u d a x) =
-        coarbitrary (TriMatrix u d a, TestVector x)
-        
-data TriMatrixMM m n e = 
-    TriMatrixMM UpLo Diag (Matrix (m,m) e) (Matrix (m,n) e) deriving (Eq, Show)
+        return $ TriMatrixMV t a x
 
-instance (Arbitrary e, BLAS1 e) => Arbitrary (TriMatrixMM m n e) where
-    arbitrary = sized $ \k ->
-        let k' = ceiling (sqrt $ fromInteger $ toInteger k :: Double)
-        in do
-            (TriMatrix u d a) <- arbitrary
-            n <- choose (0,k')
-            b <- dmatrix (numCols a, n)
-            return $ TriMatrixMM u d a b
+    coarbitrary = undefined
+
+-- | A triangular matrix, and equivalent dense matrix, and a matrix in
+-- their domain.        
+data TriMatrixMM m k n e = 
+    TriMatrixMM (Tri Matrix (m,k) e) 
+                (Matrix (m,k) e) 
+                (Matrix (k,n) e) 
+    deriving Show
+
+instance (Arbitrary e, BLAS1 e) => Arbitrary (TriMatrixMM m k n e) where
+    arbitrary = matrixSized $ \s -> do
+        (TriMatrix t a) <- arbitrary
+        n <- choose (0,s)
+        b <- dmatrix (numCols a, n)
+        return $ TriMatrixMM t a b
             
     coarbitrary = undefined
-        
-data TriMatrixSV n e = 
-    TriMatrixSV (Tri Matrix (n,n) e) (Vector n e) deriving (Show)
+
+-- | A triangular matrix and a vector in its range
+data TriMatrixSV m n e = 
+    TriMatrixSV (Tri Matrix (m,n) e) 
+                (Vector m e) 
+    deriving Show
     
-instance (Arbitrary e, BLAS3 e) => Arbitrary (TriMatrixSV n e) where
+instance (Arbitrary e, BLAS3 e) => Arbitrary (TriMatrixSV m n e) where
     arbitrary = do
-        (TriMatrix u d a) <- arbitrary
-        let t = case (u,d) of
-                    (Lower,NonUnit) -> lower  a
-                    (Lower,Unit)    -> lowerU a
-                    (Upper,NonUnit) -> upper  a
-                    (Upper,Unit)    -> upperU a
-        k <- arbitrary
-        t' <- elements [ t
-                       , k *> t
-                       ]
-        x <- dvector (numCols t')
-        let y = t' <*> x
-        return (TriMatrixSV t' y)
+        (TriMatrix t a) <- arbitrary
+        x <- dvector (numCols a)
+        let y = a <*> x
+        return (TriMatrixSV t y)
         
     coarbitrary = undefined
 
-
-data TriMatrixSM m n e = 
-    TriMatrixSM (Tri Matrix (m,m) e) (Matrix (m,n) e) 
-    deriving (Show)
+-- | A triangular matrix and a matrix in its range
+data TriMatrixSM m k n e = 
+    TriMatrixSM (Tri Matrix (m,k) e) 
+                (Matrix (m,n) e) 
+    deriving Show
     
-instance (Arbitrary e, BLAS3 e) => Arbitrary (TriMatrixSM m n e) where
-    arbitrary = sized $ \k ->
-        let k' = ceiling (sqrt $ fromInteger $ toInteger k :: Double)
-        in do
-            (TriMatrixSV t _) <- arbitrary
-            n <- choose (0, k')
-            a <- dmatrix (numCols t, n)
-            
-            let b = t <**> a
-            return (TriMatrixSM t b)
+instance (Arbitrary e, BLAS3 e) => Arbitrary (TriMatrixSM m k n e) where
+    arbitrary = matrixSized $ \s -> do
+        (TriMatrix t a) <- arbitrary
+        n <- choose (0, s)
+        b <- dmatrix (numCols a, n)
+        let c = a <**> b
+        return (TriMatrixSM t c)
         
     coarbitrary = undefined
-
     
diff --git a/Test/QuickCheck/Permutation.hs b/Test/QuickCheck/Permutation.hs
new file mode 100644
--- /dev/null
+++ b/Test/QuickCheck/Permutation.hs
@@ -0,0 +1,36 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module     : Test.QuickCheck.Permutation
+-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+
+module Test.QuickCheck.Permutation
+    where
+
+import Data.List ( sortBy )
+import Data.Ord ( comparing )
+
+import Test.QuickCheck
+
+import Data.Permutation ( Permutation )
+import qualified Data.Permutation as P
+
+permutation :: Int -> Gen Permutation
+permutation n = do
+    xs <- vector n :: Gen [Int]
+    let is = (snd . unzip) $ sortBy (comparing fst) $ zip xs [0..]
+    return $ P.permutation n is
+    
+newtype TestPermutation = TestPermutation Permutation deriving Show
+
+instance Arbitrary TestPermutation where
+    arbitrary = do
+        n <- arbitrary >>= return . abs
+        p <- permutation n
+        return $ TestPermutation p
+        
+    coarbitrary (TestPermutation p) = 
+            coarbitrary $ P.toList p
diff --git a/blas.cabal b/blas.cabal
--- a/blas.cabal
+++ b/blas.cabal
@@ -1,5 +1,5 @@
 name:            blas
-version:         0.4.1
+version:         0.5
 homepage:        http://stat.stanford.edu/~patperry/code/blas
 synopsis:        Bindings to the BLAS library
 description:
@@ -25,9 +25,15 @@
 tested-with:     GHC == 6.8.2
 
 extra-source-files:     INSTALL
+                        NEWS
                         tests/Makefile
+                        tests/Banded.hs
+                        tests/Diag.hs
                         tests/Matrix.hs
+                        tests/Perm.hs
                         tests/HermMatrix.hs
+                        tests/HermBanded.hs
+                        tests/TriBanded.hs
                         tests/TriMatrix.hs
                         tests/Vector.hs
 
@@ -87,17 +93,26 @@
                         BLAS.Types
                         BLAS.Vector
                         
-                        Data.Matrix.Banded.Internal
+                        Data.Matrix.Banded
+                            Data.Matrix.Banded.IO
+                            Data.Matrix.Banded.Internal
+                            Data.Matrix.Banded.Operations
                         
+                        Data.Matrix.Diag
+                        
                         Data.Matrix.Dense
                             Data.Matrix.Dense.IO
                             Data.Matrix.Dense.Internal
                             Data.Matrix.Dense.Operations
                         
+                        Data.Matrix.Perm
+                        
                         Data.Matrix.Herm
+                            Data.Matrix.Herm.Banded
                             Data.Matrix.Herm.Dense
                         
                         Data.Matrix.Tri
+                            Data.Matrix.Tri.Banded
                             Data.Matrix.Tri.Dense
                         
                         Data.Vector.Dense
@@ -106,11 +121,17 @@
                             Data.Vector.Dense.Operations
                         
                         Test.QuickCheck.Complex
+                        Test.QuickCheck.Permutation
                         Test.QuickCheck.Vector
                         Test.QuickCheck.Vector.Dense
                         Test.QuickCheck.Matrix
+                        Test.QuickCheck.Matrix.Banded
                         Test.QuickCheck.Matrix.Dense
+                        Test.QuickCheck.Matrix.Diag
+                        Test.QuickCheck.Matrix.Perm
+                        Test.QuickCheck.Matrix.Herm.Banded
                         Test.QuickCheck.Matrix.Herm.Dense
+                        Test.QuickCheck.Matrix.Tri.Banded
                         Test.QuickCheck.Matrix.Tri.Dense
                         
     other-modules:      BLAS.C.Double
@@ -120,7 +141,7 @@
     extensions:         BangPatterns, EmptyDataDecls, FlexibleContexts, 
                         FlexibleInstances, ForeignFunctionInterface, 
                         FunctionalDependencies, MultiParamTypeClasses
-    build-depends:      base, ieee, storable-complex, QuickCheck
+    build-depends:      base, ieee, permutation, storable-complex, QuickCheck
 
     
     if flag(atlas)
diff --git a/tests/Banded.hs b/tests/Banded.hs
new file mode 100644
--- /dev/null
+++ b/tests/Banded.hs
@@ -0,0 +1,277 @@
+{-# OPTIONS -fglasgow-exts -fno-excess-precision -cpp #-}
+-----------------------------------------------------------------------------
+-- |
+-- Copyright  : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+
+
+import qualified Data.Array as Array
+import Data.Ix   ( inRange, range )
+import Data.List ( nub, sortBy, foldl' )
+import Data.Ord  ( comparing )
+import System.Environment ( getArgs )
+import Test.QuickCheck.Parallel hiding ( vector )
+import qualified Test.QuickCheck as QC
+
+import Data.Matrix.Banded
+
+import BLAS.Elem ( BLAS1 )
+import qualified BLAS.Elem as E
+
+import Data.Complex ( Complex(..) )
+import Data.Vector.Dense hiding ( invScale )
+import Data.Matrix.Dense ( Matrix, identity )
+import qualified Data.Matrix.Dense as M
+
+import Data.AEq
+import Numeric.IEEE
+
+import Test.QuickCheck.Complex
+import Test.QuickCheck.Vector hiding ( Assocs )
+import Test.QuickCheck.Vector.Dense hiding ( Pair )
+import Test.QuickCheck.Matrix
+import Test.QuickCheck.Matrix.Banded
+        
+isUndefR x = isNaN x || isInfinite x
+isUndefC (x :+ y) = isUndefR x || isUndefR y
+        
+#ifdef COMPLEX
+field = "Complex Double"
+type E = Complex Double
+isUndef = isUndefC
+#else
+field = "Double"
+type E = Double
+isUndef = isUndefR
+#endif        
+
+type V = Vector Int E
+type M = Matrix (Int,Int) E
+type B = Banded (Int,Int) E
+
+instance (Arbitrary e, RealFloat e) => Arbitrary (Complex e) where
+    arbitrary   = arbitrary >>= \(TestComplex x) -> return x
+    coarbitrary = coarbitrary . TestComplex
+
+instance (Arbitrary e, BLAS1 e) => Arbitrary (Vector n e) where
+    arbitrary   = arbitrary >>= \(TestVector x) -> return x
+    coarbitrary = coarbitrary . TestVector
+
+instance (Arbitrary e, BLAS1 e) => Arbitrary (Banded (m,n) e) where
+    arbitrary   = arbitrary >>= \(TestBanded x) -> return x
+    coarbitrary = coarbitrary . TestBanded
+
+
+assocsEq :: (BLAS1 e, AEq e) => Banded (m,n) e -> [((Int,Int), e)] -> Bool
+assocsEq x ijes =
+    let ijs = fst $ unzip ijes
+    in filter (\(ij,e) -> ij `elem` ijs) (sortBy (comparing fst) $ assocs x) === sortBy (comparing fst) ijes
+       && (all (==0) $ map snd $ filter (\(ij,e) -> not $ ij `elem` ijs) $ assocs x)
+
+fromAssocs (Assocs mn ijes) =
+    let kl = foldl' max 0 $ map (\((i,j),_) -> (i-j)) ijes
+        ku = foldl' max 0 $ map (\((i,j),_) -> (j-i)) ijes
+    in banded mn (kl,ku) ijes 
+    
+prop_banded_shape a@(Assocs mn _) =
+    shape (fromAssocs a :: B) == mn
+    
+prop_banded_assocs a@(Assocs mn ijes) =
+    (fromAssocs a :: B) `assocsEq` ijes
+
+prop_listsBanded_shape (ListsBanded mn lu ds) =
+    shape (listsBanded mn lu ds :: B) == mn
+    
+prop_listsBanded_toLists (ListsBanded mn lu ds) =
+    toLists (listsBanded mn lu ds :: B) === (mn,lu,ds)
+
+
+
+prop_replace_elems (a :: B) (Assocs _ ijes) =
+    let ijes' = filter (\((i,j),_) -> i < numRows a 
+                                      && j < numCols a
+                                      && (i-j) <= numLower a
+                                      && (j-i) <= numUpper a) ijes
+        old = filter (\(ij,_) -> not $ ij `elem` (fst . unzip) ijes') $ assocs a
+        expected = sortBy (comparing fst) $ old ++ ijes'
+        actual   = sortBy (comparing fst) $ assocs (a // ijes')
+    in expected === actual
+
+
+prop_shape (a :: B) = 
+    shape a == (numRows a, numCols a)
+
+prop_bandwidth (a :: B) =
+    bandwidth a == ((negate . numLower) a, numUpper a)
+
+prop_size (a :: B) =
+    size a == (sum $ map diagLen (range $ bandwidth a))
+  where
+    (m,n) = shape a
+    diagLen i | i <= 0    = min (m+i) n
+              | otherwise = min m     (n-i) 
+
+prop_bounds (a :: B) =
+    bounds a == ((0,0), (numRows a - 1, numCols a - 1))
+    
+prop_at (BandedAt (a :: B) ij@(i,j)) =
+    (a!ij) === expected
+    
+  where
+    (_,(kl,ku),ds) = toLists a
+    d = ds !! (kl + (j-i))
+    expected = case undefined of
+                   _ | i - j > kl -> 0
+                   _ | j - i > ku -> 0
+                   _ | otherwise  -> d!!i
+    
+prop_row_dim (BandedAt (a :: B) (i,_)) =
+    dim (row a i) == numCols a
+prop_col_dim (BandedAt (a :: B) (_,j)) =
+    dim (col a j) == numRows a
+prop_rows_len (a :: B) =
+    length (rows a) == numRows a
+prop_cols_len (a :: B) =
+    length (cols a) == numCols a
+prop_rows_dims (a :: B) =
+    map dim (rows a) == replicate (numRows a) (numCols a)
+prop_cols_dims (a :: B) =
+    map dim (cols a) == replicate (numCols a) (numRows a)
+
+prop_indices_length (a :: B) =
+    length (indices a) == size a
+prop_indices_lower (a :: B) =
+    all (\(i,j) -> i - j <= numLower a) $ indices a
+prop_indices_upper (a :: B) =
+    all (\(i,j) -> j - i <= numUpper a) $ indices a
+
+prop_elems_length (a :: B) =
+    length (elems a) == size a
+
+prop_assocs (a :: B) = 
+    assocs a === zip (indices a) (elems a)
+prop_assocs_at (a :: B) =
+    all (\(ij,e) -> a!ij === e) $ assocs a
+
+prop_scale_elems (a :: B) k =
+    and $ zipWith (~==) (elems (k *> a)) (map (k*) (elems a))
+prop_herm_elem (BandedAt (a :: B) (i,j)) =
+    (herm a) ! (j,i) === E.conj (a!(i,j))
+prop_herm_scale (a :: B) k =
+    herm (k *> a) === (E.conj k) *> (herm a)
+
+prop_herm_shape (a :: B) =
+    shape (herm a) == (numCols a, numRows a)
+prop_herm_rows (a :: B) =
+    rows (herm a) === map conj (cols a)
+prop_herm_cols (a :: B) = 
+    cols (herm a) === map conj (rows a)
+
+prop_herm_herm (a :: B) =
+    herm (herm a) === a
+
+prop_diag_herm1 (BandedAt (a :: B) (k,_)) =
+    diag a (-k) === conj (diag (herm a) k)
+prop_diag_herm2 (BandedAt (a :: B) (_,k)) =
+    diag a k === conj (diag (herm a) (-k))
+
+
+prop_apply_basis (BandedAt (a :: B) (_,j)) =
+    a <*> (basis (numCols a) j :: V) ~== col a j
+    || (any isUndef $ elems a)
+prop_apply_herm_basis (BandedAt (a :: B) (i,_)) =
+    (herm a) <*> (basis (numRows a) i :: V) ~== conj (row a i)
+    || (any isUndef $ elems a)    
+prop_apply_scale k (BandedMV (a :: B) x) =
+    a <*> (k *> x) ~== k *> (a <*> x)
+prop_apply_linear (BandedMVPair (a :: B) x y) =
+    a <*> (x + y) ~== a <*> x + a <*> y
+
+prop_applyMat_scale_left (BandedMM (a:: B) b) k =
+    a <**> (k *> b) ~== k *> (a <**> b)    
+prop_applyMat_scale_right (BandedMM (a:: B) b) k =
+    (k *> a) <**> b ~== k *> (a <**> b)
+prop_applyMat_linear (BandedMMPair (a :: B) b c) =
+    (a <**> (b + c) ~== a <**> b + a <**> c)
+    || (any isUndef $ elems (a <**> b + a <**> c))
+prop_applyMat_cols (BandedMM (a :: B) b) =
+    M.cols (a <**> b) ~== map (a <*> ) (M.cols b)
+
+prop_scale k (a :: B) =
+    k *> a ~== amap (\e -> e * k) a
+prop_invScale k (a :: B) =
+    invScale k a ~== amap (\e -> e / k) a
+
+
+properties =
+    [ ("shape of banded"       , pDet prop_banded_shape)
+    , ("assocs of banded"      , pDet prop_banded_assocs)
+    , ("shape of listsBanded"  , pDet prop_listsBanded_shape)
+    , ("listsBanded/toLists"   , pDet prop_listsBanded_toLists)
+
+    , ("elems of replace"      , pDet prop_replace_elems)
+    
+    , ("numRows/numCols"       , pDet prop_shape)
+    , ("numLower/numUpper"     , pDet prop_bandwidth)
+    , ("size"                  , pDet prop_size)
+    , ("bounds"                , pDet prop_bounds)
+
+    , ("at"                    , pDet prop_at)
+    , ("row dim"               , pDet prop_row_dim)
+    , ("col dim"               , pDet prop_col_dim)
+    , ("rows length"           , pDet prop_rows_len)
+    , ("cols length"           , pDet prop_cols_len)
+    , ("rows dims"             , pDet prop_rows_dims)
+    , ("cols dims"             , pDet prop_cols_dims)
+
+    , ("indices length"        , pDet prop_indices_length)
+    , ("indices low bw"        , pDet prop_indices_lower)
+    , ("indices up bw"         , pDet prop_indices_upper)
+
+    , ("elems length"          , pDet prop_elems_length)
+
+    , ("assocs"                , pDet prop_assocs)
+    , ("assocs/at"             , pDet prop_assocs_at)
+
+    , ("elems of scale"        , pDet prop_scale_elems)
+    , ("elem of herm"          , pDet prop_herm_elem)
+    , ("herm/scale"            , pDet prop_herm_scale)
+                               
+    , ("shape . herm"          , pDet prop_herm_shape)
+    , ("rows . herm"           , pDet prop_herm_rows)
+    , ("cols . herm"           , pDet prop_herm_cols)
+                               
+    , ("herm . herm == id"     , pDet prop_herm_herm)
+                               
+    , ("subdiag . herm"        , pDet prop_diag_herm1)
+    , ("superdiag . herm"      , pDet prop_diag_herm2)
+                               
+    , ("apply basis"           , pDet prop_apply_basis)
+    , ("apply herm basis"      , pDet prop_apply_herm_basis)
+    , ("apply scale"           , pDet prop_apply_scale)
+    , ("apply linear"          , pDet prop_apply_linear)
+    
+    , ("applyMat scale left"   , pDet prop_applyMat_scale_left)
+    , ("applyMat scale right"  , pDet prop_applyMat_scale_right)
+    , ("applyMat linear"       , pDet prop_applyMat_linear)
+    , ("applyMat cols"         , pDet prop_applyMat_cols)
+    
+    , ("scale"                 , pDet prop_scale)
+    , ("invScale"              , pDet prop_invScale)
+    
+    ]
+
+
+main = do
+    args <- getArgs
+    n <- case args of
+             (a:_) -> readIO a
+             _     -> return 1
+    main' n
+
+main' n = do
+    putStrLn $ "Running tests for " ++ field
+    pRun n 400 properties
diff --git a/tests/Diag.hs b/tests/Diag.hs
new file mode 100644
--- /dev/null
+++ b/tests/Diag.hs
@@ -0,0 +1,101 @@
+{-# OPTIONS -fglasgow-exts -fno-excess-precision -cpp #-}
+-----------------------------------------------------------------------------
+-- |
+-- Copyright  : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+
+import System.Environment ( getArgs )
+import Test.QuickCheck.Parallel hiding ( vector )
+import qualified Test.QuickCheck as QC
+
+import Data.Complex ( Complex(..) )
+
+import qualified BLAS.Elem as E
+import Data.Vector.Dense
+import Data.Matrix.Dense
+import Data.Matrix.Diag
+
+import Data.AEq
+import Numeric.IEEE
+
+import Debug.Trace
+
+import Test.QuickCheck.Complex
+import Test.QuickCheck.Matrix.Diag
+
+isUndefR x = isNaN x || isInfinite x
+isUndefC (x :+ y) = isUndefR x || isUndefR y
+        
+#ifdef COMPLEX
+field = "Complex Double"
+type E = Complex Double
+isUndef = isUndefC
+#else
+field = "Double"
+type E = Double
+isUndef = isUndefR
+#endif        
+
+type V = Vector Int E
+type M = Matrix (Int,Int) E
+type D = Diag (Int,Int) E
+
+instance (Arbitrary e, RealFloat e) => Arbitrary (Complex e) where
+    arbitrary   = arbitrary >>= \(TestComplex x) -> return x
+    coarbitrary = coarbitrary . TestComplex
+
+prop_diag_apply (DiagMV (d :: D) a x) =
+    d <*> x ~== a <*> x
+
+prop_diag_sapply k (DiagMV (d :: D) a x) =
+    sapply k d x ~== sapply k a x
+
+prop_diag_applyMat (DiagMM (d :: D) a b) =
+    d <**> b ~== a <**> b
+
+prop_diag_sapplyMat k (DiagMM (d :: D) a b) =
+    sapplyMat k d b ~== sapplyMat k a b
+
+
+prop_diag_solve (DiagSV (d :: D) y) =
+    let x = d <\> y
+    in d <*> x ~== y || (any isUndef $ elems x)
+
+prop_diag_ssolve k (DiagSV (d :: D) y) =
+    ssolve k d y ~== d <\> (k *> y)
+
+prop_diag_solveMat (DiagSM (d :: D) b) =
+    let a = d <\\> b
+    in d <**> a ~== b || (any isUndef $ elems a)
+
+prop_diag_ssolveMat k (DiagSM (d :: D) b) =
+    ssolveMat k d b ~== d <\\> (k *> b)
+
+
+properties =
+    [ ("diag apply"             , pDet prop_diag_apply)
+    , ("diag sapply"            , pDet prop_diag_sapply)
+    , ("diag applyMat"          , pDet prop_diag_applyMat)
+    , ("diag sapplyMat"         , pDet prop_diag_sapplyMat)
+
+    , ("diag solve"             , pDet prop_diag_solve)
+    , ("diag ssolve"            , pDet prop_diag_ssolve)
+    , ("diag solveMat"          , pDet prop_diag_solveMat)
+    , ("diag ssolveMat"         , pDet prop_diag_ssolveMat)
+    
+    ]
+
+
+main = do
+    args <- getArgs
+    n <- case args of
+             (a:_) -> readIO a
+             _     -> return 1
+    main' n
+
+main' n = do
+    putStrLn $ "Running tests for " ++ field
+    pRun n 200 properties
diff --git a/tests/HermBanded.hs b/tests/HermBanded.hs
new file mode 100644
--- /dev/null
+++ b/tests/HermBanded.hs
@@ -0,0 +1,82 @@
+{-# OPTIONS -fglasgow-exts -fno-excess-precision -cpp #-}
+-----------------------------------------------------------------------------
+-- |
+-- Copyright  : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+
+import Data.AEq
+import Numeric.IEEE
+
+import System.Environment ( getArgs )
+import Test.QuickCheck.Parallel hiding ( vector )
+import qualified Test.QuickCheck as QC
+
+
+import qualified BLAS.Elem as E
+import Data.Complex ( Complex(..) )
+
+import Data.Matrix.Herm.Banded
+import Data.Matrix.Banded
+import Data.Vector.Dense
+
+import Test.QuickCheck.Complex
+import Test.QuickCheck.Matrix.Herm.Banded
+
+#ifdef COMPLEX
+field = "Complex Double"
+type E = Complex Double
+#else
+field = "Double"
+type E = Double
+#endif
+
+type V = Vector Int E
+type B = Banded (Int,Int) E
+type HB = Herm Banded (Int,Int) E
+
+instance (Arbitrary e, RealFloat e) => Arbitrary (Complex e) where
+    arbitrary   = arbitrary >>= \(TestComplex x) -> return x
+    coarbitrary = coarbitrary . TestComplex
+
+prop_herm_apply (HermBandedMV (h :: HB) a x) =
+    h <*> x ~== a <*> x
+
+prop_herm_sapply k (HermBandedMV (h :: HB) a x) =
+    sapply k h x ~== sapply k a x
+
+prop_herm_herm_apply (HermBandedMV (h :: HB) a x) =
+    herm h <*> x ~== h <*> x
+
+prop_herm_applyMat (HermBandedMM (h :: HB) a b) =
+    h <**> b ~== a <**> b
+
+prop_herm_sapplyMat k (HermBandedMM (h :: HB) a b) =
+    sapplyMat k h b ~== sapplyMat k a b
+
+prop_herm_herm_applyMat (HermBandedMM (h :: HB) _ b) =
+    herm h <**> b ~== h <**> b
+
+
+properties =
+    [ ("herm apply"            , pDet prop_herm_apply)
+    , ("herm sapply"           , pDet prop_herm_sapply)
+    , ("herm herm apply"       , pDet prop_herm_herm_apply)
+
+    , ("herm applyMat"         , pDet prop_herm_applyMat)
+    , ("herm sapplyMat"        , pDet prop_herm_sapplyMat)
+    , ("herm herm applyMat"    , pDet prop_herm_herm_applyMat)
+    ]
+
+main = do
+    args <- getArgs
+    n <- case args of
+             (a:_) -> readIO a
+             _     -> return 1
+    main' n
+
+main' n = do
+    putStrLn $ "Running tests for " ++ field
+    pRun n 400 properties
diff --git a/tests/HermMatrix.hs b/tests/HermMatrix.hs
--- a/tests/HermMatrix.hs
+++ b/tests/HermMatrix.hs
@@ -33,54 +33,40 @@
 type E = Double
 #endif
 
-type V = Vector Int E
-type M = Matrix (Int,Int) E
+type V  = Vector Int E
+type M  = Matrix (Int,Int) E
+type HM = Herm Matrix (Int,Int) E
 
 instance (Arbitrary e, RealFloat e) => Arbitrary (Complex e) where
     arbitrary   = arbitrary >>= \(TestComplex x) -> return x
     coarbitrary = coarbitrary . TestComplex
 
-prop_herm_apply (HermMatrixMV u h (a :: M) x) =
-    case u of
-        Lower -> hermL h <*> x ~== a <*> x
-        Upper -> hermU h <*> x ~== a <*> x
-
-prop_scale_herm_apply k (HermMatrixMV u h (a :: M) x) =
-    case u of
-        Lower -> (k *> hermL h) <*> x ~== (k *> (a <*> x))
-        Upper -> (k *> hermU h) <*> x ~== (k *> (a <*> x))
+prop_herm_apply (HermMatrixMV (h :: HM) a x) =
+    h <*> x ~== a <*> x
 
-prop_herm_herm_apply (HermMatrixMV u h (a :: M) x) =
-    case u of
-        Lower -> hermU (herm h) <*> x ~== a <*> x
-        Upper -> hermL (herm h) <*> x ~== a <*> x
+prop_herm_sapply k (HermMatrixMV (h :: HM) a x) =
+    sapply k h x ~== sapply k a x
 
-prop_herm_compose (HermMatrixMM u h (a :: M) b) =
-    case u of
-        Lower -> hermL h <**> b ~== a <**> b
-        Upper -> hermU h <**> b ~== a <**> b
+prop_herm_herm_apply (HermMatrixMV (h :: HM) a x) =
+    herm h <*> x ~== h <*> x
 
-prop_scale_herm_compose k (HermMatrixMM u h (a :: M) b) =
-    case u of
-        Lower -> (k *> hermL h) <**> b ~== (k *> (a <**> b))
-        Upper -> (k *> hermU h) <**> b ~== (k *> (a <**> b))
+prop_herm_applyMat (HermMatrixMM (h :: HM) a b) =
+    h <**> b ~== a <**> b
 
-prop_herm_herm_compose (HermMatrixMM u h (a :: M) b) =
-    case u of
-        Lower -> hermU (herm h) <**> b ~== a <**> b
-        Upper -> hermL (herm h) <**> b ~== a <**> b
+prop_herm_sapplyMat k (HermMatrixMM (h :: HM) a b) =
+    sapplyMat k h b ~== sapplyMat k a b
 
+prop_herm_herm_applyMat (HermMatrixMM (h :: HM) _ b) =
+    herm h <**> b ~== h <**> b
 
 properties =
-    [ 
-      ("herm apply"            , pDet prop_herm_apply)
-    , ("scale herm apply"      , pDet prop_scale_herm_apply)
+    [ ("herm apply"            , pDet prop_herm_apply)
+    , ("herm sapply"           , pDet prop_herm_sapply)
     , ("herm herm apply"       , pDet prop_herm_herm_apply)
 
-    , ("herm compose"          , pDet prop_herm_compose)
-    , ("scale herm compose"    , pDet prop_scale_herm_compose)
-    , ("herm herm compose"     , pDet prop_herm_herm_compose)
-    
+    , ("herm applyMat"         , pDet prop_herm_applyMat)
+    , ("herm sapplyMat"        , pDet prop_herm_sapplyMat)
+    , ("herm herm applyMat"    , pDet prop_herm_herm_applyMat)
     ]
 
 
diff --git a/tests/Makefile b/tests/Makefile
--- a/tests/Makefile
+++ b/tests/Makefile
@@ -32,6 +32,40 @@
 	ghc --make -fforce-recomp -O -DCOMPLEX -threaded TriMatrix.hs -o tri-matrix-cplx
 	./tri-matrix-cplx $(THREADS)
 
+	ghc --make -fforce-recomp -O -DREAL -threaded Banded.hs -o banded-real
+	./banded-real $(THREADS)
+	
+	ghc --make -fforce-recomp -O -DCOMPLEX -threaded Banded.hs -o banded-cplx
+	./banded-cplx $(THREADS)
+
+	ghc --make -fforce-recomp -O -DREAL -threaded HermBanded.hs -o herm-banded-real
+	./herm-banded-real $(THREADS)
+	
+	ghc --make -fforce-recomp -O -DCOMPLEX -threaded HermBanded.hs -o herm-banded-cplx
+	./herm-banded-cplx $(THREADS)
+
+	ghc --make -fforce-recomp -O -DREAL -threaded TriBanded.hs -o tri-banded-real
+	./tri-banded-real $(THREADS)
+	
+	ghc --make -fforce-recomp -O -DCOMPLEX -threaded TriBanded.hs -o tri-banded-cplx
+	./tri-banded-cplx $(THREADS)
+
+	ghc --make -fforce-recomp -O -DREAL -threaded Perm.hs -o perm-real
+	./perm-real $(THREADS)
+	
+	ghc --make -fforce-recomp -O -DCOMPLEX -threaded Perm.hs -o perm-cplx
+	./perm-cplx $(THREADS)
+
+	ghc --make -fforce-recomp -O -DREAL -threaded Diag.hs -o diag-real
+	./diag-real $(THREADS)
+	
+	ghc --make -fforce-recomp -O -DCOMPLEX -threaded Diag.hs -o diag-cplx
+	./diag-cplx $(THREADS)
+
+
+
 clean:
 	rm -f *~ *.hi *.o vector-real vector-cplx matrix-real matrix-cplx \
-	herm-matrix-real herm-matrix-cplx tri-matrix-real tri-matrix-cplx
+	herm-matrix-real herm-matrix-cplx tri-matrix-real tri-matrix-cplx \
+	banded-real banded-cplx herm-banded-real herm-banded-cplx \
+	tri-banded-real tri-banded-cplx perm-real perm-cplx diag-real diag-cplx
diff --git a/tests/Perm.hs b/tests/Perm.hs
new file mode 100644
--- /dev/null
+++ b/tests/Perm.hs
@@ -0,0 +1,107 @@
+{-# OPTIONS -fglasgow-exts -fno-excess-precision -cpp #-}
+-----------------------------------------------------------------------------
+-- |
+-- Copyright  : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+
+import System.Environment ( getArgs )
+import Test.QuickCheck.Parallel hiding ( vector )
+import qualified Test.QuickCheck as QC
+
+import Data.Complex ( Complex(..) )
+
+import qualified BLAS.Elem as E
+import Data.Vector.Dense
+import Data.Matrix.Dense
+import Data.Matrix.Perm
+
+import Data.Permutation ( Permutation, permutation )
+import qualified Data.Permutation as P
+
+import Data.AEq
+import Numeric.IEEE
+
+import Test.QuickCheck.Complex
+import Test.QuickCheck.Matrix.Perm
+
+isUndefR x = isNaN x || isInfinite x
+isUndefC (x :+ y) = isUndefR x || isUndefR y
+        
+#ifdef COMPLEX
+field = "Complex Double"
+type E = Complex Double
+isUndef = isUndefC
+#else
+field = "Double"
+type E = Double
+isUndef = isUndefR
+#endif        
+
+type V = Vector Int E
+type M = Matrix (Int,Int) E
+type P = Perm (Int,Int) E
+
+instance (Arbitrary e, RealFloat e) => Arbitrary (Complex e) where
+    arbitrary   = arbitrary >>= \(TestComplex x) -> return x
+    coarbitrary = coarbitrary . TestComplex
+
+prop_perm_herm (TestPerm (p :: P)) =
+    toPermutation (herm p) == P.inverse (toPermutation p)
+
+prop_perm_apply_basis (PermMBasis (p :: P) i) =
+    n > 0 ==> p <*> (basis n i) === basis n (P.apply (toPermutation p) i)
+  where
+    n = numCols p
+
+prop_perm_herm_apply (PermMV (p :: P) x) =
+    p <*> herm p <*> x === x
+
+prop_herm_perm_apply (PermMV (p :: P) x) =
+    herm p <*> p <*> x === x
+
+prop_perm_solve (PermMV (p :: P) x) =
+    p <\> x === herm p <*> x
+
+prop_perm_applyMat_cols (PermMM (p :: P) a) =
+    cols (p <**> a) === map (p <*>) (cols a)
+
+prop_perm_herm_applyMat (PermMM (p :: P) a) =
+    p <**> herm p <**> a === a
+
+prop_herm_perm_applyMat (PermMM (p :: P) a) =
+    herm p <**> p <**> a === a
+
+prop_perm_solveMat_cols (PermMM (p :: P) a) =
+    cols (p <\\> a) === map (p <\>) (cols a)
+
+prop_perm_solveMat (PermMM (p :: P) a) =
+    p <\\> a === herm p <**> a
+    
+
+properties =
+    [ ("perm herm"             , pDet prop_perm_herm)
+    , ("perm apply basis"      , pDet prop_perm_apply_basis)
+    , ("perm herm apply"       , pDet prop_perm_herm_apply)
+    , ("herm perm apply"       , pDet prop_herm_perm_apply)
+    , ("perm solve"            , pDet prop_perm_solve)
+    , ("perm applyMat cols"    , pDet prop_perm_applyMat_cols)
+    , ("perm herm applyMat"    , pDet prop_perm_herm_applyMat)
+    , ("herm perm applyMat"    , pDet prop_herm_perm_applyMat)
+    , ("perm solveMat cols"    , pDet prop_perm_solveMat_cols)
+    , ("perm solveMat"         , pDet prop_perm_solveMat)
+    ]
+
+
+main = do
+    args <- getArgs
+    n <- case args of
+             (a:_) -> readIO a
+             _     -> return 1
+    main' n
+
+main' n = do
+    putStrLn $ "Running tests for " ++ field
+    pRun n 400 properties
diff --git a/tests/TriBanded.hs b/tests/TriBanded.hs
new file mode 100644
--- /dev/null
+++ b/tests/TriBanded.hs
@@ -0,0 +1,118 @@
+{-# OPTIONS -fglasgow-exts -fno-excess-precision -cpp #-}
+-----------------------------------------------------------------------------
+-- |
+-- Copyright  : Copyright (c) 2008, Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+
+import System.Environment ( getArgs )
+import Test.QuickCheck.Parallel hiding ( vector )
+import qualified Test.QuickCheck as QC
+
+import Data.Complex ( Complex(..) )
+
+import qualified BLAS.Elem as E
+import Data.Vector.Dense
+import Data.Matrix.Dense
+import Data.Matrix.Banded
+import Data.Matrix.Tri.Banded  
+
+import Data.AEq
+import Numeric.IEEE
+
+import Test.QuickCheck.Complex
+import Test.QuickCheck.Matrix.Tri.Banded
+
+isUndefR x = isNaN x || isInfinite x
+isUndefC (x :+ y) = isUndefR x || isUndefR y
+        
+#ifdef COMPLEX
+field = "Complex Double"
+type E = Complex Double
+isUndef = isUndefC
+#else
+field = "Double"
+type E = Double
+isUndef = isUndefR
+#endif        
+
+type V = Vector Int E
+type M = Matrix (Int,Int) E
+type B = Banded (Int,Int) E
+type TB = Tri Banded (Int,Int) E
+
+instance (Arbitrary e, RealFloat e) => Arbitrary (Complex e) where
+    arbitrary   = arbitrary >>= \(TestComplex x) -> return x
+    coarbitrary = coarbitrary . TestComplex
+
+prop_tri_apply (TriBandedMV (t :: TB) a x) =
+    t <*> x ~== a <*> x
+        
+prop_herm_tri_apply (TriBandedMV (t :: TB) a x) =
+    herm t <*> x ~== herm a <*> x
+
+prop_scale_tri_apply k (TriBandedMV (t :: TB) a x) =
+    sapply k t x ~== sapply k a x
+
+prop_scale_herm_tri_apply k (TriBandedMV (t :: TB) a x) =
+    sapply k (herm t) x ~== sapply k (herm a) x
+
+
+prop_tri_compose (TriBandedMM (t :: TB) a b) =
+    t <**> b ~== a <**> b
+
+prop_herm_tri_compose (TriBandedMM (t :: TB) a b) =
+    herm t <**> b ~== herm a <**> b
+
+prop_scale_tri_compose k (TriBandedMM (t :: TB) a b) =
+    sapplyMat k t b ~== sapplyMat k a b
+
+prop_scale_herm_tri_compose k (TriBandedMM (t :: TB) a b) =
+    sapplyMat k (herm t) b ~== sapplyMat k (herm a) b
+
+
+prop_tri_solve (TriBandedSV (t :: TB) y) =
+    let x = t <\> y
+    in t <*> x ~== y || (any isUndef $ elems x)
+
+prop_tri_ssolve k (TriBandedSV (t :: TB) y) =
+    ssolve k t y ~== t <\> (k*>y)
+
+prop_tri_solveMat (TriBandedSM (t :: TB) c) =
+    let b = t <\\> c
+    in t <**> b ~== c || (any isUndef $ elems b)
+
+prop_tri_ssolveMat k (TriBandedSM (t :: TB) c) =
+    ssolveMat k t c ~== t <\\> (k*>c)
+
+properties =
+    [ ("tri apply"             , pDet prop_tri_apply)
+    , ("scale tri apply"       , pDet prop_scale_tri_apply)
+    , ("herm tri apply"        , pDet prop_herm_tri_apply)
+    , ("scale herm tri apply"  , pDet prop_scale_herm_tri_apply)
+    
+    , ("tri compose"           , pDet prop_tri_compose)
+    , ("scale tri compose"     , pDet prop_scale_tri_compose)
+    , ("herm tri compose"      , pDet prop_herm_tri_compose)
+    , ("scale herm tri compose", pDet prop_scale_herm_tri_compose)
+    
+    , ("tri solve"             , pDet prop_tri_solve)
+    , ("tri solveMat"          , pDet prop_tri_solveMat)
+    , ("tri ssolve"            , pDet prop_tri_ssolve)
+    , ("tri ssolveMat"         , pDet prop_tri_ssolveMat)
+
+    ]
+
+
+main = do
+    args <- getArgs
+    n <- case args of
+             (a:_) -> readIO a
+             _     -> return 1
+    main' n
+
+main' n = do
+    putStrLn $ "Running tests for " ++ field
+    pRun n 400 properties
diff --git a/tests/TriMatrix.hs b/tests/TriMatrix.hs
--- a/tests/TriMatrix.hs
+++ b/tests/TriMatrix.hs
@@ -45,102 +45,44 @@
     arbitrary   = arbitrary >>= \(TestComplex x) -> return x
     coarbitrary = coarbitrary . TestComplex
 
-prop_tri_apply (TriMatrixMV u d (a :: M) x) =
-    case (u,d) of
-        (Lower, NonUnit) -> (lower  a) <*> x ~== a <*> x
-        (Lower, Unit   ) -> (lowerU a) <*> x ~== a <*> x + x
-        (Upper, NonUnit) -> (upper  a) <*> x ~== a <*> x
-        (Upper, Unit   ) -> (upperU a) <*> x ~== a <*> x + x
-
-prop_herm_tri_apply (TriMatrixMV u d (a :: M) x) =
-    case (u,d) of
-        (Lower, NonUnit) -> (herm $ lower  a) <*> x ~== herm a <*> x
-        (Lower, Unit   ) -> (herm $ lowerU a) <*> x ~== herm a <*> x + x
-        (Upper, NonUnit) -> (herm $ upper  a) <*> x ~== herm a <*> x
-        (Upper, Unit   ) -> (herm $ upperU a) <*> x ~== herm a <*> x + x
-
-prop_scale_tri_apply k (TriMatrixMV u d (a :: M) x) =
-    case (u,d) of
-        (Lower, NonUnit) -> (k *> lower  a) <*> x ~== (k *> a) <*> x
-        (Lower, Unit   ) -> (k *> lowerU a) <*> x ~== (k *> a) <*> x + (k *> x)
-        (Upper, NonUnit) -> (k *> upper  a) <*> x ~== (k *> a) <*> x
-        (Upper, Unit   ) -> (k *> upperU a) <*> x ~== (k *> a) <*> x + (k *> x)
-
-prop_scale_herm_tri_apply k (TriMatrixMV u d (a :: M) x) =
-    case (u,d) of
-        (Lower, NonUnit) -> (k *> (herm $ lower  a)) <*> x ~== (k *> herm a) <*> x
-        (Lower, Unit   ) -> (k *> (herm $ lowerU a)) <*> x ~== (k *> herm a) <*> x + (k *> x)
-        (Upper, NonUnit) -> (k *> (herm $ upper  a)) <*> x ~== (k *> herm a) <*> x
-        (Upper, Unit   ) -> (k *> (herm $ upperU a)) <*> x ~== (k *> herm a) <*> x + (k *> x)
-
-prop_herm_scale_tri_apply k (TriMatrixMV u d (a :: M) x) =
-    case (u,d) of
-        (Lower, NonUnit) -> (herm $ k *> lower  a) <*> x ~== (herm $ k *> a) <*> x
-        (Lower, Unit   ) -> (herm $ k *> lowerU a) <*> x ~== (herm $ k *> a) <*> x + ((E.conj k) *> x)
-        (Upper, NonUnit) -> (herm $ k *> upper  a) <*> x ~== (herm $ k *> a) <*> x
-        (Upper, Unit   ) -> (herm $ k *> upperU a) <*> x ~== (herm $ k *> a) <*> x + ((E.conj k) *> x)
-
-
-prop_tri_compose (TriMatrixMM u d (a :: M) b) =
-    case (u,d) of
-        (Lower, NonUnit) -> (lower  a) <**> b ~== a <**> b
-        (Lower, Unit   ) -> (lowerU a) <**> b ~== a <**> b + b
-        (Upper, NonUnit) -> (upper  a) <**> b ~== a <**> b
-        (Upper, Unit   ) -> (upperU a) <**> b ~== a <**> b + b
-
-prop_herm_tri_compose (TriMatrixMM u d (a :: M) b) =
-    case (u,d) of
-        (Lower, NonUnit) -> (herm $ lower  a) <**> b ~== herm a <**> b
-        (Lower, Unit   ) -> (herm $ lowerU a) <**> b ~== herm a <**> b + b
-        (Upper, NonUnit) -> (herm $ upper  a) <**> b ~== herm a <**> b
-        (Upper, Unit   ) -> (herm $ upperU a) <**> b ~== herm a <**> b + b
+prop_tri_apply (TriMatrixMV (t :: TM) a x) =
+    t <*> x ~== a <*> x
 
-prop_scale_tri_compose k (TriMatrixMM u d (a :: M) b) =
-    case (u,d) of
-        (Lower, NonUnit) -> (k *> lower  a) <**> b ~== (k *> a) <**> b
-        (Lower, Unit   ) -> (k *> lowerU a) <**> b ~== (k *> a) <**> b + (k *> b)
-        (Upper, NonUnit) -> (k *> upper  a) <**> b ~== (k *> a) <**> b
-        (Upper, Unit   ) -> (k *> upperU a) <**> b ~== (k *> a) <**> b + (k *> b)
+prop_tri_sapply k (TriMatrixMV (t :: TM) a x) =
+    sapply k t x ~== sapply k a x
 
-prop_scale_herm_tri_compose k (TriMatrixMM u d (a :: M) b) =
-    case (u,d) of
-        (Lower, NonUnit) -> (k *> (herm $ lower  a)) <**> b ~== (k *> herm a) <**> b
-        (Lower, Unit   ) -> (k *> (herm $ lowerU a)) <**> b ~== (k *> herm a) <**> b + (k *> b)
-        (Upper, NonUnit) -> (k *> (herm $ upper  a)) <**> b ~== (k *> herm a) <**> b
-        (Upper, Unit   ) -> (k *> (herm $ upperU a)) <**> b ~== (k *> herm a) <**> b + (k *> b)
+prop_tri_applyMat (TriMatrixMM (t :: TM) a b) =
+    t <**> b ~== a <**> b
 
-prop_herm_scale_tri_compose k (TriMatrixMM u d (a :: M) b) =
-    case (u,d) of
-        (Lower, NonUnit) -> (herm $ k *> lower  a) <**> b ~== (herm $ k *> a) <**> b
-        (Lower, Unit   ) -> (herm $ k *> lowerU a) <**> b ~== (herm $ k *> a) <**> b + ((E.conj k) *> b)
-        (Upper, NonUnit) -> (herm $ k *> upper  a) <**> b ~== (herm $ k *> a) <**> b
-        (Upper, Unit   ) -> (herm $ k *> upperU a) <**> b ~== (herm $ k *> a) <**> b + ((E.conj k) *> b)
+prop_tri_sapplyMat k (TriMatrixMM (t :: TM) a b) =
+    sapplyMat k t b ~== sapplyMat k a b
 
 
 prop_tri_solve (TriMatrixSV (t :: TM) y) =
     let x = t <\> y
     in t <*> x ~== y || (any isUndef $ elems x)
 
-prop_tri_invCompose (TriMatrixSM (t :: TM) b) =
+prop_tri_ssolve k (TriMatrixSV (t :: TM) y) =
+    ssolve k t y ~== t <\> (k *> y)
+
+prop_tri_solveMat (TriMatrixSM (t :: TM) b) =
     let a = t <\\> b
     in t <**> a ~== b || (any isUndef $ elems a)
 
+prop_tri_ssolveMat k (TriMatrixSM (t :: TM) b) =
+    ssolveMat k t b ~== t <\\> (k *> b)
 
+
 properties =
     [ ("tri apply"             , pDet prop_tri_apply)
-    , ("scale tri apply"       , pDet prop_scale_tri_apply)
-    , ("herm tri apply"        , pDet prop_herm_tri_apply)
-    , ("scale herm tri apply"  , pDet prop_scale_herm_tri_apply)
-    , ("herm scale tri apply"  , pDet prop_herm_scale_tri_apply)
-    
-    , ("tri compose"           , pDet prop_tri_compose)
-    , ("scale tri compose"     , pDet prop_scale_tri_compose)
-    , ("herm tri compose"      , pDet prop_herm_tri_compose)
-    , ("scale herm tri compose", pDet prop_scale_herm_tri_compose)
-    , ("herm scale tri compose", pDet prop_herm_scale_tri_compose)
-    
+    , ("tri sapply"            , pDet prop_tri_sapply)
+    , ("tri applyMat"          , pDet prop_tri_applyMat)
+    , ("tri sapplyMat"         , pDet prop_tri_sapplyMat)
+
     , ("tri solve"             , pDet prop_tri_solve)
-    , ("tri invCompose"        , pDet prop_tri_invCompose)
+    , ("tri ssolve"            , pDet prop_tri_ssolve)
+    , ("tri solveMat"          , pDet prop_tri_solveMat)
+    , ("tri ssolveMat"         , pDet prop_tri_ssolveMat)
     
     ]
 
@@ -154,4 +96,4 @@
 
 main' n = do
     putStrLn $ "Running tests for " ++ field
-    pRun n 400 properties
+    pRun n 200 properties
diff --git a/tests/Vector.hs b/tests/Vector.hs
--- a/tests/Vector.hs
+++ b/tests/Vector.hs
@@ -7,7 +7,6 @@
 -- Stability  : experimental
 --
 
-
 import qualified Data.Array as Array
 import Data.Ix   ( inRange )
 import Data.List ( nub, sortBy )
@@ -96,11 +95,13 @@
     dim (subvectorWithStride s x o n) == n
     
 prop_subvectorWithStride_elems (SubVector s (x :: V) o n) =
-    elems (subvectorWithStride s x o n) === 
-        (map snd $ filter (\(i,_) -> (i - o >= 0) 
-                                  && ((i - o) `mod` s == 0) 
-                                  && ((i - o) `div` s <= n))
-                          (assocs x))
+    let expected = (map snd $ filter (\(i,_) -> (i - o >= 0) 
+                              && ((i - o) `mod` s == 0) 
+                              && ((i - o) `div` s < n))
+                      (assocs x))
+        actual = elems (subvectorWithStride s x o n) 
+    in expected === actual
+        
 
 prop_dim (x :: V) = 
     dim x == length (elems x)