diff --git a/lapack.cabal b/lapack.cabal
--- a/lapack.cabal
+++ b/lapack.cabal
@@ -1,5 +1,5 @@
 Name:             lapack
-Version:          0.2
+Version:          0.2.1
 License:          BSD3
 License-File:     LICENSE
 Author:           Henning Thielemann <haskell@henning-thielemann.de>
@@ -45,7 +45,7 @@
   default:     False
 
 Source-Repository this
-  Tag:         0.2
+  Tag:         0.2.1
   Type:        darcs
   Location:    http://hub.darcs.net/thielema/lapack/
 
@@ -58,7 +58,7 @@
     lapack-ffi >=0.0.1 && <0.1,
     blas-ffi >=0.0 && <0.2,
     netlib-ffi >=0.1.1 && <0.2,
-    comfort-array >=0.2 && <0.3,
+    comfort-array >=0.3 && <0.4,
     guarded-allocation >=0.0 && <0.1,
     boxes >=0.1.5 && <0.2,
     deepseq >=1.3 && <1.5,
@@ -132,6 +132,10 @@
     netlib-ffi,
     tfp,
     comfort-array,
+    data-ref >=0.0.1 && <0.1,
+    unique-logic-tf >=0.5.1 && <0.6,
+    random >=1.1 && <1.2,
+    quickcheck-transformer >=0.3 && <0.4,
     QuickCheck >=2.5 && <3,
     ChasingBottoms >=1.2.2 && <1.4,
     transformers,
@@ -157,6 +161,7 @@
     Test.Banded.Utility
     Test.Singular
     Test.Generator
+    Test.Logic
     Test.Format
     Test.Utility
 
diff --git a/src/Numeric/LAPACK/Format.hs b/src/Numeric/LAPACK/Format.hs
--- a/src/Numeric/LAPACK/Format.hs
+++ b/src/Numeric/LAPACK/Format.hs
@@ -21,7 +21,7 @@
 import qualified Type.Data.Num.Unary as Unary
 import Type.Data.Num (integralFromProxy)
 
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
 import Data.Array.Comfort.Storable (Array)
 
diff --git a/src/Numeric/LAPACK/Linear/LowerUpper.hs b/src/Numeric/LAPACK/Linear/LowerUpper.hs
--- a/src/Numeric/LAPACK/Linear/LowerUpper.hs
+++ b/src/Numeric/LAPACK/Linear/LowerUpper.hs
@@ -55,10 +55,10 @@
 import qualified Numeric.Netlib.Utility as Call
 import qualified Numeric.Netlib.Class as Class
 
-import qualified Data.Array.Comfort.Storable.Internal.Monadic as ArrayIO
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked.Monadic as ArrayIO
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 
 import Foreign.Marshal.Array (advancePtr)
 import Foreign.C.Types (CInt)
diff --git a/src/Numeric/LAPACK/Linear/Private.hs b/src/Numeric/LAPACK/Linear/Private.hs
--- a/src/Numeric/LAPACK/Linear/Private.hs
+++ b/src/Numeric/LAPACK/Linear/Private.hs
@@ -11,9 +11,9 @@
 import qualified Numeric.Netlib.Utility as Call
 import qualified Numeric.Netlib.Class as Class
 
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 
 import Foreign.Marshal.Alloc (alloca)
 import Foreign.C.Types (CInt)
diff --git a/src/Numeric/LAPACK/Matrix.hs b/src/Numeric/LAPACK/Matrix.hs
--- a/src/Numeric/LAPACK/Matrix.hs
+++ b/src/Numeric/LAPACK/Matrix.hs
@@ -64,7 +64,8 @@
           multiplyVector, multiply, multiplyVectorUnchecked)
 import Numeric.LAPACK.Matrix.Divide
          (Solve(solve), solveVector, Inverse(inverse))
-import Numeric.LAPACK.Matrix.Basic (transpose, scaleRows, scaleColumns)
+import Numeric.LAPACK.Matrix.Basic
+         (transpose, forceOrder, forceRowMajor, scaleRows, scaleColumns)
 import Numeric.LAPACK.Matrix.Private
          (Full, Tall, Wide, General, argGeneral, ZeroInt, zeroInt,
           mapExtent, fromFull, generalizeTall, generalizeWide)
@@ -79,9 +80,9 @@
 import qualified Numeric.Netlib.Class as Class
 
 import qualified Data.Array.Comfort.Boxed as BoxedArray
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 import Data.Array.Comfort.Shape ((:+:)((:+:)))
 
 import Foreign.Marshal.Array (copyArray, advancePtr, pokeArray)
@@ -468,33 +469,6 @@
    Array.mapShape
       (\shape -> (MatrixShape.fullHeight shape, MatrixShape.fullWidth shape)) .
    forceRowMajor
-
-forceRowMajor ::
-   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,
-    Class.Floating a) =>
-   Full vert horiz height width a ->
-   Full vert horiz height width a
-forceRowMajor (Array shape@(MatrixShape.Full order extent) x) =
-   case order of
-      RowMajor -> Array shape x
-      ColumnMajor ->
-         Array.unsafeCreate (MatrixShape.Full RowMajor extent) $ \yPtr ->
-         withForeignPtr x $ \xPtr -> do
-            let (height, width) = Extent.dimensions extent
-            let n = Shape.size width
-            let m = Shape.size height
-            copyTransposed n m xPtr n yPtr
-
-forceOrder ::
-   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,
-    Class.Floating a) =>
-   Order ->
-   Full vert horiz height width a ->
-   Full vert horiz height width a
-forceOrder order =
-   case order of
-      RowMajor -> forceRowMajor
-      ColumnMajor -> transpose . forceRowMajor . transpose
 
 {- |
 @adaptOrder x y@ contains the data of @y@ with the layout of @x@.
diff --git a/src/Numeric/LAPACK/Matrix/Banded/Basic.hs b/src/Numeric/LAPACK/Matrix/Banded/Basic.hs
--- a/src/Numeric/LAPACK/Matrix/Banded/Basic.hs
+++ b/src/Numeric/LAPACK/Matrix/Banded/Basic.hs
@@ -50,9 +50,9 @@
 import Type.Data.Num (integralFromProxy)
 import Type.Base.Proxy (Proxy(Proxy))
 
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 
 import Foreign.Marshal.Array (advancePtr)
 import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)
diff --git a/src/Numeric/LAPACK/Matrix/Banded/Linear.hs b/src/Numeric/LAPACK/Matrix/Banded/Linear.hs
--- a/src/Numeric/LAPACK/Matrix/Banded/Linear.hs
+++ b/src/Numeric/LAPACK/Matrix/Banded/Linear.hs
@@ -24,7 +24,7 @@
 import Type.Data.Num (integralFromProxy)
 
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 
 import System.IO.Unsafe (unsafePerformIO)
 
diff --git a/src/Numeric/LAPACK/Matrix/BandedHermitian/Basic.hs b/src/Numeric/LAPACK/Matrix/BandedHermitian/Basic.hs
--- a/src/Numeric/LAPACK/Matrix/BandedHermitian/Basic.hs
+++ b/src/Numeric/LAPACK/Matrix/BandedHermitian/Basic.hs
@@ -49,9 +49,9 @@
 import Type.Data.Num (integralFromProxy)
 import Type.Base.Proxy (Proxy(Proxy))
 
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 
 import Foreign.Marshal.Array (advancePtr)
 import Foreign.C.Types (CInt, CChar)
diff --git a/src/Numeric/LAPACK/Matrix/BandedHermitian/Eigen.hs b/src/Numeric/LAPACK/Matrix/BandedHermitian/Eigen.hs
--- a/src/Numeric/LAPACK/Matrix/BandedHermitian/Eigen.hs
+++ b/src/Numeric/LAPACK/Matrix/BandedHermitian/Eigen.hs
@@ -24,10 +24,10 @@
 import qualified Type.Data.Num.Unary as Unary
 import Type.Data.Num (integralFromProxy)
 
-import qualified Data.Array.Comfort.Storable.Internal.Monadic as ArrayIO
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked.Monadic as ArrayIO
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 
 import Foreign.C.Types (CInt, CChar)
 import Foreign.Ptr (Ptr, nullPtr)
diff --git a/src/Numeric/LAPACK/Matrix/BandedHermitianPositiveDefinite/Linear.hs b/src/Numeric/LAPACK/Matrix/BandedHermitianPositiveDefinite/Linear.hs
--- a/src/Numeric/LAPACK/Matrix/BandedHermitianPositiveDefinite/Linear.hs
+++ b/src/Numeric/LAPACK/Matrix/BandedHermitianPositiveDefinite/Linear.hs
@@ -26,9 +26,9 @@
 import Type.Data.Num (integralFromProxy)
 import Type.Base.Proxy (Proxy(Proxy))
 
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 
 import Foreign.ForeignPtr (withForeignPtr)
 
diff --git a/src/Numeric/LAPACK/Matrix/Basic.hs b/src/Numeric/LAPACK/Matrix/Basic.hs
--- a/src/Numeric/LAPACK/Matrix/Basic.hs
+++ b/src/Numeric/LAPACK/Matrix/Basic.hs
@@ -15,9 +15,9 @@
 import qualified Numeric.Netlib.Utility as Call
 import qualified Numeric.Netlib.Class as Class
 
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 
 import Foreign.ForeignPtr (withForeignPtr)
 import Foreign.Storable (poke, peek)
@@ -44,6 +44,33 @@
 flattenColumn :: General height () a -> Vector height a
 flattenColumn = Array.mapShape MatrixShape.fullHeight
 
+
+forceRowMajor ::
+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,
+    Class.Floating a) =>
+   Full vert horiz height width a ->
+   Full vert horiz height width a
+forceRowMajor (Array shape@(MatrixShape.Full order extent) x) =
+   case order of
+      RowMajor -> Array shape x
+      ColumnMajor ->
+         Array.unsafeCreate (MatrixShape.Full RowMajor extent) $ \yPtr ->
+         withForeignPtr x $ \xPtr -> do
+            let (height, width) = Extent.dimensions extent
+            let n = Shape.size width
+            let m = Shape.size height
+            Private.copyTransposed n m xPtr n yPtr
+
+forceOrder ::
+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,
+    Class.Floating a) =>
+   Order ->
+   Full vert horiz height width a ->
+   Full vert horiz height width a
+forceOrder order =
+   case order of
+      RowMajor -> forceRowMajor
+      ColumnMajor -> transpose . forceRowMajor . transpose
 
 
 scaleRows ::
