{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE EmptyDataDecls #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FlexibleContexts #-}
{-# 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)