{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE GADTs #-}
module Numeric.LAPACK.Matrix.Array.Multiply where

import qualified Numeric.LAPACK.Matrix.Full as Full
import qualified Numeric.LAPACK.Matrix.BandedHermitian as BandedHermitian
import qualified Numeric.LAPACK.Matrix.Banded as Banded
import qualified Numeric.LAPACK.Matrix.Hermitian as Hermitian
import qualified Numeric.LAPACK.Matrix.Symmetric as Symmetric
import qualified Numeric.LAPACK.Matrix.Triangular as Triangular
import qualified Numeric.LAPACK.Matrix.Mosaic.Basic as Mosaic
import qualified Numeric.LAPACK.Matrix.Banded.Basic as BandedBasic
import qualified Numeric.LAPACK.Matrix.Triangular.Basic as TriangularBasic
import qualified Numeric.LAPACK.Matrix.Square.Basic as SquareBasic
import qualified Numeric.LAPACK.Matrix.Basic as FullBasic
import qualified Numeric.LAPACK.Matrix.Type as Matrix
import qualified Numeric.LAPACK.Matrix.Array.Basic as OmniMatrix
import qualified Numeric.LAPACK.Matrix.Array.Unpacked as Unpacked
import qualified Numeric.LAPACK.Matrix.Array.Private as ArrMatrix
import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
import qualified Numeric.LAPACK.Matrix.Shape.Omni as Omni
import qualified Numeric.LAPACK.Matrix.Extent.Strict as ExtentStrict
import qualified Numeric.LAPACK.Matrix.Extent.Private as ExtentPriv
import qualified Numeric.LAPACK.Matrix.Extent as Extent
import qualified Numeric.LAPACK.Vector.Private as VectorPriv
import qualified Numeric.LAPACK.Vector as Vector
import Numeric.LAPACK.Matrix.Layout.Private
         (Filled, Empty, Bands, Packed, Unpacked, UnaryProxy, natFromProxy)
import Numeric.LAPACK.Matrix.Array.Banded (Banded)
import Numeric.LAPACK.Matrix.Array.Private (ArrayMatrix, diagTag)
import Numeric.LAPACK.Matrix.Type (Matrix)
import Numeric.LAPACK.Matrix.Modifier (Transposition(NonTransposed,Transposed))
import Numeric.LAPACK.Matrix.Shape.Omni (StripSingleton(StripBands, StripFilled))
import Numeric.LAPACK.Vector (Vector)
import Numeric.LAPACK.Shape.Private (Unchecked(Unchecked, deconsUnchecked))

import qualified Numeric.Netlib.Class as Class

import qualified Data.Array.Comfort.Storable as Array
import qualified Data.Array.Comfort.Shape as Shape

import qualified Type.Data.Num.Unary.Proof as Proof
import qualified Type.Data.Num.Unary as Unary
import qualified Type.Data.Bool as TBool
import Type.Data.Num.Unary (Succ, (:+:))
import Type.Base.Proxy (Proxy(Proxy))

import qualified Data.Stream as Stream
import Data.Stream (Stream)
import Data.Tuple.HT (double)
import Data.Function.HT (Id)



type FactorQuadratic pack property lower upper sh =
      Factor pack property lower upper
         Extent.Shape Extent.Small Extent.Small sh sh

data Factor pack property lower upper meas vert horiz height width where
   FactorUpperTriangular ::
      (Layout.Packing pack, Omni.TriDiag diag) =>
      FactorQuadratic pack diag Empty Filled sh
   FactorLowerTriangular ::
      (Layout.Packing pack, Omni.TriDiag diag) =>
      FactorQuadratic pack diag Filled Empty sh
   FactorSymmetric ::
      (Layout.Packing pack) =>
      FactorQuadratic pack Omni.Symmetric Filled Filled sh
   FactorHermitian ::
      (Layout.Packing pack, TBool.C neg, TBool.C zero, TBool.C pos) =>
      FactorQuadratic pack (Omni.Hermitian neg zero pos) Filled Filled sh
   FactorBanded ::
      (Unary.Natural sub, Unary.Natural super) =>
      Factor Packed Omni.Arbitrary (Bands sub) (Bands super)
         meas vert horiz height width
   FactorBandedHermitian ::
      (TBool.C neg, TBool.C zero, TBool.C pos, Unary.Natural offDiag) =>
      FactorQuadratic Packed (Omni.Hermitian neg zero pos)
         (Bands offDiag) (Bands offDiag) sh
   FactorUnitBandedTriangular ::
      (Omni.BandedTriangular sub super,
       Omni.BandedTriangular super sub) =>
      FactorQuadratic Packed Omni.Unit (Bands sub) (Bands super) sh
   FactorFull ::
      Factor Unpacked property lower upper meas vert horiz height width

factorSingleton ::
   (Layout.Packing pack, Omni.Property property,
    Omni.Strip lower, Omni.Strip upper,
    Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
   Omni.Omni pack property lower upper meas vert horiz height width ->
   Factor pack property lower upper meas vert horiz height width
factorSingleton shape =
   case Omni.packTag shape of
      Layout.Packed ->
         case shape of
            Omni.UpperTriangular _ -> FactorUpperTriangular
            Omni.LowerTriangular _ -> FactorLowerTriangular
            Omni.Symmetric _ -> FactorSymmetric
            Omni.Hermitian _ -> FactorHermitian
            Omni.Banded _ -> FactorBanded
            Omni.BandedHermitian _ -> FactorBandedHermitian
            Omni.UnitBandedTriangular _ -> FactorUnitBandedTriangular
      Layout.Unpacked ->
         case (Omni.extent shape, Omni.property shape, Omni.strips shape) of
            (ExtentPriv.Square _, Omni.PropArbitrary,
               (StripBands Unary.Zero, StripFilled)) ->
                  FactorUpperTriangular
            (ExtentPriv.Square _, Omni.PropArbitrary,
               (StripFilled, StripBands Unary.Zero)) ->
                  FactorLowerTriangular
            (ExtentPriv.Square _, Omni.PropUnit,
               (StripBands Unary.Zero, StripFilled)) ->
                  FactorUpperTriangular
            (ExtentPriv.Square _, Omni.PropUnit,
               (StripFilled, StripBands Unary.Zero)) ->
                  FactorLowerTriangular
            (ExtentPriv.Square _, Omni.PropSymmetric,
               (StripFilled, StripFilled)) ->
                  FactorSymmetric
            (ExtentPriv.Square _, Omni.PropHermitian,
               (StripFilled, StripFilled)) ->
                  FactorHermitian
            _ -> FactorFull

matrixVector ::
   (Layout.Packing pack) =>
   (Omni.Property property, Omni.Strip lower, Omni.Strip upper) =>
   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
   ArrayMatrix pack property lower upper meas vert horiz height width a ->
   Vector width a -> Vector height a
matrixVector a x =
   case factorSingleton $ ArrMatrix.shape a of
      FactorFull -> Full.multiplyVector (OmniMatrix.toFull a) x
      FactorUpperTriangular -> Triangular.multiplyVector a x
      FactorLowerTriangular -> Triangular.multiplyVector a x
      FactorSymmetric -> Symmetric.multiplyVector a x
      FactorHermitian -> Hermitian.multiplyVector NonTransposed a x
      FactorBanded -> Banded.multiplyVector a x
      FactorUnitBandedTriangular ->
         BandedBasic.multiplyVector (ArrMatrix.toVector a) x
      FactorBandedHermitian ->
         BandedHermitian.multiplyVector NonTransposed a x

vectorMatrix ::
   (Layout.Packing pack) =>
   (Omni.Property property, Omni.Strip lower, Omni.Strip upper) =>
   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
   (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
   Vector height a ->
   ArrayMatrix pack property lower upper meas vert horiz height width a ->
   Vector width a
vectorMatrix x a =
   case factorSingleton $ ArrMatrix.shape a of
      FactorFull ->
         Full.multiplyVector (OmniMatrix.toFull $ Unpacked.transpose a) x
      FactorUpperTriangular ->
         Triangular.multiplyVector (Triangular.transpose a) x
      FactorLowerTriangular ->
         Triangular.multiplyVector (Triangular.transpose a) x
      FactorSymmetric -> Symmetric.multiplyVector a x
      FactorHermitian -> Hermitian.multiplyVector Transposed a x
      FactorBanded -> Banded.multiplyVector (Banded.transpose a) x
      FactorUnitBandedTriangular ->
         BandedBasic.multiplyVector
            (BandedBasic.transpose $ ArrMatrix.toVector a) x
      FactorBandedHermitian ->
         BandedHermitian.multiplyVector Transposed a x


{- ToDo:
We can remove the transposition parameter
if we can transpose Hermitian in constant time.
Then we can express all variations using 'transpose'.
-}
transposableSquare ::
   (Layout.Packing pack) =>
   (Omni.Property property, Omni.Strip lower, Omni.Strip upper) =>
   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
   (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
   Transposition ->
   ArrMatrix.Quadratic pack property lower upper height a ->
   ArrMatrix.Full meas vert horiz height width a ->
   ArrMatrix.Full meas vert horiz height width a
transposableSquare trans a =
   ($a) $
   case ArrMatrix.shape a of
      Omni.Hermitian _ -> Hermitian.multiplyFull trans
      Omni.BandedHermitian _ -> BandedHermitian.multiplyFull trans
      _ ->
         case trans of
            Transposed -> squareFull . Matrix.transpose
            NonTransposed -> squareFull

extentFromSquare ::
   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
   Extent.Map Extent.Shape Extent.Small Extent.Small meas vert horiz size size
extentFromSquare = ExtentStrict.Map ExtentPriv.fromSquare

fullSquare ::
   (Layout.Packing pack) =>
   (Omni.Property property, Omni.Strip lower, Omni.Strip upper) =>
   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
   ArrMatrix.Full meas vert horiz height width a ->
   ArrMatrix.Quadratic pack property lower upper width a ->
   ArrMatrix.Full meas vert horiz height width a
fullSquare = Matrix.swapMultiply $ transposableSquare Transposed

squareFull ::
   (Layout.Packing pack) =>
   (Omni.Property property, Omni.Strip lower, Omni.Strip upper) =>
   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
   (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
   ArrMatrix.Quadratic pack property lower upper height a ->
   ArrMatrix.Full meas vert horiz height width a ->
   ArrMatrix.Full meas vert horiz height width a
squareFull a =
   case factorSingleton $ ArrMatrix.shape a of
      FactorFull -> Unpacked.multiply $ Unpacked.mapExtent extentFromSquare a
      FactorUpperTriangular -> Triangular.multiplyFull a
      FactorLowerTriangular -> Triangular.multiplyFull a
      FactorSymmetric -> Symmetric.multiplyFull a
      FactorHermitian -> Hermitian.multiplyFull NonTransposed a
      FactorBanded -> Banded.multiplyFull $ Banded.mapExtent extentFromSquare a
      FactorBandedHermitian -> BandedHermitian.multiplyFull NonTransposed a
      FactorUnitBandedTriangular ->
         ArrMatrix.lift2
            (BandedBasic.multiplyFull .
             BandedBasic.mapExtent ExtentPriv.fromSquare) a


vectorSquare :: (Shape.C sh, Class.Floating a) => Vector sh a -> Vector sh a
vectorSquare =
   VectorPriv.recheck . uncurry Vector.mul . double . VectorPriv.uncheck

square ::
   (Layout.Packing pack, Omni.Property property) =>
   (MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper) =>
   (Shape.C sh, Class.Floating a) =>
   ArrMatrix.Quadratic pack property lower upper sh a ->
   ArrMatrix.Quadratic pack property lower upper sh a
square a =
   case Omni.powerSingleton $ ArrMatrix.shape a of
      Omni.PowerIdentity -> a
      Omni.PowerUpperTriangular -> Triangular.square a
      Omni.PowerLowerTriangular -> Triangular.square a
      Omni.PowerSymmetric -> Symmetric.square a
      Omni.PowerHermitian -> Hermitian.square a
      Omni.PowerFull -> ArrMatrix.liftUnpacked1 SquareBasic.square a
      Omni.PowerDiagonal ->
         case ArrMatrix.shape a of
            Omni.Banded _ -> ArrMatrix.lift1 vectorSquare a
            -- ToDo: could be optimized by processing only the real part
            Omni.BandedHermitian _ -> ArrMatrix.lift1 vectorSquare a
            Omni.UnitBandedTriangular _ -> a
            Omni.Full _ ->
               ArrMatrix.liftUnpacked1
                  (SquareBasic.diagonalOrder (ArrMatrix.order a) .
                   vectorSquare . SquareBasic.takeDiagonal)
                  a

power ::
   (Layout.Packing pack, Omni.Property property) =>
   (MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper) =>
   (Shape.C sh, Class.Floating a) =>
   Integer ->
   ArrMatrix.Quadratic pack property lower upper sh a ->
   ArrMatrix.Quadratic pack property lower upper sh a
power n a =
   case Omni.powerSingleton $ ArrMatrix.shape a of
      Omni.PowerIdentity -> a
      Omni.PowerUpperTriangular ->
         ArrMatrix.lift1 (Mosaic.power (diagTag a) n) a
      Omni.PowerLowerTriangular ->
         ArrMatrix.lift1 (Mosaic.power (diagTag a) n) a
      Omni.PowerSymmetric ->
         ArrMatrix.lift1 (Mosaic.power Omni.Arbitrary n) a
      Omni.PowerHermitian ->
         ArrMatrix.lift1 (Mosaic.power Omni.Arbitrary n) a
      Omni.PowerFull ->
         ArrMatrix.liftUnpacked1 (SquareBasic.power n) a
      Omni.PowerDiagonal ->
         case ArrMatrix.shape a of
            Omni.Banded _ -> ArrMatrix.lift1 (Array.map (^n)) a
            Omni.BandedHermitian _ -> ArrMatrix.lift1 (Array.map (^n)) a
            Omni.UnitBandedTriangular _ -> a
            Omni.Full _ ->
               ArrMatrix.liftUnpacked1
                  (SquareBasic.diagonalOrder (ArrMatrix.order a) .
                   Array.map (^n) . SquareBasic.takeDiagonal)
                  a

powers, powers1 ::
   (Layout.Packing pack, Omni.Property property) =>
   (MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper) =>
   (Shape.C sh, Class.Floating a) =>
   ArrMatrix.Quadratic pack property lower upper sh a ->
   Stream (ArrMatrix.Quadratic pack property lower upper sh a)
powers a = Stream.Cons (OmniMatrix.identityFrom a) (powers1 a)
powers1 a =
   case Omni.powerSingleton $ ArrMatrix.shape a of
      Omni.PowerIdentity -> Stream.repeat a
      Omni.PowerUpperTriangular -> powers1Gen a Triangular.multiply
      Omni.PowerLowerTriangular -> powers1Gen a Triangular.multiply
      Omni.PowerSymmetric ->
         fmap ArrMatrix.lift0 $
         Mosaic.powers1 Omni.Arbitrary $ ArrMatrix.toVector a
      Omni.PowerHermitian ->
         fmap ArrMatrix.lift0 $
         Mosaic.powers1 Omni.Arbitrary $ ArrMatrix.toVector a
      Omni.PowerFull ->
         powers1Gen a $ ArrMatrix.liftUnpacked2 FullBasic.multiply
      Omni.PowerDiagonal ->
         case ArrMatrix.shape a of
            Omni.Banded _ ->
               fmap ArrMatrix.Array $ powers1Diagonal $ ArrMatrix.unwrap a
            Omni.BandedHermitian _ ->
               fmap ArrMatrix.Array $ powers1Diagonal $ ArrMatrix.unwrap a
            Omni.UnitBandedTriangular _ -> Stream.repeat a
            Omni.Full _ ->
               fmap (ArrMatrix.liftUnpacked0 .
                     SquareBasic.diagonalOrder (ArrMatrix.order a)) $
               powers1Diagonal $
               SquareBasic.takeDiagonal $ ArrMatrix.unpackedToVector a

powers1Diagonal ::
   (Shape.C sh, Class.Floating a) => Vector sh a -> Stream (Vector sh a)
powers1Diagonal a =
   let b = VectorPriv.uncheck a
   in fmap VectorPriv.recheck $ Stream.iterate (flip (Array.zipWith (*)) b) b

powers1Gen ::
   (Layout.Packing pack, Omni.Property property) =>
   (MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper) =>
   (Shape.C sh, Class.Floating a) =>
   ArrMatrix.Quadratic pack property lower upper sh a ->
   (ArrMatrix.Quadratic pack property lower upper (Unchecked sh) a ->
    ArrMatrix.Quadratic pack property lower upper (Unchecked sh) a ->
    ArrMatrix.Quadratic pack property lower upper (Unchecked sh) a) ->
   Stream (ArrMatrix.Quadratic pack property lower upper sh a)
powers1Gen a mul =
   let b = OmniMatrix.mapSquareSize Unchecked a
   in fmap (OmniMatrix.mapSquareSize deconsUnchecked) $
      Stream.iterate (flip mul b) b


type family MultipliedPacking packA packB
type instance MultipliedPacking Packed packB = packB
type instance MultipliedPacking Unpacked packB = Unpacked

type family PackingByStrip lower upper meas pack
type instance PackingByStrip lower upper meas Unpacked = Unpacked
type instance PackingByStrip Filled upper Extent.Size pack = Unpacked
type instance PackingByStrip Empty upper Extent.Shape pack = pack
type instance PackingByStrip (Bands k) (Bands l) meas pack = pack
type instance PackingByStrip Filled Filled meas pack = Unpacked
type instance PackingByStrip
                     Filled (Bands (Unary.Succ l)) meas pack = Unpacked


{- |
This functions allows to multiply two matrices of arbitrary special features
and returns the most special matrix type possible (almost).
At the first glance, this is handy.
At the second glance, this has some problems.
First of all, we may refine the types in future
and then multiplication may return a different, more special type than before.
Second, if you write code with polymorphic matrix types,
then 'matrixMatrix' may leave you with constraints like
@ExtentPriv.Multiply vert vert ~ vert@.
That constraint is always fulfilled but the compiler cannot infer that.
Because of these problems
you may instead consider using specialised 'Basic.multiply' functions
from the various modules for production use.
Btw. 'MultiplyVector' and 'MultiplySquare'
are much less problematic,
because the input and output are always dense vectors or dense matrices.


We require e.g. both

> Extent.Multiply vertA vertB ~ vertC

and

> Extent.Multiply vertB vertA ~ vertC

This makes the combination commutative.
At first glance this looks counterintuitive,
but this allows the type checker to infer
that "shape A" times "square" is "shape A".
-}
matrixMatrix ::
   (Layout.Packing packA, Omni.Property propertyA) =>
   (Layout.Packing packB, Omni.Property propertyB) =>
   (Layout.Packing packC, Omni.Property propertyC) =>
   (Omni.Strip lowerA, Omni.Strip upperA) =>
   (Omni.Strip lowerB, Omni.Strip upperB) =>
   (Omni.Strip lowerC, Omni.Strip upperC) =>
   (MultipliedPacking packA packB ~ pack) =>
   (MultipliedPacking packB packA ~ pack) =>
   (Omni.MultipliedStrip lowerA lowerB ~ lowerC) =>
   (Omni.MultipliedStrip lowerB lowerA ~ lowerC) =>
   (Omni.MultipliedStrip upperA upperB ~ upperC) =>
   (Omni.MultipliedStrip upperB upperA ~ upperC) =>
   (Omni.MultipliedBands lowerA lowerB ~ lowerC) =>
   (Omni.MultipliedBands lowerB lowerA ~ lowerC) =>
   (Omni.MultipliedBands upperA upperB ~ upperC) =>
   (Omni.MultipliedBands upperB upperA ~ upperC) =>
   (Omni.MultipliedProperty propertyA propertyB ~ propertyAB) =>
   (Omni.MultipliedProperty propertyB propertyA ~ propertyAB) =>
   (Omni.UnitIfTriangular lowerC upperC ~ diag) =>
   (Omni.UnitIfTriangular upperC lowerC ~ diag) =>
   (Omni.MergeUnit propertyAB diag ~ propertyC) =>
   (Omni.MergeUnit diag propertyAB ~ propertyC) =>
   (PackingByStrip lowerC upperC measC pack ~ packC) =>
   (PackingByStrip upperC lowerC measC pack ~ packC) =>
   (Extent.Measure measA, Extent.C vertA, Extent.C horizA) =>
   (Extent.Measure measB, Extent.C vertB, Extent.C horizB) =>
   (ExtentPriv.MultiplyMeasure measA measB ~ measC) =>
   (ExtentPriv.MultiplyMeasure measB measA ~ measC) =>
   (ExtentPriv.Multiply vertA  vertB  ~ vertC)  =>
   (ExtentPriv.Multiply vertB  vertA  ~ vertC)  =>
   (ExtentPriv.Multiply horizA horizB ~ horizC) =>
   (ExtentPriv.Multiply horizB horizA ~ horizC) =>
   (Shape.C height, Shape.C fuse, Eq fuse, Shape.C width, Class.Floating a) =>
   ArrayMatrix packA propertyA lowerA upperA measA vertA horizA height fuse a ->
   ArrayMatrix packB propertyB lowerB upperB measB vertB horizB fuse width a ->
   ArrayMatrix packC propertyC lowerC upperC measC vertC horizC height width a
matrixMatrix a b =
   case (ArrMatrix.shape a, ArrMatrix.shape b) of
      (shapeA@(Omni.Full _), shapeB@(Omni.Full _)) ->
         case factorSingleton shapeA of
            FactorUpperTriangular -> squareUnpacked a b
            FactorLowerTriangular -> squareUnpacked a b
            FactorSymmetric -> squareUnpacked a b
            FactorHermitian -> squareUnpacked a b
            _ ->
               case factorSingleton shapeB of
                  FactorUpperTriangular -> unpackedSquare a b
                  FactorLowerTriangular -> unpackedSquare a b
                  FactorSymmetric -> unpackedSquare a b
                  FactorHermitian -> unpackedSquare a b
                  _ -> unpackedUnpacked a b
      (Omni.Full _, Omni.UpperTriangular _) -> unpackedSquare a b
      (Omni.Full _, Omni.LowerTriangular _) -> unpackedSquare a b
      (Omni.Full _, Omni.Symmetric _) -> unpackedSquare a b
      (Omni.Full _, Omni.Hermitian _) -> unpackedSquare a b
      (Omni.Full _, Omni.Banded _) -> unpackedBanded a b
      (Omni.Full _, Omni.UnitBandedTriangular _) -> unpackedSquare a b
      (Omni.Full _, Omni.BandedHermitian _) -> unpackedSquare a b
      (Omni.UpperTriangular _, Omni.Full _) -> squareUnpacked a b
      (Omni.UpperTriangular _, Omni.UpperTriangular _) ->
         case (diagTag a, diagTag b) of
            (diagA@Omni.Unit, Omni.Unit) ->
               ArrMatrix.lift2 (TriangularBasic.multiply diagA) a b
            (diagA@Omni.Unit, Omni.Arbitrary) ->
               ArrMatrix.lift2 (TriangularBasic.multiply diagA) a b
            (diagA@Omni.Arbitrary, _) ->
               ArrMatrix.lift2 (TriangularBasic.multiply diagA) a b
      (Omni.UpperTriangular _, Omni.LowerTriangular _) ->
         squareFull a (OmniMatrix.toFull b)
      (Omni.UpperTriangular _, Omni.Symmetric _) ->
         squareFull a (Symmetric.toSquare b)
      (Omni.UpperTriangular _, Omni.Hermitian _) ->
         squareFull a (Hermitian.toSquare b)
      (Omni.UpperTriangular _, Omni.Banded _) ->
         case ArrMatrix.subBandsSingleton b of
            Unary.Succ -> unpackedBanded (OmniMatrix.toFull a) b
            Unary.Zero ->
               case Matrix.extent b of
                  ExtentPriv.Square _ ->
                     Triangular.takeUpper $ asArbitrary $
                     unpackedBanded (OmniMatrix.toFull a) b
                  ExtentPriv.Separate _ _ ->
                     unpackedBanded (OmniMatrix.toFull a) b
      (Omni.UpperTriangular _, Omni.BandedHermitian _) ->
         case ArrMatrix.subBandsSingleton b of
            Unary.Succ -> unpackedSquare (OmniMatrix.toFull a) b
            Unary.Zero ->
               ArrMatrix.lift1 TriangularBasic.takeUpper $
               fullSquare (OmniMatrix.toFull a) b
      (Omni.UpperTriangular _, Omni.UnitBandedTriangular _) ->
         case ArrMatrix.subBandsSingleton b of
            Unary.Succ -> unpackedSquare (OmniMatrix.toFull a) b
            Unary.Zero ->
               case diagTag a of
                  Omni.Unit ->
                     ArrMatrix.lift1 id $
                     Triangular.takeUpper $ fullSquare (OmniMatrix.toFull a) b
                  Omni.Arbitrary ->
                     Triangular.takeUpper $ fullSquare (OmniMatrix.toFull a) b
      (Omni.LowerTriangular _, Omni.Full _) -> squareUnpacked a b
      (Omni.LowerTriangular _, Omni.LowerTriangular _) ->
         case (diagTag a, diagTag b) of
            (diagA@Omni.Unit, Omni.Unit) ->
               ArrMatrix.lift2 (TriangularBasic.multiply diagA) a b
            (diagA@Omni.Unit, Omni.Arbitrary) ->
               ArrMatrix.lift2 (TriangularBasic.multiply diagA) a b
            (diagA@Omni.Arbitrary, _) ->
               ArrMatrix.lift2 (TriangularBasic.multiply diagA) a b
      (Omni.LowerTriangular _, Omni.UpperTriangular _) ->
         squareFull a (OmniMatrix.toFull b)
      (Omni.LowerTriangular _, Omni.Symmetric _) ->
         squareFull a (Symmetric.toSquare b)
      (Omni.LowerTriangular _, Omni.Hermitian _) ->
         squareFull a (Hermitian.toSquare b)
      (Omni.LowerTriangular _, Omni.Banded _) ->
         case ArrMatrix.superBandsSingleton b of
            Unary.Succ -> unpackedBanded (OmniMatrix.toFull a) b
            Unary.Zero ->
               case Matrix.extent b of
                  ExtentPriv.Square _ ->
                     Triangular.takeLower $ asArbitrary $
                     unpackedBanded (OmniMatrix.toFull a) b
                  ExtentPriv.Separate _ _ ->
                     unpackedBanded (OmniMatrix.toFull a) b
      (Omni.LowerTriangular _, Omni.BandedHermitian _) ->
         case ArrMatrix.superBandsSingleton b of
            Unary.Succ -> flex $ fullSquare (OmniMatrix.toFull a) b
            Unary.Zero ->
               Triangular.takeLower $ fullSquare (OmniMatrix.toFull a) b
      (Omni.LowerTriangular _, Omni.UnitBandedTriangular _) ->
         case ArrMatrix.superBandsSingleton b of
            Unary.Succ -> flex $ fullSquare (OmniMatrix.toFull a) b
            Unary.Zero ->
               case diagTag a of
                  Omni.Unit ->
                     ArrMatrix.lift1 id $
                     Triangular.takeLower $ fullSquare (OmniMatrix.toFull a) b
                  Omni.Arbitrary ->
                     Triangular.takeLower $ fullSquare (OmniMatrix.toFull a) b
      (Omni.Symmetric _, Omni.Full _) -> squareFull a (unflex b)
      (Omni.Symmetric _, Omni.LowerTriangular _) ->
         fullSquare (Symmetric.toSquare a) b
      (Omni.Symmetric _, Omni.UpperTriangular _) ->
         fullSquare (Symmetric.toSquare a) b
      (Omni.Symmetric _, Omni.Symmetric _) ->
         squareFull a (Symmetric.toSquare b)
      (Omni.Symmetric _, Omni.Hermitian _) ->
         fullSquare (Symmetric.toSquare a) b
      (Omni.Symmetric _, Omni.Banded _) ->
         unpackedBanded (Symmetric.toSquare a) b
      (Omni.Symmetric _, Omni.BandedHermitian _) ->
         fullSquare (Symmetric.toSquare a) b
      (Omni.Symmetric _, Omni.UnitBandedTriangular _) ->
         fullSquare (Symmetric.toSquare a) b
      (Omni.Hermitian _, Omni.Full _) -> squareFull a (unflex b)
      (Omni.Hermitian _, Omni.LowerTriangular _) ->
         fullSquare (Hermitian.toSquare a) b
      (Omni.Hermitian _, Omni.UpperTriangular _) ->
         fullSquare (Hermitian.toSquare a) b
      (Omni.Hermitian _, Omni.Symmetric _) ->
         squareFull a (Symmetric.toSquare b)
      (Omni.Hermitian _, Omni.Hermitian _) ->
         squareFull a (Hermitian.toSquare b)
      (Omni.Hermitian _, Omni.Banded _) ->
         unpackedBanded (Hermitian.toSquare a) b
      (Omni.Hermitian _, Omni.BandedHermitian _) ->
         fullSquare (Hermitian.toSquare a) b
      (Omni.Hermitian _, Omni.UnitBandedTriangular _) ->
         fullSquare (Hermitian.toSquare a) b
      (Omni.Banded _, Omni.Full _) -> bandedUnpacked a b
      (Omni.Banded _, Omni.Symmetric _) ->
         bandedUnpacked a (Symmetric.toSquare b)
      (Omni.Banded _, Omni.Hermitian _) ->
         bandedUnpacked a (Hermitian.toSquare b)
      (Omni.Banded _, Omni.LowerTriangular _) ->
         case ArrMatrix.superBandsSingleton a of
            Unary.Succ -> bandedUnpacked a (OmniMatrix.toFull b)
            Unary.Zero ->
               case Matrix.extent a of
                  ExtentPriv.Square _ ->
                     Triangular.takeLower $ asArbitrary $
                     bandedUnpacked a (OmniMatrix.toFull b)
                  ExtentPriv.Separate _ _ ->
                     bandedUnpacked a (OmniMatrix.toFull b)
      (Omni.Banded _, Omni.UpperTriangular _) ->
         case ArrMatrix.subBandsSingleton a of
            Unary.Succ -> bandedUnpacked a (OmniMatrix.toFull b)
            Unary.Zero ->
               case Matrix.extent a of
                  ExtentPriv.Square _ ->
                     Triangular.takeUpper $ asArbitrary $
                     bandedUnpacked a (OmniMatrix.toFull b)
                  ExtentPriv.Separate _ _ ->
                     bandedUnpacked a (OmniMatrix.toFull b)
      (Omni.UnitBandedTriangular _, Omni.Full _) -> squareUnpacked a b
      (Omni.UnitBandedTriangular _, Omni.LowerTriangular _) ->
         case ArrMatrix.superBandsSingleton a of
            Unary.Succ -> flex $ squareFull a (OmniMatrix.toFull b)
            Unary.Zero ->
               case diagTag b of
                  Omni.Unit ->
                     ArrMatrix.lift1 id $
                     Triangular.takeLower $ squareFull a (OmniMatrix.toFull b)
                  Omni.Arbitrary ->
                     Triangular.takeLower $ squareFull a (OmniMatrix.toFull b)
      (Omni.UnitBandedTriangular _, Omni.UpperTriangular _) ->
         case ArrMatrix.subBandsSingleton a of
            Unary.Succ -> flex $ squareFull a (OmniMatrix.toFull b)
            Unary.Zero ->
               case diagTag b of
                  Omni.Unit ->
                     ArrMatrix.lift1 id $
                     Triangular.takeUpper $ squareFull a (OmniMatrix.toFull b)
                  Omni.Arbitrary ->
                     Triangular.takeUpper $ squareFull a (OmniMatrix.toFull b)
      (Omni.UnitBandedTriangular _, Omni.Symmetric _) ->
         squareFull a (Symmetric.toSquare b)
      (Omni.UnitBandedTriangular _, Omni.Hermitian _) ->
         squareFull a (Hermitian.toSquare b)
      (Omni.BandedHermitian _, Omni.Full _) -> squareUnpacked a b
      (Omni.BandedHermitian _, Omni.LowerTriangular _) ->
         case ArrMatrix.superBandsSingleton a of
            Unary.Succ -> squareUnpacked a (OmniMatrix.toFull b)
            Unary.Zero ->
               Triangular.takeLower $ squareFull a (OmniMatrix.toFull b)
      (Omni.BandedHermitian _, Omni.UpperTriangular _) ->
         case ArrMatrix.subBandsSingleton a of
            Unary.Succ -> squareUnpacked a (OmniMatrix.toFull b)
            Unary.Zero ->
               Triangular.takeUpper $ squareFull a (OmniMatrix.toFull b)
      (Omni.Banded _shA, Omni.Banded _shB) -> bandedBanded a b
      (Omni.Banded _shA, Omni.UnitBandedTriangular _shB) ->
         bandedBanded a (Banded.noUnit b)
      (Omni.Banded _shA, Omni.BandedHermitian _shB) ->
         bandedBanded a (BandedHermitian.toBanded b)
      (Omni.UnitBandedTriangular _shA, Omni.Banded _shB) ->
         bandedBanded (Banded.noUnit a) b
      (Omni.UnitBandedTriangular _shA, Omni.BandedHermitian _shB) ->
         bandedBanded (Banded.noUnit a) (BandedHermitian.toBanded b)
      (Omni.UnitBandedTriangular shA, Omni.UnitBandedTriangular shB) ->
         case (Omni.bandedTriangularSingleton shA,
               Omni.bandedTriangularSingleton shB) of
            (Omni.BandedDiagonal, Omni.BandedDiagonal) -> b
            (Omni.BandedDiagonal, Omni.BandedUpper) -> b
            (Omni.BandedDiagonal, Omni.BandedLower) -> b
            (Omni.BandedUpper, Omni.BandedDiagonal) -> a
            (Omni.BandedLower, Omni.BandedDiagonal) -> a
            (Omni.BandedUpper, Omni.BandedLower) ->
               bandedBanded (Banded.noUnit a) (Banded.noUnit b)
            (Omni.BandedLower, Omni.BandedUpper) ->
               bandedBanded (Banded.noUnit a) (Banded.noUnit b)
            (Omni.BandedUpper, Omni.BandedUpper) -> bandedUpper a b
            (Omni.BandedLower, Omni.BandedLower) -> bandedLower a b
      (Omni.BandedHermitian _, Omni.Symmetric _) ->
         squareUnpacked a (Symmetric.toSquare b)
      (Omni.BandedHermitian _, Omni.Hermitian _) ->
         squareUnpacked a (Hermitian.toSquare b)
      (Omni.BandedHermitian _shA, Omni.Banded _shB) ->
         bandedBanded (BandedHermitian.toBanded a) b
      (Omni.BandedHermitian _shA, Omni.UnitBandedTriangular _shB) ->
         bandedBanded (BandedHermitian.toBanded a) (Banded.noUnit b)
      (Omni.BandedHermitian _shA, Omni.BandedHermitian _shB) ->
         bandedBanded (BandedHermitian.toBanded a) (BandedHermitian.toBanded b)

asArbitrary ::
   Id (ArrayMatrix pack
         Omni.Arbitrary Filled Filled meas vert horiz height width a)
asArbitrary = id


unpackedSquare ::
   (Layout.Packing packB) =>
   (Omni.Property propertyA, Omni.Strip lowerA, Omni.Strip upperA) =>
   (Omni.Property propertyB, Omni.Strip lowerB, Omni.Strip upperB) =>
   (Omni.Property propertyC, Omni.Strip lowerC, Omni.Strip upperC) =>
   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
   Full.Unpacked propertyA lowerA upperA meas vert horiz height width a ->
   ArrMatrix.Quadratic packB propertyB lowerB upperB width a ->
   Full.Unpacked propertyC lowerC upperC meas vert horiz height width a
unpackedSquare a b = flex $ fullSquare (unflex a) b

squareUnpacked ::
   (Layout.Packing packA) =>
   (Omni.Property propertyA, Omni.Strip lowerA, Omni.Strip upperA) =>
   (Omni.Property propertyB, Omni.Strip lowerB, Omni.Strip upperB) =>
   (Omni.Property propertyC, Omni.Strip lowerC, Omni.Strip upperC) =>
   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
   (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
   ArrMatrix.Quadratic packA propertyA lowerA upperA height a ->
   Full.Unpacked propertyB lowerB upperB meas vert horiz height width a ->
   Full.Unpacked propertyC lowerC upperC meas vert horiz height width a
squareUnpacked a b = flex $ squareFull a (unflex b)

unpackedUnpacked ::
   (Omni.Property propertyA, Omni.Strip lowerA, Omni.Strip upperA) =>
   (Omni.Property propertyB, Omni.Strip lowerB, Omni.Strip upperB) =>
   (Omni.Property propertyC, Omni.Strip lowerC, Omni.Strip upperC) =>
   (Extent.Measure measA, Extent.C vertA, Extent.C horizA) =>
   (Extent.Measure measB, Extent.C vertB, Extent.C horizB) =>
   (Shape.C height, Shape.C fuse, Eq fuse, Shape.C width, Class.Floating a) =>
   Full.Unpacked propertyA lowerA upperA measA vertA horizA height fuse a ->
   Full.Unpacked propertyB lowerB upperB measB vertB horizB fuse width a ->
   Full.Unpacked propertyC lowerC upperC
      (ExtentPriv.MultiplyMeasure measA measB)
      (ExtentPriv.Multiply vertA vertB)
      (ExtentPriv.Multiply horizA horizB)
      height width a
unpackedUnpacked a b =
   case unifyFactors a b of
      ((ExtentPriv.MeasureFact, ExtentPriv.TagFact, ExtentPriv.TagFact),
       (unifyLeft, unifyRight)) ->
         ArrMatrix.liftUnpacked2 FullBasic.multiply
            (Unpacked.mapExtent unifyLeft a)
            (Unpacked.mapExtent unifyRight b)

unpackedBanded ::
   (Omni.Property propertyA, Omni.Strip lowerA, Omni.Strip upperA) =>
   (Omni.Property propertyC, Omni.Strip lowerC, Omni.Strip upperC) =>
   (Unary.Natural sub, Unary.Natural super) =>
   (Extent.Measure measA, Extent.C vertA, Extent.C horizA) =>
   (Extent.Measure measB, Extent.C vertB, Extent.C horizB) =>
   (Shape.C height, Shape.C fuse, Eq fuse, Shape.C width, Class.Floating a) =>
   Full.Unpacked propertyA lowerA upperA measA vertA horizA height fuse a ->
   Banded sub super measB vertB horizB fuse width a ->
   Full.Unpacked propertyC lowerC upperC
      (ExtentPriv.MultiplyMeasure measA measB)
      (ExtentPriv.Multiply vertA vertB)
      (ExtentPriv.Multiply horizA horizB)
      height width a
unpackedBanded a b =
   case unifyFactors a b of
      ((ExtentPriv.MeasureFact, ExtentPriv.TagFact, ExtentPriv.TagFact),
       (unifyLeft, unifyRight)) ->
         Unpacked.swapMultiply
            (ArrMatrix.liftUnpacked1 . BandedBasic.multiplyFull .
             ArrMatrix.toVector . Banded.transpose)
            (Unpacked.mapExtent unifyLeft a)
            (Banded.mapExtent unifyRight b)

bandedUnpacked ::
   (Omni.Property propertyB, Omni.Strip lowerB, Omni.Strip upperB) =>
   (Omni.Property propertyC, Omni.Strip lowerC, Omni.Strip upperC) =>
   (Unary.Natural sub, Unary.Natural super) =>
   (Extent.Measure measA, Extent.C vertA, Extent.C horizA) =>
   (Extent.Measure measB, Extent.C vertB, Extent.C horizB) =>
   (Shape.C height, Shape.C fuse, Eq fuse, Shape.C width, Class.Floating a) =>
   Banded sub super measA vertA horizA height fuse a ->
   Full.Unpacked propertyB lowerB upperB measB vertB horizB fuse width a ->
   Full.Unpacked propertyC lowerC upperC
      (ExtentPriv.MultiplyMeasure measA measB)
      (ExtentPriv.Multiply vertA vertB)
      (ExtentPriv.Multiply horizA horizB)
      height width a
bandedUnpacked a b =
   case unifyFactors a b of
      ((ExtentPriv.MeasureFact, ExtentPriv.TagFact, ExtentPriv.TagFact),
       (unifyLeft, unifyRight)) ->
         (ArrMatrix.liftUnpacked1 . BandedBasic.multiplyFull .
          ArrMatrix.toVector)
            (Banded.mapExtent unifyLeft a)
            (Unpacked.mapExtent unifyRight b)


bandedBanded ::
   (Unary.Natural subA, Unary.Natural superA) =>
   (Unary.Natural subB, Unary.Natural superB) =>
   (Extent.Measure measA, Extent.Measure measB) =>
   (Extent.C vertA, Extent.C horizA) =>
   (Extent.C vertB, Extent.C horizB) =>
   (Shape.C height, Shape.C fuse, Eq fuse, Shape.C width, Class.Floating a) =>
   Banded subA superA measA vertA horizA height fuse a ->
   Banded subB superB measB vertB horizB fuse width a ->
   Banded
      (subA :+: subB) (superA :+: superB)
      (ExtentPriv.MultiplyMeasure measA measB)
      (ExtentPriv.Multiply vertA vertB)
      (ExtentPriv.Multiply horizA horizB)
      height width a
bandedBanded a b =
   case unifyFactors a b of
      ((ExtentPriv.MeasureFact, ExtentPriv.TagFact, ExtentPriv.TagFact),
       (unifyLeft, unifyRight)) ->
         Banded.multiply
            (Banded.mapExtent unifyLeft a)
            (Banded.mapExtent unifyRight b)

bandedUpper ::
   (Unary.Natural superA, Unary.Natural superB,
    (Unary.Succ superA :+: superB) ~ superC,
    Shape.C size, Eq size, Class.Floating a) =>
   Banded.UnitUpper (Unary.Succ superA) size a ->
   Banded.UnitUpper (Unary.Succ superB) size a ->
   Banded.UnitUpper (Unary.Succ superC) size a
bandedUpper a b =
   case (snd $ MatrixShape.bandedOffDiagonals $ ArrMatrix.shape a,
         snd $ MatrixShape.bandedOffDiagonals $ ArrMatrix.shape b) of
      (superA,superB) ->
         case Proof.addNat (natFromProxy superA) (natFromProxy $ prev superB) of
            Proof.Nat -> ArrMatrix.lift2 BandedBasic.multiply a b

bandedLower ::
   (Unary.Natural subA, Unary.Natural subB,
    (Unary.Succ subA :+: subB) ~ subC,
    Shape.C size, Eq size, Class.Floating a) =>
   Banded.UnitLower (Unary.Succ subA) size a ->
   Banded.UnitLower (Unary.Succ subB) size a ->
   Banded.UnitLower (Unary.Succ subC) size a
bandedLower a b =
   case (fst $ MatrixShape.bandedOffDiagonals $ ArrMatrix.shape a,
         fst $ MatrixShape.bandedOffDiagonals $ ArrMatrix.shape b) of
      (subA,subB) ->
         case Proof.addNat (natFromProxy subA) (natFromProxy $ prev subB) of
            Proof.Nat -> ArrMatrix.lift2 BandedBasic.multiply a b

prev :: UnaryProxy (Succ x) -> UnaryProxy x
prev Proxy = Proxy


unifyFactors ::
   (Matrix.Box typA, Matrix.Box typB) =>
   (Extent.Measure measA, Extent.C vertA, Extent.C horizA,
    Extent.Measure measB, Extent.C vertB, Extent.C horizB) =>
   (ExtentPriv.MultiplyMeasure measA measB ~ measC) =>
   (ExtentPriv.Multiply vertA vertB ~ vertC) =>
   (ExtentPriv.Multiply horizA horizB ~ horizC) =>
   Matrix typA xlA xuA lowerA upperA measA vertA horizA height fuse a ->
   Matrix typB xlB xuB lowerB upperB measB vertB horizB fuse width a ->
   ((ExtentPriv.MeasureFact measC,
     ExtentPriv.TagFact vertC, ExtentPriv.TagFact horizC),
    (Extent.Map measA vertA horizA measC vertC horizC height fuse,
     Extent.Map measB vertB horizB measC vertC horizC fuse width))
unifyFactors a b = ExtentStrict.unifiers (Matrix.extent a) (Matrix.extent b)


unflex ::
   (Omni.Property property, Omni.Strip lower, Omni.Strip upper) =>
   Full.Unpacked property lower upper meas vert horiz height width a ->
   ArrMatrix.Full meas vert horiz height width a
unflex = ArrMatrix.liftUnpacked1 id

flex ::
   (Omni.Property property, Omni.Strip lower, Omni.Strip upper) =>
   ArrMatrix.Full meas vert horiz height width a ->
   Full.Unpacked property lower upper meas vert horiz height width a
flex = ArrMatrix.liftUnpacked1 id


property ::
   (Omni.Property property) =>
   ArrayMatrix pack property lower upper meas vert horiz height width a ->
   Omni.PropertySingleton property
property _ = Omni.propertySingleton