diff --git a/Changes.md b/Changes.md
--- a/Changes.md
+++ b/Changes.md
@@ -1,5 +1,19 @@
 # Change log for the `lapack` package
 
+## 0.3.2
+
+ * `Orthogonal`: `project`, `affineKernelFromSpan`, `affineSpanFromKernel`,
+   `leastSquaresConstraint`, `gaussMarkovLinearModel`
+
+ * `Symmetric.fromHermitian`, `Hermitian.fromSymmetric`
+
+ * `instance Monoid Matrix`, especially `mempty`
+   for matrices with static shapes.
+
+ * `Extent.Dimensions`: turn from type family to data family
+
+ * Start using `doctest-extract` for simple tests
+
 ## 0.3.1
 
  * `Matrix.Symmetric`:
diff --git a/Makefile b/Makefile
--- a/Makefile
+++ b/Makefile
@@ -1,10 +1,16 @@
 all:	ReadMe.html Changes.html
 
-run-test:
+run-test:	update-test
 	runhaskell Setup configure --user --enable-tests -fbuildExamples
 	runhaskell Setup build
 	runhaskell Setup haddock
 	./dist/build/lapack-test/lapack-test
 
+run-test-cabal:	update-test
+	cabal test --show-details=streaming
+
+update-test:
+	doctest-extract -i src/ -o test/ --module-prefix DocTest --library-main=Main $$(cat test-module.list)
+
 %.html: %.md
