{-# LANGUAGE TypeFamilies #-}
module Numeric.LAPACK.Matrix.Divide where

import qualified Numeric.LAPACK.Matrix.Plain.Divide as ArrDivide
import qualified Numeric.LAPACK.Matrix.Array.Basic as Basic
import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
import qualified Numeric.LAPACK.Matrix.Permutation as PermMatrix
import qualified Numeric.LAPACK.Matrix.Multiply as Multiply
import qualified Numeric.LAPACK.Matrix.Type as Type
import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
import qualified Numeric.LAPACK.Vector as Vector
import Numeric.LAPACK.Matrix.Array (Full)
import Numeric.LAPACK.Matrix.Type (Matrix, scaleWithCheck)
import Numeric.LAPACK.Matrix.Modifier
         (Transposition(NonTransposed,Transposed),
          Inversion(Inverted))
import Numeric.LAPACK.Vector (Vector)

import qualified Numeric.Netlib.Class as Class

import qualified Data.Array.Comfort.Shape as Shape

import Data.Semigroup ((<>))


class
   (Type.Box typ, Type.HeightOf typ ~ Type.WidthOf typ) =>
      Determinant typ where
   determinant :: (Class.Floating a) => Matrix typ a -> a

class (Type.Box typ, Type.HeightOf typ ~ Type.WidthOf typ) => Solve typ where
   {-# MINIMAL solve | solveLeft,solveRight #-}
   solve ::
      (Type.HeightOf typ ~ height, Eq height, Shape.C width,
       Extent.C vert, Extent.C horiz, Class.Floating a) =>
      Transposition -> Matrix typ a ->
      Full vert horiz height width a -> Full vert horiz height width a
   solve Transposition
NonTransposed Matrix typ a
a Full vert horiz height width a
b = Matrix typ a
-> Full vert horiz height width a -> Full vert horiz height width a
forall typ height width vert horiz a.
(Solve typ, HeightOf typ ~ height, Eq height, C width, C vert,
 C horiz, Floating a) =>
Matrix typ a
-> Full vert horiz height width a -> Full vert horiz height width a
solveRight Matrix typ a
a Full vert horiz height width a
b
   solve Transposition
Transposed Matrix typ a
a Full vert horiz height width a
b = Full horiz vert width height a -> Full vert horiz height width a
forall vert horiz height width a.
(C vert, C horiz) =>
Full vert horiz height width a -> Full horiz vert width height a
Basic.transpose (Full horiz vert width height a -> Full vert horiz height width a)
-> Full horiz vert width height a -> Full vert horiz height width a
forall a b. (a -> b) -> a -> b
$ Full horiz vert width height a
-> Matrix typ a -> Full horiz vert width height a
forall typ width height vert horiz a.
(Solve typ, WidthOf typ ~ width, Eq width, C height, C vert,
 C horiz, Floating a) =>
Full vert horiz height width a
-> Matrix typ a -> Full vert horiz height width a
solveLeft (Full vert horiz height width a -> Full horiz vert width height a
forall vert horiz height width a.
(C vert, C horiz) =>
Full vert horiz height width a -> Full horiz vert width height a
Basic.transpose Full vert horiz height width a
b) Matrix typ a
a

   solveRight ::
      (Type.HeightOf typ ~ height, Eq height, Shape.C width,
       Extent.C vert, Extent.C horiz, Class.Floating a) =>
      Matrix typ a ->
      Full vert horiz height width a -> Full vert horiz height width a
   solveRight = Transposition
-> Matrix typ a
-> Full vert horiz height width a
-> Full vert horiz height width a
forall typ height width vert horiz a.
(Solve typ, HeightOf typ ~ height, Eq height, C width, C vert,
 C horiz, Floating a) =>
Transposition
-> Matrix typ a
-> Full vert horiz height width a
-> Full vert horiz height width a
solve Transposition
NonTransposed

   solveLeft ::
      (Type.WidthOf typ ~ width, Eq width, Shape.C height,
       Extent.C vert, Extent.C horiz, Class.Floating a) =>
      Full vert horiz height width a ->
      Matrix typ a -> Full vert horiz height width a
   solveLeft = (Matrix typ a
 -> Full horiz vert width height a
 -> Full horiz vert width height a)
-> Full vert horiz height width a
-> Matrix typ a
-> Full vert horiz height width a
forall vertA vertB horizA horizB matrix widthA heightA a widthB
       heightB.
(C vertA, C vertB, C horizA, C horizB) =>
(matrix
 -> Full horizA vertA widthA heightA a
 -> Full horizB vertB widthB heightB a)
-> Full vertA horizA heightA widthA a
-> matrix
-> Full vertB horizB heightB widthB a
Basic.swapMultiply ((Matrix typ a
  -> Full horiz vert width height a
  -> Full horiz vert width height a)
 -> Full vert horiz height width a
 -> Matrix typ a
 -> Full vert horiz height width a)
-> (Matrix typ a
    -> Full horiz vert width height a
    -> Full horiz vert width height a)
-> Full vert horiz height width a
-> Matrix typ a
-> Full vert horiz height width a
forall a b. (a -> b) -> a -> b
$ Transposition
-> Matrix typ a
-> Full horiz vert width height a
-> Full horiz vert width height a
forall typ height width vert horiz a.
(Solve typ, HeightOf typ ~ height, Eq height, C width, C vert,
 C horiz, Floating a) =>
Transposition
-> Matrix typ a
-> Full vert horiz height width a
-> Full vert horiz height width a
solve Transposition
Transposed

class (Solve typ, Multiply.Power typ) => Inverse typ where
   inverse :: (Class.Floating a) => Matrix typ a -> Matrix typ a

infixl 7 ##/#
infixr 7 #\##

(#\##) ::
   (Solve typ, Type.HeightOf typ ~ height, Eq height, Shape.C width,
    Extent.C vert, Extent.C horiz, Class.Floating a) =>
   Matrix typ a ->
   Full vert horiz height width a -> Full vert horiz height width a
#\## :: Matrix typ a
-> Full vert horiz height width a -> Full vert horiz height width a
(#\##) = Matrix typ a
-> Full vert horiz height width a -> Full vert horiz height width a
forall typ height width vert horiz a.
(Solve typ, HeightOf typ ~ height, Eq height, C width, C vert,
 C horiz, Floating a) =>
Matrix typ a
-> Full vert horiz height width a -> Full vert horiz height width a
solveRight

(##/#) ::
   (Solve typ, Type.WidthOf typ ~ width, Eq width, Shape.C height,
    Extent.C vert, Extent.C horiz, Class.Floating a) =>
   Full vert horiz height width a ->
   Matrix typ a -> Full vert horiz height width a
##/# :: Full vert horiz height width a
-> Matrix typ a -> Full vert horiz height width a
(##/#) = Full vert horiz height width a
-> Matrix typ a -> Full vert horiz height width a
forall typ width height vert horiz a.
(Solve typ, WidthOf typ ~ width, Eq width, C height, C vert,
 C horiz, Floating a) =>
Full vert horiz height width a
-> Matrix typ a -> Full vert horiz height width a
solveLeft


solveVector ::
   (Solve typ, Type.HeightOf typ ~ height, Eq height, Class.Floating a) =>
   Transposition -> Matrix typ a -> Vector height a -> Vector height a
solveVector :: Transposition -> Matrix typ a -> Vector height a -> Vector height a
solveVector Transposition
trans =
   Order
-> (General height () a -> General height () a)
-> Vector height a
-> Vector height a
forall height0 a height1 b.
Order
-> (General height0 () a -> General height1 () b)
-> Vector height0 a
-> Vector height1 b
ArrMatrix.unliftColumn Order
MatrixShape.ColumnMajor ((General height () a -> General height () a)
 -> Vector height a -> Vector height a)
-> (Matrix typ a -> General height () a -> General height () a)
-> Matrix typ a
-> Vector height a
-> Vector height a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Transposition
-> Matrix typ a -> General height () a -> General height () a
forall typ height width vert horiz a.
(Solve typ, HeightOf typ ~ height, Eq height, C width, C vert,
 C horiz, Floating a) =>
Transposition
-> Matrix typ a
-> Full vert horiz height width a
-> Full vert horiz height width a
solve Transposition
trans

infixl 7 -/#
infixr 7 #\|

(#\|) ::
   (Solve typ, Type.HeightOf typ ~ height, Eq height, Class.Floating a) =>
   Matrix typ a -> Vector height a -> Vector height a
#\| :: Matrix typ a -> Vector height a -> Vector height a
(#\|) = Transposition -> Matrix typ a -> Vector height a -> Vector height a
forall typ height a.
(Solve typ, HeightOf typ ~ height, Eq height, Floating a) =>
Transposition -> Matrix typ a -> Vector height a -> Vector height a
solveVector Transposition
NonTransposed

(-/#) ::
   (Solve typ, Type.HeightOf typ ~ height, Eq height, Class.Floating a) =>
   Vector height a -> Matrix typ a -> Vector height a
-/# :: Vector height a -> Matrix typ a -> Vector height a
(-/#) = (Matrix typ a -> Vector height a -> Vector height a)
-> Vector height a -> Matrix typ a -> Vector height a
forall a b c. (a -> b -> c) -> b -> a -> c
flip ((Matrix typ a -> Vector height a -> Vector height a)
 -> Vector height a -> Matrix typ a -> Vector height a)
-> (Matrix typ a -> Vector height a -> Vector height a)
-> Vector height a
-> Matrix typ a
-> Vector height a
forall a b. (a -> b) -> a -> b
$ Transposition -> Matrix typ a -> Vector height a -> Vector height a
forall typ height a.
(Solve typ, HeightOf typ ~ height, Eq height, Floating a) =>
Transposition -> Matrix typ a -> Vector height a -> Vector height a
solveVector Transposition
Transposed


instance (Shape.C shape, Eq shape) => Determinant (Type.Scale shape) where
   determinant :: Matrix (Scale shape) a -> a
determinant (Type.Scale sh a) = a
a a -> Int -> a
forall a b. (Num a, Integral b) => a -> b -> a
^ shape -> Int
forall sh. C sh => sh -> Int
Shape.size shape
sh

instance (Shape.C shape, Eq shape) => Solve (Type.Scale shape) where
   solve :: Transposition
-> Matrix (Scale shape) a
-> Full vert horiz height width a
-> Full vert horiz height width a
solve Transposition
_trans =
      String
-> (Full vert horiz height width a -> shape)
-> (a
    -> Full vert horiz height width a
    -> Full vert horiz height width a)
-> Matrix (Scale shape) a
-> Full vert horiz height width a
-> Full vert horiz height width a
forall shape b a c.
Eq shape =>
String
-> (b -> shape)
-> (a -> b -> c)
-> Matrix (Scale shape) a
-> b
-> c
scaleWithCheck String
"Matrix.Scale.solve" Full vert horiz height width a -> shape
forall typ a. Box typ => Matrix typ a -> HeightOf typ
Type.height ((a
  -> Full vert horiz height width a
  -> Full vert horiz height width a)
 -> Matrix (Scale shape) a
 -> Full vert horiz height width a
 -> Full vert horiz height width a)
-> (a
    -> Full vert horiz height width a
    -> Full vert horiz height width a)
-> Matrix (Scale shape) a
-> Full vert horiz height width a
-> Full vert horiz height width a
forall a b. (a -> b) -> a -> b
$
         (Array (Full vert horiz height width) a
 -> Array (Full vert horiz height width) a)
-> Full vert horiz height width a -> Full vert horiz height width a
forall shA a shB b.
(Array shA a -> Array shB b)
-> ArrayMatrix shA a -> ArrayMatrix shB b
ArrMatrix.lift1 ((Array (Full vert horiz height width) a
  -> Array (Full vert horiz height width) a)
 -> Full vert horiz height width a
 -> Full vert horiz height width a)
-> (a
    -> Array (Full vert horiz height width) a
    -> Array (Full vert horiz height width) a)
-> a
-> Full vert horiz height width a
-> Full vert horiz height width a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a
-> Array (Full vert horiz height width) a
-> Array (Full vert horiz height width) a
forall sh a. (C sh, Floating a) => a -> Vector sh a -> Vector sh a
Vector.scale (a
 -> Array (Full vert horiz height width) a
 -> Array (Full vert horiz height width) a)
-> (a -> a)
-> a
-> Array (Full vert horiz height width) a
-> Array (Full vert horiz height width) a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> a
forall a. Fractional a => a -> a
recip

instance (Shape.C shape, Eq shape) => Inverse (Type.Scale shape) where
   inverse :: Matrix (Scale shape) a -> Matrix (Scale shape) a
inverse (Type.Scale shape a) = shape -> a -> Matrix (Scale shape) a
forall shape a. shape -> a -> Matrix (Scale shape) a
Type.Scale shape
shape (a -> Matrix (Scale shape) a) -> a -> Matrix (Scale shape) a
forall a b. (a -> b) -> a -> b
$ a -> a
forall a. Fractional a => a -> a
recip a
a


instance (Shape.C shape) => Determinant (PermMatrix.Permutation shape) where
   determinant :: Matrix (Permutation shape) a -> a
determinant = Matrix (Permutation shape) a -> a
forall sh a. (C sh, Floating a) => Matrix (Permutation sh) a -> a
PermMatrix.determinant

instance (Shape.C shape) => Solve (PermMatrix.Permutation shape) where
   solve :: Transposition
-> Matrix (Permutation shape) a
-> Full vert horiz height width a
-> Full vert horiz height width a
solve Transposition
trans =
      Inversion
-> Matrix (Permutation height) a
-> Full vert horiz height width a
-> Full vert horiz height width a
forall vert horiz height width a.
(C vert, C horiz, C height, Eq height, C width, Floating a) =>
Inversion
-> Matrix (Permutation height) a
-> Full vert horiz height width a
-> Full vert horiz height width a
PermMatrix.multiplyFull
         (Inversion
Inverted Inversion -> Inversion -> Inversion
forall a. Semigroup a => a -> a -> a
<> Transposition -> Inversion
PermMatrix.inversionFromTransposition Transposition
trans)

instance (Shape.C shape) => Inverse (PermMatrix.Permutation shape) where
   inverse :: Matrix (Permutation shape) a -> Matrix (Permutation shape) a
inverse = Matrix (Permutation shape) a -> Matrix (Permutation shape) a
forall sh a.
C sh =>
Matrix (Permutation sh) a -> Matrix (Permutation sh) a
PermMatrix.transpose


instance
      (ArrDivide.Determinant shape) => Determinant (ArrMatrix.Array shape) where
   determinant :: Matrix (Array shape) a -> a
determinant = Array shape a -> a
forall shape a.
(Determinant shape, Floating a) =>
Array shape a -> a
ArrDivide.determinant (Array shape a -> a)
-> (Matrix (Array shape) a -> Array shape a)
-> Matrix (Array shape) a
-> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix (Array shape) a -> Array shape a
forall sh a. ArrayMatrix sh a -> Array sh a
ArrMatrix.toVector

instance (ArrDivide.Solve shape) => Solve (ArrMatrix.Array shape) where
   solve :: Transposition
-> Matrix (Array shape) a
-> Full vert horiz height width a
-> Full vert horiz height width a
solve = (Array shape a
 -> Array (Full vert horiz height width) a
 -> Array (Full vert horiz height width) a)
-> Matrix (Array shape) a
-> Full vert horiz height width a
-> Full vert horiz height width a
forall shA a shB b shC c.
(Array shA a -> Array shB b -> Array shC c)
-> ArrayMatrix shA a -> ArrayMatrix shB b -> ArrayMatrix shC c
ArrMatrix.lift2 ((Array shape a
  -> Array (Full vert horiz height width) a
  -> Array (Full vert horiz height width) a)
 -> Matrix (Array shape) a
 -> Full vert horiz height width a
 -> Full vert horiz height width a)
-> (Transposition
    -> Array shape a
    -> Array (Full vert horiz height width) a
    -> Array (Full vert horiz height width) a)
-> Transposition
-> Matrix (Array shape) a
-> Full vert horiz height width a
-> Full vert horiz height width a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Transposition
-> Array shape a
-> Array (Full vert horiz height width) a
-> Array (Full vert horiz height width) a
forall shape a height vert horiz width.
(Solve shape, Floating a, HeightOf shape ~ height, Eq height,
 C vert, C horiz, C width) =>
Transposition
-> Array shape a
-> Full vert horiz height width a
-> Full vert horiz height width a
ArrDivide.solve
   solveLeft :: Full vert horiz height width a
-> Matrix (Array shape) a -> Full vert horiz height width a
solveLeft = (Array (Full vert horiz height width) a
 -> Array shape a -> Array (Full vert horiz height width) a)
-> Full vert horiz height width a
-> Matrix (Array shape) a
-> Full vert horiz height width a
forall shA a shB b shC c.
(Array shA a -> Array shB b -> Array shC c)
-> ArrayMatrix shA a -> ArrayMatrix shB b -> ArrayMatrix shC c
ArrMatrix.lift2 Array (Full vert horiz height width) a
-> Array shape a -> Array (Full vert horiz height width) a
forall shape a width vert horiz height.
(Solve shape, Floating a, HeightOf shape ~ width, Eq width, C vert,
 C horiz, C height) =>
Full vert horiz height width a
-> Array shape a -> Full vert horiz height width a
ArrDivide.solveLeft
   solveRight :: Matrix (Array shape) a
-> Full vert horiz height width a -> Full vert horiz height width a
solveRight = (Array shape a
 -> Array (Full vert horiz height width) a
 -> Array (Full vert horiz height width) a)
-> Matrix (Array shape) a
-> Full vert horiz height width a
-> Full vert horiz height width a
forall shA a shB b shC c.
(Array shA a -> Array shB b -> Array shC c)
-> ArrayMatrix shA a -> ArrayMatrix shB b -> ArrayMatrix shC c
ArrMatrix.lift2 Array shape a
-> Array (Full vert horiz height width) a
-> Array (Full vert horiz height width) a
forall shape a height vert horiz width.
(Solve shape, Floating a, HeightOf shape ~ height, Eq height,
 C vert, C horiz, C width) =>
Array shape a
-> Full vert horiz height width a -> Full vert horiz height width a
ArrDivide.solveRight

instance (ArrDivide.Inverse shape) => Inverse (ArrMatrix.Array shape) where
   inverse :: Matrix (Array shape) a -> Matrix (Array shape) a
inverse = (Array shape a -> Array shape a)
-> Matrix (Array shape) a -> Matrix (Array shape) a
forall shA a shB b.
(Array shA a -> Array shB b)
-> ArrayMatrix shA a -> ArrayMatrix shB b
ArrMatrix.lift1 Array shape a -> Array shape a
forall shape a.
(Inverse shape, Floating a) =>
Array shape a -> Array shape a
ArrDivide.inverse