{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE FlexibleContexts #-}

-----------------------------------------------------------------------------
{-|
Module      : Math.Tensor.LinearAlgebra.Matrix
Description : Gaussian elimination subroutines.
Copyright   : (c) Nils Alex, 2020
License     : MIT
Maintainer  : nils.alex@fau.de

Gaussian elimination subroutines.
-}
-----------------------------------------------------------------------------
module Math.Tensor.LinearAlgebra.Matrix
  ( gaussianST
  , gaussianFFST
  , gaussian
  , gaussianFF
  , rrefST
  , rref
  , independentColumns
  , independentColumnsFF
  , independentColumnsRREF
  , independentColumnsVerifiedFF
  , independentColumnsMat
  , independentColumnsMatFF
  , independentColumnsMatRREF
  , pivotsU
  , pivotsUFF
  , findPivotMax
  , findPivotMaxFF
  , findRowPivot
  , isref
  , isrref
  , isrref'
  , verify
  ) where

import Numeric.LinearAlgebra
  ( Matrix
  , Vector
  , Container
  , Extractor (All, Take, Drop)
  , Z
  , toLists
  , rows
  , cols
  , find
  , (¿)
  , (??)
  , (><)
  , (===)
  , rank
  , fromZ
  )
import Numeric.LinearAlgebra.Devel
  ( STMatrix
  , RowOper (AXPY, SCAL, SWAP)
  , ColRange (FromCol)
  , RowRange (Row)
  , freezeMatrix
  , thawMatrix
  , modifyMatrix
  , readMatrix
  , rowOper
  )

import Data.List (maximumBy)

import Control.Monad (foldM)
import Control.Monad.ST
  ( ST
  , runST
  )

-- | Returns the pivot columns of an upper triangular matrix.
--
-- @
-- &#x3BB; let mat = (3 >< 4) [1, 0, 2, -3, 0, 0, 1, 0, 0, 0, 0, 0]
-- &#x3BB; mat
-- (3><4)
--  [ 1.0, 0.0, 2.0, -3.0
--  , 0.0, 0.0, 1.0,  0.0
--  , 0.0, 0.0, 0.0,  0.0 ]
-- &#x3BB; pivotsU mat
-- [0,2]
-- @
--
pivotsU :: Matrix Double -> [Int]
pivotsU :: Matrix Double -> [Int]
pivotsU Matrix Double
mat = (Int, Int) -> [Int]
go (Int
0,Int
0)
  where
    go :: (Int, Int) -> [Int]
go (Int
i,Int
j)
      = case Matrix Double -> Double -> (Int, Int) -> Maybe (Int, Int)
findPivot Matrix Double
mat Double
e (Int
i,Int
j) of
          Maybe (Int, Int)
Nothing       -> []
          Just (Int
i', Int
j') -> Int
j' Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: (Int, Int) -> [Int]
go (Int
i'Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1, Int
j'Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
    maxAbs :: Double
maxAbs = [Double] -> Double
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum ([Double] -> Double) -> [Double] -> Double
forall a b. (a -> b) -> a -> b
$ ([Double] -> Double) -> [[Double]] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> Double
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum ([Double] -> Double)
-> ([Double] -> [Double]) -> [Double] -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
abs) ([[Double]] -> [Double]) -> [[Double]] -> [Double]
forall a b. (a -> b) -> a -> b
$ Matrix Double -> [[Double]]
forall t. Element t => Matrix t -> [[t]]
toLists Matrix Double
mat
    e :: Double
e = Double
eps Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
maxAbs

pivotsUFF :: Matrix Z -> [Int]
pivotsUFF :: Matrix Z -> [Int]
pivotsUFF Matrix Z
mat = (Int, Int) -> [Int]
go (Int
0,Int
0)
  where
    go :: (Int, Int) -> [Int]
go (Int
i,Int
j)
      = case Matrix Z -> (Int, Int) -> Maybe (Int, Int)
findPivotFF Matrix Z
mat (Int
i,Int
j) of
          Maybe (Int, Int)
Nothing       -> []
          Just (Int
i', Int
j') -> Int
j' Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: (Int, Int) -> [Int]
go (Int
i'Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1, Int
j'Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)

eps :: Double
eps :: Double
eps = Double
1e-12

-- find next pivot in upper triangular matrix
findPivotFF :: Matrix Z -> (Int, Int) -> Maybe (Int, Int)
findPivotFF :: Matrix Z -> (Int, Int) -> Maybe (Int, Int)
findPivotFF Matrix Z
mat (Int
i, Int
j)
    | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j = Maybe (Int, Int)
forall a. Maybe a
Nothing
    | Int
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i = Maybe (Int, Int)
forall a. Maybe a
Nothing
    | Bool
otherwise = case [(Int, Int)]
nonZeros of
                    []           -> if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1
                                    then Maybe (Int, Int)
forall a. Maybe a
Nothing
                                    else Matrix Z -> (Int, Int) -> Maybe (Int, Int)
findPivotFF Matrix Z
mat (Int
i, Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
                    (Int
pi_, Int
pj):[(Int, Int)]
_  -> (Int, Int) -> Maybe (Int, Int)
forall a. a -> Maybe a
Just (Int
pi_, Int
pjInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
j)
    where
        m :: Int
m = Matrix Z -> Int
forall t. Matrix t -> Int
rows Matrix Z
mat
        n :: Int
n = Matrix Z -> Int
forall t. Matrix t -> Int
cols Matrix Z
mat
        col :: Matrix Z
col = Matrix Z
mat Matrix Z -> [Int] -> Matrix Z
forall t. Element t => Matrix t -> [Int] -> Matrix t
¿ [Int
j]
        nonZeros :: [(Int, Int)]
nonZeros = ((Int, Int) -> Bool) -> [(Int, Int)] -> [(Int, Int)]
forall a. (a -> Bool) -> [a] -> [a]
filter (\(Int
i', Int
_) -> Int
i' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
i) ([(Int, Int)] -> [(Int, Int)]) -> [(Int, Int)] -> [(Int, Int)]
forall a b. (a -> b) -> a -> b
$ (Z -> Bool) -> Matrix Z -> [IndexOf Matrix]
forall (c :: * -> *) e.
Container c e =>
(e -> Bool) -> c e -> [IndexOf c]
find (Z -> Z -> Bool
forall a. Eq a => a -> a -> Bool
/= Z
0) Matrix Z
col

findPivot :: Matrix Double -> Double -> (Int, Int) -> Maybe (Int, Int)
findPivot :: Matrix Double -> Double -> (Int, Int) -> Maybe (Int, Int)
findPivot Matrix Double
mat Double
e (Int
i, Int
j)
    | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j = Maybe (Int, Int)
forall a. Maybe a
Nothing
    | Int
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i = Maybe (Int, Int)
forall a. Maybe a
Nothing
    | Bool
otherwise = case [(Int, Int)]
nonZeros of
                    []           -> if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1
                                    then Maybe (Int, Int)
forall a. Maybe a
Nothing
                                    else Matrix Double -> Double -> (Int, Int) -> Maybe (Int, Int)
findPivot Matrix Double
mat Double
e (Int
i, Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
                    (Int
pi_, Int
pj):[(Int, Int)]
_  -> (Int, Int) -> Maybe (Int, Int)
forall a. a -> Maybe a
Just (Int
pi_, Int
pjInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
j)
    where
        m :: Int
m = Matrix Double -> Int
forall t. Matrix t -> Int
rows Matrix Double
mat
        n :: Int
n = Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
mat
        col :: Matrix Double
col = Matrix Double
mat Matrix Double -> [Int] -> Matrix Double
forall t. Element t => Matrix t -> [Int] -> Matrix t
¿ [Int
j]
        nonZeros :: [(Int, Int)]
nonZeros = ((Int, Int) -> Bool) -> [(Int, Int)] -> [(Int, Int)]
forall a. (a -> Bool) -> [a] -> [a]
filter (\(Int
i', Int
_) -> Int
i' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
i) ([(Int, Int)] -> [(Int, Int)]) -> [(Int, Int)] -> [(Int, Int)]
forall a b. (a -> b) -> a -> b
$ (Double -> Bool) -> Matrix Double -> [IndexOf Matrix]
forall (c :: * -> *) e.
Container c e =>
(e -> Bool) -> c e -> [IndexOf c]
find ((Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
e) (Double -> Bool) -> (Double -> Double) -> Double -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
forall a. Num a => a -> a
abs) Matrix Double
col

-- | Find pivot element below position (i, j) with greatest absolute value in the ST monad.
findPivotMaxFF :: Int -> Int -> Int -> Int -> STMatrix s Z -> ST s (Maybe (Int, Int))
findPivotMaxFF :: Int -> Int -> Int -> Int -> STMatrix s Z -> ST s (Maybe (Int, Int))
findPivotMaxFF Int
m Int
n Int
i Int
j STMatrix s Z
mat
    | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j = Maybe (Int, Int) -> ST s (Maybe (Int, Int))
forall (m :: * -> *) a. Monad m => a -> m a
return Maybe (Int, Int)
forall a. Maybe a
Nothing
    | Int
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i = Maybe (Int, Int) -> ST s (Maybe (Int, Int))
forall (m :: * -> *) a. Monad m => a -> m a
return Maybe (Int, Int)
forall a. Maybe a
Nothing
    | Bool
otherwise =
        do
          [(Int, Z)]
col      <- (Int -> ST s (Int, Z)) -> [Int] -> ST s [(Int, Z)]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (\Int
i' -> do
                                    Z
x <- STMatrix s Z -> Int -> Int -> ST s Z
forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
readMatrix STMatrix s Z
mat Int
i' Int
j
                                    (Int, Z) -> ST s (Int, Z)
forall (m :: * -> *) a. Monad m => a -> m a
return (Int
i', Z
x))
                      [Int
i..Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]
          let nonZeros :: [(Int, Z)]
nonZeros = ((Int, Z) -> Bool) -> [(Int, Z)] -> [(Int, Z)]
forall a. (a -> Bool) -> [a] -> [a]
filter ((Z -> Z -> Bool
forall a. Eq a => a -> a -> Bool
/= Z
0) (Z -> Bool) -> ((Int, Z) -> Z) -> (Int, Z) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int, Z) -> Z
forall a b. (a, b) -> b
snd) [(Int, Z)]
col
          case [(Int, Z)]
nonZeros of
            []         -> if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1
                          then Maybe (Int, Int) -> ST s (Maybe (Int, Int))
forall (m :: * -> *) a. Monad m => a -> m a
return Maybe (Int, Int)
forall a. Maybe a
Nothing
                          else Int -> Int -> Int -> Int -> STMatrix s Z -> ST s (Maybe (Int, Int))
forall s.
Int -> Int -> Int -> Int -> STMatrix s Z -> ST s (Maybe (Int, Int))
findPivotMaxFF Int
m Int
n Int
i (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) STMatrix s Z
mat
            (Int
pi_,Z
_):[(Int, Z)]
_  -> Maybe (Int, Int) -> ST s (Maybe (Int, Int))
forall (m :: * -> *) a. Monad m => a -> m a
return (Maybe (Int, Int) -> ST s (Maybe (Int, Int)))
-> Maybe (Int, Int) -> ST s (Maybe (Int, Int))
forall a b. (a -> b) -> a -> b
$ (Int, Int) -> Maybe (Int, Int)
forall a. a -> Maybe a
Just (Int
pi_, Int
j)

findPivotMax :: Int -> Int -> Int -> Int -> STMatrix s Double -> ST s (Maybe (Int, Int))
findPivotMax :: Int
-> Int
-> Int
-> Int
-> STMatrix s Double
-> ST s (Maybe (Int, Int))
findPivotMax Int
m Int
n Int
i Int
j STMatrix s Double
mat
    | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j = Maybe (Int, Int) -> ST s (Maybe (Int, Int))
forall (m :: * -> *) a. Monad m => a -> m a
return Maybe (Int, Int)
forall a. Maybe a
Nothing
    | Int
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i = Maybe (Int, Int) -> ST s (Maybe (Int, Int))
forall (m :: * -> *) a. Monad m => a -> m a
return Maybe (Int, Int)
forall a. Maybe a
Nothing
    | Bool
otherwise =
        do
          [(Int, Double)]
col      <- (Int -> ST s (Int, Double)) -> [Int] -> ST s [(Int, Double)]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (\Int
i' -> do
                                    Double
x <- STMatrix s Double -> Int -> Int -> ST s Double
forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
readMatrix STMatrix s Double
mat Int
i' Int
j
                                    (Int, Double) -> ST s (Int, Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Int
i', Double -> Double
forall a. Num a => a -> a
abs Double
x))
                      [Int
i..Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]
          let nonZeros :: [(Int, Double)]
nonZeros = ((Int, Double) -> Bool) -> [(Int, Double)] -> [(Int, Double)]
forall a. (a -> Bool) -> [a] -> [a]
filter ((Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
eps) (Double -> Bool)
-> ((Int, Double) -> Double) -> (Int, Double) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
forall a. Num a => a -> a
abs (Double -> Double)
-> ((Int, Double) -> Double) -> (Int, Double) -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int, Double) -> Double
forall a b. (a, b) -> b
snd) [(Int, Double)]
col
          let (Int
pi_, Double
_) = ((Int, Double) -> (Int, Double) -> Ordering)
-> [(Int, Double)] -> (Int, Double)
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
maximumBy (\(Int
_, Double
x) (Int
_, Double
y) -> Double
x Double -> Double -> Ordering
forall a. Ord a => a -> a -> Ordering
`compare` Double
y) [(Int, Double)]
nonZeros
          case [(Int, Double)]
nonZeros of
            [] -> if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1
                  then Maybe (Int, Int) -> ST s (Maybe (Int, Int))
forall (m :: * -> *) a. Monad m => a -> m a
return Maybe (Int, Int)
forall a. Maybe a
Nothing
                  else Int
-> Int
-> Int
-> Int
-> STMatrix s Double
-> ST s (Maybe (Int, Int))
forall s.
Int
-> Int
-> Int
-> Int
-> STMatrix s Double
-> ST s (Maybe (Int, Int))
findPivotMax Int
m Int
n Int
i (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) STMatrix s Double
mat
            [(Int, Double)]
_  -> Maybe (Int, Int) -> ST s (Maybe (Int, Int))
forall (m :: * -> *) a. Monad m => a -> m a
return (Maybe (Int, Int) -> ST s (Maybe (Int, Int)))
-> Maybe (Int, Int) -> ST s (Maybe (Int, Int))
forall a b. (a -> b) -> a -> b
$ (Int, Int) -> Maybe (Int, Int)
forall a. a -> Maybe a
Just (Int
pi_, Int
j)

findRowPivot :: Int -> Int -> Int -> Int -> STMatrix s Z -> ST s (Maybe Int)
findRowPivot :: Int -> Int -> Int -> Int -> STMatrix s Z -> ST s (Maybe Int)
findRowPivot Int
m Int
n Int
i Int
j STMatrix s Z
mat
    | Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
n       = [Char] -> ST s (Maybe Int)
forall a. HasCallStack => [Char] -> a
error [Char]
"out of bounds" -- return Nothing
    | Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
m Int
n = [Char] -> ST s (Maybe Int)
forall a. HasCallStack => [Char] -> a
error [Char]
"out of bounds" -- return Nothing
    | Bool
otherwise =
        do
         [(Int, Z)]
row <- (Int -> ST s (Int, Z)) -> [Int] -> ST s [(Int, Z)]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (\Int
j' -> do
                              Z
x <- STMatrix s Z -> Int -> Int -> ST s Z
forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
readMatrix STMatrix s Z
mat Int
i Int
j'
                              (Int, Z) -> ST s (Int, Z)
forall (m :: * -> *) a. Monad m => a -> m a
return (Int
j', Z
x))
                [Int
0 .. Int
j]
         let nonZeros :: [(Int, Z)]
nonZeros = ((Int, Z) -> Bool) -> [(Int, Z)] -> [(Int, Z)]
forall a. (a -> Bool) -> [a] -> [a]
filter ((Z -> Z -> Bool
forall a. Eq a => a -> a -> Bool
/=Z
0) (Z -> Bool) -> ((Int, Z) -> Z) -> (Int, Z) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int, Z) -> Z
forall a b. (a, b) -> b
snd) [(Int, Z)]
row
         case [(Int, Z)]
nonZeros of
           []        -> Maybe Int -> ST s (Maybe Int)
forall (m :: * -> *) a. Monad m => a -> m a
return Maybe Int
forall a. Maybe a
Nothing
           (Int
pj, Z
_):[(Int, Z)]
_ -> Maybe Int -> ST s (Maybe Int)
forall (m :: * -> *) a. Monad m => a -> m a
return (Maybe Int -> ST s (Maybe Int)) -> Maybe Int -> ST s (Maybe Int)
forall a b. (a -> b) -> a -> b
$ Int -> Maybe Int
forall a. a -> Maybe a
Just Int
pj

backwardFF' :: Int -> Int -> Int -> Int -> STMatrix s Z -> ST s ()
backwardFF' :: Int -> Int -> Int -> Int -> STMatrix s Z -> ST s ()
backwardFF' Int
m Int
n Int
i Int
j STMatrix s Z
mat
      | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = () -> ST s ()
forall (m :: * -> *) a. Monad m => a -> m a
return ()
      | Bool
otherwise = do
    Maybe Int
iPivot' <- Int -> Int -> Int -> Int -> STMatrix s Z -> ST s (Maybe Int)
forall s.
Int -> Int -> Int -> Int -> STMatrix s Z -> ST s (Maybe Int)
findRowPivot Int
m Int
n Int
i Int
j STMatrix s Z
mat
    case Maybe Int
iPivot' of
        Maybe Int
Nothing -> Int -> Int -> Int -> Int -> STMatrix s Z -> ST s ()
forall s. Int -> Int -> Int -> Int -> STMatrix s Z -> ST s ()
backwardFF' Int
m Int
n (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Int
j STMatrix s Z
mat
        Just Int
pj -> do
                    Z
pv <- STMatrix s Z -> Int -> Int -> ST s Z
forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
readMatrix STMatrix s Z
mat Int
i Int
pj
                    (Int -> ST s ()) -> [Int] -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (Z -> Int -> Int -> ST s ()
reduce Z
pv Int
pj) [Int
0 .. Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]
                    Int -> Int -> Int -> Int -> STMatrix s Z -> ST s ()
forall s. Int -> Int -> Int -> Int -> STMatrix s Z -> ST s ()
backwardFF' Int
m Int
n (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) (Int
pjInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) STMatrix s Z
mat
  where
    reduce :: Z -> Int -> Int -> ST s ()
reduce Z
pv Int
pj Int
r = do
                      Just Int
pr <- Int -> Int -> Int -> Int -> STMatrix s Z -> ST s (Maybe Int)
forall s.
Int -> Int -> Int -> Int -> STMatrix s Z -> ST s (Maybe Int)
findRowPivot Int
m Int
n Int
r Int
pj STMatrix s Z
mat
                      -- prv <- readMatrix mat r pr
                      Z
pjv <- STMatrix s Z -> Int -> Int -> ST s Z
forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
readMatrix STMatrix s Z
mat Int
r Int
pj
                      if Z
pjv Z -> Z -> Bool
forall a. Eq a => a -> a -> Bool
== Z
0
                        then () -> ST s ()
forall (m :: * -> *) a. Monad m => a -> m a
return ()
                        else
                         let op1 :: RowOper Z
op1 = Z -> RowRange -> ColRange -> RowOper Z
forall t. t -> RowRange -> ColRange -> RowOper t
SCAL Z
pv (Int -> RowRange
Row Int
r) (Int -> ColRange
FromCol Int
pr)
                             op2 :: RowOper Z
op2 = Z -> Int -> Int -> ColRange -> RowOper Z
forall t. t -> Int -> Int -> ColRange -> RowOper t
AXPY (-Z
pjv) Int
i Int
r (Int -> ColRange
FromCol Int
pj)
                         in do
                             RowOper Z -> STMatrix s Z -> ST s ()
forall t s.
(Num t, Element t) =>
RowOper t -> STMatrix s t -> ST s ()
rowOper RowOper Z
op1 STMatrix s Z
mat
                             RowOper Z -> STMatrix s Z -> ST s ()
forall t s.
(Num t, Element t) =>
RowOper t -> STMatrix s t -> ST s ()
rowOper RowOper Z
op2 STMatrix s Z
mat
                             Z
g <- (Z -> Int -> ST s Z) -> Z -> [Int] -> ST s Z
forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM (\Z
acc Int
c -> Z -> Z -> Z
forall a. Integral a => a -> a -> a
gcd Z
acc (Z -> Z) -> ST s Z -> ST s Z
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> STMatrix s Z -> Int -> Int -> ST s Z
forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
readMatrix STMatrix s Z
mat Int
r Int
c) Z
0 [Int
pr .. Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]
                             if Z
g Z -> Z -> Bool
forall a. Eq a => a -> a -> Bool
== Z
0
                               then () -> ST s ()
forall (m :: * -> *) a. Monad m => a -> m a
return()
                               else (Int -> ST s ()) -> [Int] -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
c -> STMatrix s Z -> Int -> Int -> (Z -> Z) -> ST s ()
forall t s.
Storable t =>
STMatrix s t -> Int -> Int -> (t -> t) -> ST s ()
modifyMatrix STMatrix s Z
mat Int
r Int
c (Z -> Z -> Z
forall a. Integral a => a -> a -> a
`quot` Z
g)) [Int
pr .. Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]

gaussianFF' :: Int -> Int -> Int -> Int -> STMatrix s Z -> ST s ()
gaussianFF' :: Int -> Int -> Int -> Int -> STMatrix s Z -> ST s ()
gaussianFF' Int
m Int
n Int
i Int
j STMatrix s Z
mat = do
    Maybe (Int, Int)
iPivot' <- Int -> Int -> Int -> Int -> STMatrix s Z -> ST s (Maybe (Int, Int))
forall s.
Int -> Int -> Int -> Int -> STMatrix s Z -> ST s (Maybe (Int, Int))
findPivotMaxFF Int
m Int
n Int
i Int
j STMatrix s Z
mat
    case Maybe (Int, Int)
iPivot' of
        Maybe (Int, Int)
Nothing     -> () -> ST s ()
forall (m :: * -> *) a. Monad m => a -> m a
return ()
        Just (Int
r, Int
p) -> do
                          RowOper Z -> STMatrix s Z -> ST s ()
forall t s.
(Num t, Element t) =>
RowOper t -> STMatrix s t -> ST s ()
rowOper (Int -> Int -> ColRange -> RowOper Z
forall t. Int -> Int -> ColRange -> RowOper t
SWAP Int
i Int
r (Int -> ColRange
FromCol Int
j)) STMatrix s Z
mat
                          Z
pv <- STMatrix s Z -> Int -> Int -> ST s Z
forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
readMatrix STMatrix s Z
mat Int
i Int
p
                          (Int -> ST s ()) -> [Int] -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (Z -> Int -> Int -> ST s ()
reduce Z
pv Int
p) [Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 .. Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]
                          Int -> Int -> Int -> Int -> STMatrix s Z -> ST s ()
forall s. Int -> Int -> Int -> Int -> STMatrix s Z -> ST s ()
gaussianFF' Int
m Int
n (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) (Int
pInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) STMatrix s Z
mat
  where
    reduce :: Z -> Int -> Int -> ST s ()
reduce Z
pv Int
p Int
r = do
                      Z
rv <- STMatrix s Z -> Int -> Int -> ST s Z
forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
readMatrix STMatrix s Z
mat Int
r Int
p
                      if Z
rv Z -> Z -> Bool
forall a. Eq a => a -> a -> Bool
== Z
0
                        then () -> ST s ()
forall (m :: * -> *) a. Monad m => a -> m a
return ()
                        else
                         let op1 :: RowOper Z
op1 = Z -> RowRange -> ColRange -> RowOper Z
forall t. t -> RowRange -> ColRange -> RowOper t
SCAL Z
pv (Int -> RowRange
Row Int
r) (Int -> ColRange
FromCol Int
p)
                             op2 :: RowOper Z
op2 = Z -> Int -> Int -> ColRange -> RowOper Z
forall t. t -> Int -> Int -> ColRange -> RowOper t
AXPY (-Z
rv) Int
i Int
r (Int -> ColRange
FromCol Int
p)
                         in do
                             RowOper Z -> STMatrix s Z -> ST s ()
forall t s.
(Num t, Element t) =>
RowOper t -> STMatrix s t -> ST s ()
rowOper RowOper Z
op1 STMatrix s Z
mat
                             RowOper Z -> STMatrix s Z -> ST s ()
forall t s.
(Num t, Element t) =>
RowOper t -> STMatrix s t -> ST s ()
rowOper RowOper Z
op2 STMatrix s Z
mat
                             Z
g <- (Z -> Int -> ST s Z) -> Z -> [Int] -> ST s Z
forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM (\Z
acc Int
c -> Z -> Z -> Z
forall a. Integral a => a -> a -> a
gcd Z
acc (Z -> Z) -> ST s Z -> ST s Z
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> STMatrix s Z -> Int -> Int -> ST s Z
forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
readMatrix STMatrix s Z
mat Int
r Int
c) Z
0 [Int
p .. Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]
                             if Z
g Z -> Z -> Bool
forall a. Eq a => a -> a -> Bool
== Z
0
                               then () -> ST s ()
forall (m :: * -> *) a. Monad m => a -> m a
return()
                               else (Int -> ST s ()) -> [Int] -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
c -> STMatrix s Z -> Int -> Int -> (Z -> Z) -> ST s ()
forall t s.
Storable t =>
STMatrix s t -> Int -> Int -> (t -> t) -> ST s ()
modifyMatrix STMatrix s Z
mat Int
r Int
c (Z -> Z -> Z
forall a. Integral a => a -> a -> a
`quot` Z
g)) [Int
p .. Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]

-- gaussian elimination of sub matrix below position (i, j)
gaussian' :: Int -> Int -> Int -> Int -> STMatrix s Double -> ST s ()
gaussian' :: Int -> Int -> Int -> Int -> STMatrix s Double -> ST s ()
gaussian' Int
m Int
n Int
i Int
j STMatrix s Double
mat = do
    Maybe (Int, Int)
iPivot' <- Int
-> Int
-> Int
-> Int
-> STMatrix s Double
-> ST s (Maybe (Int, Int))
forall s.
Int
-> Int
-> Int
-> Int
-> STMatrix s Double
-> ST s (Maybe (Int, Int))
findPivotMax Int
m Int
n Int
i Int
j STMatrix s Double
mat
    case Maybe (Int, Int)
iPivot' of
        Maybe (Int, Int)
Nothing     -> () -> ST s ()
forall (m :: * -> *) a. Monad m => a -> m a
return ()
        Just (Int
r, Int
p) -> do
                          RowOper Double -> STMatrix s Double -> ST s ()
forall t s.
(Num t, Element t) =>
RowOper t -> STMatrix s t -> ST s ()
rowOper (Int -> Int -> ColRange -> RowOper Double
forall t. Int -> Int -> ColRange -> RowOper t
SWAP Int
i Int
r (Int -> ColRange
FromCol Int
j)) STMatrix s Double
mat
                          Double
pv <- STMatrix s Double -> Int -> Int -> ST s Double
forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
readMatrix STMatrix s Double
mat Int
i Int
p
                          (Int -> ST s ()) -> [Int] -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (Double -> Int -> Int -> ST s ()
reduce Double
pv Int
p) [Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 .. Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]
                          Int -> Int -> Int -> Int -> STMatrix s Double -> ST s ()
forall s. Int -> Int -> Int -> Int -> STMatrix s Double -> ST s ()
gaussian' Int
m Int
n (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) (Int
pInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) STMatrix s Double
mat
  where
    reduce :: Double -> Int -> Int -> ST s ()
reduce Double
pv Int
p Int
r = do
                      Double
rv <- STMatrix s Double -> Int -> Int -> ST s Double
forall t s. Storable t => STMatrix s t -> Int -> Int -> ST s t
readMatrix STMatrix s Double
mat Int
r Int
p
                      if Double -> Double
forall a. Num a => a -> a
abs Double
rv Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
eps
                        then () -> ST s ()
forall (m :: * -> *) a. Monad m => a -> m a
return ()
                        else
                         let frac :: Double
frac = -Double
rv Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
pv
                             op :: RowOper Double
op = Double -> Int -> Int -> ColRange -> RowOper Double
forall t. t -> Int -> Int -> ColRange -> RowOper t
AXPY Double
frac Int
i Int
r (Int -> ColRange
FromCol Int
p)
                         in do
                             RowOper Double -> STMatrix s Double -> ST s ()
forall t s.
(Num t, Element t) =>
RowOper t -> STMatrix s t -> ST s ()
rowOper RowOper Double
op STMatrix s Double
mat
                             (Int -> ST s ()) -> [Int] -> ST s ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\Int
j' -> STMatrix s Double -> Int -> Int -> (Double -> Double) -> ST s ()
forall t s.
Storable t =>
STMatrix s t -> Int -> Int -> (t -> t) -> ST s ()
modifyMatrix STMatrix s Double
mat Int
r Int
j' (\Double
x -> if Double -> Double
forall a. Num a => a -> a
abs Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
eps then Double
0 else Double
x)) [Int
p..Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]