diff --git a/src/Numeric/LAPACK/Matrix/Hermitian/Basic.hs b/src/Numeric/LAPACK/Matrix/Hermitian/Basic.hs
--- a/src/Numeric/LAPACK/Matrix/Hermitian/Basic.hs
+++ b/src/Numeric/LAPACK/Matrix/Hermitian/Basic.hs
@@ -46,9 +46,9 @@
 import qualified Numeric.Netlib.Utility as Call
 import qualified Numeric.Netlib.Class as Class
 
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 
 import Foreign.C.Types (CInt, CChar)
 import Foreign.ForeignPtr (withForeignPtr)
@@ -322,7 +322,7 @@
    evalContT $ do
       let n = Shape.size sh
       aPtr <- ContT $ withForeignPtr a
-      conjPtr <- conjugateToTemp (MatrixShape.triangleSize n) a
+      conjPtr <- conjugateToTemp (Shape.triangleSize n) a
       liftIO $ do
          unpack (flipOrder order) n conjPtr bPtr -- wrong
          unpack order n aPtr bPtr
diff --git a/src/Numeric/LAPACK/Matrix/Hermitian/Eigen.hs b/src/Numeric/LAPACK/Matrix/Hermitian/Eigen.hs
--- a/src/Numeric/LAPACK/Matrix/Hermitian/Eigen.hs
+++ b/src/Numeric/LAPACK/Matrix/Hermitian/Eigen.hs
@@ -10,7 +10,7 @@
 import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
 import Numeric.LAPACK.Matrix.Hermitian.Private (TakeDiagonal(..))
 import Numeric.LAPACK.Matrix.Shape.Private
-         (Order(RowMajor,ColumnMajor), uploFromOrder, triangleSize)
+         (Order(RowMajor,ColumnMajor), uploFromOrder)
 import Numeric.LAPACK.Vector (Vector)
 import Numeric.LAPACK.Scalar (RealOf)
 import Numeric.LAPACK.Private
@@ -21,10 +21,11 @@
 import qualified Numeric.Netlib.Utility as Call
 import qualified Numeric.Netlib.Class as Class
 
-import qualified Data.Array.Comfort.Storable.Internal.Monadic as ArrayIO
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked.Monadic as ArrayIO
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
+import Data.Array.Comfort.Shape (triangleSize)
 
 import Foreign.C.Types (CInt, CChar)
 import Foreign.Ptr (Ptr, nullPtr)
diff --git a/src/Numeric/LAPACK/Matrix/Hermitian/Linear.hs b/src/Numeric/LAPACK/Matrix/Hermitian/Linear.hs
--- a/src/Numeric/LAPACK/Matrix/Hermitian/Linear.hs
+++ b/src/Numeric/LAPACK/Matrix/Hermitian/Linear.hs
@@ -17,9 +17,9 @@
 
 import qualified Numeric.Netlib.Class as Class
 
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 
 import System.IO.Unsafe (unsafePerformIO)
 
diff --git a/src/Numeric/LAPACK/Matrix/HermitianPositiveDefinite/Linear.hs b/src/Numeric/LAPACK/Matrix/HermitianPositiveDefinite/Linear.hs
--- a/src/Numeric/LAPACK/Matrix/HermitianPositiveDefinite/Linear.hs
+++ b/src/Numeric/LAPACK/Matrix/HermitianPositiveDefinite/Linear.hs
@@ -15,8 +15,7 @@
 import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
 import Numeric.LAPACK.Matrix.Hermitian.Private (Determinant(..))
 import Numeric.LAPACK.Matrix.Triangular.Private (copyTriangleToTemp)
-import Numeric.LAPACK.Matrix.Shape.Private
-         (NonUnit(NonUnit), uploFromOrder, triangleSize)
+import Numeric.LAPACK.Matrix.Shape.Private (NonUnit(NonUnit), uploFromOrder)
 import Numeric.LAPACK.Matrix.Private (Conjugation(Conjugated))
 import Numeric.LAPACK.Linear.Private (solver)
 import Numeric.LAPACK.Scalar (RealOf, realPart)
@@ -26,9 +25,10 @@
 import qualified Numeric.Netlib.Utility as Call
 import qualified Numeric.Netlib.Class as Class
 
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
+import Data.Array.Comfort.Shape (triangleSize)
 
 import Foreign.ForeignPtr (withForeignPtr)
 
diff --git a/src/Numeric/LAPACK/Matrix/Multiply.hs b/src/Numeric/LAPACK/Matrix/Multiply.hs
--- a/src/Numeric/LAPACK/Matrix/Multiply.hs
+++ b/src/Numeric/LAPACK/Matrix/Multiply.hs
@@ -32,9 +32,9 @@
 import qualified Type.Data.Num.Unary as Unary
 import Type.Data.Num.Unary ((:+:))
 
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 
 import Foreign.ForeignPtr (withForeignPtr)
 
diff --git a/src/Numeric/LAPACK/Matrix/Private.hs b/src/Numeric/LAPACK/Matrix/Private.hs
--- a/src/Numeric/LAPACK/Matrix/Private.hs
+++ b/src/Numeric/LAPACK/Matrix/Private.hs
@@ -4,9 +4,9 @@
 import qualified Numeric.LAPACK.Matrix.Extent as Extent
 import Numeric.LAPACK.Matrix.Shape.Private (Order, flipOrder)
 
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 
 import Foreign.ForeignPtr (ForeignPtr)
 
diff --git a/src/Numeric/LAPACK/Matrix/Shape/Private.hs b/src/Numeric/LAPACK/Matrix/Shape/Private.hs
--- a/src/Numeric/LAPACK/Matrix/Shape/Private.hs
+++ b/src/Numeric/LAPACK/Matrix/Shape/Private.hs
@@ -16,14 +16,14 @@
 import Type.Base.Proxy (Proxy(Proxy))
 
 import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Shape (triangleSize, triangleRoot)
 
 import Control.DeepSeq (NFData, rnf)
 import Control.Applicative (Const(Const, getConst))
 
-import qualified Data.NonEmpty as NonEmpty
 import Data.Functor.Identity (Identity(Identity), runIdentity)
 import Data.List (tails)
