* Add getSum/sum.  Put in re-write rules from scale and add to axpy.  Here
 is the old code and re-write rules:
 
    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 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)

    {-# RULES
    "scale/plus"   forall k l x y. plus (scale k x) (scale l y) = add k x l y
    "scale1/plus"  forall k x y.   plus (scale k x) y = add k x 1 y
    "scale2/plus"  forall k x y.   plus x (scale k y) = add 1 x k y

    "scale/minus"  forall k l x y. minus (scale k x) (scale l y) = 
                                       add k x (-l) y
    "scale1/minus" forall k x y.   minus (scale k x) y = add k x (-1) y
    "scale2/minus" forall k x y.   minus x (scale k y) = add 1 x (-k) y
      #-}

* Matrix creation functions are inconsistent in phantom types.  Some have
  pairs, some do not.

* Add getSubmatrix function to MMatrix.

* Add getDense/toDense to MMatrix/IMatrix

* Put better checks in swapElem for WriteTensor

* STDiag and IODiag instead of current type?

* Add "Mult" class and get rid of <**> ?

* Add "Solv" class and get rid of <\\> ?

* lazy getRows/getCols needs to be re-implemented

* Optimization: Put 'isReal' in Elem class and specialize swap, copy, axpy, dot.

* Optimization: better "swapElem" for vector and matrix

* Optimization: Put vector/matrix fields in the typeclass and inline them.  
Implementation for IOVector:

    fptrOfVector (DV f _ _ _ _) = f
    {-# INLINE fptrOfVector #-}

    offsetOfVector (DV _ o _ _ _) = o
    {-# INLINE offsetOfVector #-}
    
    dim (DV _ _ n _ _) = n
    {-# INLINE dim #-}
    
    stride (DV _ _ _ s _) = s
    {-# INLINE stride #-}
    
    isConj (DV _ _ _ _ c) = c
    {-# INLINE isConj #-}

* Optimization: The ReadTensor/WriteTensor functions for Banded matrices are
really terrible.  Re-write these to be more efficient.