gaussianFFST :: Int -> Int -> STMatrix s Z -> ST s ()
gaussianFFST :: Int -> Int -> STMatrix s Z -> ST s ()
gaussianFFST Int
m Int
n = Int -> Int -> Int -> Int -> STMatrix s Z -> ST s ()
forall s. Int -> Int -> Int -> Int -> STMatrix s Z -> ST s ()
gaussianFF' Int
m Int
n Int
0 Int
0

-- | Gaussian elimination perfomed in-place in the @'ST'@ monad.
gaussianST :: Int -> Int -> STMatrix s Double -> ST s ()
gaussianST :: Int -> Int -> STMatrix s Double -> ST s ()
gaussianST Int
m Int
n = Int -> Int -> Int -> Int -> STMatrix s Double -> ST s ()
forall s. Int -> Int -> Int -> Int -> STMatrix s Double -> ST s ()
gaussian' Int
m Int
n Int
0 Int
0

rrefST :: Int -> Int -> STMatrix s Z -> ST s ()
rrefST :: Int -> Int -> STMatrix s Z -> ST s ()
rrefST Int
m Int
n STMatrix s Z
mat = do
                    Int -> Int -> Int -> Int -> STMatrix s Z -> ST s ()
forall s. Int -> Int -> Int -> Int -> STMatrix s Z -> ST s ()
gaussianFF' Int
m Int
n Int
0 Int
0 STMatrix s Z
mat
                    Int -> Int -> Int -> Int -> STMatrix s Z -> ST s ()