-import Data.Tuple.HT (mapPair, swap, double)
+import Data.Tuple.HT (mapSnd, swap, double)
 import Data.Bool.HT (if')
 
 
@@ -100,11 +100,12 @@
    sizeOffset (Full RowMajor extent) =
       Shape.sizeOffset (Extent.dimensions extent)
    sizeOffset (Full ColumnMajor extent) =
-      Shape.sizeOffset (swap $ Extent.dimensions extent) . swap
+      mapSnd (.swap) $ Shape.sizeOffset (swap $ Extent.dimensions extent)
    uncheckedSizeOffset (Full RowMajor extent) =
       Shape.uncheckedSizeOffset (Extent.dimensions extent)
    uncheckedSizeOffset (Full ColumnMajor extent) =
-      Shape.uncheckedSizeOffset (swap $ Extent.dimensions extent) . swap
+      mapSnd (.swap) $
+      Shape.uncheckedSizeOffset (swap $ Extent.dimensions extent)
 
    inBounds (Full _ extent) = Shape.inBounds (Extent.dimensions extent)
 
@@ -295,18 +296,25 @@
    uncheckedOffset (Split _ ColumnMajor extent) =
       Shape.uncheckedOffset (swap $ Extent.dimensions extent) . swap . snd
 
-   sizeOffset sh@(Split _ order extent) (part,ix) =
-      if part == splitPart sh ix
-        then
-            case order of
-               RowMajor -> Shape.sizeOffset (Extent.dimensions extent) ix
-               ColumnMajor ->
-                  Shape.sizeOffset (swap $ Extent.dimensions extent) (swap ix)
-        else error "Shape.Split.sizeOffset: wrong matrix part"
+   sizeOffset sh@(Split _ order extent) =
+      let check (part,ix) a =
+            if part == splitPart sh ix
+              then a
+              else error "Shape.Split.sizeOffset: wrong matrix part"
+      in case order of
+            RowMajor ->
+               mapSnd (\getOffset (part,ix) -> check (part,ix) $ getOffset ix) $
+               Shape.sizeOffset (Extent.dimensions extent)
+            ColumnMajor ->
+               mapSnd
+                  (\getOffset (part,ix) ->
+                     check (part,ix) $ getOffset $ swap ix) $
+               Shape.sizeOffset (swap $ Extent.dimensions extent)
    uncheckedSizeOffset (Split _ RowMajor extent) =
-      Shape.uncheckedSizeOffset (Extent.dimensions extent) . snd
+      mapSnd (.snd) $ Shape.uncheckedSizeOffset (Extent.dimensions extent)
    uncheckedSizeOffset (Split _ ColumnMajor extent) =
-      Shape.uncheckedSizeOffset (swap $ Extent.dimensions extent) . swap . snd
+      mapSnd (.swap.snd) $
+      Shape.uncheckedSizeOffset (swap $ Extent.dimensions extent)
 
    inBounds sh@(Split _ _ extent) (part,ix) =
       Shape.inBounds (Extent.dimensions extent) ix
@@ -351,21 +359,13 @@
 
    indices (Hermitian order size) = triangleIndices order size
 
-   uncheckedOffset sh ix =
-      snd $ Shape.uncheckedSizeOffset sh ix
-
-   sizeOffset sh ix =
-      if Shape.inBounds sh ix
-        then Shape.uncheckedSizeOffset sh ix
-        else error "Shape.Hermitian.sizeOffset: wrong matrix part"
+   offset (Hermitian order size) = triangleOffset order size
+   uncheckedOffset (Hermitian order size) =
+      triangleUncheckedOffset order size
 
-   uncheckedSizeOffset (Hermitian order size) (rs,cs) =
-      let (n,r) = Shape.uncheckedSizeOffset size rs
-          c = Shape.uncheckedOffset size cs
-      in (triangleSize n,
-          case order of
-            RowMajor -> triangleOffset n (r,c)
-            ColumnMajor -> triangleSize c + r)
+   sizeOffset (Hermitian order size) = triangleSizeOffset order size
+   uncheckedSizeOffset (Hermitian order size) =
+      triangleUncheckedSizeOffset order size
 
    inBounds (Hermitian _ size) ix@(r,c) =
       Shape.inBounds (size,size) ix
@@ -570,29 +570,31 @@
          (map swap $ triangleIndices (flipOrder order) size)
          (triangleIndices order size)
 
-   uncheckedOffset sh ix =
-      snd $ Shape.uncheckedSizeOffset sh ix
-
-   sizeOffset sh ix =
-      if Shape.inBounds sh ix
-        then Shape.uncheckedSizeOffset sh ix
-        else error "Shape.Triangular.sizeOffset: wrong matrix part"
+   offset (Triangular _diag uplo order size) =
+      caseDiagUpLoSym uplo
+         (Shape.offset size . snd)
+         (triangleOffset order size)
+         (triangleOffset (flipOrder order) size . swap)
+         (triangleOffset order size)
+   uncheckedOffset (Triangular _diag uplo order size) =
+      caseDiagUpLoSym uplo
+         (Shape.offset size . snd)
+         (triangleUncheckedOffset order size)
+         (triangleUncheckedOffset (flipOrder order) size . swap)
+         (triangleUncheckedOffset order size)
 
-   uncheckedSizeOffset (Triangular _diag uplo order size) (rs,cs) =
-      let (n,r) = Shape.uncheckedSizeOffset size rs
-          c = Shape.uncheckedOffset size cs
-          triSize = triangleSize n
-      in case order of
-            RowMajor ->
-               caseDiagUpLoSym uplo (n,c)
-                  (triSize, triangleOffset n (r,c))
-                  (triSize, triangleSize r + c)
-                  (triSize, triangleOffset n (r,c))
-            ColumnMajor ->
-               caseDiagUpLoSym uplo (n,c)
-                  (triSize, triangleSize c + r)
-                  (triSize, triangleOffset n (c,r))
-                  (triSize, triangleSize c + r)
+   sizeOffset (Triangular _diag uplo order size) =
+      caseDiagUpLoSym uplo
+         (mapSnd (.snd) $ Shape.sizeOffset size)
+         (triangleSizeOffset order size)
+         (mapSnd (.swap) $ triangleSizeOffset (flipOrder order) size)
+         (triangleSizeOffset order size)
+   uncheckedSizeOffset (Triangular _diag uplo order size) =
+      caseDiagUpLoSym uplo
+         (mapSnd (.snd) $ Shape.uncheckedSizeOffset size)
+         (triangleUncheckedSizeOffset order size)
+         (mapSnd (.swap) $ triangleUncheckedSizeOffset (flipOrder order) size)
+         (triangleUncheckedSizeOffset order size)
 
    inBounds (Triangular _diag uplo _ size) ix@(r,c) =
       Shape.inBounds (size,size) ix
@@ -615,16 +617,6 @@
          (triangleIndexFromOffset order size k)
 
 
-triangleSize :: Int -> Int
-triangleSize n = div (n*(n+1)) 2
-
-triangleOffset :: Int -> (Int,Int) -> Int
-triangleOffset s (r,c) =
-   triangleSize s - triangleSize (s-r) + c-r
-
-triangleRoot :: Floating a => a -> a
-triangleRoot size = (sqrt (8*size+1)-1)/2
-
 triangleRootDouble :: Int -> Double
 triangleRootDouble = triangleRoot . fromIntegral
 
@@ -637,30 +629,48 @@
 
 triangleIndices ::
    (Shape.Indexed sh) => Order -> sh -> [(Shape.Index sh, Shape.Index sh)]
-triangleIndices RowMajor size =
-   let ixs = Shape.indices size
-   in  concat $ zipWith (\r cs -> map ((,) r) cs) ixs $ tails ixs
-triangleIndices ColumnMajor size =
-   let ixs = Shape.indices size
-   in concat $
-      zipWith (\rs c -> map (flip (,) c) rs)
-         (NonEmpty.tail $ NonEmpty.inits ixs) ixs
+triangleIndices RowMajor = Shape.indices . Shape.upperTriangular
+triangleIndices ColumnMajor = map swap . Shape.indices . Shape.lowerTriangular
 
+triangleOffset ::
+   (Shape.Indexed sh) => Order -> sh -> (Shape.Index sh, Shape.Index sh) -> Int
+triangleOffset order size =
+   case order of
+      RowMajor -> Shape.offset (Shape.upperTriangular size)
+      ColumnMajor -> Shape.offset (Shape.lowerTriangular size) . swap
+
+triangleUncheckedOffset ::
+   (Shape.Indexed sh) => Order -> sh -> (Shape.Index sh, Shape.Index sh) -> Int
+triangleUncheckedOffset order size =
+   case order of
+      RowMajor -> Shape.uncheckedOffset (Shape.upperTriangular size)
+      ColumnMajor -> Shape.uncheckedOffset (Shape.lowerTriangular size) . swap
+
+triangleSizeOffset ::
+   (Shape.Indexed sh) =>
+   Order -> sh -> (Int, (Shape.Index sh, Shape.Index sh) -> Int)
+triangleSizeOffset order size =
+   case order of
+      RowMajor -> Shape.sizeOffset (Shape.upperTriangular size)
+      ColumnMajor ->
+         mapSnd (.swap) $ Shape.sizeOffset (Shape.lowerTriangular size)
+
+triangleUncheckedSizeOffset ::
+   (Shape.Indexed sh) =>
+   Order -> sh -> (Int, (Shape.Index sh, Shape.Index sh) -> Int)
+triangleUncheckedSizeOffset order size =
+   case order of
+      RowMajor -> Shape.uncheckedSizeOffset (Shape.upperTriangular size)
+      ColumnMajor ->
+         mapSnd (.swap) $ Shape.uncheckedSizeOffset (Shape.lowerTriangular size)
+
 triangleIndexFromOffset ::
    (Shape.InvIndexed sh) =>
    Order -> sh -> Int -> (Shape.Index sh, Shape.Index sh)
-triangleIndexFromOffset order size k =
-   mapPair (Shape.indexFromOffset size, Shape.indexFromOffset size) $
+triangleIndexFromOffset order size =
    case order of
-      RowMajor ->
-         let n = Shape.size size
-             triSize = triangleSize n
-             rr = ceiling (triangleRootDouble (triSize-k))
-             r = n - rr
-         in (r, k+r - (triSize - triangleSize rr))
-      ColumnMajor ->
-         let c = floor (triangleRootDouble k)
-         in (k - triangleSize c, c)
+      RowMajor -> Shape.indexFromOffset (Shape.upperTriangular size)
+      ColumnMajor -> swap . Shape.indexFromOffset (Shape.lowerTriangular size)
 
 
 type UnaryProxy a = Proxy (Unary.Un a)
diff --git a/src/Numeric/LAPACK/Matrix/Square/Basic.hs b/src/Numeric/LAPACK/Matrix/Square/Basic.hs
--- a/src/Numeric/LAPACK/Matrix/Square/Basic.hs
+++ b/src/Numeric/LAPACK/Matrix/Square/Basic.hs
@@ -46,9 +46,9 @@
 import qualified Numeric.Netlib.Utility as Call
 import qualified Numeric.Netlib.Class as Class
 
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 
 import Foreign.ForeignPtr (withForeignPtr)
 import Foreign.Storable (Storable, peek, poke)
diff --git a/src/Numeric/LAPACK/Matrix/Square/Eigen.hs b/src/Numeric/LAPACK/Matrix/Square/Eigen.hs
--- a/src/Numeric/LAPACK/Matrix/Square/Eigen.hs
+++ b/src/Numeric/LAPACK/Matrix/Square/Eigen.hs
@@ -21,8 +21,8 @@
 import qualified Numeric.Netlib.Utility as Call
 import qualified Numeric.Netlib.Class as Class
 
-import qualified Data.Array.Comfort.Storable.Internal.Monadic as ArrayIO
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked.Monadic as ArrayIO
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
 import Data.Array.Comfort.Storable (Array)
 
diff --git a/src/Numeric/LAPACK/Matrix/Square/Linear.hs b/src/Numeric/LAPACK/Matrix/Square/Linear.hs
--- a/src/Numeric/LAPACK/Matrix/Square/Linear.hs
+++ b/src/Numeric/LAPACK/Matrix/Square/Linear.hs
@@ -21,9 +21,9 @@
 import qualified Numeric.Netlib.Utility as Call
 import qualified Numeric.Netlib.Class as Class
 
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 
 import System.IO.Unsafe (unsafePerformIO)
 
diff --git a/src/Numeric/LAPACK/Matrix/Symmetric/Private.hs b/src/Numeric/LAPACK/Matrix/Symmetric/Private.hs
--- a/src/Numeric/LAPACK/Matrix/Symmetric/Private.hs
+++ b/src/Numeric/LAPACK/Matrix/Symmetric/Private.hs
@@ -5,7 +5,7 @@
          (diagonalPointerPairs, columnMajorPointers, rowMajorPointers,
           forPointers, pack, unpackToTemp, copyTriangleToTemp)
 import Numeric.LAPACK.Matrix.Shape.Private
-         (Order(RowMajor,ColumnMajor), uploFromOrder, triangleSize)
+         (Order(RowMajor,ColumnMajor), uploFromOrder)
 import Numeric.LAPACK.Matrix.Private
          (Full, Conjugation(NonConjugated, Conjugated))
 import Numeric.LAPACK.Linear.Private (solver, withDeterminantInfo, withInfo)
@@ -18,6 +18,7 @@
 import qualified Numeric.Netlib.Class as Class
 
 import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Shape (triangleSize)
 
 import Foreign.Marshal.Array (advancePtr)
 import Foreign.C.Types (CInt)
diff --git a/src/Numeric/LAPACK/Matrix/Triangular/Basic.hs b/src/Numeric/LAPACK/Matrix/Triangular/Basic.hs
--- a/src/Numeric/LAPACK/Matrix/Triangular/Basic.hs
+++ b/src/Numeric/LAPACK/Matrix/Triangular/Basic.hs
@@ -23,6 +23,12 @@
    takeUpper,
    takeLower,
 
+   fromLowerRowMajor, toLowerRowMajor,
+   fromUpperRowMajor, toUpperRowMajor,
+   forceOrder, adaptOrder,
+
+   add, sub,
+
    Tri.PowerDiag,
    multiplyVector,
    square, squareGeneric,
@@ -55,9 +61,9 @@
 import qualified Numeric.Netlib.Utility as Call
 import qualified Numeric.Netlib.Class as Class
 
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 
 import Foreign.C.Types (CChar, CInt)
 import Foreign.ForeignPtr (withForeignPtr)
@@ -201,6 +207,100 @@
          (MatrixShape.Triangular MatrixShape.NonUnit
             MatrixShape.upper order width) $ \bPtr ->
       withForeignPtr a $ \aPtr -> packRect order n k aPtr bPtr
+
+
+fromLowerRowMajor ::
+   (Shape.C sh, Class.Floating a) =>
+   Array (Shape.Triangular Shape.Lower sh) a -> Lower sh a
+fromLowerRowMajor =
+   Array.mapShape
+      (MatrixShape.Triangular MatrixShape.NonUnit MatrixShape.lower RowMajor .
+       Shape.triangularSize)
+
+fromUpperRowMajor ::
+   (Shape.C sh, Class.Floating a) =>
+   Array (Shape.Triangular Shape.Upper sh) a -> Upper sh a
+fromUpperRowMajor =
+   Array.mapShape
+      (MatrixShape.Triangular MatrixShape.NonUnit MatrixShape.upper RowMajor .
+       Shape.triangularSize)
+
+toLowerRowMajor ::
+   (Shape.C sh, Class.Floating a) =>
+   Lower sh a -> Array (Shape.Triangular Shape.Lower sh) a
+toLowerRowMajor =
+   Array.mapShape (Shape.Triangular Shape.Lower . MatrixShape.triangularSize)
+   .
+   forceOrder MatrixShape.RowMajor
+
+toUpperRowMajor ::
+   (Shape.C sh, Class.Floating a) =>
+   Upper sh a -> Array (Shape.Triangular Shape.Upper sh) a
+toUpperRowMajor =
+   Array.mapShape (Shape.Triangular Shape.Upper . MatrixShape.triangularSize)
+   .
+   forceOrder MatrixShape.RowMajor
+
+{-
+This is not maximally efficient.
+It fills up a whole square.
+This wastes memory but enables more regular memory access patterns.
+Additionally, it fills unused parts of the square with zero or mirrored values.
+-}
+forceOrder ::
+   (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag,
+    Shape.C sh, Class.Floating a) =>
+   Order -> Triangular lo diag up sh a -> Triangular lo diag up sh a
+forceOrder newOrder =
+   Tri.getMap $
+   MatrixShape.switchDiagUpLoSym
+      (Tri.Map $
+       Array.mapShape (\sh -> sh{MatrixShape.triangularOrder = newOrder}))
+      (forceOrderMap newOrder takeUpper)
+      (forceOrderMap newOrder takeLower)
+      (forceOrderMap newOrder $
+         Array.mapShape
+            (\sh -> sh{MatrixShape.triangularUplo = MatrixShape.autoUplo})
+         .
+         takeUpper)
+
+forceOrderMap ::
+   (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag,
+    Shape.C sh, Class.Floating a) =>
+   Order ->
+   (Square sh a -> Triangular lo NonUnit up sh a) ->
+   Tri.Map diag sh a lo up
+forceOrderMap newOrder f = Tri.Map $ \a ->
+   if MatrixShape.triangularOrder (Array.shape a) == newOrder
+      then a
+      else uncheckedRelaxNonUnitDiagonal $
+           f $ Basic.forceOrder newOrder $ toSquare a
+
+uncheckedRelaxNonUnitDiagonal ::
+   (MatrixShape.TriDiag diag) =>
+   Triangular lo NonUnit up sh a -> Triangular lo diag up sh a
+uncheckedRelaxNonUnitDiagonal =
+   Array.mapShape (\sh -> sh{MatrixShape.triangularDiag = MatrixShape.autoDiag})
+
+{- |
+@adaptOrder x y@ contains the data of @y@ with the layout of @x@.
+-}
+adaptOrder ::
+   (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag,
+    Shape.C sh, Class.Floating a) =>
+   Triangular lo diag up sh a ->
+   Triangular lo diag up sh a ->
+   Triangular lo diag up sh a
+adaptOrder x = forceOrder (MatrixShape.triangularOrder $ Array.shape x)
+
+add, sub ::
+   (MatrixShape.Content lo, MatrixShape.Content up,
+    Eq lo, Eq up, Eq sh, Shape.C sh, Class.Floating a) =>
+   Triangular lo NonUnit up sh a ->
+   Triangular lo NonUnit up sh a ->
+   Triangular lo NonUnit up sh a
+add x y = Vector.add (adaptOrder y x) y
+sub x y = Vector.sub (adaptOrder y x) y
 
 
 identity ::
diff --git a/src/Numeric/LAPACK/Matrix/Triangular/Eigen.hs b/src/Numeric/LAPACK/Matrix/Triangular/Eigen.hs
--- a/src/Numeric/LAPACK/Matrix/Triangular/Eigen.hs
+++ b/src/Numeric/LAPACK/Matrix/Triangular/Eigen.hs
@@ -12,8 +12,7 @@
           forPointers, rowMajorPointers)
 import Numeric.LAPACK.Matrix.Triangular.Basic (Triangular)
 import Numeric.LAPACK.Matrix.Shape.Private
