{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE EmptyDataDecls #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE UndecidableInstances #-}
module Numeric.LAPACK.Matrix.Inverse where

import qualified Numeric.LAPACK.Matrix.Type.Private as Matrix
import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
import qualified Numeric.LAPACK.Matrix.Shape.Omni as Omni
import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
import qualified Numeric.LAPACK.Matrix.Class as MatrixClass
import qualified Numeric.LAPACK.Matrix.Divide as Divide
import qualified Numeric.LAPACK.Matrix.Multiply as Multiply
import Numeric.LAPACK.Matrix.Divide ((#\|), (-/#))
import Numeric.LAPACK.Matrix.Type.Private (Matrix)

import qualified Type.Data.Num.Unary as Unary


{- |
You may wrap non-diagonal band matrices in 'Wrapper.FillStrips' first
in order to meet the 'Omni.PowerStrip' constraint.
-}
data Inverse typ
data instance
   Matrix (Inverse typ)
      extraLower extraUpper lowerf upperf meas vert horiz height width a
         where
      Inverse ::
         Matrix.QuadraticMeas typ xl xu upper lower meas width height a ->
         Matrix.QuadraticMeas (Inverse typ) (xl,lower) (xu,upper)
            lower upper meas height width a

type family Fill offDiag
type instance Fill (Layout.Bands Unary.Zero) = Layout.Bands Unary.Zero
type instance Fill (Layout.Bands (Unary.Succ k)) = Layout.Filled
type instance Fill Layout.Filled = Layout.Filled

type family InverseExtra xl
type instance InverseExtra (xl,lower) = xl
type family InverseStrip xl
type instance InverseStrip (xl,lower) = lower

data PowerStripFact c = (Omni.PowerStrip c) => PowerStripFact

filledPowerStripFact ::
   (Omni.Strip c) => Omni.StripSingleton c -> PowerStripFact (Fill c)
filledPowerStripFact w =
   case w of
      Omni.StripFilled -> PowerStripFact
      Omni.StripBands Unary.Zero -> PowerStripFact
      Omni.StripBands Unary.Succ -> PowerStripFact



instance (Matrix.Transpose typ) => Matrix.Transpose (Inverse typ) where
   type TransposeExtra (Inverse typ) extra =
         Matrix.TransposeExtra typ (InverseExtra extra)
   transpose (Inverse a) = Inverse $ Matrix.transpose a

instance (Matrix.MultiplySame typ) => Matrix.MultiplySame (Inverse typ) where
   type MultiplySameExtra (Inverse typ) extra =
         (Matrix.MultiplySameExtra typ (InverseExtra extra),
          MatrixShape.PowerStrip (InverseStrip extra))
   multiplySame (Inverse a) (Inverse b) = Inverse $ Matrix.multiplySame b a

instance (Matrix.Box typ) => Matrix.Box (Inverse typ) where
   type BoxExtra (Inverse typ) extra = Matrix.BoxExtra typ (InverseExtra extra)
   extent (Inverse m) = Extent.transpose $ Matrix.extent m
   height (Inverse m) = Matrix.width m
   width (Inverse m) = Matrix.height m

instance (Matrix.ToQuadratic typ) => Matrix.ToQuadratic (Inverse typ) where
   heightToQuadratic (Inverse m) = Inverse $ Matrix.widthToQuadratic m
   widthToQuadratic (Inverse m) = Inverse $ Matrix.heightToQuadratic m

instance (MatrixClass.Complex typ) => MatrixClass.Complex (Inverse typ) where
   conjugate (Inverse m) = Inverse $ MatrixClass.conjugate m
   fromReal (Inverse m) = Inverse $ MatrixClass.fromReal m
   toComplex (Inverse m) = Inverse $ MatrixClass.toComplex m


instance
   (Divide.Solve typ, Matrix.ToQuadratic typ) =>
      Multiply.MultiplyVector (Inverse typ) where
   type MultiplyVectorExtra (Inverse typ) extra =
         (Multiply.MultiplyVectorExtra typ (InverseExtra extra),
          Divide.SolveExtra typ (InverseExtra extra),
          Matrix.BoxExtra typ (InverseExtra extra),
          Omni.Strip (InverseStrip extra))
   matrixVector (Inverse a) x = a#\|x
   vectorMatrix x (Inverse a) = x-/#a

instance
   (Divide.Solve typ, Matrix.ToQuadratic typ) =>
      Multiply.MultiplySquare (Inverse typ) where
   type MultiplySquareExtra (Inverse typ) extra =
         (Multiply.MultiplySquareExtra typ (InverseExtra extra),
          Divide.SolveExtra typ (InverseExtra extra),
          Matrix.BoxExtra typ (InverseExtra extra),
          Omni.Strip (InverseStrip extra))
   transposableSquare trans (Inverse a) = Divide.solve trans a
   squareFull (Inverse a) b = Divide.solveRight a b
   fullSquare b (Inverse a) = Divide.solveLeft b a

instance (Multiply.Power typ) => Multiply.Power (Inverse typ) where
   type PowerExtra (Inverse typ) extra =
         (Multiply.PowerExtra typ (InverseExtra extra),
          MatrixShape.PowerStrip (InverseStrip extra))
   square (Inverse a) = Inverse $ Multiply.square a
   power n (Inverse a) = Inverse $ Multiply.power n a
   powers1 (Inverse a) = fmap Inverse $ Multiply.powers1 a


instance (Divide.Determinant typ) => Divide.Determinant (Inverse typ) where
   type DeterminantExtra (Inverse typ) extra =
         (Divide.DeterminantExtra typ (InverseExtra extra),
          MatrixShape.PowerStrip (InverseStrip extra))
   determinant (Inverse a) = recip $ Divide.determinant a

instance
   (Multiply.MultiplySquare typ, Matrix.ToQuadratic typ) =>
      Divide.Solve (Inverse typ) where
   type SolveExtra (Inverse typ) extra =
         (Divide.SolveExtra typ (InverseExtra extra),
          Multiply.MultiplySquareExtra typ (InverseExtra extra),
          MatrixShape.PowerStrip (InverseStrip extra))
   solve trans (Inverse a) = Multiply.transposableSquare trans a
   solveRight (Inverse a) b = Multiply.squareFull a b
   solveLeft b (Inverse a) = Multiply.fullSquare b a

instance
   (Divide.Inverse typ, Multiply.MultiplySquare typ, Matrix.ToQuadratic typ) =>
      Divide.Inverse (Inverse typ) where
   type InverseExtra (Inverse typ) extra =
         (Divide.InverseExtra typ (InverseExtra extra),
          Multiply.MultiplySquareExtra typ (InverseExtra extra),
          MatrixShape.PowerStrip (InverseStrip extra))
   inverse (Inverse a) = Inverse $ Divide.inverse a