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

import qualified Numeric.LAPACK.Matrix.HermitianPositiveDefinite as HermitianPD
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.Square.Linear as SquareLin
import qualified Numeric.LAPACK.Matrix.Square as Square
import qualified Numeric.LAPACK.Matrix.Banded.Linear as BandedLin
import qualified Numeric.LAPACK.Matrix.BandedHermitianPositiveDefinite
                  as BandedHermitianPD
import qualified Numeric.LAPACK.Matrix.BandedHermitian as BandedHermitian
import qualified Numeric.LAPACK.Matrix.Banded as Banded

import qualified Numeric.LAPACK.Matrix.Array.Basic as OmniMatrix
import qualified Numeric.LAPACK.Matrix.Array.Private as ArrMatrix
import qualified Numeric.LAPACK.Matrix.Shape.Omni as Omni
import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
import qualified Numeric.LAPACK.Matrix.Extent as Extent
import qualified Numeric.LAPACK.Vector as Vector
import qualified Numeric.LAPACK.Scalar as Scalar
import Numeric.LAPACK.Matrix.Array.Multiply (Factor(..), factorSingleton)
import Numeric.LAPACK.Matrix.Array.Private (Full, Quadratic, QuadraticMeas)

import qualified Numeric.Netlib.Class as Class

import qualified Type.Data.Bool as TBool

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


determinant ::
   (Layout.Packing pack, Omni.Property property,
    Omni.Strip lower, Omni.Strip upper,
    Shape.C sh, Class.Floating a) =>
   Quadratic pack property lower upper sh a -> a
determinant a =
   let shape = ArrMatrix.shape a in
   case factorSingleton shape of
      FactorFull -> Square.determinant $ OmniMatrix.toFull a
      FactorUpperTriangular ->
         case ArrMatrix.diagTag a of
            MatrixShape.Unit -> Scalar.one
            MatrixShape.Arbitrary -> Triangular.determinant a
      FactorLowerTriangular ->
         case ArrMatrix.diagTag a of
            MatrixShape.Unit -> Scalar.one
            MatrixShape.Arbitrary -> Triangular.determinant a

      FactorSymmetric -> Symmetric.determinant a
      FactorHermitian -> Scalar.fromReal $ Hermitian.determinant a

      FactorBanded -> Banded.determinant a
      FactorUnitBandedTriangular -> Scalar.one
      FactorBandedHermitian ->
         case Omni.hermitianSet shape of
            (TBool.False, TBool.False, TBool.True) ->
               Scalar.fromReal $ BandedHermitianPD.determinant a
            (TBool.True, TBool.False, TBool.False) ->
               (OmniMatrix.signNegativeDeterminant shape *) $ Scalar.fromReal $
               BandedHermitianPD.determinant $ Hermitian.negate a
            _ -> Banded.determinant $ BandedHermitian.toBanded a


{- |
We use the solvers for positive definite Hermitian matrices
also for negative definite matrices.
We even use them for semidefinite matrices,
because semidefinite matrices with full rank are definite.
-}
solve ::
   (Layout.Packing pack, Omni.Property property,
    Omni.Strip lower, Omni.Strip upper,
    Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C sh, Eq sh, Shape.C nrhs, Class.Floating a) =>
   Quadratic pack property lower upper sh a ->
   Full meas vert horiz sh nrhs a -> Full meas vert horiz sh nrhs a
solve a =
   case factorSingleton $ ArrMatrix.shape a of
      FactorFull -> Square.solve $ OmniMatrix.toFull a
      FactorUpperTriangular -> Triangular.solve a
      FactorLowerTriangular -> Triangular.solve a
      FactorSymmetric -> Symmetric.solve a
      FactorHermitian -> Hermitian.solve a
      FactorBanded -> Banded.solve a
      FactorUnitBandedTriangular ->
         ArrMatrix.lift2 (BandedLin.solveTriangular MatrixShape.Unit) a
      FactorBandedHermitian ->
         case Omni.hermitianSet $ ArrMatrix.shape a of
            (TBool.False, _, TBool.True) ->
               BandedHermitianPD.solve (HermitianPD.assureFullRank a)
            (TBool.True, _, TBool.False) ->
               ArrMatrix.negate .
               BandedHermitianPD.solve
                  (Hermitian.negate $ HermitianPD.assureFullRank a)
            _ -> Banded.solve $ BandedHermitian.toBanded a


inverse ::
   (Layout.Packing pack, Omni.Property property,
    MatrixShape.PowerStrip lower, MatrixShape.PowerStrip upper,
    Extent.Measure meas, Shape.C height, Shape.C width, Class.Floating a) =>
   QuadraticMeas pack property lower upper meas height width a ->
   QuadraticMeas pack property lower upper meas width height a
inverse a =
   case Omni.powerSingleton $ ArrMatrix.shape a of
      Omni.PowerIdentity -> a
      Omni.PowerUpperTriangular -> Triangular.inverse a
      Omni.PowerLowerTriangular -> Triangular.inverse a
      Omni.PowerSymmetric -> Symmetric.inverse a
      Omni.PowerHermitian -> Hermitian.inverse a
      Omni.PowerFull -> ArrMatrix.liftUnpacked1 SquareLin.inverse a
      Omni.PowerDiagonal ->
         case ArrMatrix.shape a of
            Omni.Banded _ ->
               ArrMatrix.lift1
                  (Array.mapShape Layout.diagonalInverse . Vector.recip) a
            Omni.BandedHermitian _ -> ArrMatrix.lift1 Vector.recip a
            Omni.UnitBandedTriangular _ -> a
            Omni.Full _ -> ArrMatrix.liftUnpacked1 SquareLin.inverse a