-         (Order(ColumnMajor,RowMajor), caseLoUp, uploOrder,
-          NonUnit(NonUnit), triangleSize)
+         (Order(ColumnMajor,RowMajor), caseLoUp, uploOrder, NonUnit(NonUnit))
 import Numeric.LAPACK.Vector (Vector)
 import Numeric.LAPACK.Scalar (zero)
 import Numeric.LAPACK.Private (lacgv, withInfo, errorCodeMsg)
@@ -24,10 +23,11 @@
 import qualified Numeric.Netlib.Utility as Call
 import qualified Numeric.Netlib.Class as Class
 
-import qualified Data.Array.Comfort.Storable.Internal.Monadic as ArrayIO
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked.Monadic as ArrayIO
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
+import Data.Array.Comfort.Shape (triangleSize)
 
 import Foreign.C.Types (CInt, CChar)
 import Foreign.Ptr (Ptr, nullPtr)
diff --git a/src/Numeric/LAPACK/Matrix/Triangular/Linear.hs b/src/Numeric/LAPACK/Matrix/Triangular/Linear.hs
--- a/src/Numeric/LAPACK/Matrix/Triangular/Linear.hs
+++ b/src/Numeric/LAPACK/Matrix/Triangular/Linear.hs
@@ -19,8 +19,7 @@
 import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
 import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
 import Numeric.LAPACK.Matrix.Shape.Private
-         (transposeFromOrder, uploFromOrder, uploOrder,
-          charFromTriDiag, triangleSize)
+         (transposeFromOrder, uploFromOrder, uploOrder, charFromTriDiag)
 import Numeric.LAPACK.Matrix.Private (Conjugation(NonConjugated))
 import Numeric.LAPACK.Private (copyBlock, copyToTemp)
 
@@ -28,9 +27,10 @@
 import qualified Numeric.Netlib.Utility as Call
 import qualified Numeric.Netlib.Class as Class
 
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
+import Data.Array.Comfort.Shape (triangleSize)
 
 import System.IO.Unsafe (unsafePerformIO)
 
diff --git a/src/Numeric/LAPACK/Matrix/Triangular/Private.hs b/src/Numeric/LAPACK/Matrix/Triangular/Private.hs
--- a/src/Numeric/LAPACK/Matrix/Triangular/Private.hs
+++ b/src/Numeric/LAPACK/Matrix/Triangular/Private.hs
@@ -16,9 +16,9 @@
 import qualified Numeric.Netlib.Utility as Call
 import qualified Numeric.Netlib.Class as Class
 
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 
 import Foreign.Marshal.Alloc (alloca)
 import Foreign.Marshal.Array (advancePtr)
diff --git a/src/Numeric/LAPACK/Orthogonal.hs b/src/Numeric/LAPACK/Orthogonal.hs
--- a/src/Numeric/LAPACK/Orthogonal.hs
+++ b/src/Numeric/LAPACK/Orthogonal.hs
@@ -37,9 +37,9 @@
 import qualified Numeric.Netlib.Utility as Call
 import qualified Numeric.Netlib.Class as Class
 
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 
 import System.IO.Unsafe (unsafePerformIO)
 
diff --git a/src/Numeric/LAPACK/Orthogonal/Private.hs b/src/Numeric/LAPACK/Orthogonal/Private.hs
--- a/src/Numeric/LAPACK/Orthogonal/Private.hs
+++ b/src/Numeric/LAPACK/Orthogonal/Private.hs
@@ -26,10 +26,10 @@
 import qualified Numeric.Netlib.Utility as Call
 import qualified Numeric.Netlib.Class as Class
 
-import qualified Data.Array.Comfort.Storable.Internal.Monadic as ArrayIO
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked.Monadic as ArrayIO
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 
 import Foreign.Marshal.Array (advancePtr)
 import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)
diff --git a/src/Numeric/LAPACK/Permutation/Private.hs b/src/Numeric/LAPACK/Permutation/Private.hs
--- a/src/Numeric/LAPACK/Permutation/Private.hs
+++ b/src/Numeric/LAPACK/Permutation/Private.hs
@@ -15,9 +15,9 @@
 import qualified Numeric.Netlib.Utility as Call
 import qualified Numeric.Netlib.Class as Class
 
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 
 import qualified Text.PrettyPrint.Boxes as TextBox
 
diff --git a/src/Numeric/LAPACK/Private.hs b/src/Numeric/LAPACK/Private.hs
--- a/src/Numeric/LAPACK/Private.hs
+++ b/src/Numeric/LAPACK/Private.hs
@@ -29,7 +29,7 @@
 import Control.Monad (when, foldM)
 import Control.Applicative ((<$>))
 
-import qualified Data.Array.Comfort.Storable.Internal.Monadic as ArrayIO
+import qualified Data.Array.Comfort.Storable.Unchecked.Monadic as ArrayIO
 import qualified Data.Array.Comfort.Shape as Shape
 import Data.Array.Comfort.Storable (Array)
 
diff --git a/src/Numeric/LAPACK/Singular.hs b/src/Numeric/LAPACK/Singular.hs
--- a/src/Numeric/LAPACK/Singular.hs
+++ b/src/Numeric/LAPACK/Singular.hs
@@ -36,10 +36,10 @@
 import qualified Numeric.Netlib.Utility as Call
 import qualified Numeric.Netlib.Class as Class
 
-import qualified Data.Array.Comfort.Storable.Internal.Monadic as ArrayIO
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked.Monadic as ArrayIO
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 
 import System.IO.Unsafe (unsafePerformIO)
 
diff --git a/src/Numeric/LAPACK/Split.hs b/src/Numeric/LAPACK/Split.hs
--- a/src/Numeric/LAPACK/Split.hs
+++ b/src/Numeric/LAPACK/Split.hs
@@ -24,9 +24,9 @@
 import qualified Numeric.Netlib.Utility as Call
 import qualified Numeric.Netlib.Class as Class
 
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 
 import System.IO.Unsafe (unsafePerformIO)
 
