{-# 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
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