forall s. Int -> Int -> Int -> Int -> STMatrix s Z -> ST s ()
backwardFF' Int
m Int
n (Int
r'Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) STMatrix s Z
mat
    where
        r' :: Int
r' = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
m Int
n

isref :: (Num a, Ord a, Container Vector a) => Matrix a -> Bool
isref :: Matrix a -> Bool
isref Matrix a
mat = case [(Int, Int)]
pivot of
              []      -> Bool
True
              (Int
r,Int
p):[(Int, Int)]
_ -> (Int
r Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0)
                           Bool -> Bool -> Bool
&&
                             (let leftMat :: Matrix a
leftMat  = Matrix a
mat Matrix a -> (Extractor, Extractor) -> Matrix a
forall t.
Element t =>
Matrix t -> (Extractor, Extractor) -> Matrix t
?? (Int -> Extractor
Drop Int
1, Int -> Extractor
Take (Int
pInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1))
                                  rightMat :: Matrix a
rightMat = Matrix a
mat Matrix a -> (Extractor, Extractor) -> Matrix a
forall t.
Element t =>
Matrix t -> (Extractor, Extractor) -> Matrix t
?? (Int -> Extractor
Drop Int
1, Int -> Extractor
Drop (Int
pInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1))
                                  leftZero :: Bool
leftZero = [(Int, Int)] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null ([(Int, Int)] -> Bool) -> [(Int, Int)] -> Bool
forall a b. (a -> b) -> a -> b
$ (a -> Bool) -> Matrix a -> [IndexOf Matrix]
forall (c :: * -> *) e.
Container c e =>
(e -> Bool) -> c e -> [IndexOf c]
find (a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/=a
0) Matrix a
leftMat
                                  rightRef :: Bool
rightRef = Matrix a -> Bool
forall a. (Num a, Ord a, Container Vector a) => Matrix a -> Bool
isref Matrix a
rightMat
                              in Bool
leftZero Bool -> Bool -> Bool
&& Bool
rightRef)
    where
        pivot :: [IndexOf Matrix]
pivot = (a -> Bool) -> Matrix a -> [IndexOf Matrix]
forall (c :: * -> *) e.
Container c e =>
(e -> Bool) -> c e -> [IndexOf c]
find (a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/=a
0) Matrix a
mat

isrref' :: (Num a, Ord a, Container Vector a) => Int -> Matrix a -> Bool
isrref' :: Int -> Matrix a -> Bool
isrref' Int
r Matrix a
mat = case [(Int, Int)]
pivot of
              []       -> Bool
True
              (Int
r',Int
p):[(Int, Int)]
_ -> (Int
r' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0)
                           Bool -> Bool -> Bool
&& (let leftMat :: Matrix a
leftMat  = Matrix a
subMat Matrix a -> (Extractor, Extractor) -> Matrix a
forall t.
Element t =>
Matrix t -> (Extractor, Extractor) -> Matrix t
?? (Int -> Extractor
Drop Int
1, Int -> Extractor
Take (Int
pInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1))
                                   col :: Matrix a
col      = Matrix a
mat Matrix a -> [Int] -> Matrix a
forall t. Element t => Matrix t -> [Int] -> Matrix t
¿ [Int
p]
                                   colSingleton :: Bool
colSingleton = case (a -> Bool) -> Matrix a -> [IndexOf Matrix]
forall (c :: * -> *) e.
Container c e =>
(e -> Bool) -> c e -> [IndexOf c]
find (a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/=a
0) Matrix a
col of
                                                    [IndexOf Matrix
_] -> Bool
True
                                                    [IndexOf Matrix]
_   -> Bool
False
                                   leftZero :: Bool
leftZero = [(Int, Int)] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null ([(Int, Int)] -> Bool) -> [(Int, Int)] -> Bool
forall a b. (a -> b) -> a -> b
$ (a -> Bool) -> Matrix a -> [IndexOf Matrix]
forall (c :: * -> *) e.
Container c e =>
(e -> Bool) -> c e -> [IndexOf c]
find (a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/=a
0) Matrix a
leftMat
                                   nextRref :: Bool
nextRref = Int -> Matrix a -> Bool
forall a.
(Num a, Ord a, Container Vector a) =>
Int -> Matrix a -> Bool
isrref' (Int
rInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Matrix a
mat
                               in Bool
leftZero Bool -> Bool -> Bool
&& Bool
colSingleton Bool -> Bool -> Bool
&& Bool
nextRref)
    where
        subMat :: Matrix a
subMat = Matrix a
mat Matrix a -> (Extractor, Extractor) -> Matrix a
forall t.
Element t =>
Matrix t -> (Extractor, Extractor) -> Matrix t
?? (Int -> Extractor
Drop Int
r, Extractor
All)
        pivot :: [IndexOf Matrix]
pivot  = (a -> Bool) -> Matrix a -> [IndexOf Matrix]
forall (c :: * -> *) e.
Container c e =>
(e -> Bool) -> c e -> [IndexOf c]
find (a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/=a
0) Matrix a
subMat

isrref :: (Num a, Ord a, Container Vector a) => Matrix a -> Bool
isrref :: Matrix a -> Bool
isrref = Int -> Matrix a -> Bool
forall a.
(Num a, Ord a, Container Vector a) =>
Int -> Matrix a -> Bool
isrref' Int
0

rref :: Matrix Z -> Matrix Z
rref :: Matrix Z -> Matrix Z
rref Matrix Z
mat = (forall s. ST s (Matrix Z)) -> Matrix Z
forall a. (forall s. ST s a) -> a
runST ((forall s. ST s (Matrix Z)) -> Matrix Z)
-> (forall s. ST s (Matrix Z)) -> Matrix Z
forall a b. (a -> b) -> a -> b
$ do
    STMatrix s Z
matST <- Matrix Z -> ST s (STMatrix s Z)
forall t s. Element t => Matrix t -> ST s (STMatrix s t)
thawMatrix Matrix Z
mat
    Int -> Int -> STMatrix s Z -> ST s ()
forall s. Int -> Int -> STMatrix s Z -> ST s ()
rrefST Int
m Int
n STMatrix s Z
matST
    STMatrix s Z -> ST s (Matrix Z)
forall t s. Element t => STMatrix s t -> ST s (Matrix t)
freezeMatrix STMatrix s Z
matST
  where
    m :: Int
m = Matrix Z -> Int
forall t. Matrix t -> Int
rows Matrix Z
mat
    n :: Int
n = Matrix Z -> Int
forall t. Matrix t -> Int
cols Matrix Z
mat

gaussianFF :: Matrix Z -> Matrix Z
gaussianFF :: Matrix Z -> Matrix Z
gaussianFF Matrix Z
mat = (forall s. ST s (Matrix Z)) -> Matrix Z
forall a. (forall s. ST s a) -> a
runST ((forall s. ST s (Matrix Z)) -> Matrix Z)
-> (forall s. ST s (Matrix Z)) -> Matrix Z
forall a b. (a -> b) -> a -> b
$ do
    STMatrix s Z
matST <- Matrix Z -> ST s (STMatrix s Z)
forall t s. Element t => Matrix t -> ST s (STMatrix s t)
thawMatrix Matrix Z
mat
    Int -> Int -> STMatrix s Z -> ST s ()
forall s. Int -> Int -> STMatrix s Z -> ST s ()
gaussianFFST Int
m Int
n STMatrix s Z
matST
    STMatrix s Z -> ST s (Matrix Z)
forall t s. Element t => STMatrix s t -> ST s (Matrix t)
freezeMatrix STMatrix s Z
matST
  where
    m :: Int
m = Matrix Z -> Int
forall t. Matrix t -> Int
rows Matrix Z
mat
    n :: Int
n = Matrix Z -> Int
forall t. Matrix t -> Int
cols Matrix Z
mat

-- | Gaussian elimination as pure function. Involves a copy of the input matrix.
--
-- @
-- &#x3BB; let mat = (3 >< 4) [1, 1, -2, 0, 0, 2, -6, -4, 3, 0, 3, 1]
-- &#x3BB; mat
-- (3><4)
--  [ 1.0, 1.0, -2.0,  0.0
--  , 0.0, 2.0, -6.0, -4.0
--  , 3.0, 0.0,  3.0,  1.0 ]
-- &#x3BB; gaussian mat
-- (3><4)
--  [ 3.0, 0.0,  3.0,                1.0
--  , 0.0, 2.0, -6.0,               -4.0
--  , 0.0, 0.0,  0.0, 1.6666666666666667 ]
-- @
--
gaussian :: Matrix Double -> Matrix Double
gaussian :: Matrix Double -> Matrix Double
gaussian Matrix Double
mat = (forall s. ST s (Matrix Double)) -> Matrix Double
forall a. (forall s. ST s a) -> a
runST ((forall s. ST s (Matrix Double)) -> Matrix Double)
-> (forall s. ST s (Matrix Double)) -> Matrix Double
forall a b. (a -> b) -> a -> b
$ do
    STMatrix s Double
matST <- Matrix Double -> ST s (STMatrix s Double)
forall t s. Element t => Matrix t -> ST s (STMatrix s t)
thawMatrix Matrix Double
mat
    Int -> Int -> STMatrix s Double -> ST s ()
forall s. Int -> Int -> STMatrix s Double -> ST s ()
gaussianST Int
m Int
n STMatrix s Double
matST
    STMatrix s Double -> ST s (Matrix Double)
forall t s. Element t => STMatrix s t -> ST s (Matrix t)
freezeMatrix STMatrix s Double
matST
  where
    m :: Int
m = Matrix Double -> Int
forall t. Matrix t -> Int
rows Matrix Double
mat
    n :: Int
n = Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
mat

independentColumnsRREF :: Matrix Z -> [Int]
independentColumnsRREF :: Matrix Z -> [Int]
independentColumnsRREF Matrix Z
mat = Matrix Z -> [Int]
pivotsUFF Matrix Z
mat'
    where
        mat' :: Matrix Z
mat' = Matrix Z -> Matrix Z
rref Matrix Z
mat

independentColumnsFF :: Matrix Z -> [Int]
independentColumnsFF :: Matrix Z -> [Int]
independentColumnsFF Matrix Z
mat = Matrix Z -> [Int]
pivotsUFF Matrix Z
mat'
    where
        mat' :: Matrix Z
mat' = Matrix Z -> Matrix Z
gaussianFF Matrix Z
mat

independentColumnsVerifiedFF :: Matrix Z -> [Int]
independentColumnsVerifiedFF :: Matrix Z -> [Int]
independentColumnsVerifiedFF Matrix Z
mat
        | Matrix Z -> Bool
forall a. (Num a, Ord a, Container Vector a) => Matrix a -> Bool
isref Matrix Z
mat' Bool -> Bool -> Bool
&& Matrix Z -> Matrix Z -> Bool
verify Matrix Z
mat Matrix Z
mat'
                     = Matrix Z -> [Int]
pivotsUFF Matrix Z
mat'
        | Bool
otherwise  = [Char] -> [Int]
forall a. HasCallStack => [Char] -> a
error [Char]
"could not verify row echelon form"
    where
        mat' :: Matrix Z
mat' = Matrix Z -> Matrix Z
gaussianFF Matrix Z
mat

-- | Returns the indices of a maximal linearly independent subset of the columns
--   in the matrix.
--
-- @
-- &#x3BB; let mat = (3 >< 4) [1, 1, -2, 0, 0, 2, -6, -4, 3, 0, 3, 1]
-- &#x3BB; mat
-- (3><4)
--  [ 1.0, 1.0, -2.0,  0.0
--  , 0.0, 2.0, -6.0, -4.0
--  , 3.0, 0.0,  3.0,  1.0 ]
-- &#x3BB; independentColumns mat
-- [0,1,3]
-- @
--
independentColumns :: Matrix Double -> [Int]
independentColumns :: Matrix Double -> [Int]
independentColumns Matrix Double
mat = Matrix Double -> [Int]
pivotsU Matrix Double
mat'
    where
        mat' :: Matrix Double
mat' = Matrix Double -> Matrix Double
gaussian Matrix Double
mat

verify :: Matrix Z -> Matrix Z -> Bool
verify :: Matrix Z -> Matrix Z -> Bool
verify Matrix Z
mat Matrix Z
ref = Int
rank1 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
rank2 Bool -> Bool -> Bool
&& Int
rank1 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
rank3
    where
        matD :: Matrix Double
matD = Matrix Z -> Matrix Double
forall (c :: * -> *) e. Container c e => c Z -> c e
fromZ Matrix Z
mat :: Matrix Double
        refD :: Matrix Double
refD = Matrix Z -> Matrix Double
forall (c :: * -> *) e. Container c e => c Z -> c e
fromZ Matrix Z
ref :: Matrix Double
        rank1 :: Int
rank1 = Matrix Double -> Int
forall t. Field t => Matrix t -> Int
rank Matrix Double
matD
        rank2 :: Int
rank2 = Matrix Double -> Int
forall t. Field t => Matrix t -> Int
rank Matrix Double
refD
        rank3 :: Int
rank3 = Matrix Double -> Int
forall t. Field t => Matrix t -> Int
rank (Matrix Double -> Int) -> Matrix Double -> Int
forall a b. (a -> b) -> a -> b
$ Matrix Double
matD Matrix Double -> Matrix Double -> Matrix Double
forall t. Element t => Matrix t -> Matrix t -> Matrix t
=== Matrix Double
refD

independentColumnsMatRREF :: Matrix Z -> Matrix Z
independentColumnsMatRREF :: Matrix Z -> Matrix Z
independentColumnsMatRREF Matrix Z
mat =
  case Matrix Z -> [Int]
independentColumnsRREF Matrix Z
mat of
    [] -> (Matrix Z -> Int
forall t. Matrix t -> Int
rows Matrix Z
mat Int -> Int -> [Z] -> Matrix Z
forall a. Storable a => Int -> Int -> [a] -> Matrix a
>< Int
1) ([Z] -> Matrix Z) -> [Z] -> Matrix Z
forall a b. (a -> b) -> a -> b
$ Z -> [Z]
forall a. a -> [a]
repeat Z
0
    [Int]
cs -> Matrix Z
mat Matrix Z -> [Int] -> Matrix Z
forall t. Element t => Matrix t -> [Int] -> Matrix t
¿ [Int]
cs

independentColumnsMatFF :: Matrix Z -> Matrix Z
independentColumnsMatFF :: Matrix Z -> Matrix Z
independentColumnsMatFF Matrix Z
mat =
  case Matrix Z -> [Int]
independentColumnsFF Matrix Z
mat of
    [] -> (Matrix Z -> Int
forall t. Matrix t -> Int
rows Matrix Z
mat Int -> Int -> [Z] -> Matrix Z
forall a. Storable a => Int -> Int -> [a] -> Matrix a
>< Int
1) ([Z] -> Matrix Z) -> [Z] -> Matrix Z
forall a b. (a -> b) -> a -> b
$ Z -> [Z]
forall a. a -> [a]
repeat Z
0
    [Int]
cs -> Matrix Z
mat Matrix Z -> [Int] -> Matrix Z
forall t. Element t => Matrix t -> [Int] -> Matrix t
¿ [Int]
cs

-- | Returns a sub matrix containing a maximal linearly independent subset of
--   the columns in the matrix.
--
-- @
-- &#x3BB; let mat = (3 >< 4) [1, 1, -2, 0, 0, 2, -6, -4, 3, 0, 3, 1]
-- &#x3BB; mat
-- (3><4)
--  [ 1.0, 1.0, -2.0,  0.0
--  , 0.0, 2.0, -6.0, -4.0
--  , 3.0, 0.0,  3.0,  1.0 ]
-- &#x3BB; independentColumnsMat mat
-- (3><3)
--  [ 1.0, 1.0,  0.0
--  , 0.0, 2.0, -4.0
--  , 3.0, 0.0,  1.0 ]
-- @
--
independentColumnsMat :: Matrix Double -> Matrix Double
independentColumnsMat :: Matrix Double -> Matrix Double
independentColumnsMat Matrix Double
mat =
  case Matrix Double -> [Int]
independentColumns Matrix Double
mat of
    [] -> (Matrix Double -> Int
forall t. Matrix t -> Int
rows Matrix Double
mat Int -> Int -> [Double] -> Matrix Double
forall a. Storable a => Int -> Int -> [a] -> Matrix a
>< Int
1) ([Double] -> Matrix Double) -> [Double] -> Matrix Double
forall a b. (a -> b) -> a -> b
$ Double -> [Double]
forall a. a -> [a]
repeat Double
0
    [Int]
cs -> Matrix Double
mat Matrix Double -> [Int] -> Matrix Double
forall t. Element t => Matrix t -> [Int] -> Matrix t
¿ [Int]
cs