diff --git a/src/Numeric/LAPACK/Vector.hs b/src/Numeric/LAPACK/Vector.hs
--- a/src/Numeric/LAPACK/Vector.hs
+++ b/src/Numeric/LAPACK/Vector.hs
@@ -65,9 +65,9 @@
 import Control.Monad.IO.Class (liftIO)
 import Control.Applicative (Const(Const,getConst), liftA3, (<$>))
 
-import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 import Data.Array.Comfort.Shape ((:+:)((:+:)))
 
 import Data.Complex (Complex)
diff --git a/test/Test/Banded.hs b/test/Test/Banded.hs
--- a/test/Test/Banded.hs
+++ b/test/Test/Banded.hs
@@ -221,8 +221,8 @@
 
 
 checkForAll ::
-   (Show a, QC.Testable test, Gen.Required required) =>
-   Gen.T tag required actual a -> (a -> test) -> Tagged tag QC.Property
+   (Show a, QC.Testable test) =>
+   Gen.T tag dim a -> (a -> test) -> Tagged tag QC.Property
 checkForAll gen = Util.checkForAll (Gen.run gen 10 5)
 
 
diff --git a/test/Test/BandedHermitian.hs b/test/Test/BandedHermitian.hs
--- a/test/Test/BandedHermitian.hs
+++ b/test/Test/BandedHermitian.hs
@@ -343,13 +343,13 @@
 
 
 checkForAll ::
-   (Show a, QC.Testable test, Gen.Required required) =>
-   Gen.T tag required actual a -> (a -> test) -> Tagged tag QC.Property
+   (Show a, QC.Testable test) =>
+   Gen.T tag dim a -> (a -> test) -> Tagged tag QC.Property
 checkForAll gen = Util.checkForAll (Gen.run gen 6 5)
 
 checkForAllExtra ::
-   (Show a, Show b, QC.Testable test, Gen.Required required) =>
-   QC.Gen a -> Gen.T tag required actual b ->
+   (Show a, Show b, QC.Testable test) =>
+   QC.Gen a -> Gen.T tag dim b ->
    (a -> b -> test) -> Tagged tag QC.Property
 checkForAllExtra = Gen.withExtra checkForAll
 
diff --git a/test/Test/Generator.hs b/test/Test/Generator.hs
--- a/test/Test/Generator.hs
+++ b/test/Test/Generator.hs
@@ -1,10 +1,16 @@
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE Rank2Types #-}
+{-# LANGUAGE EmptyDataDecls #-}
 module Test.Generator where
 
+import qualified Test.Logic as Logic
 import qualified Test.Utility as Util
-import Test.Utility (Match(Match,Mismatch))
+import Test.Logic (Dim, MatchMode(DontForceMatch,ForceMatch))
+import Test.Utility (Match)
 
+import qualified UniqueLogic.ST.TF.System.Simple as Sys
+
 import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
 import qualified Numeric.LAPACK.Matrix.Triangular as Triangular
 import qualified Numeric.LAPACK.Matrix.Square as Square
@@ -19,127 +25,125 @@
 import qualified Data.Array.Comfort.Storable as Array
 import qualified Data.Array.Comfort.Shape as Shape
 import Data.Array.Comfort.Storable (Array)
-import Data.Array.Comfort.Shape ((:+:)((:+:)))
+import Data.Array.Comfort.Shape ((:+:))
 
 import qualified Control.Monad.Trans.RWS as MRWS
-import qualified Control.Monad.Trans.Class as MT
-import qualified Control.Functor.HT as FuncHT
-import Control.Applicative (liftA2, (<$>))
+import qualified Control.Applicative.HT as AppHT
+import Control.Applicative (liftA2, (<*>), (<$>))
 
 import Data.Traversable (for)
-import Data.Tuple.HT (mapFst, mapSnd, mapPair, swap)
+import Data.Tuple.HT (swap, mapFst3, mapThd3)
 
 import qualified Test.QuickCheck as QC
 
 
 
 {- |
-@Cons generator@ with @generator maxElem maxDim fixedDims@.
-@generator@ constructs an array with possibly fixed height or width
-and returns its actual dimensions.
-Non-fixed dimensions will be choosen arbitrarily from the range @(0,maxDim)@.
+@Cons generator@ with @generator maxElem@.
+@generator@ constructs an array and maintains relations between the dimensions.
+Dimensions will be choosen arbitrarily from the range @(0,maxDim)@.
 Elements are choosen from the range @(-maxElem,maxElem)@.
--}
-newtype T tag required actual array = Cons (required -> ExtGen (array,actual))
 
-instance Functor (T tag required actual) where
-   fmap f (Cons gen) = Cons $ \fixed -> mapFst f <$> gen fixed
+I moved to the 's' tag to within the 'Cons' constructor
+and furthermore defined 'TaggedVariables' to strip the 's' tag
+from the Variables in 'dim'.
+This way, we can easily define 'checkForAll' in the test modules.
+Otherwise there would not be a way to quantify 'dim' while containing 's' tags.
+That is, we would have to reset 'dim' to () before every call to 'checkForAll'.
+-}
+newtype T tag dim array =
+   Cons {
+      decons :: forall s.
+         Integer ->
+         Logic.M s (TaggedVariables s dim, Logic.System s, Logic.M s array)
+   }
 
-type ExtGen = MRWS.RWST (Integer,Int,MatchMode) Match () QC.Gen
+data Variable dim
 
-data MatchMode = DontForceMatch | ForceMatch
-   deriving (Eq, Show)
+type family TaggedVariables s tuple
+type instance TaggedVariables s (Variable dim) = Logic.Variable s dim
+type instance TaggedVariables s () = ()
+type instance TaggedVariables s (a,b) =
+                  (TaggedVariables s a, TaggedVariables s b)
 
-class Required required where nothingRequired :: required
-instance Required () where nothingRequired = ()
-instance Required (Maybe a) where nothingRequired = Nothing
-instance (Required a, Required b) => Required (a,b) where
-   nothingRequired = (nothingRequired,nothingRequired)
+instance Functor (T tag dim) where
+   fmap f (Cons gen) = Cons $ \mapElem -> mapThd3 (fmap f) <$> gen mapElem
 
 run ::
-   (Required required) =>
-   T tag required actual array -> Integer -> Int ->
+   T tag dim array -> Integer -> Int ->
    Util.TaggedGen tag (array, Match)
-run (Cons gen) maxElem maxDim =
-   Util.Tagged $ do
-      forceMatch <- QC.elements [DontForceMatch, ForceMatch]
-      ((array, _actualDim), match) <-
-         MRWS.evalRWST (gen nothingRequired) (maxElem, maxDim, forceMatch) ()
-      return (array, match)
+run gen maxElem maxDim =
+   Util.Tagged $
+      QC.elements [DontForceMatch, ForceMatch] >>=
+      Logic.runSTInGen
+         (do (_dim, sys, queries) <- decons gen maxElem
+             Sys.solve sys
+             queries)
+         maxDim
 
 withExtra ::
-   (T tag required actual (a,b) -> ((a,b) -> c) -> io) ->
-   QC.Gen a -> T tag required actual b -> (a -> b -> c) -> io
+   (T tag dim (a,b) -> ((a,b) -> c) -> io) ->
+   QC.Gen a -> T tag dim b -> (a -> b -> c) -> io
 withExtra checkForAll genA genB test =
    checkForAll (mapGen (\_ b -> flip (,) b <$> genA) genB) (uncurry test)
 
 
 mapGen ::
    (Integer -> a -> QC.Gen b) ->
-   T tag required actual a -> T tag required actual b
+   T tag dim a -> T tag dim b
 mapGen f (Cons gen) =
-   Cons $ \fixed -> do
-      (maxElem, _maxDim, _match) <- MRWS.ask
-      MT.lift . FuncHT.mapFst (f maxElem) =<< gen fixed
+   Cons $ \maxElem ->
+      mapThd3 (Logic.liftGen . f maxElem =<<) <$> gen maxElem
 
 mapGenDim ::
    (Integer -> Int -> a -> QC.Gen b) ->
-   T tag required actual a -> T tag required actual b
+   T tag dim a -> T tag dim b
 mapGenDim f (Cons gen) =
-   Cons $ \fixed -> do
-      (maxElem, maxDim, _match) <- MRWS.ask
-      MT.lift . FuncHT.mapFst (f maxElem maxDim) =<< gen fixed
-
-
-chooseDimMin :: Int -> ExtGen Int
-chooseDimMin k = do
-   (_maxElem, maxDim, _match) <- MRWS.ask
-   MT.lift $ QC.choose (k,maxDim)
-
-
-class Dim dim where chooseDim :: ExtGen dim
-instance Dim Int where chooseDim = chooseDimMin 0
-instance (Dim dimA, Dim dimB) => Dim (dimA:+:dimB) where
-   chooseDim = liftA2 (:+:) chooseDim chooseDim
+   Cons $ \maxElem -> do
+      (maxDim, _matchMode) <- Logic.M MRWS.ask
+      mapThd3 (Logic.liftGen . f maxElem maxDim =<<) <$> gen maxElem
 
 
-matchDim :: (Dim i, Eq i) => i -> ExtGen i
-matchDim size = do
-   (_maxElem, _maxDim, match) <- MRWS.ask
-   case match of
-      ForceMatch -> return size
-      DontForceMatch -> do
-         newSize <- chooseDim
-         MRWS.tell $ if newSize==size then Match else Mismatch
-         return newSize
+combine ::
+   (forall s.
+    TaggedVariables s dimF -> TaggedVariables s dimA ->
+    (TaggedVariables s dimB, Logic.System s)) ->
+   T tag dimF (a -> b) ->
+   T tag dimA a ->
+   T tag dimB b
+combine combineDim (Cons genF) (Cons genA) =
+   Cons $ \maxElem -> do
+      (dimF,sysF,f) <- genF maxElem
+      (dimA,sysA,a) <- genA maxElem
+      let (dimB, constraint) = combineDim dimF dimA
+      return (dimB, sysF >> sysA >> constraint, f <*> a)
 
 
-type Scalar tag = T tag () ()
+type Scalar tag = T tag ()
 
 scalar :: (Class.Floating a) => Scalar a a
 scalar =