-	pandoc --output=$@ $<
+	pandoc --standalone --metadata pagetitle="$*" --output=$@ $<
diff --git a/ReadMe.md b/ReadMe.md
--- a/ReadMe.md
+++ b/ReadMe.md
@@ -99,7 +99,7 @@
 Vectors are `Storable.Array`s from the
 [comfort-array](https://hackage.haskell.org/package/comfort-array) package.
 An array can have a fancy shape
-like a shape defined by an Haskell enumeration type,
+like a shape defined by an enumeration type,
 the shape of two appended arrays,
 a shape that is compatible to a Haskell container type,
 a rectangular or triangular shape.
@@ -157,7 +157,7 @@
 ~~~~
 
 The relations are defined using two type tags
-in order to support matrix transposition flawlessly.
+in order to support matrix transposition without hassle.
 Using `Small Small` for square matrices
 and `Big Big` for general matrices appears to be arbitrary,
 but is chosen such that altering `Small` to `Big`
@@ -180,6 +180,8 @@
 that always have a unit diagonal by construction.
 The property of a unit diagonal is preserved by some operations
 and enables some optimizations by LAPACK.
+E.g. solving with a unit triangular matrix does not require division
+and thus cannot fail due to division by zero.
 `NonUnit` is a bit of a misnomer.
 A `NonUnit` matrix can still have a unit diagonal,
 but in general it has not and no optimizations will take place.
diff --git a/bug/GGGLM.hs b/bug/GGGLM.hs
new file mode 100644
--- /dev/null
+++ b/bug/GGGLM.hs
@@ -0,0 +1,19 @@
+module Main where
+
+import qualified Numeric.LAPACK.Orthogonal as Ortho
+import qualified Numeric.LAPACK.Matrix as Matrix
+import qualified Numeric.LAPACK.Vector as Vector
+import Numeric.LAPACK.Matrix (shapeInt)
+
+import Data.Array.Comfort.Shape ((:+:)((:+:)))
+
+
+main :: IO ()
+main = do
+   let ab = Matrix.fromList (shapeInt 0) (shapeInt 0 :+: shapeInt 3) ([]::[Double])
+   let d = Vector.fromList (shapeInt 0) []
+   let xy = Vector.fromList (shapeInt 0 :+: shapeInt 3) [-3,1,2]
+   let a = Matrix.tallFromGeneral $ Matrix.takeLeft ab
+   let b = Matrix.takeRight ab
+   print $ Ortho.gaussMarkovLinearModel a b d
+   print $ Ortho.project (Matrix.wideFromGeneral ab) d xy
diff --git a/lapack.cabal b/lapack.cabal
--- a/lapack.cabal
+++ b/lapack.cabal
@@ -1,10 +1,11 @@
+Cabal-Version:    2.2
 Name:             lapack
-Version:          0.3.1
-License:          BSD3
+Version:          0.3.2
+License:          BSD-3-Clause
 License-File:     LICENSE
 Author:           Henning Thielemann <haskell@henning-thielemann.de>
 Maintainer:       Henning Thielemann <haskell@henning-thielemann.de>
-Homepage:         http://hub.darcs.net/thielema/lapack/
+Homepage:         https://hub.darcs.net/thielema/lapack/
 Category:         Data Structures
 Synopsis:         Numerical Linear Algebra using LAPACK
 Description:
@@ -17,7 +18,8 @@
   .
   * Based on @comfort-array@:
     Allows to precisely express one-column or one-row matrices,
-    as well as dense, square, triangular, banded, symmetric and block matrices.
+    as well as dense, square, triangular, banded, symmetric, Hermitian,
+    banded Hermitian, blocked and LU or QR decomposed matrices.
   .
   * Support all data types that are supported by LAPACK,
     i.e. Float, Double, Complex Float, Complex Double
@@ -40,22 +42,26 @@
   .
   See also: @hmatrix@.
 Tested-With:      GHC==7.4.2, GHC==7.8.4, GHC==8.2.2
-Cabal-Version:    >=1.8
 Build-Type:       Simple
 Extra-Source-Files:
   ReadMe.md
   Changes.md
   Makefile
+  test-module.list
 
 Source-Repository this
-  Tag:         0.3.1
+  Tag:         0.3.2
   Type:        darcs
-  Location:    http://hub.darcs.net/thielema/lapack/
+  Location:    https://hub.darcs.net/thielema/lapack/
 
 Source-Repository head
   Type:        darcs
-  Location:    http://hub.darcs.net/thielema/lapack/
+  Location:    https://hub.darcs.net/thielema/lapack/
 
+Flag debug
+  Description: Compile programs for demonstrating bugs in the LAPACK implementation.
+  Default: False
+
 Library
   Build-Depends:
     lapack-ffi >=0.0.1 && <0.1,
@@ -77,6 +83,7 @@
     utility-ht >=0.0.10 && <0.1,
     base >=4.5 && <5
 
+  Default-Language: Haskell98
   GHC-Options:      -Wall
   Hs-Source-Dirs:   src
   Exposed-Modules:
@@ -107,6 +114,7 @@
     Numeric.LAPACK.Format
     Numeric.LAPACK.Output
     Numeric.LAPACK.Example.EconomicAllocation
+    Numeric.LAPACK.Example.DividedDifference
   Other-Modules:
     Numeric.LAPACK.Singular.Plain
     Numeric.LAPACK.Orthogonal.Plain
@@ -136,7 +144,6 @@
     Numeric.LAPACK.Matrix.BandedHermitianPositiveDefinite.Linear
     Numeric.LAPACK.Matrix.Shape.Private
     Numeric.LAPACK.Matrix.Extent.Private
-    Numeric.LAPACK.Matrix.Extent.Kind
     Numeric.LAPACK.Matrix.Divide
     Numeric.LAPACK.Matrix.Multiply
     Numeric.LAPACK.Matrix.Class
@@ -171,6 +178,7 @@
     data-ref >=0.0.1 && <0.1,
     unique-logic-tf >=0.5.1 && <0.6,
     random >=1.1 && <1.2,
+    doctest-exitcode-stdio >=0.0 && <0.1,
     quickcheck-transformer >=0.3 && <0.4,
     QuickCheck >=2.5 && <3,
     ChasingBottoms >=1.2.2 && <1.4,
@@ -181,12 +189,16 @@
     utility-ht,
     base
 
+  Default-Language: Haskell98
   GHC-Options:      -Wall
   Hs-Source-Dirs:   test
   Main-Is:          Main.hs
   Other-Modules:
+    DocTest.Numeric.LAPACK.Example.EconomicAllocation
+    DocTest.Numeric.LAPACK.Example.DividedDifference
+    DocTest.Numeric.LAPACK.Permutation.Private
+    DocTest.Main
     Test.Shape
-    Test.Example
     Test.Indexed
     Test.Divide
     Test.Multiply
@@ -208,3 +220,16 @@
     Test.Logic
     Test.Format
     Test.Utility
+
+Executable lapack-debug
+  Default-Language: Haskell98
+  GHC-Options:      -Wall
+  Hs-Source-Dirs:   bug
+  Main-Is:          GGGLM.hs
+  If flag(debug)
+    Build-Depends:
+      lapack,
+      comfort-array,
+      base
+  Else
+    Buildable: False
diff --git a/src/Numeric/LAPACK/Example/DividedDifference.hs b/src/Numeric/LAPACK/Example/DividedDifference.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Example/DividedDifference.hs
@@ -0,0 +1,116 @@
+{- |
+This module demonstrates triangular matrices.
+
+It verifies that the divided difference scheme
+nicely fits into a triangular matrix,
+where function addition is mapped to matrix addition
+and function multiplication is mapped to matrix multiplication.
+
+<http://en.wikipedia.org/wiki/Divided_difference>
+-}
+module Numeric.LAPACK.Example.DividedDifference where
+
+import qualified Numeric.LAPACK.Matrix.Triangular as Triangular
+import qualified Numeric.LAPACK.Matrix.Square as Square
+import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
+import qualified Numeric.LAPACK.Matrix as Matrix
+import qualified Numeric.LAPACK.Vector as Vector
+import Numeric.LAPACK.Matrix (ShapeInt, (#+#), (#-#))
+import Numeric.LAPACK.Vector (Vector, (|+|), (|-|))
+import Numeric.LAPACK.Format ((##))
+
+import qualified Data.Array.Comfort.Shape as Shape
+import qualified Data.Array.Comfort.Storable as Array
+
+import qualified Data.List as List
+import Data.Semigroup ((<>))
+
+
+{- $setup
+>>> import qualified Test.Utility as Util
+>>> import Test.Utility (approxArray)
+>>>
+>>> import qualified Numeric.LAPACK.Vector as Vector
+>>> import Numeric.LAPACK.Example.DividedDifference (dividedDifferencesMatrix)
+>>> import Numeric.LAPACK.Matrix (ShapeInt, (#+#))
+>>> import Numeric.LAPACK.Vector ((|+|))
+>>>
+>>> import qualified Data.Array.Comfort.Storable as Array
+>>>
+>>> import qualified Test.QuickCheck as QC
+>>>
+>>> import Control.Monad (liftM2)
+>>> import Data.Tuple.HT (mapPair)
+>>> import Data.Semigroup ((<>))
+>>>
+>>> type Vector = Vector.Vector ShapeInt Float
+>>>
+>>> genDD :: QC.Gen (Vector, (Vector, Vector))
+>>> genDD = do
+>>>    (ys0,ys1) <-
+>>>       fmap (mapPair (Vector.autoFromList, Vector.autoFromList) .
+>>>             unzip . take 10) $
+>>>       QC.listOf $ liftM2 (,) (Util.genElement 10) (Util.genElement 10)
+>>>    xs <- Util.genDistinct 10 10 $ Array.shape ys0
+>>>    return (xs,(ys0,ys1))
+-}
+
+
+size :: Vector ShapeInt a -> Int
+size = Shape.zeroBasedSize . Array.shape
+
+subSlices :: Int -> Vector ShapeInt Float -> Vector ShapeInt Float
+subSlices k xs = Vector.drop k xs |-| Vector.take (size xs - k) xs
+
+parameterDifferences :: Vector ShapeInt Float -> [Vector ShapeInt Float]
+parameterDifferences xs = map (flip subSlices xs) [1 .. size xs - 1]
+
+dividedDifferences ::
+   Vector ShapeInt Float -> Vector ShapeInt Float -> [Vector ShapeInt Float]
+dividedDifferences xs ys =
+   scanl
+      (\ddys dxs -> Vector.divide (subSlices 1 ddys) dxs)
+      ys
+      (parameterDifferences xs)
+
+{- |
+prop> QC.forAll genDD $ \(xs, (ys0,ys1)) -> approxArray (dividedDifferencesMatrix xs (ys0|+|ys1)) (dividedDifferencesMatrix xs ys0 #+# dividedDifferencesMatrix xs ys1)
+prop> QC.forAll genDD $ \(xs, (ys0,ys1)) -> approxArray (dividedDifferencesMatrix xs (Vector.mul ys0 ys1)) (dividedDifferencesMatrix xs ys0 <> dividedDifferencesMatrix xs ys1)
+-}
+dividedDifferencesMatrix ::
+   Vector ShapeInt Float -> Vector ShapeInt Float ->
+   Triangular.Upper ShapeInt Float
+dividedDifferencesMatrix xs ys =
+   Triangular.upperFromList MatrixShape.RowMajor (Array.shape xs) $
+   concat $ List.transpose $ map Vector.toList $ dividedDifferences xs ys
+
+
+parameterDifferencesMatrix ::
+   Vector ShapeInt Float -> Triangular.Upper ShapeInt Float
+parameterDifferencesMatrix xs =
+   let ones = Vector.one $ Array.shape xs
+       tp = Matrix.tensorProduct MatrixShape.RowMajor
+   in Triangular.takeUpper $ Square.fromGeneral $ tp ones xs #-# tp xs ones
+
+
+main :: IO ()
+main = do
+   let xs  = Vector.autoFromList [0,1,4,9,16]
+   let ys0 = Vector.autoFromList [3,1,4,1,5]
+   let ys1 = Vector.autoFromList [2,7,1,8,1]
+
+   mapM_ (## "%.4f") $ parameterDifferences xs
+   parameterDifferencesMatrix xs ## "%.4f"
+
+   let ddys0 = dividedDifferencesMatrix xs ys0
+   let ddys1 = dividedDifferencesMatrix xs ys1
+   ddys0 ## "%.4f"
+   ddys1 ## "%.4f"
+   putStrLn ""
+
+   dividedDifferencesMatrix xs (ys0|+|ys1) ## "%.4f"
+   ddys0 #+# ddys1 ## "%.4f"
+   putStrLn ""
+
+   dividedDifferencesMatrix xs (Vector.mul ys0 ys1) ## "%.4f"
+   ddys0 <> ddys1 ## "%.4f"
diff --git a/src/Numeric/LAPACK/Example/EconomicAllocation.hs b/src/Numeric/LAPACK/Example/EconomicAllocation.hs
--- a/src/Numeric/LAPACK/Example/EconomicAllocation.hs
+++ b/src/Numeric/LAPACK/Example/EconomicAllocation.hs
@@ -17,6 +17,17 @@
 import Data.Array.Comfort.Shape ((:+:)((:+:)))
 
 
+{- $setup
+>>> import Numeric.LAPACK.Example.EconomicAllocation
+>>> import Test.Utility (approxVector)
+>>>
+>>> import qualified Numeric.LAPACK.Vector as Vector
+>>> import Numeric.LAPACK.Vector ((+++))
+>>>
+>>> import qualified Data.Array.Comfort.Storable as Array
+-}
+
+
 type ZeroInt2 = ShapeInt:+:ShapeInt
 type Vector sh = Vector.Vector sh Double
 type Matrix height width = Matrix.General height width Double
@@ -75,6 +86,9 @@
    in x |-| p #*| (k #-# Square.identityFrom k) #\| y
 
 
+{- |
+prop> let result = iterated expenses0 balances0 in approxVector result $ compensated expenses0 balances0 +++ Vector.zero (Array.shape $ Vector.takeRight result)
+-}
 main :: IO ()
 main = do
    iterated expenses0 balances0 ## "%10.2f"
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
@@ -55,7 +55,8 @@
 import Numeric.LAPACK.Matrix.Private (Full)
 import Numeric.LAPACK.Linear.Private (solver, withInfo)
 import Numeric.LAPACK.Vector (Vector)
-import Numeric.LAPACK.Private (copyBlock, copyTransposed, copyToColumnMajor)
+import Numeric.LAPACK.Private
+         (copyBlock, copyTransposed, copyToColumnMajor, copyToColumnMajorTemp)
 
 import qualified Numeric.LAPACK.FFI.Generic as LapackGen
 import qualified Numeric.BLAS.FFI.Generic as BlasGen
@@ -166,11 +167,7 @@
             Transposed -> 'T'
       aPtr <-
          case orderLU of
-            RowMajor -> do
-               aPtr <- ContT $ withForeignPtr lu
-               atmpPtr <- Call.allocaArray (n*n)
-               liftIO $ copyToColumnMajor orderLU n n aPtr atmpPtr
-               return atmpPtr
+            RowMajor -> copyToColumnMajorTemp orderLU n n lu
             ColumnMajor -> ContT $ withForeignPtr lu
       ldaPtr <- Call.leadingDim lda
       ipivPtr <- fmap Perm.deconsElementPtr $ ContT $ withForeignPtr ipiv
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
@@ -50,6 +50,7 @@
    Basic.OrderBias, leftBias, rightBias, contiguousBias,
    (|||), beside,
    (===), above,
+   stack,
 
    (|*-),
    tensorProduct,
@@ -605,6 +606,17 @@
 contiguousBias = Basic.ContiguousBias
 
 
+stack ::
+   (Extent.C vert, Extent.C horiz,
+    Shape.C heightA, Eq heightA, Shape.C heightB, Eq heightB,
+    Shape.C widthA, Eq widthA, Shape.C widthB, Eq widthB, Class.Floating a) =>
+   Full vert horiz heightA widthA a -> General heightA widthB a ->
+   General heightB widthA a -> Full vert horiz heightB widthB a ->
+   Full vert horiz (heightA:+:heightB) (widthA:+:widthB) a
+stack = ArrMatrix.lift4 Basic.stack
+
+
+
 rowSums ::
    (Extent.C vert, Extent.C horiz,
     Shape.C height, Shape.C width, Class.Floating a) =>
@@ -745,7 +757,7 @@
 multiplySquare ::
    (Multiply.MultiplySquare typ,
     Type.HeightOf typ ~ height, Eq height, Shape.C width,
-    Extent.C horiz, Extent.C vert, Class.Floating a) =>
+    Extent.C vert, Extent.C horiz, Class.Floating a) =>
    Mod.Transposition -> Type.Matrix typ a ->
    Full vert horiz height width a -> Full vert horiz height width a
 multiplySquare = Multiply.transposableSquare
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,5 +1,6 @@
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE EmptyDataDecls #-}
+{-# LANGUAGE UndecidableInstances #-}
 module Numeric.LAPACK.Matrix.Array (
    Matrix(Array),
    ArrayMatrix,
@@ -178,6 +179,13 @@
 
 instance (Multiply.MultiplySame sh) => Type.MultiplySame (Array sh) where
    multiplySame = lift2 Multiply.same
+
+-- requires UndecidableInstances
+instance
+   (Plain.SquareShape sh, Box.WidthOf sh ~ width, Shape.Static width) =>
+      Type.StaticIdentity (Array sh) where
+   staticIdentity =
+      lift0 $ Plain.identityOrder MatrixShape.RowMajor Shape.static
 
 
 zero ::
diff --git a/src/Numeric/LAPACK/Matrix/Basic.hs b/src/Numeric/LAPACK/Matrix/Basic.hs
--- a/src/Numeric/LAPACK/Matrix/Basic.hs
+++ b/src/Numeric/LAPACK/Matrix/Basic.hs
@@ -11,7 +11,7 @@
          (Order(RowMajor, ColumnMajor), transposeFromOrder, flipOrder)
 import Numeric.LAPACK.Matrix.Modifier (Conjugation(NonConjugated))
 import Numeric.LAPACK.Matrix.Private
-         (Full, Tall, Wide, General, ShapeInt, revealOrder)
+         (Full, Tall, Wide, General, fromFull, ShapeInt, revealOrder)
 import Numeric.LAPACK.Vector (Vector)
 import Numeric.LAPACK.Scalar (RealOf, zero, one)
 import Numeric.LAPACK.Shape.Private (Unchecked(Unchecked))
@@ -347,6 +347,24 @@
    Full Extent.Big horizC (heightA:+:heightB) width a
 above orderBias appendMode a b =
    transpose $ beside orderBias appendMode (transpose a) (transpose b)
+
+stack ::
+   (Extent.C vert, Extent.C horiz,
+    Shape.C heightA, Eq heightA, Shape.C heightB, Eq heightB,
+    Shape.C widthA, Eq widthA, Shape.C widthB, Eq widthB, Class.Floating a) =>
+   Full vert horiz heightA widthA a -> General heightA widthB a ->
+   General heightB widthA a -> Full vert horiz heightB widthB a ->
+   Full vert horiz (heightA:+:heightB) (widthA:+:widthB) a
+stack a b c d =
+   Array.mapShape
+      (\(MatrixShape.Full order _) ->
+         MatrixShape.Full order $
+         Extent.stack
+            (MatrixShape.fullExtent $ Array.shape a)
+            (MatrixShape.fullExtent $ Array.shape d)) $
+   above RightBias Extent.appendAny
+      (beside RightBias Extent.appendAny (fromFull a) b)
+      (beside RightBias Extent.appendAny c (fromFull d))
 
 
 liftRowMajor ::
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
@@ -33,7 +33,7 @@
    {-# MINIMAL solve | solveLeft,solveRight #-}
    solve ::
       (Type.HeightOf typ ~ height, Eq height, Shape.C width,
-       Extent.C horiz, Extent.C vert, Class.Floating a) =>
+       Extent.C vert, Extent.C horiz, Class.Floating a) =>
       Transposition -> Matrix typ a ->
       Full vert horiz height width a -> Full vert horiz height width a
    solve NonTransposed a b = solveRight a b
@@ -41,14 +41,14 @@
 
    solveRight ::
       (Type.HeightOf typ ~ height, Eq height, Shape.C width,
-       Extent.C horiz, Extent.C vert, Class.Floating a) =>
+       Extent.C vert, Extent.C horiz, Class.Floating a) =>
       Matrix typ a ->
       Full vert horiz height width a -> Full vert horiz height width a
    solveRight = solve NonTransposed
 
    solveLeft ::
       (Type.WidthOf typ ~ width, Eq width, Shape.C height,
-       Extent.C horiz, Extent.C vert, Class.Floating a) =>
+       Extent.C vert, Extent.C horiz, Class.Floating a) =>
       Full vert horiz height width a ->
       Matrix typ a -> Full vert horiz height width a
    solveLeft = Basic.swapMultiply $ solve Transposed
@@ -61,14 +61,14 @@
 
 (#\##) ::
    (Solve typ, Type.HeightOf typ ~ height, Eq height, Shape.C width,
-    Extent.C horiz, Extent.C vert, Class.Floating a) =>
+    Extent.C vert, Extent.C horiz, Class.Floating a) =>
    Matrix typ a ->
    Full vert horiz height width a -> Full vert horiz height width a
 (#\##) = solveRight
 
 (##/#) ::
    (Solve typ, Type.WidthOf typ ~ width, Eq width, Shape.C height,
-    Extent.C horiz, Extent.C vert, Class.Floating a) =>
+    Extent.C vert, Extent.C horiz, Class.Floating a) =>
    Full vert horiz height width a ->
    Matrix typ a -> Full vert horiz height width a
 (##/#) = solveLeft
diff --git a/src/Numeric/LAPACK/Matrix/Extent/Kind.hs b/src/Numeric/LAPACK/Matrix/Extent/Kind.hs
deleted file mode 100644
--- a/src/Numeric/LAPACK/Matrix/Extent/Kind.hs
+++ /dev/null
@@ -1,35 +0,0 @@
-{-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE GADTs #-}
-module Numeric.LAPACK.Matrix.Extent.Kind where
-
-
-data General height width =
-   General {
-      generalHeight :: height,
-      generalWidth :: width
-   } deriving (Eq, Show)
-
-data Tall height width =
-   Tall {
-      tallHeight :: height,
-      tallWidth :: width
-   } deriving (Eq, Show)
-
-data Wide height width =
-   Wide {
-      wideHeight :: height,
-      wideWidth :: width
-   } deriving (Eq, Show)
-
-data Square height width =
-   (height ~ width) =>
-   Square {
-      squareSize :: height
-   }
-
-instance (Eq height, Eq width) => Eq (Square height width) where
-   Square a == Square b  =  a==b
-
-instance (Show height, Show width) => Show (Square height width) where
-   showsPrec p (Square s) =
-      showParen (p>10) (showString "Square " . showsPrec 11 s)
diff --git a/src/Numeric/LAPACK/Matrix/Extent/Private.hs b/src/Numeric/LAPACK/Matrix/Extent/Private.hs
--- a/src/Numeric/LAPACK/Matrix/Extent/Private.hs
+++ b/src/Numeric/LAPACK/Matrix/Extent/Private.hs
@@ -4,7 +4,6 @@
 {-# LANGUAGE GADTs #-}
 module Numeric.LAPACK.Matrix.Extent.Private where
 
-import qualified Numeric.LAPACK.Matrix.Extent.Kind as EK
 import Numeric.LAPACK.Shape.Private (Unchecked(deconsUnchecked))
 import Numeric.LAPACK.Wrapper (Flip(Flip, getFlip))
 
@@ -13,15 +12,10 @@
 import Control.DeepSeq (NFData, rnf)
 
 import Data.Maybe.HT (toMaybe)
-import Data.Tuple.HT (swap)
 import Data.Eq.HT (equating)
 
 
-data Extent vertical horizontal height width =
-   Extent {
-      extentDir :: (vertical,horizontal),
-      extentDim :: Dimensions vertical horizontal height width
-   }
+data family Extent vertical horizontal :: * -> * -> *
 
 instance
    (C vertical, C horizontal, NFData height, NFData width) =>
@@ -29,10 +23,10 @@
    rnf =
       getAccessor $
       switchTagPair
-         (Accessor $ \(Extent o (EK.Square s)) -> rnf (o,s))
-         (Accessor $ \(Extent o (EK.Wide h w)) -> rnf (o,(h,w)))
-         (Accessor $ \(Extent o (EK.Tall h w)) -> rnf (o,(h,w)))
-         (Accessor $ \(Extent o (EK.General h w)) -> rnf (o,(h,w)))
+         (Accessor $ \(Square s) -> rnf s)
+         (Accessor $ \(Wide h w) -> rnf (h,w))
+         (Accessor $ \(Tall h w) -> rnf (h,w))
+         (Accessor $ \(General h w) -> rnf (h,w))
 
 
 data Big = Big deriving (Eq,Show)
@@ -47,25 +41,42 @@
 type Square sh = Extent Small Small sh sh
 
 
-type family Dimensions vertical horizontal :: * -> * -> *
+data instance Extent Big Big height width =
+   General {
+      generalHeight :: height,
+      generalWidth :: width
+   }
 
-type instance Dimensions Big Big = EK.General
-type instance Dimensions Big Small = EK.Tall
-type instance Dimensions Small Big = EK.Wide
-type instance Dimensions Small Small = EK.Square
+data instance Extent Big Small height width =
+   Tall {
+      tallHeight :: height,
+      tallWidth :: width
+   }
 
+data instance Extent Small Big height width =
+   Wide {
+      wideHeight :: height,
+      wideWidth :: width
+   }
 
+data instance Extent Small Small height width =
+   (height ~ width) =>
+   Square {
+      squareSize :: height
+   }
+
+
 general :: height -> width -> General height width
-general h w = Extent (Big,Big) $ EK.General h w
+general = General
 
 tall :: height -> width -> Tall height width
-tall h w = Extent (Big,Small) $ EK.Tall h w
+tall = Tall
 
 wide :: height -> width -> Wide height width
-wide h w = Extent (Small,Big) $ EK.Wide h w
+wide = Wide
 
 square :: sh -> Square sh
-square sh = Extent (Small,Small) $ EK.Square sh
+square = Square
 
 
 newtype Map vertA horizA vertB horizB height width =
@@ -106,10 +117,10 @@
 caseTallWide ge =
    getCaseTallWide $
    switchTagPair
-      (CaseTallWide $ \(Extent _ (EK.Square sh)) -> Left $ tall sh sh)
+      (CaseTallWide $ \(Square sh) -> Left $ tall sh sh)
       (CaseTallWide Right)
       (CaseTallWide Left)
-      (CaseTallWide $ \(Extent _ (EK.General h w)) ->
+      (CaseTallWide $ \(General h w) ->
          if ge h w
             then Left $ tall h w
             else Right $ wide h w)
@@ -138,8 +149,8 @@
 generalizeTall =
    getGenTall $
    switchTagPair
-      (GenTall id) (GenTall $ \(Extent _ (EK.Square s)) -> wide s s)
-      (GenTall id) (GenTall $ \(Extent _ (EK.Tall h w)) -> general h w)
+      (GenTall id) (GenTall $ \(Square s) -> wide s s)
+      (GenTall id) (GenTall $ \(Tall h w) -> general h w)
 
 newtype GenWide height width vert horiz =
    GenWide {
@@ -154,8 +165,8 @@
    switchTagPair
       (GenWide id)
       (GenWide id)
-      (GenWide $ \(Extent _ (EK.Square s)) -> tall s s)
-      (GenWide $ \(Extent _ (EK.Wide h w)) -> general h w)
+      (GenWide $ \(Square s) -> tall s s)
+      (GenWide $ \(Wide h w) -> general h w)
 
 
 newtype GenToTall height width vert horiz =
@@ -169,8 +180,8 @@
 genToTall =
    getGenToTall $
    switchTagPair
-      (GenToTall $ \(Extent _ (EK.Square s)) -> tall s s)
-      (GenToTall $ \(Extent _ (EK.Wide h w)) -> general h w)
+      (GenToTall $ \(Square s) -> tall s s)
+      (GenToTall $ \(Wide h w) -> general h w)
       (GenToTall id)
       (GenToTall id)
 
@@ -186,16 +197,12 @@
 genToWide =
    getGenToWide $
    switchTagPair
-      (GenToWide $ \(Extent _ (EK.Square s)) -> wide s s)
+      (GenToWide $ \(Square s) -> wide s s)
       (GenToWide id)
-      (GenToWide $ \(Extent _ (EK.Tall h w)) -> general h w)
+      (GenToWide $ \(Tall h w) -> general h w)
       (GenToWide id)
 
 
-squareSize :: Square sh -> sh
-squareSize (Extent (Small,Small) (EK.Square sh)) = sh
-
-
 newtype Accessor a height width vert horiz =
    Accessor {getAccessor :: Extent vert horiz height width -> a}
 
@@ -203,19 +210,19 @@
 height =
    getAccessor $
    switchTagPair
-      (Accessor (\(Extent _ (EK.Square s)) -> s))
-      (Accessor (EK.wideHeight . extentDim))
-      (Accessor (EK.tallHeight . extentDim))
-      (Accessor (EK.generalHeight . extentDim))
+      (Accessor squareSize)
+      (Accessor wideHeight)
+      (Accessor tallHeight)
+      (Accessor generalHeight)
 
 width :: (C vert, C horiz) => Extent vert horiz height width -> width
 width =
    getAccessor $
    switchTagPair
-      (Accessor (\(Extent _ (EK.Square s)) -> s))
-      (Accessor (EK.wideWidth . extentDim))
-      (Accessor (EK.tallWidth . extentDim))
-      (Accessor (EK.generalWidth . extentDim))
+      (Accessor (\(Square s) -> s))
+      (Accessor wideWidth)
+      (Accessor tallWidth)
+      (Accessor generalWidth)
 
 
 dimensions ::
@@ -232,7 +239,7 @@
 
 fromSquareLiberal :: (C vert, C horiz) =>
    Extent Small Small height width -> Extent vert horiz height width
-fromSquareLiberal x@(Extent _ (EK.Square _)) = genSquare $ height x
+fromSquareLiberal (Square s) = genSquare s
 
 squareFromGeneral ::
    (C vert, C horiz, Eq size) =>
@@ -258,10 +265,10 @@
 transpose =
    getTranspose $
    switchTagPair
-      (Transpose $ \(Extent o (EK.Square s)) -> Extent o (EK.Square s))
-      (Transpose $ \(Extent o (EK.Wide h w)) -> Extent (swap o) (EK.Tall w h))
-      (Transpose $ \(Extent o (EK.Tall h w)) -> Extent (swap o) (EK.Wide w h))
-      (Transpose $ \(Extent o (EK.General h w)) -> Extent o (EK.General w h))
+      (Transpose $ \(Square s) -> Square s)
+      (Transpose $ \(Wide h w) -> Tall w h)
+      (Transpose $ \(Tall h w) -> Wide w h)
+      (Transpose $ \(General h w) -> General w h)
 
 
 newtype Equal height width vert horiz =
@@ -277,10 +284,11 @@
    (==) =
       getEqual $
       switchTagPair
-         (Equal $ equating extentDim)
-         (Equal $ equating extentDim)
-         (Equal $ equating extentDim)
-         (Equal $ equating extentDim)
+         (Equal $ \(Square a) (Square b) -> a==b)
+         (Equal $ \a b -> equating wideHeight a b && equating wideWidth a b)
+         (Equal $ \a b -> equating tallHeight a b && equating tallWidth a b)
+         (Equal $ \a b ->
+            equating generalHeight a b && equating generalWidth a b)
 
 
 instance
@@ -324,8 +332,8 @@
 widen w =
    getWiden $
    switchTag
-      (Widen (\(Extent o x) -> Extent o (x{EK.wideWidth = w})))
-      (Widen (\(Extent o x) -> Extent o (x{EK.generalWidth = w})))
+      (Widen (\x -> x{wideWidth = w}))
+      (Widen (\x -> x{generalWidth = w}))
 
 reduceWideHeight ::
    (C vert) =>
@@ -333,8 +341,8 @@
 reduceWideHeight h =
    getWiden $
    switchTag
-      (Widen (\(Extent o x) -> Extent o (x{EK.wideHeight = h})))
-      (Widen (\(Extent o x) -> Extent o (x{EK.generalHeight = h})))
+      (Widen (\x -> x{wideHeight = h}))
+      (Widen (\x -> x{generalHeight = h}))
 
 
 newtype Adapt heightA widthA heightB widthB vert horiz =
@@ -351,10 +359,10 @@
 reduceConsistent h w =
    getAdapt $
    switchTagPair
-      (Adapt $ \(Extent o (EK.Square _)) -> Extent o (EK.Square h))
-      (Adapt $ \(Extent o (EK.Wide _ _)) -> Extent o (EK.Wide h w))
-      (Adapt $ \(Extent o (EK.Tall _ _)) -> Extent o (EK.Tall h w))
-      (Adapt $ \(Extent o (EK.General _ _)) -> Extent o (EK.General h w))
+      (Adapt $ \(Square _) -> Square h)
+      (Adapt $ \(Wide _ _) -> Wide h w)
+      (Adapt $ \(Tall _ _) -> Tall h w)
+      (Adapt $ \(General _ _) -> General h w)
 
 
 class (C vert, C horiz) => GeneralTallWide vert horiz where
@@ -371,9 +379,9 @@
 mapHeight f =
    getAdapt $
    switchTagGTW
-      (Adapt $ \(Extent o (EK.Wide h w)) -> Extent o (EK.Wide (f h) w))
-      (Adapt $ \(Extent o (EK.Tall h w)) -> Extent o (EK.Tall (f h) w))
-      (Adapt $ \(Extent o (EK.General h w)) -> Extent o (EK.General (f h) w))
+      (Adapt $ \(Wide h w) -> Wide (f h) w)
+      (Adapt $ \(Tall h w) -> Tall (f h) w)
+      (Adapt $ \(General h w) -> General (f h) w)
 
 mapWidth ::
    (GeneralTallWide vert horiz) =>
@@ -382,12 +390,12 @@
 mapWidth f =
    getAdapt $
    switchTagGTW
-      (Adapt $ \(Extent o (EK.Wide h w)) -> Extent o (EK.Wide h (f w)))
-      (Adapt $ \(Extent o (EK.Tall h w)) -> Extent o (EK.Tall h (f w)))
-      (Adapt $ \(Extent o (EK.General h w)) -> Extent o (EK.General h (f w)))
+      (Adapt $ \(Wide h w) -> Wide h (f w))
+      (Adapt $ \(Tall h w) -> Tall h (f w))
+      (Adapt $ \(General h w) -> General h (f w))
 
 mapSquareSize :: (shA -> shB) -> Square shA -> Square shB
-mapSquareSize f (Extent o (EK.Square s)) = Extent o (EK.Square (f s))
+mapSquareSize f (Square s) = Square (f s)
 
 
 mapWrap ::
@@ -399,11 +407,10 @@
 mapWrap fh fw =
    getAdapt $
    switchTagPair
-      (Adapt $ \(Extent o (EK.Square h)) -> Extent o (EK.Square (fh h)))
-      (Adapt $ \(Extent o (EK.Wide h w)) -> Extent o (EK.Wide (fh h) (fw w)))
-      (Adapt $ \(Extent o (EK.Tall h w)) -> Extent o (EK.Tall (fh h) (fw w)))
-      (Adapt $ \(Extent o (EK.General h w)) ->
-                  Extent o (EK.General (fh h) (fw w)))
+      (Adapt $ \(Square h) -> Square (fh h))
+      (Adapt $ \(Wide h w) -> Wide (fh h) (fw w))
+      (Adapt $ \(Tall h w) -> Tall (fh h) (fw w))
+      (Adapt $ \(General h w) -> General (fh h) (fw w))
 
 {- only admissible since GHC-7.8
 mapUnwrap ::
@@ -415,11 +422,10 @@
 mapUnwrap fh fw =
    getAdapt $
    switchTagPair
-      (Adapt $ \(Extent o (EK.Square h)) -> Extent o (EK.Square (fh h)))
-      (Adapt $ \(Extent o (EK.Wide h w)) -> Extent o (EK.Wide (fh h) (fw w)))
-      (Adapt $ \(Extent o (EK.Tall h w)) -> Extent o (EK.Tall (fh h) (fw w)))
-      (Adapt $ \(Extent o (EK.General h w)) ->
-                  Extent o (EK.General (fh h) (fw w)))
+      (Adapt $ \(Square h) -> Square (fh h))
+      (Adapt $ \(Wide h w) -> Wide (fh h) (fw w))
+      (Adapt $ \(Tall h w) -> Tall (fh h) (fw w))
+      (Adapt $ \(General h w) -> General (fh h) (fw w))
 -}
 
 recheck ::
@@ -429,14 +435,11 @@
 recheck =
    getAdapt $
    switchTagPair
-      (Adapt $ \(Extent o (EK.Square h)) ->
-         Extent o (EK.Square (deconsUnchecked h)))
-      (Adapt $ \(Extent o (EK.Wide h w)) ->
-         Extent o (EK.Wide (deconsUnchecked h) (deconsUnchecked w)))
-      (Adapt $ \(Extent o (EK.Tall h w)) ->
-         Extent o (EK.Tall (deconsUnchecked h) (deconsUnchecked w)))
-      (Adapt $ \(Extent o (EK.General h w)) ->
-         Extent o (EK.General (deconsUnchecked h) (deconsUnchecked w)))
+      (Adapt $ \(Square h) -> Square (deconsUnchecked h))
+      (Adapt $ \(Wide h w) -> Wide (deconsUnchecked h) (deconsUnchecked w))
+      (Adapt $ \(Tall h w) -> Tall (deconsUnchecked h) (deconsUnchecked w))
+      (Adapt $ \(General h w) ->
+                              General (deconsUnchecked h) (deconsUnchecked w))
 
 
 
@@ -456,48 +459,18 @@
 fuse =
    getFuse $
    switchTagPair
-      (Fuse $
-       \(Extent o (EK.Square s0)) (Extent _ (EK.Square s1)) ->
-         toMaybe (s0==s1) $ Extent o (EK.Square s0))
-      (Fuse $
-       \(Extent o (EK.Wide h f0)) (Extent _ (EK.Wide f1 w)) ->
-         toMaybe (f0==f1) $ Extent o (EK.Wide h w))
-      (Fuse $
-       \(Extent o (EK.Tall h f0)) (Extent _ (EK.Tall f1 w)) ->
-         toMaybe (f0==f1) $ Extent o (EK.Tall h w))
-      (Fuse $
-       \(Extent o (EK.General h f0)) (Extent _ (EK.General f1 w)) ->
-         toMaybe (f0==f1) $ Extent o (EK.General h w))
+      (Fuse $ \(Square s0) (Square s1) -> toMaybe (s0==s1) $ Square s0)
+      (Fuse $ \(Wide h f0) (Wide f1 w) -> toMaybe (f0==f1) $ Wide h w)
+      (Fuse $ \(Tall h f0) (Tall f1 w) -> toMaybe (f0==f1) $ Tall h w)
+      (Fuse $ \(General h f0) (General f1 w) -> toMaybe (f0==f1) $ General h w)
 
 
-newtype Kronecker heightA widthA heightB widthB vert horiz =
-   Kronecker {
-      getKronecker ::
-         Extent vert horiz heightA widthA ->
-         Extent vert horiz heightB widthB ->
-         Extent vert horiz (heightA,heightB) (widthA,widthB)
-   }
-
 kronecker ::
    (C vert, C horiz) =>
    Extent vert horiz heightA widthA ->
    Extent vert horiz heightB widthB ->
    Extent vert horiz (heightA,heightB) (widthA,widthB)
-kronecker =
-   getKronecker $
-   switchTagPair
-      (Kronecker $
-       \(Extent o (EK.Square s0)) (Extent _ (EK.Square s1)) ->
-         Extent o (EK.Square (s0,s1)))
-      (Kronecker $
-       \(Extent o (EK.Wide h0 w0)) (Extent _ (EK.Wide h1 w1)) ->
-         Extent o (EK.Wide (h0,h1) (w0,w1)))
-      (Kronecker $
-       \(Extent o (EK.Tall h0 w0)) (Extent _ (EK.Tall h1 w1)) ->
-         Extent o (EK.Tall (h0,h1) (w0,w1)))
-      (Kronecker $
-       \(Extent o (EK.General h0 w0)) (Extent _ (EK.General h1 w1)) ->
-         Extent o (EK.General (h0,h1) (w0,w1)))
+kronecker = stackGen (,) (,)
 
 
 
@@ -561,6 +534,41 @@
    AppendMode vertA vertB (Append vertA vertB) height widthA widthB
 appendAny =
    getAppendAny $ switchTag (AppendAny appendLeftAux) (AppendAny appendRight)
+
+
+stack ::
+   (C vert, C horiz) =>
+   Extent vert horiz heightA widthA ->
+   Extent vert horiz heightB widthB ->
+   Extent vert horiz (heightA:+:heightB) (widthA:+:widthB)
+stack = stackGen (:+:) (:+:)
+
+newtype Stack f heightA widthA heightB widthB vert horiz =
+   Stack {
+      getStack ::
+         Extent vert horiz heightA widthA ->
+         Extent vert horiz heightB widthB ->
+         Extent vert horiz (f heightA heightB) (f widthA widthB)
+   }
+
+stackGen ::
+   (C vert, C horiz) =>
+   (heightA -> heightB -> f heightA heightB) ->
+   (widthA -> widthB -> f widthA widthB) ->
+   Extent vert horiz heightA widthA ->
+   Extent vert horiz heightB widthB ->
+   Extent vert horiz (f heightA heightB) (f widthA widthB)
+stackGen fh fw =
+   getStack $
+   switchTagPair
+      (Stack $ \(Square sa) (Square sb) ->
+         Square (fh sa sb))
+      (Stack $ \(Wide ha wa) (Wide hb wb) ->
+         Wide (fh ha hb) (fw wa wb))
+      (Stack $ \(Tall ha wa) (Tall hb wb) ->
+         Tall (fh ha hb) (fw wa wb))
+      (Stack $ \(General ha wa) (General hb wb) ->
+         General (fh ha hb) (fw wa wb))
 
 
 
diff --git a/src/Numeric/LAPACK/Matrix/Hermitian.hs b/src/Numeric/LAPACK/Matrix/Hermitian.hs
--- a/src/Numeric/LAPACK/Matrix/Hermitian.hs
+++ b/src/Numeric/LAPACK/Matrix/Hermitian.hs
@@ -25,7 +25,7 @@
    sumRank1, sumRank1NonEmpty,
    sumRank2, sumRank2NonEmpty,
 
-   toSquare,
+   toSquare, fromSymmetric,
    gramian,            gramianAdjoint,
    congruenceDiagonal, congruenceDiagonalAdjoint,
    congruence,         congruenceAdjoint,
@@ -46,7 +46,7 @@
 import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
 import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
 import qualified Numeric.LAPACK.Matrix.Extent as Extent
-import Numeric.LAPACK.Matrix.Array.Triangular (Hermitian)
+import Numeric.LAPACK.Matrix.Array.Triangular (Hermitian, Symmetric)
 import Numeric.LAPACK.Matrix.Array (Full, General, Square)
 import Numeric.LAPACK.Matrix.Shape.Private (Order)
 import Numeric.LAPACK.Matrix.Modifier (Transposition(NonTransposed, Transposed))
@@ -186,9 +186,15 @@
 sumRank2NonEmpty order (NonEmpty.Cons xy xys) =
    sumRank2 order (Array.shape $ fst $ snd xy) (xy:xys)
 
+
 toSquare ::
    (Shape.C sh, Class.Floating a) => Hermitian sh a -> Square sh a
 toSquare = ArrMatrix.lift1 Basic.toSquare
+
+fromSymmetric :: (Shape.C sh, Class.Real a) => Symmetric sh a -> Hermitian sh a
+fromSymmetric =
+   ArrMatrix.lift1 $ Array.mapShape MatrixShape.hermitianFromSymmetric
+
 
 {- |
 gramian A = A^H * A
diff --git a/src/Numeric/LAPACK/Matrix/Hermitian/Basic.hs b/src/Numeric/LAPACK/Matrix/Hermitian/Basic.hs
--- a/src/Numeric/LAPACK/Matrix/Hermitian/Basic.hs
+++ b/src/Numeric/LAPACK/Matrix/Hermitian/Basic.hs
@@ -428,7 +428,7 @@
    \bPtr -> gramianIO order a bPtr $ gramianParameters order extent
 
 gramianParameters ::
-   (Extent.C horiz, Extent.C vert, Shape.C height, Shape.C width) =>
+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>
    Order ->
    Extent.Extent vert horiz height width ->
    ((Int, Int), (Char, Char, Int))
@@ -450,7 +450,7 @@
    \bPtr -> gramianIO order a bPtr $ gramianAdjointParameters order extent
 
 gramianAdjointParameters ::
-   (Extent.C horiz, Extent.C vert, Shape.C height, Shape.C width) =>
+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>
    Order ->
    Extent.Extent vert horiz height width ->
    ((Int, Int), (Char, Char, Int))
diff --git a/src/Numeric/LAPACK/Matrix/HermitianPositiveDefinite/Linear.hs b/src/Numeric/LAPACK/Matrix/HermitianPositiveDefinite/Linear.hs
--- a/src/Numeric/LAPACK/Matrix/HermitianPositiveDefinite/Linear.hs
+++ b/src/Numeric/LAPACK/Matrix/HermitianPositiveDefinite/Linear.hs
@@ -61,8 +61,7 @@
             LapackGen.pptrs uploPtr nPtr nrhsPtr apPtr xPtr ldxPtr
 
 
-inverse ::
-   (Shape.C sh, Class.Floating a) => Hermitian sh a -> Hermitian sh a
+inverse :: (Shape.C sh, Class.Floating a) => Hermitian sh a -> Hermitian sh a
 inverse
    (Array shape@(MatrixShape.Hermitian order sh) a) =
       Array.unsafeCreateWithSize shape $ \triSize bPtr -> do
@@ -75,8 +74,7 @@
          withInfo definiteMsg "pptrf" $ LapackGen.pptrf uploPtr nPtr bPtr
          withInfo rankMsg "pptri" $ LapackGen.pptri uploPtr nPtr bPtr
 
-decompose ::
-   (Shape.C sh, Class.Floating a) => Hermitian sh a -> Upper sh a
+decompose :: (Shape.C sh, Class.Floating a) => Hermitian sh a -> Upper sh a
 decompose
    (Array (MatrixShape.Hermitian order sh) a) =
       Array.unsafeCreateWithSize
@@ -91,8 +89,7 @@
          withInfo definiteMsg "pptrf" $ LapackGen.pptrf uploPtr nPtr bPtr
 
 
-determinant ::
-   (Shape.C sh, Class.Floating a) => Hermitian sh a -> RealOf a
+determinant :: (Shape.C sh, Class.Floating a) => Hermitian sh a -> RealOf a
 determinant =
    getDeterminant $
    Class.switchFloating
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
@@ -75,7 +75,7 @@
    {-# MINIMAL transposableSquare | fullSquare,squareFull #-}
    transposableSquare ::
       (Type.HeightOf typ ~ height, Eq height, Shape.C width,
-       Extent.C horiz, Extent.C vert, Class.Floating a) =>
+       Extent.C vert, Extent.C horiz, Class.Floating a) =>
       Transposition -> Matrix typ a ->
       Full vert horiz height width a -> Full vert horiz height width a
    transposableSquare NonTransposed a b = squareFull a b
@@ -84,14 +84,14 @@
 
    squareFull ::
       (Type.HeightOf typ ~ height, Eq height, Shape.C width,
-       Extent.C horiz, Extent.C vert, Class.Floating a) =>
+       Extent.C vert, Extent.C horiz, Class.Floating a) =>
       Matrix typ a ->
       Full vert horiz height width a -> Full vert horiz height width a
    squareFull = transposableSquare NonTransposed
 
    fullSquare ::
       (Type.WidthOf typ ~ width, Eq width, Shape.C height,
-       Extent.C horiz, Extent.C vert, Class.Floating a) =>
+       Extent.C vert, Extent.C horiz, Class.Floating a) =>
       Full vert horiz height width a ->
       Matrix typ a -> Full vert horiz height width a
    fullSquare b a =
@@ -102,14 +102,14 @@
 
 (#*##) ::
    (MultiplySquare typ, Type.HeightOf typ ~ height, Eq height, Shape.C width,
-    Extent.C horiz, Extent.C vert, Class.Floating a) =>
+    Extent.C vert, Extent.C horiz, Class.Floating a) =>
    Matrix typ a ->
    Full vert horiz height width a -> Full vert horiz height width a
 (#*##) = squareFull
 
 (##*#) ::
    (MultiplySquare typ, Type.WidthOf typ ~ width, Eq width, Shape.C height,
-    Extent.C horiz, Extent.C vert, Class.Floating a) =>
+    Extent.C vert, Extent.C horiz, Class.Floating a) =>
    Full vert horiz height width a ->
    Matrix typ a -> Full vert horiz height width a
 (##*#) = fullSquare
diff --git a/src/Numeric/LAPACK/Matrix/Plain/Divide.hs b/src/Numeric/LAPACK/Matrix/Plain/Divide.hs
--- a/src/Numeric/LAPACK/Matrix/Plain/Divide.hs
+++ b/src/Numeric/LAPACK/Matrix/Plain/Divide.hs
@@ -50,7 +50,7 @@
    {-# MINIMAL solve | solveLeft,solveRight #-}
    solve ::
       (Class.Floating a, Box.HeightOf shape ~ height, Eq height,
-       Extent.C horiz, Extent.C vert, Shape.C width) =>
+       Extent.C vert, Extent.C horiz, Shape.C width) =>
       Transposition -> Array shape a ->
       Full vert horiz height width a -> Full vert horiz height width a
    solve NonTransposed a b = solveRight a b
@@ -58,14 +58,14 @@
 
    solveRight ::
       (Class.Floating a, Box.HeightOf shape ~ height, Eq height,
-       Extent.C horiz, Extent.C vert, Shape.C width) =>
+       Extent.C vert, Extent.C horiz, Shape.C width) =>
       Array shape a ->
       Full vert horiz height width a -> Full vert horiz height width a
    solveRight = solve NonTransposed
 
    solveLeft ::
       (Class.Floating a, Box.HeightOf shape ~ width, Eq width,
-       Extent.C horiz, Extent.C vert, Shape.C height) =>
+       Extent.C vert, Extent.C horiz, Shape.C height) =>
       Full vert horiz height width a ->
       Array shape a ->
       Full vert horiz height width a
diff --git a/src/Numeric/LAPACK/Matrix/Plain/Multiply.hs b/src/Numeric/LAPACK/Matrix/Plain/Multiply.hs
--- a/src/Numeric/LAPACK/Matrix/Plain/Multiply.hs
+++ b/src/Numeric/LAPACK/Matrix/Plain/Multiply.hs
@@ -104,7 +104,7 @@
    {-# MINIMAL transposableSquare | fullSquare,squareFull #-}
    transposableSquare ::
       (Box.HeightOf shape ~ height, Eq height, Shape.C width,
-       Extent.C horiz, Extent.C vert, Class.Floating a) =>
+       Extent.C vert, Extent.C horiz, Class.Floating a) =>
       Transposition -> Array shape a ->
       Full vert horiz height width a -> Full vert horiz height width a
    transposableSquare NonTransposed a b = squareFull a b
@@ -112,14 +112,14 @@
 
    squareFull ::
       (Box.HeightOf shape ~ height, Eq height, Shape.C width,
-       Extent.C horiz, Extent.C vert, Class.Floating a) =>
+       Extent.C vert, Extent.C horiz, Class.Floating a) =>
       Array shape a ->
       Full vert horiz height width a -> Full vert horiz height width a
    squareFull = transposableSquare NonTransposed
 
    fullSquare ::
       (Box.WidthOf shape ~ width, Eq width, Shape.C height,
-       Extent.C horiz, Extent.C vert, Class.Floating a) =>
+       Extent.C vert, Extent.C horiz, Class.Floating a) =>
       Full vert horiz height width a ->
       Array shape a -> Full vert horiz height width a
    fullSquare = swapMultiply $ transposableSquare Transposed
diff --git a/src/Numeric/LAPACK/Matrix/Private.hs b/src/Numeric/LAPACK/Matrix/Private.hs
--- a/src/Numeric/LAPACK/Matrix/Private.hs
+++ b/src/Numeric/LAPACK/Matrix/Private.hs
@@ -69,7 +69,7 @@
 
 
 revealOrder ::
-   (Extent.C horiz, Extent.C vert) =>
+   (Extent.C vert, Extent.C horiz) =>
    Full vert horiz height width a ->
    Either (Array (height,width) a) (Array (width,height) a)
 revealOrder (Array (MatrixShape.Full order extent) a) =
diff --git a/src/Numeric/LAPACK/Matrix/Shape/Private.hs b/src/Numeric/LAPACK/Matrix/Shape/Private.hs
--- a/src/Numeric/LAPACK/Matrix/Shape/Private.hs
+++ b/src/Numeric/LAPACK/Matrix/Shape/Private.hs
@@ -363,6 +363,10 @@
    height = hermitianSize
    width = hermitianSize
 
+hermitianFromSymmetric :: Symmetric size -> Hermitian size
+hermitianFromSymmetric (Triangular _diag _uplo order size) =
+   Hermitian order size
+
 uploFromOrder :: Order -> Char
 uploFromOrder RowMajor = 'L'
 uploFromOrder ColumnMajor = 'U'
@@ -413,6 +417,10 @@
    type WidthOf (Triangular lo diag up size) = size
    height = triangularSize
    width = triangularSize
+
+symmetricFromHermitian :: Hermitian size -> Symmetric size
+symmetricFromHermitian (Hermitian order size) =
+   Triangular NonUnit autoUplo order size
 
 
 data Unit = Unit deriving (Eq, Show)
diff --git a/src/Numeric/LAPACK/Matrix/Square/Basic.hs b/src/Numeric/LAPACK/Matrix/Square/Basic.hs
--- a/src/Numeric/LAPACK/Matrix/Square/Basic.hs
+++ b/src/Numeric/LAPACK/Matrix/Square/Basic.hs
@@ -196,13 +196,7 @@
    Square sizeA a -> Full vert horiz sizeA sizeB a ->
    Full horiz vert sizeB sizeA a -> Square sizeB a ->
    Square (sizeA:+:sizeB) a
-stack a b c d =
-   fromGeneral $
-   Basic.above Basic.RightBias Extent.appendAny
-      (Basic.beside Basic.RightBias Extent.appendAny
-         (Matrix.fromFull a) (Matrix.fromFull b))
-      (Basic.beside Basic.RightBias Extent.appendAny
-         (Matrix.fromFull c) (Matrix.fromFull d))
+stack a b c d = Basic.stack a (Matrix.fromFull b) (Matrix.fromFull c) d
 
 
 square :: (Shape.C sh, Class.Floating a) => Square sh a -> Square sh a
diff --git a/src/Numeric/LAPACK/Matrix/Square/Linear.hs b/src/Numeric/LAPACK/Matrix/Square/Linear.hs
--- a/src/Numeric/LAPACK/Matrix/Square/Linear.hs
+++ b/src/Numeric/LAPACK/Matrix/Square/Linear.hs
@@ -14,7 +14,7 @@
 import Numeric.LAPACK.Matrix.Shape.Private (transposeFromOrder)
 import Numeric.LAPACK.Matrix.Private (Full, Square, argSquare)
 import Numeric.LAPACK.Private
-         (withAutoWorkspaceInfo, copyBlock, copyToTemp, copyToColumnMajor)
+         (withAutoWorkspaceInfo, copyBlock, copyToTemp, copyToColumnMajorTemp)
 
 import qualified Numeric.LAPACK.FFI.Generic as LapackGen
 import qualified Numeric.Netlib.Utility as Call
@@ -55,14 +55,12 @@
 _solve =
    argSquare $ \orderA shA a ->
    solver "Square.solve" shA $ \n nPtr nrhsPtr xPtr ldxPtr -> do
-      aPtr <- ContT $ withForeignPtr a
-      atmpPtr <- Call.allocaArray (n*n)
+      aPtr <- copyToColumnMajorTemp orderA n n a
       ldaPtr <- Call.leadingDim n
       ipivPtr <- Call.allocaArray n
       liftIO $ do
-         copyToColumnMajor orderA n n aPtr atmpPtr
          withInfo "gesv" $
-            LapackGen.gesv nPtr nrhsPtr atmpPtr ldaPtr ipivPtr xPtr ldxPtr
+            LapackGen.gesv nPtr nrhsPtr aPtr ldaPtr ipivPtr xPtr ldxPtr
 
 
 inverse :: (Shape.C sh, Class.Floating a) => Square sh a -> Square sh a
diff --git a/src/Numeric/LAPACK/Matrix/Symmetric.hs b/src/Numeric/LAPACK/Matrix/Symmetric.hs
--- a/src/Numeric/LAPACK/Matrix/Symmetric.hs
+++ b/src/Numeric/LAPACK/Matrix/Symmetric.hs
@@ -13,6 +13,7 @@
    split,
 
    toSquare,
+   fromHermitian,
 
    gramian,            gramianTransposed,
    congruenceDiagonal, congruenceDiagonalTransposed,
@@ -26,7 +27,7 @@
 import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
 import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
 import qualified Numeric.LAPACK.Matrix.Extent as Extent
-import Numeric.LAPACK.Matrix.Array.Triangular (Symmetric)
+import Numeric.LAPACK.Matrix.Array.Triangular (Symmetric, Hermitian)
 import Numeric.LAPACK.Matrix.Array (Full, General, Square)
 import Numeric.LAPACK.Matrix.Shape.Private (Order)
 import Numeric.LAPACK.Matrix.Private (ShapeInt)
@@ -35,6 +36,7 @@
 
 import qualified Numeric.Netlib.Class as Class
 
+import qualified Data.Array.Comfort.Storable as Array
 import qualified Data.Array.Comfort.Shape as Shape
 import Data.Array.Comfort.Shape ((:+:))
 
@@ -60,6 +62,10 @@
 
 toSquare :: (Shape.C sh, Class.Floating a) => Symmetric sh a -> Square sh a
 toSquare = Triangular.toSquare
+
+fromHermitian :: (Shape.C sh, Class.Real a) => Hermitian sh a -> Symmetric sh a
+fromHermitian =
+   ArrMatrix.lift1 $ Array.mapShape MatrixShape.symmetricFromHermitian
 
 
 identity :: (Shape.C sh, Class.Floating a) => Order -> sh -> Symmetric sh a
diff --git a/src/Numeric/LAPACK/Matrix/Symmetric/Basic.hs b/src/Numeric/LAPACK/Matrix/Symmetric/Basic.hs
--- a/src/Numeric/LAPACK/Matrix/Symmetric/Basic.hs
+++ b/src/Numeric/LAPACK/Matrix/Symmetric/Basic.hs
@@ -51,7 +51,7 @@
    \bPtr -> gramianIO order a bPtr $ gramianParameters order extent
 
 gramianParameters ::
-   (Extent.C horiz, Extent.C vert, Shape.C height, Shape.C width) =>
+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>
    Order ->
    Extent.Extent vert horiz height width ->
    ((Int, Int), (Char, Char, Int))
@@ -73,7 +73,7 @@
    \bPtr -> gramianIO order a bPtr $ gramianTransposedParameters order extent
 
 gramianTransposedParameters ::
-   (Extent.C horiz, Extent.C vert, Shape.C height, Shape.C width) =>
+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>
    Order ->
    Extent.Extent vert horiz height width ->
    ((Int, Int), (Char, Char, Int))
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
@@ -16,6 +16,7 @@
 
 import qualified Data.Array.Comfort.Shape as Shape
 
+import Data.Monoid (Monoid, mempty, mappend)
 import Data.Semigroup (Semigroup, (<>))
 
 
@@ -80,6 +81,22 @@
 instance (Shape.C sh, Eq sh) => MultiplySame (Perm.Permutation sh) where
    multiplySame (Permutation a) (Permutation b) =
       Permutation $ Perm.multiply b a
+
+
+instance
+   (MultiplySame typ, StaticIdentity typ, Class.Floating a) =>
+      Monoid (Matrix typ a) where
+   mappend = (<>)
+   mempty = staticIdentity
+
+class StaticIdentity typ where
+   staticIdentity :: (Class.Floating a) => Matrix typ a
+
+instance (Shape.Static shape) => StaticIdentity (Scale shape) where
+   staticIdentity = Scale Shape.static 1
+
+instance (Shape.Static sh) => StaticIdentity (Perm.Permutation sh) where
+   staticIdentity = Permutation $ Perm.identity Shape.static
 
 
 scaleWithCheck :: (Eq shape) =>
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
@@ -5,9 +5,15 @@
    leastSquaresMinimumNormRCond,
    pseudoInverseRCond,
 
+   project,
+   leastSquaresConstraint,
+   gaussMarkovLinearModel,
+
    determinant,
    determinantAbsolute,
    complement,
+   affineSpanFromKernel,
+   affineKernelFromSpan,
 
    householder,
    householderTall,
@@ -16,12 +22,17 @@
 import qualified Numeric.LAPACK.Orthogonal.Householder as HH
 import qualified Numeric.LAPACK.Orthogonal.Plain as Plain
 
+import qualified Numeric.LAPACK.Matrix.Multiply as Multiply
 import qualified Numeric.LAPACK.Matrix.Triangular as Triangular
+import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
 import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
 import qualified Numeric.LAPACK.Matrix.Basic as Basic
+import qualified Numeric.LAPACK.Matrix.Type as Matrix
 import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
-import Numeric.LAPACK.Matrix.Array (Full, Tall, Square)
+import qualified Numeric.LAPACK.Vector as Vector
+import Numeric.LAPACK.Matrix.Array (Full, General, Tall, Wide, Square)
 import Numeric.LAPACK.Matrix.Private (ShapeInt)
+import Numeric.LAPACK.Vector (Vector, (|-|))
 import Numeric.LAPACK.Scalar (RealOf)
 
 import qualified Numeric.Netlib.Class as Class
@@ -94,6 +105,65 @@
    Basic.uncheck . ArrMatrix.toVector
 
 
+{- |
+@project b d x@ projects @x@ on the plane described by @B*x = d@.
+
+@b@ must have full rank.
+-}
+project ::
+   (Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>
+   Wide height width a -> Vector height a ->
+   Vector width a -> Vector width a
+project b d x =
+   x
+   |-|
+   ArrMatrix.unliftColumn MatrixShape.ColumnMajor
+      (minimumNorm b) (Multiply.matrixVector b x |-| d)
+
+
+{- |
+@leastSquaresConstraint a c b d@ computes @x@
+with minimal @|| c - A*x ||_2@ and constraint @B*x = d@.
+
+@b@ must be wide and @a===b@ must be tall
+and both matrices must have full rank.
+-}
+leastSquaresConstraint ::
+   (Shape.C height, Eq height,
+    Shape.C width, Eq width,
+    Shape.C constraints, Eq constraints, Class.Floating a) =>
+   General height width a -> Vector height a ->
+   Wide constraints width a -> Vector constraints a ->
+   Vector width a
+leastSquaresConstraint a c b d =
+   Plain.leastSquaresConstraint
+      (ArrMatrix.toVector a) c
+      (ArrMatrix.toVector b) d
+
+{- |
+@gaussMarkovLinearModel a b d@ computes @(x,y)@
+with minimal @|| y ||_2@ and constraint @d = A*x + B*y@.
+
+@a@ must be tall and @a|||b@ must be wide
+and both matrices must have full rank.
+-}
+gaussMarkovLinearModel ::
+   (Shape.C height, Eq height,
+    Shape.C width, Eq width,
+    Shape.C opt, Eq opt, Class.Floating a) =>
+   Tall height width a -> General height opt a -> Vector height a ->
+   (Vector width a, Vector opt a)
+gaussMarkovLinearModel a b d =
+   {-
+   Fortran-LAPACK and OpenBLAS would leave Y uninitalized
+   instead of setting Y to zeros.
+   -}
+   if Shape.size (Matrix.height a) == 0
+      then (Vector.zero (Matrix.width a), Vector.zero (Matrix.width b))
+      else Plain.gaussMarkovLinearModel
+               (ArrMatrix.toVector a) (ArrMatrix.toVector b) d
+
+
 {-
 @(q,r) = householder a@
 means that @q@ is unitary and @r@ is upper triangular and @a = multiply q r@.
@@ -145,3 +215,41 @@
    (Shape.C height, Shape.C width, Class.Floating a) =>
    Tall height width a -> Tall height ShapeInt a
 complement = ArrMatrix.lift1 Plain.complement
+
+
+{- |
+> affineSpanFromKernel a b == (c,d)
+
+Means:
+An affine subspace is given implicitly by {x : a#*|x == b}.
+Convert this into an explicit representation {c#*|y|+|d : y}.
+Matrix @a@ must have full rank,
+otherwise the explicit representation will miss dimensions
+and we cannot easily determine the origin @d@ as a minimum norm solution.
+-}
+affineSpanFromKernel ::
+   (Shape.C width, Eq width, Shape.C height, Eq height, Class.Floating a) =>
+   Wide height width a -> Vector height a ->
+   (Tall width ShapeInt a, Vector width a)
+affineSpanFromKernel a b =
+   let qr = HH.fromMatrix $ ArrMatrix.lift1 Basic.adjoint a
+   in (ArrMatrix.lift0 $ Plain.extractComplement qr,
+       ArrMatrix.unliftColumn MatrixShape.ColumnMajor (HH.minimumNorm qr) b)
+
+{- |
+This conversion is somehow inverse to 'affineSpanFromKernel'.
+However, it is not precisely inverse in either direction.
+This is because both 'affineSpanFromKernel' and 'affineKernelFromSpan'
+accept non-orthogonal matrices but always return orthogonal ones.
+
+In @affineKernelFromSpan c d@,
+matrix @c@ should have full rank,
+otherwise the implicit representation will miss dimensions.
+-}
+affineKernelFromSpan ::
+   (Shape.C width, Eq width, Shape.C height, Eq height, Class.Floating a) =>
+   Tall height width a -> Vector height a ->
+   (Wide ShapeInt height a, Vector ShapeInt a)
+affineKernelFromSpan c d =
+   let a = Basic.adjoint $ Plain.complement $ ArrMatrix.toVector c
+   in (ArrMatrix.lift0 a, Basic.multiplyVector a d)
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
@@ -5,8 +5,12 @@
    leastSquaresMinimumNormRCond,
    pseudoInverseRCond,
 
+   leastSquaresConstraint,
+   gaussMarkovLinearModel,
+
    determinantAbsolute,
    complement,
+   extractComplement,
    ) where
 
 import qualified Numeric.LAPACK.Orthogonal.Private as HH
@@ -17,11 +21,14 @@
 import qualified Numeric.LAPACK.Matrix.Basic as Basic
 import qualified Numeric.LAPACK.Vector as Vector
 import Numeric.LAPACK.Matrix.Shape.Private (Order(RowMajor,ColumnMajor))
-import Numeric.LAPACK.Matrix.Private (Full, Tall, ShapeInt, shapeInt)
+import Numeric.LAPACK.Matrix.Private
+         (Full, General, Tall, Wide, ShapeInt, shapeInt)
+import Numeric.LAPACK.Vector (Vector)
 import Numeric.LAPACK.Scalar (RealOf, zero, absolute)
 import Numeric.LAPACK.Private
          (lacgv, peekCInt,
-          copySubMatrix, copyToTemp, copyToColumnMajor, copyToSubColumnMajor,
+          copySubMatrix, copyToTemp,
+          copyToColumnMajorTemp, copyToSubColumnMajor,
           withAutoWorkspaceInfo, rankMsg, errorCodeMsg, createHigherArray)
 
 import qualified Numeric.LAPACK.FFI.Generic as LapackGen
@@ -30,6 +37,7 @@
 import qualified Numeric.Netlib.Utility as Call
 import qualified Numeric.Netlib.Class as Class
 
+import qualified Data.Array.Comfort.Storable.Unchecked.Monadic as ArrayIO
 import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Shape as Shape
 import Data.Array.Comfort.Storable.Unchecked (Array(Array))
@@ -56,7 +64,7 @@
    Full vert horiz width nrhs a
 leastSquares
    (Array shapeA@(MatrixShape.Full orderA extentA) a)
-   (Array shapeB@(MatrixShape.Full orderB extentB) b) =
+   (Array        (MatrixShape.Full orderB extentB) b) =
 
  case Extent.fuse (Extent.generalizeWide $ Extent.transpose extentA) extentB of
   Nothing -> error "leastSquares: height shapes mismatch"
@@ -76,15 +84,12 @@
       nrhsPtr <- Call.cint nrhs
       (transPtr,aPtr) <- adjointA orderA (m*n) a
       ldaPtr <- Call.leadingDim lda
-      bPtr <- ContT $ withForeignPtr b
       ldbPtr <- Call.leadingDim ldb
-      let bSize = Shape.size shapeB
-      btmpPtr <- Call.allocaArray bSize
-      liftIO $ copyToColumnMajor orderB ldb nrhs bPtr btmpPtr
+      bPtr <- copyToColumnMajorTemp orderB ldb nrhs b
       liftIO $ withAutoWorkspaceInfo rankMsg "gels" $
          LapackGen.gels transPtr
-            mPtr nPtr nrhsPtr aPtr ldaPtr btmpPtr ldbPtr
-      liftIO $ copySubMatrix ldx nrhs ldb btmpPtr ldx xPtr
+            mPtr nPtr nrhsPtr aPtr ldaPtr bPtr ldbPtr
+      liftIO $ copySubMatrix ldx nrhs ldb bPtr ldx xPtr
 
 minimumNorm ::
    (Extent.C vert, Extent.C horiz,
@@ -184,12 +189,9 @@
 leastSquaresMinimumNormIO rcond shapeX orderA a orderB b m n nrhs =
    createHigherArray shapeX m n nrhs $ \(tmpPtr,ldtmp) -> do
 
-   let aSize = m*n
    let lda = m
    evalContT $ do
-      aPtr <- ContT $ withForeignPtr a
-      atmpPtr <- Call.allocaArray aSize
-      liftIO $ copyToColumnMajor orderA m n aPtr atmpPtr
+      aPtr <- copyToColumnMajorTemp orderA m n a
       ldaPtr <- Call.leadingDim lda
       ldtmpPtr <- Call.leadingDim ldtmp
       bPtr <- ContT $ withForeignPtr b
@@ -197,7 +199,7 @@
       jpvtPtr <- Call.allocaArray n
       liftIO $ pokeArray jpvtPtr (replicate n 0)
       rankPtr <- Call.alloca
-      gelsy m n nrhs atmpPtr ldaPtr tmpPtr ldtmpPtr jpvtPtr rcond rankPtr
+      gelsy m n nrhs aPtr ldaPtr tmpPtr ldtmpPtr jpvtPtr rcond rankPtr
       liftIO $ peekCInt rankPtr
 
 
@@ -259,6 +261,82 @@
          MatrixShape.fullHeight $ Array.shape a
 
 
+leastSquaresConstraint ::
+   (Shape.C height, Eq height,
+    Shape.C width, Eq width,
+    Shape.C constraints, Eq constraints, Class.Floating a) =>
+   General height width a -> Vector height a ->
+   Wide constraints width a -> Vector constraints a ->
+   Vector width a
+leastSquaresConstraint
+   (Array (MatrixShape.Full orderA extentA) a) c
+   (Array (MatrixShape.Full orderB extentB) b) d =
+
+ let sameShape name shape0 shape1 =
+      if shape0 == shape1
+         then shape0
+         else error $ "leastSquaresConstraint: " ++ name ++ " shapes mismatch"
+     width = sameShape "width" (Extent.width extentA) (Extent.width extentB)
+ in
+   Array.unsafeCreate width $ \xPtr -> do
+
+   let height = sameShape "height" (Extent.height extentA) (Array.shape c)
+   let constraints =
+         sameShape "constraints" (Extent.height extentB) (Array.shape d)
+   let m = Shape.size height
+   let n = Shape.size width
+   let p = Shape.size constraints
+   evalContT $ do
+      mPtr <- Call.cint m
+      nPtr <- Call.cint n
+      pPtr <- Call.cint p
+      aPtr <- copyToColumnMajorTemp orderA m n a
+      ldaPtr <- Call.leadingDim m
+      bPtr <- copyToColumnMajorTemp orderB p n b
+      ldbPtr <- Call.leadingDim p
+      cPtr <- copyToTemp m (Array.buffer c)
+      dPtr <- copyToTemp p (Array.buffer d)
+      liftIO $ withAutoWorkspaceInfo rankMsg "gglse" $
+         LapackGen.gglse
+            mPtr nPtr pPtr aPtr ldaPtr bPtr ldbPtr cPtr dPtr xPtr
+
+gaussMarkovLinearModel ::
+   (Shape.C height, Eq height,
+    Shape.C width, Eq width,
+    Shape.C opt, Eq opt, Class.Floating a) =>
+   Tall height width a -> General height opt a -> Vector height a ->
+   (Vector width a, Vector opt a)
+gaussMarkovLinearModel
+   (Array (MatrixShape.Full orderA extentA) a)
+   (Array (MatrixShape.Full orderB extentB) b) d =
+
+   let width = Extent.width extentA in
+   let opt = Extent.width extentB in
+   Array.unsafeCreateWithSizeAndResult width $ \m xPtr -> do
+   ArrayIO.unsafeCreateWithSize opt $ \p yPtr -> do
+
+   let sameHeight shape0 shape1 =
+         if shape0 == shape1
+            then shape0
+            else error $ "gaussMarkovLinearModel: height shapes mismatch"
+       height =
+         sameHeight (Array.shape d) $
+         sameHeight (Extent.height extentA) (Extent.height extentB)
+   let n = Shape.size height
+   evalContT $ do
+      mPtr <- Call.cint m
+      nPtr <- Call.cint n
+      pPtr <- Call.cint p
+      aPtr <- copyToColumnMajorTemp orderA n m a
+      ldaPtr <- Call.leadingDim n
+      bPtr <- copyToColumnMajorTemp orderB n p b
+      ldbPtr <- Call.leadingDim n
+      dPtr <- copyToTemp n (Array.buffer d)
+      liftIO $ withAutoWorkspaceInfo rankMsg "ggglm" $
+         LapackGen.ggglm
+            nPtr mPtr pPtr aPtr ldaPtr bPtr ldbPtr dPtr xPtr yPtr
+
+
 determinantAbsolute ::
    (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width,
     Class.Floating a) =>
@@ -272,7 +350,13 @@
 complement ::
    (Shape.C height, Shape.C width, Class.Floating a) =>
    Tall height width a -> Tall height ShapeInt a
-complement a =
-   Basic.dropColumns (Shape.size $ MatrixShape.fullWidth $ Array.shape a) $
+complement = extractComplement . HH.fromMatrix
+
+extractComplement ::
+   (Shape.C height, Shape.C width, Class.Floating a) =>
+   HH.Tall height width a -> Tall height ShapeInt a
+extractComplement qr =
+   Basic.dropColumns
+      (Shape.size $ MatrixShape.splitWidth $ Array.shape $ HH.split_ qr) $
    Basic.mapWidth (shapeInt . Shape.size) $ Square.toFull $
-   HH.extractQ $ HH.fromMatrix a
+   HH.extractQ qr
diff --git a/src/Numeric/LAPACK/Orthogonal/Private.hs b/src/Numeric/LAPACK/Orthogonal/Private.hs
--- a/src/Numeric/LAPACK/Orthogonal/Private.hs
+++ b/src/Numeric/LAPACK/Orthogonal/Private.hs
@@ -474,9 +474,9 @@
    (vert ~ Extent.Small, horiz ~ Extent.Small,
     Shape.C height, height ~ width) =>
       Divide.Solve (Hh vert horiz height width) where
-   solveRight = ArrMatrix.lift1 . leastSquares . mapExtent Extent.generalizeWide
+   solveRight = ArrMatrix.lift1 . leastSquares . mapExtent Extent.fromSquare
    solveLeft =
       flip $ \a -> ArrMatrix.lift1 $
          Basic.adjoint .
-         minimumNorm (mapExtent Extent.generalizeWide a) .
+         minimumNorm (mapExtent Extent.fromSquare a) .
          Basic.adjoint
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
@@ -44,6 +44,27 @@
 import Prelude hiding (odd)
 
 
+{- $setup
+>>> import qualified Test.QuickCheck as QC
+>>> import Test.Permutation (genPerm, genPivots)
+>>>
+>>> import qualified Numeric.LAPACK.Permutation as Perm
+>>> import Numeric.LAPACK.Permutation (Permutation, Inversion(NonInverted), determinant, multiply, transpose)
+>>> import Numeric.LAPACK.Matrix (ShapeInt)
+>>>
+>>> import qualified Data.Array.Comfort.Storable as Array
+>>> import Data.Eq.HT (equating)
+>>> import Data.Semigroup ((<>))
+>>>
+>>> import Control.Applicative (liftA2)
+>>>
+>>> genPerm2 :: QC.Gen (Permutation ShapeInt, Permutation ShapeInt)
+>>> genPerm2 = do
+>>>    nat <- QC.arbitrary
+>>>    liftA2 (,) (genPerm nat) (genPerm nat)
+-}
+
+
 newtype Permutation sh = Permutation (Vector (Shape sh) (Element sh))
    deriving (Show)
 
@@ -62,6 +83,9 @@
 identity :: (Shape.C sh) => sh -> Permutation sh
 identity shape = Permutation $ CheckedArray.sample (Shape shape) id
 
+{- |
+prop> QC.forAll QC.arbitraryBoundedEnum $ \inv -> QC.forAll (QC.arbitrary >>= genPivots) $ \xs -> Array.toList xs == Array.toList (Perm.toPivots inv (Perm.fromPivots inv xs))
+-}
 fromPivots ::
    (Shape.C sh) =>
    Inversion -> Vector (Shape sh) (Element sh) -> Permutation sh
@@ -140,6 +164,9 @@
 {-
 We could also count the cycles of even number. This might be a little faster.
 -}
+{- |
+prop> QC.forAll genPerm2 $ \(p0,p1) -> determinant (multiply p0 p1) == determinant p0 <> determinant p1
+-}
 determinant :: (Shape.C sh) => Permutation sh -> Sign
 determinant =
    (\oddp -> if oddp then Negative else Positive) .
@@ -166,6 +193,9 @@
 dropEven xs = xs
 
 
+{- |
+prop> QC.forAll genPerm2 $ \(p0,p1) -> equating (Array.toList . Perm.toPivots NonInverted) (transpose $ multiply p0 p1) (multiply (transpose p1) (transpose p0))
+-}
 transpose :: (Shape.C sh) => Permutation sh -> Permutation sh
 transpose (Permutation perm) =
    Permutation $ runST (MutArray.unsafeFreeze =<< transposeToMutable perm)
diff --git a/src/Numeric/LAPACK/Private.hs b/src/Numeric/LAPACK/Private.hs
--- a/src/Numeric/LAPACK/Private.hs
+++ b/src/Numeric/LAPACK/Private.hs
@@ -195,6 +195,16 @@
            else copySubMatrix m n m aPtr ldb bPtr
 
 
+copyToColumnMajorTemp ::
+   (Class.Floating a) =>
+   Order -> Int -> Int -> ForeignPtr a -> ContT r IO (Ptr a)
+copyToColumnMajorTemp order m n fptr = do
+   ptr <- ContT $ withForeignPtr fptr
+   tmpPtr <- Call.allocaArray (m*n)
+   liftIO $ copyToColumnMajor order m n ptr tmpPtr
+   return tmpPtr
+
+
 pointerSeq :: (Storable a) => Int -> Ptr a -> [Ptr a]
 pointerSeq k ptr = iterate (flip advancePtr k) ptr
 
diff --git a/src/Numeric/LAPACK/Singular/Plain.hs b/src/Numeric/LAPACK/Singular/Plain.hs
--- a/src/Numeric/LAPACK/Singular/Plain.hs
+++ b/src/Numeric/LAPACK/Singular/Plain.hs
@@ -30,7 +30,7 @@
 import Numeric.LAPACK.Scalar (RealOf, zero)
 import Numeric.LAPACK.Private
          (withAutoWorkspace, peekCInt, createHigherArray,
-          copyToTemp, copyToColumnMajor, copyToSubColumnMajor)
+          copyToTemp, copyToColumnMajorTemp, copyToSubColumnMajor)
 
 import qualified Numeric.LAPACK.FFI.Complex as LapackComplex
 import qualified Numeric.LAPACK.FFI.Real as LapackReal
@@ -325,15 +325,12 @@
    createHigherArray shapeX m n nrhs $ \(tmpPtr,ldtmp) -> do
 
    let mn = min m n
-   let aSize = m*n
    let lda = m
    evalContT $ do
       mPtr <- Call.cint m
       nPtr <- Call.cint n
       nrhsPtr <- Call.cint nrhs
-      aPtr <- Call.allocaArray aSize
-      liftIO $ withForeignPtr a $ \asrcPtr ->
-         copyToColumnMajor orderA m n asrcPtr aPtr
+      aPtr <- copyToColumnMajorTemp orderA m n a
       ldaPtr <- Call.leadingDim lda
       ldtmpPtr <- Call.leadingDim ldtmp
       liftIO $ withForeignPtr b $ \bPtr ->
diff --git a/test-module.list b/test-module.list
new file mode 100644
--- /dev/null
+++ b/test-module.list
@@ -0,0 +1,3 @@
+Numeric.LAPACK.Permutation.Private
+Numeric.LAPACK.Example.DividedDifference
+Numeric.LAPACK.Example.EconomicAllocation
diff --git a/test/DocTest/Main.hs b/test/DocTest/Main.hs
new file mode 100644
--- /dev/null
+++ b/test/DocTest/Main.hs
@@ -0,0 +1,14 @@
+-- Do not edit! Automatically created with doctest-extract.
+module DocTest.Main where
+
+import qualified DocTest.Numeric.LAPACK.Permutation.Private
+import qualified DocTest.Numeric.LAPACK.Example.DividedDifference
+import qualified DocTest.Numeric.LAPACK.Example.EconomicAllocation
+
+import qualified Test.DocTest.Driver as DocTest
+
+main :: DocTest.T ()
+main = do
+    DocTest.Numeric.LAPACK.Permutation.Private.test
+    DocTest.Numeric.LAPACK.Example.DividedDifference.test
+    DocTest.Numeric.LAPACK.Example.EconomicAllocation.test
diff --git a/test/DocTest/Numeric/LAPACK/Example/DividedDifference.hs b/test/DocTest/Numeric/LAPACK/Example/DividedDifference.hs
new file mode 100644
--- /dev/null
+++ b/test/DocTest/Numeric/LAPACK/Example/DividedDifference.hs
@@ -0,0 +1,47 @@
+-- Do not edit! Automatically created with doctest-extract from src/Numeric/LAPACK/Example/DividedDifference.hs
+{-# LINE 29 "src/Numeric/LAPACK/Example/DividedDifference.hs" #-}
+
+module DocTest.Numeric.LAPACK.Example.DividedDifference where
+
+import qualified Test.DocTest.Driver as DocTest
+
+{-# LINE 30 "src/Numeric/LAPACK/Example/DividedDifference.hs" #-}
+import     qualified Test.Utility as Util
+import     Test.Utility (approxArray)
+
+import     qualified Numeric.LAPACK.Vector as Vector
+import     Numeric.LAPACK.Example.DividedDifference (dividedDifferencesMatrix)
+import     Numeric.LAPACK.Matrix (ShapeInt, (#+#))
+import     Numeric.LAPACK.Vector ((|+|))
+
+import     qualified Data.Array.Comfort.Storable as Array
+
+import     qualified Test.QuickCheck as QC
+
+import     Control.Monad (liftM2)
+import     Data.Tuple.HT (mapPair)
+import     Data.Semigroup ((<>))
+
+type     Vector = Vector.Vector ShapeInt Float
+
+genDD     :: QC.Gen (Vector, (Vector, Vector))
+genDD     = do
+       (ys0,ys1) <-
+          fmap (mapPair (Vector.autoFromList, Vector.autoFromList) .
+                unzip . take 10) $
+          QC.listOf $ liftM2 (,) (Util.genElement 10) (Util.genElement 10)
+       xs <- Util.genDistinct 10 10 $ Array.shape ys0
+       return (xs,(ys0,ys1))
+
+test :: DocTest.T ()
+test = do
+ DocTest.printPrefix "Numeric.LAPACK.Example.DividedDifference:77: "
+{-# LINE 77 "src/Numeric/LAPACK/Example/DividedDifference.hs" #-}
+ DocTest.property
+{-# LINE 77 "src/Numeric/LAPACK/Example/DividedDifference.hs" #-}
+     (QC.forAll genDD $ \(xs, (ys0,ys1)) -> approxArray (dividedDifferencesMatrix xs (ys0|+|ys1)) (dividedDifferencesMatrix xs ys0 #+# dividedDifferencesMatrix xs ys1))
+ DocTest.printPrefix "Numeric.LAPACK.Example.DividedDifference:78: "
+{-# LINE 78 "src/Numeric/LAPACK/Example/DividedDifference.hs" #-}
+ DocTest.property
+{-# LINE 78 "src/Numeric/LAPACK/Example/DividedDifference.hs" #-}
+     (QC.forAll genDD $ \(xs, (ys0,ys1)) -> approxArray (dividedDifferencesMatrix xs (Vector.mul ys0 ys1)) (dividedDifferencesMatrix xs ys0 <> dividedDifferencesMatrix xs ys1))
diff --git a/test/DocTest/Numeric/LAPACK/Example/EconomicAllocation.hs b/test/DocTest/Numeric/LAPACK/Example/EconomicAllocation.hs
new file mode 100644
--- /dev/null
+++ b/test/DocTest/Numeric/LAPACK/Example/EconomicAllocation.hs
@@ -0,0 +1,23 @@
+-- Do not edit! Automatically created with doctest-extract from src/Numeric/LAPACK/Example/EconomicAllocation.hs
+{-# LINE 20 "src/Numeric/LAPACK/Example/EconomicAllocation.hs" #-}
+
+module DocTest.Numeric.LAPACK.Example.EconomicAllocation where
+
+import qualified Test.DocTest.Driver as DocTest
+
+{-# LINE 21 "src/Numeric/LAPACK/Example/EconomicAllocation.hs" #-}
+import     Numeric.LAPACK.Example.EconomicAllocation
+import     Test.Utility (approxVector)
+
+import     qualified Numeric.LAPACK.Vector as Vector
+import     Numeric.LAPACK.Vector ((+++))
+
+import     qualified Data.Array.Comfort.Storable as Array
+
+test :: DocTest.T ()
+test = do
+ DocTest.printPrefix "Numeric.LAPACK.Example.EconomicAllocation:90: "
+{-# LINE 90 "src/Numeric/LAPACK/Example/EconomicAllocation.hs" #-}
+ DocTest.property
+{-# LINE 90 "src/Numeric/LAPACK/Example/EconomicAllocation.hs" #-}
+     (let result = iterated expenses0 balances0 in approxVector result $ compensated expenses0 balances0 +++ Vector.zero (Array.shape $ Vector.takeRight result))
diff --git a/test/DocTest/Numeric/LAPACK/Permutation/Private.hs b/test/DocTest/Numeric/LAPACK/Permutation/Private.hs
new file mode 100644
--- /dev/null
+++ b/test/DocTest/Numeric/LAPACK/Permutation/Private.hs
@@ -0,0 +1,43 @@
+-- Do not edit! Automatically created with doctest-extract from 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 48 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
+import     qualified Test.QuickCheck as QC
+import     Test.Permutation (genPerm, genPivots)
+
+import     qualified Numeric.LAPACK.Permutation as Perm
+import     Numeric.LAPACK.Permutation (Permutation, Inversion(NonInverted), determinant, multiply, transpose)
+import     Numeric.LAPACK.Matrix (ShapeInt)
+
+import     qualified Data.Array.Comfort.Storable as Array
+import     Data.Eq.HT (equating)
+import     Data.Semigroup ((<>))
+
+import     Control.Applicative (liftA2)
+
+genPerm2     :: QC.Gen (Permutation ShapeInt, Permutation ShapeInt)
+genPerm2     = do
+       nat <- QC.arbitrary
+       liftA2 (,) (genPerm nat) (genPerm nat)
+
+test :: DocTest.T ()
+test = do
+ DocTest.printPrefix "Numeric.LAPACK.Permutation.Private:87: "
+{-# LINE 87 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
+ DocTest.property
+{-# LINE 87 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
+     (QC.forAll QC.arbitraryBoundedEnum $ \inv -> QC.forAll (QC.arbitrary >>= genPivots) $ \xs -> Array.toList xs == Array.toList (Perm.toPivots inv (Perm.fromPivots inv xs)))
+ DocTest.printPrefix "Numeric.LAPACK.Permutation.Private:168: "
+{-# LINE 168 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
+ DocTest.property
+{-# LINE 168 "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:197: "
+{-# LINE 197 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
+ DocTest.property
+{-# LINE 197 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
+     (QC.forAll genPerm2 $ \(p0,p1) -> equating (Array.toList . Perm.toPivots NonInverted) (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
@@ -14,11 +14,12 @@
 import qualified Test.Singular as Singular
 import qualified Test.Shape as Shape
 import qualified Test.Permutation as Permutation
-import qualified Test.Example as Example
+import qualified DocTest.Main as DocTestMain
 import Test.Format ()
 import Test.Utility (Tagged(Tagged), prefix)
 
 import qualified Test.QuickCheck as QC
+import qualified Test.DocTest.Driver as DocTest
 
 import Numeric.LAPACK.Scalar (RealOf)
 
@@ -26,6 +27,15 @@
 
 import Type.Base.Proxy (Proxy(Proxy))
 
+import System.Exit (exitFailure)
+
+import Text.Printf (printf)
+
+import qualified Control.Monad.Trans.Writer.Strict as MW
+import qualified Control.Monad.Trans.Reader as MR
+import Control.Monad.IO.Class (liftIO)
+import Control.Monad (when, void)
+
 import qualified Data.List as List
 import Data.Complex (Complex)
 import Data.Tuple.HT (mapSnd)
@@ -68,14 +78,28 @@
 simpleTests :: [(String, QC.Property)]
 simpleTests =
    prefix "Shape" Shape.tests ++
-   prefix "Permutation" Permutation.tests ++
-   prefix "Example" Example.tests ++
    []
 
+
+
+run :: MW.WriterT DocTest.Count IO () -> IO ()
+run act = do
+   count <- MW.execWriterT act
+   putStrLn ""
+   void $ printf "Total: %d\n" $ DocTest.numTotal count
+   void $ printf "Failures: %d\n" $ DocTest.numFailures count
+   when (DocTest.numFailures count > 0) exitFailure
+
 main :: IO ()
 main =
-   mapM_ (\(name,act) -> putStr (name ++ ": ") >> act) $
+   run $
+   (>> MR.runReaderT DocTestMain.main QC.stdArgs) $
 
-   map (mapSnd (QC.quickCheckWith (QC.stdArgs {QC.maxSuccess=200}))) tests
+   mapM_
+      (\(name,(args,act)) -> do
+         liftIO (putStr (name ++ ": "))
+         MR.runReaderT (DocTest.property act) args) $
+
+   map (mapSnd ((,) (QC.stdArgs {QC.maxSuccess=200}))) tests
    ++
-   map (mapSnd QC.quickCheck) simpleTests
+   map (mapSnd ((,) QC.stdArgs)) simpleTests
diff --git a/test/Test/Example.hs b/test/Test/Example.hs
deleted file mode 100644
--- a/test/Test/Example.hs
+++ /dev/null
@@ -1,21 +0,0 @@
-module Test.Example where
-
-import qualified Test.Utility as Util
-
-import qualified Numeric.LAPACK.Example.EconomicAllocation as Eco
-import qualified Numeric.LAPACK.Vector as Vector
-
-import qualified Data.Array.Comfort.Storable as Array
-
-import qualified Test.QuickCheck as QC
-
-
-tests :: [(String, QC.Property)]
-tests =
-   ("economicAllocation",
-      QC.property $
-         let iterated = Eco.iterated Eco.expenses0 Eco.balances0
-         in Util.approxVector iterated
-               (Vector.append (Eco.compensated Eco.expenses0 Eco.balances0) $
-                Vector.zero $ Array.shape $ Vector.takeRight iterated)) :
-   []
diff --git a/test/Test/Generator.hs b/test/Test/Generator.hs
--- a/test/Test/Generator.hs
+++ b/test/Test/Generator.hs
@@ -429,6 +429,24 @@
       Util.genArrayExtraDiag maxElem shape
          (const $ fromReal <$> Util.genReal maxElem)
 
+lscStack ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   Matrix (ShapeInt:+:ShapeInt) ShapeInt a
+      (Matrix.Tall (ShapeInt:+:ShapeInt) ShapeInt a)
+lscStack =
+   condition
+      (Util.fullRankTall . Matrix.transpose .
+       Matrix.wideFromGeneral . Matrix.takeTop . Matrix.fromFull) $
+   condition Util.fullRankTall $
+   mapGen Util.genArray $
+   shapeFromDims (uncurry . MatrixShape.tall) $
+   Base $ do
+      height <- newVariable
+      width <- newVariableWith $ Logic.between height
+      return
+         ((height,width),
+          liftA2 (,) (Logic.query height) (Logic.query width))
+
 
 {-
 There cannot be a pure/point function.
diff --git a/test/Test/Logic.hs b/test/Test/Logic.hs
--- a/test/Test/Logic.hs
+++ b/test/Test/Logic.hs
@@ -87,6 +87,31 @@
    assignmentM (liftZeroBased $ \a -> choose (\ maxk -> (a,maxk))) va vb
    assignmentM (liftZeroBased $ \b -> choose (\_maxk -> (0,b))) vb va
 
+{-
+We cannot split this into something like
+
+a <!= c <> c <!= a+b
+
+because the solver might choose in this order
+
+a+b=6:+:4
+c=5
+
+and then it is left with a>c from which it cannot recover.
+-}
+between :: Variable s (ShapeInt:+:ShapeInt) -> Variable s ShapeInt -> System s
+between vab vc  =  AppMn.Cons $ do
+   assignmentM
+      (\(Shape.ZeroBased a :+: Shape.ZeroBased b) ->
+         Shape.ZeroBased <$> choose (\_maxk -> (a,a+b)))
+      vab vc
+   assignmentM
+      (\(Shape.ZeroBased c) ->
+         liftA2 (:+:)
+            (Shape.ZeroBased <$> choose (\_maxk -> (0,c)))
+            (Shape.ZeroBased <$> choose (\ maxk -> (c,maxk))))
+      vc vab
+
 
 class (Shape.C dim) => Dim dim where chooseDim :: M s dim
 instance (i ~ Int) => Dim (Shape.ZeroBased i) where
diff --git a/test/Test/Orthogonal.hs b/test/Test/Orthogonal.hs
--- a/test/Test/Orthogonal.hs
+++ b/test/Test/Orthogonal.hs
@@ -1,13 +1,14 @@
 {-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
 {-# LANGUAGE FlexibleContexts #-}
 module Test.Orthogonal (testsVar) where
 
 import qualified Test.Divide as Divide
 import qualified Test.Generator as Gen
 import qualified Test.Utility as Util
-import Test.Generator ((<#*#>), (<#\#>))
+import Test.Generator ((<#*#>), (<#\#>), (<-*#>), (<#*|>), (<|=|>))
 import Test.Utility
-         (approx, approxReal, approxArrayTol, approxMatrix,
+         (approx, approxReal, approxArrayTol, approxMatrix, approxVectorTol,
           Tagged, isIdentity, isUnitary, maybeConjugate)
 
 import qualified Numeric.LAPACK.Orthogonal.Householder as HH
@@ -20,13 +21,16 @@
 import qualified Numeric.LAPACK.Matrix as Matrix
 import qualified Numeric.LAPACK.Vector as Vector
 import Numeric.LAPACK.Matrix.Square (Square)
-import Numeric.LAPACK.Matrix (General, ShapeInt, (#*#), (##*#), (#*##), (#\##))
+import Numeric.LAPACK.Matrix
+         (General, ShapeInt, (#*#), (##*#), (#*##), (#\##), (#*|))
+import Numeric.LAPACK.Vector (Vector, (|+|), (|-|))
 import Numeric.LAPACK.Scalar (RealOf, absolute, selectReal)
 
 import qualified Numeric.Netlib.Class as Class
 
 import qualified Data.Array.Comfort.Storable as Array
 import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Shape ((:+:))
 
 import Control.Applicative (liftA2, (<$>))
 import Data.Semigroup ((<>))
@@ -240,6 +244,166 @@
 complementOrthogonal = isUnitary (selectReal 1e-3 1e-7) . Ortho.complement
 
 
+affineSpanFromKernel ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   (Matrix.Wide ShapeInt ShapeInt a, Vector ShapeInt a) -> Bool
+affineSpanFromKernel (a, by) =
+   let b = Vector.take (Shape.size $ Matrix.height a) by
+       y = Vector.drop (Shape.size $ Matrix.height a) by
+       (c,d) = Ortho.affineSpanFromKernel a b
+   in approxVectorTol
+         (selectReal 1e-3 1e-7)
+         b
+         (a#*|(c#*|y|+|d))
+
+affineKernelFromSpan ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   (Matrix.Tall ShapeInt ShapeInt a, Vector ShapeInt a, Vector ShapeInt a) ->
+   Bool
+affineKernelFromSpan (c,y,d) =
+   let (a,b) = Ortho.affineKernelFromSpan c d
+   in approxVectorTol
+         (selectReal 1e-3 1e-7)
+         b
+         (a#*|(c#*|y|+|d))
+
+
+projectHit ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   (Matrix.Wide ShapeInt ShapeInt a,
+    Vector ShapeInt a,
+    Vector ShapeInt a) ->
+   Bool
+projectHit (b,x,d) =
+   approxVectorTol
+      (selectReal 1e-3 1e-9)
+      d
+      (b #*| Ortho.project b d x)
+
+
+leastSquaresNoConstraint ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   (Matrix.Tall ShapeInt ShapeInt a, Vector ShapeInt a) -> Bool
+leastSquaresNoConstraint (a, b) =
+   approxVectorTol
+      (selectReal 0.1 1e-7)
+      (ArrMatrix.unliftColumn MatrixShape.ColumnMajor (Ortho.leastSquares a) b)
+      (Ortho.leastSquaresConstraint
+         (Matrix.fromFull a) b
+         (Matrix.zero $
+          MatrixShape.wide MatrixShape.ColumnMajor Shape.Zero (Matrix.width a))
+         (Vector.zero Shape.Zero))
+
+leastSquaresConstraintUnique ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   (Vector ShapeInt a,
+    Matrix.General ShapeInt ShapeInt a, Matrix.Square ShapeInt a,
+    Vector ShapeInt a) ->
+   Bool
+leastSquaresConstraintUnique (c,a,b,d) =
+   approxVectorTol
+      (selectReal 0.1 1e-7)
+      d
+      (b #*| Ortho.leastSquaresConstraint a c (Matrix.generalizeTall b) d)
+
+splitLSCStack ::
+   (Shape.C height, Shape.C constraints, Shape.C width, Class.Floating a) =>
+   Matrix.Tall (constraints:+:height) width a ->
+   Vector (constraints:+:height) a ->
+   ((Matrix.General height width a, Vector height a),
+    (Matrix.Wide constraints width a, Vector constraints a))
+splitLSCStack baTall dc =
+   let ba = Matrix.fromFull baTall
+       b = Matrix.wideFromGeneral $ Matrix.takeTop ba
+       a = Matrix.takeBottom ba
+       (d,c) = (Vector.takeLeft dc, Vector.takeRight dc)
+   in ((a,c),(b,d))
+
+leastSquaresConstraintAdmissible ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   (Matrix.Tall (ShapeInt:+:ShapeInt) ShapeInt a,
+    Vector (ShapeInt:+:ShapeInt) a) ->
+   Bool
+leastSquaresConstraintAdmissible (ba,dc) =
+   let ((a,c),(b,d)) = splitLSCStack ba dc
+   in approxVectorTol
+         (selectReal 0.1 1e-7)
+         d
+         (b #*| Ortho.leastSquaresConstraint a c b d)
+
+leastSquaresConstraintMinimal ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   (Matrix.Tall (ShapeInt:+:ShapeInt) ShapeInt a,
+    Vector ShapeInt a,
+    Vector (ShapeInt:+:ShapeInt) a) ->
+   Bool
+leastSquaresConstraintMinimal (ba,x,dc) =
+   let ((a,c),(b,d)) = splitLSCStack ba dc
+   in Vector.norm2 (c |-| a #*| Ortho.leastSquaresConstraint a c b d)
+      <=
+      Vector.norm2 (c |-| a #*| Ortho.project b d x) + selectReal 1e-1 1e-10
+
+
+gaussMarkovLinearModelMinimumNorm ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar, Eq a) =>
+   (Matrix.Wide ShapeInt ShapeInt a, Vector ShapeInt a) -> Bool
+gaussMarkovLinearModelMinimumNorm (b, d) =
+   let (x,y) =
+         Ortho.gaussMarkovLinearModel
+            (Matrix.zero $
+             MatrixShape.tall MatrixShape.ColumnMajor
+               (Matrix.height b) Shape.Zero)
+            (Matrix.fromFull b) d
+   in x == Vector.zero Shape.Zero
+      &&
+      approxVectorTol
+         (selectReal 0.1 1e-7)
+         y
+         (ArrMatrix.unliftColumn MatrixShape.ColumnMajor
+            (Ortho.minimumNorm b) d)
+
+gaussMarkovLinearModelUnique ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar, Eq a) =>
+   (Matrix.Square ShapeInt a, Vector ShapeInt a) -> Bool
+gaussMarkovLinearModelUnique (a,d) =
+   let (x,y) =
+         Ortho.gaussMarkovLinearModel
+            (Matrix.generalizeWide a)
+            (Matrix.zero $
+             MatrixShape.general MatrixShape.ColumnMajor
+               (Matrix.height a) Shape.Zero)
+            d
+   in y == Vector.zero Shape.Zero
+      &&
+      approxVectorTol (selectReal 0.1 1e-7) d (a #*| x)
+
+gaussMarkovLinearModelAdmissible ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   (Matrix.Wide ShapeInt (ShapeInt:+:ShapeInt) a, Vector ShapeInt a) -> Bool
+gaussMarkovLinearModelAdmissible (abWide,d) =
+   let ab = Matrix.fromFull abWide
+       a = Matrix.tallFromGeneral $ Matrix.takeLeft ab
+       b = Matrix.takeRight ab
+   in approxVectorTol
+         (selectReal 0.1 1e-7)
+         d
+         (ab #*| uncurry Vector.append (Ortho.gaussMarkovLinearModel a b d))
+
+gaussMarkovLinearModelMinimal ::
+   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   (Matrix.Wide ShapeInt (ShapeInt:+:ShapeInt) a,
+    Vector ShapeInt a,
+    Vector (ShapeInt:+:ShapeInt) a) -> Bool
+gaussMarkovLinearModelMinimal (abWide,d,xy) =
+   let ab = Matrix.fromFull abWide
+       a = Matrix.tallFromGeneral $ Matrix.takeLeft ab
+       b = Matrix.takeRight ab
+   in Vector.norm2 (snd $ Ortho.gaussMarkovLinearModel a b d)
+      <=
+      Vector.norm2 (Vector.takeRight $ Ortho.project abWide d xy)
+         + selectReal 1e-3 1e-10
+
+
 householderReconstruction ::
    (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
    Matrix.General ShapeInt ShapeInt a -> Bool
@@ -392,6 +556,53 @@
       checkForAll Gen.tall complementBiorthogonal) :
    ("complementOrthogonal",
       checkForAll Gen.tall complementOrthogonal) :
+
+   ("affineSpanFromKernel",
+      checkForAll
+         ((,) <$> Gen.fullRankWide <#*|> Gen.vector) affineSpanFromKernel) :
+   ("affineKernelFromSpan",
+      checkForAll
+         ((,,) <$> Gen.tall <#*|> Gen.vector <|=|> Gen.vector)
+         affineKernelFromSpan) :
+
+   ("projectHit",
+      checkForAll
+         ((,,) <$> Gen.fullRankWide <#*|> Gen.vector <|=|> Gen.vector)
+         projectHit) :
+   ("leastSquaresNoConstraint",
+      checkForAll
+         ((,) <$> Gen.transpose Gen.fullRankTall <#*|> Gen.vector)
+         leastSquaresNoConstraint) :
+   ("leastSquaresConstraintUnique",
+      checkForAll
+         ((,,,) <$>
+            Gen.vector <-*#> Gen.matrix <-*#> Gen.invertible <|=|> Gen.vector)
+         leastSquaresConstraintUnique) :
+   ("leastSquaresConstraintAdmissible",
+      checkForAll
+         ((,) <$> Gen.transpose Gen.lscStack <#*|> Gen.vector)
+         leastSquaresConstraintAdmissible) :
+   ("leastSquaresConstraintMinimal",
+      checkForAll
+         ((,,) <$> Gen.lscStack <#*|> Gen.vector <|=|> Gen.vector)
+         leastSquaresConstraintMinimal) :
+   ("gaussMarkovLinearModelMinimumNorm",
+      checkForAll
+         ((,) <$> Gen.transpose Gen.fullRankWide <#*|> Gen.vector)
+         gaussMarkovLinearModelMinimumNorm) :
+   ("gaussMarkovLinearModelUnique",
+      checkForAll
+         ((,) <$> Gen.invertible <#*|> Gen.vector)
+         gaussMarkovLinearModelUnique) :
+   ("gaussMarkovLinearModelAdmissible",
+      checkForAll
+         ((,) <$> fmap Matrix.transpose Gen.lscStack <#*|> Gen.vector)
+         gaussMarkovLinearModelAdmissible) :
+   ("gaussMarkovLinearModelMinimal",
+      checkForAll
+         ((,,) <$> fmap Matrix.transpose Gen.lscStack
+               <#*|> Gen.vector <|=|> Gen.vector)
+         gaussMarkovLinearModelMinimal) :
 
    ("triangularLeastSquares",
       checkForAll genFullRankTallRHS triangularLeastSquares) :
diff --git a/test/Test/Permutation.hs b/test/Test/Permutation.hs
--- a/test/Test/Permutation.hs
+++ b/test/Test/Permutation.hs
@@ -20,13 +20,10 @@
 
 import qualified Numeric.Netlib.Class as Class
 
-import qualified Data.Array.Comfort.Storable as Array
 import qualified Data.Array.Comfort.Shape as Shape
 
 import Control.Monad (forM)
-import Control.Applicative (liftA2, (<$>))
-
-import Data.Semigroup ((<>))
+import Control.Applicative ((<$>))
 
 import qualified Test.QuickCheck as QC
 
@@ -44,27 +41,6 @@
 genPerm = fmap (Perm.fromPivots NonInverted) . genPivots
 
 
-permutationPivots :: Inversion -> Pivots -> Bool
-permutationPivots inv xs =
-   Array.toList (Perm.toPivots inv (Perm.fromPivots inv xs))
-   ==
-   Array.toList xs
-
-determinantMultiply :: (Permutation ShapeInt, Permutation ShapeInt) -> Bool
-determinantMultiply (p0,p1) =
-   Perm.determinant (Perm.multiply p0 p1)
-   ==
-   Perm.determinant p0 <> Perm.determinant p1
-
-transposeMultiply :: (Permutation ShapeInt, Permutation ShapeInt) -> Bool
-transposeMultiply (p0,p1) =
-   (Array.toList $ Perm.toPivots NonInverted $
-    Perm.transpose (Perm.multiply p0 p1))
-   ==
-   (Array.toList $ Perm.toPivots NonInverted $
-    Perm.multiply (Perm.transpose p1) (Perm.transpose p0))
-
-
 genPermutation :: (Dim sh) => Gen.Matrix sh sh a (Permutation sh)
 genPermutation =
    flip Gen.mapGen Gen.squareDim $ \_maxElem sh ->
@@ -117,28 +93,6 @@
    (Class.Floating a, Eq a) => Matrix (Permutation ShapeInt) a -> Bool
 determinantNumber p =
    PermMatrix.determinant p == Square.determinant (PermMatrix.toMatrix p)
-
-
-tests :: [(String, QC.Property)]
-tests =
-   ("permutationPivots",
-      QC.property $
-         QC.forAll QC.arbitraryBoundedEnum $ \inv ->
-         QC.forAll (QC.arbitrary >>= genPivots) $ \pivot ->
-            permutationPivots inv pivot) :
-   ("determinantMultiply",
-      QC.property $
-         QC.forAll
-            (do nat <- QC.arbitrary
-                liftA2 (,) (genPerm nat) (genPerm nat))
-            determinantMultiply) :
-   ("transposeMultiply",
-      QC.property $
-         QC.forAll
-            (do nat <- QC.arbitrary
-                liftA2 (,) (genPerm nat) (genPerm nat))
-            transposeMultiply) :
-   []
 
 
 
diff --git a/test/Test/Shape.hs b/test/Test/Shape.hs
--- a/test/Test/Shape.hs
+++ b/test/Test/Shape.hs
@@ -99,19 +99,19 @@
    Banded (MatrixShape.Banded sub super vert horiz height width)
 
 instance
-   (Extent.C horiz, Extent.C vert,
+   (Extent.C vert, Extent.C horiz,
     Show height, Show width, Shape.C height, Shape.C width) =>
       Show (Banded vert horiz height width) where
    showsPrec p (Banded sh) = showsPrec p sh
 
 instance
-   (Extent.C horiz, Extent.C vert, Shape.C height, Shape.C width) =>
+   (Extent.C vert, Extent.C horiz, Shape.C height, Shape.C width) =>
       Shape.C (Banded vert horiz height width) where
    size (Banded sh) = Shape.size sh
    uncheckedSize (Banded sh) = Shape.uncheckedSize sh
 
 instance
-   (Extent.C horiz, Extent.C vert,
+   (Extent.C vert, Extent.C horiz,
     Shape.Indexed height, Shape.Indexed width) =>
       Shape.Indexed (Banded vert horiz height width) where
    type Index (Banded vert horiz height width) =
@@ -125,7 +125,7 @@
    uncheckedSizeOffset (Banded sh) = Shape.uncheckedSizeOffset sh
 
 instance
-   (Extent.C horiz, Extent.C vert,
+   (Extent.C vert, Extent.C horiz,
     Shape.InvIndexed height, Shape.InvIndexed width) =>
       Shape.InvIndexed (Banded vert horiz height width) where