Inverse.hs

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

import qualified Numeric.LAPACK.Matrix.Type 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 (Matrix)

import qualified Type.Data.Num.Unary as Unary


data Inverse typ
data instance
   Matrix (Inverse typ)
      extraLower extraUpper lowerf upperf meas vert horiz height width a
         where
      Inverse ::
         (Omni.Strip lower, Fill lower ~ lowerf, Omni.PowerStrip lowerf,
          Omni.Strip upper, Fill upper ~ upperf, Omni.PowerStrip upperf) =>
         Matrix.QuadraticMeas typ xl xu upper lower meas width height a ->
         Matrix.QuadraticMeas (Inverse typ) (xl,lower) (xu,upper)
            lowerf upperf meas height width a

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

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
   transpose (Inverse a) = Inverse $ Matrix.transpose a

instance
   (Matrix.MultiplySame typ xl xu,
    MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper) =>
      Matrix.MultiplySame (Inverse typ) (xl,lower) (xu,upper) where
   multiplySame (Inverse a) (Inverse b) = Inverse $ Matrix.multiplySame b a

instance (Matrix.Box typ) => Matrix.Box (Inverse typ) where
   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 xl xu, Matrix.ToQuadratic typ,
    Omni.Strip lower, Omni.Strip upper) =>
      Multiply.MultiplyVector (Inverse typ) (xl,lower) (xu,upper) where
   matrixVector (Inverse a) x = a#\|x
   vectorMatrix x (Inverse a) = x-/#a

instance
   (Divide.Solve typ xl xu, Matrix.ToQuadratic typ,
    Omni.Strip lower, Omni.Strip upper) =>
      Multiply.MultiplySquare (Inverse typ) (xl,lower) (xu,upper) where
   transposableSquare trans (Inverse a) = Divide.solve trans a
   squareFull (Inverse a) b = Divide.solveRight a b
   fullSquare b (Inverse a) = Divide.solveLeft b a

instance
   (Multiply.Power typ xl xu,
    MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper) =>
      Multiply.Power (Inverse typ) (xl,lower) (xu,upper) where
   square (Inverse a) = Inverse $ Multiply.square a
   power n (Inverse a) = Inverse $ Multiply.power n a
   powers1 (Inverse a) = fmap Inverse $ Multiply.powers1 a


instance
   (Divide.Determinant typ xl xu,
    MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper) =>
      Divide.Determinant (Inverse typ) (xl,lower) (xu,upper) where
   determinant (Inverse a) = recip $ Divide.determinant a

instance
   (Multiply.MultiplySquare typ xl xu, Matrix.ToQuadratic typ,
    MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper) =>
      Divide.Solve (Inverse typ) (xl,lower) (xu,upper) where
   solve trans (Inverse a) = Multiply.transposableSquare trans a
   solveRight (Inverse a) b = Multiply.squareFull a b
   solveLeft b (Inverse a) = Multiply.fullSquare b a

instance
   (Divide.Inverse typ xl xu, Multiply.MultiplySquare typ xl xu,
    Matrix.ToQuadratic typ, MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper) =>
      Divide.Inverse (Inverse typ) (xl,lower) (xu,upper) where
   inverse (Inverse a) = Inverse $ Divide.inverse a