-   Cons $ \ _fixed -> do
-      (maxElem, _maxDim, _match) <- MRWS.ask
-      MT.lift $ flip (,) () <$> Util.genElement maxElem
+   Cons $ \ maxElem ->
+      return ((), return (), Logic.liftGen $ Util.genElement maxElem)
 
 (<.*.>) ::
+   (Dim size, Eq size) =>
    Vector tag size (a -> b) ->
    Vector tag size a ->
    Scalar tag b
-(<.*.>) (Cons genA) (Cons genB) =
-   Cons $ \() -> do
-      (f,size) <- genA Nothing
-      (a,_) <- genB $ Just size
-      return (f a, ())
+(<.*.>) = combine (\dimF dimA -> ((), Logic.ruleEqualDim dimF dimA))
 
 
-type Vector tag size = T tag (Maybe size) size
+queryZeroInt :: Logic.Variable s Int -> Logic.M s ZeroInt
+queryZeroInt var = zeroInt <$> Logic.query var
 
+type Vector tag size = T tag (Variable size)
+
 vectorDim :: (Class.Floating a) => Vector a Int ZeroInt
 vectorDim =
-   Cons $ \ fixed -> do
-      dims <- maybe chooseDim return fixed
-      return (zeroInt dims, dims)
+   Cons $ \ _maxElem -> do
+      dim <- Sys.globalVariable
+      return (dim, return (), queryZeroInt dim)
 
 vector :: (Class.Floating a) => Vector a Int (Vector.Vector ZeroInt a)
 vector = mapGen Util.genArray vectorDim
@@ -154,48 +158,34 @@
    Vector tag height (a -> b) ->
    Matrix tag height width a ->
    Vector tag width b
-(<.*|>) (Cons genA) (Cons genB) =
-   Cons $ \fixed -> do
-      (a,(height,width)) <- genB $ Right <$> fixed
-      (f,_) <- genA . Just =<< matchDim height
-      return (f a, width)
+(<.*|>) =
+   combine (\dim (height,width) -> (width, Logic.ruleEqualDim dim height))
 
 (<|*.>) ::
    (Dim width, Eq width) =>
    Matrix tag height width (a -> b) ->
    Vector tag width a ->
    Vector tag height b
-(<|*.>) (Cons genA) (Cons genB) =
-   Cons $ \fixed -> do
-      (f,(height,width)) <- genA $ Left <$> fixed
-      (a,_) <- genB . Just =<< matchDim width
-      return (f a, height)
+(<|*.>) =
+   combine (\(height,width) dim -> (height, Logic.ruleEqualDim width dim))
 
 (<.=.>) ::
    (Dim size, Eq size) =>
    Vector tag size (a -> b) ->
    Vector tag size a ->
    Vector tag size b
-(<.=.>) (Cons genA) (Cons genB) =
-   Cons $ \fixed -> do
-      (f,size) <- genA fixed
-      (a,_) <- genB . Just =<< matchDim size
-      return (f a, size)
+(<.=.>) =
+   combine (\sizeF sizeA -> (sizeF, Logic.ruleEqualDim sizeF sizeA))
 
 
-type Matrix tag height width =
-      T tag (Maybe (Either height width)) (height,width)
+type Matrix tag height width = T tag (Variable height, Variable width)
 
 matrixDims ::
    (Class.Floating a) => Matrix a Int Int (ZeroInt, ZeroInt)
 matrixDims =
-   Cons $ \ fixed -> do
-      dims <-
-         case fixed of
-            Nothing -> liftA2 (,) chooseDim chooseDim
-            Just (Left h) -> (,) h <$> chooseDim
-            Just (Right w) -> flip (,) w <$> chooseDim
-      return (mapPair (zeroInt,zeroInt) dims, dims)
+   Cons $ \ _maxElem -> do
+      dims <- liftA2 (,) Sys.globalVariable Sys.globalVariable
+      return (dims, return (), AppHT.mapPair (queryZeroInt,queryZeroInt) dims)
 
 matrix ::
    (Class.Floating a) => Matrix a Int Int (Matrix.General ZeroInt ZeroInt a)
@@ -207,13 +197,9 @@
 
 squareDim :: (Class.Floating a) => Matrix a Int Int ZeroInt
 squareDim =
-   Cons $ \ fixed -> do
-      size <-
-         case fixed of
-            Nothing -> chooseDim
-            Just (Left h) -> return h
-            Just (Right w) -> return w
-      return (zeroInt size, (size,size))
+   Cons $ \ _maxElem -> do
+      dim <- Sys.globalVariable
+      return ((dim,dim), return (), queryZeroInt dim)
 
 squareShaped ::
    (Shape.C sh, Class.Floating a) =>
@@ -302,20 +288,11 @@
 
 tallDims :: (Class.Floating a) => Matrix a Int Int (ZeroInt, ZeroInt)
 tallDims =
-   Cons $ \ fixed -> do
-      dims <-
-         case fixed of
-            Nothing -> do
-               h <- chooseDim
-               w <- MT.lift $ QC.choose (0,h)
-               return (h,w)
-            Just (Left h) -> do
-               w <- MT.lift $ QC.choose (0,h)
-               return (h,w)
-            Just (Right w) -> do
-               h <- chooseDimMin w
-               return (h,w)
-      return (mapPair (zeroInt,zeroInt) dims, dims)
+   Cons $ \ _maxElem -> do
+      height <- Sys.globalVariable
+      width  <- Sys.globalVariable
+      return ((height,width), Logic.ruleLessOrEqual width height,
+              liftA2 (,) (queryZeroInt height) (queryZeroInt width))
 
 tall ::
    (Class.Floating a) =>
@@ -335,40 +312,15 @@
          `QC.suchThat` Util.fullRankTall
 
 
-wideDims :: (Class.Floating a) => Matrix a Int Int (ZeroInt, ZeroInt)
-wideDims =
-   Cons $ \ fixed -> do
-      dims <-
-         case fixed of
-            Nothing -> do
-               w <- chooseDim
-               h <- MT.lift $ QC.choose (0,w)
-               return (h,w)
-            Just (Left h) -> do
-               w <- chooseDimMin h
-               return (h,w)
-            Just (Right w) -> do
-               h <- MT.lift $ QC.choose (0,w)
-               return (h,w)
-      return (mapPair (zeroInt,zeroInt) dims, dims)
-
 wide ::
    (Class.Floating a) =>
    Matrix a Int Int (Matrix.Wide ZeroInt ZeroInt a)
