diff --git a/Changes.md b/Changes.md
--- a/Changes.md
+++ b/Changes.md
@@ -1,5 +1,14 @@
 # Change log for the `lapack` package
 
+## 0.5
+
+ * `Matrix.Block` for Block matrices.
+   Add `*Extra` constraint families to many type classes
+   in order to handle the data stored in the extra type parameters of `Matrix`.
+
+ * `Format.format` now uses custom type `Config`
+   instead of a plain format `String`.
+
 ## 0.4
 
  * Unified `Matrix` type that provides the same type parameters
diff --git a/lapack.cabal b/lapack.cabal
--- a/lapack.cabal
+++ b/lapack.cabal
@@ -1,6 +1,6 @@
 Cabal-Version:    2.2
 Name:             lapack
-Version:          0.4.1
+Version:          0.5
 License:          BSD-3-Clause
 License-File:     LICENSE
 Author:           Henning Thielemann <haskell@henning-thielemann.de>
@@ -53,7 +53,7 @@
   test-module.list
 
 Source-Repository this
-  Tag:         0.4.1
+  Tag:         0.5
   Type:        darcs
   Location:    https://hub.darcs.net/thielema/lapack/
 
@@ -95,6 +95,9 @@
   Exposed-Modules:
     Numeric.LAPACK.Matrix
     Numeric.LAPACK.Matrix.Special
+    Numeric.LAPACK.Matrix.Block
+    Numeric.LAPACK.Matrix.Block.Type
+    Numeric.LAPACK.Matrix.Type
     Numeric.LAPACK.Matrix.Superscript
     Numeric.LAPACK.Matrix.Array
     Numeric.LAPACK.Matrix.Extent
@@ -134,7 +137,7 @@
     Numeric.LAPACK.Linear.Private
     Numeric.LAPACK.Split
     Numeric.LAPACK.Permutation.Private
-    Numeric.LAPACK.Matrix.Block
+    Numeric.LAPACK.Matrix.Block.Private
     Numeric.LAPACK.Matrix.Square.Basic
     Numeric.LAPACK.Matrix.Square.Linear
     Numeric.LAPACK.Matrix.Square.Eigen
@@ -170,7 +173,8 @@
     Numeric.LAPACK.Matrix.Plain.Format
     Numeric.LAPACK.Matrix.Plain
     Numeric.LAPACK.Matrix.Indexed
-    Numeric.LAPACK.Matrix.Type
+    Numeric.LAPACK.Matrix.Type.Private
+    Numeric.LAPACK.Matrix.Wrapper
     Numeric.LAPACK.Matrix.Inverse
     Numeric.LAPACK.Matrix.Array.Indexed
     Numeric.LAPACK.Matrix.Array.Banded
@@ -224,6 +228,7 @@
     DocTest.Main
     Test.Shape
     Test.Indexed
+    Test.Block
     Test.Function
     Test.Divide
     Test.Multiply
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
@@ -1,21 +1,25 @@
 {-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE ConstraintKinds #-}
+{-# LANGUAGE UndecidableInstances #-}
 module Numeric.LAPACK.Format (
    (##),
    print,
    hyper,
    Format(format),
    FormatArray(formatArray),
-   FormatMatrix(formatMatrix),
    ArrFormat.deflt,
+   defltConfig,
+   Config(..),
    ) where
 
+import qualified Numeric.LAPACK.Matrix.Type.Private as Matrix
 import qualified Numeric.LAPACK.Matrix.Plain.Format as ArrFormat
 import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
 import qualified Numeric.LAPACK.Output as Output
 import qualified Numeric.LAPACK.Permutation.Private as Perm
-import Numeric.LAPACK.Matrix.Type (Matrix, FormatMatrix(formatMatrix))
+import Numeric.LAPACK.Matrix.Type.Private (Matrix)
 import Numeric.LAPACK.Matrix.Plain.Format
-         (FormatArray, formatArray, printfComplex)
+         (Config(..), defltConfig, FormatArray, formatArray, printfFloating)
 import Numeric.LAPACK.Output (Output, (/+/))
 
 import qualified Numeric.Netlib.Class as Class
@@ -37,52 +41,53 @@
 
 infix 0 ##
 
-print :: (Format a) => String -> a -> IO ()
-print fmt a = putStr $ trim $ TextBox.render $ format fmt a
+print :: (Format a) => Config -> a -> IO ()
+print cfg a = putStr $ trim $ TextBox.render $ format cfg a
 
 (##) :: (Format a) => a -> String -> IO ()
-(##) = flip print
+a##cfg = print (defltConfig {configFormat = cfg}) a
 
 trim :: String -> String
 trim = unlines . map (ListRev.dropWhile isSpace) . lines
 
-hyper :: (Format a) => String -> a -> Hyper.Graphic
-hyper fmt = Output.hyper . format fmt
+hyper :: (Format a) => Config -> a -> Hyper.Graphic
+hyper cfg = Output.hyper . format cfg
 
 
 class Format a where
-   format :: (Output out) => String -> a -> out
+   format :: (Output out) => Config -> a -> out
 
 instance Format Int where
-   format _fmt = Output.text . show
+   format _cfg = Output.text . show
 
 instance Format Float where
-   format fmt = Output.text . printf fmt
+   format cfg = Output.text . printf (configFormat cfg)
 
 instance Format Double where
-   format fmt = Output.text . printf fmt
+   format cfg = Output.text . printf (configFormat cfg)
 
 instance (Class.Real a) => Format (Complex a) where
-   format fmt = Output.text . fold . printfComplex fmt
+   format cfg = Output.text . fold . printfFloating cfg
 
 instance (Format a) => Format [a] where
-   format fmt = Output.formatColumn . map (format fmt)
+   format cfg = Output.formatColumn . map (format cfg)
 
 instance (Format a, Format b) => Format (a,b) where
-   format fmt (a,b) = format fmt a /+/ format fmt b
+   format cfg (a,b) = format cfg a /+/ format cfg b
 
 instance (Format a, Format b, Format c) => Format (a,b,c) where
-   format fmt (a,b,c) = format fmt a /+/ format fmt b /+/ format fmt c
+   format cfg (a,b,c) = format cfg a /+/ format cfg b /+/ format cfg c
 
 instance (Shape.C sh) => Format (Perm.Permutation sh) where
-   format _fmt = Perm.format
+   format _cfg = Perm.format
 
 instance (FormatArray sh, Class.Floating a) => Format (Array sh a) where
    format = formatArray
 
 instance
-   (FormatMatrix typ, Extent.Measure meas, Extent.C vert, Extent.C horiz,
-    Shape.C width, Shape.C height, Class.Floating a) =>
+   (Matrix.Format typ, Matrix.FormatExtra typ xl, Matrix.FormatExtra typ xu,
+    Extent.Measure meas, Extent.C vert, Extent.C horiz,
+    Shape.C height, Shape.C width, Class.Floating a) =>
    Format
       (Matrix typ xl xu lower upper meas vert horiz height width a) where
-   format = formatMatrix
+   format = Matrix.format
diff --git a/src/Numeric/LAPACK/Linear/Plain.hs b/src/Numeric/LAPACK/Linear/Plain.hs
--- a/src/Numeric/LAPACK/Linear/Plain.hs
+++ b/src/Numeric/LAPACK/Linear/Plain.hs
@@ -3,6 +3,7 @@
 {-# LANGUAGE EmptyDataDecls #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE ConstraintKinds #-}
 {-# LANGUAGE StandaloneDeriving #-}
 module Numeric.LAPACK.Linear.Plain (
    LowerUpper,
@@ -36,10 +37,11 @@
 
 import qualified Numeric.LAPACK.Matrix.Divide as Divide
 import qualified Numeric.LAPACK.Matrix.Multiply as Multiply
-import qualified Numeric.LAPACK.Matrix.Type as Matrix
+import qualified Numeric.LAPACK.Matrix.Type.Private as Matrix
 import qualified Numeric.LAPACK.Matrix.Banded.Basic as Banded
 import qualified Numeric.LAPACK.Matrix.Basic as Basic
 import qualified Numeric.LAPACK.Matrix.Array.Private as ArrMatrix
+import qualified Numeric.LAPACK.Matrix.Plain.Format as ArrFormat
 import qualified Numeric.LAPACK.Matrix.Layout as LayoutPub
 import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
 import qualified Numeric.LAPACK.Matrix.Extent.Strict as ExtentStrict
@@ -50,7 +52,7 @@
 import qualified Numeric.LAPACK.Split as Split
 import Numeric.LAPACK.Output ((/+/))
 import Numeric.LAPACK.Matrix.Plain.Format (formatArray)
-import Numeric.LAPACK.Matrix.Type (Matrix, FormatMatrix(formatMatrix))
+import Numeric.LAPACK.Matrix.Type.Private (Matrix)
 import Numeric.LAPACK.Matrix.Triangular.Basic (Lower, Upper)
 import Numeric.LAPACK.Matrix.Layout.Private
          (Order(RowMajor, ColumnMajor), Triangle(Triangle))
@@ -121,13 +123,20 @@
 type SquareMeas meas height width =
          LowerUpper meas Extent.Small Extent.Small height width
 
-instance FormatMatrix LU where
-   formatMatrix fmt lu@(LowerUpper _ipiv m) =
+instance Matrix.Format LU where
+   type FormatExtra LU extra = extra ~ ()
+   format fmt lu@(LowerUpper _ipiv m) =
       Perm.format (extractP NonInverted lu)
       /+/
       formatArray fmt m
 
+instance Matrix.Layout LU where
+   type LayoutExtra LU extra = extra ~ ()
+   layout (LowerUpper _ipiv m) =
+      ArrFormat.splitArrayFromList2 (Layout.splitExtent $ Array.shape m) $
+      ArrFormat.layoutSplit m
 
+
 mapExtent ::
    (Extent.C vertA, Extent.C horizA) =>
    (Extent.C vertB, Extent.C horizB) =>
@@ -575,6 +584,7 @@
 
 
 instance Matrix.Box LU where
+   type BoxExtra LU extra = extra ~ ()
    extent = Layout.splitExtent . Array.shape . split_
 
 instance Matrix.ToQuadratic LU where
@@ -591,10 +601,13 @@
 layoutPivotSquare :: sh -> Layout.Diagonal sh
 layoutPivotSquare = LayoutPub.diagonal Layout.ColumnMajor
 
-instance (xl ~ (), xu ~ ()) => Matrix.MapExtent LU xl xu lower upper where
+instance Matrix.MapExtent LU where
+   type MapExtentExtra LU extra = extra ~ ()
+   type MapExtentStrip LU strip = ()
    mapExtent = mapExtent
 
-instance (xl ~ (), xu ~ ()) => Multiply.MultiplyVector LU xl xu where
+instance Multiply.MultiplyVector LU where
+   type MultiplyVectorExtra LU extra = extra ~ ()
    matrixVector lu =
       Basic.unliftColumn Layout.ColumnMajor
          (multiplyFull (mapExtent Extent.toGeneral lu))
@@ -605,7 +618,8 @@
             Array.mapShape deconsUnchecked .
             transMultiplyVector (mapWidth Unchecked lu)
 
-instance (xl ~ (), xu ~ ()) => Multiply.MultiplySquare LU xl xu where
+instance Multiply.MultiplySquare LU where
+   type MultiplySquareExtra LU extra = extra ~ ()
    squareFull lu =
       ArrMatrix.lift1 $
          multiplyP NonInverted lu .
@@ -620,8 +634,10 @@
          multiplyP Inverted lu .
          Basic.transpose
 
-instance (xl ~ (), xu ~ ()) => Divide.Determinant LU xl xu where
+instance Divide.Determinant LU where
+   type DeterminantExtra LU extra = extra ~ ()
    determinant = determinant
 
-instance (xl ~ (), xu ~ ()) => Divide.Solve LU xl xu where
+instance Divide.Solve LU where
+   type SolveExtra LU extra = extra ~ ()
    solve trans = ArrMatrix.lift1 . solveTrans trans
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
@@ -1,5 +1,6 @@
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE ConstraintKinds #-}
 module Numeric.LAPACK.Matrix (
    Matrix.Matrix,
    Full,
@@ -17,20 +18,26 @@
    BandedHermitian.BandedHermitian,
    PermMatrix.Permutation,
 
+   Matrix.Format(Matrix.format), Matrix.FormatExtra, Matrix.formatWithLayout,
+   Matrix.Layout(Matrix.layout), Matrix.LayoutExtra,
+
    ShapeInt, shapeInt,
-   Matrix.Transpose, Matrix.transpose, MatrixClass.adjoint,
-   Matrix.height, Matrix.width,
-   Matrix.Box, Matrix.indices,
+   Matrix.Transpose, Matrix.TransposeExtra,
+   Matrix.transpose, MatrixClass.adjoint,
+   Matrix.height, Matrix.width, Matrix.extent,
+   Matrix.Box, Matrix.BoxExtra, Matrix.indices,
    ArrMatrix.reshape,
    ArrMatrix.mapShape,
    caseTallWide,
    fromScalar, toScalar,
    fromList,
-   MatrixClass.Unpack, MatrixClass.toFull, OmniMatrix.unpack,
+   MatrixClass.Unpack, MatrixClass.UnpackExtra,
+   MatrixClass.toFull, MatrixClass.unpack,
    Matrix.mapExtent, fromFull,
    asGeneral, asTall, asWide,
    tallFromGeneral, wideFromGeneral,
    generalizeTall, generalizeWide,
+   MatrixClass.MapSquareSize, MatrixClass.MapSize,
    MatrixClass.mapHeight, MatrixClass.mapWidth,
    MatrixClass.mapSquareSize,
    Quadratic.identity,
@@ -70,7 +77,7 @@
    map,
    MatrixClass.Complex, MatrixClass.conjugate,
    MatrixClass.fromReal, MatrixClass.toComplex,
-   MatrixClass.SquareShape, MatrixClass.toSquare,
+   MatrixClass.SquareShape, MatrixClass.SquareShapeExtra, MatrixClass.toSquare,
    OmniMatrix.identityFromShape,
    MatrixClass.identityFrom,
    MatrixClass.takeDiagonal, MatrixClass.trace,
@@ -85,26 +92,33 @@
    Full.multiply,
    Full.multiplyVector,
 
-   Matrix.ToQuadratic,
+   Matrix.ToQuadratic(..),
 
-   ArrMatrix.zero, ArrMatrix.negate,
-   ArrMatrix.scale, ArrMatrix.scaleReal, ArrMatrix.scaleRealReal,
-   (ArrMatrix..*#),
-   ArrMatrix.add, ArrMatrix.sub,
-   (ArrMatrix.#+#), (ArrMatrix.#-#),
+   MatrixClass.Homogeneous, MatrixClass.HomogeneousExtra,
+   ArrMatrix.zero, MatrixClass.zeroFrom, MatrixClass.negate,
+   MatrixClass.Scale, MatrixClass.ScaleExtra, MatrixClass.scale,
+   MatrixClass.scaleReal, MatrixClass.scaleRealReal,
+   (MatrixClass..*#),
+   MatrixClass.Additive, MatrixClass.Subtractive,
+   MatrixClass.AdditiveExtra, MatrixClass.SubtractiveExtra,
+   MatrixClass.add, MatrixClass.sub,
+   (MatrixClass.#+#), (MatrixClass.#-#),
    Multiply.Multiply, (Multiply.#*#),
-   Multiply.MultiplyVector, (Multiply.#*|), (Multiply.-*#),
-   Multiply.MultiplySquare, multiplySquare,
-   Multiply.Power, Multiply.square,
+   Multiply.MultiplyVector, Multiply.MultiplyVectorExtra,
+   (Multiply.#*|), (Multiply.-*#),
+   Multiply.MultiplySquare, Multiply.MultiplySquareExtra, multiplySquare,
+   Matrix.MultiplySame(multiplySame), Matrix.MultiplySameExtra,
+   Multiply.Power, Multiply.PowerExtra, Multiply.square,
    Multiply.power, Multiply.powers, Multiply.powers1,
    (Multiply.##*#), (Multiply.#*##),
    Indexed.Indexed, (Indexed.#!),
 
-   Divide.Determinant, Divide.determinant,
-   Divide.Solve, Divide.solve, Divide.solveLeft, Divide.solveRight,
+   Divide.Determinant, Divide.DeterminantExtra, Divide.determinant,
+   Divide.Solve, Divide.SolveExtra,
+   Divide.solve, Divide.solveLeft, Divide.solveRight,
    (Divide.##/#), (Divide.#\##),
    Divide.solveVector, (Divide.-/#), (Divide.#\|),
-   Divide.Inverse, Divide.inverse,
+   Divide.Inverse, Divide.InverseExtra, Divide.inverse,
    Mod.Transposition(..),
    ) where
 
@@ -120,7 +134,7 @@
 import qualified Numeric.LAPACK.Matrix.Basic as Basic
 import qualified Numeric.LAPACK.Matrix.Array.Basic as OmniMatrix
 import qualified Numeric.LAPACK.Matrix.Array.Private as ArrMatrix
-import qualified Numeric.LAPACK.Matrix.Type as Matrix
+import qualified Numeric.LAPACK.Matrix.Type.Private as Matrix
 import qualified Numeric.LAPACK.Matrix.Plain as Plain
 import qualified Numeric.LAPACK.Matrix.Modifier as Mod
 import qualified Numeric.LAPACK.Matrix.Divide as Divide
@@ -720,7 +734,8 @@
 
 
 multiplySquare ::
-   (Multiply.MultiplySquare typ xl xu) =>
+   (Multiply.MultiplySquare typ) =>
+   (Multiply.MultiplySquareExtra typ xl, Multiply.MultiplySquareExtra typ xu) =>
    (Omni.Strip lower, Omni.Strip upper) =>
    (Extent.Measure meas, Extent.C vert, Extent.C horiz,
     Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
diff --git a/src/Numeric/LAPACK/Matrix/Array.hs b/src/Numeric/LAPACK/Matrix/Array.hs
--- a/src/Numeric/LAPACK/Matrix/Array.hs
+++ b/src/Numeric/LAPACK/Matrix/Array.hs
@@ -1,6 +1,7 @@
 module Numeric.LAPACK.Matrix.Array (
    module Numeric.LAPACK.Matrix.Array.Private,
    OmniMatrix.identityFromShape,
+   OmniMatrix.unpack,
    ) where
 
 import qualified Numeric.LAPACK.Matrix.Array.Basic as OmniMatrix
diff --git a/src/Numeric/LAPACK/Matrix/Array/Multiply.hs b/src/Numeric/LAPACK/Matrix/Array/Multiply.hs
--- a/src/Numeric/LAPACK/Matrix/Array/Multiply.hs
+++ b/src/Numeric/LAPACK/Matrix/Array/Multiply.hs
@@ -1,5 +1,6 @@
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE ConstraintKinds #-}
 {-# LANGUAGE GADTs #-}
 module Numeric.LAPACK.Matrix.Array.Multiply where
 
@@ -14,7 +15,7 @@
 import qualified Numeric.LAPACK.Matrix.Triangular.Basic as TriangularBasic
 import qualified Numeric.LAPACK.Matrix.Square.Basic as SquareBasic
 import qualified Numeric.LAPACK.Matrix.Basic as FullBasic
-import qualified Numeric.LAPACK.Matrix.Type as Matrix
+import qualified Numeric.LAPACK.Matrix.Type.Private as Matrix
 import qualified Numeric.LAPACK.Matrix.Array.Basic as OmniMatrix
 import qualified Numeric.LAPACK.Matrix.Array.Unpacked as Unpacked
 import qualified Numeric.LAPACK.Matrix.Array.Private as ArrMatrix
@@ -30,7 +31,7 @@
          (Filled, Empty, Bands, Packed, Unpacked, UnaryProxy, natFromProxy)
 import Numeric.LAPACK.Matrix.Array.Banded (Banded)
 import Numeric.LAPACK.Matrix.Array.Private (ArrayMatrix, diagTag)
-import Numeric.LAPACK.Matrix.Type (Matrix)
+import Numeric.LAPACK.Matrix.Type.Private (Matrix)
 import Numeric.LAPACK.Matrix.Modifier (Transposition(NonTransposed,Transposed))
 import Numeric.LAPACK.Matrix.Shape.Omni (StripSingleton(StripBands, StripFilled))
 import Numeric.LAPACK.Vector (Vector)
@@ -816,6 +817,8 @@
 
 unifyFactors ::
    (Matrix.Box typA, Matrix.Box typB) =>
+   (Matrix.BoxExtra typA xlA, Matrix.BoxExtra typA xuA) =>
+   (Matrix.BoxExtra typB xlB, Matrix.BoxExtra typB xuB) =>
    (Extent.Measure measA, Extent.C vertA, Extent.C horizA,
     Extent.Measure measB, Extent.C vertB, Extent.C horizB) =>
    (ExtentPriv.MultiplyMeasure measA measB ~ measC) =>
diff --git a/src/Numeric/LAPACK/Matrix/Array/Private.hs b/src/Numeric/LAPACK/Matrix/Array/Private.hs
--- a/src/Numeric/LAPACK/Matrix/Array/Private.hs
+++ b/src/Numeric/LAPACK/Matrix/Array/Private.hs
@@ -3,6 +3,8 @@
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE GADTs #-}
 {-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE ConstraintKinds #-}
+{-# LANGUAGE UndecidableInstances #-}
 {-# LANGUAGE StandaloneDeriving #-}
 module Numeric.LAPACK.Matrix.Array.Private (
    Matrix(Array),
@@ -56,20 +58,21 @@
    liftOmni1,
    liftOmni2,
 
-   Homogeneous, Scale, zero, negate, scaleReal, scale, scaleRealReal, (.*#),
+   Homogeneous, Scale, zero, negate, scaleReal, scale, scaleRealReal,
    order, forceOrder, adaptOrder,
-   Additive, add, (#+#),
-   Subtractive, sub, (#-#),
+   Additive, add,
+   Subtractive, sub,
 
    MapExtent, mapExtent,
    ) where
 
 import qualified Numeric.LAPACK.Matrix.Plain.Class as ArrClass
-import qualified Numeric.LAPACK.Matrix.Type as Matrix
+import qualified Numeric.LAPACK.Matrix.Type.Private as Matrix
 import qualified Numeric.LAPACK.Matrix.BandedHermitian.Basic as BandedHermitian
 import qualified Numeric.LAPACK.Matrix.Banded.Basic as Banded
 import qualified Numeric.LAPACK.Matrix.Triangular.Basic as Triangular
 import qualified Numeric.LAPACK.Matrix.Mosaic.Basic as Mosaic
+import qualified Numeric.LAPACK.Matrix.Plain.Format as ArrFormat
 import qualified Numeric.LAPACK.Matrix.Plain as Plain
 import qualified Numeric.LAPACK.Matrix.Basic as Basic
 import qualified Numeric.LAPACK.Matrix.Shape.Omni as Omni
@@ -78,13 +81,13 @@
 import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
 import qualified Numeric.LAPACK.Matrix.Extent.Strict as ExtentStrict
 import qualified Numeric.LAPACK.Vector as Vector
-import Numeric.LAPACK.Matrix.Plain.Format (formatArray)
 import Numeric.LAPACK.Matrix.Shape.Omni (Omni, Arbitrary)
-import Numeric.LAPACK.Matrix.Layout.Private (Filled, Bands, Packed, Unpacked)
+import Numeric.LAPACK.Matrix.Layout.Private
+         (Filled, Bands, GetBands, Packed, Unpacked)
 import Numeric.LAPACK.Matrix.Extent.Private (Extent, Shape, Size, Big, Small)
-import Numeric.LAPACK.Matrix.Type (Matrix)
+import Numeric.LAPACK.Matrix.Type.Private (Matrix)
 import Numeric.LAPACK.Vector (Vector)
-import Numeric.LAPACK.Scalar (RealOf)
+import Numeric.LAPACK.Scalar (RealOf, conjugate)
 
 import qualified Type.Data.Num.Unary as Unary
 import qualified Type.Data.Bool as TBool
@@ -101,6 +104,8 @@
 
 import Foreign.Storable (Storable)
 
+import GHC.Exts (Constraint)
+
 import Prelude hiding (negate)
 
 
@@ -147,9 +152,11 @@
    rnf (Array arr) = DeepSeq.rnf arr
 
 instance Matrix.Box (Array pack property) where
+   type BoxExtra (Array pack property) extra = extra ~ ()
    extent (Array arr) = Omni.extent $ Array.shape arr
 
 instance Matrix.Transpose (Array pack property) where
+   type TransposeExtra (Array pack property) extra = extra ~ ()
    transpose a@(Array _) =
       case shape a of
          Omni.Full _ -> liftUnpacked1 Basic.transpose a
@@ -197,14 +204,14 @@
 subBandsSingleton ::
    (Unary.Natural sub) =>
    ArrayMatrix pack property
-      (Layout.Bands sub) upper meas vert horiz height width a ->
+      (Bands sub) upper meas vert horiz height width a ->
    Unary.HeadSingleton sub
 subBandsSingleton _ = Unary.headSingleton
 
 superBandsSingleton ::
    (Unary.Natural super) =>
    ArrayMatrix pack property
-      lower (Layout.Bands super) meas vert horiz height width a ->
+      lower (Bands super) meas vert horiz height width a ->
    Unary.HeadSingleton super
 superBandsSingleton _ = Unary.headSingleton
 
@@ -418,18 +425,38 @@
 unliftColumn order_ = Basic.unliftColumn order_ . unlift1
 
 
-instance Matrix.FormatMatrix (Array pack property) where
-   formatMatrix fmt a@(Array _) =
+instance Matrix.Format (Array pack property) where
+   type FormatExtra (Array pack property) extra = ()
+   format = Matrix.formatWithLayout
+
+instance Matrix.Layout (Array pack property) where
+   type LayoutExtra (Array pack property) extra = ()
+   layout a@(Array _) =
       case shape a of
          Omni.Full fullShape ->
-            formatArray fmt $ Array.reshape fullShape $ unwrap a
-         Omni.UpperTriangular _ -> formatArray fmt $ toVector a
-         Omni.LowerTriangular _ -> formatArray fmt $ toVector a
-         Omni.Symmetric _ -> formatArray fmt $ toVector a
-         Omni.Hermitian _ -> formatArray fmt $ toVector a
-         Omni.Banded _ -> formatArray fmt $ toVector a
-         Omni.UnitBandedTriangular _ -> formatArray fmt $ toVector a
-         Omni.BandedHermitian _ -> formatArray fmt $ toVector a
+            ArrFormat.arrayFromList2 (Matrix.extent a) $
+            ArrFormat.layoutFull $ Array.reshape fullShape $ unwrap a
+         Omni.UpperTriangular _ ->
+            ArrFormat.incompleteArrayFromList2 (Matrix.extent a) $
+            ArrFormat.layoutTriangular $ toVector a
+         Omni.LowerTriangular _ ->
+            ArrFormat.incompleteArrayFromList2 (Matrix.extent a) $
+            ArrFormat.layoutTriangular $ toVector a
+         Omni.Symmetric _ ->
+            ArrFormat.splitArrayFromList2 (Matrix.extent a) $
+            ArrFormat.layoutMirrored id $ toVector a
+         Omni.Hermitian _ ->
+            ArrFormat.splitArrayFromList2 (Matrix.extent a) $
+            ArrFormat.layoutMirrored conjugate $ toVector a
+         Omni.Banded _ ->
+            ArrFormat.incompleteArrayFromList2 (Matrix.extent a) $
+            ArrFormat.layoutBanded $ toVector a
+         Omni.UnitBandedTriangular _ ->
+            ArrFormat.incompleteArrayFromList2 (Matrix.extent a) $
+            ArrFormat.layoutBanded $ toVector a
+         Omni.BandedHermitian _ ->
+            ArrFormat.incompleteSplitArrayFromList2 (Matrix.extent a) $
+            ArrFormat.layoutBandedHermitian $ toVector a
 
 
 packTag ::
@@ -445,8 +472,9 @@
 diagTag _ = Omni.autoDiag
 
 instance
-   (Layout.Packing pack, Omni.TriDiag diag, xl ~ (), xu ~ ()) =>
-      Matrix.MultiplySame (Array pack diag) xl xu where
+   (Layout.Packing pack, Omni.TriDiag diag) =>
+      Matrix.MultiplySame (Array pack diag) where
+   type MultiplySameExtra (Array pack diag) extra = extra ~ ()
    multiplySame a =
       case Omni.powerSingleton $ shape a of
          Omni.PowerIdentity -> \b ->
@@ -533,18 +561,15 @@
 instance Scale Arbitrary where
 instance Scale Omni.Symmetric where
 
-scale, (.*#) ::
+scale ::
    (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
    (Scale property, Shape.C height, Shape.C width, Class.Floating a) =>
    a ->
    ArrayMatrix pack property lower upper meas vert horiz height width a ->
    ArrayMatrix pack property lower upper meas vert horiz height width a
 scale a (Array v) = Array $ Vector.scale a v
-(.*#) = scale
 
-infixl 7 .*#
 
-
 order ::
    ArrayMatrix pack property lower upper meas vert horiz height width a ->
    Layout.Order
@@ -591,9 +616,9 @@
 instance (TBool.C neg, TBool.C zero, neg ~ pos) =>
             Subtractive (Omni.Hermitian neg zero pos) where
 
-infixl 6 #+#, #-#, `add`, `sub`
+infixl 6 `add`, `sub`
 
-add, (#+#) ::
+add ::
    (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
    (Additive property,
     Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>
@@ -602,7 +627,7 @@
    ArrayMatrix pack property lower upper meas vert horiz height width a
 add a b = liftOmni2 Vector.add (adaptOrder b a) b
 
-sub, (#-#) ::
+sub ::
    (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
    (Subtractive property,
     Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>
@@ -611,10 +636,7 @@
    ArrayMatrix pack property lower upper meas vert horiz height width a
 sub a b = liftOmni2 Vector.sub (adaptOrder b a) b
 
-(#+#) = add
-(#-#) = sub
 
-
 liftOmni1 ::
    (OmniArray packA propA lowerA upperA measA vertA horizA heightA widthA a ->
     OmniArray packB propB lowerB upperB measB vertB horizB heightB widthB b)
@@ -665,22 +687,29 @@
 
 
 instance
-   (MapExtent pack property lower upper, xl ~ (), xu ~ ()) =>
-      Matrix.MapExtent (Array pack property) xl xu lower upper where
+   (MapExtent pack, property ~ Arbitrary) =>
+      Matrix.MapExtent (Array pack property) where
+   type MapExtentExtra (Array pack property) extra = extra ~ ()
+   type MapExtentStrip (Array pack property) strip = MapExtentStrip pack strip
    mapExtent = mapExtent
 
-class MapExtent pack property lower upper where
+class MapExtent pack where
+   type MapExtentStrip pack strip :: Constraint
    mapExtent ::
+      (property ~ Arbitrary) =>
       (Extent.Measure measA, Extent.C vertA, Extent.C horizA) =>
       (Extent.Measure measB, Extent.C vertB, Extent.C horizB) =>
+      (MapExtentStrip pack lower, MapExtentStrip pack upper) =>
       ExtentStrict.Map measA vertA horizA measB vertB horizB height width ->
       ArrayMatrix pack property lower upper measA vertA horizA height width a ->
       ArrayMatrix pack property lower upper measB vertB horizB height width a
 
-instance MapExtent Unpacked Arbitrary Filled Filled where
+instance MapExtent Unpacked where
+   type MapExtentStrip Unpacked strip = strip ~ Filled
    mapExtent = lift1 . Plain.mapExtent . ExtentStrict.apply
 
-instance
-   (Unary.Natural sub, Unary.Natural super) =>
-      MapExtent Packed Arbitrary (Bands sub) (Bands super) where
+instance MapExtent Packed where
+   type MapExtentStrip Packed strip =
+            (strip ~ Bands (GetBands strip),
+             Unary.Natural (GetBands strip))
    mapExtent = lift1 . Banded.mapExtent . ExtentStrict.apply
diff --git a/src/Numeric/LAPACK/Matrix/Block.hs b/src/Numeric/LAPACK/Matrix/Block.hs
--- a/src/Numeric/LAPACK/Matrix/Block.hs
+++ b/src/Numeric/LAPACK/Matrix/Block.hs
@@ -1,232 +1,79 @@
-{-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE TypeOperators #-}
-{-# LANGUAGE GADTs #-}
-{-# LANGUAGE EmptyDataDecls #-}
-{-# LANGUAGE FlexibleInstances #-}
-{-# LANGUAGE UndecidableInstances #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
-module Numeric.LAPACK.Matrix.Block where
-
-import qualified Numeric.LAPACK.Matrix as Matrix
-import qualified Numeric.LAPACK.Matrix.Class as MatrixClass
-import qualified Numeric.LAPACK.Matrix.Divide as Divide
-import qualified Numeric.LAPACK.Matrix.Multiply as Multiply
-import qualified Numeric.LAPACK.Matrix.Square as Square
-import qualified Numeric.LAPACK.Matrix.Array.Unpacked as Unpacked
-import qualified Numeric.LAPACK.Matrix.Array.Private as ArrMatrix
-import qualified Numeric.LAPACK.Matrix.Shape.Omni as Omni
-import qualified Numeric.LAPACK.Matrix.Layout as Layout
-import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
-import qualified Numeric.Netlib.Class as Class
-import Numeric.LAPACK.Matrix.Divide (determinant)
-import Numeric.LAPACK.Matrix.Type (Matrix, Quadratic, extent, squareSize)
-import Numeric.LAPACK.Matrix.Layout.Private (Filled)
-import Numeric.LAPACK.Matrix.Extent.Private (Size, Big)
-import Numeric.LAPACK.Matrix ((#*#), (#-#), (===), (|||))
-import Numeric.LAPACK.Vector ((|+|))
-import Numeric.LAPACK.Shape.Private (Unchecked(Unchecked), deconsUnchecked)
-
-import qualified Data.Array.Comfort.Storable as Array
-import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Shape ((::+)((::+)))
+{- |
+Matrices that are assembled from smaller matrices.
 
+We can nest block matrices,
+but we still not have appropriate type classes for their multiplications.
+E.g. a Square matrix with more than 2x2 blocks would have
+the top-level structure:
 
+> Quadratic   Block.Beside
+> Block.Above Block.Square
 
-data Square typ00 typ01 typ10 typ11
-data instance
-   Matrix (Square typ00 typ01 typ10 typ11) xl xu
-      lower upper meas vert horiz height width a where
-   Square ::
-      (Extent.Measure measOff, Extent.C vertOff, Extent.C horizOff,
-       Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1) =>
-      Quadratic typ00 xl00 xu00 Filled Filled sh0 a ->
-      Matrix typ01 xl01 xu01 Filled Filled measOff vertOff horizOff sh0 sh1 a ->
-      Matrix typ10 xl10 xu10 Filled Filled measOff horizOff vertOff sh1 sh0 a ->
-      Quadratic typ11 xl11 xu11 Filled Filled sh1 a ->
-      Quadratic
-         (Square typ00 typ01 typ10 typ11)
-         (xl00,xl10,xl01,xl11) (xu00,xu01,xu10,xu11)
-         Filled Filled (sh0::+sh1) a
+Thus we would need
+e.g. multiplication Beside times Above with Quadratic result
+and multiplication Beside times Square with Above result.
+-}
+module Numeric.LAPACK.Matrix.Block (
+   Matrix(Diagonal, Above, Beside, Square, Upper, Lower, Symmetric),
+   Diagonal,
+   Above, Beside, Block.aboveFromFull, Block.besideFromFull,
+   Square,
+   LowerTriangular, UpperTriangular,
+   Symmetric, Block.squareFromSymmetric, Block.schurComplement,
+   ) where
 
-instance
-   (Matrix.Box typ00, Matrix.Box typ11) =>
-      Matrix.Box (Square typ00 typ01 typ10 typ11) where
-   extent (Square a _b _c d) = Extent.square (squareSize a ::+ squareSize d)
+import qualified Numeric.LAPACK.Matrix.Block.Private as Block
+import Numeric.LAPACK.Matrix.Type.Private (Matrix, Quadratic)
+import Numeric.LAPACK.Matrix.Layout.Private (Filled)
+import Numeric.LAPACK.Matrix.Extent.Private (Size, Big)
 
-transposeSquare ::
-   (Matrix.Transpose typ00, Matrix.Transpose typ01) =>
-   (Matrix.Transpose typ10, Matrix.Transpose typ11) =>
-   (Class.Floating a) =>
-   Matrix (Square typ00 typ01 typ10 typ11) xl xu
-      lower upper meas vert horiz height width a ->
-   Matrix (Square typ00 typ10 typ01 typ11) xu xl
-      upper lower meas horiz vert width height a
-transposeSquare (Square a b c d) =
-   Square
-      (Matrix.transpose a) (Matrix.transpose c)
-      (Matrix.transpose b) (Matrix.transpose d)
+import Data.Array.Comfort.Shape ((::+))
 
-instance
-   (xl ~ (xl00,xl10,xl01,xl11),
-    xu ~ (xu00,xu01,xu10,xu11),
-    Multiply.MultiplyVector typ00 xl00 xu00,
-    Multiply.MultiplyVector typ01 xl01 xu01,
-    Multiply.MultiplyVector typ10 xl10 xu10,
-    Multiply.MultiplyVector typ11 xl11 xu11) =>
-      Multiply.MultiplyVector (Square typ00 typ01 typ10 typ11) xl xu where
-   matrixVector (Square a b c d) x =
-      let (x0,x1) = Array.split x
-      in Array.append
-            (Multiply.matrixVector a x0 |+| Multiply.matrixVector b x1)
-            (Multiply.matrixVector c x0 |+| Multiply.matrixVector d x1)
-   vectorMatrix x (Square a b c d) =
-      let (x0,x1) = Array.split x
-      in Array.append
-            (Multiply.vectorMatrix x0 a |+| Multiply.vectorMatrix x1 c)
-            (Multiply.vectorMatrix x0 b |+| Multiply.vectorMatrix x1 d)
+import qualified Type.Data.Bool as TBool
 
 
-type TypeFull = ArrMatrix.Array Layout.Unpacked Omni.Arbitrary
-
-schurComplement ::
-   (Divide.Solve typ11 xl11 xu11,
-    Omni.Strip lower, Omni.Strip upper, Class.Floating a) =>
-   Quadratic
-      (Square TypeFull TypeFull TypeFull typ11)
-      ((),(),(),xl11) ((),(),(),xu11)
-      lower upper (sh0::+sh1) a ->
-   Square.Square sh0 a
-schurComplement (Square a b c d) =
-   Unpacked.fillBoth a
-   #-#
-   Square.fromFull
-      (Matrix.fromFull b #*# Matrix.fromFull (Divide.solveRight d c))
-
-{- |
-Requires that the right bottom sub-matrix is invertible.
--}
-instance
-   (xl ~ ((),(),(),xl11),
-    xu ~ ((),(),(),xu11),
-    typ00 ~ TypeFull,
-    typ01 ~ TypeFull,
-    typ10 ~ TypeFull,
-    Divide.Solve typ11 xl11 xu11,
-    Divide.Determinant typ11 xl11 xu11) =>
-      Divide.Determinant (Square typ00 typ01 typ10 typ11) xl xu where
-   determinant sq@(Square _a _b _c d) =
-      determinant d * determinant (schurComplement sq)
-
+type Diagonal typ0 xl0 xu0 typ1 xl1 xu1 lower upper sh0 sh1 =
+      Quadratic
+         (Block.Diagonal typ0 typ1) (xl0,xl1) (xu0,xu1)
+         lower upper (sh0::+sh1)
 
-withoutHeightCheck ::
-   (MatrixClass.MapSize typ0, Extent.C vert0, Extent.C horiz0) =>
-   (MatrixClass.MapSize typ1, Extent.C vert1, Extent.C horiz1) =>
-   (MatrixClass.MapSize typ2, Extent.C vert2, Extent.C horiz2) =>
-   (Matrix typ0 xl0 xu0 lower0 upper0 Size vert0 horiz0 ~ matrix0) =>
-   (Matrix typ1 xl1 xu1 lower1 upper1 Size vert1 horiz1 ~ matrix1) =>
-   (Matrix typ2 xl2 xu2 lower2 upper2 Size vert2 horiz2 ~ matrix2) =>
-   (Shape.C height0, Shape.C width0) =>
-   (Shape.C height1, Shape.C width1) =>
-   (Shape.C height2, Shape.C width2) =>
-   (matrix0 (Unchecked height0) width0 a0 ->
-    matrix1 (Unchecked height1) width1 a1 ->
-    matrix2 (Unchecked height2) width2 a2) ->
-   matrix0 height0 width0 a0 ->
-   matrix1 height1 width1 a1 ->
-   matrix2 height2 width2 a2
-withoutHeightCheck op a b =
-   Matrix.mapHeight deconsUnchecked $
-   Matrix.mapHeight Unchecked a `op` Matrix.mapHeight Unchecked b
+type Square typ00 xl00 xu00 typ01 xl01 xu01 typ10 xl10 xu10 typ11 xl11 xu11
+         measOff vertOff horizOff sh0 sh1 =
+      Quadratic
+         (Block.Square typ00 measOff vertOff horizOff typ11)
+         (xl00,xl11,(typ10,xl10,xu10))
+         (xu00,xu11,(typ01,xu01,xl01))
+         Filled Filled (sh0::+sh1)
 
-withoutWidthCheck ::
-   (MatrixClass.MapSize typ0, Extent.C vert0, Extent.C horiz0) =>
-   (MatrixClass.MapSize typ1, Extent.C vert1, Extent.C horiz1) =>
-   (MatrixClass.MapSize typ2, Extent.C vert2, Extent.C horiz2) =>
-   (Matrix typ0 xl0 xu0 lower0 upper0 Size vert0 horiz0 ~ matrix0) =>
-   (Matrix typ1 xl1 xu1 lower1 upper1 Size vert1 horiz1 ~ matrix1) =>
-   (Matrix typ2 xl2 xu2 lower2 upper2 Size vert2 horiz2 ~ matrix2) =>
-   (Shape.C height0, Shape.C width0) =>
-   (Shape.C height1, Shape.C width1) =>
-   (Shape.C height2, Shape.C width2) =>
-   (matrix0 height0 (Unchecked width0) a0 ->
-    matrix1 height1 (Unchecked width1) a1 ->
-    matrix2 height2 (Unchecked width2) a2) ->
-   matrix0 height0 width0 a0 ->
-   matrix1 height1 width1 a1 ->
-   matrix2 height2 width2 a2
-withoutWidthCheck op a b =
-   Matrix.mapWidth deconsUnchecked $
-   Matrix.mapWidth Unchecked a `op` Matrix.mapWidth Unchecked b
+type Above typ0 xl0 xu0 typ1 xl1 xu1 horiz height0 height1 width =
+      Matrix
+         (Block.Append typ0 typ1 height0 height1)
+         (xl0,xl1,TBool.False) (xu0,xu1,TBool.True)
+         Filled Filled Size Big horiz (height0::+height1) width
 
-{- |
-Requires that the right bottom sub-matrix is invertible.
--}
-instance
-   (xl ~ ((),(),(),xl11),
-    xu ~ ((),(),(),xu11),
-    typ00 ~ TypeFull,
-    typ01 ~ TypeFull,
-    typ10 ~ TypeFull,
-    Divide.Solve typ11 xl11 xu11) =>
-      Divide.Solve (Square typ00 typ01 typ10 typ11) xl xu where
-   solveRight sq@(Square _a b c d) x =
-      let x0 = Matrix.takeTop    $ Matrix.fromFull x
-          x1 = Matrix.takeBottom $ Matrix.fromFull x
-          xComplement =
-            withoutWidthCheck (#-#) x0 $
-               Matrix.fromFull b #*# Divide.solveRight d x1
-          y = Divide.solveRight (schurComplement sq) xComplement
-      in -- ToDo: does it always has correct order?
-         ArrMatrix.reshape (ArrMatrix.shape x) $
-         withoutWidthCheck (===) y
-            (Divide.solveRight d $
-             withoutWidthCheck (#-#) x1 $ Matrix.fromFull c #*# y)
-   solveLeft x sq@(Square _a b c d) =
-      let x0 = Matrix.takeLeft  $ Matrix.fromFull x
-          x1 = Matrix.takeRight $ Matrix.fromFull x
-          xComplement =
-            withoutHeightCheck (#-#) x0 $
-               Divide.solveLeft x1 d #*# Matrix.fromFull c
-          y = Divide.solveLeft xComplement (schurComplement sq)
-      in -- ToDo: does it always has correct order?
-         ArrMatrix.reshape (ArrMatrix.shape x) $
-         withoutHeightCheck (|||) y
-            (Divide.solveLeft
-               (withoutHeightCheck (#-#) x1 $ y #*# Matrix.fromFull b)
-               d)
+type Beside typ0 xl0 xu0 typ1 xl1 xu1 vert height width0 width1 =
+      Matrix
+         (Block.Append typ0 typ1 width0 width1)
+         (xl0,xl1,TBool.True) (xu0,xu1,TBool.False)
+         Filled Filled Size vert Big height (width0::+width1)
 
-{-
-instance
-   (xl ~ ((),(),(),xl11),
-    xu ~ ((),(),(),xu11),
-    typ00 ~ TypeFull,
-    typ01 ~ TypeFull,
-    typ10 ~ TypeFull,
-    Divide.Inverse typ11 xl11 xu11) =>
-      Divide.Inverse (Square typ00 typ01 typ10 typ11) xl xu where
-   inverse (Square a b c d) =
-      let as = a #-# b #*# Divide.solveRight d c
-          bdinv = Divide.solveLeft b d
-          dinvc = Divide.solveRight d c
-          br = Divide.solveRight as bdinv
-          cr = Divide.solveLeft dinvc as
-      in Square
-            (inverse as)
-            (Matrix.negate br)
-            (Matrix.negate cr)
-            (inverse d #+# br #*# bdinv)
--}
+type UpperTriangular typ0 xl0 xu0 typOff xlOff xuOff typ1 xl1 xu1
+         lower sh0 sh1 =
+      Quadratic
+         (Block.Triangular typ0 typOff typ1)
+         (xl0,xlOff,xl1,TBool.False) (xu0,xuOff,xu1,TBool.True)
+         lower Filled (sh0::+sh1)
 
+type LowerTriangular typ0 xl0 xu0 typOff xlOff xuOff typ1 xl1 xu1
+         upper sh0 sh1 =
+      Quadratic
+         (Block.Triangular typ0 typOff typ1)
+         (xl0,xlOff,xl1,TBool.True) (xu0,xuOff,xu1,TBool.False)
+         Filled upper (sh0::+sh1)
 
-data Mosaic typ0 typOff typ1
-data instance
-   Matrix (Mosaic typ0 typOff typ1) xl xu
-      lower upper meas vert horiz height width a where
-   Mosaic ::
-      Quadratic typ0 xl0 xu0 lower upper sh0 a ->
-      Matrix typOff xlOff xuOff Filled Filled Size Big Big sh0 sh1 a ->
-      Quadratic typ1 xl1 xu1 lower upper sh1 a ->
+type Symmetric typ0 xl0 xu0 typOff xlOff xuOff typ1 xl1 xu1 sh0 sh1 =
       Quadratic
-         (Mosaic typ0 typOff typ1)
-         (xl0,xlOff,xl1) (xu0,xuOff,xu1)
-         lower upper (sh0::+sh1) a
+         (Block.Symmetric typ0 typOff xlOff xuOff typ1)
+         (xl0,xl1) (xu0,xu1)
+         Filled Filled (sh0::+sh1)
diff --git a/src/Numeric/LAPACK/Matrix/Block/Private.hs b/src/Numeric/LAPACK/Matrix/Block/Private.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Matrix/Block/Private.hs
@@ -0,0 +1,1604 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE EmptyDataDecls #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE ConstraintKinds #-}
+{-# LANGUAGE UndecidableInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE StandaloneDeriving #-}
+module Numeric.LAPACK.Matrix.Block.Private (
+   Matrix(Diagonal, Above, Beside, Square, Upper, Lower, Symmetric),
+   Diagonal,
+   Append, aboveFromFull, besideFromFull,
+   Square,
+   Triangular,
+   Symmetric, squareFromSymmetric, schurComplement,
+   ) where
+
+import qualified Numeric.LAPACK.Matrix as Matrix
+import qualified Numeric.LAPACK.Matrix.Class as MatrixClass
+import qualified Numeric.LAPACK.Matrix.Divide as Divide
+import qualified Numeric.LAPACK.Matrix.Multiply as Multiply
+import qualified Numeric.LAPACK.Matrix.Square as Square
+import qualified Numeric.LAPACK.Matrix.Array.Private as ArrMatrix
+import qualified Numeric.LAPACK.Matrix.Plain.Format as ArrFormat
+import qualified Numeric.LAPACK.Matrix.Shape.Omni as Omni
+import qualified Numeric.LAPACK.Matrix.Layout as Layout
+import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
+import qualified Numeric.LAPACK.Vector.Private as Vector
+import qualified Numeric.LAPACK.Output as Output
+import qualified Numeric.Netlib.Class as Class
+import Numeric.LAPACK.Matrix.Divide (determinant)
+import Numeric.LAPACK.Matrix.Type.Private
+         (Matrix, Layout(layout), LayoutExtra, Quadratic, squareSize)
+import Numeric.LAPACK.Matrix.Layout.Private (Filled)
+import Numeric.LAPACK.Matrix.Extent.Private (Size, Big)
+import Numeric.LAPACK.Matrix ((#*#), (#*##), (##*#), (#+#), (#-#), (===), (|||))
+import Numeric.LAPACK.Vector (Vector, (|+|))
+import Numeric.LAPACK.Shape.Private (Unchecked(Unchecked), deconsUnchecked)
+
+import qualified Data.Array.Comfort.Storable as Array
+import qualified Data.Array.Comfort.Boxed as BoxedArray
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Shape ((::+)((::+)))
+
+import qualified Data.Stream as Stream
+import Data.Function.HT (powerAssociative)
+import Data.Tuple.HT (swap)
+
+import qualified Type.Data.Bool as TBool
+
+
+type family ShapeHead sh; type instance ShapeHead (sh0::+sh1) = sh0
+type family ShapeTail sh; type instance ShapeTail (sh0::+sh1) = sh1
+
+
+data Diagonal typ0 typ1
+data instance
+   Matrix (Diagonal typ0 typ1) xl xu
+      lower upper meas vert horiz height width a where
+   Diagonal ::
+      (Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1) =>
+      Quadratic typ0 xl0 xu0 lower upper sh0 a ->
+      Quadratic typ1 xl1 xu1 lower upper sh1 a ->
+      Quadratic
+         (Diagonal typ0 typ1) (xl0,xl1) (xu0,xu1)
+         lower upper (sh0::+sh1) a
+
+type family Diagonal0 extra; type instance Diagonal0 (x0,x1) = x0
+type family Diagonal1 extra; type instance Diagonal1 (x0,x1) = x1
+
+deriving instance
+   (Show (Quadratic typ0 (Diagonal0 xl) (Diagonal0 xu)
+            lower upper (ShapeHead height) a),
+    Show (Quadratic typ1 (Diagonal1 xl) (Diagonal1 xu)
+            lower upper (ShapeTail height) a),
+    Shape.C height, Shape.C width, Show height, Show width, Show a) =>
+   Show (Matrix (Diagonal typ0 typ1) xl xu
+            lower upper meas vert horiz height width a)
+
+instance
+   (Matrix.Box typ0, Matrix.Box typ1) =>
+      Matrix.Box (Diagonal typ0 typ1) where
+   type BoxExtra (Diagonal typ0 typ1) extra =
+         (Matrix.BoxExtra typ0 (Diagonal0 extra),
+          Matrix.BoxExtra typ1 (Diagonal1 extra))
+   extent (Diagonal a b) = Extent.square (squareSize a ::+ squareSize b)
+
+instance
+   (Matrix.Transpose typ0, Matrix.Transpose typ1) =>
+      Matrix.Transpose (Diagonal typ0 typ1) where
+   type TransposeExtra (Diagonal typ0 typ1) extra =
+         (Matrix.TransposeExtra typ0 (Diagonal0 extra),
+          Matrix.TransposeExtra typ1 (Diagonal1 extra))
+   transpose (Diagonal a b) =
+      Diagonal (Matrix.transpose a) (Matrix.transpose b)
+
+instance
+   (Matrix.Box typ0, Matrix.Box typ1) =>
+      Matrix.ToQuadratic (Diagonal typ0 typ1) where
+   heightToQuadratic a@(Diagonal _ _) = a
+   widthToQuadratic a@(Diagonal _ _) = a
+
+instance (Layout typ0, Layout typ1) => Layout (Diagonal typ0 typ1) where
+   type LayoutExtra (Diagonal typ0 typ1) extra =
+         (Matrix.BoxExtra typ0 (Diagonal0 extra),
+          Matrix.BoxExtra typ1 (Diagonal1 extra),
+          LayoutExtra typ0 (Diagonal0 extra),
+          LayoutExtra typ1 (Diagonal1 extra))
+   layout (Diagonal a b) =
+      finalizeLayout $
+      let sha = squareSize a in
+      let shb = squareSize b in
+      toRows a |||# layoutEmpty sha shb
+      ===#
+      layoutEmpty shb sha |||# toRows b
+
+instance (Layout typ0, Layout typ1) => Matrix.Format (Diagonal typ0 typ1) where
+   type FormatExtra (Diagonal typ0 typ1) extra =
+         (Matrix.BoxExtra typ0 (Diagonal0 extra),
+          Matrix.BoxExtra typ1 (Diagonal1 extra),
+          LayoutExtra typ0 (Diagonal0 extra),
+          LayoutExtra typ1 (Diagonal1 extra))
+   format = Matrix.formatWithLayout
+
+instance
+   (Matrix.SquareShape typ0, Matrix.SquareShape typ1) =>
+      Matrix.SquareShape (Diagonal typ0 typ1) where
+   type SquareShapeExtra (Diagonal typ0 typ1) extra =
+         (Matrix.SquareShapeExtra typ0 (Diagonal0 extra),
+          Matrix.SquareShapeExtra typ1 (Diagonal1 extra))
+   takeDiagonal (Diagonal a b) =
+      Array.append (MatrixClass.takeDiagonal a) (MatrixClass.takeDiagonal b)
+   identityFrom (Diagonal a b) =
+      Diagonal (MatrixClass.identityFrom a) (MatrixClass.identityFrom b)
+
+instance
+   (Matrix.Unpack typ0, Matrix.Unpack typ1) =>
+      Matrix.Unpack (Diagonal typ0 typ1) where
+   type UnpackExtra (Diagonal typ0 typ1) extra =
+         (Matrix.UnpackExtra typ0 (Diagonal0 extra),
+          Matrix.UnpackExtra typ1 (Diagonal1 extra))
+   unpack (Diagonal a0 b0) =
+      let a = MatrixClass.toFull a0 in
+      let b = MatrixClass.toFull b0 in
+      let zero shape = Matrix.asGeneral $ Matrix.zero shape in
+      let order = ArrMatrix.order b in
+      let dims = (squareSize a, squareSize b) in
+      ArrMatrix.liftUnpacked1 id $
+      Square.stack
+         a             (zero $ Omni.cons order dims)
+         (zero $ Omni.cons order $ swap dims)      b
+
+instance
+   (Matrix.Homogeneous typ0, Matrix.Homogeneous typ1) =>
+      Matrix.Homogeneous (Diagonal typ0 typ1) where
+   type HomogeneousExtra (Diagonal typ0 typ1) extra =
+         (Matrix.HomogeneousExtra typ0 (Diagonal0 extra),
+          Matrix.HomogeneousExtra typ1 (Diagonal1 extra))
+   zeroFrom (Diagonal a b) = Diagonal (Matrix.zeroFrom a) (Matrix.zeroFrom b)
+   negate (Diagonal a b) = Diagonal (Matrix.negate a) (Matrix.negate b)
+   scaleReal k (Diagonal a b) =
+      Diagonal (Matrix.scaleReal k a) (Matrix.scaleReal k b)
+
+instance
+   (Matrix.Scale typ0, Matrix.Scale typ1) =>
+      Matrix.Scale (Diagonal typ0 typ1) where
+   type ScaleExtra (Diagonal typ0 typ1) extra =
+         (Matrix.ScaleExtra typ0 (Diagonal0 extra),
+          Matrix.ScaleExtra typ1 (Diagonal1 extra))
+   scale k (Diagonal a b) = Diagonal (Matrix.scale k a) (Matrix.scale k b)
+
+instance
+   (Matrix.Additive typ0, Matrix.Additive typ1) =>
+      Matrix.Additive (Diagonal typ0 typ1) where
+   type AdditiveExtra (Diagonal typ0 typ1) extra =
+         (Matrix.AdditiveExtra typ0 (Diagonal0 extra),
+          Matrix.AdditiveExtra typ1 (Diagonal1 extra))
+   add (Diagonal a0 b0) (Diagonal a1 b1) =
+      Diagonal (Matrix.add a0 a1) (Matrix.add b0 b1)
+
+instance
+   (Matrix.Subtractive typ0, Matrix.Subtractive typ1) =>
+      Matrix.Subtractive (Diagonal typ0 typ1) where
+   type SubtractiveExtra (Diagonal typ0 typ1) extra =
+         (Matrix.SubtractiveExtra typ0 (Diagonal0 extra),
+          Matrix.SubtractiveExtra typ1 (Diagonal1 extra))
+   sub (Diagonal a0 b0) (Diagonal a1 b1) =
+      Diagonal (Matrix.sub a0 a1) (Matrix.sub b0 b1)
+
+instance
+   (Multiply.MultiplyVector typ0, Multiply.MultiplyVector typ1) =>
+      Multiply.MultiplyVector (Diagonal typ0 typ1) where
+   type MultiplyVectorExtra (Diagonal typ0 typ1) extra =
+         (Multiply.MultiplyVectorExtra typ0 (Diagonal0 extra),
+          Multiply.MultiplyVectorExtra typ1 (Diagonal1 extra))
+   matrixVector (Diagonal a b) x =
+      let (x0,x1) = Array.split x
+      in Array.append (Multiply.matrixVector a x0) (Multiply.matrixVector b x1)
+   vectorMatrix x (Diagonal a b) =
+      let (x0,x1) = Array.split x
+      in Array.append (Multiply.vectorMatrix x0 a) (Multiply.vectorMatrix x1 b)
+
+instance
+   (Multiply.MultiplySquare typ0, Multiply.MultiplySquare typ1) =>
+      Multiply.MultiplySquare (Diagonal typ0 typ1) where
+   type MultiplySquareExtra (Diagonal typ0 typ1) extra =
+         (Multiply.MultiplySquareExtra typ0 (Diagonal0 extra),
+          Multiply.MultiplySquareExtra typ1 (Diagonal1 extra))
+   squareFull (Diagonal a b) =
+      withVerticalSplit $ \(Above top bottom) ->
+         Above (Multiply.squareFull a top) (Multiply.squareFull b bottom)
+   fullSquare x (Diagonal a b) =
+      withHorizontalSplit x $ \(Beside left right) ->
+         Beside (Multiply.fullSquare left a) (Multiply.fullSquare right b)
+
+instance
+   (Matrix.MultiplySame typ0, Matrix.MultiplySame typ1) =>
+      Matrix.MultiplySame (Diagonal typ0 typ1) where
+   type MultiplySameExtra (Diagonal typ0 typ1) extra =
+         (Matrix.MultiplySameExtra typ0 (Diagonal0 extra),
+          Matrix.MultiplySameExtra typ1 (Diagonal1 extra))
+   multiplySame (Diagonal a0 b0) (Diagonal a1 b1) =
+      Diagonal (Matrix.multiplySame a0 a1) (Matrix.multiplySame b0 b1)
+
+instance
+   (Multiply.Power typ0, Multiply.Power typ1) =>
+      Multiply.Power (Diagonal typ0 typ1) where
+   type PowerExtra (Diagonal typ0 typ1) extra =
+         (Multiply.PowerExtra typ0 (Diagonal0 extra),
+          Multiply.PowerExtra typ1 (Diagonal1 extra))
+   square (Diagonal a b) = Diagonal (Multiply.square a) (Multiply.square b)
+   power n (Diagonal a b) = Diagonal (Multiply.power n a) (Multiply.power n b)
+   powers1 (Diagonal a b) =
+      Stream.zipWith Diagonal (Multiply.powers1 a) (Multiply.powers1 b)
+
+instance
+   (Divide.Determinant typ0, Divide.Determinant typ1) =>
+      Divide.Determinant (Diagonal typ0 typ1) where
+   type DeterminantExtra (Diagonal typ0 typ1) extra =
+         (Matrix.DeterminantExtra typ0 (Diagonal0 extra),
+          Matrix.DeterminantExtra typ1 (Diagonal1 extra))
+   determinant (Diagonal a b) = determinant a * determinant b
+
+instance
+   (Divide.Solve typ0, Divide.Solve typ1) =>
+      Divide.Solve (Diagonal typ0 typ1) where
+   type SolveExtra (Diagonal typ0 typ1) extra =
+         (Divide.SolveExtra typ0 (Diagonal0 extra),
+          Divide.SolveExtra typ1 (Diagonal1 extra))
+   solveRight (Diagonal a b) =
+      withVerticalSplit $ \(Above top bottom) ->
+         Above (Divide.solveRight a top) (Divide.solveRight b bottom)
+   solveLeft x (Diagonal a b) =
+      withHorizontalSplit x $ \(Beside left right) ->
+         Beside (Divide.solveLeft left a) (Divide.solveLeft right b)
+
+instance
+   (Divide.Inverse typ0, Divide.Inverse typ1) =>
+      Divide.Inverse (Diagonal typ0 typ1) where
+   type InverseExtra (Diagonal typ0 typ1) extra =
+         (Divide.InverseExtra typ0 (Diagonal0 extra),
+          Divide.InverseExtra typ1 (Diagonal1 extra))
+   inverse (Diagonal a b) = Diagonal (Divide.inverse a) (Divide.inverse b)
+
+
+data Square typ00 measOff vertOff horizOff typ11
+data instance
+   Matrix (Square typ00 measOff vertOff horizOff typ11) xl xu
+      lower upper meas vert horiz height width a where
+   Square ::
+      (Extent.Measure measOff, Extent.C vertOff, Extent.C horizOff,
+       Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1) =>
+      Quadratic typ00 xl00 xu00 Filled Filled sh0 a ->
+      Matrix typ01 xl01 xu01 Filled Filled measOff vertOff horizOff sh0 sh1 a ->
+      Matrix typ10 xl10 xu10 Filled Filled measOff horizOff vertOff sh1 sh0 a ->
+      Quadratic typ11 xl11 xu11 Filled Filled sh1 a ->
+      Quadratic
+         (Square typ00 measOff vertOff horizOff typ11)
+         (xl00,xl11,(typ10,xl10,xu10))
+         (xu00,xu11,(typ01,xu01,xl01))
+         Filled Filled (sh0::+sh1) a
+
+type family SquareType extra
+type instance SquareType (x00,x11,(typ,x10,x01)) = typ
+type family Square00 extra
+type family Square01 extra
+type family Square10 extra
+type family Square11 extra
+type instance Square00 (x00,x11,(typ,x10,x01)) = x00
+type instance Square01 (x00,x11,(typ,x10,x01)) = x01
+type instance Square10 (x00,x11,(typ,x10,x01)) = x10
+type instance Square11 (x00,x11,(typ,x10,x01)) = x11
+
+deriving instance
+   (Show (Quadratic typ00 (Square00 xl) (Square00 xu)
+            Filled Filled (ShapeHead height) a),
+    Show (Quadratic typ11 (Square11 xl) (Square11 xu)
+            Filled Filled (ShapeTail height) a),
+    Show (Matrix (SquareType xu) (Square01 xu) (Square10 xu)
+            Filled Filled measOff vertOff horizOff
+            (ShapeHead height) (ShapeTail width) a),
+    Show (Matrix (SquareType xl) (Square10 xl) (Square01 xl)
+            Filled Filled measOff horizOff vertOff
+            (ShapeTail height) (ShapeHead width) a),
+    Shape.C height, Shape.C width, Show height, Show width, Show a) =>
+   Show (Matrix (Square typ00 measOff vertOff horizOff typ11) xl xu
+            lower upper meas vert horiz height width a)
+
+instance
+   (Matrix.Box typ00, Matrix.Box typ11) =>
+      Matrix.Box (Square typ00 measOff vertOff horizOff typ11) where
+   type BoxExtra (Square typ00 measOff vertOff horizOff typ11) extra =
+         (Matrix.BoxExtra typ00 (Square00 extra),
+          Matrix.BoxExtra typ11 (Square11 extra))
+   extent (Square a _b _c d) = Extent.square (squareSize a ::+ squareSize d)
+
+instance
+   (Matrix.Transpose typ00, Matrix.Transpose typ11) =>
+      Matrix.Transpose (Square typ00 measOff vertOff horizOff typ11) where
+   type TransposeExtra (Square typ00 measOff vertOff horizOff typ11) extra =
+         (Matrix.Transpose (SquareType extra),
+          Matrix.TransposeExtra typ00 (Square00 extra),
+          Matrix.TransposeExtra (SquareType extra) (Square01 extra),
+          Matrix.TransposeExtra (SquareType extra) (Square10 extra),
+          Matrix.TransposeExtra typ11 (Square11 extra))
+   transpose (Square a b c d) =
+      Square
+         (Matrix.transpose a) (Matrix.transpose c)
+         (Matrix.transpose b) (Matrix.transpose d)
+
+instance
+   (Matrix.Box typ00, Matrix.Box typ11) =>
+      Matrix.ToQuadratic (Square typ00 measOff vertOff horizOff typ11) where
+   heightToQuadratic a@(Square _ _ _ _) = a
+   widthToQuadratic a@(Square _ _ _ _) = a
+
+instance
+   (Layout typ00, Layout typ11) =>
+      Layout (Square typ00 measOff vertOff horizOff typ11) where
+   type LayoutExtra (Square typ00 measOff vertOff horizOff typ11) extra =
+         (Matrix.BoxExtra typ00 (Square00 extra),
+          Matrix.BoxExtra typ11 (Square11 extra),
+          Matrix.BoxExtra (SquareType extra) (Square01 extra),
+          Matrix.BoxExtra (SquareType extra) (Square10 extra),
+          Layout (SquareType extra),
+          LayoutExtra typ00 (Square00 extra),
+          LayoutExtra typ11 (Square11 extra),
+          LayoutExtra (SquareType extra) (Square01 extra),
+          LayoutExtra (SquareType extra) (Square10 extra))
+   layout (Square a b c d) = finalizeLayout $
+      toRows a |||# toRows b
+      ===#
+      toRows c |||# toRows d
+
+instance
+   (Layout typ00, Layout typ11) =>
+      Matrix.Format (Square typ00 measOff vertOff horizOff typ11) where
+   type FormatExtra (Square typ00 measOff vertOff horizOff typ11) extra =
+         (Matrix.BoxExtra typ00 (Square00 extra),
+          Matrix.BoxExtra typ11 (Square11 extra),
+          Matrix.BoxExtra (SquareType extra) (Square01 extra),
+          Matrix.BoxExtra (SquareType extra) (Square10 extra),
+          Layout (SquareType extra),
+          LayoutExtra typ00 (Square00 extra),
+          LayoutExtra typ11 (Square11 extra),
+          LayoutExtra (SquareType extra) (Square01 extra),
+          LayoutExtra (SquareType extra) (Square10 extra))
+   format = Matrix.formatWithLayout
+
+instance
+   (Matrix.SquareShape typ00, Matrix.SquareShape typ11) =>
+      Matrix.SquareShape (Square typ00 measOff vertOff horizOff typ11) where
+   type SquareShapeExtra (Square typ00 measOff vertOff horizOff typ11) extra =
+         (Matrix.SquareShape (SquareType extra),
+          Matrix.SquareShapeExtra typ00 (Square00 extra),
+          Matrix.SquareShapeExtra typ11 (Square11 extra),
+          Matrix.SquareShapeExtra (SquareType extra) (Square01 extra),
+          Matrix.SquareShapeExtra (SquareType extra) (Square10 extra),
+          Matrix.Homogeneous (SquareType extra),
+          Matrix.HomogeneousExtra (SquareType extra) (Square01 extra),
+          Matrix.HomogeneousExtra (SquareType extra) (Square10 extra))
+   takeDiagonal (Square a _b _c d) =
+      Array.append (MatrixClass.takeDiagonal a) (MatrixClass.takeDiagonal d)
+   identityFrom (Square a b c d) =
+      Square
+         (MatrixClass.identityFrom a)
+         (MatrixClass.zeroFrom b)
+         (MatrixClass.zeroFrom c)
+         (MatrixClass.identityFrom d)
+
+instance
+   (Matrix.Unpack typ00, Matrix.Unpack typ11) =>
+      Matrix.Unpack (Square typ00 measOff vertOff horizOff typ11) where
+   type UnpackExtra (Square typ00 measOff vertOff horizOff typ11) extra =
+         (Matrix.Unpack (SquareType extra),
+          Matrix.UnpackExtra typ00 (Square00 extra),
+          Matrix.UnpackExtra typ11 (Square11 extra),
+          Matrix.UnpackExtra (SquareType extra) (Square01 extra),
+          Matrix.UnpackExtra (SquareType extra) (Square10 extra))
+   unpack (Square a b c d) =
+      Square.stack
+         (MatrixClass.toFull a) (MatrixClass.toFull b)
+         (MatrixClass.toFull c) (MatrixClass.toFull d)
+
+instance
+   (Matrix.Homogeneous typ00, Matrix.Homogeneous typ11) =>
+      Matrix.Homogeneous (Square typ00 measOff vertOff horizOff typ11) where
+   type HomogeneousExtra (Square typ00 measOff vertOff horizOff typ11) extra =
+         (Matrix.Homogeneous (SquareType extra),
+          Matrix.HomogeneousExtra typ00 (Square00 extra),
+          Matrix.HomogeneousExtra typ11 (Square11 extra),
+          Matrix.HomogeneousExtra (SquareType extra) (Square01 extra),
+          Matrix.HomogeneousExtra (SquareType extra) (Square10 extra))
+   zeroFrom (Square a b c d) =
+      Square
+         (MatrixClass.zeroFrom a) (MatrixClass.zeroFrom b)
+         (MatrixClass.zeroFrom c) (MatrixClass.zeroFrom d)
+   negate (Square a b c d) =
+      Square
+         (MatrixClass.negate a) (MatrixClass.negate b)
+         (MatrixClass.negate c) (MatrixClass.negate d)
+   scaleReal k (Square a b c d) =
+      Square
+         (MatrixClass.scaleReal k a) (MatrixClass.scaleReal k b)
+         (MatrixClass.scaleReal k c) (MatrixClass.scaleReal k d)
+
+instance
+   (Matrix.Scale typ00, Matrix.Scale typ11) =>
+      Matrix.Scale (Square typ00 measOff vertOff horizOff typ11) where
+   type ScaleExtra (Square typ00 measOff vertOff horizOff typ11) extra =
+         (Matrix.Scale (SquareType extra),
+          Matrix.ScaleExtra typ00 (Square00 extra),
+          Matrix.ScaleExtra typ11 (Square11 extra),
+          Matrix.ScaleExtra (SquareType extra) (Square01 extra),
+          Matrix.ScaleExtra (SquareType extra) (Square10 extra))
+   scale k (Square a b c d) =
+      Square
+         (MatrixClass.scale k a) (MatrixClass.scale k b)
+         (MatrixClass.scale k c) (MatrixClass.scale k d)
+
+instance
+   (Matrix.Additive typ00, Matrix.Additive typ11) =>
+      Matrix.Additive (Square typ00 measOff vertOff horizOff typ11) where
+   type AdditiveExtra (Square typ00 measOff vertOff horizOff typ11) extra =
+         (Matrix.Additive (SquareType extra),
+          Matrix.AdditiveExtra typ00 (Square00 extra),
+          Matrix.AdditiveExtra typ11 (Square11 extra),
+          Matrix.AdditiveExtra (SquareType extra) (Square01 extra),
+          Matrix.AdditiveExtra (SquareType extra) (Square10 extra))
+   add (Square a0 b0 c0 d0) (Square a1 b1 c1 d1) =
+      Square
+         (MatrixClass.add a0 a1) (MatrixClass.add b0 b1)
+         (MatrixClass.add c0 c1) (MatrixClass.add d0 d1)
+
+instance
+   (Matrix.Subtractive typ00, Matrix.Subtractive typ11) =>
+      Matrix.Subtractive (Square typ00 measOff vertOff horizOff typ11) where
+   type SubtractiveExtra (Square typ00 measOff vertOff horizOff typ11) extra =
+         (Matrix.Subtractive (SquareType extra),
+          Matrix.SubtractiveExtra typ00 (Square00 extra),
+          Matrix.SubtractiveExtra typ11 (Square11 extra),
+          Matrix.SubtractiveExtra (SquareType extra) (Square01 extra),
+          Matrix.SubtractiveExtra (SquareType extra) (Square10 extra))
+   sub (Square a0 b0 c0 d0) (Square a1 b1 c1 d1) =
+      Square
+         (MatrixClass.sub a0 a1) (MatrixClass.sub b0 b1)
+         (MatrixClass.sub c0 c1) (MatrixClass.sub d0 d1)
+
+instance
+   (Multiply.MultiplyVector typ00, Multiply.MultiplyVector typ11) =>
+      Multiply.MultiplyVector (Square typ00 measOff vertOff horizOff typ11)
+         where
+   type
+      MultiplyVectorExtra (Square typ00 measOff vertOff horizOff typ11) extra =
+         (Multiply.MultiplyVector (SquareType extra),
+          Multiply.MultiplyVectorExtra typ00 (Square00 extra),
+          Multiply.MultiplyVectorExtra (SquareType extra) (Square01 extra),
+          Multiply.MultiplyVectorExtra (SquareType extra) (Square10 extra),
+          Multiply.MultiplyVectorExtra typ11 (Square11 extra))
+   matrixVector (Square a b c d) x =
+      let (x0,x1) = Array.split x
+      in Array.append
+            (Multiply.matrixVector a x0 |+| Multiply.matrixVector b x1)
+            (Multiply.matrixVector c x0 |+| Multiply.matrixVector d x1)
+   vectorMatrix x (Square a b c d) =
+      let (x0,x1) = Array.split x
+      in Array.append
+            (Multiply.vectorMatrix x0 a |+| Multiply.vectorMatrix x1 c)
+            (Multiply.vectorMatrix x0 b |+| Multiply.vectorMatrix x1 d)
+
+instance
+   (Multiply.MultiplySquare typ00, Multiply.MultiplySquare typ11) =>
+      Multiply.MultiplySquare (Square typ00 measOff vertOff horizOff typ11)
+         where
+   type
+      MultiplySquareExtra (Square typ00 measOff vertOff horizOff typ11) extra =
+         (SquareType extra ~ TypeFull,
+          Square01 extra ~ (), Square10 extra ~ (),
+          Multiply.MultiplySquareExtra typ00 (Square00 extra),
+          Multiply.MultiplySquareExtra typ11 (Square11 extra))
+   squareFull (Square a b c d) =
+      withVerticalSplit $ \(Above x0 x1) ->
+         Above
+            (Multiply.squareFull a x0 #+# Matrix.fromFull b #*# x1)
+            (Matrix.fromFull c #*# x0 #+# Multiply.squareFull d x1)
+   fullSquare x (Square a b c d) =
+      withHorizontalSplit x $ \(Beside x0 x1) ->
+         Beside
+            (Multiply.fullSquare x0 a #+# x1 #*# Matrix.fromFull c)
+            (x0 #*# Matrix.fromFull b #+# Multiply.fullSquare x1 d)
+
+infixl 7 ##*##
+
+(##*##) ::
+   (Shape.C height, Shape.C width, Eq height, Eq width, Class.Floating a,
+    Extent.Measure meas0, Extent.C vert0, Extent.C horiz0,
+    Extent.Measure meas1, Extent.C vert1, Extent.C horiz1) =>
+   ArrMatrix.Full meas0 vert0 horiz0 height width a ->
+   ArrMatrix.Full meas1 vert1 horiz1 width height a ->
+   ArrMatrix.Square height a
+a ##*## b = Square.fromFull (Matrix.fromFull a #*# Matrix.fromFull b)
+
+instance
+   (typ00 ~ TypeFull, typ11 ~ TypeFull) =>
+      Matrix.MultiplySame (Square typ00 measOff vertOff horizOff typ11) where
+   type MultiplySameExtra (Square typ00 measOff vertOff horizOff typ11) extra =
+         (SquareType extra ~ TypeFull,
+          Square00 extra ~ (), Square11 extra ~ (),
+          Square01 extra ~ (), Square10 extra ~ ())
+   multiplySame (Square a0 b0 c0 d0) (Square a1 b1 c1 d1) =
+      Square
+         (a0#*#a1 #+# b0##*##c1) (a0#*##b1 #+# b0##*#d1)
+         (c0##*#a1 #+# d0#*##c1) (c0##*##b1 #+# d0#*#d1)
+
+instance
+   (typ00 ~ TypeFull, typ11 ~ TypeFull) =>
+      Multiply.Power (Square typ00 measOff vertOff horizOff typ11) where
+   type PowerExtra (Square typ00 measOff vertOff horizOff typ11) extra =
+         (SquareType extra ~ TypeFull,
+          Square00 extra ~ (), Square11 extra ~ (),
+          Square01 extra ~ (), Square10 extra ~ ())
+   square a@(Square _ _ _ _) = Matrix.multiplySame a a
+   power n m@(Square a b c d) =
+      powerAssociative Matrix.multiplySame
+         (Square
+            (Matrix.identityFrom a) (Matrix.zeroFrom b)
+            (Matrix.zeroFrom c) (Matrix.identityFrom d))
+         m n
+   powers1 a@(Square _ _ _ _) = Stream.iterate (flip Matrix.multiplySame a) a
+
+
+schurComplement ::
+   (Divide.Solve typ11,
+    Divide.SolveExtra typ11 xl11, Divide.SolveExtra typ11 xu11,
+    Class.Floating a) =>
+   Quadratic
+      (Square TypeFull measOff vertOff horizOff typ11)
+      ((),xl11,(TypeFull,(),())) ((),xu11,(TypeFull,(),()))
+      Filled Filled (sh0::+sh1) a ->
+   Square.Square sh0 a
+schurComplement (Square a b c d) = a #-# b ##*## Divide.solveRight d c
+
+{- |
+Requires that the bottom right sub-matrix is regular.
+
+The order is chosen such that no nested Schur complements are necessary.
+However, in some common examples like the resistor network
+and Lagrange multiplicators we have a zero bottom right sub-matrix
+and the top left matrix is regular.
+-}
+instance
+   (typ00 ~ TypeFull, Divide.Solve typ11, Divide.Determinant typ11) =>
+      Divide.Determinant (Square typ00 measOff vertOff horizOff typ11) where
+   type DeterminantExtra (Square typ00 measOff vertOff horizOff typ11) extra =
+         (SquareType extra ~ TypeFull,
+          Square01 extra ~ (), Square10 extra ~ (),
+          Matrix.DeterminantExtra typ00 (Square00 extra),
+          Matrix.DeterminantExtra typ11 (Square11 extra),
+          Divide.SolveExtra typ11 (Square11 extra))
+   determinant sq@(Square _a _b _c d) =
+      determinant d * determinant (schurComplement sq)
+
+{- |
+Requires that the bottom right sub-matrix is regular.
+-}
+instance
+   (typ00 ~ TypeFull, Divide.Solve typ11) =>
+      Divide.Solve (Square typ00 measOff vertOff horizOff typ11) where
+   type SolveExtra (Square typ00 measOff vertOff horizOff typ11) extra =
+         (SquareType extra ~ TypeFull,
+          Square01 extra ~ (), Square10 extra ~ (),
+          Divide.SolveExtra typ00 (Square00 extra),
+          Divide.SolveExtra typ11 (Square11 extra))
+   solveRight sq@(Square _a b c d) =
+      withVerticalSplit $ \(Above x0 x1) ->
+      let xComplement = x0 #-# Matrix.fromFull b #*# Divide.solveRight d x1
+          y = Divide.solveRight (schurComplement sq) xComplement
+      in Above y (Divide.solveRight d $ x1 #-# Matrix.fromFull c #*# y)
+   solveLeft x sq@(Square _a b c d) =
+      withHorizontalSplit x $ \(Beside x0 x1) ->
+      let xComplement = x0 #-# Divide.solveLeft x1 d #*# Matrix.fromFull c
+          y = Divide.solveLeft xComplement (schurComplement sq)
+      in Beside y (Divide.solveLeft (x1 #-# y #*# Matrix.fromFull b) d)
+
+instance
+   (typ00 ~ TypeFull, typ11 ~ TypeFull) =>
+      Divide.Inverse (Square typ00 measOff vertOff horizOff typ11) where
+   type InverseExtra (Square typ00 measOff vertOff horizOff typ11) extra =
+         (SquareType extra ~ TypeFull,
+          Square01 extra ~ (), Square10 extra ~ (),
+          Divide.SolveExtra typ00 (Square00 extra),
+          Divide.SolveExtra typ11 (Square11 extra))
+   inverse sq@(Square a b c d) =
+      if False
+         then
+            let s = schurComplement (Square d c b a)
+                cainv = Divide.solveLeft c a
+                ainvb = Divide.solveRight a b
+                br = Divide.solveLeft ainvb s
+                cr = Divide.solveRight s cainv
+            in Square
+                  (Matrix.inverse a #+# br ##*## cainv)
+                  (Matrix.negate br)
+                  (Matrix.negate cr)
+                  (Matrix.inverse s)
+
+         else
+            let s = schurComplement sq
+                bdinv = Divide.solveLeft b d
+                dinvc = Divide.solveRight d c
+                cr = Divide.solveLeft dinvc s
+                br = Divide.solveRight s bdinv
+            in Square
+                  (Matrix.inverse s)
+                  (Matrix.negate br)
+                  (Matrix.negate cr)
+                  (Matrix.inverse d #+# cr ##*## bdinv)
+
+
+data Append typ0 typ1 sh0 sh1
+data instance
+   Matrix (Append typ0 typ1 sh0 sh1) xl xu
+      lower upper meas vert horiz height width a where
+   Above ::
+      (Shape.C sh0, Shape.C sh1, Eq width) =>
+      Matrix typ0 xl0 xu0 Filled Filled Size Big horiz sh0 width a ->
+      Matrix typ1 xl1 xu1 Filled Filled Size Big horiz sh1 width a ->
+      Matrix
+         (Append typ0 typ1 sh0 sh1)
+         (xl0,xl1,TBool.False) (xu0,xu1,TBool.True)
+         Filled Filled Size Big horiz (sh0::+sh1) width a
+   Beside ::
+      (Shape.C sh0, Shape.C sh1, Eq height) =>
+      Matrix typ0 xl0 xu0 Filled Filled Size vert Big height sh0 a ->
+      Matrix typ1 xl1 xu1 Filled Filled Size vert Big height sh1 a ->
+      Matrix
+         (Append typ0 typ1 sh0 sh1)
+         (xl0,xl1,TBool.True) (xu0,xu1,TBool.False)
+         Filled Filled Size vert Big height (sh0::+sh1) a
+
+type family Append0 extra; type instance Append0 (x0,x1,b) = x0
+type family Append1 extra; type instance Append1 (x0,x1,b) = x1
+
+type family AppendSelectShape xl sh part
+type instance AppendSelectShape (xl0,xl1,TBool.False) sh part = part
+type instance AppendSelectShape (xl0,xl1,TBool.True) sh part = sh
+
+type AboveGeneral height0 height1 width =
+   Matrix
+      (Append TypeFull TypeFull height0 height1)
+      ((),(),TBool.False) ((),(),TBool.True)
+      Filled Filled Size Big Big (height0::+height1) width
+
+type BesideGeneral height width0 width1 =
+   Matrix
+      (Append TypeFull TypeFull width0 width1)
+      ((),(),TBool.True) ((),(),TBool.False)
+      Filled Filled Size Big Big height (width0::+width1)
+
+aboveFromFull ::
+   (Shape.C height0, Eq height0, Shape.C height1, Eq height1,
+    Shape.C width, Eq width, Class.Floating a) =>
+   Matrix.General (height0 ::+ height1) width a ->
+   AboveGeneral height0 height1 width a
+aboveFromFull a0 =
+   let a = Matrix.toFull a0
+   in Above (Matrix.takeTop a) (Matrix.takeBottom a)
+
+besideFromFull ::
+   (Shape.C height, Eq height,
+    Shape.C width0, Eq width0, Shape.C width1, Eq width1, Class.Floating a) =>
+   Matrix.General height (width0 ::+ width1) a ->
+   BesideGeneral height width0 width1 a
+besideFromFull a0 =
+   let a = Matrix.toFull a0
+   in Beside (Matrix.takeLeft a) (Matrix.takeRight a)
+
+deriving instance
+   (Show (Matrix typ0 (Append0 xl) (Append0 xu)
+            lower upper meas vert horiz
+            (AppendSelectShape xl height sh0)
+            (AppendSelectShape xu width sh0) a),
+    Show (Matrix typ1 (Append1 xl) (Append1 xu)
+            lower upper meas vert horiz
+            (AppendSelectShape xl height sh1)
+            (AppendSelectShape xu width sh1) a),
+    Shape.C height, Shape.C width, Show height, Show width, Show a) =>
+   Show (Matrix (Append typ0 typ1 sh0 sh1) xl xu
+            lower upper meas vert horiz height width a)
+
+instance
+   (Matrix.Box typ0, Matrix.Box typ1) =>
+      Matrix.Box (Append typ0 typ1 sh0 sh1) where
+   type BoxExtra (Append typ0 typ1 sh0 sh1) extra =
+         (Matrix.BoxExtra typ0 (Append0 extra),
+          Matrix.BoxExtra typ1 (Append1 extra))
+   extent (Above a b) =
+      case Extent.appendSame of
+         Extent.AppendMode append ->
+            Extent.transpose $
+            append
+               (Extent.transpose $ Matrix.extent a)
+               (Extent.transpose $ Matrix.extent b)
+   extent (Beside a b) =
+      case Extent.appendSame of
+         Extent.AppendMode append ->
+            append (Matrix.extent a) (Matrix.extent b)
+
+instance
+   (Matrix.Transpose typ0, Matrix.Transpose typ1) =>
+      Matrix.Transpose (Append typ0 typ1 sh0 sh1) where
+   type TransposeExtra (Append typ0 typ1 sh0 sh1) extra =
+         (Matrix.TransposeExtra typ0 (Append0 extra),
+          Matrix.TransposeExtra typ1 (Append1 extra))
+   transpose (Above a b) = Beside (Matrix.transpose a) (Matrix.transpose b)
+   transpose (Beside a b) = Above (Matrix.transpose a) (Matrix.transpose b)
+
+instance
+   (Layout typ0, Layout typ1) =>
+      Layout (Append typ0 typ1 sh0 sh1) where
+   type LayoutExtra (Append typ0 typ1 sh0 sh1) extra =
+         (Matrix.BoxExtra typ0 (Append0 extra),
+          Matrix.BoxExtra typ1 (Append1 extra),
+          LayoutExtra typ0 (Append0 extra),
+          LayoutExtra typ1 (Append1 extra))
+   layout (Above a b) = finalizeLayout $ toRows a ===# toRows b
+   layout (Beside a b) = finalizeLayout $ toRows a |||# toRows b
+
+instance
+   (Layout typ0, Layout typ1) =>
+      Matrix.Format (Append typ0 typ1 sh0 sh1) where
+   type FormatExtra (Append typ0 typ1 sh0 sh1) extra =
+         (Matrix.BoxExtra typ0 (Append0 extra),
+          Matrix.BoxExtra typ1 (Append1 extra),
+          LayoutExtra typ0 (Append0 extra),
+          LayoutExtra typ1 (Append1 extra))
+   format = Matrix.formatWithLayout
+
+instance
+   (Matrix.Unpack typ0, Matrix.Unpack typ1) =>
+      Matrix.Unpack (Append typ0 typ1 sh0 sh1) where
+   type UnpackExtra (Append typ0 typ1 sh0 sh1) extra =
+         (Matrix.UnpackExtra typ0 (Append0 extra),
+          Matrix.UnpackExtra typ1 (Append1 extra))
+   unpack (Above a b) =
+      Matrix.above Matrix.rightBias Extent.appendSame
+         (MatrixClass.toFull a)
+         (MatrixClass.toFull b)
+   unpack (Beside a b) =
+      Matrix.beside Matrix.rightBias Extent.appendSame
+         (MatrixClass.toFull a)
+         (MatrixClass.toFull b)
+
+instance
+   (Matrix.Homogeneous typ0, Matrix.Homogeneous typ1) =>
+      Matrix.Homogeneous (Append typ0 typ1 sh0 sh1) where
+   type HomogeneousExtra (Append typ0 typ1 sh0 sh1) extra =
+         (Matrix.HomogeneousExtra typ0 (Append0 extra),
+          Matrix.HomogeneousExtra typ1 (Append1 extra))
+   zeroFrom (Above a b) = Above (Matrix.zeroFrom a) (Matrix.zeroFrom b)
+   zeroFrom (Beside a b) = Beside (Matrix.zeroFrom a) (Matrix.zeroFrom b)
+   negate (Above a b) = Above (Matrix.negate a) (Matrix.negate b)
+   negate (Beside a b) = Beside (Matrix.negate a) (Matrix.negate b)
+   scaleReal k (Above a b) =
+      Above (Matrix.scaleReal k a) (Matrix.scaleReal k b)
+   scaleReal k (Beside a b) =
+      Beside (Matrix.scaleReal k a) (Matrix.scaleReal k b)
+
+instance
+   (Matrix.Scale typ0, Matrix.Scale typ1) =>
+      Matrix.Scale (Append typ0 typ1 sh0 sh1) where
+   type ScaleExtra (Append typ0 typ1 sh0 sh1) extra =
+         (Matrix.ScaleExtra typ0 (Append0 extra),
+          Matrix.ScaleExtra typ1 (Append1 extra))
+   scale k (Above a b) = Above (Matrix.scale k a) (Matrix.scale k b)
+   scale k (Beside a b) = Beside (Matrix.scale k a) (Matrix.scale k b)
+
+instance
+   (Matrix.Additive typ0, Matrix.Additive typ1, Eq sh0, Eq sh1) =>
+      Matrix.Additive (Append typ0 typ1 sh0 sh1) where
+   type AdditiveExtra (Append typ0 typ1 sh0 sh1) extra =
+         (Matrix.AdditiveExtra typ0 (Append0 extra),
+          Matrix.AdditiveExtra typ1 (Append1 extra))
+   add (Above a0 b0) (Above a1 b1) =
+      Above (Matrix.add a0 a1) (Matrix.add b0 b1)
+   add (Beside a0 b0) (Beside a1 b1) =
+      Beside (Matrix.add a0 a1) (Matrix.add b0 b1)
+
+instance
+   (Matrix.Subtractive typ0, Matrix.Subtractive typ1, Eq sh0, Eq sh1) =>
+      Matrix.Subtractive (Append typ0 typ1 sh0 sh1) where
+   type SubtractiveExtra (Append typ0 typ1 sh0 sh1) extra =
+         (Matrix.SubtractiveExtra typ0 (Append0 extra),
+          Matrix.SubtractiveExtra typ1 (Append1 extra))
+   sub (Above a0 b0) (Above a1 b1) =
+      Above (Matrix.sub a0 a1) (Matrix.sub b0 b1)
+   sub (Beside a0 b0) (Beside a1 b1) =
+      Beside (Matrix.sub a0 a1) (Matrix.sub b0 b1)
+
+instance
+   (Multiply.MultiplyVector typ0, Multiply.MultiplyVector typ1,
+    Eq sh0, Eq sh1) =>
+      Multiply.MultiplyVector (Append typ0 typ1 sh0 sh1) where
+   type MultiplyVectorExtra (Append typ0 typ1 sh0 sh1) extra =
+         (Multiply.MultiplyVectorExtra typ0 (Append0 extra),
+          Multiply.MultiplyVectorExtra typ1 (Append1 extra))
+   matrixVector (Above a b) x =
+      Array.append (Multiply.matrixVector a x) (Multiply.matrixVector b x)
+   matrixVector (Beside a b) x =
+      let (x0,x1) = Array.split x
+      in withoutSizeCheck (|+|)
+            (Multiply.matrixVector a x0)
+            (Multiply.matrixVector b x1)
+   vectorMatrix x (Above a b) =
+      let (x0,x1) = Array.split x
+      in withoutSizeCheck (|+|)
+            (Multiply.vectorMatrix x0 a)
+            (Multiply.vectorMatrix x1 b)
+   vectorMatrix x (Beside a b) =
+      Array.append (Multiply.vectorMatrix x a) (Multiply.vectorMatrix x b)
+
+
+data Triangular typ0 typOff typ1
+data instance
+   Matrix (Triangular typ0 typOff typ1) xl xu
+      lower upper meas vert horiz height width a where
+   Upper ::
+      (Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1) =>
+      Quadratic typ0 xl0 xu0 lower Filled sh0 a ->
+      Matrix typOff xlOff xuOff Filled Filled Size Big Big sh0 sh1 a ->
+      Quadratic typ1 xl1 xu1 lower Filled sh1 a ->
+      Quadratic
+         (Triangular typ0 typOff typ1)
+         (xl0,xlOff,xl1,TBool.False) (xu0,xuOff,xu1,TBool.True)
+         lower Filled (sh0::+sh1) a
+   Lower ::
+      (Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1) =>
+      Quadratic typ0 xl0 xu0 Filled upper sh0 a ->
+      Matrix typOff xlOff xuOff Filled Filled Size Big Big sh1 sh0 a ->
+      Quadratic typ1 xl1 xu1 Filled upper sh1 a ->
+      Quadratic
+         (Triangular typ0 typOff typ1)
+         (xl0,xlOff,xl1,TBool.True) (xu0,xuOff,xu1,TBool.False)
+         Filled upper (sh0::+sh1) a
+
+type family Triangular0 extra; type instance Triangular0 (x0,xo,x1,b) = x0
+type family Triangular1 extra; type instance Triangular1 (x0,xo,x1,b) = x1
+type family TriangularOff extra; type instance TriangularOff (x0,xo,x1,b) = xo
+
+type family TriangularFstShape xl xu sh
+type family TriangularSndShape xl xu sh
+type instance
+   TriangularFstShape
+      (xl0,xlo,xl1,TBool.False) (xu0,xuo,xu1,TBool.True) (sh0::+sh1) = sh0
+type instance
+   TriangularSndShape
+      (xl0,xlo,xl1,TBool.False) (xu0,xuo,xu1,TBool.True) (sh0::+sh1) = sh1
+type instance
+   TriangularFstShape
+      (xl0,xlo,xl1,TBool.True) (xu0,xuo,xu1,TBool.False) (sh0::+sh1) = sh1
+type instance
+   TriangularSndShape
+      (xl0,xlo,xl1,TBool.True) (xu0,xuo,xu1,TBool.False) (sh0::+sh1) = sh0
+
+deriving instance
+   (Show (Quadratic typ0 (Triangular0 xl) (Triangular0 xu)
+            lower upper (ShapeHead height) a),
+    Show (Quadratic typ1 (Triangular1 xl) (Triangular1 xu)
+            lower upper (ShapeTail height) a),
+    Show (Matrix typOff (TriangularOff xl) (TriangularOff xu)
+            Filled Filled Size Big Big
+            (TriangularFstShape xl xu height)
+            (TriangularSndShape xl xu height)
+            a),
+    Shape.C height, Shape.C width, Show height, Show width, Show a) =>
+   Show (Matrix (Triangular typ0 typOff typ1) xl xu
+            lower upper meas vert horiz height width a)
+
+instance
+   (Matrix.Box typ0, Matrix.Box typOff, Matrix.Box typ1) =>
+      Matrix.Box (Triangular typ0 typOff typ1) where
+   type BoxExtra (Triangular typ0 typOff typ1) extra =
+         (Matrix.BoxExtra typ0 (Triangular0 extra),
+          Matrix.BoxExtra typ1 (Triangular1 extra),
+          Matrix.BoxExtra typOff (TriangularOff extra))
+   extent (Upper a _ b) = Extent.square (squareSize a ::+ squareSize b)
+   extent (Lower a _ b) = Extent.square (squareSize a ::+ squareSize b)
+
+instance
+   (Matrix.Box typ0, Matrix.Box typOff, Matrix.Box typ1) =>
+      Matrix.ToQuadratic (Triangular typ0 typOff typ1) where
+   heightToQuadratic a@(Upper _ _ _) = a
+   heightToQuadratic a@(Lower _ _ _) = a
+   widthToQuadratic a@(Upper _ _ _) = a
+   widthToQuadratic a@(Lower _ _ _) = a
+
+instance
+   (Matrix.Transpose typ0, Matrix.Transpose typOff, Matrix.Transpose typ1) =>
+      Matrix.Transpose (Triangular typ0 typOff typ1) where
+   type TransposeExtra (Triangular typ0 typOff typ1) extra =
+         (Matrix.TransposeExtra typ0 (Triangular0 extra),
+          Matrix.TransposeExtra typ1 (Triangular1 extra),
+          Matrix.TransposeExtra typOff (TriangularOff extra))
+   transpose (Upper a o b) =
+      Lower (Matrix.transpose a) (Matrix.transpose o) (Matrix.transpose b)
+   transpose (Lower a o b) =
+      Upper (Matrix.transpose a) (Matrix.transpose o) (Matrix.transpose b)
+
+instance
+   (Layout typ0, Layout typOff, Layout typ1) =>
+      Layout (Triangular typ0 typOff typ1) where
+   type LayoutExtra (Triangular typ0 typOff typ1) extra =
+         (Matrix.BoxExtra typ0 (Triangular0 extra),
+          Matrix.BoxExtra typ1 (Triangular1 extra),
+          Matrix.BoxExtra typOff (TriangularOff extra),
+          LayoutExtra typ0 (Triangular0 extra),
+          LayoutExtra typ1 (Triangular1 extra),
+          LayoutExtra typOff (TriangularOff extra))
+   layout (Upper a o b) = finalizeLayout $
+      toRows a |||# toRows o
+      ===#
+      layoutEmpty (squareSize b) (squareSize a) |||# toRows b
+   layout (Lower a o b) = finalizeLayout $
+      toRows a |||# layoutEmpty (squareSize a) (squareSize b)
+      ===#
+      toRows o |||# toRows b
+
+instance
+   (Layout typ0, Layout typOff, Layout typ1) =>
+      Matrix.Format (Triangular typ0 typOff typ1) where
+   type FormatExtra (Triangular typ0 typOff typ1) extra =
+         (Matrix.BoxExtra typ0 (Triangular0 extra),
+          Matrix.BoxExtra typ1 (Triangular1 extra),
+          Matrix.BoxExtra typOff (TriangularOff extra),
+          LayoutExtra typ0 (Triangular0 extra),
+          LayoutExtra typ1 (Triangular1 extra),
+          LayoutExtra typOff (TriangularOff extra))
+   format = Matrix.formatWithLayout
+
+instance
+   (Matrix.SquareShape typ0, Matrix.SquareShape typ1,
+    Matrix.Box typOff, Matrix.Homogeneous typOff) =>
+      Matrix.SquareShape (Triangular typ0 typOff typ1) where
+   type SquareShapeExtra (Triangular typ0 typOff typ1) extra =
+         (Matrix.SquareShapeExtra typ0 (Triangular0 extra),
+          Matrix.SquareShapeExtra typ1 (Triangular1 extra),
+          Matrix.SquareShapeExtra typOff (TriangularOff extra),
+          Matrix.HomogeneousExtra typOff (TriangularOff extra))
+   takeDiagonal (Upper a _b c) =
+      Array.append (MatrixClass.takeDiagonal a) (MatrixClass.takeDiagonal c)
+   takeDiagonal (Lower a _b c) =
+      Array.append (MatrixClass.takeDiagonal a) (MatrixClass.takeDiagonal c)
+   identityFrom (Upper a b c) =
+      Upper
+         (MatrixClass.identityFrom a)
+         (MatrixClass.zeroFrom b)
+         (MatrixClass.identityFrom c)
+   identityFrom (Lower a b c) =
+      Lower
+         (MatrixClass.identityFrom a)
+         (MatrixClass.zeroFrom b)
+         (MatrixClass.identityFrom c)
+
+instance
+   (Matrix.Unpack typ0, Matrix.Unpack typOff, Matrix.Unpack typ1) =>
+      Matrix.Unpack (Triangular typ0 typOff typ1) where
+   type UnpackExtra (Triangular typ0 typOff typ1) extra =
+         (Matrix.UnpackExtra typ0 (Triangular0 extra),
+          Matrix.UnpackExtra typ1 (Triangular1 extra),
+          Matrix.UnpackExtra typOff (TriangularOff extra))
+   unpack (Upper a0 o0 b0) =
+      let a = MatrixClass.toFull a0 in
+      let b = MatrixClass.toFull b0 in
+      let o = MatrixClass.toFull o0 in
+      let zero shape = Matrix.asGeneral $ Matrix.zero shape in
+      let order = ArrMatrix.order b in
+      let dims = (squareSize b, squareSize a) in
+      ArrMatrix.liftUnpacked1 id $
+      Square.stack
+         a                              o
+         (zero $ Omni.cons order dims)  b
+   unpack (Lower a0 o0 b0) =
+      let a = MatrixClass.toFull a0 in
+      let b = MatrixClass.toFull b0 in
+      let o = MatrixClass.toFull o0 in
+      let zero shape = Matrix.asGeneral $ Matrix.zero shape in
+      let order = ArrMatrix.order b in
+      let dims = (squareSize a, squareSize b) in
+      ArrMatrix.liftUnpacked1 id $
+      Square.stack
+         a  (zero $ Omni.cons order dims)
+         o  b
+
+instance
+   (Matrix.Homogeneous typ0, Matrix.Homogeneous typOff,
+    Matrix.Homogeneous typ1) =>
+      Matrix.Homogeneous (Triangular typ0 typOff typ1) where
+   type HomogeneousExtra (Triangular typ0 typOff typ1) extra =
+         (Matrix.HomogeneousExtra typ0 (Triangular0 extra),
+          Matrix.HomogeneousExtra typ1 (Triangular1 extra),
+          Matrix.HomogeneousExtra typOff (TriangularOff extra))
+   zeroFrom (Upper a o b) =
+      Upper (Matrix.zeroFrom a) (Matrix.zeroFrom o) (Matrix.zeroFrom b)
+   zeroFrom (Lower a o b) =
+      Lower (Matrix.zeroFrom a) (Matrix.zeroFrom o) (Matrix.zeroFrom b)
+   negate (Upper a o b) =
+      Upper (Matrix.negate a) (Matrix.negate o) (Matrix.negate b)
+   negate (Lower a o b) =
+      Lower (Matrix.negate a) (Matrix.negate o) (Matrix.negate b)
+   scaleReal k (Upper a o b) =
+      Upper (Matrix.scaleReal k a) (Matrix.scaleReal k o) (Matrix.scaleReal k b)
+   scaleReal k (Lower a o b) =
+      Lower (Matrix.scaleReal k a) (Matrix.scaleReal k o) (Matrix.scaleReal k b)
+
+instance
+   (Matrix.Scale typ0, Matrix.Scale typOff, Matrix.Scale typ1) =>
+      Matrix.Scale (Triangular typ0 typOff typ1) where
+   type ScaleExtra (Triangular typ0 typOff typ1) extra =
+         (Matrix.ScaleExtra typ0 (Triangular0 extra),
+          Matrix.ScaleExtra typ1 (Triangular1 extra),
+          Matrix.ScaleExtra typOff (TriangularOff extra))
+   scale k (Upper a o b) =
+      Upper (Matrix.scale k a) (Matrix.scale k o) (Matrix.scale k b)
+   scale k (Lower a o b) =
+      Lower (Matrix.scale k a) (Matrix.scale k o) (Matrix.scale k b)
+
+instance
+   (Matrix.Additive typ0, Matrix.Additive typOff, Matrix.Additive typ1) =>
+      Matrix.Additive (Triangular typ0 typOff typ1) where
+   type AdditiveExtra (Triangular typ0 typOff typ1) extra =
+         (Matrix.AdditiveExtra typ0 (Triangular0 extra),
+          Matrix.AdditiveExtra typ1 (Triangular1 extra),
+          Matrix.AdditiveExtra typOff (TriangularOff extra))
+   add (Upper a0 o0 b0) (Upper a1 o1 b1) =
+      Upper (Matrix.add a0 a1) (Matrix.add o0 o1) (Matrix.add b0 b1)
+   add (Lower a0 o0 b0) (Lower a1 o1 b1) =
+      Lower (Matrix.add a0 a1) (Matrix.add o0 o1) (Matrix.add b0 b1)
+
+instance
+   (Matrix.Subtractive typ0, Matrix.Subtractive typOff,
+    Matrix.Subtractive typ1) =>
+      Matrix.Subtractive (Triangular typ0 typOff typ1) where
+   type SubtractiveExtra (Triangular typ0 typOff typ1) extra =
+         (Matrix.SubtractiveExtra typ0 (Triangular0 extra),
+          Matrix.SubtractiveExtra typ1 (Triangular1 extra),
+          Matrix.SubtractiveExtra typOff (TriangularOff extra))
+   sub (Upper a0 o0 b0) (Upper a1 o1 b1) =
+      Upper (Matrix.sub a0 a1) (Matrix.sub o0 o1) (Matrix.sub b0 b1)
+   sub (Lower a0 o0 b0) (Lower a1 o1 b1) =
+      Lower (Matrix.sub a0 a1) (Matrix.sub o0 o1) (Matrix.sub b0 b1)
+
+instance
+   (Multiply.MultiplyVector typ0, Multiply.MultiplyVector typ1,
+    Multiply.MultiplyVector typOff) =>
+      Multiply.MultiplyVector (Triangular typ0 typOff typ1) where
+   type MultiplyVectorExtra (Triangular typ0 typOff typ1) extra =
+         (Multiply.MultiplyVectorExtra typ0 (Triangular0 extra),
+          Multiply.MultiplyVectorExtra typ1 (Triangular1 extra),
+          Multiply.MultiplyVectorExtra typOff (TriangularOff extra))
+   matrixVector (Upper a o b) x =
+      let (x0,x1) = Array.split x
+      in Array.append
+            (Multiply.matrixVector a x0 |+| Multiply.matrixVector o x1)
+            (Multiply.matrixVector b x1)
+   matrixVector (Lower a o b) x =
+      let (x0,x1) = Array.split x
+      in Array.append
+            (Multiply.matrixVector a x0)
+            (Multiply.matrixVector o x0 |+| Multiply.matrixVector b x1)
+   vectorMatrix x (Upper a o b) =
+      let (x0,x1) = Array.split x
+      in Array.append
+            (Multiply.vectorMatrix x0 a)
+            (Multiply.vectorMatrix x0 o |+| Multiply.vectorMatrix x1 b)
+   vectorMatrix x (Lower a o b) =
+      let (x0,x1) = Array.split x
+      in Array.append
+            (Multiply.vectorMatrix x0 a |+| Multiply.vectorMatrix x1 o)
+            (Multiply.vectorMatrix x1 b)
+
+instance
+   (Multiply.MultiplySquare typ0, Multiply.MultiplySquare typ1,
+    typOff ~ TypeFull) =>
+      Multiply.MultiplySquare (Triangular typ0 typOff typ1) where
+   type MultiplySquareExtra (Triangular typ0 typOff typ1) extra =
+         (Matrix.MultiplySquareExtra typ0 (Triangular0 extra),
+          Matrix.MultiplySquareExtra typ1 (Triangular1 extra),
+          TriangularOff extra ~ ())
+
+   squareFull (Upper a o b) =
+      withVerticalSplit $ \(Above top bottom) ->
+         Above
+            (Multiply.squareFull a top #+# o#*#bottom)
+            (Multiply.squareFull b bottom)
+   squareFull (Lower a o b) =
+      withVerticalSplit $ \(Above top bottom) ->
+         Above
+            (Multiply.squareFull a top)
+            (o#*#top #+# Multiply.squareFull b bottom)
+   fullSquare x (Lower a o b) =
+      withHorizontalSplit x $ \(Beside left right) ->
+         Beside
+            (Multiply.fullSquare left a #+# right#*#o)
+            (Multiply.fullSquare right b)
+   fullSquare x (Upper a o b) =
+      withHorizontalSplit x $ \(Beside left right) ->
+         Beside
+            (Multiply.fullSquare left a)
+            (left#*#o #+# Multiply.fullSquare right b)
+
+instance
+   (typ0 ~ TypeFull, typOff ~ TypeFull, typ1 ~ TypeFull) =>
+      Matrix.MultiplySame (Triangular typ0 typOff typ1) where
+   type MultiplySameExtra (Triangular typ0 typOff typ1) extra =
+         (Triangular0 extra ~ (), Triangular1 extra ~ (),
+          TriangularOff extra ~ ())
+   multiplySame (Upper a0 b0 c0) (Upper a1 b1 c1) =
+      case Omni.powerStrips $ ArrMatrix.shape a0 of
+         (Omni.Empty,  _) -> Upper (a0#*#a1) (a0#*##b1 #+# b0##*#c1) (c0#*#c1)
+         (Omni.Filled, _) -> Upper (a0#*#a1) (a0#*##b1 #+# b0##*#c1) (c0#*#c1)
+   multiplySame (Lower a0 b0 c0) (Lower a1 b1 c1) =
+      case Omni.powerStrips $ ArrMatrix.shape a0 of
+         (_, Omni.Empty ) -> Lower (a0#*#a1) (b0##*#a1 #+# c0#*##b1) (c0#*#c1)
+         (_, Omni.Filled) -> Lower (a0#*#a1) (b0##*#a1 #+# c0#*##b1) (c0#*#c1)
+
+instance
+   (typ0 ~ TypeFull, typOff ~ TypeFull, typ1 ~ TypeFull) =>
+      Multiply.Power (Triangular typ0 typOff typ1) where
+   type PowerExtra (Triangular typ0 typOff typ1) extra =
+         (Triangular0 extra ~ (), Triangular1 extra ~ (),
+          TriangularOff extra ~ ())
+   square a@(Upper _ _ _) = Matrix.multiplySame a a
+   square a@(Lower _ _ _) = Matrix.multiplySame a a
+   power n m@(Upper a b c) =
+      powerAssociative Matrix.multiplySame
+         (Upper (Matrix.identityFrom a)
+            (Matrix.zeroFrom b) (Matrix.identityFrom c))
+         m n
+   power n m@(Lower a b c) =
+      powerAssociative Matrix.multiplySame
+         (Lower (Matrix.identityFrom a)
+            (Matrix.zeroFrom b) (Matrix.identityFrom c))
+         m n
+   powers1 a@(Upper _ _ _) = Stream.iterate (flip Matrix.multiplySame a) a
+   powers1 a@(Lower _ _ _) = Stream.iterate (flip Matrix.multiplySame a) a
+
+instance
+   (Divide.Determinant typ0, Matrix.Box typOff, Divide.Determinant typ1) =>
+      Divide.Determinant (Triangular typ0 typOff typ1) where
+   type DeterminantExtra (Triangular typ0 typOff typ1) extra =
+         (Matrix.DeterminantExtra typ0 (Triangular0 extra),
+          Matrix.DeterminantExtra typ1 (Triangular1 extra),
+          Matrix.BoxExtra typOff (TriangularOff extra))
+   determinant (Upper a _ b) = determinant a * determinant b
+   determinant (Lower a _ b) = determinant a * determinant b
+
+instance
+   (Divide.Solve typ0, Divide.Solve typ1, typOff ~ TypeFull) =>
+      Divide.Solve (Triangular typ0 typOff typ1) where
+   type SolveExtra (Triangular typ0 typOff typ1) extra =
+         (Matrix.SolveExtra typ0 (Triangular0 extra),
+          Matrix.SolveExtra typ1 (Triangular1 extra),
+          TriangularOff extra ~ ())
+
+   solveRight (Upper a o b) =
+      withVerticalSplit $ \(Above rhsTop rhsBottom) ->
+      let xBottom = Divide.solveRight b rhsBottom
+      in Above (Divide.solveRight a $ rhsTop #-# o #*# xBottom) xBottom
+   solveRight (Lower a o b) =
+      withVerticalSplit $ \(Above rhsTop rhsBottom) ->
+      let xTop = Divide.solveRight a rhsTop
+      in Above xTop (Divide.solveRight b $ rhsBottom #-# o #*# xTop)
+
+   solveLeft rhs (Lower a o b) =
+      withHorizontalSplit rhs $ \(Beside rhsLeft rhsRight) ->
+      let xRight = Divide.solveLeft rhsRight b
+      in Beside (Divide.solveLeft (rhsLeft #-# xRight #*# o) a) xRight
+   solveLeft rhs (Upper a o b) =
+      withHorizontalSplit rhs $ \(Beside rhsLeft rhsRight) ->
+      let xLeft = Divide.solveLeft rhsLeft a
+      in Beside xLeft (Divide.solveLeft (rhsRight #-# xLeft #*# o) b)
+
+instance
+   (Divide.Inverse typ0, Divide.Inverse typ1, typOff ~ TypeFull) =>
+      Divide.Inverse (Triangular typ0 typOff typ1) where
+   type InverseExtra (Triangular typ0 typOff typ1) extra =
+         (Matrix.InverseExtra typ0 (Triangular0 extra),
+          Matrix.InverseExtra typ1 (Triangular1 extra),
+          Matrix.SolveExtra typ0 (Triangular0 extra),
+          Matrix.SolveExtra typ1 (Triangular1 extra),
+          TriangularOff extra ~ ())
+   inverse (Upper a o b) =
+      Upper
+         (Divide.inverse a)
+         (Matrix.negate $ Divide.solveRight a $ Divide.solveLeft o b)
+         (Divide.inverse b)
+   inverse (Lower a o b) =
+      Lower
+         (Divide.inverse a)
+         (Matrix.negate $ Divide.solveRight b $ Divide.solveLeft o a)
+         (Divide.inverse b)
+
+
+{- |
+It is not necessary that the square blocks on the diagonal are symmetric.
+But if they are, according specialised algorithms are used.
+
+We have no algorithms that benefit from this structure.
+Using this data type merely documents the structure
+and saves a field in the record
+(which comes in handy for nested block matrices).
+-}
+data Symmetric typ0 typOff xlOff xuOff typ1
+data instance
+   Matrix (Symmetric typ0 typOff xlOff xuOff typ1) xl xu
+      lower upper meas vert horiz height width a where
+   Symmetric ::
+      (Shape.C sh0, Eq sh0, Shape.C sh1, Eq sh1) =>
+      Quadratic typ0 xl0 xu0 Filled Filled sh0 a ->
+      Matrix typOff xlOff xuOff Filled Filled Size Big Big sh0 sh1 a ->
+      Quadratic typ1 xl1 xu1 Filled Filled sh1 a ->
+      Quadratic
+         (Symmetric typ0 typOff xlOff xuOff typ1)
+         (xl0,xl1) (xu0,xu1)
+         Filled Filled (sh0::+sh1) a
+
+squareFromSymmetric ::
+   (Matrix.Transpose typOff,
+    Matrix.TransposeExtra typOff xlOff,
+    Matrix.TransposeExtra typOff xuOff,
+    Class.Floating a) =>
+   Quadratic
+      (Symmetric typ0 typOff xlOff xuOff typ1)
+      (xl0,xl1) (xu0,xu1)
+      lower upper sh a ->
+   Quadratic
+      (Square typ0 Size Big Big typ1)
+      (xl0,xl1,(typOff,xuOff,xlOff))
+      (xu0,xu1,(typOff,xuOff,xlOff))
+      lower upper sh a
+squareFromSymmetric (Symmetric a o b) = Square a o (Matrix.transpose o) b
+
+type family Symmetric0 extra; type instance Symmetric0 (x0,x1) = x0
+type family Symmetric1 extra; type instance Symmetric1 (x0,x1) = x1
+
+deriving instance
+   (Show (Quadratic typ0 (Symmetric0 xl) (Symmetric0 xu)
+            lower upper (ShapeHead height) a),
+    Show (Quadratic typ1 (Symmetric1 xl) (Symmetric1 xu)
+            lower upper (ShapeTail height) a),
+    Show (Matrix typOff xlOff xuOff Filled Filled Size Big Big
+            (ShapeHead height) (ShapeTail height) a),
+    Shape.C height, Shape.C width, Show height, Show width, Show a) =>
+   Show (Matrix (Symmetric typ0 typOff xlOff xuOff typ1) xl xu
+            lower upper meas vert horiz height width a)
+
+instance
+   (Matrix.Box typ0, Matrix.Box typ1, Matrix.Box typOff,
+    Matrix.BoxExtra typOff xlOff, Matrix.BoxExtra typOff xuOff) =>
+      Matrix.Box (Symmetric typ0 typOff xlOff xuOff typ1) where
+   type BoxExtra (Symmetric typ0 typOff xlOff xuOff typ1) extra =
+         (Matrix.BoxExtra typ0 (Symmetric0 extra),
+          Matrix.BoxExtra typ1 (Symmetric1 extra))
+   extent (Symmetric a _ b) = Extent.square (squareSize a ::+ squareSize b)
+
+instance
+   (Matrix.Transpose typ0, Matrix.Transpose typ1, Matrix.Transpose typOff,
+    Matrix.BoxExtra typOff xlOff, Matrix.BoxExtra typOff xuOff,
+    Matrix.TransposeExtra typOff xlOff, Matrix.TransposeExtra typOff xuOff) =>
+      Matrix.Transpose (Symmetric typ0 typOff xlOff xuOff typ1) where
+   type TransposeExtra (Symmetric typ0 typOff xlOff xuOff typ1) extra =
+         (Matrix.TransposeExtra typ0 (Symmetric0 extra),
+          Matrix.TransposeExtra typ1 (Symmetric1 extra))
+   transpose (Symmetric a o b) =
+      Symmetric (Matrix.transpose a) o (Matrix.transpose b)
+
+instance
+   (Matrix.Box typ0, Matrix.Box typOff, Matrix.Box typ1,
+    Matrix.BoxExtra typOff xlOff, Matrix.BoxExtra typOff xuOff) =>
+      Matrix.ToQuadratic (Symmetric typ0 typOff xlOff xuOff typ1) where
+   heightToQuadratic a@(Symmetric _ _ _) = a
+   widthToQuadratic a@(Symmetric _ _ _) = a
+
+instance
+   (Matrix.BoxExtra typOff xlOff, Matrix.BoxExtra typOff xuOff,
+    Layout typ0, Layout typ1,
+    Layout typOff, LayoutExtra typOff xlOff, LayoutExtra typOff xuOff) =>
+      Layout (Symmetric typ0 typOff xlOff xuOff typ1) where
+   type LayoutExtra (Symmetric typ0 typOff xlOff xuOff typ1) extra =
+         (Matrix.BoxExtra typ0 (Symmetric0 extra),
+          Matrix.BoxExtra typ1 (Symmetric1 extra),
+          LayoutExtra typ0 (Symmetric0 extra),
+          LayoutExtra typ1 (Symmetric1 extra))
+   layout (Symmetric a o b) = finalizeLayout $
+      toRows a |||# toRows o
+      ===#
+      toColumns o |||# toRows b
+
+instance
+   (Matrix.BoxExtra typOff xlOff, Matrix.BoxExtra typOff xuOff,
+    Layout typ0, Layout typ1,
+    Layout typOff, LayoutExtra typOff xlOff, LayoutExtra typOff xuOff) =>
+      Matrix.Format (Symmetric typ0 typOff xlOff xuOff typ1) where
+   type FormatExtra (Symmetric typ0 typOff xlOff xuOff typ1) extra =
+         (Matrix.BoxExtra typ0 (Symmetric0 extra),
+          Matrix.BoxExtra typ1 (Symmetric1 extra),
+          LayoutExtra typ0 (Symmetric0 extra),
+          LayoutExtra typ1 (Symmetric1 extra))
+   format = Matrix.formatWithLayout
+
+instance
+   (Matrix.SquareShape typ0, Matrix.SquareShape typ1,
+    Matrix.Box typOff, Matrix.Homogeneous typOff,
+    Matrix.BoxExtra typOff xlOff, Matrix.BoxExtra typOff xuOff,
+    MatrixClass.HomogeneousExtra typOff xlOff,
+    MatrixClass.HomogeneousExtra typOff xuOff) =>
+      Matrix.SquareShape (Symmetric typ0 typOff xlOff xuOff typ1) where
+   type SquareShapeExtra (Symmetric typ0 typOff xlOff xuOff typ1) extra =
+         (Matrix.SquareShapeExtra typ0 (Symmetric0 extra),
+          Matrix.SquareShapeExtra typ1 (Symmetric1 extra))
+   takeDiagonal (Symmetric a _b c) =
+      Array.append (MatrixClass.takeDiagonal a) (MatrixClass.takeDiagonal c)
+   identityFrom (Symmetric a b c) =
+      Symmetric
+         (MatrixClass.identityFrom a)
+         (MatrixClass.zeroFrom b)
+         (MatrixClass.identityFrom c)
+
+instance
+   (Matrix.Unpack typ0, Matrix.Unpack typOff, Matrix.Unpack typ1,
+    MatrixClass.UnpackExtra typOff xlOff,
+    MatrixClass.UnpackExtra typOff xuOff) =>
+      Matrix.Unpack (Symmetric typ0 typOff xlOff xuOff typ1) where
+   type UnpackExtra (Symmetric typ0 typOff xlOff xuOff typ1) extra =
+         (Matrix.UnpackExtra typ0 (Symmetric0 extra),
+          Matrix.UnpackExtra typ1 (Symmetric1 extra))
+   unpack (Symmetric a0 o0 b0) =
+      let o = MatrixClass.toFull o0 in
+      Square.stack
+         (MatrixClass.toFull a0)  o
+         (Matrix.transpose o)     (MatrixClass.toFull b0)
+
+instance
+   (Matrix.Homogeneous typ0, Matrix.Homogeneous typOff, Matrix.Homogeneous typ1,
+    MatrixClass.HomogeneousExtra typOff xlOff,
+    MatrixClass.HomogeneousExtra typOff xuOff) =>
+      Matrix.Homogeneous (Symmetric typ0 typOff xlOff xuOff typ1) where
+   type HomogeneousExtra (Symmetric typ0 typOff xlOff xuOff typ1) extra =
+         (Matrix.HomogeneousExtra typ0 (Symmetric0 extra),
+          Matrix.HomogeneousExtra typ1 (Symmetric1 extra))
+   zeroFrom (Symmetric a o b) =
+      Symmetric (Matrix.zeroFrom a) (Matrix.zeroFrom o) (Matrix.zeroFrom b)
+   negate (Symmetric a o b) =
+      Symmetric (Matrix.negate a) (Matrix.negate o) (Matrix.negate b)
+   scaleReal k (Symmetric a o b) =
+      Symmetric
+         (Matrix.scaleReal k a) (Matrix.scaleReal k o) (Matrix.scaleReal k b)
+
+instance
+   (Matrix.Scale typ0, Matrix.Scale typOff, Matrix.Scale typ1,
+    MatrixClass.ScaleExtra typOff xlOff,
+    MatrixClass.ScaleExtra typOff xuOff,
+    MatrixClass.HomogeneousExtra typOff xlOff,
+    MatrixClass.HomogeneousExtra typOff xuOff) =>
+      Matrix.Scale (Symmetric typ0 typOff xlOff xuOff typ1) where
+   type ScaleExtra (Symmetric typ0 typOff xlOff xuOff typ1) extra =
+         (Matrix.ScaleExtra typ0 (Symmetric0 extra),
+          Matrix.ScaleExtra typ1 (Symmetric1 extra))
+   scale k (Symmetric a o b) =
+      Symmetric (Matrix.scale k a) (Matrix.scale k o) (Matrix.scale k b)
+
+instance
+   (Matrix.Additive typ0, Matrix.Additive typOff, Matrix.Additive typ1,
+    MatrixClass.AdditiveExtra typOff xlOff,
+    MatrixClass.AdditiveExtra typOff xuOff) =>
+      Matrix.Additive (Symmetric typ0 typOff xlOff xuOff typ1) where
+   type AdditiveExtra (Symmetric typ0 typOff xlOff xuOff typ1) extra =
+         (Matrix.AdditiveExtra typ0 (Symmetric0 extra),
+          Matrix.AdditiveExtra typ1 (Symmetric1 extra))
+   add (Symmetric a0 o0 b0) (Symmetric a1 o1 b1) =
+      Symmetric (Matrix.add a0 a1) (Matrix.add o0 o1) (Matrix.add b0 b1)
+
+instance
+   (Matrix.Subtractive typ0, Matrix.Subtractive typOff, Matrix.Subtractive typ1,
+    MatrixClass.SubtractiveExtra typOff xlOff,
+    MatrixClass.SubtractiveExtra typOff xuOff,
+    MatrixClass.AdditiveExtra typOff xlOff,
+    MatrixClass.AdditiveExtra typOff xuOff) =>
+      Matrix.Subtractive (Symmetric typ0 typOff xlOff xuOff typ1) where
+   type SubtractiveExtra (Symmetric typ0 typOff xlOff xuOff typ1) extra =
+         (Matrix.SubtractiveExtra typ0 (Symmetric0 extra),
+          Matrix.SubtractiveExtra typ1 (Symmetric1 extra))
+   sub (Symmetric a0 o0 b0) (Symmetric a1 o1 b1) =
+      Symmetric (Matrix.sub a0 a1) (Matrix.sub o0 o1) (Matrix.sub b0 b1)
+
+instance
+   (Multiply.MultiplyVector typ0, Multiply.MultiplyVector typ1,
+    Matrix.BoxExtra typOff xlOff, Matrix.BoxExtra typOff xuOff,
+    Multiply.MultiplyVector typOff,
+    Multiply.MultiplyVectorExtra typOff xlOff,
+    Multiply.MultiplyVectorExtra typOff xuOff) =>
+      Multiply.MultiplyVector (Symmetric typ0 typOff xlOff xuOff typ1) where
+   type MultiplyVectorExtra (Symmetric typ0 typOff xlOff xuOff typ1) extra =
+         (Multiply.MultiplyVectorExtra typ0 (Symmetric0 extra),
+          Multiply.MultiplyVectorExtra typ1 (Symmetric1 extra))
+   matrixVector (Symmetric a o b) x =
+      let (x0,x1) = Array.split x
+      in Array.append
+            (Multiply.matrixVector a x0 |+| Multiply.matrixVector o x1)
+            (Multiply.vectorMatrix x0 o |+| Multiply.matrixVector b x1)
+   vectorMatrix x (Symmetric a o b) =
+      let (x0,x1) = Array.split x
+      in Array.append
+            (Multiply.vectorMatrix x0 a |+| Multiply.matrixVector o x1)
+            (Multiply.vectorMatrix x0 o |+| Multiply.vectorMatrix x1 b)
+
+instance
+   (Multiply.MultiplySquare typ0, Multiply.MultiplySquare typ1,
+    typOff ~ TypeFull, xlOff ~ (), xuOff ~ ()) =>
+      Multiply.MultiplySquare (Symmetric typ0 typOff xlOff xuOff typ1) where
+   type MultiplySquareExtra (Symmetric typ0 typOff xlOff xuOff typ1) extra =
+         (Matrix.MultiplySquareExtra typ0 (Symmetric0 extra),
+          Matrix.MultiplySquareExtra typ1 (Symmetric1 extra))
+
+   squareFull (Symmetric a o b) =
+      withVerticalSplit $ \(Above top bottom) ->
+         Above
+            (Multiply.squareFull a top #+# o#*#bottom)
+            (Matrix.transpose o #*# top #+# Multiply.squareFull b bottom)
+   fullSquare x (Symmetric a o b) =
+      withHorizontalSplit x $ \(Beside left right) ->
+         Beside
+            (Multiply.fullSquare left a #+# right #*# Matrix.transpose o)
+            (left#*#o #+# Multiply.fullSquare right b)
+
+instance
+   (typ0 ~ TypeFull, Divide.Solve typ1,
+    typOff ~ TypeFull, xlOff ~ (), xuOff ~ ()) =>
+      Divide.Determinant (Symmetric typ0 typOff xlOff xuOff typ1) where
+   type DeterminantExtra (Symmetric typ0 typOff xlOff xuOff typ1) extra =
+         (Matrix.DeterminantExtra typ1 (Symmetric1 extra),
+          Matrix.Determinant typ1,
+          Divide.SolveExtra typ1 (Symmetric1 extra),
+          Symmetric0 extra ~ ())
+   determinant a@(Symmetric _ _ _) = determinant $ squareFromSymmetric a
+
+instance
+   (typ0 ~ TypeFull, Divide.Solve typ1,
+    typOff ~ TypeFull, xlOff ~ (), xuOff ~ ()) =>
+      Divide.Solve (Symmetric typ0 typOff xlOff xuOff typ1) where
+   type SolveExtra (Symmetric typ0 typOff xlOff xuOff typ1) extra =
+         (Matrix.SolveExtra typ0 (Symmetric0 extra),
+          Matrix.SolveExtra typ1 (Symmetric1 extra))
+
+   solveRight a@(Symmetric _ _ _) =
+      Divide.solveRight $ squareFromSymmetric a
+   solveLeft rhs a@(Symmetric _ _ _) =
+      Divide.solveLeft rhs $ squareFromSymmetric a
+
+
+
+-- * Helper types and functions
+
+type TypeFull = ArrMatrix.Array Layout.Unpacked Omni.Arbitrary
+
+
+withoutSizeCheck ::
+   (Vector (Unchecked sh0) a0 ->
+    Vector (Unchecked sh1) a1 ->
+    Vector (Unchecked sh2) a2) ->
+   Vector sh0 a0 -> Vector sh1 a1 -> Vector sh2 a2
+withoutSizeCheck op a b =
+   Vector.recheck $ Vector.uncheck a `op` Vector.uncheck b
+
+
+withVerticalSplit ::
+   (Extent.Measure meas, Extent.C vert, Extent.C horiz,
+    Shape.C height0, Shape.C height1, Shape.C width, Class.Floating a) =>
+   (AboveGeneral height0 height1 (Unchecked width) a ->
+    AboveGeneral height0 height1 (Unchecked width) a) ->
+   Matrix.Full meas vert horiz (height0::+height1) width a ->
+   Matrix.Full meas vert horiz (height0::+height1) width a
+withVerticalSplit f x =
+   restoreShape x $ Matrix.mapWidth deconsUnchecked $
+      (\(Above a b) -> a===b) $ f $ splitVertically x
+
+withHorizontalSplit ::
+   (Extent.Measure meas, Extent.C vert, Extent.C horiz,
+    Shape.C height, Shape.C width0, Shape.C width1, Class.Floating a) =>
+   Matrix.Full meas vert horiz height (width0::+width1) a ->
+   (BesideGeneral (Unchecked height) width0 width1 a ->
+    BesideGeneral (Unchecked height) width0 width1 a) ->
+   Matrix.Full meas vert horiz height (width0::+width1) a
+withHorizontalSplit x f =
+   restoreShape x $ Matrix.mapHeight deconsUnchecked $
+      (\(Beside a b) -> a|||b) $ f $ splitHorizontally x
+
+splitVertically ::
+   (Extent.Measure meas, Extent.C vert, Extent.C horiz,
+    Shape.C height0, Shape.C height1, Shape.C width, Class.Floating a) =>
+   Matrix.Full meas vert horiz (height0 ::+ height1) width a ->
+   AboveGeneral height0 height1 (Unchecked width) a
+splitVertically a0 =
+   let a = Matrix.mapWidth Unchecked $ Matrix.fromFull a0
+   in Above (Matrix.takeTop a) (Matrix.takeBottom a)
+
+splitHorizontally ::
+   (Extent.Measure meas, Extent.C vert, Extent.C horiz,
+    Shape.C height, Shape.C width0, Shape.C width1, Class.Floating a) =>
+   Matrix.Full meas vert horiz height (width0 ::+ width1) a ->
+   BesideGeneral (Unchecked height) width0 width1 a
+splitHorizontally a0 =
+   let a = Matrix.mapHeight Unchecked $ Matrix.fromFull a0
+   in Beside (Matrix.takeLeft a) (Matrix.takeRight a)
+
+restoreShape ::
+   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+   (Shape.C height, Shape.C width) =>
+   Matrix.Full meas vert horiz height width a ->
+   Matrix.General height width a ->
+   Matrix.Full meas vert horiz height width a
+restoreShape a b =
+   ArrMatrix.reshape
+      (case ArrMatrix.shape a of
+         Omni.Full shape ->
+            Omni.Full (shape{Layout.fullOrder = ArrMatrix.order b}))
+      b
+
+
+
+
+finalizeLayout :: (Shape.C shape) => (shape, [[a]]) -> BoxedArray.Array shape a
+finalizeLayout (shape, a) = BoxedArray.fromList shape $ concat a
+
+infixr 3 |||#
+infixr 2 ===#
+
+(|||#) ::
+   (Eq height) =>
+   ((height, widthA), [[a]]) ->
+   ((height, widthB), [[a]]) ->
+   ((height, widthA ::+ widthB), [[a]])
+((heightA,widthA),a) |||# ((heightB,widthB),b) =
+   if heightA == heightB
+      then ((heightA,widthA::+widthB), zipWith (++) a b)
+      else error "Block.beside: mismatching heights"
+
+(===#) ::
+   (Eq width) =>
+   ((heightA, width), [a]) ->
+   ((heightB, width), [a]) ->
+   ((heightA ::+ heightB, width), [a])
+((heightA,widthA),a) ===# ((heightB,widthB),b) =
+   if widthA == widthB
+      then ((heightA::+heightB,widthA), a ++ b)
+      else error "Block.above: mismatching widths"
+
+toRows ::
+   (Matrix.BoxExtra typ xl, Matrix.BoxExtra typ xu,
+    Layout typ, LayoutExtra typ xl, LayoutExtra typ xu) =>
+   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+   (Shape.C height, Shape.C width, Class.Floating a) =>
+   Matrix typ xl xu lower upper meas vert horiz height width a ->
+   ((height, width), [[(Output.Separator, Maybe (Output.Style, a))]])
+toRows a = ((Matrix.height a, Matrix.width a), ArrFormat.toRows (layout a))
+
+toColumns ::
+   (Matrix.BoxExtra typ xl, Matrix.BoxExtra typ xu,
+    Layout typ, LayoutExtra typ xl, LayoutExtra typ xu) =>
+   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+   (Shape.C height, Shape.C width, Class.Floating a) =>
+   Matrix typ xl xu lower upper meas vert horiz height width a ->
+   ((width, height), [[(Output.Separator, Maybe (Output.Style, a))]])
+toColumns a =
+   ((Matrix.width a, Matrix.height a), ArrFormat.toColumns (layout a))
+
+layoutEmpty ::
+   (Shape.C height, Shape.C width) =>
+   height -> width -> ((height, width), [[(Output.Separator, Maybe a)]])
+layoutEmpty height width =
+   (,) (height, width) $
+   replicate (Shape.size height) $
+   replicate (Shape.size width) (Output.Space, Nothing)
diff --git a/src/Numeric/LAPACK/Matrix/Block/Type.hs b/src/Numeric/LAPACK/Matrix/Block/Type.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Matrix/Block/Type.hs
@@ -0,0 +1,11 @@
+module Numeric.LAPACK.Matrix.Block.Type (
+   Matrix(Diagonal, Above, Beside, Square, Upper, Lower, Symmetric),
+   Block.Diagonal,
+   Block.Append, Block.aboveFromFull, Block.besideFromFull,
+   Block.Square,
+   Block.Triangular,
+   Block.Symmetric, Block.squareFromSymmetric, Block.schurComplement,
+   ) where
+
+import qualified Numeric.LAPACK.Matrix.Block.Private as Block
+import Numeric.LAPACK.Matrix.Type.Private (Matrix)
diff --git a/src/Numeric/LAPACK/Matrix/Class.hs b/src/Numeric/LAPACK/Matrix/Class.hs
--- a/src/Numeric/LAPACK/Matrix/Class.hs
+++ b/src/Numeric/LAPACK/Matrix/Class.hs
@@ -1,17 +1,24 @@
 {-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE ConstraintKinds #-}
 {-# LANGUAGE GADTs #-}
 module Numeric.LAPACK.Matrix.Class (
-   SquareShape(toSquare, takeDiagonal, mapSquareSize, identityFrom),
+   SquareShape(takeDiagonal, identityFrom), SquareShapeExtra, toSquare,
+   MapSquareSize(mapSquareSize),
    MapSize(mapHeight, mapWidth),
    trace,
    Complex(conjugate, fromReal, toComplex),
    adjoint,
-   Unpack(unpack), toFull,
+   Unpack(unpack), UnpackExtra, toFull,
+
+   Homogeneous, HomogeneousExtra, Scale, ScaleExtra,
+   zeroFrom, negate, scaleReal, scale, scaleRealReal, (.*#),
+   Additive, AdditiveExtra, add, (#+#),
+   Subtractive, SubtractiveExtra, sub, (#-#),
    ) where
 
 import qualified Numeric.LAPACK.Matrix.Array.Basic as OmniMatrix
 import qualified Numeric.LAPACK.Matrix.Array.Private as ArrMatrix
-import qualified Numeric.LAPACK.Matrix.Type as Matrix
+import qualified Numeric.LAPACK.Matrix.Type.Private as Matrix
 import qualified Numeric.LAPACK.Matrix.Banded.Basic as Banded
 import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
 import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
@@ -21,7 +28,7 @@
 import qualified Numeric.LAPACK.Permutation as PermPub
 import qualified Numeric.LAPACK.Vector as Vector
 import qualified Numeric.LAPACK.Scalar as Scalar
-import Numeric.LAPACK.Matrix.Type (Matrix)
+import Numeric.LAPACK.Matrix.Type.Private (Matrix)
 import Numeric.LAPACK.Vector (Vector)
 import Numeric.LAPACK.Scalar (RealOf, ComplexOf)
 
@@ -29,7 +36,11 @@
 
 import qualified Data.Array.Comfort.Shape as Shape
 
+import Prelude hiding (negate)
 
+import GHC.Exts (Constraint)
+
+
 class Complex typ where
    conjugate ::
       (Matrix typ xl xu lower upper meas vert horiz height width ~ matrix,
@@ -64,6 +75,7 @@
 
 adjoint ::
    (Matrix.Transpose typ, Complex typ) =>
+   (Matrix.TransposeExtra typ xl, Matrix.TransposeExtra typ xu) =>
    (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
    (Shape.C height, Shape.C width, Class.Floating a) =>
    Matrix typ xl xu lower upper meas vert horiz height width a ->
@@ -72,58 +84,68 @@
 
 
 class (Matrix.Box typ) => SquareShape typ where
-   toSquare ::
-      (Shape.C sh, Class.Floating a) =>
-      Matrix.Quadratic typ xl xu lower upper sh a -> ArrMatrix.Square sh a
+   type SquareShapeExtra typ extra :: Constraint
    takeDiagonal ::
+      (SquareShapeExtra typ xl, SquareShapeExtra typ xu) =>
       (Shape.C sh, Class.Floating a) =>
       Matrix.Quadratic typ xl xu lower upper sh a -> Vector sh a
-   {- |
-   The number of rows and columns
-   must be maintained by the shape mapping function.
-   -}
-   mapSquareSize ::
-      (Shape.C shA, Shape.C shB) =>
-      (shA -> shB) ->
-      Matrix.Quadratic typ xl xu lower upper shA a ->
-      Matrix.Quadratic typ xl xu lower upper shB a
    identityFrom ::
+      (SquareShapeExtra typ xl, SquareShapeExtra typ xu) =>
       (Shape.C sh, Class.Floating a) =>
       Matrix.Quadratic typ xl xu lower upper sh a ->
       Matrix.Quadratic typ xl xu lower upper sh a
 
 instance SquareShape (ArrMatrix.Array pack property) where
-   toSquare a@(ArrMatrix.Array _) = OmniMatrix.toFull a
+   type SquareShapeExtra (ArrMatrix.Array pack property) extra = ()
    takeDiagonal a@(ArrMatrix.Array _) = OmniMatrix.takeDiagonal a
-   mapSquareSize f a@(ArrMatrix.Array _) = OmniMatrix.mapSquareSize f a
    identityFrom a@(ArrMatrix.Array _) = OmniMatrix.identityFrom a
 
 instance SquareShape Matrix.Scale where
-   toSquare (Matrix.Scale sh a) =
-      ArrMatrix.lift0 $ Banded.toFull $
-      Banded.diagonal Layout.RowMajor $ Vector.constant sh a
+   type SquareShapeExtra Matrix.Scale extra = ()
    takeDiagonal (Matrix.Scale sh a) = Vector.constant sh a
-   mapSquareSize f (Matrix.Scale sh a) =
-      Matrix.Scale (Layout.mapChecked "Scale.mapSquareSize" f sh) a
    identityFrom (Matrix.Scale sh _a) = Matrix.Scale sh Scalar.one
 
 instance SquareShape Matrix.Permutation where
-   toSquare (Matrix.Permutation perm) = PermPub.toMatrix perm
+   type SquareShapeExtra Matrix.Permutation extra = ()
    takeDiagonal a@(Matrix.Permutation _) =
       Perm.takeDiagonal . Permutation.toPermutation $ a
-   mapSquareSize f (Matrix.Permutation perm) =
-      Matrix.Permutation $ Perm.mapSize f perm
    identityFrom (Matrix.Permutation perm) =
       Matrix.Permutation $ Perm.identity $ Perm.size perm
 
 
 trace ::
-   (SquareShape typ, Shape.C sh, Class.Floating a) =>
+   (SquareShape typ, SquareShapeExtra typ xl, SquareShapeExtra typ xu) =>
+   (Shape.C sh, Class.Floating a) =>
    Matrix.Quadratic typ xl xu lower upper sh a -> a
 trace = Vector.sum . takeDiagonal
 
 
 
+class (SquareShape typ) => MapSquareSize typ where
+   {- |
+   The number of rows and columns
+   must be maintained by the shape mapping function.
+
+   Not available for `Block` matrices.
+   -}
+   mapSquareSize ::
+      (Shape.C shA, Shape.C shB) =>
+      (shA -> shB) ->
+      Matrix.Quadratic typ xl xu lower upper shA a ->
+      Matrix.Quadratic typ xl xu lower upper shB a
+
+instance MapSquareSize (ArrMatrix.Array pack property) where
+   mapSquareSize f a@(ArrMatrix.Array _) = OmniMatrix.mapSquareSize f a
+
+instance MapSquareSize Matrix.Scale where
+   mapSquareSize f (Matrix.Scale sh a) =
+      Matrix.Scale (Layout.mapChecked "Scale.mapSquareSize" f sh) a
+
+instance MapSquareSize Matrix.Permutation where
+   mapSquareSize f (Matrix.Permutation perm) =
+      Matrix.Permutation $ Perm.mapSize f perm
+
+
 class (Matrix.Box typ) => MapSize typ where
    {- |
    The number of rows and columns
@@ -152,8 +174,10 @@
 
 
 class Unpack typ where
+   type UnpackExtra typ extra :: Constraint
    -- In contrast to OmniMatrix.unpack it cannot maintain the matrix property.
    unpack ::
+      (UnpackExtra typ xl, UnpackExtra typ xu) =>
       (Omni.Strip lower, Omni.Strip upper) =>
       (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
       (Shape.C height, Shape.C width, Class.Floating a) =>
@@ -162,23 +186,182 @@
          lower upper meas vert horiz height width a
 
 instance (Omni.Property prop) => Unpack (ArrMatrix.Array pack prop) where
+   type UnpackExtra (ArrMatrix.Array pack prop) extra = extra ~ ()
    unpack a@(ArrMatrix.Array _) =
       ArrMatrix.liftUnpacked1 id $ OmniMatrix.unpack a
 
 instance Unpack Matrix.Scale where
+   type UnpackExtra Matrix.Scale extra = extra ~ ()
    unpack (Matrix.Scale sh a) =
       ArrMatrix.liftUnpacked0 $ Banded.toFull $
       Banded.diagonal Layout.RowMajor $ Vector.constant sh a
 
 instance Unpack Matrix.Permutation where
+   type UnpackExtra Matrix.Permutation extra = extra ~ ()
    unpack (Matrix.Permutation perm) =
       ArrMatrix.liftUnpacked1 id $ PermPub.toMatrix perm
 
 toFull ::
-   (Unpack typ) =>
+   (Unpack typ, UnpackExtra typ xl, UnpackExtra typ xu) =>
    (Omni.Strip lower, Omni.Strip upper) =>
    (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
    (Shape.C height, Shape.C width, Class.Floating a) =>
    Matrix typ xl xu lower upper meas vert horiz height width a ->
    ArrMatrix.Full meas vert horiz height width a
 toFull = OmniMatrix.toFull . unpack
+
+toSquare ::
+   (Unpack typ, UnpackExtra typ xl, UnpackExtra typ xu) =>
+   (Omni.Strip lower, Omni.Strip upper) =>
+   (Shape.C sh, Class.Floating a) =>
+   Matrix.Quadratic typ xl xu lower upper sh a -> ArrMatrix.Square sh a
+toSquare = toFull
+
+
+
+class Homogeneous typ where
+   type HomogeneousExtra typ extra :: Constraint
+   zeroFrom ::
+      (HomogeneousExtra typ xl, HomogeneousExtra typ xu) =>
+      (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+      (Shape.C height, Shape.C width, Class.Floating a) =>
+      Matrix typ xl xu lower upper meas vert horiz height width a ->
+      Matrix typ xl xu lower upper meas vert horiz height width a
+   negate ::
+      (HomogeneousExtra typ xl, HomogeneousExtra typ xu) =>
+      (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+      (Shape.C height, Shape.C width, Class.Floating a) =>
+      Matrix typ xl xu lower upper meas vert horiz height width a ->
+      Matrix typ xl xu lower upper meas vert horiz height width a
+   scaleReal ::
+      (HomogeneousExtra typ xl, HomogeneousExtra typ xu) =>
+      (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+      (Shape.C height, Shape.C width, Class.Floating a) =>
+      RealOf a ->
+      Matrix typ xl xu lower upper meas vert horiz height width a ->
+      Matrix typ xl xu lower upper meas vert horiz height width a
+
+instance
+   (ArrMatrix.Homogeneous property) =>
+      Homogeneous (ArrMatrix.Array pack property) where
+   type HomogeneousExtra (ArrMatrix.Array pack property) extra = extra ~ ()
+   zeroFrom = ArrMatrix.zero . ArrMatrix.shape
+   negate = ArrMatrix.negate
+   scaleReal = ArrMatrix.scaleReal
+
+instance Homogeneous Matrix.Scale where
+   type HomogeneousExtra Matrix.Scale extra = extra ~ ()
+   zeroFrom (Matrix.Scale sh _a) = Matrix.Scale sh Scalar.zero
+   negate (Matrix.Scale sh a) = Matrix.Scale sh (-a)
+   scaleReal c (Matrix.Scale sh a) = Matrix.Scale sh (Scalar.fromReal c*a)
+
+newtype ScaleReal f a = ScaleReal {getScaleReal :: a -> f a -> f a}
+
+scaleRealReal ::
+   (Homogeneous typ, HomogeneousExtra typ xl, HomogeneousExtra typ xu) =>
+   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+   (Shape.C height, Shape.C width, Class.Real a) =>
+   a ->
+   Matrix typ xl xu lower upper meas vert horiz height width a ->
+   Matrix typ xl xu lower upper meas vert horiz height width a
+scaleRealReal =
+   getScaleReal $ Class.switchReal (ScaleReal scaleReal) (ScaleReal scaleReal)
+
+
+class (Homogeneous typ) => Scale typ where
+   type ScaleExtra typ extra :: Constraint
+   scale ::
+      (ScaleExtra typ xl, ScaleExtra typ xu) =>
+      (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+      (Shape.C height, Shape.C width, Class.Floating a) =>
+      a ->
+      Matrix typ xl xu lower upper meas vert horiz height width a ->
+      Matrix typ xl xu lower upper meas vert horiz height width a
+
+instance
+   (ArrMatrix.Scale property) =>
+      Scale (ArrMatrix.Array pack property) where
+   type ScaleExtra (ArrMatrix.Array pack property) extra = extra ~ ()
+   scale = ArrMatrix.scale
+
+instance Scale Matrix.Scale where
+   type ScaleExtra Matrix.Scale extra = extra ~ ()
+   scale c (Matrix.Scale sh a) = Matrix.Scale sh (c*a)
+
+(.*#) ::
+   (Scale typ, ScaleExtra typ xl, ScaleExtra typ xu) =>
+   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+   (Shape.C height, Shape.C width, Class.Floating a) =>
+   a ->
+   Matrix typ xl xu lower upper meas vert horiz height width a ->
+   Matrix typ xl xu lower upper meas vert horiz height width a
+(.*#) = scale
+
+infixl 7 .*#
+
+
+class Additive typ where
+   type AdditiveExtra typ extra :: Constraint
+   add ::
+      (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+      (AdditiveExtra typ xl, AdditiveExtra typ xu,
+       Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>
+      Matrix typ xl xu lower upper meas vert horiz height width a ->
+      Matrix typ xl xu lower upper meas vert horiz height width a ->
+      Matrix typ xl xu lower upper meas vert horiz height width a
+
+instance
+   (ArrMatrix.Additive property) =>
+      Additive (ArrMatrix.Array pack property) where
+   type AdditiveExtra (ArrMatrix.Array pack property) extra = extra ~ ()
+   add = ArrMatrix.add
+
+instance Additive Matrix.Scale where
+   type AdditiveExtra Matrix.Scale extra = extra ~ ()
+   add (Matrix.Scale sha a) (Matrix.Scale shb b) =
+      if sha == shb
+         then Matrix.Scale sha (a+b)
+         else error "Matrix.add Scale: dimensions mismatch"
+
+class (Additive typ) => Subtractive typ where
+   type SubtractiveExtra typ extra :: Constraint
+   sub ::
+      (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+      (SubtractiveExtra typ xl, SubtractiveExtra typ xu,
+       Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>
+      Matrix typ xl xu lower upper meas vert horiz height width a ->
+      Matrix typ xl xu lower upper meas vert horiz height width a ->
+      Matrix typ xl xu lower upper meas vert horiz height width a
+
+instance
+   (ArrMatrix.Subtractive property) =>
+      Subtractive (ArrMatrix.Array pack property) where
+   type SubtractiveExtra (ArrMatrix.Array pack property) extra = extra ~ ()
+   sub = ArrMatrix.sub
+
+instance Subtractive Matrix.Scale where
+   type SubtractiveExtra Matrix.Scale extra = extra ~ ()
+   sub (Matrix.Scale sha a) (Matrix.Scale shb b) =
+      if sha == shb
+         then Matrix.Scale sha (a-b)
+         else error "Matrix.sub Scale: dimensions mismatch"
+
+infixl 6 #+#, #-#, `add`, `sub`
+
+(#+#) ::
+   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+   (Additive typ, AdditiveExtra typ xl, AdditiveExtra typ xu,
+    Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>
+   Matrix typ xl xu lower upper meas vert horiz height width a ->
+   Matrix typ xl xu lower upper meas vert horiz height width a ->
+   Matrix typ xl xu lower upper meas vert horiz height width a
+(#+#) = add
+
+(#-#) ::
+   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+   (Subtractive typ, SubtractiveExtra typ xl, SubtractiveExtra typ xu,
+    Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>
+   Matrix typ xl xu lower upper meas vert horiz height width a ->
+   Matrix typ xl xu lower upper meas vert horiz height width a ->
+   Matrix typ xl xu lower upper meas vert horiz height width a
+(#-#) = sub
diff --git a/src/Numeric/LAPACK/Matrix/Divide.hs b/src/Numeric/LAPACK/Matrix/Divide.hs
--- a/src/Numeric/LAPACK/Matrix/Divide.hs
+++ b/src/Numeric/LAPACK/Matrix/Divide.hs
@@ -2,7 +2,8 @@
 {-# LANGUAGE GADTs #-}
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE FlexibleContexts #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE ConstraintKinds #-}
+{-# LANGUAGE UndecidableInstances #-}
 module Numeric.LAPACK.Matrix.Divide where
 
 import qualified Numeric.LAPACK.Matrix.Multiply as Multiply
@@ -10,14 +11,14 @@
 import qualified Numeric.LAPACK.Matrix.Array.Unpacked as Unpacked
 import qualified Numeric.LAPACK.Matrix.Array.Private as ArrMatrix
 import qualified Numeric.LAPACK.Matrix.Permutation as PermMatrix
-import qualified Numeric.LAPACK.Matrix.Type as Matrix
+import qualified Numeric.LAPACK.Matrix.Type.Private as Matrix
 import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
 import qualified Numeric.LAPACK.Matrix.Shape.Omni as Omni
 import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
 import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
 import qualified Numeric.LAPACK.Vector as Vector
 import Numeric.LAPACK.Matrix.Array.Private (Full)
-import Numeric.LAPACK.Matrix.Type (scaleWithCheck)
+import Numeric.LAPACK.Matrix.Type.Private (scaleWithCheck)
 import Numeric.LAPACK.Matrix.Modifier
          (Transposition(NonTransposed,Transposed),
           Inversion(Inverted))
@@ -30,16 +31,22 @@
 
 import Data.Semigroup ((<>))
 
+import GHC.Exts (Constraint)
 
-class (Matrix.Box typ) => Determinant typ xl xu where
+
+class (Matrix.Box typ) => Determinant typ where
+   type DeterminantExtra typ extra :: Constraint
    determinant ::
+      (DeterminantExtra typ xl, DeterminantExtra typ xu) =>
       (Omni.Strip lower, Omni.Strip upper) =>
       (Shape.C sh, Class.Floating a) =>
       Matrix.Quadratic typ xl xu lower upper sh a -> a
 
-class (Matrix.Box typ) => Solve typ xl xu where
+class (Matrix.Box typ) => Solve typ where
+   type SolveExtra typ extra :: Constraint
    {-# MINIMAL solve | solveLeft,solveRight #-}
    solve ::
+      (SolveExtra typ xl, SolveExtra typ xu) =>
       (Omni.Strip lower, Omni.Strip upper) =>
       (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
       (Shape.C height, Shape.C width, Eq height, Class.Floating a) =>
@@ -52,6 +59,7 @@
       Unpacked.transpose $ solveLeft (Unpacked.transpose b) a
 
    solveRight ::
+      (SolveExtra typ xl, SolveExtra typ xu) =>
       (Omni.Strip lower, Omni.Strip upper) =>
       (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
       (Shape.C height, Shape.C width, Eq height, Class.Floating a) =>
@@ -61,6 +69,7 @@
    solveRight = solve NonTransposed
 
    solveLeft ::
+      (SolveExtra typ xl, SolveExtra typ xu) =>
       (Omni.Strip lower, Omni.Strip upper) =>
       (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
       (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
@@ -69,8 +78,10 @@
       Full meas vert horiz height width a
    solveLeft = Unpacked.swapMultiply $ solve Transposed
 
-class (Solve typ xl xu) => Inverse typ xl xu where
+class (Solve typ) => Inverse typ where
+   type InverseExtra typ extra :: Constraint
    inverse ::
+      (InverseExtra typ xl, InverseExtra typ xu) =>
       (MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper,
        Extent.Measure meas, Shape.C height, Shape.C width, Class.Floating a) =>
       Matrix.QuadraticMeas typ xl xu lower upper meas height width a ->
@@ -80,7 +91,9 @@
 infixr 7 #\##
 
 (#\##) ::
-   (Solve typ xl xu, Matrix.ToQuadratic typ,
+   (Solve typ, Matrix.ToQuadratic typ,
+    SolveExtra typ xl, SolveExtra typ xu,
+    Matrix.BoxExtra typ xl, Matrix.BoxExtra typ xu,
     Omni.Strip lower, Omni.Strip upper,
     Shape.C height, Eq height, Shape.C width, Shape.C nrhs,
     Extent.Measure measA, Extent.Measure measB, Extent.Measure measC,
@@ -94,7 +107,9 @@
          Multiply.reshapeHeight (Matrix.transpose ident) (solveRight q b)
 
 (##/#) ::
-   (Solve typ xl xu, Matrix.ToQuadratic typ,
+   (Solve typ, Matrix.ToQuadratic typ,
+    SolveExtra typ xl, SolveExtra typ xu,
+    Matrix.BoxExtra typ xl, Matrix.BoxExtra typ xu,
     Omni.Strip lower, Omni.Strip upper,
     Shape.C height, Shape.C width, Eq width, Shape.C nrhs,
     Extent.Measure measA, Extent.Measure measB, Extent.Measure measC,
@@ -110,7 +125,8 @@
 
 
 solveVector ::
-   (Solve typ xl xu, Omni.Strip lower, Omni.Strip upper,
+   (Solve typ, SolveExtra typ xl, SolveExtra typ xu,
+    Omni.Strip lower, Omni.Strip upper,
     Shape.C sh, Eq sh, Class.Floating a) =>
    Transposition ->
    Matrix.Quadratic typ xl xu lower upper sh a ->
@@ -123,7 +139,9 @@
 infixr 7 #\|
 
 (#\|) ::
-   (Solve typ xl xu, Matrix.ToQuadratic typ,
+   (Solve typ, Matrix.ToQuadratic typ,
+    SolveExtra typ xl, SolveExtra typ xu,
+    Matrix.BoxExtra typ xl, Matrix.BoxExtra typ xu,
     Omni.Strip lower, Omni.Strip upper, Extent.Measure meas,
     Shape.C height, Shape.C width, Eq height, Class.Floating a) =>
    Matrix.QuadraticMeas typ xl xu lower upper meas height width a ->
@@ -134,7 +152,9 @@
          reshapeVector (Matrix.transpose ident) . solveVector NonTransposed q
 
 (-/#) ::
-   (Solve typ xl xu, Matrix.ToQuadratic typ,
+   (Solve typ, Matrix.ToQuadratic typ,
+    SolveExtra typ xl, SolveExtra typ xu,
+    Matrix.BoxExtra typ xl, Matrix.BoxExtra typ xu,
     Omni.Strip lower, Omni.Strip upper, Extent.Measure meas,
     Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
    Vector width a ->
@@ -152,44 +172,53 @@
 reshapeVector (Matrix.Identity extent) = Array.reshape (Extent.height extent)
 
 
-instance (xl ~ (), xu ~ ()) => Determinant Matrix.Scale xl xu where
+instance Determinant Matrix.Scale where
+   type DeterminantExtra Matrix.Scale extra = extra ~ ()
    determinant (Matrix.Scale sh a) = a ^ Shape.size sh
 
-instance (xl ~ (), xu ~ ()) => Solve Matrix.Scale xl xu where
+instance Solve Matrix.Scale where
+   type SolveExtra Matrix.Scale extra = extra ~ ()
    solve _trans =
       scaleWithCheck "Matrix.Scale.solve" Matrix.height $
          ArrMatrix.lift1 . Vector.scale . recip
 
-instance (xl ~ (), xu ~ ()) => Inverse Matrix.Scale xl xu where
+instance Inverse Matrix.Scale where
+   type InverseExtra Matrix.Scale extra = extra ~ ()
    inverse (Matrix.Scale shape a) = Matrix.Scale shape $ recip a
 
 
-instance (xl ~ (), xu ~ ()) => Determinant Matrix.Permutation xl xu where
+instance Determinant Matrix.Permutation where
+   type DeterminantExtra Matrix.Permutation extra = extra ~ ()
    determinant = PermMatrix.determinant
 
-instance (xl ~ (), xu ~ ()) => Solve Matrix.Permutation xl xu where
+instance Solve Matrix.Permutation where
+   type SolveExtra Matrix.Permutation extra = extra ~ ()
    solve trans =
       PermMatrix.multiplyFull
          (Inverted <> PermMatrix.inversionFromTransposition trans)
 
-instance (xl ~ (), xu ~ ()) => Inverse Matrix.Permutation xl xu where
+instance Inverse Matrix.Permutation where
+   type InverseExtra Matrix.Permutation extra = extra ~ ()
    inverse a@(Matrix.Permutation _) =
       case Matrix.powerStrips a of
          (MatrixShape.Filled, MatrixShape.Filled) -> PermMatrix.transpose a
          _ -> a -- identity matrix
 
 instance
-   (Layout.Packing pack, Omni.Property property, xl ~ (), xu ~ ()) =>
-      Determinant (ArrMatrix.Array pack property) xl xu where
+   (Layout.Packing pack, Omni.Property property) =>
+      Determinant (ArrMatrix.Array pack property) where
+   type DeterminantExtra (ArrMatrix.Array pack property) extra = extra ~ ()
    determinant = ArrDivide.determinant
 
 instance
-   (Layout.Packing pack, Omni.Property property, xl ~ (), xu ~ ()) =>
-      Solve (ArrMatrix.Array pack property) xl xu where
+   (Layout.Packing pack, Omni.Property property) =>
+      Solve (ArrMatrix.Array pack property) where
+   type SolveExtra (ArrMatrix.Array pack property) extra = extra ~ ()
    solveRight = ArrDivide.solve
    solveLeft = Matrix.swapMultiply $  ArrDivide.solve . Matrix.transpose
 
 instance
-   (Layout.Packing pack, Omni.Property property, xl ~ (), xu ~ ()) =>
-      Inverse (ArrMatrix.Array pack property) xl xu where
+   (Layout.Packing pack, Omni.Property property) =>
+      Inverse (ArrMatrix.Array pack property) where
+   type InverseExtra (ArrMatrix.Array pack property) extra = extra ~ ()
    inverse = ArrDivide.inverse
diff --git a/src/Numeric/LAPACK/Matrix/Extent/Strict.hs b/src/Numeric/LAPACK/Matrix/Extent/Strict.hs
--- a/src/Numeric/LAPACK/Matrix/Extent/Strict.hs
+++ b/src/Numeric/LAPACK/Matrix/Extent/Strict.hs
@@ -15,6 +15,13 @@
 newtype Map measA vertA horizA measB vertB horizB height width =
    Map {apply :: Extent.Map measA vertA horizA measB vertB horizB height width}
 
+transpose ::
+   (Measure measA, C vertA, C horizA) =>
+   (Measure measB, C vertB, C horizB) =>
+   Map measA vertA horizA measB vertB horizB height width ->
+   Map measA horizA vertA measB horizB vertB width height
+transpose (Map m) = Map (Extent.transpose . m . Extent.transpose)
+
 
 {- |
 Admissible tag combinations are:
diff --git a/src/Numeric/LAPACK/Matrix/Indexed.hs b/src/Numeric/LAPACK/Matrix/Indexed.hs
--- a/src/Numeric/LAPACK/Matrix/Indexed.hs
+++ b/src/Numeric/LAPACK/Matrix/Indexed.hs
@@ -5,9 +5,9 @@
 import qualified Numeric.LAPACK.Matrix.Array.Indexed as ArrIndexed
 import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
 import qualified Numeric.LAPACK.Matrix.Array.Private as ArrMatrix
-import qualified Numeric.LAPACK.Matrix.Type as Matrix
+import qualified Numeric.LAPACK.Matrix.Type.Private as Matrix
 import qualified Numeric.LAPACK.Permutation as Perm
-import Numeric.LAPACK.Matrix.Type (Matrix)
+import Numeric.LAPACK.Matrix.Type.Private (Matrix)
 import Numeric.LAPACK.Scalar (one, zero)
 import Numeric.LAPACK.Permutation.Private (Permutation(Permutation))
 
diff --git a/src/Numeric/LAPACK/Matrix/Inverse.hs b/src/Numeric/LAPACK/Matrix/Inverse.hs
--- a/src/Numeric/LAPACK/Matrix/Inverse.hs
+++ b/src/Numeric/LAPACK/Matrix/Inverse.hs
@@ -2,10 +2,11 @@
 {-# LANGUAGE GADTs #-}
 {-# LANGUAGE EmptyDataDecls #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE ConstraintKinds #-}
 {-# LANGUAGE UndecidableInstances #-}
 module Numeric.LAPACK.Matrix.Inverse where
 
-import qualified Numeric.LAPACK.Matrix.Type as Matrix
+import qualified Numeric.LAPACK.Matrix.Type.Private as Matrix
 import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
 import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
 import qualified Numeric.LAPACK.Matrix.Shape.Omni as Omni
@@ -14,28 +15,35 @@
 import qualified Numeric.LAPACK.Matrix.Divide as Divide
 import qualified Numeric.LAPACK.Matrix.Multiply as Multiply
 import Numeric.LAPACK.Matrix.Divide ((#\|), (-/#))
-import Numeric.LAPACK.Matrix.Type (Matrix)
+import Numeric.LAPACK.Matrix.Type.Private (Matrix)
 
 import qualified Type.Data.Num.Unary as Unary
 
 
+{- |
+You may wrap non-diagonal band matrices in 'Wrapper.FillStrips' first
+in order to meet the 'Omni.PowerStrip' constraint.
+-}
 data Inverse typ
 data instance
    Matrix (Inverse typ)
       extraLower extraUpper lowerf upperf meas vert horiz height width a
          where
       Inverse ::
-         (Omni.Strip lower, Fill lower ~ lowerf, Omni.PowerStrip lowerf,
-          Omni.Strip upper, Fill upper ~ upperf, Omni.PowerStrip upperf) =>
          Matrix.QuadraticMeas typ xl xu upper lower meas width height a ->
          Matrix.QuadraticMeas (Inverse typ) (xl,lower) (xu,upper)
-            lowerf upperf meas height width a
+            lower upper meas height width a
 
 type family Fill offDiag
 type instance Fill (Layout.Bands Unary.Zero) = Layout.Bands Unary.Zero
 type instance Fill (Layout.Bands (Unary.Succ k)) = Layout.Filled
 type instance Fill Layout.Filled = Layout.Filled
 
+type family InverseExtra xl
+type instance InverseExtra (xl,lower) = xl
+type family InverseStrip xl
+type instance InverseStrip (xl,lower) = lower
+
 data PowerStripFact c = (Omni.PowerStrip c) => PowerStripFact
 
 filledPowerStripFact ::
@@ -49,15 +57,18 @@
 
 
 instance (Matrix.Transpose typ) => Matrix.Transpose (Inverse typ) where
+   type TransposeExtra (Inverse typ) extra =
+         Matrix.TransposeExtra typ (InverseExtra extra)
    transpose (Inverse a) = Inverse $ Matrix.transpose a
 
-instance
-   (Matrix.MultiplySame typ xl xu,
-    MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper) =>
-      Matrix.MultiplySame (Inverse typ) (xl,lower) (xu,upper) where
+instance (Matrix.MultiplySame typ) => Matrix.MultiplySame (Inverse typ) where
+   type MultiplySameExtra (Inverse typ) extra =
+         (Matrix.MultiplySameExtra typ (InverseExtra extra),
+          MatrixShape.PowerStrip (InverseStrip extra))
    multiplySame (Inverse a) (Inverse b) = Inverse $ Matrix.multiplySame b a
 
 instance (Matrix.Box typ) => Matrix.Box (Inverse typ) where
+   type BoxExtra (Inverse typ) extra = Matrix.BoxExtra typ (InverseExtra extra)
    extent (Inverse m) = Extent.transpose $ Matrix.extent m
    height (Inverse m) = Matrix.width m
    width (Inverse m) = Matrix.height m
@@ -73,45 +84,59 @@
 
 
 instance
-   (Divide.Solve typ xl xu, Matrix.ToQuadratic typ,
-    Omni.Strip lower, Omni.Strip upper) =>
-      Multiply.MultiplyVector (Inverse typ) (xl,lower) (xu,upper) where
+   (Divide.Solve typ, Matrix.ToQuadratic typ) =>
+      Multiply.MultiplyVector (Inverse typ) where
+   type MultiplyVectorExtra (Inverse typ) extra =
+         (Multiply.MultiplyVectorExtra typ (InverseExtra extra),
+          Divide.SolveExtra typ (InverseExtra extra),
+          Matrix.BoxExtra typ (InverseExtra extra),
+          Omni.Strip (InverseStrip extra))
    matrixVector (Inverse a) x = a#\|x
    vectorMatrix x (Inverse a) = x-/#a
 
 instance
-   (Divide.Solve typ xl xu, Matrix.ToQuadratic typ,
-    Omni.Strip lower, Omni.Strip upper) =>
-      Multiply.MultiplySquare (Inverse typ) (xl,lower) (xu,upper) where
+   (Divide.Solve typ, Matrix.ToQuadratic typ) =>
+      Multiply.MultiplySquare (Inverse typ) where
+   type MultiplySquareExtra (Inverse typ) extra =
+         (Multiply.MultiplySquareExtra typ (InverseExtra extra),
+          Divide.SolveExtra typ (InverseExtra extra),
+          Matrix.BoxExtra typ (InverseExtra extra),
+          Omni.Strip (InverseStrip extra))
    transposableSquare trans (Inverse a) = Divide.solve trans a
    squareFull (Inverse a) b = Divide.solveRight a b
    fullSquare b (Inverse a) = Divide.solveLeft b a
 
-instance
-   (Multiply.Power typ xl xu,
-    MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper) =>
-      Multiply.Power (Inverse typ) (xl,lower) (xu,upper) where
+instance (Multiply.Power typ) => Multiply.Power (Inverse typ) where
+   type PowerExtra (Inverse typ) extra =
+         (Multiply.PowerExtra typ (InverseExtra extra),
+          MatrixShape.PowerStrip (InverseStrip extra))
    square (Inverse a) = Inverse $ Multiply.square a
    power n (Inverse a) = Inverse $ Multiply.power n a
    powers1 (Inverse a) = fmap Inverse $ Multiply.powers1 a
 
 
-instance
-   (Divide.Determinant typ xl xu,
-    MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper) =>
-      Divide.Determinant (Inverse typ) (xl,lower) (xu,upper) where
+instance (Divide.Determinant typ) => Divide.Determinant (Inverse typ) where
+   type DeterminantExtra (Inverse typ) extra =
+         (Divide.DeterminantExtra typ (InverseExtra extra),
+          MatrixShape.PowerStrip (InverseStrip extra))
    determinant (Inverse a) = recip $ Divide.determinant a
 
 instance
-   (Multiply.MultiplySquare typ xl xu, Matrix.ToQuadratic typ,
-    MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper) =>
-      Divide.Solve (Inverse typ) (xl,lower) (xu,upper) where
+   (Multiply.MultiplySquare typ, Matrix.ToQuadratic typ) =>
+      Divide.Solve (Inverse typ) where
+   type SolveExtra (Inverse typ) extra =
+         (Divide.SolveExtra typ (InverseExtra extra),
+          Multiply.MultiplySquareExtra typ (InverseExtra extra),
+          MatrixShape.PowerStrip (InverseStrip extra))
    solve trans (Inverse a) = Multiply.transposableSquare trans a
    solveRight (Inverse a) b = Multiply.squareFull a b
    solveLeft b (Inverse a) = Multiply.fullSquare b a
 
 instance
-   (Divide.Inverse typ xl xu, Multiply.MultiplySquare typ xl xu,
-    Matrix.ToQuadratic typ, MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper) =>
-      Divide.Inverse (Inverse typ) (xl,lower) (xu,upper) where
+   (Divide.Inverse typ, Multiply.MultiplySquare typ, Matrix.ToQuadratic typ) =>
+      Divide.Inverse (Inverse typ) where
+   type InverseExtra (Inverse typ) extra =
+         (Divide.InverseExtra typ (InverseExtra extra),
+          Multiply.MultiplySquareExtra typ (InverseExtra extra),
+          MatrixShape.PowerStrip (InverseStrip extra))
    inverse (Inverse a) = Inverse $ Divide.inverse a
diff --git a/src/Numeric/LAPACK/Matrix/Layout/Private.hs b/src/Numeric/LAPACK/Matrix/Layout/Private.hs
--- a/src/Numeric/LAPACK/Matrix/Layout/Private.hs
+++ b/src/Numeric/LAPACK/Matrix/Layout/Private.hs
@@ -474,6 +474,10 @@
 
 
 newtype Bands offDiag = Bands (UnaryProxy offDiag) deriving (Eq, Show)
+
+type family GetBands strip
+type instance GetBands (Bands offDiag) = offDiag
+
 type Empty = Bands TypeNum.U0
 data Filled = Filled deriving (Eq, Show)
 
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
@@ -2,6 +2,7 @@
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE UndecidableInstances #-}
 {-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE ConstraintKinds #-}
 {-# LANGUAGE GADTs #-}
 module Numeric.LAPACK.Matrix.Multiply where
 
@@ -10,7 +11,7 @@
 import qualified Numeric.LAPACK.Matrix.Array.Private as ArrMatrix
 import qualified Numeric.LAPACK.Matrix.Permutation as PermMatrix
 import qualified Numeric.LAPACK.Matrix.Class as MatrixClass
-import qualified Numeric.LAPACK.Matrix.Type as Matrix
+import qualified Numeric.LAPACK.Matrix.Type.Private as Matrix
 import qualified Numeric.LAPACK.Matrix.Basic as FullBasic
 import qualified Numeric.LAPACK.Matrix.Modifier as Mod
 import qualified Numeric.LAPACK.Matrix.Shape.Omni as Omni
@@ -21,7 +22,7 @@
 import qualified Numeric.LAPACK.Permutation.Private as Perm
 import qualified Numeric.LAPACK.Vector as Vector
 import Numeric.LAPACK.Matrix.Array.Private (Full)
-import Numeric.LAPACK.Matrix.Type (Matrix, scaleWithCheck)
+import Numeric.LAPACK.Matrix.Type.Private (Matrix, scaleWithCheck)
 import Numeric.LAPACK.Matrix.Modifier (Transposition(NonTransposed,Transposed))
 import Numeric.LAPACK.Vector (Vector)
 
@@ -34,13 +35,16 @@
 import Data.Stream (Stream)
 import Data.Maybe (fromMaybe)
 
+import GHC.Exts (Constraint)
 
 
+
 infixl 7 -*#
 infixr 7 #*|
 
 (#*|) ::
-   (MultiplyVector typ xl xu, Omni.Strip lower, Omni.Strip upper) =>
+   (MultiplyVector typ, Omni.Strip lower, Omni.Strip upper) =>
+   (MultiplyVectorExtra typ xl, MultiplyVectorExtra typ xu) =>
    (Extent.Measure meas, Extent.C vert, Extent.C horiz,
     Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
    Matrix typ xl xu lower upper meas vert horiz height width a ->
@@ -48,7 +52,8 @@
 (#*|) = matrixVector
 
 (-*#) ::
-   (MultiplyVector typ xl xu, Omni.Strip lower, Omni.Strip upper) =>
+   (MultiplyVector typ, Omni.Strip lower, Omni.Strip upper) =>
+   (MultiplyVectorExtra typ xl, MultiplyVectorExtra typ xu) =>
    (Extent.Measure meas, Extent.C vert, Extent.C horiz,
     Shape.C height, Shape.C width, Eq height, Class.Floating a) =>
    Vector height a ->
@@ -57,8 +62,10 @@
 (-*#) = vectorMatrix
 
 
-class (Matrix.Box typ) => MultiplyVector typ xl xu where
+class (Matrix.Box typ) => MultiplyVector typ where
+   type MultiplyVectorExtra typ extra :: Constraint
    matrixVector ::
+      (MultiplyVectorExtra typ xl, MultiplyVectorExtra typ xu) =>
       (Omni.Strip lower, Omni.Strip upper) =>
       (Extent.Measure meas, Extent.C vert, Extent.C horiz,
        Shape.C height, Shape.C width) =>
@@ -66,6 +73,7 @@
       Matrix typ xl xu lower upper meas vert horiz height width a ->
       Vector width a -> Vector height a
    vectorMatrix ::
+      (MultiplyVectorExtra typ xl, MultiplyVectorExtra typ xu) =>
       (Omni.Strip lower, Omni.Strip upper) =>
       (Extent.Measure meas, Extent.C vert, Extent.C horiz,
        Shape.C height, Shape.C width) =>
@@ -74,7 +82,8 @@
       Matrix typ xl xu lower upper meas vert horiz height width a ->
       Vector width a
 
-instance (xl ~ (), xu ~ ()) => MultiplyVector Matrix.Scale xl xu where
+instance MultiplyVector Matrix.Scale where
+   type MultiplyVectorExtra Matrix.Scale extra = extra ~ ()
    matrixVector a@(Matrix.Scale _ _) =
       scaleWithCheck "Matrix.Multiply.matrixVector Scale"
          Array.shape Vector.scale a
@@ -82,23 +91,27 @@
       scaleWithCheck "Matrix.Multiply.vectorMatrix Scale"
          Array.shape Vector.scale a v
 
-instance (xl ~ (), xu ~ ()) => MultiplyVector Matrix.Permutation xl xu where
+instance MultiplyVector Matrix.Permutation where
+   type MultiplyVectorExtra Matrix.Permutation extra = extra ~ ()
    matrixVector a@(Matrix.Permutation _) =
       PermMatrix.multiplyVector Mod.NonInverted a
    vectorMatrix v a@(Matrix.Permutation _) =
       PermMatrix.multiplyVector Mod.Inverted a v
 
 instance
-   (Layout.Packing pack, Omni.Property property, xl ~ (), xu ~ ()) =>
-      MultiplyVector (ArrMatrix.Array pack property) xl xu where
+   (Layout.Packing pack, Omni.Property property) =>
+      MultiplyVector (ArrMatrix.Array pack property) where
+   type MultiplyVectorExtra (ArrMatrix.Array pack property) extra = extra ~ ()
    matrixVector = Multiply.matrixVector
    vectorMatrix = Multiply.vectorMatrix
 
 
 
-class (Matrix.Box typ) => MultiplySquare typ xl xu where
+class (Matrix.Box typ) => MultiplySquare typ where
+   type MultiplySquareExtra typ extra :: Constraint
    {-# MINIMAL transposableSquare | fullSquare,squareFull #-}
    transposableSquare ::
+      (MultiplySquareExtra typ xl, MultiplySquareExtra typ xu) =>
       (Omni.Strip lower, Omni.Strip upper) =>
       (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
       (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
@@ -111,6 +124,7 @@
       Unpacked.transpose $ fullSquare (Unpacked.transpose b) a
 
    squareFull ::
+      (MultiplySquareExtra typ xl, MultiplySquareExtra typ xu) =>
       (Omni.Strip lower, Omni.Strip upper) =>
       (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
       (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
@@ -120,6 +134,7 @@
    squareFull = transposableSquare NonTransposed
 
    fullSquare ::
+      (MultiplySquareExtra typ xl, MultiplySquareExtra typ xu) =>
       (Omni.Strip lower, Omni.Strip upper) =>
       (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
       (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
@@ -139,7 +154,9 @@
 infixr 7 #*##
 
 (#*##) ::
-   (MultiplySquare typ xl xu, Matrix.ToQuadratic typ) =>
+   (MultiplySquare typ, Matrix.ToQuadratic typ) =>
+   (MultiplySquareExtra typ xl, MultiplySquareExtra typ xu) =>
+   (Matrix.BoxExtra typ xl, Matrix.BoxExtra typ xu) =>
    (Omni.Strip lowerA, Omni.Strip upperA) =>
    (Omni.Strip lowerB, Omni.Strip upperB) =>
    (Omni.Strip lowerC, Omni.Strip upperC) =>
@@ -161,7 +178,9 @@
          reshapeHeight ident (squareFull q (ArrMatrix.liftUnpacked1 id b))
 
 (##*#) ::
-   (MultiplySquare typ xl xu, Matrix.ToQuadratic typ) =>
+   (MultiplySquare typ, Matrix.ToQuadratic typ) =>
+   (MultiplySquareExtra typ xl, MultiplySquareExtra typ xu) =>
+   (Matrix.BoxExtra typ xl, Matrix.BoxExtra typ xu) =>
    (Omni.Strip lowerA, Omni.Strip upperA) =>
    (Omni.Strip lowerB, Omni.Strip upperB) =>
    (Omni.Strip lowerC, Omni.Strip upperC) =>
@@ -189,6 +208,7 @@
 
 factorIdentityLeft ::
    (Matrix.ToQuadratic typ, Extent.Measure meas) =>
+   (Matrix.BoxExtra typ xl, Matrix.BoxExtra typ xu) =>
    Matrix.QuadraticMeas typ xl xu lower upper meas height width a ->
    (IdentityMaes meas height width a,
     Matrix.Quadratic typ xl xu lower upper width a)
@@ -197,6 +217,7 @@
 
 factorIdentityRight ::
    (Matrix.ToQuadratic typ, Extent.Measure meas) =>
+   (Matrix.BoxExtra typ xl, Matrix.BoxExtra typ xu) =>
    Matrix.QuadraticMeas typ xl xu lower upper meas height width a ->
    (Matrix.Quadratic typ xl xu lower upper height a,
     IdentityMaes meas height width a)
@@ -234,62 +255,75 @@
             (Extent.transpose extentA) (Extent.transpose extentB))
 
 
-instance (xl ~ (), xu ~ ()) => MultiplySquare Matrix.Scale xl xu where
+instance MultiplySquare Matrix.Scale where
+   type MultiplySquareExtra Matrix.Scale extra = extra ~ ()
    transposableSquare _trans =
       scaleWithCheck "Matrix.Multiply.transposableSquare" Matrix.height $
          ArrMatrix.lift1 . Vector.scale
 
-instance (xl ~ (), xu ~ ()) => MultiplySquare Matrix.Permutation xl xu where
+instance MultiplySquare Matrix.Permutation where
+   type MultiplySquareExtra Matrix.Permutation extra = extra ~ ()
    transposableSquare =
       PermMatrix.multiplyFull . Perm.inversionFromTransposition
 
 instance
-   (Layout.Packing pack, Omni.Property property, xl ~ (), xu ~ ()) =>
-      MultiplySquare (ArrMatrix.Array pack property) xl xu where
+   (Layout.Packing pack, Omni.Property property) =>
+      MultiplySquare (ArrMatrix.Array pack property) where
+   type MultiplySquareExtra (ArrMatrix.Array pack property) extra = extra ~ ()
    transposableSquare = Multiply.transposableSquare
    fullSquare = Multiply.fullSquare
    squareFull = Multiply.squareFull
 
 
-class (Matrix.Box typ) => Power typ xl xu where
+class (Matrix.Box typ) => Power typ where
+   type PowerExtra typ extra :: Constraint
    square ::
+      (PowerExtra typ xl, PowerExtra typ xu) =>
       (MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper) =>
       (Shape.C sh, Class.Floating a) =>
       Matrix.Quadratic typ xl xu lower upper sh a ->
       Matrix.Quadratic typ xl xu lower upper sh a
    power ::
+      (PowerExtra typ xl, PowerExtra typ xu) =>
       (MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper) =>
       (Shape.C sh, Class.Floating a) =>
       Integer ->
       Matrix.Quadratic typ xl xu lower upper sh a ->
       Matrix.Quadratic typ xl xu lower upper sh a
    powers1 ::
+      (PowerExtra typ xl, PowerExtra typ xu) =>
       (MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper) =>
       (Shape.C sh, Class.Floating a) =>
       Matrix.Quadratic typ xl xu lower upper sh a ->
       Stream (Matrix.Quadratic typ xl xu lower upper sh a)
 
-instance (xl ~ (), xu ~ ()) => Power Matrix.Scale xl xu where
+instance Power Matrix.Scale where
+   type PowerExtra Matrix.Scale extra = extra ~ ()
    square (Matrix.Scale sh a) = Matrix.Scale sh (a*a)
    power n (Matrix.Scale sh a) = Matrix.Scale sh (a^n)
    powers1 (Matrix.Scale sh a) =
       fmap (Matrix.Scale sh) $ Stream.iterate (*a) a
 
-instance (xl ~ (), xu ~ ()) => Power Matrix.Permutation xl xu where
+instance Power Matrix.Permutation where
+   type PowerExtra Matrix.Permutation extra = extra ~ ()
    square (Matrix.Permutation p) = Matrix.Permutation $ Perm.square p
    power n (Matrix.Permutation p) = Matrix.Permutation $ Perm.power n p
    powers1 (Matrix.Permutation p) =
       fmap Matrix.Permutation $ Stream.iterate (flip Perm.multiplyUnchecked p) p
 
 instance
-   (Layout.Packing pack, Omni.Property property, xl ~ (), xu ~ ()) =>
-      Power (ArrMatrix.Array pack property) xl xu where
+   (Layout.Packing pack, Omni.Property property) =>
+      Power (ArrMatrix.Array pack property) where
+   type PowerExtra (ArrMatrix.Array pack property) extra = extra ~ ()
    square = Multiply.square
    power = Multiply.power
    powers1 = Multiply.powers1
 
 powers ::
-   (Power typ xl xu, MatrixClass.SquareShape typ) =>
+   (Power typ, PowerExtra typ xl, PowerExtra typ xu) =>
+   (MatrixClass.SquareShape typ,
+    MatrixClass.SquareShapeExtra typ xl,
+    MatrixClass.SquareShapeExtra typ xu) =>
    (MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper) =>
    (Shape.C sh, Class.Floating a) =>
    Matrix.Quadratic typ xl xu lower upper sh a ->
diff --git a/src/Numeric/LAPACK/Matrix/Permutation.hs b/src/Numeric/LAPACK/Matrix/Permutation.hs
--- a/src/Numeric/LAPACK/Matrix/Permutation.hs
+++ b/src/Numeric/LAPACK/Matrix/Permutation.hs
@@ -16,14 +16,14 @@
    ) where
 
 import qualified Numeric.LAPACK.Matrix.Array.Private as ArrMatrix
-import qualified Numeric.LAPACK.Matrix.Type as Matrix
+import qualified Numeric.LAPACK.Matrix.Type.Private as Matrix
 import qualified Numeric.LAPACK.Matrix.Shape.Omni as Omni
 import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
 import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
 import qualified Numeric.LAPACK.Matrix.Modifier as Mod
 import qualified Numeric.LAPACK.Permutation.Private as Plain
 import qualified Numeric.LAPACK.Permutation as Perm
-import Numeric.LAPACK.Matrix.Type (Matrix(Permutation))
+import Numeric.LAPACK.Matrix.Type.Private (Matrix(Permutation))
 import Numeric.LAPACK.Vector (Vector)
 
 import qualified Numeric.Netlib.Class as Class
diff --git a/src/Numeric/LAPACK/Matrix/Plain/Format.hs b/src/Numeric/LAPACK/Matrix/Plain/Format.hs
--- a/src/Numeric/LAPACK/Matrix/Plain/Format.hs
+++ b/src/Numeric/LAPACK/Matrix/Plain/Format.hs
@@ -8,10 +8,12 @@
 import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
 import qualified Numeric.LAPACK.Vector as Vector
 import qualified Numeric.LAPACK.Output as Output
+import qualified Numeric.LAPACK.Shape as ExtShape
 import Numeric.LAPACK.Output (Output, formatRow, (<+>))
 import Numeric.LAPACK.Matrix.Layout.Private
          (Order(RowMajor, ColumnMajor), UnaryProxy)
 import Numeric.LAPACK.Matrix.Private (Full)
+import Numeric.LAPACK.Matrix.Extent.Private (Extent)
 import Numeric.LAPACK.Scalar (conjugate)
 import Numeric.LAPACK.Private (caseRealComplexFunc)
 
@@ -31,23 +33,38 @@
 
 import qualified Data.List.Match as Match
 import qualified Data.List.HT as ListHT
-import Data.Functor.Compose (Compose(Compose, getCompose))
 import Data.Foldable (Foldable, fold)
 import Data.List (mapAccumL, transpose)
 import Data.Complex (Complex((:+)))
 import Data.Ix (Ix)
 
+import Control.Applicative ((<$>))
 
+
 deflt :: String
 deflt = "%.4g"
 
 
+data Config =
+   Config {
+      configFormat :: String,
+      configEmpty :: String
+   }
+
+defltConfig :: Config
+defltConfig =
+   Config {
+      configFormat = deflt,
+      configEmpty = ""
+   }
+
+
 class (Shape.C sh) => FormatArray sh where
    {-
    We use constraint @(Class.Floating a)@ and not @(Format a)@
    because it allows us to align the components of complex numbers.
    -}
-   formatArray :: (Class.Floating a, Output out) => String -> Array sh a -> out
+   formatArray :: (Class.Floating a, Output out) => Config -> Array sh a -> out
 
 instance (Integral i) => FormatArray (Shape.ZeroBased i) where
    formatArray = formatVector
@@ -65,49 +82,51 @@
    formatArray = formatVector
 
 instance (FormatArray sh) => FormatArray (Shape.Deferred sh) where
-   formatArray fmt =
-      formatArray fmt . Array.mapShape (\(Shape.Deferred sh) -> sh)
+   formatArray cfg =
+      formatArray cfg . Array.mapShape (\(Shape.Deferred sh) -> sh)
 
 instance (FormatArray sh0, FormatArray sh1) => FormatArray (sh0::+sh1) where
-   formatArray fmt v =
-      formatArray fmt (Vector.takeLeft v)
+   formatArray cfg v =
+      formatArray cfg (Vector.takeLeft v)
       <+>
-      formatArray fmt (Vector.takeRight v)
+      formatArray cfg (Vector.takeRight v)
 
+instance (Shape.C sh) => FormatArray (ExtShape.IntIndexed sh) where
+   formatArray = formatVector
+
 formatVector ::
-   (Shape.C sh, Class.Floating a, Output out) => String -> Array sh a -> out
-formatVector fmt =
-   formatRow . map (Output.text . fold . printfFloating fmt) . Array.toList
+   (Shape.C sh, Class.Floating a, Output out) => Config -> Array sh a -> out
+formatVector cfg =
+   formatRow . map (Output.text . fold . printfFloating cfg) . Array.toList
 
 instance
    (Extent.Measure meas, Extent.C vert, Extent.C horiz,
     Shape.C height, Shape.C width) =>
       FormatArray (Layout.Full meas vert horiz height width) where
-   formatArray = formatFull
+   formatArray cfg = formatAligned (printfFloating cfg) . layoutFull
 
-formatFull ::
+layoutFull ::
    (Extent.Measure meas, Extent.C vert, Extent.C horiz,
-    Shape.C height, Shape.C width, Output out, Class.Floating a) =>
-   String -> Full meas vert horiz height width a -> out
-formatFull fmt m =
+    Shape.C height, Shape.C width, Class.Floating a) =>
+   Full meas vert horiz height width a -> [[a]]
+layoutFull m =
    let Layout.Full order extent = Array.shape m
-   in  formatAligned (printfFloating fmt) $
-       splitRows order (Extent.dimensions extent) $ Array.toList m
+   in  splitRows order (Extent.dimensions extent) $ Array.toList m
 
 instance
    (Eq lower, Extent.Measure meas, Extent.C vert, Extent.C horiz,
     Shape.C height, Shape.C width) =>
       FormatArray (Layout.Split lower meas vert horiz height width) where
-   formatArray = formatHouseholder
+   formatArray cfg =
+      formatSeparateTriangle (printfFloating cfg) . layoutSplit
 
-formatHouseholder ::
+layoutSplit ::
    (Eq lower, Extent.Measure meas, Extent.C vert, Extent.C horiz,
-    Shape.C height, Shape.C width, Class.Floating a, Output out) =>
-   String -> Array (Layout.Split lower meas vert horiz height width) a -> out
-formatHouseholder fmt m =
+    Shape.C height, Shape.C width, Class.Floating a) =>
+   Array (Layout.Split lower meas vert horiz height width) a -> [[a]]
+layoutSplit m =
    let Layout.Split _ order extent = Array.shape m
-   in formatSeparateTriangle (printfFloating fmt) $
-      splitRows order (Extent.dimensions extent) $ Array.toList m
+   in splitRows order (Extent.dimensions extent) $ Array.toList m
 
 instance (Shape.C size) => FormatArray (Layout.Hermitian size) where
    formatArray = formatMirrored conjugate
@@ -118,15 +137,21 @@
 formatMirrored ::
    (Shape.C size, Class.Floating a, Output out) =>
    (a -> a) ->
-   String ->
+   Config ->
    Array (Layout.Mosaic Layout.Packed mirror Shape.Upper size) a ->
    out
-formatMirrored adapt fmt m =
-   let Layout.Mosaic Layout.Packed _mirror Layout.Upper
-            order size
+formatMirrored adapt cfg m =
+   formatSeparateTriangle (printfFloating cfg) $ layoutMirrored adapt m
+
+layoutMirrored ::
+   (Shape.C size, Class.Floating a) =>
+   (a -> a) ->
+   Array (Layout.Mosaic Layout.Packed mirror Shape.Upper size) a ->
+   [[a]]
+layoutMirrored adapt m =
+   let Layout.Mosaic Layout.Packed _mirror Layout.Upper order size
          = Array.shape m
-   in  formatSeparateTriangle (printfFloating fmt) $
-       complementTriangle adapt order (Shape.size size) $ Array.toList m
+   in  complementTriangle adapt order (Shape.size size) $ Array.toList m
 
 complementTriangle ::
    (Class.Floating a) => (a -> a) -> Order -> Int -> [a] -> [[a]]
@@ -145,27 +170,25 @@
 
 formatDiagonal ::
    (Shape.C size, Class.Floating a, Output out) =>
-   String -> Order -> size -> [a] -> out
-formatDiagonal fmt order size xs =
+   Config -> Order -> size -> [a] -> out
+formatDiagonal cfg order size xs =
    let n0 = Unary.unary TypeNum.u0
-   in formatAligned (printfFloatingMaybe fmt) $
+   in formatAligned (printfFloatingMaybe cfg) $
       padBanded (n0,n0) order (size,size) xs
 
 
 instance
    (Layout.UpLo uplo, Shape.C size) =>
       FormatArray (Layout.Triangular uplo size) where
-   formatArray = formatTriangular
+   formatArray cfg = formatAligned (printfFloatingMaybe cfg) . layoutTriangular
 
-formatTriangular ::
-   (Layout.UpLo uplo, Shape.C size, Class.Floating a, Output out) =>
-   String -> Array (Layout.Triangular uplo size) a -> out
-formatTriangular fmt m =
-   let Layout.Mosaic Layout.Packed Layout.NoMirror
-            uplo order size
+layoutTriangular ::
+   (Layout.UpLo uplo, Shape.C size, Class.Floating a) =>
+   Array (Layout.Triangular uplo size) a -> [[Maybe a]]
+layoutTriangular m =
+   let Layout.Mosaic Layout.Packed Layout.NoMirror uplo order size
          = Array.shape m
-   in formatAligned (printfFloatingMaybe fmt) $
-      padTriangle uplo order (Shape.size size) $ Array.toList m
+   in padTriangle uplo order (Shape.size size) $ Array.toList m
 
 padTriangle ::
    Layout.UpLoSingleton uplo -> Order -> Int -> [a] -> [[Maybe a]]
@@ -206,12 +229,17 @@
     Extent.Measure meas, Extent.C vert, Extent.C horiz,
     Shape.C height, Shape.C width) =>
       FormatArray (Layout.Banded sub super meas vert horiz height width) where
-   formatArray fmt m =
-      let Layout.Banded offDiag order extent = Array.shape m
-      in  formatAligned (printfFloatingMaybe fmt) $
-          padBanded offDiag order (Extent.dimensions extent) $
-          Array.toList m
+   formatArray cfg = formatAligned (printfFloatingMaybe cfg) . layoutBanded
 
+layoutBanded ::
+   (Unary.Natural sub, Unary.Natural super,
+    Extent.Measure meas, Extent.C vert, Extent.C horiz,
+    Shape.C height, Shape.C width, Class.Floating a) =>
+   Array (Layout.Banded sub super meas vert horiz height width) a -> [[Maybe a]]
+layoutBanded m =
+   let Layout.Banded offDiag order extent = Array.shape m
+   in  padBanded offDiag order (Extent.dimensions extent) $ Array.toList m
+
 padBanded ::
    (Shape.C height, Shape.C width, Unary.Natural sub, Unary.Natural super) =>
    (UnaryProxy sub, UnaryProxy super) -> Order ->
@@ -240,11 +268,16 @@
 instance
    (Unary.Natural offDiag, Shape.C size) =>
       FormatArray (Layout.BandedHermitian offDiag size) where
-   formatArray fmt m =
-      let Layout.BandedHermitian offDiag order size = Array.shape m
-      in  formatSeparateTriangle (printfFloatingMaybe fmt) $
-          padBandedHermitian offDiag order size $ Array.toList m
+   formatArray cfg =
+      formatSeparateTriangle (printfFloatingMaybe cfg) . layoutBandedHermitian
 
+layoutBandedHermitian ::
+   (Unary.Natural offDiag, Shape.C sh, Class.Floating a) =>
+   Array (Layout.BandedHermitian offDiag sh) a -> [[Maybe a]]
+layoutBandedHermitian m =
+   let Layout.BandedHermitian offDiag order size = Array.shape m
+   in padBandedHermitian offDiag order size $ Array.toList m
+
 padBandedHermitian ::
    (Unary.Natural offDiag, Shape.C size, Class.Floating a) =>
    UnaryProxy offDiag -> Order -> size -> [a] -> [[Maybe a]]
@@ -285,17 +318,78 @@
 
 
 formatAligned ::
-   (Functor f, Foldable f) => Output out => (a -> f String) -> [[a]] -> out
+   (Functor f, Foldable f, Output out) =>
+   (a -> f String) -> [[a]] -> out
 formatAligned printFmt =
    Output.formatAligned . map (map (fmap Output.text . printFmt))
 
 formatSeparateTriangle ::
-   (Functor f, Foldable f) => Output out => (a -> f String) -> [[a]] -> out
+   (Functor f, Foldable f, Output out) =>
+   (a -> f String) -> [[a]] -> out
 formatSeparateTriangle printFmt =
    Output.formatSeparateTriangle . map (map (fmap Output.text . printFmt))
 
 
+arrayFromList2 ::
+   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+   (Shape.C height, Shape.C width) =>
+   Extent meas vert horiz height width -> [[a]] ->
+   BoxedArray.Array (height, width) (Output.Separator, Maybe (Output.Style, a))
+arrayFromList2 extent =
+   incompleteArrayFromList2 extent . map (map Just)
 
+incompleteArrayFromList2 ::
+   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+   (Shape.C height, Shape.C width) =>
+   Extent meas vert horiz height width -> [[Maybe a]] ->
+   BoxedArray.Array (height, width) (Output.Separator, Maybe (Output.Style, a))
+incompleteArrayFromList2 extent =
+   let (height, width) = Extent.dimensions extent in
+   let n = Shape.size width in
+   BoxedArray.fromList (height, width) . concat .
+   map (map ((,) Output.Space) .
+        ListHT.padRight Nothing n . map (fmap ((,) Output.Stored)))
+
+splitArrayFromList2 ::
+   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+   (Shape.C height, Shape.C width) =>
+   Extent meas vert horiz height width -> [[a]] ->
+   BoxedArray.Array (height, width) (Output.Separator, Maybe (Output.Style, a))
+splitArrayFromList2 extent =
+   incompleteSplitArrayFromList2 extent . map (map Just)
+
+incompleteSplitArrayFromList2 ::
+   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+   (Shape.C height, Shape.C width) =>
+   Extent meas vert horiz height width -> [[Maybe a]] ->
+   BoxedArray.Array (height, width) (Output.Separator, Maybe (Output.Style, a))
+incompleteSplitArrayFromList2 extent =
+   let (height, width) = Extent.dimensions extent in
+   let n = Shape.size width in
+   BoxedArray.fromList (height, width) . concat .
+   Output.decorateTriangle
+      (\seps styles row ->
+         zip seps $ ListHT.padRight Nothing n $
+         zipWith (\s r -> (,) s <$> r) styles row)
+
+-- ToDo: could be in BoxedArray
+toRows ::
+   (Shape.C height, Shape.C width) =>
+   BoxedArray.Array (height,width) a -> [[a]]
+toRows arr =
+   case BoxedArray.toList arr of
+      [] -> replicate (Shape.size $ fst $ BoxedArray.shape arr) []
+      xs -> ListHT.sliceVertical (Shape.size $ snd $ BoxedArray.shape arr) xs
+
+toColumns ::
+   (Shape.C height, Shape.C width) =>
+   BoxedArray.Array (height,width) a -> [[a]]
+toColumns arr =
+   ListHT.sliceHorizontal (Shape.size $ snd $ BoxedArray.shape arr) $
+   BoxedArray.toList arr
+
+
+
 data TupleShape a = TupleShape
 
 instance (Class.Floating a) => Shape.C (TupleShape a) where
@@ -309,8 +403,8 @@
 
 newtype ToTuple a = ToTuple {getToTuple :: a -> Tuple a String}
 
-printfFloating :: (Class.Floating a) => String -> a -> Tuple a String
-printfFloating fmt =
+printfFloatingPlain :: (Class.Floating a) => String -> a -> Tuple a String
+printfFloatingPlain fmt =
    getToTuple $
    Class.switchFloating
       (ToTuple $ fillTuple . printf fmt)
@@ -318,22 +412,18 @@
       (ToTuple $ printfComplex fmt)
       (ToTuple $ printfComplex fmt)
 
+printfFloating :: (Class.Floating a) => Config -> a -> Tuple a String
+printfFloating cfg = printfFloatingPlain (configFormat cfg)
+
 printfFloatingMaybe ::
-   (Class.Floating a) => String -> Maybe a -> Tuple a String
-printfFloatingMaybe = maybe (fillTuple "") . printfFloating
+   (Class.Floating a) => Config -> Maybe a -> Tuple a String
+printfFloatingMaybe cfg =
+   maybe (fillTuple $ configEmpty cfg) (printfFloating cfg)
 
 printfComplex ::
-   (Class.Real a) => String -> Complex a -> Tuple (Complex a) String
-printfComplex fmt =
-   getToTuple $ getCompose $
-   Class.switchReal
-      (Compose $ ToTuple $ printfComplexAux fmt)
-      (Compose $ ToTuple $ printfComplexAux fmt)
-
-printfComplexAux ::
    (PrintfArg a, Class.Real a) =>
    String -> Complex a -> Tuple (Complex a) String
-printfComplexAux fmt (r:+i) =
+printfComplex fmt (r:+i) =
    if i<0 || isNegativeZero i
      then complexTuple (printf (fmt ++ "-") r) (printf (fmt ++ "i") (-i))
      else complexTuple (printf (fmt ++ "+") r) (printf (fmt ++ "i") i)
diff --git a/src/Numeric/LAPACK/Matrix/Quadratic.hs b/src/Numeric/LAPACK/Matrix/Quadratic.hs
--- a/src/Numeric/LAPACK/Matrix/Quadratic.hs
+++ b/src/Numeric/LAPACK/Matrix/Quadratic.hs
@@ -26,7 +26,7 @@
 import qualified Numeric.LAPACK.Matrix.Square.Basic as Square
 import qualified Numeric.LAPACK.Matrix.Basic as Full
 
-import qualified Numeric.LAPACK.Matrix.Type as Matrix
+import qualified Numeric.LAPACK.Matrix.Type.Private as Matrix
 import qualified Numeric.LAPACK.Matrix.Class as MatrixClass
 import qualified Numeric.LAPACK.Matrix.Array.Basic as OmniMatrix
 import qualified Numeric.LAPACK.Matrix.Array.Private as ArrMatrix
diff --git a/src/Numeric/LAPACK/Matrix/Shape/Omni.hs b/src/Numeric/LAPACK/Matrix/Shape/Omni.hs
--- a/src/Numeric/LAPACK/Matrix/Shape/Omni.hs
+++ b/src/Numeric/LAPACK/Matrix/Shape/Omni.hs
@@ -43,6 +43,7 @@
    mapHeight, mapWidth, mapSquareSize,
    order,
    transpose,
+   Cons(cons),
    Plain,
    ToPlain, toPlain,
    FromPlain, fromPlain,
diff --git a/src/Numeric/LAPACK/Matrix/Special.hs b/src/Numeric/LAPACK/Matrix/Special.hs
--- a/src/Numeric/LAPACK/Matrix/Special.hs
+++ b/src/Numeric/LAPACK/Matrix/Special.hs
@@ -4,24 +4,18 @@
    ) where
 
 import qualified Numeric.LAPACK.Matrix.Inverse as Inverse
-import qualified Numeric.LAPACK.Matrix.Type as Matrix
+import qualified Numeric.LAPACK.Matrix.Type.Private as Matrix
 import qualified Numeric.LAPACK.Matrix.Shape.Omni as Omni
 import Numeric.LAPACK.Matrix.Layout.Private (Empty, Filled)
 
-import Data.Tuple.HT (mapPair)
 
-
 type Scale sh = Matrix.Quadratic Matrix.Scale () () Empty Empty sh
 type Inverse typ lower upper sh =
       Matrix.Quadratic (Inverse.Inverse typ) lower upper Filled Filled sh
 
 inverse ::
-   (Omni.Strip lower, Inverse.Fill lower ~ lowerf,
-    Omni.Strip upper, Inverse.Fill upper ~ upperf) =>
+   (Omni.PowerStrip lower, Omni.PowerStrip upper) =>
    Matrix.QuadraticMeas typ xl xu upper lower meas width height a ->
    Matrix.QuadraticMeas (Inverse.Inverse typ) (xl,lower) (xu,upper)
-      lowerf upperf meas height width a
-inverse a =
-   case mapPair (Inverse.filledPowerStripFact, Inverse.filledPowerStripFact) $
-        Matrix.strips a of
-      (Inverse.PowerStripFact, Inverse.PowerStripFact) -> Inverse.Inverse a
+      lower upper meas height width a
+inverse a = Inverse.Inverse a
diff --git a/src/Numeric/LAPACK/Matrix/Superscript.hs b/src/Numeric/LAPACK/Matrix/Superscript.hs
--- a/src/Numeric/LAPACK/Matrix/Superscript.hs
+++ b/src/Numeric/LAPACK/Matrix/Superscript.hs
@@ -1,5 +1,6 @@
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE GADTs #-}
+{-# LANGUAGE ConstraintKinds #-}
 {-# LANGUAGE EmptyDataDecls #-}
 {- |
 This module provides an infix operator and some data constructors
@@ -67,7 +68,8 @@
       typA xlA xuA lowerA upperA measA vertA horizA heightA widthA
       typB xlB xuB lowerB upperB measB vertB horizB heightB widthB a where
          T ::
-            (Matrix.Transpose typ) =>
+            (Matrix.Transpose typ,
+             Matrix.TransposeExtra typ xl, Matrix.TransposeExtra typ xu) =>
             Exponent Transpose
                typ xl xu lower upper meas vert horiz height width
                typ xu xl upper lower meas horiz vert width height a
@@ -80,7 +82,8 @@
       typA xlA xuA lowerA upperA measA vertA horizA heightA widthA
       typB xlB xuB lowerB upperB measB vertB horizB heightB widthB a where
          A ::
-            (Matrix.Transpose typ, Matrix.Complex typ) =>
+            (Matrix.Transpose typ, Matrix.Complex typ,
+             Matrix.TransposeExtra typ xl, Matrix.TransposeExtra typ xu) =>
             Exponent Adjoint
                typ xl xu lower upper meas vert horiz height width
                typ xu xl upper lower meas horiz vert width height a
@@ -107,7 +110,8 @@
       typB xlB xuB lowerB upperB measB vertB horizB heightB widthB a where
          -- We could relax lowerB and upperB using Filled type function
          Inv ::
-            (Matrix.Inverse typ xl xu,
+            (Matrix.Inverse typ,
+             Matrix.InverseExtra typ xl, Matrix.InverseExtra typ xu,
              Omni.PowerStrip lower, Omni.PowerStrip upper) =>
             Exponent Inverse
                typ xl xu lower upper meas Small Small height width
@@ -136,7 +140,8 @@
       typB xlB xuB lowerB upperB measB vertB horizB heightB widthB a where
          -- We could relax lowerB and upperB using Filled type function
          Exp ::
-            (Matrix.Power typ xl xu,
+            (Matrix.Power typ,
+             Matrix.PowerExtra typ xl, Matrix.PowerExtra typ xu,
              Omni.PowerStrip lower, Omni.PowerStrip upper) =>
             Integer ->
             Exponent Power
diff --git a/src/Numeric/LAPACK/Matrix/Symmetric/Unified.hs b/src/Numeric/LAPACK/Matrix/Symmetric/Unified.hs
--- a/src/Numeric/LAPACK/Matrix/Symmetric/Unified.hs
+++ b/src/Numeric/LAPACK/Matrix/Symmetric/Unified.hs
@@ -610,7 +610,7 @@
 
 solve ::
    (Mirror mirror, Extent.Measure meas, Extent.C vert, Extent.C horiz,
-    Shape.C width, Shape.C height, Eq height, Class.Floating a) =>
+    Eq height, Shape.C height, Shape.C width, Class.Floating a) =>
    Mosaic pack mirror Shape.Upper height a ->
    Full meas vert horiz height width a ->
    Full meas vert horiz height width a
diff --git a/src/Numeric/LAPACK/Matrix/Triangular.hs b/src/Numeric/LAPACK/Matrix/Triangular.hs
--- a/src/Numeric/LAPACK/Matrix/Triangular.hs
+++ b/src/Numeric/LAPACK/Matrix/Triangular.hs
@@ -64,7 +64,7 @@
 import qualified Numeric.LAPACK.Matrix.Array.Unpacked as ArrUnpacked
 import qualified Numeric.LAPACK.Matrix.Array.Basic as OmniMatrix
 import qualified Numeric.LAPACK.Matrix.Array.Private as ArrMatrix
-import qualified Numeric.LAPACK.Matrix.Type as Matrix
+import qualified Numeric.LAPACK.Matrix.Type.Private as Matrix
 import qualified Numeric.LAPACK.Matrix.Class as MatrixClass
 import qualified Numeric.LAPACK.Matrix.Shape.Omni as Omni
 import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
diff --git a/src/Numeric/LAPACK/Matrix/Type.hs b/src/Numeric/LAPACK/Matrix/Type.hs
--- a/src/Numeric/LAPACK/Matrix/Type.hs
+++ b/src/Numeric/LAPACK/Matrix/Type.hs
@@ -1,337 +1,13 @@
-{-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE GADTs #-}
-{-# LANGUAGE EmptyDataDecls #-}
-{-# LANGUAGE FlexibleInstances #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE StandaloneDeriving #-}
-module Numeric.LAPACK.Matrix.Type where
-
-import qualified Numeric.LAPACK.Matrix.Plain.Format as ArrFormat
-import qualified Numeric.LAPACK.Output as Output
-import qualified Numeric.LAPACK.Permutation.Private as Perm
-import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
-import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
-import qualified Numeric.LAPACK.Matrix.Shape.Omni as Omni
-import qualified Numeric.LAPACK.Matrix.Extent.Strict as ExtentStrict
-import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
-import Numeric.LAPACK.Matrix.Layout.Private (Empty, Filled)
-import Numeric.LAPACK.Matrix.Extent.Private (Extent, Shape, Small)
-import Numeric.LAPACK.Output (Output)
-
-import qualified Numeric.Netlib.Class as Class
-
-import qualified Hyper
-
-import qualified Control.DeepSeq as DeepSeq
-
-import qualified Data.Array.Comfort.Shape as Shape
-
-import Data.Function.HT (Id)
-import Data.Monoid (Monoid, mempty, mappend)
-import Data.Semigroup (Semigroup, (<>))
-
-
-
-data family
-   Matrix typ extraLower extraUpper lower upper meas vert horiz height width a
-
-type Quadratic typ extraLower extraUpper lower upper sh =
-      QuadraticMeas typ extraLower extraUpper lower upper Shape sh sh
-type QuadraticMeas typ extraLower extraUpper lower upper meas =
-      Matrix typ extraLower extraUpper lower upper meas Small Small
-
-
-asQuadratic ::
-   Id (QuadraticMeas typ extraLower extraUpper lower upper meas height width a)
-asQuadratic = id
-
-
-data Product fuse
-data instance
-   Matrix (Product fuse) xl xu lower upper meas vert horiz height width a where
-      Product ::
-         (Omni.MultipliedBands lowerA lowerB ~ lowerC,
-          Omni.MultipliedBands lowerB lowerA ~ lowerC,
-          Omni.MultipliedBands upperA upperB ~ upperC,
-          Omni.MultipliedBands upperB upperA ~ upperC) =>
-         Matrix typA xlA xuA lowerA upperA meas vert horiz height fuse a ->
-         Matrix typB xlB xuB lowerB upperB meas vert horiz fuse width a ->
-         Matrix (Product fuse)
-            (typA,lowerA,lowerB,xlA,xlB) (typB,upperB,upperA,xuB,xuA)
-            lowerC upperC meas vert horiz height width a
-
-
-data Scale
-data instance
-   Matrix Scale xl xu lower upper meas vert horiz height width a where
-      Scale :: sh -> a -> Quadratic Scale () () Empty Empty sh a
-
-deriving instance
-   (Shape.C height, Show height, Show a) =>
-   Show (Matrix Scale xl xu lower upper meas vert horiz height width a)
-
-
-data Identity
-data instance
-   Matrix Identity xl xu lower upper meas vert horiz height width a where
-      Identity ::
-         (Extent.Measure meas) =>
-         Extent meas Small Small height width ->
-         QuadraticMeas Identity () () Empty Empty meas height width a
-
-
-data Permutation
-data instance
-   Matrix Permutation xl xu lower upper meas vert horiz height width a where
-   Permutation ::
-      Perm.Permutation sh -> Quadratic Permutation () () lower upper sh a
-
-deriving instance
-   (Shape.C height, Show height) =>
-   Show (Matrix Permutation xl xu lower upper meas vert horiz height width a)
-
-deriving instance
-   (Shape.C height, Eq height) =>
-   Eq (Matrix Permutation xl xu lower upper meas vert horiz height width a)
-
-
-instance
-   (NFData typ,
-    Extent.Measure meas, Extent.C vert, Extent.C horiz,
-    DeepSeq.NFData height, DeepSeq.NFData width, DeepSeq.NFData a) =>
-   DeepSeq.NFData (Matrix typ xl xu lower upper meas vert horiz height width a)
-      where
-   rnf = rnf
-
-class NFData typ where
-   rnf ::
-      (Extent.Measure meas, Extent.C vert, Extent.C horiz,
-       DeepSeq.NFData height, DeepSeq.NFData width, DeepSeq.NFData a) =>
-      Matrix typ xl xu lower upper meas vert horiz height width a -> ()
-
-
-
-instance
-   (FormatMatrix typ, Extent.Measure meas, Extent.C vert, Extent.C horiz,
-    Shape.C width, Shape.C height, Class.Floating a) =>
-      Hyper.Display
-         (Matrix typ xl xu lower upper meas vert horiz height width a) where
-   display = Output.hyper . formatMatrix ArrFormat.deflt
-
-
-class FormatMatrix typ where
-   {-
-   We use constraint @(Class.Floating a)@ and not @(Format a)@
-   because it allows us to align the components of complex numbers.
-   -}
-   formatMatrix ::
-      (Extent.Measure meas, Extent.C vert, Extent.C horiz,
-       Shape.C width, Shape.C height, Class.Floating a, Output out) =>
-      String ->
-      Matrix typ xl xu lower upper meas vert horiz height width a ->
-      out
-
-instance FormatMatrix Scale where
-   formatMatrix fmt (Scale shape a) =
-      ArrFormat.formatDiagonal fmt Layout.RowMajor shape $
-      replicate (Shape.size shape) a
-
-instance FormatMatrix Permutation where
-   formatMatrix _fmt (Permutation perm) = Perm.format perm
-
-
-
-instance
-   (MultiplySame typ xl xu,
-    MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper,
-    Extent.Measure meas, Extent.C vert, Extent.C horiz,
-    Shape.C height, Eq height, height ~ width, Class.Floating a) =>
-      Semigroup (Matrix typ xl xu lower upper meas vert horiz height width a)
-         where
-   (<>) = multiplySame
-
-class (Box typ) => MultiplySame typ xl xu where
-   multiplySame ::
-      (matrix ~ Matrix typ xl xu lower upper meas vert horiz sh sh a,
-       MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper,
-       Extent.Measure meas, Extent.C vert, Extent.C horiz,
-       Shape.C sh, Eq sh, Class.Floating a) =>
-      matrix -> matrix -> matrix
-
-instance (xl ~ (), xu ~ ()) => MultiplySame Scale xl xu where
-   multiplySame =
-      scaleWithCheck "Scale.multiplySame" height
-         (\a (Scale shape b) -> Scale shape $ a*b)
-
-instance (xl ~ (), xu ~ ()) => MultiplySame Permutation xl xu where
-   multiplySame (Permutation a) (Permutation b) =
-      Permutation $ Perm.multiply b a
-
-
-instance
-   (MultiplySame typ xl xu, StaticIdentity typ xl xu lower upper,
-    MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper,
-    meas ~ Shape, vert ~ Small, horiz ~ Small,
-    Shape.Static height, Eq height, height ~ width, Class.Floating a) =>
-      Monoid (Matrix typ xl xu lower upper meas vert horiz height width a) where
-   mappend = (<>)
-   mempty = staticIdentity
-
-class StaticIdentity typ xl xu lower upper where
-   staticIdentity ::
-      (Shape.Static sh, Class.Floating a) =>
-      Quadratic typ xl xu lower upper sh a
-
-instance
-   (xl ~ (), xu ~ (), lower ~ Empty, upper ~ Empty) =>
-      StaticIdentity Scale xl xu lower upper where
-   staticIdentity = Scale Shape.static 1
-
-instance
-   (xl ~ (), xu ~ (), lower ~ Filled, upper ~ Filled) =>
-      StaticIdentity Permutation xl xu lower upper where
-   staticIdentity = Permutation $ Perm.identity Shape.static
-
-
-scaleWithCheck :: (Eq shape) =>
-   String -> (b -> shape) ->
-   (a -> b -> c) ->
-   Matrix Scale xl xu lower upper meas vert horiz shape shape a ->
-   b -> c
-scaleWithCheck name getSize f (Scale shape a) b =
-   if shape == getSize b
-      then f a b
-      else error $ name ++ ": dimensions mismatch"
-
-
-class Box typ where
-   extent ::
-      (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
-      Matrix typ xl xu lower upper meas vert horiz height width a ->
-      Extent.Extent meas vert horiz height width
-   height ::
-      (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
-      Matrix typ xl xu lower upper meas vert horiz height width a -> height
-   height = Extent.height . extent
-   width ::
-      (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
-      Matrix typ xl xu lower upper meas vert horiz height width a -> width
-   width = Extent.width . extent
-
-instance Box Scale where
-   extent (Scale shape _) = Extent.square shape
-   height (Scale shape _) = shape
-   width (Scale shape _) = shape
-
-instance Box Identity where
-   extent (Identity extent_) = extent_
-
-instance Box Permutation where
-   extent (Permutation perm) = Extent.square $ Perm.size perm
-   height (Permutation perm) = Perm.size perm
-   width (Permutation perm) = Perm.size perm
-
-{- ToDo: requires parameters xl and xu for Box class
-
-instance (Eq fuse) => Box (Product fuse) where
-   extent (Product a b) =
-      fromMaybe (error "Matrix.Product: shapes mismatch") $
-      Extent.fuse (extent a) (extent b)
--}
-
-squareSize :: (Box typ) => Quadratic typ xl xu lower upper sh a -> sh
-squareSize = height
-
-indices ::
-   (Box typ, Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
-   (Shape.Indexed height, Shape.Indexed width) =>
-   Matrix typ xl xu lower upper meas vert horiz height width a ->
-   [(Shape.Index height, Shape.Index width)]
-indices sh = Shape.indices (height sh, width sh)
-
-
-class (Box typ) => ToQuadratic typ where
-   heightToQuadratic ::
-      (Extent.Measure meas) =>
-      QuadraticMeas typ xl xu lower upper meas height width a ->
-      Quadratic typ xl xu lower upper height a
-   widthToQuadratic ::
-      (Extent.Measure meas) =>
-      QuadraticMeas typ xl xu lower upper meas height width a ->
-      Quadratic typ xl xu lower upper width a
-
-instance ToQuadratic Scale where
-   heightToQuadratic (Scale shape a) = Scale shape a
-   widthToQuadratic (Scale shape a) = Scale shape a
-
-instance ToQuadratic Identity where
-   heightToQuadratic (Identity extent_) =
-      Identity $ Extent.square $ Extent.height extent_
-   widthToQuadratic (Identity extent_) =
-      Identity $ Extent.square $ Extent.width extent_
-
-instance ToQuadratic Permutation where
-   heightToQuadratic (Permutation perm) = Permutation perm
-   widthToQuadratic (Permutation perm) = Permutation perm
-
-
-class (Box typ) => MapExtent typ xl xu lower upper where
-   mapExtent ::
-      (Extent.Measure measA, Extent.C vertA, Extent.C horizA) =>
-      (Extent.Measure measB, Extent.C vertB, Extent.C horizB) =>
-      ExtentStrict.Map measA vertA horizA measB vertB horizB height width ->
-      Matrix typ xl xu lower upper measA vertA horizA height width a ->
-      Matrix typ xl xu lower upper measB vertB horizB height width a
-
-
-class (Box typ) => Transpose typ where
-   transpose ::
-      (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
-      (Shape.C width, Shape.C height, Class.Floating a) =>
-      Matrix typ xl xu lower upper meas vert horiz height width a ->
-      Matrix typ xu xl upper lower meas horiz vert width height a
-
-instance Transpose Scale where
-   transpose (Scale shape a) = Scale shape a
-
-instance Transpose Identity where
-   transpose (Identity extent_) = Identity $ Extent.transpose extent_
-
-instance Transpose Permutation where
-   transpose (Permutation perm) = Permutation $ Perm.transpose perm
-
-{-
-instance (Shape.C fuse, Eq fuse) => Transpose (Product fuse) where
-   transpose (Product a b) = Product (transpose b) (transpose a)
--}
-
-
-swapMultiply ::
-   (Transpose typA, Transpose typB) =>
-   (Extent.Measure measA, Extent.C vertA, Extent.C horizA) =>
-   (Extent.Measure measB, Extent.C vertB, Extent.C horizB) =>
-   (Shape.C heightA, Shape.C widthA) =>
-   (Shape.C heightB, Shape.C widthB) =>
-   (Class.Floating a) =>
-   (matrix ->
-    Matrix typA xuA xlA upperA lowerA measA horizA vertA widthA heightA a ->
-    Matrix typB xuB xlB upperB lowerB measB horizB vertB widthB heightB a) ->
-   Matrix typA xlA xuA lowerA upperA measA vertA horizA heightA widthA a ->
-   matrix ->
-   Matrix typB xlB xuB lowerB upperB measB vertB horizB heightB widthB a
-swapMultiply multiplyTrans a b =
-   transpose $ multiplyTrans b $ transpose a
-
-powerStrips ::
-   (MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper) =>
-   Matrix typ xl xu lower upper meas vert horiz height width a ->
-   (MatrixShape.PowerStripSingleton lower,
-    MatrixShape.PowerStripSingleton upper)
-powerStrips _ =
-   (MatrixShape.powerStripSingleton, MatrixShape.powerStripSingleton)
+module Numeric.LAPACK.Matrix.Type (
+   Type.Scale,
+   Type.Identity,
+   Type.Permutation,
+   Type.Product,
+   Inverse.Inverse,
+   Wrapper.FillStrips,
+   Wrapper.MapExtent,
+   ) where
 
-strips ::
-   (MatrixShape.Strip lower, MatrixShape.Strip upper) =>
-   Matrix typ xl xu lower upper meas vert horiz height width a ->
-   (MatrixShape.StripSingleton lower, MatrixShape.StripSingleton upper)
-strips _ = (MatrixShape.stripSingleton, MatrixShape.stripSingleton)
+import qualified Numeric.LAPACK.Matrix.Type.Private as Type
+import qualified Numeric.LAPACK.Matrix.Inverse as Inverse
+import qualified Numeric.LAPACK.Matrix.Wrapper as Wrapper
diff --git a/src/Numeric/LAPACK/Matrix/Type/Private.hs b/src/Numeric/LAPACK/Matrix/Type/Private.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Matrix/Type/Private.hs
@@ -0,0 +1,469 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE EmptyDataDecls #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE ConstraintKinds #-}
+{-# LANGUAGE UndecidableInstances #-}
+{-# LANGUAGE StandaloneDeriving #-}
+module Numeric.LAPACK.Matrix.Type.Private where
+
+import qualified Numeric.LAPACK.Matrix.Plain.Format as ArrFormat
+import qualified Numeric.LAPACK.Output as Output
+import qualified Numeric.LAPACK.Permutation.Private as Perm
+import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
+import qualified Numeric.LAPACK.Matrix.Shape.Omni as Omni
+import qualified Numeric.LAPACK.Matrix.Extent.Strict as ExtentStrict
+import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
+import Numeric.LAPACK.Matrix.Layout.Private (Empty, Filled)
+import Numeric.LAPACK.Matrix.Extent.Private (Extent, Shape, Small)
+import Numeric.LAPACK.Output (Output)
+
+import qualified Numeric.Netlib.Class as Class
+
+import qualified Hyper
+
+import qualified Control.DeepSeq as DeepSeq
+import Control.Applicative ((<$>))
+
+import qualified Data.Array.Comfort.Boxed as BoxedArray
+import qualified Data.Array.Comfort.Shape as Shape
+
+import qualified Data.Foldable as Fold
+import Data.Function.HT (Id)
+import Data.Monoid (Monoid, mempty, mappend)
+import Data.Semigroup (Semigroup, (<>))
+import Data.Foldable (Foldable)
+import Data.Maybe (fromMaybe)
+import Data.Tuple.HT (mapSnd)
+
+import GHC.Exts (Constraint)
+
+
+
+data family
+   Matrix typ extraLower extraUpper lower upper meas vert horiz height width a
+
+type Quadratic typ extraLower extraUpper lower upper sh =
+      QuadraticMeas typ extraLower extraUpper lower upper Shape sh sh
+type QuadraticMeas typ extraLower extraUpper lower upper meas =
+      Matrix typ extraLower extraUpper lower upper meas Small Small
+
+
+asQuadratic ::
+   Id (QuadraticMeas typ extraLower extraUpper lower upper meas height width a)
+asQuadratic = id
+
+
+data Product fuse
+data instance
+   Matrix (Product fuse) xl xu lower upper meas vert horiz height width a where
+      Product ::
+         (Omni.MultipliedBands lowerA lowerB ~ lowerC,
+          Omni.MultipliedBands lowerB lowerA ~ lowerC,
+          Omni.MultipliedBands upperA upperB ~ upperC,
+          Omni.MultipliedBands upperB upperA ~ upperC) =>
+         Matrix typA xlA xuA lowerA upperA meas vert horiz height fuse a ->
+         Matrix typB xlB xuB lowerB upperB meas vert horiz fuse width a ->
+         Matrix (Product fuse)
+            (typA,xlA,xuA,lowerA,upperA) (typB,xuB,xlB,upperB,lowerB)
+            lowerC upperC meas vert horiz height width a
+
+type family ProductType   extra
+type family ProductExtraL extra
+type family ProductExtraU extra
+type family ProductLower  extra
+type family ProductUpper  extra
+type instance ProductType   (typ,xl,xu,lower,upper) = typ
+type instance ProductExtraL (typ,xl,xu,lower,upper) = xl
+type instance ProductExtraU (typ,xl,xu,lower,upper) = xu
+type instance ProductLower  (typ,xl,xu,lower,upper) = lower
+type instance ProductUpper  (typ,xl,xu,lower,upper) = upper
+
+
+data Scale
+data instance
+   Matrix Scale xl xu lower upper meas vert horiz height width a where
+      Scale :: sh -> a -> Quadratic Scale () () Empty Empty sh a
+
+deriving instance
+   (Shape.C height, Show height, Show a) =>
+   Show (Matrix Scale xl xu lower upper meas vert horiz height width a)
+
+
+data Identity
+data instance
+   Matrix Identity xl xu lower upper meas vert horiz height width a where
+      Identity ::
+         (Extent.Measure meas) =>
+         Extent meas Small Small height width ->
+         QuadraticMeas Identity () () Empty Empty meas height width a
+
+
+data Permutation
+data instance
+   Matrix Permutation xl xu lower upper meas vert horiz height width a where
+   Permutation ::
+      Perm.Permutation sh -> Quadratic Permutation () () lower upper sh a
+
+deriving instance
+   (Shape.C height, Show height) =>
+   Show (Matrix Permutation xl xu lower upper meas vert horiz height width a)
+
+deriving instance
+   (Shape.C height, Eq height) =>
+   Eq (Matrix Permutation xl xu lower upper meas vert horiz height width a)
+
+
+instance
+   (NFData typ,
+    Extent.Measure meas, Extent.C vert, Extent.C horiz,
+    DeepSeq.NFData height, DeepSeq.NFData width, DeepSeq.NFData a) =>
+   DeepSeq.NFData (Matrix typ xl xu lower upper meas vert horiz height width a)
+      where
+   rnf = rnf
+
+class NFData typ where
+   rnf ::
+      (Extent.Measure meas, Extent.C vert, Extent.C horiz,
+       DeepSeq.NFData height, DeepSeq.NFData width, DeepSeq.NFData a) =>
+      Matrix typ xl xu lower upper meas vert horiz height width a -> ()
+
+
+
+instance
+   (Format typ, FormatExtra typ xl, FormatExtra typ xu,
+    Extent.Measure meas, Extent.C vert, Extent.C horiz,
+    Shape.C height, Shape.C width, Class.Floating a) =>
+      Hyper.Display
+         (Matrix typ xl xu lower upper meas vert horiz height width a) where
+   display = Output.hyper . format ArrFormat.defltConfig
+
+
+class Format typ where
+   type FormatExtra typ extra :: Constraint
+   format ::
+      (FormatExtra typ xl, FormatExtra typ xu,
+       Extent.Measure meas, Extent.C vert, Extent.C horiz,
+       Shape.C height, Shape.C width, Class.Floating a, Output out) =>
+      ArrFormat.Config ->
+      Matrix typ xl xu lower upper meas vert horiz height width a ->
+      out
+
+{- |
+Default implementation of 'format'.
+Some matrices need more than one array for display,
+e.g. @Householder@ and @LowerUpper@.
+'Layout' class is still needed for @Block@ matrices.
+-}
+formatWithLayout ::
+   (Layout typ, LayoutExtra typ xl, LayoutExtra typ xu,
+    Extent.Measure meas, Extent.C vert, Extent.C horiz,
+    Shape.C height, Shape.C width, Class.Floating a, Output out) =>
+   ArrFormat.Config ->
+   Matrix typ xl xu lower upper meas vert horiz height width a ->
+   out
+formatWithLayout fmt =
+   Output.formatTable . map (concatMap (uncurry attachSeparators)) .
+   ArrFormat.toRows . fmap (mapSnd (formatCell fmt)) . layout
+
+attachSeparators :: (Foldable f) =>
+   Output.Separator -> f (style,str) -> [(Output.Separator, style, str)]
+attachSeparators sep0 =
+   map (\(sep,(style,x)) -> (sep,style,x)) .
+   zip (sep0 : repeat Output.Empty) . Fold.toList
+
+formatCell ::
+   (Class.Floating a, Output out) =>
+   ArrFormat.Config -> Maybe (Output.Style, a) ->
+   ArrFormat.Tuple a (Output.Style, out)
+formatCell fmt =
+   maybe
+      (ArrFormat.fillTuple
+         (Output.Stored, Output.text $ ArrFormat.configEmpty fmt))
+      (\(style,a) ->
+         (,) style . Output.text <$> ArrFormat.printfFloating fmt a)
+
+
+instance Format Scale where
+   type FormatExtra Scale extra = ()
+   format = formatWithLayout
+
+instance Format Permutation where
+   type FormatExtra Permutation extra = ()
+   format _fmt (Permutation perm) = Perm.format perm
+
+
+{- |
+Layout matrix elements for use in formatting a block matrix.
+Optimally its implementation is reused in 'format' via 'formatWithLayout',
+but sometimes that is not possible.
+-}
+class (Box typ) => Layout typ where
+   type LayoutExtra typ extra :: Constraint
+   {-
+   We use constraint @(Class.Floating a)@ and not @(Format a)@
+   because it allows us to align the components of complex numbers.
+
+   We use a BoxedArray instead of a nested list,
+   although both the underlying formatters
+   and the frontend use a nested list.
+   This gives us a little more type safety for block matrices.
+   -}
+   layout ::
+      (LayoutExtra typ xl, LayoutExtra typ xu,
+       Extent.Measure meas, Extent.C vert, Extent.C horiz,
+       Shape.C height, Shape.C width, Class.Floating a) =>
+      Matrix typ xl xu lower upper meas vert horiz height width a ->
+      BoxedArray.Array (height, width)
+         (Output.Separator, Maybe (Output.Style, a))
+
+instance Layout Scale where
+   type LayoutExtra Scale extra = ()
+   layout (Scale shape a) =
+      let n = Shape.size shape in
+      -- ToDo: 'take' no longer needed when BoxedArray.fromList includes it
+      BoxedArray.fromList (shape,shape) $ take (n*n) $
+      cycle $
+      (Output.Space, Just (Output.Stored, a))
+      :
+      replicate n (Output.Space, Nothing)
+
+instance Layout Permutation where
+   type LayoutExtra Permutation extra = ()
+   layout (Permutation perm) =
+      let sh = Perm.size perm in
+      BoxedArray.fromList (sh,sh) $ concat $
+      map (map ((,) Output.Space . fmap ((,) Output.Stored))) $
+      Perm.layout perm
+
+
+
+instance
+   (MultiplySame typ, MultiplySameExtra typ xl, MultiplySameExtra typ xu,
+    MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper,
+    Extent.Measure meas, Extent.C vert, Extent.C horiz,
+    Shape.C height, Eq height, height ~ width, Class.Floating a) =>
+      Semigroup (Matrix typ xl xu lower upper meas vert horiz height width a)
+         where
+   (<>) = multiplySame
+
+class (Box typ) => MultiplySame typ where
+   type MultiplySameExtra typ extra :: Constraint
+   multiplySame ::
+      (matrix ~ Matrix typ xl xu lower upper meas vert horiz sh sh a,
+       MultiplySameExtra typ xl, MultiplySameExtra typ xu,
+       MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper,
+       Extent.Measure meas, Extent.C vert, Extent.C horiz,
+       Shape.C sh, Eq sh, Class.Floating a) =>
+      matrix -> matrix -> matrix
+
+instance MultiplySame Scale where
+   type MultiplySameExtra Scale extra = extra ~ ()
+   multiplySame =
+      scaleWithCheck "Scale.multiplySame" height
+         (\a (Scale shape b) -> Scale shape $ a*b)
+
+instance MultiplySame Permutation where
+   type MultiplySameExtra Permutation extra = extra ~ ()
+   multiplySame (Permutation a) (Permutation b) =
+      Permutation $ Perm.multiply b a
+
+
+instance
+   (MultiplySame typ, StaticIdentity typ,
+    MultiplySameExtra typ xl, MultiplySameExtra typ xu,
+    StaticIdentityExtra typ xl, StaticIdentityStrip typ lower,
+    StaticIdentityExtra typ xu, StaticIdentityStrip typ upper,
+    MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper,
+    meas ~ Shape, vert ~ Small, horiz ~ Small,
+    Shape.Static height, Eq height, height ~ width, Class.Floating a) =>
+      Monoid (Matrix typ xl xu lower upper meas vert horiz height width a) where
+   mappend = (<>)
+   mempty = staticIdentity
+
+class StaticIdentity typ where
+   type StaticIdentityExtra typ extra :: Constraint
+   type StaticIdentityStrip typ strip :: Constraint
+   staticIdentity ::
+      (StaticIdentityExtra typ xl, StaticIdentityStrip typ lower) =>
+      (StaticIdentityExtra typ xu, StaticIdentityStrip typ upper) =>
+      (Shape.Static sh, Class.Floating a) =>
+      Quadratic typ xl xu lower upper sh a
+
+instance StaticIdentity Scale where
+   type StaticIdentityExtra Scale extra = extra ~ ()
+   type StaticIdentityStrip Scale strip = strip ~ Empty
+   staticIdentity = Scale Shape.static 1
+
+instance StaticIdentity Permutation where
+   type StaticIdentityExtra Permutation extra = extra ~ ()
+   type StaticIdentityStrip Permutation strip = strip ~ Filled
+   staticIdentity = Permutation $ Perm.identity Shape.static
+
+
+scaleWithCheck :: (Eq shape) =>
+   String -> (b -> shape) ->
+   (a -> b -> c) ->
+   Matrix Scale xl xu lower upper meas vert horiz shape shape a ->
+   b -> c
+scaleWithCheck name getSize f (Scale shape a) b =
+   if shape == getSize b
+      then f a b
+      else error $ name ++ ": dimensions mismatch"
+
+
+class Box typ where
+   type BoxExtra typ extra :: Constraint
+   extent ::
+      (BoxExtra typ xl, BoxExtra typ xu) =>
+      (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+      Matrix typ xl xu lower upper meas vert horiz height width a ->
+      Extent.Extent meas vert horiz height width
+   height ::
+      (BoxExtra typ xl, BoxExtra typ xu) =>
+      (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+      Matrix typ xl xu lower upper meas vert horiz height width a -> height
+   height = Extent.height . extent
+   width ::
+      (BoxExtra typ xl, BoxExtra typ xu) =>
+      (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+      Matrix typ xl xu lower upper meas vert horiz height width a -> width
+   width = Extent.width . extent
+
+instance Box Scale where
+   type BoxExtra Scale extra = ()
+   extent (Scale shape _) = Extent.square shape
+   height (Scale shape _) = shape
+   width (Scale shape _) = shape
+
+instance Box Identity where
+   type BoxExtra Identity extra = ()
+   extent (Identity extent_) = extent_
+
+instance Box Permutation where
+   type BoxExtra Permutation extra = ()
+   extent (Permutation perm) = Extent.square $ Perm.size perm
+   height (Permutation perm) = Perm.size perm
+   width (Permutation perm) = Perm.size perm
+
+instance (Eq fuse) => Box (Product fuse) where
+   type BoxExtra (Product fuse) extra =
+         (Box (ProductType extra),
+          BoxExtra (ProductType extra) (ProductExtraL extra),
+          BoxExtra (ProductType extra) (ProductExtraU extra))
+   extent (Product a b) =
+      fromMaybe (error "Matrix.Product: shapes mismatch") $
+      Extent.fuse (extent a) (extent b)
+
+squareSize ::
+   (Box typ, BoxExtra typ xl, BoxExtra typ xu) =>
+   Quadratic typ xl xu lower upper sh a -> sh
+squareSize = height
+
+indices ::
+   (Box typ, BoxExtra typ xl, BoxExtra typ xu,
+    Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+   (Shape.Indexed height, Shape.Indexed width) =>
+   Matrix typ xl xu lower upper meas vert horiz height width a ->
+   [(Shape.Index height, Shape.Index width)]
+indices sh = Shape.indices (height sh, width sh)
+
+
+class (Box typ) => ToQuadratic typ where
+   heightToQuadratic ::
+      (Extent.Measure meas) =>
+      QuadraticMeas typ xl xu lower upper meas height width a ->
+      Quadratic typ xl xu lower upper height a
+   widthToQuadratic ::
+      (Extent.Measure meas) =>
+      QuadraticMeas typ xl xu lower upper meas height width a ->
+      Quadratic typ xl xu lower upper width a
+
+instance ToQuadratic Scale where
+   heightToQuadratic (Scale shape a) = Scale shape a
+   widthToQuadratic (Scale shape a) = Scale shape a
+
+instance ToQuadratic Identity where
+   heightToQuadratic (Identity extent_) =
+      Identity $ Extent.square $ Extent.height extent_
+   widthToQuadratic (Identity extent_) =
+      Identity $ Extent.square $ Extent.width extent_
+
+instance ToQuadratic Permutation where
+   heightToQuadratic (Permutation perm) = Permutation perm
+   widthToQuadratic (Permutation perm) = Permutation perm
+
+
+class (Box typ) => MapExtent typ where
+   type MapExtentExtra typ extra :: Constraint
+   type MapExtentStrip typ strip :: Constraint
+   mapExtent ::
+      (MapExtentExtra typ xl, MapExtentStrip typ lower) =>
+      (MapExtentExtra typ xu, MapExtentStrip typ upper) =>
+      (Extent.Measure measA, Extent.C vertA, Extent.C horizA) =>
+      (Extent.Measure measB, Extent.C vertB, Extent.C horizB) =>
+      ExtentStrict.Map measA vertA horizA measB vertB horizB height width ->
+      Matrix typ xl xu lower upper measA vertA horizA height width a ->
+      Matrix typ xl xu lower upper measB vertB horizB height width a
+
+
+class (Box typ) => Transpose typ where
+   type TransposeExtra typ extra :: Constraint
+   transpose ::
+      (TransposeExtra typ xl, TransposeExtra typ xu) =>
+      (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+      (Shape.C height, Shape.C width, Class.Floating a) =>
+      Matrix typ xl xu lower upper meas vert horiz height width a ->
+      Matrix typ xu xl upper lower meas horiz vert width height a
+
+instance Transpose Scale where
+   type TransposeExtra Scale extra = ()
+   transpose (Scale shape a) = Scale shape a
+
+instance Transpose Identity where
+   type TransposeExtra Identity extra = ()
+   transpose (Identity extent_) = Identity $ Extent.transpose extent_
+
+instance Transpose Permutation where
+   type TransposeExtra Permutation extra = ()
+   transpose (Permutation perm) = Permutation $ Perm.transpose perm
+
+instance (Shape.C fuse, Eq fuse) => Transpose (Product fuse) where
+   type TransposeExtra (Product fuse) extra =
+         (Transpose (ProductType extra),
+          TransposeExtra (ProductType extra) (ProductExtraL extra),
+          TransposeExtra (ProductType extra) (ProductExtraU extra))
+   transpose (Product a b) = Product (transpose b) (transpose a)
+
+
+swapMultiply ::
+   (Transpose typA, Transpose typB) =>
+   (TransposeExtra typA xlA, TransposeExtra typA xuA) =>
+   (TransposeExtra typB xlB, TransposeExtra typB xuB) =>
+   (Extent.Measure measA, Extent.C vertA, Extent.C horizA) =>
+   (Extent.Measure measB, Extent.C vertB, Extent.C horizB) =>
+   (Shape.C heightA, Shape.C widthA) =>
+   (Shape.C heightB, Shape.C widthB) =>
+   (Class.Floating a) =>
+   (matrix ->
+    Matrix typA xuA xlA upperA lowerA measA horizA vertA widthA heightA a ->
+    Matrix typB xuB xlB upperB lowerB measB horizB vertB widthB heightB a) ->
+   Matrix typA xlA xuA lowerA upperA measA vertA horizA heightA widthA a ->
+   matrix ->
+   Matrix typB xlB xuB lowerB upperB measB vertB horizB heightB widthB a
+swapMultiply multiplyTrans a b =
+   transpose $ multiplyTrans b $ transpose a
+
+powerStrips ::
+   (MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper) =>
+   Matrix typ xl xu lower upper meas vert horiz height width a ->
+   (MatrixShape.PowerStripSingleton lower,
+    MatrixShape.PowerStripSingleton upper)
+powerStrips _ =
+   (MatrixShape.powerStripSingleton, MatrixShape.powerStripSingleton)
+
+strips ::
+   (MatrixShape.Strip lower, MatrixShape.Strip upper) =>
+   Matrix typ xl xu lower upper meas vert horiz height width a ->
+   (MatrixShape.StripSingleton lower, MatrixShape.StripSingleton upper)
+strips _ = (MatrixShape.stripSingleton, MatrixShape.stripSingleton)
diff --git a/src/Numeric/LAPACK/Matrix/Wrapper.hs b/src/Numeric/LAPACK/Matrix/Wrapper.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Matrix/Wrapper.hs
@@ -0,0 +1,367 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE EmptyDataDecls #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE ConstraintKinds #-}
+{-# LANGUAGE UndecidableInstances #-}
+module Numeric.LAPACK.Matrix.Wrapper where
+
+import qualified Numeric.LAPACK.Matrix.Type.Private as Matrix
+import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
+import qualified Numeric.LAPACK.Matrix.Extent.Strict as ExtentStrict
+import qualified Numeric.LAPACK.Matrix.Extent as Extent
+import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
+import qualified Numeric.LAPACK.Matrix.Shape.Omni as Omni
+import qualified Numeric.LAPACK.Matrix.Array.Unpacked as Unpacked
+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
+import qualified Numeric.LAPACK.Matrix.Class as MatrixClass
+import qualified Numeric.LAPACK.Matrix.Divide as Divide
+import qualified Numeric.LAPACK.Matrix.Multiply as Multiply
+import Numeric.LAPACK.Matrix.Type.Private (Matrix)
+import Numeric.LAPACK.Matrix.Shape (Filled)
+
+import qualified Type.Data.Num.Unary as Unary
+
+import Data.Tuple.HT (mapPair)
+
+
+
+data MapExtent typ meas
+data instance
+   Matrix (MapExtent typ meas)
+      extraLower extraUpper lower upper meas1 vert1 horiz1 height width a
+         where
+      MapExtent ::
+         (Extent.C vert0, Extent.C horiz0) =>
+         Extent.Map meas0 vert0 horiz0 meas1 vert1 horiz1 height width ->
+         Matrix typ xl xu lower upper meas0 vert0 horiz0 height width a ->
+         Matrix (MapExtent typ meas0) (xl,vert0) (xu,horiz0)
+            lower upper meas1 vert1 horiz1 height width a
+
+type family MapExtentExtra xl
+type instance MapExtentExtra (xl,ex) = xl
+type family MapExtentExtent xl
+type instance MapExtentExtent (xl,ex) = ex
+
+
+instance
+   (Matrix.Box typ, Extent.Measure meas) =>
+      Matrix.Box (MapExtent typ meas) where
+   type BoxExtra (MapExtent typ meas) extra =
+         (Matrix.BoxExtra typ (MapExtentExtra extra))
+   extent (MapExtent m a) = ExtentStrict.apply m $ Matrix.extent a
+
+instance
+   (Matrix.Transpose typ, Extent.Measure meas) =>
+      Matrix.Transpose (MapExtent typ meas) where
+   type TransposeExtra (MapExtent typ meas) extra =
+         (Matrix.TransposeExtra typ (MapExtentExtra extra))
+   transpose (MapExtent m a) =
+      MapExtent (ExtentStrict.transpose m) (Matrix.transpose a)
+
+instance
+   (Matrix.Layout typ, Extent.Measure meas) =>
+      Matrix.Layout (MapExtent typ meas) where
+   type LayoutExtra (MapExtent typ meas) extra =
+         (Matrix.LayoutExtra typ (MapExtentExtra extra))
+   layout (MapExtent _ a) = Matrix.layout a
+
+instance
+   (Matrix.Format typ, Extent.Measure meas) =>
+      Matrix.Format (MapExtent typ meas) where
+   type FormatExtra (MapExtent typ meas) extra =
+         (Matrix.FormatExtra typ (MapExtentExtra extra))
+   format fmt (MapExtent _ a) = Matrix.format fmt a
+
+instance
+   (MatrixClass.Unpack typ, Extent.Measure meas) =>
+      MatrixClass.Unpack (MapExtent typ meas) where
+   type UnpackExtra (MapExtent typ meas) extra =
+         (MatrixClass.UnpackExtra typ (MapExtentExtra extra))
+   unpack (MapExtent m a) =
+      ArrMatrix.liftUnpacked1 id $ ArrMatrix.mapExtent m $ MatrixClass.toFull a
+
+instance
+   (Matrix.MultiplySame typ, Extent.Measure meas) =>
+      Matrix.MultiplySame (MapExtent typ meas) where
+   type MultiplySameExtra (MapExtent typ meas) extra =
+         (Matrix.MultiplySameExtra typ (MapExtentExtra extra))
+   multiplySame (MapExtent m a) (MapExtent _ b) =
+      MapExtent m $ Matrix.multiplySame a b
+
+instance
+   (MatrixClass.Complex typ, Extent.Measure meas) =>
+      MatrixClass.Complex (MapExtent typ meas) where
+   conjugate (MapExtent m a) = MapExtent m $ MatrixClass.conjugate a
+   fromReal (MapExtent m a) = MapExtent m $ MatrixClass.fromReal a
+   toComplex (MapExtent m a) = MapExtent m $ MatrixClass.toComplex a
+
+
+instance
+   (Multiply.MultiplyVector typ, Matrix.ToQuadratic typ, Extent.Measure meas) =>
+      Multiply.MultiplyVector (MapExtent typ meas) where
+   type MultiplyVectorExtra (MapExtent typ meas) extra =
+         (Multiply.MultiplyVectorExtra typ (MapExtentExtra extra),
+          Matrix.BoxExtra typ (MapExtentExtra extra))
+   matrixVector (MapExtent _ a) x = Multiply.matrixVector a x
+   vectorMatrix x (MapExtent _ a) = Multiply.vectorMatrix x a
+
+
+
+data FillStrips typ
+data instance
+   Matrix (FillStrips typ)
+      extraLower extraUpper lower upper meas vert horiz height width a
+         where
+      FillStrips ::
+         (Omni.Strip lower, Omni.Strip upper) =>
+         Matrix typ xl xu lower upper meas vert horiz height width a ->
+         Matrix (FillStrips typ) (xl,lower) (xu,upper)
+            Filled Filled meas vert horiz height width a
+
+type family FillStripsExtra xl
+type instance FillStripsExtra (xl,lower) = xl
+type family FillStripsStrip xl
+type instance FillStripsStrip (xl,lower) = lower
+
+
+instance (Matrix.Box typ) => Matrix.Box (FillStrips typ) where
+   type BoxExtra (FillStrips typ) extra =
+         Matrix.BoxExtra typ (FillStripsExtra extra)
+   extent (FillStrips m) = Matrix.extent m
+
+instance (Matrix.Transpose typ) => Matrix.Transpose (FillStrips typ) where
+   type TransposeExtra (FillStrips typ) extra =
+         Matrix.TransposeExtra typ (FillStripsExtra extra)
+   transpose (FillStrips a) = FillStrips $ Matrix.transpose a
+
+instance (Matrix.Layout typ) => Matrix.Layout (FillStrips typ) where
+   type LayoutExtra (FillStrips typ) extra =
+         Matrix.LayoutExtra typ (FillStripsExtra extra)
+   layout (FillStrips m) = Matrix.layout m
+
+instance (Matrix.Format typ) => Matrix.Format (FillStrips typ) where
+   type FormatExtra (FillStrips typ) extra =
+         Matrix.FormatExtra typ (FillStripsExtra extra)
+   format fmt (FillStrips a) = Matrix.format fmt a
+
+instance (Matrix.ToQuadratic typ) => Matrix.ToQuadratic (FillStrips typ) where
+   heightToQuadratic (FillStrips m) = FillStrips $ Matrix.heightToQuadratic m
+   widthToQuadratic (FillStrips m) = FillStrips $ Matrix.widthToQuadratic m
+
+instance
+   (MatrixClass.Unpack typ) =>
+      MatrixClass.Unpack (FillStrips typ) where
+   type UnpackExtra (FillStrips typ) extra =
+         MatrixClass.UnpackExtra typ (FillStripsExtra extra)
+   unpack (FillStrips m) = Unpacked.fillBoth $ MatrixClass.unpack m
+
+instance (MatrixClass.Complex typ) => MatrixClass.Complex (FillStrips typ) where
+   conjugate (FillStrips m) = FillStrips $ MatrixClass.conjugate m
+   fromReal (FillStrips m) = FillStrips $ MatrixClass.fromReal m
+   toComplex (FillStrips m) = FillStrips $ MatrixClass.toComplex m
+
+
+instance (Matrix.MultiplySame typ) => Matrix.MultiplySame (FillStrips typ) where
+   type MultiplySameExtra (FillStrips typ) extra =
+         (Matrix.MultiplySameExtra typ (FillStripsExtra extra),
+          MatrixShape.PowerStrip (FillStripsStrip extra))
+   multiplySame (FillStrips a) (FillStrips b) =
+      FillStrips $ Matrix.multiplySame a b
+
+instance
+   (Multiply.MultiplyVector typ, Matrix.ToQuadratic typ) =>
+      Multiply.MultiplyVector (FillStrips typ) where
+   type MultiplyVectorExtra (FillStrips typ) extra =
+         (Multiply.MultiplyVectorExtra typ (FillStripsExtra extra),
+          Matrix.BoxExtra typ (FillStripsExtra extra),
+          Omni.Strip (FillStripsStrip extra))
+   matrixVector (FillStrips a) x = Multiply.matrixVector a x
+   vectorMatrix x (FillStrips a) = Multiply.vectorMatrix x a
+
+instance
+   (Multiply.MultiplySquare typ, Matrix.ToQuadratic typ) =>
+      Multiply.MultiplySquare (FillStrips typ) where
+   type MultiplySquareExtra (FillStrips typ) extra =
+         (Multiply.MultiplySquareExtra typ (FillStripsExtra extra),
+          Matrix.BoxExtra typ (FillStripsExtra extra),
+          Omni.Strip (FillStripsStrip extra))
+   transposableSquare trans (FillStrips a) = Multiply.transposableSquare trans a
+   squareFull (FillStrips a) b = Multiply.squareFull a b
+   fullSquare b (FillStrips a) = Multiply.fullSquare b a
+
+instance (Multiply.Power typ) => Multiply.Power (FillStrips typ) where
+   type PowerExtra (FillStrips typ) extra =
+         (Multiply.PowerExtra typ (FillStripsExtra extra),
+          MatrixShape.PowerStrip (FillStripsStrip extra))
+   square (FillStrips a) = FillStrips $ Multiply.square a
+   power n (FillStrips a) = FillStrips $ Multiply.power n a
+   powers1 (FillStrips a) = fmap FillStrips $ Multiply.powers1 a
+
+
+instance (Divide.Determinant typ) => Divide.Determinant (FillStrips typ) where
+   type DeterminantExtra (FillStrips typ) extra =
+         (Divide.DeterminantExtra typ (FillStripsExtra extra))
+   determinant (FillStrips a) = Divide.determinant a
+
+instance
+   (Divide.Solve typ, Matrix.ToQuadratic typ) =>
+      Divide.Solve (FillStrips typ) where
+   type SolveExtra (FillStrips typ) extra =
+         (Divide.SolveExtra typ (FillStripsExtra extra))
+   solve trans (FillStrips a) = Divide.solve trans a
+   solveRight (FillStrips a) b = Divide.solveRight a b
+   solveLeft b (FillStrips a) = Divide.solveLeft b a
+
+instance
+   (Divide.Inverse typ, Matrix.ToQuadratic typ) =>
+      Divide.Inverse (FillStrips typ) where
+   type InverseExtra (FillStrips typ) extra =
+         (Divide.InverseExtra typ (FillStripsExtra extra),
+          MatrixShape.PowerStrip (FillStripsStrip extra))
+   inverse (FillStrips a) = FillStrips $ Divide.inverse a
+
+
+
+{- |
+I do not know, if you will ever need this.
+For diagonal matrices you may not need a wrapper at all
+and for other matrices you may use 'FillStrips'.
+-}
+data PowerStrips typ
+data instance
+   Matrix (PowerStrips typ)
+      extraLower extraUpper lowerf upperf meas vert horiz height width a
+         where
+      PowerStrips ::
+         (Omni.Strip lower, Fill lower ~ lowerf, Omni.PowerStrip lowerf,
+          Omni.Strip upper, Fill upper ~ upperf, Omni.PowerStrip upperf) =>
+         Matrix typ xl xu lower upper meas vert horiz height width a ->
+         Matrix (PowerStrips typ) (xl,lower) (xu,upper)
+            lowerf upperf meas vert horiz height width a
+
+powerStrips ::
+   (Omni.Strip lower, Omni.Strip upper) =>
+   Matrix.QuadraticMeas typ xl xu lower upper meas height width a ->
+   Matrix.QuadraticMeas (PowerStrips typ) (xl,lower) (xu,upper)
+      (Fill lower) (Fill upper) meas height width a
+powerStrips a =
+   case mapPair (filledPowerStripFact, filledPowerStripFact) $
+        Matrix.strips a of
+      (PowerStripFact, PowerStripFact) -> PowerStrips a
+
+type family Fill offDiag
+type instance Fill (Layout.Bands Unary.Zero) = Layout.Bands Unary.Zero
+type instance Fill (Layout.Bands (Unary.Succ k)) = Layout.Filled
+type instance Fill Layout.Filled = Layout.Filled
+
+type family PowerStripsExtra xl
+type instance PowerStripsExtra (xl,lower) = xl
+type family PowerStripsStrip xl
+type instance PowerStripsStrip (xl,lower) = lower
+
+data PowerStripFact c = (Omni.PowerStrip c) => PowerStripFact
+
+filledPowerStripFact ::
+   (Omni.Strip c) => Omni.StripSingleton c -> PowerStripFact (Fill c)
+filledPowerStripFact w =
+   case w of
+      Omni.StripFilled -> PowerStripFact
+      Omni.StripBands Unary.Zero -> PowerStripFact
+      Omni.StripBands Unary.Succ -> PowerStripFact
+
+
+instance (Matrix.Box typ) => Matrix.Box (PowerStrips typ) where
+   type BoxExtra (PowerStrips typ) extra =
+         Matrix.BoxExtra typ (PowerStripsExtra extra)
+   extent (PowerStrips m) = Matrix.extent m
+
+instance (Matrix.Transpose typ) => Matrix.Transpose (PowerStrips typ) where
+   type TransposeExtra (PowerStrips typ) extra =
+         Matrix.TransposeExtra typ (PowerStripsExtra extra)
+   transpose (PowerStrips a) = PowerStrips $ Matrix.transpose a
+
+instance (Matrix.Layout typ) => Matrix.Layout (PowerStrips typ) where
+   type LayoutExtra (PowerStrips typ) extra =
+         Matrix.LayoutExtra typ (PowerStripsExtra extra)
+   layout (PowerStrips m) = Matrix.layout m
+
+instance (Matrix.Format typ) => Matrix.Format (PowerStrips typ) where
+   type FormatExtra (PowerStrips typ) extra =
+         Matrix.FormatExtra typ (PowerStripsExtra extra)
+   format fmt (PowerStrips a) = Matrix.format fmt a
+
+instance
+      (Matrix.MultiplySame typ) => Matrix.MultiplySame (PowerStrips typ) where
+   type MultiplySameExtra (PowerStrips typ) extra =
+         (Matrix.MultiplySameExtra typ (PowerStripsExtra extra),
+          MatrixShape.PowerStrip (PowerStripsStrip extra))
+   multiplySame (PowerStrips a) (PowerStrips b) =
+      PowerStrips $ Matrix.multiplySame a b
+
+instance (Matrix.ToQuadratic typ) => Matrix.ToQuadratic (PowerStrips typ) where
+   heightToQuadratic (PowerStrips m) = PowerStrips $ Matrix.heightToQuadratic m
+   widthToQuadratic (PowerStrips m) = PowerStrips $ Matrix.widthToQuadratic m
+
+instance
+   (MatrixClass.Complex typ) =>
+      MatrixClass.Complex (PowerStrips typ) where
+   conjugate (PowerStrips m) = PowerStrips $ MatrixClass.conjugate m
+   fromReal (PowerStrips m) = PowerStrips $ MatrixClass.fromReal m
+   toComplex (PowerStrips m) = PowerStrips $ MatrixClass.toComplex m
+
+
+instance
+   (Multiply.MultiplyVector typ, Matrix.ToQuadratic typ) =>
+      Multiply.MultiplyVector (PowerStrips typ) where
+   type MultiplyVectorExtra (PowerStrips typ) extra =
+         (Multiply.MultiplyVectorExtra typ (PowerStripsExtra extra),
+          Divide.SolveExtra typ (PowerStripsExtra extra),
+          Matrix.BoxExtra typ (PowerStripsExtra extra),
+          Omni.Strip (PowerStripsStrip extra))
+   matrixVector (PowerStrips a) x = Multiply.matrixVector a x
+   vectorMatrix x (PowerStrips a) = Multiply.vectorMatrix x a
+
+instance
+   (Multiply.MultiplySquare typ, Matrix.ToQuadratic typ) =>
+      Multiply.MultiplySquare (PowerStrips typ) where
+   type MultiplySquareExtra (PowerStrips typ) extra =
+         (Multiply.MultiplySquareExtra typ (PowerStripsExtra extra),
+          Divide.SolveExtra typ (PowerStripsExtra extra),
+          Matrix.BoxExtra typ (PowerStripsExtra extra),
+          Omni.Strip (PowerStripsStrip extra))
+   transposableSquare trans (PowerStrips a) =
+      Multiply.transposableSquare trans a
+   squareFull (PowerStrips a) b = Multiply.squareFull a b
+   fullSquare b (PowerStrips a) = Multiply.fullSquare b a
+
+instance (Multiply.Power typ) => Multiply.Power (PowerStrips typ) where
+   type PowerExtra (PowerStrips typ) extra =
+         (Multiply.PowerExtra typ (PowerStripsExtra extra),
+          MatrixShape.PowerStrip (PowerStripsStrip extra))
+   square (PowerStrips a) = PowerStrips $ Multiply.square a
+   power n (PowerStrips a) = PowerStrips $ Multiply.power n a
+   powers1 (PowerStrips a) = fmap PowerStrips $ Multiply.powers1 a
+
+
+instance (Divide.Determinant typ) => Divide.Determinant (PowerStrips typ) where
+   type DeterminantExtra (PowerStrips typ) extra =
+         (Divide.DeterminantExtra typ (PowerStripsExtra extra))
+   determinant (PowerStrips a) = Divide.determinant a
+
+instance
+   (Divide.Solve typ, Matrix.ToQuadratic typ) =>
+      Divide.Solve (PowerStrips typ) where
+   type SolveExtra (PowerStrips typ) extra =
+         (Divide.SolveExtra typ (PowerStripsExtra extra))
+   solve trans (PowerStrips a) = Divide.solve trans a
+   solveRight (PowerStrips a) b = Divide.solveRight a b
+   solveLeft b (PowerStrips a) = Divide.solveLeft b a
+
+instance
+   (Divide.Inverse typ, Matrix.ToQuadratic typ) =>
+      Divide.Inverse (PowerStrips typ) where
+   type InverseExtra (PowerStrips typ) extra =
+         (Divide.InverseExtra typ (PowerStripsExtra extra),
+          MatrixShape.PowerStrip (PowerStripsStrip extra))
+   inverse (PowerStrips a) = PowerStrips $ Divide.inverse a
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
@@ -28,7 +28,7 @@
 import qualified Numeric.LAPACK.Matrix.Array.Unpacked as Unpacked
 import qualified Numeric.LAPACK.Matrix.Array.Private as ArrMatrix
 import qualified Numeric.LAPACK.Matrix.Basic as Basic
-import qualified Numeric.LAPACK.Matrix.Type as Matrix
+import qualified Numeric.LAPACK.Matrix.Type.Private as Matrix
 import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
 import qualified Numeric.LAPACK.Vector as Vector
 import qualified Numeric.LAPACK.Shape as ExtShape
diff --git a/src/Numeric/LAPACK/Orthogonal/Plain.hs b/src/Numeric/LAPACK/Orthogonal/Plain.hs
--- a/src/Numeric/LAPACK/Orthogonal/Plain.hs
+++ b/src/Numeric/LAPACK/Orthogonal/Plain.hs
@@ -3,16 +3,18 @@
 {-# LANGUAGE EmptyDataDecls #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE ConstraintKinds #-}
 {-# LANGUAGE StandaloneDeriving #-}
 module Numeric.LAPACK.Orthogonal.Plain where
 
 import qualified Numeric.LAPACK.Matrix.Divide as Divide
 import qualified Numeric.LAPACK.Matrix.Multiply as Multiply
-import qualified Numeric.LAPACK.Matrix.Type as Matrix
+import qualified Numeric.LAPACK.Matrix.Type.Private as Matrix
 import qualified Numeric.LAPACK.Matrix.Array.Private as ArrMatrix
 import qualified Numeric.LAPACK.Matrix.Banded.Basic as Banded
 import qualified Numeric.LAPACK.Matrix.Basic as Basic
 import qualified Numeric.LAPACK.Matrix.Private as MatrixPriv
+import qualified Numeric.LAPACK.Matrix.Plain.Format as ArrFormat
 import qualified Numeric.LAPACK.Matrix.Layout as LayoutPub
 import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
 import qualified Numeric.LAPACK.Matrix.Extent.Strict as ExtentStrict
@@ -21,7 +23,7 @@
 import qualified Numeric.LAPACK.Split as Split
 import Numeric.LAPACK.Output ((/+/))
 import Numeric.LAPACK.Matrix.Plain.Format (formatArray)
-import Numeric.LAPACK.Matrix.Type (Matrix, FormatMatrix(formatMatrix))
+import Numeric.LAPACK.Matrix.Type.Private (Matrix)
 import Numeric.LAPACK.Matrix.Triangular.Basic (Upper)
 import Numeric.LAPACK.Matrix.Layout.Private
          (Order(RowMajor, ColumnMajor), sideSwapFromOrder)
@@ -147,12 +149,19 @@
       Layout.caseTallWideSplit shape
 
 
-instance FormatMatrix Hh where
-   formatMatrix fmt (Householder tau m) =
+instance Matrix.Format Hh where
+   type FormatExtra Hh extra = extra ~ ()
+   format fmt (Householder tau m) =
       formatArray fmt (Array.mapShape (Shape.ZeroBased . Shape.size) tau)
       /+/
       formatArray fmt m
 
+instance Matrix.Layout Hh where
+   type LayoutExtra Hh extra = extra ~ ()
+   layout (Householder _tau m) =
+      ArrFormat.splitArrayFromList2 (Layout.splitExtent $ Array.shape m) $
+      ArrFormat.layoutSplit m
+
 fromMatrix ::
    (Extent.Measure meas, Extent.C vert, Extent.C horiz,
     Shape.C height, Shape.C width, Class.Floating a) =>
@@ -516,6 +525,7 @@
 
 
 instance Matrix.Box Hh where
+   type BoxExtra Hh extra = extra ~ ()
    extent = Layout.splitExtent . Array.shape . split_
 
 instance Matrix.ToQuadratic Hh where
@@ -531,10 +541,13 @@
 layoutTauSquare :: sh -> Layout.Diagonal sh
 layoutTauSquare = LayoutPub.diagonal Layout.ColumnMajor
 
-instance (xl ~ (), xu ~ ()) => Matrix.MapExtent Hh xl xu lower upper where
+instance Matrix.MapExtent Hh where
+   type MapExtentExtra Hh extra = extra ~ ()
+   type MapExtentStrip Hh strip = ()
    mapExtent = mapExtent . ExtentStrict.apply
 
-instance (xl ~ (), xu ~ ()) => Multiply.MultiplyVector Hh xl xu where
+instance Multiply.MultiplyVector Hh where
+   type MultiplyVectorExtra Hh extra = extra ~ ()
    matrixVector qr x =
       Array.mapShape deconsUnchecked $
       Basic.unliftColumn Layout.ColumnMajor
@@ -546,7 +559,8 @@
       Basic.unliftColumn Layout.ColumnMajor
          (multiplyQ Transposed NonConjugated qr) x
 
-instance (xl ~ (), xu ~ ()) => Multiply.MultiplySquare Hh xl xu where
+instance Multiply.MultiplySquare Hh where
+   type MultiplySquareExtra Hh extra = extra ~ ()
    squareFull qr =
       ArrMatrix.lift1 $
          multiplyQ NonTransposed NonConjugated qr .
@@ -559,10 +573,12 @@
          multiplyQ Transposed NonConjugated qr .
          Basic.transpose
 
-instance (xl ~ (), xu ~ ()) => Divide.Determinant Hh xl xu where
+instance Divide.Determinant Hh where
+   type DeterminantExtra Hh extra = extra ~ ()
    determinant = determinant
 
-instance (xl ~ (), xu ~ ()) => Divide.Solve Hh xl xu where
+instance Divide.Solve Hh where
+   type SolveExtra Hh extra = extra ~ ()
    solveRight = ArrMatrix.lift1 . leastSquares . mapExtent ExtentPriv.fromSquare
    solveLeft =
       flip $ \a -> ArrMatrix.lift1 $
diff --git a/src/Numeric/LAPACK/Output.hs b/src/Numeric/LAPACK/Output.hs
--- a/src/Numeric/LAPACK/Output.hs
+++ b/src/Numeric/LAPACK/Output.hs
@@ -1,7 +1,7 @@
 module Numeric.LAPACK.Output (
-   Output
-      (text, above, beside, formatRow, formatColumn,
-       formatAligned, formatSeparateTriangle),
+   Output (text, above, beside, formatRow, formatColumn, formatTable),
+   formatAligned, formatSeparateTriangle, decorateTriangle,
+   Separator(..), Style(..),
 
    (/+/),
    (<+>),
@@ -32,9 +32,11 @@
    above :: out -> out -> out
    beside :: out -> out -> out
    formatRow, formatColumn :: [out] -> out
-   formatAligned :: (Foldable f) => [[f out]] -> out
-   formatSeparateTriangle :: (Foldable f) => [[f out]] -> out
+   formatTable :: [[(Separator, Style, out)]] -> out
 
+data Style = Stored | Derived deriving (Eq, Enum)
+
+
 (/+/) :: (Output out) => out -> out -> out
 (/+/) = above
 
@@ -53,15 +55,11 @@
    beside (Html a) (Html b) = Html $ a >> Html.string " " >> b
    formatRow = Html . Html.table . Html.tr . mapM_ (td . unHtml)
    formatColumn = Html . Html.table . mapM_ (Html.tr . td . unHtml)
-   formatAligned =
-      Html . Html.table .
-      mapM_ (Html.tr . mapM_ (td . unHtml) . concatMap Fold.toList)
-   formatSeparateTriangle =
+   formatTable =
+      let applyStyle style = case style of Stored -> id; Derived -> Html.i in
       Html . Html.table .
-      mapM_ (Html.tr . mapM_ td . concat) .
-      zipWith
-         (zipWith $ \it -> map (it . unHtml) . Fold.toList)
-         (iterate (Html.i:) (repeat id))
+      mapM_
+         (Html.tr . mapM_ (\(_sep,style,x) -> td $ applyStyle style $ unHtml x))
 
 td :: Html.Html -> Html.Html
 td = Html.td ! Attr.align (fromString "right")
@@ -73,16 +71,26 @@
    beside = (TextBox.<+>)
    formatRow = TextBox.hsep 1 TextBox.right
    formatColumn = TextBox.vsep 1 TextBox.right
-   formatAligned = alignSeparated . map (concatMap (attachSeparators Space))
-   formatSeparateTriangle =
-      alignSeparated . map concat .
-      zipWith
-         (zipWith attachSeparators)
-         (ListHT.outerProduct
-            (\row col -> if row==col then Bar else Space)
-            [(0::Int)..] [0..])
+   formatTable = alignSeparated . map (map (\(sep,_style,x) -> (sep,x)))
 
 
+formatAligned :: (Foldable f, Output out) => [[f out]] -> out
+formatAligned = formatTable . map (concatMap (plainCells Space Stored))
+
+formatSeparateTriangle :: (Foldable f, Output out) => [[f out]] -> out
+formatSeparateTriangle =
+   formatTable . decorateTriangle (((concat.).) . zipWith3 plainCells)
+
+decorateTriangle :: ([Separator] -> [Style] -> f a -> f b) -> [f a] -> [f b]
+decorateTriangle f =
+   zipWith3 f
+      (iterate (Space:) (Bar : repeat Space))
+      (iterate (Derived:) (repeat Stored))
+
+plainCells :: (Foldable f, Output c) => a -> b -> f c -> [(a, b, c)]
+plainCells sep style = map ((,,) sep style) . Fold.toList
+
+
 data Separator = Empty | Space | Bar
    deriving (Eq, Ord, Show)
 
@@ -102,6 +110,3 @@
 
 formatSeparator :: Separator -> String
 formatSeparator sep = case sep of Empty -> ""; Space -> " "; Bar -> "|"
-
-attachSeparators :: (Foldable f) => Separator -> f str -> [(Separator, str)]
-attachSeparators sep = zip (sep:repeat Empty) . Fold.toList
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
@@ -36,6 +36,7 @@
 import Control.Applicative (liftA2, (<$>))
 
 import qualified Data.Tuple.HT as Tuple
+import Data.Functor.Identity (Identity(Identity))
 import Data.Function.HT (powerAssociative)
 import Data.Monoid (Monoid, mempty, mappend)
 import Data.Semigroup (Semigroup, (<>))
@@ -69,9 +70,17 @@
 format :: (Shape.C sh, Output out) => Permutation sh -> out
 format (Permutation perm) =
    let n = Shape.size $ Array.shape perm
+       s0 = Output.text "."
+       s1 = Output.text "1"
    in formatAligned $
-      map (map ((:[]) . Output.text . (:""))) $
-      map (\k -> (replicate (k-1) '.' ++ '1' : replicate (n-k) '.')) $
+      map (\k -> map Identity $ replicate (k-1) s0 ++ s1 : replicate (n-k) s0) $
+      map (fromIntegral . deconsElement) $ Array.toList perm
+
+layout :: (Shape.C sh, Class.Floating a) => Permutation sh -> [[Maybe a]]
+layout (Permutation perm) =
+   let n = Shape.size $ Array.shape perm
+       z = Nothing
+   in map (\k -> replicate (k-1) z ++ Just one : replicate (n-k) z) $
       map (fromIntegral . deconsElement) $ Array.toList perm
 
 
diff --git a/test/DocTest/Numeric/LAPACK/Permutation/Private.hs b/test/DocTest/Numeric/LAPACK/Permutation/Private.hs
--- a/test/DocTest/Numeric/LAPACK/Permutation/Private.hs
+++ b/test/DocTest/Numeric/LAPACK/Permutation/Private.hs
@@ -1,11 +1,11 @@
 -- Do not edit! Automatically created with doctest-extract from src/Numeric/LAPACK/Permutation/Private.hs
-{-# LINE 46 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
+{-# LINE 47 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
 
 module DocTest.Numeric.LAPACK.Permutation.Private where
 
 import qualified Test.DocTest.Driver as DocTest
 
-{-# LINE 47 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
+{-# LINE 48 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
 import     qualified Test.QuickCheck as QC
 import     Test.Permutation (genPerm, genPivots)
 
@@ -25,18 +25,18 @@
 
 test :: DocTest.T ()
 test = do
- DocTest.printPrefix "Numeric.LAPACK.Permutation.Private:100: "
-{-# LINE 100 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
+ DocTest.printPrefix "Numeric.LAPACK.Permutation.Private:109: "
+{-# LINE 109 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
  DocTest.property
-{-# LINE 100 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
+{-# LINE 109 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
      (QC.forAll QC.arbitraryBoundedEnum $ \inv -> QC.forAll (QC.arbitrary >>= genPivots) $ \xs -> xs == Perm.toPivots inv (Perm.fromPivots inv xs))
- DocTest.printPrefix "Numeric.LAPACK.Permutation.Private:183: "
-{-# LINE 183 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
+ DocTest.printPrefix "Numeric.LAPACK.Permutation.Private:192: "
+{-# LINE 192 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
  DocTest.property
-{-# LINE 183 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
+{-# LINE 192 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
      (QC.forAll genPerm2 $ \(p0,p1) -> determinant (multiply p0 p1) == determinant p0 <> determinant p1)
- DocTest.printPrefix "Numeric.LAPACK.Permutation.Private:212: "
-{-# LINE 212 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
+ DocTest.printPrefix "Numeric.LAPACK.Permutation.Private:221: "
+{-# LINE 221 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
  DocTest.property
-{-# LINE 212 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
+{-# LINE 221 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
      (QC.forAll genPerm2 $ \(p0,p1) -> transpose (multiply p0 p1) == multiply (transpose p1) (transpose p0))
diff --git a/test/Main.hs b/test/Main.hs
--- a/test/Main.hs
+++ b/test/Main.hs
@@ -13,6 +13,7 @@
 import qualified Test.Orthogonal as Orthogonal
 import qualified Test.Singular as Singular
 import qualified Test.Function as Function
+import qualified Test.Block as Block
 import qualified Test.Shape as Shape
 import qualified Test.Permutation as Permutation
 import qualified DocTest.Main as DocTestMain
@@ -60,6 +61,7 @@
    prefix "Orthogonal" Orthogonal.testsVar ++
    prefix "Singular" Singular.testsVar ++
    prefix "Function" Function.testsVar ++
+   prefix "Block" Block.testsVar ++
    []
 
 testsReal ::
diff --git a/test/Test/Banded.hs b/test/Test/Banded.hs
--- a/test/Test/Banded.hs
+++ b/test/Test/Banded.hs
@@ -59,8 +59,8 @@
    Banded (Banded.General sub super height width a)
 
 instance
-   (Show width, Show height, Show a,
-    Shape.C width, Shape.C height, Storable a) =>
+   (Show height, Show width, Show a,
+    Shape.C height, Shape.C width, Storable a) =>
       Show (Banded height width a) where
    showsPrec p (Banded a) = showsPrec p a
 
@@ -103,8 +103,8 @@
       (Banded.General sub super height width a)
 
 instance
-   (Show width, Show height, Show a,
-    Shape.C width, Shape.C height, Storable a) =>
+   (Show height, Show width, Show a,
+    Shape.C height, Shape.C width, Storable a) =>
       Show (Banded2 height width a) where
    showsPrec p (Banded2 a b) =
       showParen True $ showsPrec p a . showString ", " . showsPrec p b
diff --git a/test/Test/BandedHermitian.hs b/test/Test/BandedHermitian.hs
--- a/test/Test/BandedHermitian.hs
+++ b/test/Test/BandedHermitian.hs
@@ -127,8 +127,7 @@
 genSquareShape = Gen.mapGen (const return) Gen.squareShape
 
 genBandedHermitian2 ::
-   (Dim sh, Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix, Eq ix,
-    Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   (ShapeInt ~ sh, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
    Gen.Square sh a (BandedHermitian2 sh a)
 genBandedHermitian2 =
    flip Gen.mapQCDim ((,) <$> genSquareShape <#=#> genSquareShape) $
diff --git a/test/Test/Block.hs b/test/Test/Block.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Block.hs
@@ -0,0 +1,276 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE GADTs #-}
+module Test.Block (testsVar) where
+
+import qualified Test.Divide as Divide
+import qualified Test.Multiply as Multiply
+import qualified Test.Generic as Generic
+import qualified Test.Generator as Gen
+import qualified Test.Logic as Logic
+import qualified Test.Utility as Util
+import Test.Generator ((<|||>), (<===>))
+import Test.Utility (Tagged, prefix)
+
+import qualified Numeric.LAPACK.Matrix.Block as Block
+import qualified Numeric.LAPACK.Matrix.Square as Square
+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
+import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
+import qualified Numeric.LAPACK.Matrix.Layout as Layout
+import qualified Numeric.LAPACK.Matrix.Extent as Extent
+import qualified Numeric.LAPACK.Matrix as Matrix
+import Numeric.LAPACK.Matrix.Shape (Filled)
+import Numeric.LAPACK.Matrix.Extent (Big)
+import Numeric.LAPACK.Matrix (ShapeInt)
+import Numeric.LAPACK.Scalar (RealOf)
+
+import qualified Numeric.Netlib.Class as Class
+
+import Data.Array.Comfort.Shape ((::+))
+
+import Control.Applicative ((<$>))
+
+import qualified Test.QuickCheck as QC
+
+
+type TypeFull = ArrMatrix.Array Layout.Unpacked MatrixShape.Arbitrary
+
+genIdentity ::
+   (Logic.Dim sh, Class.Floating a) => Gen.Square sh a (Matrix.Square sh a)
+genIdentity =
+   Gen.fromBase (return . Matrix.identityFromShape <$> Gen.squareShape)
+
+
+type Diagonal =
+      Block.Diagonal TypeFull () () TypeFull () ()
+         Filled Filled ShapeInt ShapeInt
+
+diagonalTransposeUnpack :: (Class.Floating a) => Diagonal a -> Bool
+diagonalTransposeUnpack a =
+   Util.equalArray
+      (Matrix.toFull (Matrix.transpose a))
+      (Matrix.transpose (Matrix.toFull a))
+
+
+genDiagonal ::
+   (Class.Floating a) => Gen.Square (ShapeInt::+ShapeInt) a (Diagonal a)
+genDiagonal =
+   uncurry Block.Diagonal <$> Gen.stackDiagonal Gen.square Gen.square
+
+genIdentityDiagonal ::
+   (Class.Floating a) => Gen.Square (ShapeInt::+ShapeInt) a (Diagonal a)
+genIdentityDiagonal =
+   uncurry Block.Diagonal <$> Gen.stackDiagonal genIdentity genIdentity
+
+genInvertibleDiagonal ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   Gen.Square (ShapeInt::+ShapeInt) a (Diagonal a)
+genInvertibleDiagonal =
+   uncurry Block.Diagonal <$> Gen.stackDiagonal Gen.invertible Gen.invertible
+
+
+type Square =
+      Block.Square
+         TypeFull () () TypeFull () () TypeFull () () TypeFull () ()
+         Extent.Size Extent.Big Extent.Big ShapeInt ShapeInt
+
+squareTransposeUnpack :: (Class.Floating a) => Square a -> Bool
+squareTransposeUnpack a =
+   Util.equalArray
+      (Matrix.toFull (Matrix.transpose a))
+      (Matrix.transpose (Matrix.toFull a))
+
+
+_genSquareCompose, genSquare ::
+   (Class.Floating a) => Gen.Square (ShapeInt::+ShapeInt) a (Square a)
+_genSquareCompose =
+   (\(a,b) (c,d) -> Block.Square a b c d) <$>
+   ((,) <$> Gen.square <|||> Gen.matrix)
+   <===>
+   ((,) <$> Gen.matrix <|||> Gen.square)
+genSquare =
+   (\m -> case Square.split m of (a,b,c,d) -> Block.Square a b c d)
+   <$>
+   Gen.square
+
+genIdentitySquare ::
+   (Class.Floating a) => Gen.Square (ShapeInt::+ShapeInt) a (Square a)
+genIdentitySquare =
+   (\m -> case Square.split m of (a,b,c,d) -> Block.Square a b c d)
+   <$>
+   genIdentity
+
+genInvertibleSquare ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   Gen.Square (ShapeInt::+ShapeInt) a (Square a)
+genInvertibleSquare =
+   Gen.condition
+      (\m@(Block.Square _a _b _c d) ->
+         Util.invertible d && Util.invertible (Block.schurComplement m))
+      genSquare
+
+
+type Above =
+      Block.Above TypeFull () () TypeFull () () Big ShapeInt ShapeInt ShapeInt
+
+aboveTransposeUnpack :: (Class.Floating a) => Above a -> Bool
+aboveTransposeUnpack a =
+   Util.equalArray
+      (Matrix.toFull (Matrix.transpose a))
+      (Matrix.transpose (Matrix.toFull a))
+
+genAbove ::
+   (Class.Floating a) => Gen.Matrix (ShapeInt::+ShapeInt) ShapeInt a (Above a)
+genAbove = Block.Above <$> Gen.matrix <===> Gen.matrix
+
+
+type Beside =
+      Block.Beside TypeFull () () TypeFull () () Big ShapeInt ShapeInt ShapeInt
+
+besideTransposeUnpack :: (Class.Floating a) => Beside a -> Bool
+besideTransposeUnpack a =
+   Util.equalArray
+      (Matrix.toFull (Matrix.transpose a))
+      (Matrix.transpose (Matrix.toFull a))
+
+genBeside ::
+   (Class.Floating a) => Gen.Matrix ShapeInt (ShapeInt::+ShapeInt) a (Beside a)
+genBeside = Block.Beside <$> Gen.matrix <|||> Gen.matrix
+
+
+type Upper =
+      Block.UpperTriangular TypeFull () () TypeFull () () TypeFull () ()
+         MatrixShape.Filled ShapeInt ShapeInt
+
+upperTransposeUnpack :: (Class.Floating a) => Upper a -> Bool
+upperTransposeUnpack a =
+   Util.equalArray
+      (Matrix.toFull (Matrix.transpose a))
+      (Matrix.transpose (Matrix.toFull a))
+
+
+genUpper ::
+   (Class.Floating a) => Gen.Square (ShapeInt::+ShapeInt) a (Upper a)
+genUpper =
+   (\m -> case Square.split m of (a,b,_c,d) -> Block.Upper a b d)
+   <$>
+   Gen.square
+
+genIdentityUpper ::
+   (Class.Floating a) => Gen.Square (ShapeInt::+ShapeInt) a (Upper a)
+genIdentityUpper =
+   (\m -> case Square.split m of (a,b,_c,d) -> Block.Upper a b d)
+   <$>
+   genIdentity
+
+genInvertibleUpper ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   Gen.Square (ShapeInt::+ShapeInt) a (Upper a)
+genInvertibleUpper =
+   Gen.condition
+      (\(Block.Upper a _ d) -> Util.invertible a && Util.invertible d)
+      genUpper
+
+
+type Lower =
+      Block.LowerTriangular TypeFull () () TypeFull () () TypeFull () ()
+         MatrixShape.Filled ShapeInt ShapeInt
+
+lowerTransposeUnpack :: (Class.Floating a) => Lower a -> Bool
+lowerTransposeUnpack a =
+   Util.equalArray
+      (Matrix.toFull (Matrix.transpose a))
+      (Matrix.transpose (Matrix.toFull a))
+
+
+genLower ::
+   (Class.Floating a) => Gen.Square (ShapeInt::+ShapeInt) a (Lower a)
+genLower =
+   (\m -> case Square.split m of (a,_b,c,d) -> Block.Lower a c d)
+   <$>
+   Gen.square
+
+genIdentityLower ::
+   (Class.Floating a) => Gen.Square (ShapeInt::+ShapeInt) a (Lower a)
+genIdentityLower =
+   (\m -> case Square.split m of (a,_b,c,d) -> Block.Lower a c d)
+   <$>
+   genIdentity
+
+genInvertibleLower ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   Gen.Square (ShapeInt::+ShapeInt) a (Lower a)
+genInvertibleLower =
+   Gen.condition
+      (\(Block.Lower a _ d) -> Util.invertible a && Util.invertible d)
+      genLower
+
+
+type Symmetric =
+      Block.Symmetric
+         TypeFull () () TypeFull () () TypeFull () ()
+         ShapeInt ShapeInt
+
+symmetricTransposeUnpack :: (Class.Floating a) => Symmetric a -> Bool
+symmetricTransposeUnpack a =
+   Util.equalArray
+      (Matrix.toFull (Matrix.transpose a))
+      (Matrix.transpose (Matrix.toFull a))
+
+
+genSymmetric ::
+   (Class.Floating a) => Gen.Square (ShapeInt::+ShapeInt) a (Symmetric a)
+genSymmetric =
+   (\m -> case Square.split m of (a,b,_c,d) -> Block.Symmetric a b d)
+   <$>
+   Gen.square
+
+
+checkForAll ::
+   (Show a, QC.Testable test) =>
+   Gen.T dim tag a -> (a -> test) -> Tagged tag QC.Property
+checkForAll gen = Util.checkForAll (Gen.run gen 3 5)
+
+
+testsVar ::
+   (Show a, Class.Floating a, Eq a, RealOf a ~ ar, Class.Real ar, Show ar) =>
+   [(String, Tagged a QC.Property)]
+testsVar =
+   prefix "Diagonal"
+      (("transposeUnpack", checkForAll genDiagonal diagonalTransposeUnpack) :
+       Generic.testsDistributive genDiagonal ++
+       Multiply.testsVar genIdentityDiagonal genDiagonal ++
+       Divide.testsVar genInvertibleDiagonal) ++
+
+   prefix "Square"
+      (("transposeUnpack", checkForAll genSquare squareTransposeUnpack) :
+       Generic.testsDistributive genSquare ++
+       Multiply.testsVar genIdentitySquare genSquare ++
+       Divide.testsVar genInvertibleSquare) ++
+
+   prefix "Above"
+      (("transposeUnpack", checkForAll genAbove aboveTransposeUnpack) :
+       Generic.testsDistributive genAbove ++
+       Multiply.testsGeneralVar genAbove) ++
+
+   prefix "Beside"
+      (("transposeUnpack", checkForAll genBeside besideTransposeUnpack) :
+       Generic.testsDistributive genBeside ++
+       Multiply.testsGeneralVar genBeside) ++
+
+   prefix "Upper"
+      (("transposeUnpack", checkForAll genUpper upperTransposeUnpack) :
+       Generic.testsDistributive genUpper ++
+       Multiply.testsVar genIdentityUpper genUpper ++
+       Divide.testsVar genInvertibleUpper) ++
+
+   prefix "Lower"
+      (("transposeUnpack", checkForAll genLower lowerTransposeUnpack) :
+       Generic.testsDistributive genLower ++
+       Multiply.testsVar genIdentityLower genLower ++
+       Divide.testsVar genInvertibleLower) ++
+
+   prefix "Symmetric"
+      (("transposeUnpack", checkForAll genSymmetric symmetricTransposeUnpack) :
+       Generic.testsDistributive genSymmetric) ++
+   []
diff --git a/test/Test/Divide.hs b/test/Test/Divide.hs
--- a/test/Test/Divide.hs
+++ b/test/Test/Divide.hs
@@ -1,5 +1,6 @@
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE GADTs #-}
+{-# LANGUAGE ConstraintKinds #-}
 {-# LANGUAGE ExistentialQuantification #-}
 module Test.Divide (
    testsVar,
@@ -11,6 +12,7 @@
    ) where
 
 import qualified Test.Generator as Gen
+import qualified Test.Logic as Logic
 import qualified Test.Utility as Util
 import Test.Generator ((<#\#>), (<#/#>))
 import Test.Utility (Tagged, approxMatrix)
@@ -99,17 +101,13 @@
    (Shape.C size, Eq size, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
    (SquareMatrix size a, Matrix.General size ShapeInt a) -> Bool
 multiplyInverseRight (SquareMatrix a, b) =
-   let aInv = Special.inverse a
-   in case aInv of
-         Special.Inverse _ -> approxRelMatrix b (a #*## (aInv #*## b))
+   approxRelMatrix b (a #*## (Special.Inverse a #*## b))
 
 multiplyInverseLeft ::
    (Shape.C size, Eq size, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
    (Matrix.General ShapeInt size a, SquareMatrix size a) -> Bool
 multiplyInverseLeft (b, SquareMatrix a) =
-   let aInv = Special.inverse a
-   in case aInv of
-         Special.Inverse _ -> approxRelMatrix b ((b ##*# aInv) ##*# a)
+   approxRelMatrix b ((b ##*# Special.Inverse a) ##*# a)
 
 
 solve ::
@@ -153,7 +151,10 @@
 
 data SquareMatrix size a =
    forall typ xl xu lower upper matrix.
-   (Matrix.Solve typ xl xu, Matrix.MultiplySquare typ xl xu,
+   (Matrix.BoxExtra typ xl, Matrix.BoxExtra typ xu,
+    Matrix.Solve typ, Matrix.SolveExtra typ xl, Matrix.SolveExtra typ xu,
+    Matrix.MultiplySquare typ,
+    Matrix.MultiplySquareExtra typ xl, Matrix.MultiplySquareExtra typ xu,
     Matrix.ToQuadratic typ,
     MatrixShape.Strip lower, MatrixShape.Strip upper,
     Matrix.Quadratic typ xl xu lower upper size a ~ matrix, Show matrix) =>
@@ -164,8 +165,9 @@
 
 
 testsVarAny ::
+   (Logic.Dim sh, Shape.C sh, Show sh, Eq sh) =>
    (Show a, Class.Floating a, Eq a, RealOf a ~ ar, Class.Real ar) =>
-   [(String, Gen.MatrixInt a (SquareMatrix ShapeInt a) -> Tagged a QC.Property)]
+   [(String, Gen.Square sh a (SquareMatrix sh a) -> Tagged a QC.Property)]
 testsVarAny =
    ("multiplySolveTrans",
       \gen ->
@@ -186,12 +188,16 @@
    []
 
 testsVar ::
-   (Matrix.Solve typ xl xu, Matrix.MultiplySquare typ xl xu,
+   (Matrix.BoxExtra typ xl, Matrix.BoxExtra typ xu,
+    Matrix.Solve typ, Matrix.SolveExtra typ xl, Matrix.SolveExtra typ xu,
+    Matrix.MultiplySquare typ,
+    Matrix.MultiplySquareExtra typ xl, Matrix.MultiplySquareExtra typ xu,
     Matrix.ToQuadratic typ,
     MatrixShape.Strip lower, MatrixShape.Strip upper,
-    Matrix.Quadratic typ xl xu lower upper ShapeInt a ~ matrix, Show matrix,
+    Logic.Dim sh, Shape.C sh, Show sh, Eq sh,
+    Matrix.Quadratic typ xl xu lower upper sh a ~ matrix, Show matrix,
     Show a, Class.Floating a, Eq a, RealOf a ~ ar, Class.Real ar) =>
-   Gen.MatrixInt a (Matrix.Quadratic typ xl xu lower upper ShapeInt a) ->
+   Gen.Square sh a (Matrix.Quadratic typ xl xu lower upper sh a) ->
    [(String, Tagged a QC.Property)]
 testsVar gen =
    map (mapSnd ($ (SquareMatrix <$> gen))) testsVarAny
diff --git a/test/Test/Generator.hs b/test/Test/Generator.hs
--- a/test/Test/Generator.hs
+++ b/test/Test/Generator.hs
@@ -354,8 +354,11 @@
 squareShape :: (Dim sh) => SquareBase sh (MatrixShape.Square sh)
 squareShape = shapeFromDims MatrixShape.square squareDim
 
-square :: (Class.Floating a) => MatrixInt a (Square.Square ShapeInt a)
+square :: (Dim sh, Class.Floating a) => Square sh a (Square.Square sh a)
 square = mapGen Util.genArray squareShape
+
+squareInt :: (Class.Floating a) => MatrixInt a (Square.Square ShapeInt a)
+squareInt = square
 
 invertible ::
    (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
diff --git a/test/Test/Generic.hs b/test/Test/Generic.hs
--- a/test/Test/Generic.hs
+++ b/test/Test/Generic.hs
@@ -1,17 +1,19 @@
 {-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE ConstraintKinds #-}
 module Test.Generic where
 
 import qualified Test.Generator as Gen
 import qualified Test.Logic as Logic
+import qualified Test.Utility as Util
 import Test.Generator ((<#*|>), (<#=#>))
-import Test.Utility (approxVector, equalArray)
+import Test.Utility (Tagged, approxVector, equalArray)
 
 import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
 import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
 import qualified Numeric.LAPACK.Matrix.Extent as Extent
 import qualified Numeric.LAPACK.Matrix as Matrix
 import Numeric.LAPACK.Matrix.Array (ArrayMatrix)
-import Numeric.LAPACK.Matrix (ShapeInt, (#*|), (#+#), (#-#))
+import Numeric.LAPACK.Matrix (Matrix, ShapeInt, (#*|), (#+#), (#-#))
 import Numeric.LAPACK.Vector (Vector, (|+|), (|-|))
 import Numeric.LAPACK.Scalar (RealOf)
 
@@ -21,7 +23,9 @@
 
 import Control.Applicative ((<$>))
 
+import qualified Test.QuickCheck as QC
 
+
 forceOrder ::
    (MatrixShape.Packing pack, MatrixShape.Property property,
     MatrixShape.Strip lower, MatrixShape.Strip upper,
@@ -63,13 +67,42 @@
 genDistribution gen = genDistribution2 ((,) <$> gen <#=#> gen)
 
 addDistributive, subDistributive ::
-   (MatrixShape.Packing pack, ArrMatrix.Subtractive property,
-    MatrixShape.Strip lower, MatrixShape.Strip upper,
+   (Matrix.MultiplyVector typ,
+    Matrix.MultiplyVectorExtra typ xl, Matrix.MultiplyVectorExtra typ xu) =>
+   (Matrix.Subtractive typ,
+    Matrix.SubtractiveExtra typ xl, Matrix.SubtractiveExtra typ xu) =>
+   (Matrix.AdditiveExtra typ xl, Matrix.AdditiveExtra typ xu) =>
+   (MatrixShape.Strip lower, MatrixShape.Strip upper,
     Extent.Measure meas, Extent.C vert, Extent.C horiz,
-    Shape.C height, Eq height, width ~ ShapeInt,
+    Shape.C height, Eq height, Shape.C width, Eq width,
     Class.Floating a, RealOf a ~ ar, Class.Real ar,
-    ArrayMatrix pack property lower upper meas vert horiz height width a ~
-      matrix) =>
+    Matrix typ xl xu lower upper meas vert horiz height width a ~ matrix) =>
    ((matrix, matrix), Vector width a) -> Bool
 addDistributive ((a,b),x) = approxVector ((a#+#b) #*| x) (a#*|x |+| b#*|x)
 subDistributive ((a,b),x) = approxVector ((a#-#b) #*| x) (a#*|x |-| b#*|x)
+
+
+checkForAll ::
+   (Show a, QC.Testable test) =>
+   Gen.T dim tag a -> (a -> test) -> Tagged tag QC.Property
+checkForAll gen = Util.checkForAll (Gen.run gen 5 10)
+
+testsDistributive ::
+   (Matrix.MultiplyVector typ,
+    Matrix.MultiplyVectorExtra typ xl, Matrix.MultiplyVectorExtra typ xu) =>
+   (Matrix.Subtractive typ,
+    Matrix.SubtractiveExtra typ xl, Matrix.SubtractiveExtra typ xu) =>
+   (Matrix.AdditiveExtra typ xl, Matrix.AdditiveExtra typ xu) =>
+   (MatrixShape.Strip lower, MatrixShape.Strip upper,
+    Extent.Measure meas, Extent.C vert, Extent.C horiz,
+    Logic.Dim height, Show height, Shape.C height, Eq height,
+    Logic.Dim width, Show width, Shape.C width, Eq width,
+    Class.Floating a, RealOf a ~ ar, Class.Real ar, Show a,
+    Matrix typ xl xu lower upper meas vert horiz height width a ~ matrix,
+    Show matrix) =>
+   Gen.Matrix height width a matrix ->
+   [(String, Tagged a QC.Property)]
+testsDistributive gen =
+   ("addDistributive", checkForAll (genDistribution gen) addDistributive) :
+   ("subDistributive", checkForAll (genDistribution gen) subDistributive) :
+   []
diff --git a/test/Test/Hermitian.hs b/test/Test/Hermitian.hs
--- a/test/Test/Hermitian.hs
+++ b/test/Test/Hermitian.hs
@@ -466,14 +466,7 @@
          ((,) <$> genHermitian p <#*|> Gen.vector) Generic.forceOrder) :
    ("forceOrderInverse",
       checkForAll (genHermitian p) Generic.forceOrderInverse) :
-   ("addDistributive",
-      checkForAll
-         (Generic.genDistribution (genHermitian p))
-         Generic.addDistributive) :
-   ("subDistributive",
-      checkForAll
-         (Generic.genDistribution (genHermitian p))
-         Generic.subDistributive) :
+   Generic.testsDistributive (Gen.asMatrixInt $ genHermitian p) ++
 
    ("stack",
       checkForAll
diff --git a/test/Test/Indexed.hs b/test/Test/Indexed.hs
--- a/test/Test/Indexed.hs
+++ b/test/Test/Indexed.hs
@@ -1,4 +1,5 @@
 {-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE ConstraintKinds #-}
 module Test.Indexed where
 
 import qualified Test.Generator as Gen
@@ -32,6 +33,7 @@
 
 genMatrixIndex ::
    (Matrix.Indexed typ,
+    Matrix.BoxExtra typ xl, Matrix.BoxExtra typ xu,
     Extent.Measure meas, Extent.C vert, Extent.C horiz,
     Shape.Indexed height, Shape.Indexed width, Class.Floating a) =>
    GenMatrixNonEmpty a
@@ -42,7 +44,9 @@
 genMatrixIndex = genMatrixIndexGen Matrix.indices
 
 unitDot ::
-   (Matrix.Indexed typ, Matrix.MultiplyVector typ xl xu,
+   (Matrix.BoxExtra typ xl, Matrix.BoxExtra typ xu,
+    Matrix.Indexed typ, Matrix.MultiplyVector typ,
+    Matrix.MultiplyVectorExtra typ xl, Matrix.MultiplyVectorExtra typ xu,
     MatrixShape.Strip lower, MatrixShape.Strip upper,
     Extent.Measure meas, Extent.C vert, Extent.C horiz,
     Shape.Indexed height, Eq height,
diff --git a/test/Test/Matrix.hs b/test/Test/Matrix.hs
--- a/test/Test/Matrix.hs
+++ b/test/Test/Matrix.hs
@@ -585,14 +585,9 @@
          ((,) <$> Gen.matrixInt <#*|> Gen.vector) Generic.forceOrder) :
    ("forceOrderInverse",
       checkForAll Gen.matrixInt Generic.forceOrderInverse) :
-   ("addDistributive",
-      checkForAll
-         (Generic.genDistribution Gen.matrixInt)
-         Generic.addDistributive) :
-   ("subDistributive",
-      checkForAll
-         (Generic.genDistribution Gen.matrixInt)
-         Generic.subDistributive) :
+
+   Generic.testsDistributive Gen.matrixInt ++
+
    ("rowArgAbsMaximums",
       checkForAll Gen.matrix rowArgAbsMaximums) :
 
diff --git a/test/Test/Mosaic.hs b/test/Test/Mosaic.hs
--- a/test/Test/Mosaic.hs
+++ b/test/Test/Mosaic.hs
@@ -10,16 +10,12 @@
 
 import qualified Test.Multiply as Multiply
 import qualified Test.Generator as Gen
-import qualified Test.Logic as Logic
-import qualified Test.Utility as Util
-import Test.Generator ((<-*#>), (<#*|>), (<#*#>))
 import Test.Utility (Tagged)
 
 import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
 import qualified Numeric.LAPACK.Matrix.Shape.Omni as Omni
 import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
 import qualified Numeric.LAPACK.Matrix.Layout as Layout
-import qualified Numeric.LAPACK.Matrix as Matrix
 import Numeric.LAPACK.Matrix (ShapeInt)
 import Numeric.LAPACK.Scalar (RealOf)
 
@@ -50,8 +46,7 @@
       Property diag MatrixShape.Empty MatrixShape.Filled
 
 genMosaic ::
-   (Logic.Dim sh, Shape.Indexed sh, Shape.Index sh ~ ix, Eq ix,
-    Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   (sh ~ ShapeInt, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
    Property prop lo up ->
    Layout.PackingSingleton pack ->
    Gen.Square sh a (ArrMatrix.Quadratic pack prop lo up sh a)
@@ -77,18 +72,6 @@
       Upper     -> repack p <$> Gen.identity `asTypeOf` Gen.triangular
 
 
-checkForAll ::
-   (Show a, QC.Testable test) =>
-   Gen.T dim tag a -> (a -> test) -> Tagged tag QC.Property
-checkForAll gen = Util.checkForAll (Gen.run gen 3 5)
-
-checkForAllExtra ::
-   (Show a, Show b, QC.Testable test) =>
-   QC.Gen a -> Gen.T dim tag b ->
-   (a -> b -> test) -> Tagged tag QC.Property
-checkForAllExtra = Gen.withExtra checkForAll
-
-
 repack ::
    (Shape.C sh, Class.Floating a) =>
    Layout.PackingSingleton pack ->
@@ -97,7 +80,7 @@
 repack pack m =
    case pack of
       Layout.Packed -> m
-      Layout.Unpacked -> Matrix.unpack m
+      Layout.Unpacked -> ArrMatrix.unpack m
 
 testsVar ::
    (MatrixShape.Property prop,
@@ -108,41 +91,4 @@
    Layout.PackingSingleton pack ->
    [(String, Tagged a QC.Property)]
 testsVar prop p =
-   ("multiplySquare",
-      checkForAll (genMosaic prop p) Multiply.multiplySquare) :
-   ("squareSquare",
-      checkForAll (genMosaic prop p) Multiply.squareSquare) :
-   ("power",
-      checkForAllExtra (QC.choose (0,10)) (genMosaic prop p) Multiply.power) :
-
-   ("multiplyIdentityVector",
-      checkForAll
-         ((,) <$> genIdentity prop p <#*|> Gen.vector)
-         Multiply.multiplyIdentityVector) :
-   ("multiplyIdentityFull",
-      checkForAll
-         ((,) <$> genIdentity prop p <#*#> Gen.matrix)
-         Multiply.multiplyIdentityFull) :
-   ("multiplyVector",
-      checkForAll ((,) <$> genMosaic prop p <#*|> Gen.vector)
-         Multiply.multiplyVector) :
-   ("multiplyFull",
-      checkForAll ((,) <$> genMosaic prop p <#*#> Gen.matrix)
-         Multiply.multiplyFull) :
-   ("multiplyVectorLeft",
-      checkForAll
-         ((,) <$> Gen.vector <-*#> genMosaic prop p)
-         Multiply.multiplyVectorLeft) :
-   ("multiplyVectorRight",
-      checkForAll
-         ((,) <$> genMosaic prop p <#*|> Gen.vector)
-         Multiply.multiplyVectorRight) :
-   ("multiplyLeft",
-      checkForAll
-         ((,) <$> Gen.matrix <#*#> genMosaic prop p)
-         Multiply.multiplyLeft) :
-   ("multiplyRight",
-      checkForAll
-         ((,) <$> genMosaic prop p <#*#> Gen.matrix)
-         Multiply.multiplyRight) :
-   []
+   Multiply.testsVar (genIdentity prop p) (genMosaic prop p)
diff --git a/test/Test/Multiply.hs b/test/Test/Multiply.hs
--- a/test/Test/Multiply.hs
+++ b/test/Test/Multiply.hs
@@ -1,4 +1,5 @@
 {-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE ConstraintKinds #-}
 module Test.Multiply (
    multiplySquare,
    squareSquare,
@@ -12,125 +13,264 @@
    multiplyVectorRight,
    multiplyLeft,
    multiplyRight,
+
+   testsGeneralVar,
+   testsVar,
    ) where
 
+import qualified Test.Generator as Gen
+import qualified Test.Logic as Logic
 import qualified Test.Utility as Util
-import Test.Utility (approxArray, approxMatrix, approxVector)
+import Test.Generator ((<-*#>), (<#*|>), (<#*#>))
+import Test.Utility (Tagged, approxArray, approxMatrix, approxVector)
 
 import qualified Numeric.LAPACK.Matrix.Square as Square
 import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
 import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
+import qualified Numeric.LAPACK.Matrix.Extent as Extent
 import qualified Numeric.LAPACK.Matrix as Matrix
 import qualified Numeric.LAPACK.Vector as Vector
-import Numeric.LAPACK.Matrix (ShapeInt, (-*#), (##*#), (#*##), (#*|))
+import Numeric.LAPACK.Matrix (Matrix, ShapeInt, (-*#), (##*#), (#*##), (#*|))
 import Numeric.LAPACK.Vector (Vector)
 import Numeric.LAPACK.Scalar (RealOf)
 
 import qualified Numeric.Netlib.Class as Class
 
+import qualified Data.Array.Comfort.Shape as Shape
 
+import Control.Applicative ((<$>))
 
+import qualified Test.QuickCheck as QC
+
+
+
 multiplySquare ::
-   (Matrix.Power typ xl xu, Matrix.MultiplySquare typ xl xu,
-    Matrix.SquareShape typ, Matrix.ToQuadratic typ,
+   (Matrix.BoxExtra typ xl, Matrix.BoxExtra typ xu,
+    Matrix.Power typ, Matrix.PowerExtra typ xl, Matrix.PowerExtra typ xu,
+    Matrix.MultiplySquare typ,
+    Matrix.MultiplySquareExtra typ xl, Matrix.MultiplySquareExtra typ xu,
+    Matrix.Unpack typ, Matrix.ToQuadratic typ,
+    Matrix.UnpackExtra typ xl, Matrix.UnpackExtra typ xu,
     MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper,
-    Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
-   Matrix.Quadratic typ xl xu lower upper ShapeInt a -> Bool
+    Shape.C sh, Eq sh, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   Matrix.Quadratic typ xl xu lower upper sh a -> Bool
 multiplySquare a =
    Util.approxArray -- (Scalar.selectReal 1e-1 1e-5)
-      (Matrix.toSquare $ Matrix.square a)
-      (a #*## Matrix.toSquare a)
+      (Matrix.toFull $ Matrix.square a)
+      (a #*## Matrix.toFull a)
 
 squareSquare ::
-   (Matrix.Power typ xl xu, Matrix.MultiplySquare typ xl xu,
-    Matrix.SquareShape typ, Matrix.ToQuadratic typ,
+   (Matrix.Power typ, Matrix.PowerExtra typ xl, Matrix.PowerExtra typ xu,
+    Matrix.MultiplySquare typ,
+    Matrix.MultiplySquareExtra typ xl, Matrix.MultiplySquareExtra typ xu,
+    Matrix.Unpack typ, Matrix.ToQuadratic typ,
+    Matrix.UnpackExtra typ xl, Matrix.UnpackExtra typ xu,
     MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper,
-    Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
-   Matrix.Quadratic typ xl xu lower upper ShapeInt a -> Bool
+    Shape.C sh, Eq sh, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   Matrix.Quadratic typ xl xu lower upper sh a -> Bool
 squareSquare a =
    Util.approxArray -- (Scalar.selectReal 1e-1 1e-5)
-      (Matrix.toSquare $ Matrix.square a)
-      (Square.square $ Matrix.toSquare a)
+      (Matrix.toFull $ Matrix.square a)
+      (Square.square $ Matrix.toFull a)
 
 power ::
-   (Matrix.Power typ xl xu, Matrix.MultiplySquare typ xl xu,
-    Matrix.SquareShape typ, Matrix.ToQuadratic typ,
+   (Matrix.BoxExtra typ xl, Matrix.BoxExtra typ xu,
+    Matrix.Power typ, Matrix.PowerExtra typ xl, Matrix.PowerExtra typ xu,
+    Matrix.MultiplySquare typ,
+    Matrix.MultiplySquareExtra typ xl, Matrix.MultiplySquareExtra typ xu,
+    Matrix.Unpack typ, Matrix.ToQuadratic typ,
+    Matrix.UnpackExtra typ xl, Matrix.UnpackExtra typ xu,
     MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper,
-    Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
-   Int -> Matrix.Quadratic typ xl xu lower upper ShapeInt a -> Bool
+    Shape.C sh, Eq sh, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   Int -> Matrix.Quadratic typ xl xu lower upper sh a -> Bool
 power n0 a =
    let n = fromIntegral n0
-       b = Matrix.toSquare (Matrix.power (n+1) a)
-       c = a #*## Matrix.toSquare (Matrix.power n a)
+       b = Matrix.toFull (Matrix.power (n+1) a)
+       c = a #*## Matrix.toFull (Matrix.power n a)
        normInf1 = Vector.normInf1 . ArrMatrix.toVector
    in Util.approxArrayTol (1e-6 * (normInf1 b + normInf1 c)) b c
 
 
 
 multiplyIdentityVector ::
-   (MatrixShape.Packing pack) =>
-   (MatrixShape.Property prop, MatrixShape.PowerStrip lo, MatrixShape.PowerStrip up,
-    Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
-   (ArrMatrix.Quadratic pack prop lo up ShapeInt a, Vector ShapeInt a) -> Bool
+   (Matrix.Unpack typ, Matrix.MultiplyVector typ) =>
+   (Matrix.UnpackExtra typ xl, Matrix.UnpackExtra typ xu) =>
+   (Matrix.MultiplyVectorExtra typ xl, Matrix.MultiplyVectorExtra typ xu) =>
+   (MatrixShape.Strip lo, MatrixShape.Strip up,
+    Shape.C sh, Eq sh, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   (Matrix.Quadratic typ xl xu lo up sh a, Vector sh a) -> Bool
 multiplyIdentityVector (eye,a) = approxVector a (eye #*| a)
 
 multiplyIdentityFull ::
-   (MatrixShape.Packing pack) =>
-   (MatrixShape.Property prop, MatrixShape.PowerStrip lo, MatrixShape.PowerStrip up,
-    Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
-   (ArrMatrix.Quadratic pack prop lo up ShapeInt a,
-    Matrix.General ShapeInt ShapeInt a) ->
+   (Matrix.MultiplySquare typ, Matrix.ToQuadratic typ) =>
+   (Matrix.MultiplySquareExtra typ xl, Matrix.MultiplySquareExtra typ xu) =>
+   (Matrix.BoxExtra typ xl, Matrix.BoxExtra typ xu) =>
+   (MatrixShape.Strip lo, MatrixShape.Strip up,
+    Shape.C sh, Eq sh, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   (Matrix.Quadratic typ xl xu lo up sh a, Matrix.General sh ShapeInt a) ->
    Bool
 multiplyIdentityFull (eye,a) = approxArray a (eye #*## a)
 
 multiplyVector ::
-   (MatrixShape.Packing pack) =>
-   (MatrixShape.Property prop, MatrixShape.PowerStrip lo, MatrixShape.PowerStrip up,
-    Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
-   (ArrMatrix.Quadratic pack prop lo up ShapeInt a, Vector ShapeInt a) -> Bool
-multiplyVector (a,x) = approxVector (Matrix.toSquare a #*| x) (a #*| x)
+   (Matrix.Unpack typ, Matrix.MultiplyVector typ) =>
+   (Matrix.UnpackExtra typ xl, Matrix.UnpackExtra typ xu) =>
+   (Matrix.MultiplyVectorExtra typ xl, Matrix.MultiplyVectorExtra typ xu) =>
+   (MatrixShape.Strip lo, MatrixShape.Strip up,
+    Shape.C sh, Eq sh, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   (Matrix.Quadratic typ xl xu lo up sh a, Vector sh a) -> Bool
+multiplyVector (a,x) = approxVector (Matrix.toFull a #*| x) (a #*| x)
 
 multiplyFull ::
-   (MatrixShape.Packing pack) =>
-   (MatrixShape.Property prop, MatrixShape.PowerStrip lo, MatrixShape.PowerStrip up,
-    Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
-   (ArrMatrix.Quadratic pack prop lo up ShapeInt a,
-    Matrix.General ShapeInt ShapeInt a) ->
+   (Matrix.ToQuadratic typ) =>
+   (Matrix.BoxExtra typ xl, Matrix.BoxExtra typ xu) =>
+   (Matrix.Unpack typ, Matrix.MultiplySquare typ) =>
+   (Matrix.UnpackExtra typ xl, Matrix.UnpackExtra typ xu) =>
+   (Matrix.MultiplySquareExtra typ xl, Matrix.MultiplySquareExtra typ xu) =>
+   (MatrixShape.Strip lo, MatrixShape.Strip up,
+    Shape.C sh, Eq sh, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   (Matrix.Quadratic typ xl xu lo up sh a,
+    Matrix.General sh ShapeInt a) ->
    Bool
-multiplyFull (a,b) = approxArray (Matrix.toSquare a #*## b) (a #*## b)
+multiplyFull (a,b) = approxArray (Matrix.toFull a #*## b) (a #*## b)
 
 
 multiplyVectorLeft ::
-   (MatrixShape.Packing pack) =>
-   (MatrixShape.Property prop, MatrixShape.PowerStrip lo, MatrixShape.PowerStrip up,
+   (Matrix.Unpack typ, Matrix.MultiplyVector typ) =>
+   (Matrix.UnpackExtra typ xl, Matrix.UnpackExtra typ xu) =>
+   (Matrix.MultiplyVectorExtra typ xl, Matrix.MultiplyVectorExtra typ xu) =>
+   (MatrixShape.Strip lo, MatrixShape.Strip up,
+    Extent.Measure meas, Extent.C vert, Extent.C horiz,
+    Shape.C height, Eq height, Shape.C width, Eq width,
     Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
-   (Vector ShapeInt a, ArrMatrix.Quadratic pack prop lo up ShapeInt a) -> Bool
+   (Vector height a, Matrix typ xl xu lo up meas vert horiz height width a) ->
+   Bool
 multiplyVectorLeft (x,a) =
-   approxVector (x -*# Matrix.toSquare a) (x -*# a)
+   approxVector (x -*# Matrix.toFull a) (x -*# a)
 
 multiplyVectorRight ::
-   (MatrixShape.Packing pack) =>
-   (MatrixShape.Property prop, MatrixShape.PowerStrip lo, MatrixShape.PowerStrip up,
+   (Matrix.Unpack typ, Matrix.MultiplyVector typ) =>
+   (Matrix.UnpackExtra typ xl, Matrix.UnpackExtra typ xu) =>
+   (Matrix.MultiplyVectorExtra typ xl, Matrix.MultiplyVectorExtra typ xu) =>
+   (MatrixShape.Strip lo, MatrixShape.Strip up,
+    Extent.Measure meas, Extent.C vert, Extent.C horiz,
+    Shape.C height, Eq height, Shape.C width, Eq width,
     Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
-   (ArrMatrix.Quadratic pack prop lo up ShapeInt a, Vector ShapeInt a) -> Bool
+   (Matrix typ xl xu lo up meas vert horiz height width a, Vector width a) ->
+   Bool
 multiplyVectorRight (a,x) =
-   approxVector (Matrix.toSquare a #*| x) (a #*| x)
+   approxVector (Matrix.toFull a #*| x) (a #*| x)
 
 
 multiplyLeft ::
-   (MatrixShape.Packing pack) =>
-   (MatrixShape.Property prop, MatrixShape.PowerStrip lo, MatrixShape.PowerStrip up,
-    Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
-   (Matrix.General ShapeInt ShapeInt a,
-    ArrMatrix.Quadratic pack prop lo up ShapeInt a) -> Bool
+   (Matrix.ToQuadratic typ) =>
+   (Matrix.BoxExtra typ xl, Matrix.BoxExtra typ xu) =>
+   (Matrix.Unpack typ) =>
+   (Matrix.UnpackExtra typ xl, Matrix.UnpackExtra typ xu) =>
+   (Matrix.MultiplySquare typ) =>
+   (Matrix.MultiplySquareExtra typ xl, Matrix.MultiplySquareExtra typ xu) =>
+   (MatrixShape.Strip lo, MatrixShape.Strip up) =>
+   (Shape.C sh, Eq sh, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   (Matrix.General ShapeInt sh a,
+    Matrix.Quadratic typ xl xu lo up sh a) -> Bool
 multiplyLeft (a,b) =
-   approxMatrix 1e-5 (a ##*# Matrix.toSquare b) (a ##*# b)
+   approxMatrix 1e-5 (a ##*# Matrix.toFull b) (a ##*# b)
 
 multiplyRight ::
-   (MatrixShape.Packing pack) =>
-   (MatrixShape.Property prop, MatrixShape.PowerStrip lo, MatrixShape.PowerStrip up,
-    Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
-   (ArrMatrix.Quadratic pack prop lo up ShapeInt a,
-    Matrix.General ShapeInt ShapeInt a) -> Bool
+   (Matrix.ToQuadratic typ) =>
+   (Matrix.BoxExtra typ xl, Matrix.BoxExtra typ xu) =>
+   (Matrix.Unpack typ) =>
+   (Matrix.UnpackExtra typ xl, Matrix.UnpackExtra typ xu) =>
+   (Matrix.MultiplySquare typ) =>
+   (Matrix.MultiplySquareExtra typ xl, Matrix.MultiplySquareExtra typ xu) =>
+   (MatrixShape.Strip lo, MatrixShape.Strip up) =>
+   (Shape.C sh, Eq sh, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   (Matrix.Quadratic typ xl xu lo up sh a,
+    Matrix.General sh ShapeInt a) -> Bool
 multiplyRight (a,b) =
-   approxArray (Matrix.toSquare a #*## b) (a #*## b)
+   approxArray (Matrix.toFull a #*## b) (a #*## b)
+
+
+
+checkForAll ::
+   (Show a, QC.Testable test) =>
+   Gen.T dim tag a -> (a -> test) -> Tagged tag QC.Property
+checkForAll gen = Util.checkForAll (Gen.run gen 3 5)
+
+checkForAllExtra ::
+   (Show a, Show b, QC.Testable test) =>
+   QC.Gen a -> Gen.T dim tag b ->
+   (a -> b -> test) -> Tagged tag QC.Property
+checkForAllExtra = Gen.withExtra checkForAll
+
+
+testsGeneralVar ::
+   (Matrix.Box typ, Matrix.BoxExtra typ xl, Matrix.BoxExtra typ xu,
+    Matrix.Unpack typ, Matrix.UnpackExtra typ xl, Matrix.UnpackExtra typ xu,
+    Matrix.MultiplyVector typ,
+    Matrix.MultiplyVectorExtra typ xl, Matrix.MultiplyVectorExtra typ xu) =>
+   (MatrixShape.Strip lo, MatrixShape.Strip up) =>
+   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+   (Matrix typ xl xu lo up meas vert horiz height width a ~ matrix) =>
+   (Show matrix) =>
+   (Show height, Eq height, Logic.Dim height,
+    Show width, Eq width, Logic.Dim width,
+    Show a, Class.Floating a, Eq a, RealOf a ~ ar, Class.Real ar) =>
+   Gen.Matrix height width a matrix ->
+   [(String, Tagged a QC.Property)]
+testsGeneralVar genMatrix =
+   ("multiplyVectorLeft",
+      checkForAll ((,) <$> Gen.vector <-*#> genMatrix) multiplyVectorLeft) :
+   ("multiplyVectorRight",
+      checkForAll ((,) <$> genMatrix <#*|> Gen.vector) multiplyVectorRight) :
+   []
+
+testsVar ::
+   (Matrix.ToQuadratic typ) =>
+   (Matrix.Box typ, Matrix.BoxExtra typ xl, Matrix.BoxExtra typ xu,
+    Matrix.Unpack typ, Matrix.UnpackExtra typ xl, Matrix.UnpackExtra typ xu,
+    Matrix.Power typ, Matrix.PowerExtra typ xl, Matrix.PowerExtra typ xu,
+    Matrix.MultiplySquare typ,
+    Matrix.MultiplySquareExtra typ xl, Matrix.MultiplySquareExtra typ xu) =>
+   (Matrix.MultiplyVector typ,
+    Matrix.MultiplyVectorExtra typ xl, Matrix.MultiplyVectorExtra typ xu) =>
+   (MatrixShape.PowerStrip lo, MatrixShape.PowerStrip up) =>
+   (Matrix.Quadratic typ xl xu lo up sh a ~ matrix) =>
+   (Show matrix) =>
+   (Show sh, Eq sh, Logic.Dim sh,
+    Show a, Class.Floating a, Eq a, RealOf a ~ ar, Class.Real ar) =>
+   Gen.Square sh a matrix ->
+   Gen.Square sh a matrix ->
+   [(String, Tagged a QC.Property)]
+testsVar genIdentity genMatrix =
+   testsGeneralVar genMatrix ++
+
+   ("multiplyIdentityVector",
+      checkForAll
+         ((,) <$> genIdentity <#*|> Gen.vector)
+         multiplyIdentityVector) :
+   ("multiplyIdentityFull",
+      checkForAll
+         ((,) <$> genIdentity <#*#> Gen.matrix)
+         multiplyIdentityFull) :
+   ("multiplyVector",
+      checkForAll ((,) <$> genMatrix <#*|> Gen.vector)
+         multiplyVector) :
+   ("multiplyFull",
+      checkForAll ((,) <$> genMatrix <#*#> Gen.matrix)
+         multiplyFull) :
+   ("multiplyLeft",
+      checkForAll
+         ((,) <$> Gen.matrix <#*#> genMatrix)
+         multiplyLeft) :
+   ("multiplyRight",
+      checkForAll
+         ((,) <$> genMatrix <#*#> Gen.matrix)
+         multiplyRight) :
+
+   ("multiplySquare",
+      checkForAll genMatrix multiplySquare) :
+   ("squareSquare",
+      checkForAll genMatrix squareSquare) :
+   ("power",
+      checkForAllExtra (QC.choose (0,10)) genMatrix power) :
+   []
diff --git a/test/Test/Square.hs b/test/Test/Square.hs
--- a/test/Test/Square.hs
+++ b/test/Test/Square.hs
@@ -160,12 +160,12 @@
    [(String, Tagged a QC.Property)]
 testsVar =
    ("multiplySquare",
-      checkForAll Gen.square Multiply.multiplySquare) :
+      checkForAll Gen.squareInt Multiply.multiplySquare) :
    ("squareSquare",
-      checkForAll Gen.square Multiply.squareSquare) :
+      checkForAll Gen.squareInt Multiply.squareSquare) :
    ("power",
       Gen.withExtra checkForAll
-         (QC.choose (0,10)) Gen.square Multiply.power) :
+         (QC.choose (0,10)) Gen.squareInt Multiply.power) :
    ("congruence",
       checkForAll ((,) <$> Gen.square <#*#> Gen.matrix) congruence) :
    ("congruenceAdjoint",
diff --git a/test/Test/Symmetric.hs b/test/Test/Symmetric.hs
--- a/test/Test/Symmetric.hs
+++ b/test/Test/Symmetric.hs
@@ -305,14 +305,7 @@
          ((,) <$> genSymmetric p <#*|> Gen.vector) Generic.forceOrder) :
    ("forceOrderInverse",
       checkForAll (genSymmetric p) Generic.forceOrderInverse) :
-   ("addDistributive",
-      checkForAll
-         (Generic.genDistribution (genSymmetric p))
-         Generic.addDistributive) :
-   ("subDistributive",
-      checkForAll
-         (Generic.genDistribution (genSymmetric p))
-         Generic.subDistributive) :
+   Generic.testsDistributive (Gen.asMatrixInt $ genSymmetric p) ++
 
    ("stack",
       checkForAll
diff --git a/test/Test/Triangular.hs b/test/Test/Triangular.hs
--- a/test/Test/Triangular.hs
+++ b/test/Test/Triangular.hs
@@ -81,7 +81,7 @@
 expandTriangle = Matrix.fromFull . Matrix.toFull
 
 transposedZero ::
-   (Class.Floating a, Shape.C width, Shape.C height) =>
+   (Shape.C height, Shape.C width, Class.Floating a) =>
    General height width a -> General width height a
 transposedZero = ArrMatrix.zero . ArrMatrix.shape . Matrix.transpose
 
@@ -303,12 +303,8 @@
    Util.prefix "Unit" (testsVarExt cont MatrixShape.Unit p) ++
    Util.prefix "Arbitrary" (testsVarExt cont MatrixShape.Arbitrary p) ++
 
-   ("addDistributive",
-      let gen = genTriangular cont MatrixShape.Arbitrary p
-      in checkForAll (Generic.genDistribution gen) Generic.addDistributive) :
-   ("subDistributive",
-      let gen = genTriangular cont MatrixShape.Arbitrary p
-      in checkForAll (Generic.genDistribution gen) Generic.subDistributive) :
+   Generic.testsDistributive
+      (Gen.asMatrixInt $ genTriangular cont MatrixShape.Arbitrary p) ++
 
    ("eigensystem",
       checkForAllExtra Util.genOrder
diff --git a/test/Test/Utility.hs b/test/Test/Utility.hs
--- a/test/Test/Utility.hs
+++ b/test/Test/Utility.hs
@@ -264,7 +264,8 @@
 
 
 invertible ::
-   (Matrix.Determinant typ xl xu,
+   (Matrix.Determinant typ,
+    Matrix.DeterminantExtra typ xl, Matrix.DeterminantExtra typ xu,
     MatrixShape.Strip lower, MatrixShape.Strip upper,
     Shape.C sh, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
    Matrix.Quadratic typ xl xu lower upper sh a -> Bool