module AERN2.Linear.Matrix.Inverse
(inverse)
where
import MixedTypesNumPrelude
import qualified Numeric.CollectErrors as CN
import AERN2.Linear.Matrix.Type
import qualified AERN2.Linear.Vector.Type as V
import Data.Maybe
import AERN2.MP.Float
import AERN2.MP.Dyadic
import AERN2.MP
inverse :: Matrix (CN MPBall) -> Maybe (Matrix (CN MPBall))
inverse :: Matrix (CN MPBall) -> Maybe (Matrix (CN MPBall))
inverse Matrix (CN MPBall)
m =
if Matrix (CN MPBall) -> Bool
forall es. CanTestErrorsPresent es => es -> Bool
CN.hasError Matrix (CN MPBall)
m
Bool -> Bool -> AndOrType Bool Bool
forall a b. CanAndOrAsymmetric a b => a -> b -> AndOrType a b
|| (Maybe (Matrix MPFloat) -> Bool
forall a. Maybe a -> Bool
isNothing Maybe (Matrix MPFloat)
maybeY)
Bool -> Bool -> AndOrType Bool Bool
forall a b. CanAndOrAsymmetric a b => a -> b -> AndOrType a b
|| (Bool -> NegType Bool
forall t. CanNeg t => t -> NegType t
not (Bool -> NegType Bool) -> Bool -> NegType Bool
forall a b. (a -> b) -> a -> b
$ CN MPBall
nerr CN MPBall -> Integer -> Bool
forall a b. HasOrderCertainlyAsymmetric a b => a -> b -> Bool
!<! Integer
1) then
Maybe (Matrix (CN MPBall))
forall a. Maybe a
Nothing
else
Matrix (CN MPBall) -> Maybe (Matrix (CN MPBall))
forall a. a -> Maybe a
Just (Matrix (CN MPBall) -> Maybe (Matrix (CN MPBall)))
-> Matrix (CN MPBall) -> Maybe (Matrix (CN MPBall))
forall a b. (a -> b) -> a -> b
$ Matrix (CN MPBall) -> Matrix (CN MPBall)
aux Matrix (CN MPBall)
x0
where
cm :: Matrix MPFloat
cm = Integer -> Vector MPFloat -> Matrix MPFloat
forall a. Integer -> Vector a -> Matrix a
Matrix (Matrix (CN MPBall) -> Integer
forall a. Matrix a -> Integer
width Matrix (CN MPBall)
m) ((CN MPBall -> MPFloat) -> Vector (CN MPBall) -> Vector MPFloat
forall a b. (a -> b) -> Vector a -> Vector b
V.map (\CN MPBall
x -> Dyadic -> MPFloat
forall t. CanBeMPFloat t => t -> MPFloat
mpFloat (Dyadic -> MPFloat) -> Dyadic -> MPFloat
forall a b. (a -> b) -> a -> b
$ MPBall -> CentreType MPBall
forall t. IsBall t => t -> CentreType t
centre (MPBall -> CentreType MPBall) -> MPBall -> CentreType MPBall
forall a b. (a -> b) -> a -> b
$ CN MPBall -> MPBall
forall p. CN p -> p
unCN CN MPBall
x) (Vector (CN MPBall) -> Vector MPFloat)
-> Vector (CN MPBall) -> Vector MPFloat
forall a b. (a -> b) -> a -> b
$ Matrix (CN MPBall) -> Vector (CN MPBall)
forall a. Matrix a -> Vector a
entries Matrix (CN MPBall)
m)
idM :: Matrix (CN MPBall)
idM = Integer -> Integer -> Matrix (CN MPBall)
forall a. HasIntegers a => Integer -> Integer -> Matrix a
identity (Matrix (CN MPBall) -> Integer
forall a. Matrix a -> Integer
width Matrix (CN MPBall)
m) (Matrix (CN MPBall) -> Integer
forall a. Matrix a -> Integer
width Matrix (CN MPBall)
m) :: Matrix (CN MPBall)
maybeY :: Maybe (Matrix MPFloat)
maybeY = Matrix MPFloat -> Maybe (Matrix MPFloat)
inverseGauss Matrix MPFloat
cm
y :: Matrix (CN MPBall)
y = (MPFloat -> CN MPBall) -> Matrix MPFloat -> Matrix (CN MPBall)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (MPBall -> CN MPBall
forall v. v -> CN v
cn (MPBall -> CN MPBall)
-> (MPFloat -> MPBall) -> MPFloat -> CN MPBall
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Dyadic -> MPBall
forall t. CanBeMPBall t => t -> MPBall
mpBall (Dyadic -> MPBall) -> (MPFloat -> Dyadic) -> MPFloat -> MPBall
forall b c a. (b -> c) -> (a -> b) -> a -> c
. MPFloat -> Dyadic
forall t. CanBeDyadic t => t -> Dyadic
dyadic) (Matrix MPFloat -> Matrix (CN MPBall))
-> Matrix MPFloat -> Matrix (CN MPBall)
forall a b. (a -> b) -> a -> b
$ Maybe (Matrix MPFloat) -> Matrix MPFloat
forall a. HasCallStack => Maybe a -> a
fromJust Maybe (Matrix MPFloat)
maybeY
yid :: MulType (Matrix (CN MPBall)) (Matrix (CN MPBall))
yid = Matrix (CN MPBall)
y Matrix (CN MPBall)
-> Matrix (CN MPBall)
-> MulType (Matrix (CN MPBall)) (Matrix (CN MPBall))
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
* Matrix (CN MPBall)
idM
nye :: Vector (CN MPBall)
nye = (Vector (CN MPBall) -> CN MPBall)
-> Vector (Vector (CN MPBall)) -> Vector (CN MPBall)
forall a b. (a -> b) -> Vector a -> Vector b
V.map Vector (CN MPBall) -> CN MPBall
forall a. (HasIntegers a, CanMinMaxSameType a) => Vector a -> a
V.inftyNorm (Matrix (CN MPBall) -> Vector (Vector (CN MPBall))
forall a. Matrix a -> Vector (Vector a)
columns MulType (Matrix (CN MPBall)) (Matrix (CN MPBall))
Matrix (CN MPBall)
yid)
err :: SubType (Matrix (CN MPBall)) (Matrix (CN MPBall))
err = Matrix (CN MPBall)
idM Matrix (CN MPBall)
-> Matrix (CN MPBall)
-> SubType (Matrix (CN MPBall)) (Matrix (CN MPBall))
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
- Matrix (CN MPBall)
y Matrix (CN MPBall)
-> Matrix (CN MPBall)
-> MulType (Matrix (CN MPBall)) (Matrix (CN MPBall))
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
* Matrix (CN MPBall)
m
nerr :: CN MPBall
nerr = Matrix (CN MPBall) -> CN MPBall
forall a.
(CanAddSameType a, CanSubSameType a, CanAbsSameType a,
HasIntegers a, CanMinMaxSameType a) =>
Matrix a -> a
inftyNorm SubType (Matrix (CN MPBall)) (Matrix (CN MPBall))
Matrix (CN MPBall)
err
ui :: CN MPBall
ui = MPBall -> CN MPBall
forall v. v -> CN v
cn (MPBall -> CN MPBall) -> MPBall -> CN MPBall
forall a b. (a -> b) -> a -> b
$ MPBall -> MPBall -> MPBall
forall i.
(IsInterval i, CanMinMaxSameType (IntervalEndpoint i)) =>
i -> i -> i
fromEndpointsAsIntervals (Integer -> MPBall
forall t. CanBeMPBall t => t -> MPBall
mpBall (Integer -> MPBall) -> Integer -> MPBall
forall a b. (a -> b) -> a -> b
$ -Integer
1) (Integer -> MPBall
forall t. CanBeMPBall t => t -> MPBall
mpBall Integer
1)
x0 :: Matrix (CN MPBall)
x0 = Integer
-> Integer
-> (Integer -> Integer -> CN MPBall)
-> Matrix (CN MPBall)
forall a.
Integer -> Integer -> (Integer -> Integer -> a) -> Matrix a
create (Matrix (CN MPBall) -> Integer
forall a. Matrix a -> Integer
width Matrix (CN MPBall)
m) (Matrix (CN MPBall) -> Integer
forall a. Matrix a -> Integer
width Matrix (CN MPBall)
m) (\Integer
_ Integer
j -> (CN MPBall
ui CN MPBall -> CN MPBall -> MulType (CN MPBall) (CN MPBall)
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
* (Vector (CN MPBall)
nye Vector (CN MPBall) -> Integer -> CN MPBall
forall a. Vector a -> Integer -> a
V.! Integer
j)) CN MPBall -> CN MPBall -> DivType (CN MPBall) (CN MPBall)
forall t1 t2. CanDiv t1 t2 => t1 -> t2 -> DivType t1 t2
/ (Integer
1 Integer -> CN MPBall -> SubType Integer (CN MPBall)
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
- CN MPBall
nerr))
aux :: Matrix (CN MPBall) -> Matrix (CN MPBall)
aux Matrix (CN MPBall)
x = let x' :: IntersectionType (Matrix (CN MPBall)) (Matrix (CN MPBall))
x' = Matrix (CN MPBall)
-> IntersectionType (Matrix (CN MPBall)) (Matrix (CN MPBall))
it Matrix (CN MPBall)
x in if Matrix (CN MPBall) -> Accuracy
forall a. HasAccuracy a => a -> Accuracy
getAccuracy IntersectionType (Matrix (CN MPBall)) (Matrix (CN MPBall))
Matrix (CN MPBall)
x' Accuracy -> Accuracy -> OrderCompareType Accuracy Accuracy
forall a b.
HasOrderAsymmetric a b =>
a -> b -> OrderCompareType a b
<= Matrix (CN MPBall) -> Accuracy
forall a. HasAccuracy a => a -> Accuracy
getAccuracy Matrix (CN MPBall)
x then IntersectionType (Matrix (CN MPBall)) (Matrix (CN MPBall))
Matrix (CN MPBall)
x' else Matrix (CN MPBall) -> Matrix (CN MPBall)
aux IntersectionType (Matrix (CN MPBall)) (Matrix (CN MPBall))
Matrix (CN MPBall)
x'
it :: Matrix (CN MPBall)
-> IntersectionType (Matrix (CN MPBall)) (Matrix (CN MPBall))
it Matrix (CN MPBall)
x = Matrix (CN MPBall)
-> Matrix (CN MPBall)
-> IntersectionType (Matrix (CN MPBall)) (Matrix (CN MPBall))
forall e1 e2.
CanIntersectAsymmetric e1 e2 =>
e1 -> e2 -> IntersectionType e1 e2
intersect Matrix (CN MPBall)
x (Matrix (CN MPBall)
-> IntersectionType (Matrix (CN MPBall)) (Matrix (CN MPBall)))
-> Matrix (CN MPBall)
-> IntersectionType (Matrix (CN MPBall)) (Matrix (CN MPBall))
forall a b. (a -> b) -> a -> b
$ MulType (Matrix (CN MPBall)) (Matrix (CN MPBall))
Matrix (CN MPBall)
yid Matrix (CN MPBall)
-> Matrix (CN MPBall)
-> AddType (Matrix (CN MPBall)) (Matrix (CN MPBall))
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ SubType (Matrix (CN MPBall)) (Matrix (CN MPBall))
Matrix (CN MPBall)
err Matrix (CN MPBall)
-> Matrix (CN MPBall)
-> MulType (Matrix (CN MPBall)) (Matrix (CN MPBall))
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
* Matrix (CN MPBall)
x
inverseGauss :: Matrix MPFloat -> Maybe (Matrix MPFloat)
inverseGauss :: Matrix MPFloat -> Maybe (Matrix MPFloat)
inverseGauss Matrix MPFloat
m =
Integer
-> Matrix MPFloat -> Matrix MPFloat -> Maybe (Matrix MPFloat)
aux Integer
0 Matrix MPFloat
m Matrix MPFloat
inv0
where
inv0 :: Matrix MPFloat
inv0 = Integer -> Integer -> Matrix MPFloat
forall a. HasIntegers a => Integer -> Integer -> Matrix a
identity (Matrix MPFloat -> Integer
forall a. Matrix a -> Integer
width Matrix MPFloat
m) (Matrix MPFloat -> Integer
forall a. Matrix a -> Integer
width Matrix MPFloat
m) :: Matrix MPFloat
aux :: Integer
-> Matrix MPFloat -> Matrix MPFloat -> Maybe (Matrix MPFloat)
aux Integer
j Matrix MPFloat
n Matrix MPFloat
inv =
if Integer
j Integer -> Integer -> EqCompareType Integer Integer
forall a b. HasEqAsymmetric a b => a -> b -> EqCompareType a b
== Matrix MPFloat -> Integer
forall a. Matrix a -> Integer
width Matrix MPFloat
m then
Matrix MPFloat -> Maybe (Matrix MPFloat)
forall a. a -> Maybe a
Just Matrix MPFloat
inv
else
let
n' :: Matrix MPFloat
n' = Matrix MPFloat -> Integer -> Matrix MPFloat -> Matrix MPFloat
pivot Matrix MPFloat
n Integer
j Matrix MPFloat
n
inv' :: Matrix MPFloat
inv' = Matrix MPFloat -> Integer -> Matrix MPFloat -> Matrix MPFloat
pivot Matrix MPFloat
n Integer
j Matrix MPFloat
inv
a :: MPFloat
a = Matrix MPFloat -> Integer -> Integer -> MPFloat
forall a. Matrix a -> Integer -> Integer -> a
get Matrix MPFloat
n' Integer
j Integer
j
ia :: MPFloat
ia = ((Integer -> MPFloat
forall t. CanBeMPFloat t => t -> MPFloat
mpFloat Integer
1)MPFloat -> MPFloat -> MPFloat
/.MPFloat
a)
n'' :: Matrix MPFloat
n'' = Matrix MPFloat -> Integer -> MPFloat -> Matrix MPFloat
multiplyRow Matrix MPFloat
n' Integer
j MPFloat
ia
inv'' :: Matrix MPFloat
inv'' = Matrix MPFloat -> Integer -> MPFloat -> Matrix MPFloat
multiplyRow Matrix MPFloat
inv' Integer
j MPFloat
ia
elim :: Integer -> Matrix MPFloat -> IfThenElseType Bool (Matrix MPFloat)
elim Integer
i Matrix MPFloat
l =
if Integer
i Integer -> Integer -> EqCompareType Integer Integer
forall a b. HasEqAsymmetric a b => a -> b -> EqCompareType a b
/= Integer
j then
let
b :: MPFloat
b = Matrix MPFloat -> Integer -> Integer -> MPFloat
forall a. Matrix a -> Integer -> Integer -> a
get Matrix MPFloat
l Integer
i Integer
j
in
Matrix MPFloat -> MPFloat -> Integer -> Integer -> Matrix MPFloat
addRows Matrix MPFloat
l (-MPFloat
b) Integer
j Integer
i
else
Matrix MPFloat
l
n''' :: Matrix MPFloat
n''' =
(Integer -> Matrix MPFloat -> Matrix MPFloat)
-> Matrix MPFloat -> [Integer] -> Matrix MPFloat
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr Integer -> Matrix MPFloat -> IfThenElseType Bool (Matrix MPFloat)
Integer -> Matrix MPFloat -> Matrix MPFloat
elim Matrix MPFloat
n'' [Integer
0 .. (Matrix MPFloat -> Integer
forall a. Matrix a -> Integer
height Matrix MPFloat
m) Integer -> Integer -> SubType Integer Integer
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
- Integer
1]
inv''' :: Matrix MPFloat
inv''' =
(Integer -> Matrix MPFloat -> Matrix MPFloat)
-> Matrix MPFloat -> [Integer] -> Matrix MPFloat
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr Integer -> Matrix MPFloat -> IfThenElseType Bool (Matrix MPFloat)
Integer -> Matrix MPFloat -> Matrix MPFloat
elim Matrix MPFloat
inv'' [Integer
0 .. (Matrix MPFloat -> Integer
forall a. Matrix a -> Integer
height Matrix MPFloat
m) Integer -> Integer -> SubType Integer Integer
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
- Integer
1]
in
if MPFloat
a MPFloat -> MPFloat -> EqCompareType MPFloat MPFloat
forall a b. HasEqAsymmetric a b => a -> b -> EqCompareType a b
/= Integer -> MPFloat
forall t. CanBeMPFloat t => t -> MPFloat
mpFloat Integer
0 then
Integer
-> Matrix MPFloat -> Matrix MPFloat -> Maybe (Matrix MPFloat)
aux (Integer
j Integer -> Integer -> AddType Integer Integer
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Integer
1) Matrix MPFloat
n''' Matrix MPFloat
inv'''
else
Maybe (Matrix MPFloat)
forall a. Maybe a
Nothing
largestRowIndex :: Matrix MPFloat -> Integer -> Integer
largestRowIndex :: Matrix MPFloat -> Integer -> Integer
largestRowIndex Matrix MPFloat
m Integer
j =
Integer -> Integer
forall t. CanBeInteger t => t -> Integer
integer (Integer -> Integer) -> Integer -> Integer
forall a b. (a -> b) -> a -> b
$ Integer -> Integer -> MPFloat -> Integer
aux (Integer
j Integer -> Integer -> AddType Integer Integer
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Integer
1) Integer
j (Matrix MPFloat -> Integer -> Integer -> MPFloat
forall a. Matrix a -> Integer -> Integer -> a
get Matrix MPFloat
m Integer
j Integer
j)
where
h :: Integer
h = Matrix MPFloat -> Integer
forall a. Matrix a -> Integer
height Matrix MPFloat
m
aux :: Integer -> Integer -> MPFloat -> Integer
aux Integer
k Integer
i MPFloat
x =
if Integer
k Integer -> Integer -> EqCompareType Integer Integer
forall a b. HasEqAsymmetric a b => a -> b -> EqCompareType a b
== Integer
h then
Integer
i
else
let
x' :: MPFloat
x' = Matrix MPFloat -> Integer -> Integer -> MPFloat
forall a. Matrix a -> Integer -> Integer -> a
get Matrix MPFloat
m Integer
i Integer
j
in
if MPFloat
x' MPFloat -> MPFloat -> OrderCompareType MPFloat MPFloat
forall a b.
HasOrderAsymmetric a b =>
a -> b -> OrderCompareType a b
> MPFloat
x then
Integer -> Integer -> MPFloat -> Integer
aux (Integer
k Integer -> Integer -> AddType Integer Integer
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Integer
1) Integer
k MPFloat
x'
else
Integer -> Integer -> MPFloat -> Integer
aux (Integer
k Integer -> Integer -> AddType Integer Integer
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Integer
1) Integer
i MPFloat
x
pivot :: Matrix MPFloat -> Integer -> Matrix MPFloat -> Matrix MPFloat
pivot :: Matrix MPFloat -> Integer -> Matrix MPFloat -> Matrix MPFloat
pivot Matrix MPFloat
m Integer
j Matrix MPFloat
n =
Matrix MPFloat -> Integer -> Integer -> Matrix MPFloat
swapRows Matrix MPFloat
n Integer
j Integer
i
where
i :: Integer
i = Matrix MPFloat -> Integer -> Integer
largestRowIndex Matrix MPFloat
m Integer
j
multiplyRow :: Matrix MPFloat -> Integer -> MPFloat -> Matrix MPFloat
multiplyRow :: Matrix MPFloat -> Integer -> MPFloat -> Matrix MPFloat
multiplyRow Matrix MPFloat
m Integer
i MPFloat
a =
(Integer -> Integer -> MPFloat -> MPFloat)
-> Matrix MPFloat -> Matrix MPFloat
forall a. (Integer -> Integer -> a -> a) -> Matrix a -> Matrix a
imap Integer -> Integer -> MPFloat -> MPFloat
Integer -> Integer -> MPFloat -> IfThenElseType Bool MPFloat
h Matrix MPFloat
m
where
h :: Integer -> Integer -> MPFloat -> IfThenElseType Bool MPFloat
h Integer
k Integer
_ MPFloat
x =
if Integer
k Integer -> Integer -> EqCompareType Integer Integer
forall a b. HasEqAsymmetric a b => a -> b -> EqCompareType a b
== Integer
i then
MPFloat
a MPFloat -> MPFloat -> MPFloat
*. MPFloat
x
else
MPFloat
x
addRows :: Matrix MPFloat -> MPFloat -> Integer -> Integer -> Matrix MPFloat
addRows :: Matrix MPFloat -> MPFloat -> Integer -> Integer -> Matrix MPFloat
addRows Matrix MPFloat
m MPFloat
a Integer
i0 Integer
i1 =
(Integer -> Integer -> MPFloat -> MPFloat)
-> Matrix MPFloat -> Matrix MPFloat
forall a. (Integer -> Integer -> a -> a) -> Matrix a -> Matrix a
imap Integer -> Integer -> MPFloat -> MPFloat
Integer -> Integer -> MPFloat -> IfThenElseType Bool MPFloat
h Matrix MPFloat
m
where
h :: Integer -> Integer -> MPFloat -> IfThenElseType Bool MPFloat
h Integer
k Integer
l MPFloat
x =
if Integer
k Integer -> Integer -> EqCompareType Integer Integer
forall a b. HasEqAsymmetric a b => a -> b -> EqCompareType a b
== Integer
i1 then
MPFloat
x MPFloat -> MPFloat -> MPFloat
+. MPFloat
a MPFloat -> MPFloat -> MPFloat
*. Matrix MPFloat -> Integer -> Integer -> MPFloat
forall a. Matrix a -> Integer -> Integer -> a
get Matrix MPFloat
m Integer
i0 Integer
l
else
MPFloat
x
swapRows :: Matrix MPFloat -> Integer -> Integer -> Matrix MPFloat
swapRows :: Matrix MPFloat -> Integer -> Integer -> Matrix MPFloat
swapRows Matrix MPFloat
m Integer
i0 Integer
i1 =
(Integer -> Integer -> MPFloat -> MPFloat)
-> Matrix MPFloat -> Matrix MPFloat
forall a. (Integer -> Integer -> a -> a) -> Matrix a -> Matrix a
imap Integer -> Integer -> MPFloat -> MPFloat
Integer -> Integer -> MPFloat -> IfThenElseType Bool MPFloat
h Matrix MPFloat
m
where
h :: Integer -> Integer -> MPFloat -> IfThenElseType Bool MPFloat
h Integer
k Integer
l MPFloat
x =
if Integer
k Integer -> Integer -> EqCompareType Integer Integer
forall a b. HasEqAsymmetric a b => a -> b -> EqCompareType a b
== Integer
i0 then
Matrix MPFloat -> Integer -> Integer -> MPFloat
forall a. Matrix a -> Integer -> Integer -> a
get Matrix MPFloat
m Integer
i1 Integer
l
else if Integer
k Integer -> Integer -> EqCompareType Integer Integer
forall a b. HasEqAsymmetric a b => a -> b -> EqCompareType a b
== Integer
i1 then
Matrix MPFloat -> Integer -> Integer -> MPFloat
forall a. Matrix a -> Integer -> Integer -> a
get Matrix MPFloat
m Integer
i0 Integer
l
else
MPFloat
x