-wide =
-   flip mapGen wideDims $ \maxElem dims -> do
-      order <- Util.genOrder
-      Util.genArray maxElem $ uncurry (MatrixShape.wide order) dims
+wide = Matrix.transpose <$> transpose tall
 
 fullRankWide ::
    (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
    Matrix a Int Int (Matrix.Wide ZeroInt ZeroInt a)
-fullRankWide =
-   flip mapGen wideDims $ \maxElem dims -> do
-      order <- Util.genOrder
-      fmap Matrix.transpose $
-         Util.genArray maxElem (uncurry (MatrixShape.tall order) (swap dims))
-            `QC.suchThat` Util.fullRankTall
+fullRankWide = Matrix.transpose <$> transpose fullRankTall
 
 
 hermitian ::
@@ -400,27 +352,15 @@
    Matrix tag height fuse (a -> b) ->
    Matrix tag fuse width a ->
    Matrix tag height width b
-(<|*|>) (Cons genA) (Cons genB) =
-   Cons $ \fixed ->
-      case fixed of
-         Just (Right width) -> do
-            (a,(fuse,_)) <- genB $ Just $ Right width
-            (f,(height,_)) <- genA . Just . Right =<< matchDim fuse
-            return (f a, (height,width))
-         Just (Left height) -> do
-            (f,(_,fuse)) <- genA $ Just $ Left height
-            (a,(_,width)) <- genB . Just . Left =<< matchDim fuse
-            return (f a, (height,width))
-         Nothing -> do
-            (f,(height,fuse)) <- genA Nothing
-            (a,(_,width)) <- genB . Just . Left =<< matchDim fuse
-            return (f a, (height,width))
+(<|*|>) =
+   combine
+      (\(height,fuseF) (fuseA,width) ->
+         ((height,width), Logic.ruleEqualDim fuseF fuseA))
 
 transpose ::
    Matrix tag height width a ->
    Matrix tag width height a
-transpose (Cons gen) =
-   Cons $ fmap (mapSnd swap) . gen . fmap (either Right Left)
+transpose (Cons gen) = Cons $ \maxElem -> mapFst3 swap <$> gen maxElem
 
 (<|\|>) ::
    (Dim height, Eq height) =>
@@ -433,65 +373,43 @@
    Vector tag height (a -> b) ->
    Vector tag width a ->
    Matrix tag height width b
-(<***>) (Cons genA) (Cons genB) =
-   Cons $ \fixed -> do
-      (f,height) <- genA $ either Just (const Nothing) =<< fixed
-      (a,width) <- genB $ either (const Nothing) Just =<< fixed
-      return (f a, (height,width))
+(<***>) = combine (\height width -> ((height,width), return ()))
 
 
-{-
-We need this type because the test stackRowsColumnsCommutative
-requires to fix both height and width of the bottom right matrix.
-
-Conversely, we cannot use the type e.g. for Square matrices,
-because Square does not allow independent choice of height and width.
--}
-type Matrix2 tag height width =
-      T tag (Maybe height, Maybe width) (height,width)
-
-matrix2Dims :: (Class.Floating a) => Matrix2 a Int Int (ZeroInt, ZeroInt)
-matrix2Dims =
-   Cons $ \ (fixedHeight,fixedWidth) -> do
-      let maybeChooseDim = maybe chooseDim return
-      dims <-
-         liftA2 (,) (maybeChooseDim fixedHeight) (maybeChooseDim fixedWidth)
-      return (mapPair (zeroInt,zeroInt) dims, dims)
-
-matrix2 ::
-   (Class.Floating a) => Matrix2 a Int Int (Matrix.General ZeroInt ZeroInt a)
-matrix2 =
-   flip mapGen matrix2Dims $ \maxElem dims -> do
-      order <- Util.genOrder
-      Util.genArray maxElem $ uncurry (MatrixShape.general order) dims
+(<|=|>) ::
+   (Dim height, Eq height) =>
+   (Dim width, Eq width) =>
+   Matrix tag height width (a -> b) ->
+   Matrix tag height width a ->
+   Matrix tag height width b
+(<|=|>) =
+   combine $ \(heightF,widthF) (heightA,widthA) ->
+      ((heightF,widthF),
+       Logic.ruleEqualDim heightF heightA >> Logic.ruleEqualDim widthF widthA)
 
 (<===>) ::
    (Dim width, Eq width) =>
-   Matrix2 tag heightA width (a -> b) ->
-   Matrix2 tag heightB width a ->
-   Matrix2 tag (heightA:+:heightB) width b
-(<===>) (Cons genA) (Cons genB) =
-   Cons $ \(fixedHeight,fixedWidth) -> do
-      (f,(heightA,width)) <-
-         genA ((\(heightA:+:_) -> heightA) <$> fixedHeight, fixedWidth)
-      matchingWidth <- matchDim width
-      (a,(heightB,_)) <-
-         genB ((\(_:+:heightB) -> heightB) <$> fixedHeight, Just matchingWidth)
-      return (f a, (heightA:+:heightB, width))
+   Matrix tag heightA width (a -> b) ->
+   Matrix tag heightB width a ->
+   Matrix tag (heightA:+:heightB) width b
+(<===>) (Cons genF) (Cons genA) =
+   Cons $ \maxElem -> do
+      heightB <- Sys.globalVariable
+      ((heightF,widthF),sysF,f) <- genF maxElem
+      ((heightA,widthA),sysA,a) <- genA maxElem
+      return
+         ((heightB,widthF),
+          do sysF >> sysA
+             Logic.ruleEqualDim widthF widthA
+             Logic.ruleAppendDim heightF heightA heightB,
+          f <*> a)
 
 (<|||>) ::
    (Dim height, Eq height) =>
-   Matrix2 tag height widthA (a -> b) ->
-   Matrix2 tag height widthB a ->
-   Matrix2 tag height (widthA:+:widthB) b
-(<|||>) (Cons genA) (Cons genB) =
-   Cons $ \(fixedHeight,fixedWidth) -> do
-      (f,(height,widthA)) <-
-         genA (fixedHeight, (\(widthA:+:_) -> widthA) <$> fixedWidth)
-      matchingHeight <- matchDim height
-      (a,(_,widthB)) <-
-         genB (Just matchingHeight, (\(_:+:widthB) -> widthB) <$> fixedWidth)
-      return (f a, (height, widthA:+:widthB))
+   Matrix tag height widthA (a -> b) ->
+   Matrix tag height widthB a ->
+   Matrix tag height (widthA:+:widthB) b
+(<|||>) f a = transpose $ transpose f <===> transpose a
 
 
-infixl 4 <.*.>, <.*|>, <|*.>, <|*|>, <|\|>, <***>, <.=.>, <===>, <|||>
+infixl 4 <.*.>, <.*|>, <|*.>, <|*|>, <|\|>, <***>, <.=.>, <|=|>, <===>, <|||>
diff --git a/test/Test/Hermitian.hs b/test/Test/Hermitian.hs
--- a/test/Test/Hermitian.hs
+++ b/test/Test/Hermitian.hs
@@ -301,13 +301,13 @@
 
 
 checkForAll ::
-   (Show a, QC.Testable test, Gen.Required required) =>
-   Gen.T tag required actual a -> (a -> test) -> Tagged tag QC.Property
+   (Show a, QC.Testable test) =>
+   Gen.T tag dim a -> (a -> test) -> Tagged tag QC.Property
 checkForAll gen = Util.checkForAll (Gen.run gen 3 5)
 
 checkForAllExtra ::
-   (Show a, Show b, QC.Testable test, Gen.Required required) =>
-   QC.Gen a -> Gen.T tag required actual b ->
+   (Show a, Show b, QC.Testable test) =>
+   QC.Gen a -> Gen.T tag dim b ->
    (a -> b -> test) -> Tagged tag QC.Property
 checkForAllExtra = Gen.withExtra checkForAll
 
diff --git a/test/Test/Logic.hs b/test/Test/Logic.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Logic.hs
@@ -0,0 +1,125 @@
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE Rank2Types #-}
+module Test.Logic where
+
+import Test.Utility (Match(Match,Mismatch))
+
+import qualified UniqueLogic.ST.TF.Rule as Rule
+import qualified UniqueLogic.ST.TF.System.Simple as Sys
+import qualified UniqueLogic.ST.TF.System as Sys (runApplyM)
+
+import qualified Data.Ref as Ref
+import Data.STRef (newSTRef, writeSTRef, readSTRef)
+
+import Data.Array.Comfort.Shape ((:+:)((:+:)))
+
+import qualified Control.Monad.Trans.Class as MT
+import qualified Control.Monad.Trans.RWS as MRWS
+import Control.Monad.ST (ST, runST)
+import Control.Applicative (Applicative, liftA2, (<$>))
+
+import qualified Test.QuickCheck as QC
+import qualified Test.QuickCheck.GenT as GenT
+import Test.QuickCheck.Gen (Gen(MkGen))
+import Test.QuickCheck.GenT (GenT)
+
+import System.Random (Random)
+
+
+data MatchMode = DontForceMatch | ForceMatch
+   deriving (Eq, Show)
+
+newtype M s a = M {runM :: MRWS.RWST (Int, MatchMode) Match () (GenT (ST s)) a}
+   deriving (Functor, Applicative, Monad)
+
+instance Ref.C (M s) where
+   new =
+      Ref.newCons
+         (liftST . newSTRef)
+         ((liftST .) . writeSTRef)
+         (liftST . readSTRef)
+
+liftST :: ST s a -> M s a
+liftST = M . MT.lift . MT.lift
+
+liftGen :: QC.Gen a -> M s a
+liftGen = M . MT.lift . GenT.liftGen
+
+
+type Variable s a = Sys.Variable (M s) a
+type System s = Sys.T (M s) ()
+
+
+example :: Int -> MatchMode -> QC.Gen ([Int], Match)
+example =
+   runSTInGen (do
+      a <- Sys.globalVariable
+      b <- Sys.globalVariable
+      c <- Sys.globalVariable
+      d <- Sys.globalVariable
+      e <- Sys.globalVariable
+      f <- Sys.globalVariable
+      Sys.solve $ do
+         ruleLessOrEqual a b
+         ruleEqualDim d f
+         ruleEqualDim c d
+      mapM query [a,b,c,d,e,f])
+
+
+choose :: (Random a) => (Int -> (a,a)) -> M s a
+choose f = liftGen . QC.choose . f . fst =<< M MRWS.ask
+
+ruleLessOrEqual :: Variable s Int -> Variable s Int -> System s
+ruleLessOrEqual va vb = do
+   assignmentM (\a -> choose (\ maxk -> (a,maxk))) va vb
+   assignmentM (\b -> choose (\_maxk -> (0,b))) vb va
+
+
+class Dim dim where chooseDim :: M s dim
+instance Dim Int where chooseDim = choose ((,) 0)
+instance (Dim dimA, Dim dimB) => Dim (dimA:+:dimB) where
+   chooseDim = liftA2 (:+:) chooseDim chooseDim
+
+ruleEqualDim ::
+   (Dim dim, Eq dim) => Variable s dim -> Variable s dim -> System s
+ruleEqualDim va vb = do
+   let equalM =
+         assignmentM $ \x -> do
+            matchMode <- M $ MRWS.asks snd
+            case matchMode of
+               ForceMatch -> return x
+               DontForceMatch -> do
+                  y <- chooseDim
+                  M $ MRWS.tell $ if x==y then Match else Mismatch
+                  return y
+   equalM va vb
+   equalM vb va
+
+assignmentM ::
+   (Ref.C s) => (a -> s b) -> Sys.Variable s a -> Sys.Variable s b -> Sys.T s ()
+assignmentM f vx vy = Sys.runApplyM (f <$> Sys.arg vx) vy
+
+
+ruleAppendDim ::
+   Variable s dimA -> Variable s dimB -> Variable s (dimA:+:dimB) -> System s
+ruleAppendDim va vb vab = do
+   Sys.assignment3 (:+:) va vb vab
+   Sys.assignment2 (\(a:+:_) -> a) vab va
+   Sys.assignment2 (\(_:+:b) -> b) vab vb
+
+
+runSTInGen :: (forall s. M s b) -> Int -> MatchMode -> QC.Gen (b, Match)
+runSTInGen m =
+   \maxDim matchMode -> MkGen $ \r n ->
+      runST (GenT.unGenT (MRWS.evalRWST (runM m) (maxDim,matchMode) ()) r n)
+
+query :: (Dim dim) => Variable s dim -> M s dim
+query v = do
+   mk <- Sys.query v
+   case mk of
+      Just k -> return k
+      Nothing -> do
+         k <- chooseDim
+         Sys.solve $ Rule.equ v =<< Sys.constant k
+         return k
diff --git a/test/Test/Matrix.hs b/test/Test/Matrix.hs
--- a/test/Test/Matrix.hs
+++ b/test/Test/Matrix.hs
@@ -4,7 +4,8 @@
 
 import qualified Test.Generator as Gen
 import qualified Test.Utility as Util
-import Test.Generator ((<|*|>), (<|*.>), (<.*.>), (<***>), (<|||>), (<===>))
+import Test.Generator
+         ((<|*|>), (<|*.>), (<.*.>), (<***>), (<|=|>), (<|||>), (<===>))
 import Test.Utility
          (approx, approxArray, approxMatrix,
           genOrder, Tagged(Tagged), TaggedGen)
@@ -159,7 +160,8 @@
 
 outerInner ::
    (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
-   MatrixShape.Order -> (Vector ZeroInt a, Vector ZeroInt a, Vector ZeroInt a) -> Bool
+   MatrixShape.Order ->
+   (Vector ZeroInt a, Vector ZeroInt a, Vector ZeroInt a) -> Bool
 outerInner order (x,y,z) =
    approxArray (Matrix.outer order x y #> z) (Vector.scale (Vector.inner y z) x)
 
@@ -200,15 +202,15 @@
 
 genMatrix2EqHeight ::
    (Class.Floating a) =>
-   Gen.Matrix2 a Int (Int:+:Int)
+   Gen.Matrix a Int (Int:+:Int)
       (General ZeroInt ZeroInt a, General ZeroInt ZeroInt a)
-genMatrix2EqHeight = (,) <$> Gen.matrix2 <|||> Gen.matrix2
+genMatrix2EqHeight = (,) <$> Gen.matrix <|||> Gen.matrix
 
 genMatrix2EqWidth ::
    (Class.Floating a) =>
-   Gen.Matrix2 a (Int:+:Int) Int
+   Gen.Matrix a (Int:+:Int) Int
       (General ZeroInt ZeroInt a, General ZeroInt ZeroInt a)
-genMatrix2EqWidth = (,) <$> Gen.matrix2 <===> Gen.matrix2
+genMatrix2EqWidth = (,) <$> Gen.matrix <===> Gen.matrix
 
 reverseRows :: (Class.Floating a) => General ZeroInt ZeroInt a -> Bool
 reverseRows x =
@@ -375,6 +377,32 @@
       (Matrix.toRowMajor $ (x===z)|||(y===w))
 
 
+forceOrder ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   MatrixShape.Order ->
+   (General ZeroInt ZeroInt a, Vector ZeroInt a) ->
+   Bool
+forceOrder order (a,x) =
+   let ao = Matrix.forceOrder order a
+   in MatrixShape.fullOrder (Array.shape ao) == order
+      &&
+      approxArray (a #> x) (ao #> x)
+
+addDistributive ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   ((General ZeroInt ZeroInt a, General ZeroInt ZeroInt a), Vector ZeroInt a) ->
+   Bool
+addDistributive ((a,b),x) =
+   approxArray (Matrix.add a b #> x) (Vector.add (a#>x) (b#>x))
+
+subDistributive ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   ((General ZeroInt ZeroInt a, General ZeroInt ZeroInt a), Vector ZeroInt a) ->
+   Bool
+subDistributive ((a,b),x) =
+   approxArray (Matrix.sub a b #> x) (Vector.sub (a#>x) (b#>x))
+
+
 multiplyDiagonalMatrix ::
    (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
    (Triangular.Diagonal ZeroInt a, General ZeroInt ZeroInt a) -> Bool
@@ -390,13 +418,13 @@
 
 
 checkForAll ::
-   (Show a, QC.Testable test, Gen.Required required) =>
-   Gen.T tag required actual a -> (a -> test) -> Tagged tag QC.Property
+   (Show a, QC.Testable test) =>
+   Gen.T tag dim a -> (a -> test) -> Tagged tag QC.Property
 checkForAll gen = Util.checkForAll (Gen.run gen 10 5)
 
 checkForAllExtra ::
-   (Show a, Show b, QC.Testable test, Gen.Required required) =>
-   QC.Gen a -> Gen.T tag required actual b ->
+   (Show a, Show b, QC.Testable test) =>
+   QC.Gen a -> Gen.T tag dim b ->
    (a -> b -> test) -> Tagged tag QC.Property
 checkForAllExtra = Gen.withExtra checkForAll
 
@@ -483,16 +511,34 @@
       checkForAll genMatrix2EqHeight dropStackColumns) :
    ("stackRowsAssociative",
       checkForAll
-         ((,,) <$> Gen.matrix2 <===> Gen.matrix2 <===> Gen.matrix2)
+         ((,,) <$> Gen.matrix <===> Gen.matrix <===> Gen.matrix)
          stackRowsAssociative) :
    ("stackColumnsAssociative",
       checkForAll
-         ((,,) <$> Gen.matrix2 <|||> Gen.matrix2 <|||> Gen.matrix2)
+         ((,,) <$> Gen.matrix <|||> Gen.matrix <|||> Gen.matrix)
          stackColumnsAssociative) :
    ("stackRowsColumnsCommutative",
       checkForAll
          ((,) <$> genMatrix2EqHeight <===> genMatrix2EqHeight)
          stackRowsColumnsCommutative) :
+
+   ("forceOrder",
+      checkForAllExtra genOrder
+         ((,) <$> Gen.matrix <|*.> Gen.vector) forceOrder) :
+   ("addDistributive",
+      checkForAll
+         ((,) <$>
+            ((,) <$> Gen.matrix <|=|> Gen.matrix)
+            <|*.>
+            Gen.vector)
+         addDistributive) :
+   ("subDistributive",
+      checkForAll
+         ((,) <$>
+            ((,) <$> Gen.matrix <|=|> Gen.matrix)
+            <|*.>
+            Gen.vector)
+         subDistributive) :
 
    ("multiplyDiagonalMatrix",
       checkForAll
diff --git a/test/Test/Orthogonal.hs b/test/Test/Orthogonal.hs
--- a/test/Test/Orthogonal.hs
+++ b/test/Test/Orthogonal.hs
@@ -334,8 +334,8 @@
 
 
 checkForAll ::
-   (Show a, QC.Testable test, Gen.Required required) =>
-   Gen.T tag required actual a -> (a -> test) -> Tagged tag QC.Property
+   (Show a, QC.Testable test) =>
+   Gen.T tag dim a -> (a -> test) -> Tagged tag QC.Property
 checkForAll gen = Util.checkForAll (Gen.run gen 3 5)
 
 
diff --git a/test/Test/Singular.hs b/test/Test/Singular.hs
--- a/test/Test/Singular.hs
+++ b/test/Test/Singular.hs
@@ -116,8 +116,8 @@
 
 
 checkForAll ::
-   (Show a, QC.Testable test, Gen.Required required) =>
-   Gen.T tag required actual a -> (a -> test) -> Tagged tag QC.Property
+   (Show a, QC.Testable test) =>
+   Gen.T tag dim a -> (a -> test) -> Tagged tag QC.Property
 checkForAll gen = Util.checkForAll (Gen.run gen 3 5)
 
 testsVar ::
diff --git a/test/Test/Square.hs b/test/Test/Square.hs
--- a/test/Test/Square.hs
+++ b/test/Test/Square.hs
@@ -135,8 +135,8 @@
 
 
 checkForAll ::
-   (Show a, QC.Testable test, Gen.Required required) =>
-   Gen.T tag required actual a -> (a -> test) -> Tagged tag QC.Property
+   (Show a, QC.Testable test) =>
+   Gen.T tag dim a -> (a -> test) -> Tagged tag QC.Property
 checkForAll gen = Util.checkForAll (Gen.run gen 3 5)
 
 
diff --git a/test/Test/Triangular.hs b/test/Test/Triangular.hs
--- a/test/Test/Triangular.hs
+++ b/test/Test/Triangular.hs
@@ -5,7 +5,7 @@
 
 import qualified Test.Generator as Gen
 import qualified Test.Utility as Util
-import Test.Generator ((<.*|>), (<|*.>), (<|*|>), (<|\|>))
+import Test.Generator ((<.*|>), (<|*.>), (<|*|>), (<|\|>), (<|=|>))
 import Test.Utility (approx, approxArray, approxArrayTol, approxMatrix, Tagged)
 
 import qualified Numeric.LAPACK.Matrix.Triangular as Triangular
@@ -32,6 +32,46 @@
 import qualified Test.QuickCheck as QC
 
 
+forceOrder ::
+   (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag,
+    Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   MatrixShape.Order ->
+   (Triangular lo diag up ZeroInt a, Vector ZeroInt a) -> Bool
+forceOrder order (a,x) =
+   let ao = Triangular.forceOrder order a
+   in MatrixShape.triangularOrder (Array.shape ao) == order
+      &&
+      approxArray (a #> x) (ao #> x)
+
+addDistributive ::
+   (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag,
+    Eq lo, Eq up, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   ((Triangular lo diag up ZeroInt a,
+     Triangular lo diag up ZeroInt a),
+    Vector ZeroInt a) ->
+   Bool
+addDistributive ((a,b),x) =
+   approxArray
+      (Triangular.add
+         (Triangular.strictNonUnitDiagonal a)
+         (Triangular.strictNonUnitDiagonal b) #> x)
+      (Vector.add (a#>x) (b#>x))
+
+subDistributive ::
+   (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag,
+    Eq lo, Eq up, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   ((Triangular lo diag up ZeroInt a,
+     Triangular lo diag up ZeroInt a),
+    Vector ZeroInt a) ->
+   Bool
+subDistributive ((a,b),x) =
+   approxArray
+      (Triangular.sub
+         (Triangular.strictNonUnitDiagonal a)
+         (Triangular.strictNonUnitDiagonal b) #> x)
+      (Vector.sub (a#>x) (b#>x))
+
+
 multiplyIdentityVector ::
    (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag,
     Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
@@ -245,13 +285,13 @@
 
 
 checkForAll ::
-   (Show a, QC.Testable test, Gen.Required required) =>
-   Gen.T tag required actual a -> (a -> test) -> Tagged tag QC.Property
+   (Show a, QC.Testable test) =>
+   Gen.T tag dim a -> (a -> test) -> Tagged tag QC.Property
 checkForAll gen = Util.checkForAll (Gen.run gen 3 5)
 
 checkForAllExtra ::
-   (Show a, Show b, QC.Testable test, Gen.Required required) =>
-   QC.Gen a -> Gen.T tag required actual b ->
+   (Show a, Show b, QC.Testable test) =>
+   QC.Gen a -> Gen.T tag dim b ->
    (a -> b -> test) -> Tagged tag QC.Property
 checkForAllExtra = Gen.withExtra checkForAll
 
@@ -341,6 +381,24 @@
    (Show a, Show ar, Class.Floating a, Eq a, RealOf a ~ ar, Class.Real ar) =>
    [(String, Tagged a QC.Property)]
 testsVar =
+   checkAny "forceOrder"
+      (\gen ->
+         checkForAllExtra Util.genOrder
+            ((,) <$> gen <|*.> Gen.vector) forceOrder)
+      Gen.triangular ++
+   checkAny "addDistributive"
+      (\gen ->
+         checkForAll
+            ((,) <$> ((,) <$> gen <|=|> gen) <|*.> Gen.vector)
+            addDistributive)
+      Gen.triangular ++
+   checkAny "subDistributive"
+      (\gen ->
+         checkForAll
+            ((,) <$> ((,) <$> gen <|=|> gen) <|*.> Gen.vector)
+            subDistributive)
+      Gen.triangular ++
+
    checkAny "multiplyIdentityVector"
       (\gen -> checkForAll ((,) <$> gen <|*.> Gen.vector) multiplyIdentityVector)
       (Triangular.relaxUnitDiagonal <$> Gen.identity) ++
diff --git a/test/Test/Vector.hs b/test/Test/Vector.hs
--- a/test/Test/Vector.hs
+++ b/test/Test/Vector.hs
@@ -101,8 +101,8 @@
 
 
 checkForAll ::
-   (Show a, QC.Testable test, Gen.Required required) =>
-   Gen.T tag required actual a -> (a -> test) -> Tagged tag QC.Property
+   (Show a, QC.Testable test) =>
+   Gen.T tag dim a -> (a -> test) -> Tagged tag QC.Property
 checkForAll gen = Util.checkForAll (Gen.run gen 10 5)