-- SPDX-License-Identifier: MPL-2.0

{- | Pure reference decompositions with explicit shape and tolerance contracts.
Cholesky uses a symmetric positive-definite matrix; unpivoted LDL^T requires
nonzero scalar pivots. Bidiagonalization supports rectangular real matrices.
-}
module HMath.Matrix.Decomposition (
    gaussianElimination,
    gaussJordan,
    cholesky,
    ldlt,
    givensQR,
    Bidiagonal (..),
    bidiagonalize,
    tridiagonalize,
    Iteration (..),
    powerIteration,
    inverseIteration,
    rayleighIteration,
    Arnoldi (..),
    arnoldi,
    lanczos,
    qrIteration,
    jordanChains,
    jordanForm,
) where

import Control.Monad (foldM)
import Data.List (findIndex, maximumBy, nub)
import Data.Ord (comparing)
import HMath.Error
import HMath.Matrix
import HMath.Matrix.Internal
import HMath.Matrix.LAPACK (Hessenberg (..), hessenberg)
import HMath.Matrix.LU (inverse, luDecomposeWithTolerance, solve, solveLU)

{- | Given finite nonnegative relative tolerance then finite matrix, return normalized
row-echelon form and zero-based pivot columns. Partial row pivoting skips pivots below
tolerance times the global entry scale. Invalid tolerance/non-finite arithmetic fails.
-}
gaussianElimination :: (RealFloat a) => a -> Matrix a -> Either MathError (Matrix a, [Int])
gaussianElimination :: forall a.
RealFloat a =>
a -> Matrix a -> Either MathError (Matrix a, [Int])
gaussianElimination a
tol Matrix a
a = Bool -> a -> Matrix a -> Either MathError (Matrix a, [Int])
forall a.
RealFloat a =>
Bool -> a -> Matrix a -> Either MathError (Matrix a, [Int])
eliminate Bool
False a
tol Matrix a
a

{- | Given relative tolerance then finite matrix, return reduced row-echelon form and
zero-based pivot columns using row pivoting. Invalid tolerance/non-finite arithmetic
fails.
-}
gaussJordan :: (RealFloat a) => a -> Matrix a -> Either MathError (Matrix a, [Int])
gaussJordan :: forall a.
RealFloat a =>
a -> Matrix a -> Either MathError (Matrix a, [Int])
gaussJordan a
tol Matrix a
a = Bool -> a -> Matrix a -> Either MathError (Matrix a, [Int])
forall a.
RealFloat a =>
Bool -> a -> Matrix a -> Either MathError (Matrix a, [Int])
eliminate Bool
True a
tol Matrix a
a

eliminate :: (RealFloat a) => Bool -> a -> Matrix a -> Either MathError (Matrix a, [Int])
eliminate :: forall a.
RealFloat a =>
Bool -> a -> Matrix a -> Either MathError (Matrix a, [Int])
eliminate Bool
reduced a
tol Matrix a
a
    | a
tol a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0 Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
tol Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
tol = MathError -> Either MathError (Matrix a, [Int])
forall a b. a -> Either a b
Left MathError
InvalidTolerance
    | Bool
otherwise = do
        Matrix a -> Either MathError ()
forall a. RealFloat a => Matrix a -> Either MathError ()
finiteMatrix Matrix a
a
        (rows, pivots) <- Int -> Int -> [[a]] -> [Int] -> Either MathError ([[a]], [Int])
go Int
0 Int
0 (Matrix a -> [[a]]
forall a. Matrix a -> [[a]]
entries Matrix a
a) []
        result <- fromListChecked m n (concat rows)
        finiteMatrix result
        Right (result, pivots)
  where
    (Int
m, Int
n) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a
    scale :: a
scale = [a] -> a
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (a
0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map a -> a
forall a. Num a => a -> a
abs (Matrix a -> [a]
forall a. Matrix a -> [a]
toList Matrix a
a))
    go :: Int -> Int -> [[a]] -> [Int] -> Either MathError ([[a]], [Int])
go Int
r Int
c [[a]]
rows [Int]
pivots
        | Int
r Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
m Bool -> Bool -> Bool
|| Int
c Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
n = ([[a]], [Int]) -> Either MathError ([[a]], [Int])
forall a b. b -> Either a b
Right ([[a]]
rows, [Int]
pivots)
        | Bool
otherwise = do
            let pivot :: Int
pivot = (Int, a) -> Int
forall a b. (a, b) -> a
fst (((Int, a) -> (Int, a) -> Ordering) -> [(Int, a)] -> (Int, a)
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
maximumBy (((Int, a) -> a) -> (Int, a) -> (Int, a) -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing (Int, a) -> a
forall a b. (a, b) -> b
snd) [(Int
i, a -> a
forall a. Num a => a -> a
abs ([[a]]
rows [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
c)) | Int
i <- [Int
r .. Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]])
                value :: a
value = [[a]]
rows [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
pivot [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
c
            if a
value a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 Bool -> Bool -> Bool
|| a -> a
forall a. Num a => a -> a
abs a
value a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
tol a -> a -> a
forall a. Num a => a -> a -> a
* a
scale
                then Int -> Int -> [[a]] -> [Int] -> Either MathError ([[a]], [Int])
go Int
r (Int
c Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [[a]]
rows [Int]
pivots
                else do
                    let swapped :: [[a]]
swapped = Int -> Int -> [[a]] -> [[a]]
forall a. Int -> Int -> [a] -> [a]
swap Int
r Int
pivot [[a]]
rows
                        normalized :: [a]
normalized = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
value) ([[a]]
swapped [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
r)
                        update :: Int -> [a] -> [a]
update Int
i [a]
row
                            | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
r = [a]
normalized
                            | (Bool
reduced Bool -> Bool -> Bool
&& Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
r) Bool -> Bool -> Bool
|| Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
r = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [a]
row ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (([a]
row [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
c) a -> a -> a
forall a. Num a => a -> a -> a
*) [a]
normalized)
                            | Bool
otherwise = [a]
row
                        next :: [[a]]
next = (Int -> [a] -> [a]) -> [Int] -> [[a]] -> [[a]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Int -> [a] -> [a]
update [Int
0 ..] [[a]]
swapped
                    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues ([[a]] -> [a]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[a]]
next)
                    Int -> Int -> [[a]] -> [Int] -> Either MathError ([[a]], [Int])
go (Int
r Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) (Int
c Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [[a]]
next ([Int]
pivots [Int] -> [Int] -> [Int]
forall a. [a] -> [a] -> [a]
++ [Int
c])
checkSymmetric :: (RealFloat a) => Matrix a -> Either MathError Int
checkSymmetric :: forall a. RealFloat a => Matrix a -> Either MathError Int
checkSymmetric Matrix a
a
    | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
m = MathError -> Either MathError Int
forall a b. a -> Either a b
Left ((Int, Int) -> MathError
NonSquareMatrix (Int
n, Int
m))
    | ((a, a) -> Bool) -> [(a, a)] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any
        (\(a
x, a
y) -> a -> a
forall a. Num a => a -> a
abs (a
x a -> a -> a
forall a. Num a => a -> a -> a
- a
y) a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> Matrix a -> a
forall a. RealFloat a => Matrix a -> a
defaultTolerance Matrix a
a a -> a -> a
forall a. Num a => a -> a -> a
* a -> a -> a
forall a. Ord a => a -> a -> a
max a
1 (a -> a -> a
forall a. Ord a => a -> a -> a
max (a -> a
forall a. Num a => a -> a
abs a
x) (a -> a
forall a. Num a => a -> a
abs a
y)))
        ([a] -> [a] -> [(a, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip (Matrix a -> [a]
forall a. Matrix a -> [a]
toList Matrix a
a) (Matrix a -> [a]
forall a. Matrix a -> [a]
toList (Matrix a -> Matrix a
forall a. Matrix a -> Matrix a
transposeMatrix Matrix a
a))) =
        MathError -> Either MathError Int
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"symmetric matrix required")
    | Bool
otherwise = Matrix a -> Either MathError ()
forall a. RealFloat a => Matrix a -> Either MathError ()
finiteMatrix Matrix a
a Either MathError () -> Either MathError Int -> Either MathError Int
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Int -> Either MathError Int
forall a b. b -> Either a b
Right Int
n
  where
    (Int
n, Int
m) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a

{- | Factor a finite symmetric positive-definite matrix as L*transpose(L), returning
lower-triangular L. Non-square/asymmetric/non-finite input or a nonpositive diagonal
pivot returns Left MathError.
-}
cholesky :: (RealFloat a) => Matrix a -> Either MathError (Matrix a)
cholesky :: forall a. RealFloat a => Matrix a -> Either MathError (Matrix a)
cholesky Matrix a
a = do
    n <- Matrix a -> Either MathError Int
forall a. RealFloat a => Matrix a -> Either MathError Int
checkSymmetric Matrix a
a
    rows <- foldM (step n) [] [0 .. n - 1]
    fromListChecked n n (concat rows)
  where
    step :: Int -> [[a]] -> Int -> Either MathError [[a]]
step Int
n [[a]]
rows Int
i = do
        prefix <- ([a] -> Int -> Either MathError [a])
-> [a] -> [Int] -> Either MathError [a]
forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM ([[a]] -> Int -> [a] -> Int -> Either MathError [a]
entry [[a]]
rows Int
i) [] [Int
0 .. Int
i]
        Right (rows ++ [prefix ++ replicate (n - i - 1) 0])
    entry :: [[a]] -> Int -> [a] -> Int -> Either MathError [a]
entry [[a]]
rows Int
i [a]
xs Int
j = do
        let value :: a
value = Matrix a -> Int -> Int -> a
forall a. Matrix a -> Int -> Int -> a
get Matrix a
a Int
i Int
j a -> a -> a
forall a. Num a => a -> a -> a
- [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [[a]
xs [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
k a -> a -> a
forall a. Num a => a -> a -> a
* (if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then [a]
xs [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
k else [[a]]
rows [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
k) | Int
k <- [Int
0 .. Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
        if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j Bool -> Bool -> Bool
&& a
value a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
0
            then MathError -> Either MathError [a]
forall a b. a -> Either a b
Left (Int -> MathError
NotPositiveDefinite Int
i)
            else do
                let result :: a
result = if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then a -> a
forall a. Floating a => a -> a
sqrt a
value else a
value a -> a -> a
forall a. Fractional a => a -> a -> a
/ ([[a]]
rows [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
j)
                [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues [a
result]
                [a] -> Either MathError [a]
forall a b. b -> Either a b
Right ([a]
xs [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a
result])

{- | Factor a finite symmetric matrix as L*diag(d)*transpose(L), returning unit lower L and
diagonal list d. This unpivoted method requires nonzero scalar pivots; it can reject
invertible indefinite matrices needing block pivots.
-}
ldlt :: (RealFloat a) => Matrix a -> Either MathError (Matrix a, [a])
ldlt :: forall a.
RealFloat a =>
Matrix a -> Either MathError (Matrix a, [a])
ldlt Matrix a
a = do
    n <- Matrix a -> Either MathError Int
forall a. RealFloat a => Matrix a -> Either MathError Int
checkSymmetric Matrix a
a
    (rows, ds) <- foldM (step n) ([], []) [0 .. n - 1]
    lower <- fromListChecked n n (concat rows)
    Right (lower, ds)
  where
    step :: Int -> ([[a]], [a]) -> Int -> Either MathError ([[a]], [a])
step Int
n ([[a]]
rows, [a]
ds) Int
i = do
        prefix <-
            ([a] -> Int -> Either MathError [a])
-> [a] -> [Int] -> Either MathError [a]
forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM
                ( \[a]
xs Int
j -> do
                    let v :: a
v = (Matrix a -> Int -> Int -> a
forall a. Matrix a -> Int -> Int -> a
get Matrix a
a Int
i Int
j a -> a -> a
forall a. Num a => a -> a -> a
- [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [[a]
xs [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
k a -> a -> a
forall a. Num a => a -> a -> a
* [[a]]
rows [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
k a -> a -> a
forall a. Num a => a -> a -> a
* [a]
ds [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
k | Int
k <- [Int
0 .. Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]) a -> a -> a
forall a. Fractional a => a -> a -> a
/ ([a]
ds [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
j)
                    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues [a
v]
                    [a] -> Either MathError [a]
forall a b. b -> Either a b
Right ([a]
xs [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a
v])
                )
                []
                [Int
0 .. Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
        let d = Matrix a -> Int -> Int -> a
forall a. Matrix a -> Int -> Int -> a
get Matrix a
a Int
i Int
i a -> a -> a
forall a. Num a => a -> a -> a
- [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [[a]
prefix [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
k a -> a -> a
forall a. Num a => a -> a -> a
* [a]
prefix [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
k a -> a -> a
forall a. Num a => a -> a -> a
* [a]
ds [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
k | Int
k <- [Int
0 .. Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
        if d == 0
            then Left (SingularMatrix i)
            else do
                finiteValues [d]
                Right (rows ++ [prefix ++ 1 : replicate (n - i - 1) 0], ds ++ [d])

{- | Factor finite m-by-n matrix as Q*R using Givens rotations, returning square m-by-m Q
and m-by-n R without column pivoting. Non-finite input/intermediate results fail.
-}
givensQR :: (RealFloat a) => Matrix a -> Either MathError (Matrix a, Matrix a)
givensQR :: forall a.
RealFloat a =>
Matrix a -> Either MathError (Matrix a, Matrix a)
givensQR Matrix a
a = do
    Matrix a -> Either MathError ()
forall a. RealFloat a => Matrix a -> Either MathError ()
finiteMatrix Matrix a
a
    q0 <- Int -> Either MathError (Matrix a)
forall a. Num a => Int -> Either MathError (Matrix a)
identityChecked Int
m
    let rotatePair ([[a]]
q, [[a]]
r) (Int
j, Int
i) = do
            let x :: a
x = [[a]]
r [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
j; y :: a
y = [[a]]
r [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
j; radius :: a
radius = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a
x, a
y]
            if a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
radius Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
radius
                then MathError -> Either MathError ([[a]], [[a]])
forall a b. a -> Either a b
Left MathError
NonFiniteResult
                else
                    if a
radius a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0
                        then ([[a]], [[a]]) -> Either MathError ([[a]], [[a]])
forall a b. b -> Either a b
Right ([[a]]
q, [[a]]
r)
                        else do
                            let c :: a
c = a
x a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
radius
                                s :: a
s = a
y a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
radius
                                rowJ :: [a]
rowJ = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\a
u a
v -> a
c a -> a -> a
forall a. Num a => a -> a -> a
* a
u a -> a -> a
forall a. Num a => a -> a -> a
+ a
s a -> a -> a
forall a. Num a => a -> a -> a
* a
v) ([[a]]
r [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) ([[a]]
r [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i)
                                rowI :: [a]
rowI = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\a
u a
v -> -a
s a -> a -> a
forall a. Num a => a -> a -> a
* a
u a -> a -> a
forall a. Num a => a -> a -> a
+ a
c a -> a -> a
forall a. Num a => a -> a -> a
* a
v) ([[a]]
r [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) ([[a]]
r [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i)
                                rnew :: [[a]]
rnew = Int -> [a] -> [[a]] -> [[a]]
forall a. Int -> a -> [a] -> [a]
replace Int
i [a]
rowI (Int -> [a] -> [[a]] -> [[a]]
forall a. Int -> a -> [a] -> [a]
replace Int
j [a]
rowJ [[a]]
r)
                                qnew :: [[a]]
qnew =
                                    ([a] -> [a]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map
                                        ( \[a]
row -> Int -> a -> [a] -> [a]
forall a. Int -> a -> [a] -> [a]
replace Int
i (-a
s a -> a -> a
forall a. Num a => a -> a -> a
* ([a]
row [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) a -> a -> a
forall a. Num a => a -> a -> a
+ a
c a -> a -> a
forall a. Num a => a -> a -> a
* ([a]
row [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
i)) (Int -> a -> [a] -> [a]
forall a. Int -> a -> [a] -> [a]
replace Int
j (a
c a -> a -> a
forall a. Num a => a -> a -> a
* ([a]
row [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) a -> a -> a
forall a. Num a => a -> a -> a
+ a
s a -> a -> a
forall a. Num a => a -> a -> a
* ([a]
row [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
i)) [a]
row)
                                        )
                                        [[a]]
q
                            [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues ([[a]] -> [a]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[a]]
rnew [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [[a]] -> [a]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[a]]
qnew)
                            ([[a]], [[a]]) -> Either MathError ([[a]], [[a]])
forall a b. b -> Either a b
Right ([[a]]
qnew, [[a]]
rnew)
    (q, r) <-
        foldM
            rotatePair
            (entries q0, entries a)
            [(j, i) | j <- [0 .. min m n - 1], i <- reverse [j + 1 .. m - 1]]
    (,) <$> fromListChecked m m (concat q) <*> fromListChecked m n (concat r)
  where
    (Int
m, Int
n) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a

-- | Orthogonal factors A=U*B*transpose(V), retaining a rectangular bidiagonal B.
data Bidiagonal a = Bidiagonal
    { forall a. Bidiagonal a -> Matrix a
bidiagonalU :: Matrix a
    -- ^ Left orthogonal factor in A=U*B*transpose(V).
    , forall a. Bidiagonal a -> Matrix a
bidiagonalB :: Matrix a
    -- ^ Bidiagonal matrix in the factorization, retaining input shape.
    , forall a. Bidiagonal a -> Matrix a
bidiagonalV :: Matrix a
    -- ^ Right orthogonal factor in A=U*B*transpose(V).
    }
    deriving (Bidiagonal a -> Bidiagonal a -> Bool
(Bidiagonal a -> Bidiagonal a -> Bool)
-> (Bidiagonal a -> Bidiagonal a -> Bool) -> Eq (Bidiagonal a)
forall a. Eq a => Bidiagonal a -> Bidiagonal a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall a. Eq a => Bidiagonal a -> Bidiagonal a -> Bool
== :: Bidiagonal a -> Bidiagonal a -> Bool
$c/= :: forall a. Eq a => Bidiagonal a -> Bidiagonal a -> Bool
/= :: Bidiagonal a -> Bidiagonal a -> Bool
Eq, Int -> Bidiagonal a -> ShowS
[Bidiagonal a] -> ShowS
Bidiagonal a -> String
(Int -> Bidiagonal a -> ShowS)
-> (Bidiagonal a -> String)
-> ([Bidiagonal a] -> ShowS)
-> Show (Bidiagonal a)
forall a. Show a => Int -> Bidiagonal a -> ShowS
forall a. Show a => [Bidiagonal a] -> ShowS
forall a. Show a => Bidiagonal a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> Bidiagonal a -> ShowS
showsPrec :: Int -> Bidiagonal a -> ShowS
$cshow :: forall a. Show a => Bidiagonal a -> String
show :: Bidiagonal a -> String
$cshowList :: forall a. Show a => [Bidiagonal a] -> ShowS
showList :: [Bidiagonal a] -> ShowS
Show)

{- | Reduce a finite rectangular matrix to bidiagonal form with orthogonal factors,
returning Bidiagonal fields. Input/numeric validation errors propagate; inspect the
factor types and module reconstruction example.
-}
bidiagonalize :: (RealFloat a) => Matrix a -> Either MathError (Bidiagonal a)
bidiagonalize :: forall a.
RealFloat a =>
Matrix a -> Either MathError (Bidiagonal a)
bidiagonalize Matrix a
a = do
    Matrix a -> Either MathError ()
forall a. RealFloat a => Matrix a -> Either MathError ()
finiteMatrix Matrix a
a
    u0 <- Int -> Either MathError (Matrix a)
forall a. Num a => Int -> Either MathError (Matrix a)
identityChecked Int
m
    v0 <- identityChecked n
    (u, b, v) <- foldM step (u0, a, v0) [0 .. min m n - 1]
    Right (Bidiagonal u b v)
  where
    (Int
m, Int
n) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a
    reflector :: Int -> Int -> [a] -> Either MathError (Matrix a)
reflector Int
size Int
start [a]
xs = do
        let s :: a
s = [a] -> a
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (a
0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map a -> a
forall a. Num a => a -> a
abs [a]
xs)
        if a
s a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0
            then Int -> Either MathError (Matrix a)
forall a. Num a => Int -> Either MathError (Matrix a)
identityChecked Int
size
            else do
                let y :: [a]
y = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
s) [a]
xs
                    alpha :: a
alpha = if [a]
y [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
0 a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
0 then a -> a
forall a. Num a => a -> a
negate ([a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
y) else [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
y
                    z :: [a]
z = Int -> a -> [a] -> [a]
forall a. Int -> a -> [a] -> [a]
replace Int
0 ([a]
y [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
0 a -> a -> a
forall a. Num a => a -> a -> a
- a
alpha) [a]
y
                    unit :: [a]
unit = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
forall a. Fractional a => a -> a -> a
/ [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
z) [a]
z
                    vv :: [a]
vv = Int -> a -> [a]
forall a. Int -> a -> [a]
replicate Int
start a
0 [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a]
unit
                Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
forall a.
Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
build Int
size Int
size (\Int
i Int
j -> (if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then a
1 else a
0) a -> a -> a
forall a. Num a => a -> a -> a
- a
2 a -> a -> a
forall a. Num a => a -> a -> a
* ([a]
vv [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) a -> a -> a
forall a. Num a => a -> a -> a
* ([a]
vv [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
j))
    step :: (Matrix a, Matrix a, Matrix a)
-> Int -> Either MathError (Matrix a, Matrix a, Matrix a)
step (Matrix a
u, Matrix a
b, Matrix a
v) Int
k = do
        h <- Int -> Int -> [a] -> Either MathError (Matrix a)
forall {a}.
RealFloat a =>
Int -> Int -> [a] -> Either MathError (Matrix a)
reflector Int
m Int
k [Matrix a -> Int -> Int -> a
forall a. Matrix a -> Int -> Int -> a
get Matrix a
b Int
i Int
k | Int
i <- [Int
k .. Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
        b1 <- multiplyChecked h b
        u1 <- multiplyChecked u h
        finiteValues (toList b1 ++ toList u1)
        if k + 1 >= n
            then Right (u1, b1, v)
            else do
                g <- reflector n (k + 1) [get b1 k j | j <- [k + 1 .. n - 1]]
                b2 <- multiplyChecked b1 g
                v1 <- multiplyChecked v g
                finiteValues (toList b2 ++ toList v1)
                Right (u1, b2, v1)

{- | Reduce a symmetric finite square matrix to @(Q,T)@, satisfying
@A = Q * T * transpose Q@ with tridiagonal @T@. Non-square/asymmetric input
or numerical failure returns Left MathError.
-}
tridiagonalize :: Matrix Double -> Either MathError (Matrix Double, Matrix Double)
tridiagonalize :: Matrix Double -> Either MathError (Matrix Double, Matrix Double)
tridiagonalize Matrix Double
a = do
    _ <- Matrix Double -> Either MathError Int
forall a. RealFloat a => Matrix a -> Either MathError Int
checkSymmetric Matrix Double
a
    f <- hessenberg a
    Right (hessenbergQ f, hessenbergH f)

{- | Accepted eigenpair estimate with residual and iteration count. See the iterative
solver's tolerance and failure contract.
-}
data Iteration a = Iteration
    { forall a. Iteration a -> a
iterationValue :: a
    -- ^ Returned real eigenvalue estimate.
    , forall a. Iteration a -> [a]
iterationVector :: [a]
    -- ^ Returned normalized eigenvector estimate.
    , forall a. Iteration a -> a
iterationResidual :: a
    -- ^ Eigenpair residual norm used to assess convergence.
    , forall a. Iteration a -> Int
iterationCount :: Int
    -- ^ Number of eigen-iteration updates performed.
    }
    deriving (Iteration a -> Iteration a -> Bool
(Iteration a -> Iteration a -> Bool)
-> (Iteration a -> Iteration a -> Bool) -> Eq (Iteration a)
forall a. Eq a => Iteration a -> Iteration a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall a. Eq a => Iteration a -> Iteration a -> Bool
== :: Iteration a -> Iteration a -> Bool
$c/= :: forall a. Eq a => Iteration a -> Iteration a -> Bool
/= :: Iteration a -> Iteration a -> Bool
Eq, Int -> Iteration a -> ShowS
[Iteration a] -> ShowS
Iteration a -> String
(Int -> Iteration a -> ShowS)
-> (Iteration a -> String)
-> ([Iteration a] -> ShowS)
-> Show (Iteration a)
forall a. Show a => Int -> Iteration a -> ShowS
forall a. Show a => [Iteration a] -> ShowS
forall a. Show a => Iteration a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> Iteration a -> ShowS
showsPrec :: Int -> Iteration a -> ShowS
$cshow :: forall a. Show a => Iteration a -> String
show :: Iteration a -> String
$cshowList :: forall a. Show a => [Iteration a] -> ShowS
showList :: [Iteration a] -> ShowS
Show)

iterateEigen ::
    (RealFloat a) =>
    String ->
    a ->
    Int ->
    (a -> [a] -> Either MathError [a]) ->
    Matrix a ->
    [a] ->
    Either MathError (Iteration a)
iterateEigen :: forall a.
RealFloat a =>
String
-> a
-> Int
-> (a -> [a] -> Either MathError [a])
-> Matrix a
-> [a]
-> Either MathError (Iteration a)
iterateEigen String
name a
tol Int
budget a -> [a] -> Either MathError [a]
action Matrix a
a [a]
initial
    | Int
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n = MathError -> Either MathError (Iteration a)
forall a b. a -> Either a b
Left ((Int, Int) -> MathError
NonSquareMatrix (Int
m, Int
n))
    | [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
initial Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n = MathError -> Either MathError (Iteration a)
forall a b. a -> Either a b
Left (Int -> Int -> MathError
VectorLengthMismatch Int
n ([a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
initial))
    | a
tol a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
0 Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
tol Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
tol = MathError -> Either MathError (Iteration a)
forall a b. a -> Either a b
Left MathError
InvalidTolerance
    | Int
budget Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
1 = MathError -> Either MathError (Iteration a)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"eigen iteration: positive budget required")
    | Bool
otherwise = do
        Matrix a -> Either MathError ()
forall a. RealFloat a => Matrix a -> Either MathError ()
finiteMatrix Matrix a
a
        [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues [a]
initial
        let norm :: a
norm = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
initial
        if a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
norm Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
norm Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
matrixScale Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
matrixScale
            then MathError -> Either MathError (Iteration a)
forall a b. a -> Either a b
Left MathError
NonFiniteResult
            else
                if a
norm a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0
                    then MathError -> Either MathError (Iteration a)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"eigen iteration: initial vector is zero")
                    else Int -> [a] -> Either MathError (Iteration a)
go Int
0 ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
norm) [a]
initial)
  where
    (Int
m, Int
n) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a
    go :: Int -> [a] -> Either MathError (Iteration a)
go Int
k [a]
x = do
        let ax :: [a]
ax = [[a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
dotV [a]
row [a]
x | [a]
row <- Matrix a -> [[a]]
forall a. Matrix a -> [[a]]
entries Matrix a
a]
            lambda :: a
lambda = [a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
dotV [a]
x [a]
ax
            residual :: a
residual = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\a
u a
v -> a
u a -> a -> a
forall a. Num a => a -> a -> a
- a
lambda a -> a -> a
forall a. Num a => a -> a -> a
* a
v) [a]
ax [a]
x)
        [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues (a
lambda a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a
residual a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a]
ax)
        if a
residual a -> a -> a
forall a. Fractional a => a -> a -> a
/ a -> a -> a
forall a. Ord a => a -> a -> a
max a
1 a
matrixScale a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
tol
            then Iteration a -> Either MathError (Iteration a)
forall a b. b -> Either a b
Right (a -> [a] -> a -> Int -> Iteration a
forall a. a -> [a] -> a -> Int -> Iteration a
Iteration a
lambda [a]
x a
residual Int
k)
            else
                if Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
budget
                    then MathError -> Either MathError (Iteration a)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
name Int
budget)
                    else do
                        y <- a -> [a] -> Either MathError [a]
action a
lambda [a]
x
                        let norm = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
y
                        if isNaN norm || isInfinite norm
                            then Left NonFiniteResult
                            else
                                if norm == 0 then Left (DomainError (name ++ ": zero iterate")) else go (k + 1) (map (/ norm) y)
    matrixScale :: a
matrixScale = [a] -> a
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (a
0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: ([a] -> a) -> [[a]] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map ([a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([a] -> a) -> ([a] -> [a]) -> [a] -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map a -> a
forall a. Num a => a -> a
abs) (Matrix a -> [[a]]
forall a. Matrix a -> [[a]]
entries Matrix a
a))

{- | Given positive finite tolerance,positive budget,square finite matrix,then nonzero
starting vector, estimate a dominant eigenpair. Return Iteration with
value/vector/residual/count; invalid dimensions/options, non-finite arithmetic or
exhausted budget fails.
-}
powerIteration :: (RealFloat a) => a -> Int -> Matrix a -> [a] -> Either MathError (Iteration a)
powerIteration :: forall a.
RealFloat a =>
a -> Int -> Matrix a -> [a] -> Either MathError (Iteration a)
powerIteration a
tol Int
budget Matrix a
a = String
-> a
-> Int
-> (a -> [a] -> Either MathError [a])
-> Matrix a
-> [a]
-> Either MathError (Iteration a)
forall a.
RealFloat a =>
String
-> a
-> Int
-> (a -> [a] -> Either MathError [a])
-> Matrix a
-> [a]
-> Either MathError (Iteration a)
iterateEigen String
"powerIteration" a
tol Int
budget (\a
_ [a]
x -> [a] -> Either MathError [a]
forall a b. b -> Either a b
Right [[a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
dotV [a]
row [a]
x | [a]
row <- Matrix a -> [[a]]
forall a. Matrix a -> [[a]]
entries Matrix a
a]) Matrix a
a

{- | Given tolerance,budget,shift,square matrix,then starting vector, estimate an eigenpair
near shift by repeated shifted solves. Require valid finite inputs and nonzero start;
singular shifted systems or convergence/numerical failures return Left MathError.
-}
inverseIteration ::
    (RealFloat a) => a -> Int -> a -> Matrix a -> [a] -> Either MathError (Iteration a)
inverseIteration :: forall a.
RealFloat a =>
a -> Int -> a -> Matrix a -> [a] -> Either MathError (Iteration a)
inverseIteration a
tol Int
budget a
shift Matrix a
a [a]
initial = do
    shifted <- Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
forall a.
Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
build Int
n Int
n (\Int
i Int
j -> Matrix a -> Int -> Int -> a
forall a. Matrix a -> Int -> Int -> a
get Matrix a
a Int
i Int
j a -> a -> a
forall a. Num a => a -> a -> a
- if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then a
shift else a
0)
    iterateEigen "inverseIteration" tol budget (\a
_ -> Matrix a -> [a] -> Either MathError [a]
forall a. RealFloat a => Matrix a -> [a] -> Either MathError [a]
solve Matrix a
shifted) a initial
  where
    n :: Int
n = (Int, Int) -> Int
forall a b. (a, b) -> a
fst (Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a)

{- | Given tolerance,budget,square matrix,then nonzero initial vector, refine an eigenpair
with Rayleigh-quotient shifts. Shape/options, singular shifted solves or
convergence/numerical failures return Left MathError.
-}
rayleighIteration :: (RealFloat a) => a -> Int -> Matrix a -> [a] -> Either MathError (Iteration a)
rayleighIteration :: forall a.
RealFloat a =>
a -> Int -> Matrix a -> [a] -> Either MathError (Iteration a)
rayleighIteration a
tol Int
budget Matrix a
a = String
-> a
-> Int
-> (a -> [a] -> Either MathError [a])
-> Matrix a
-> [a]
-> Either MathError (Iteration a)
forall a.
RealFloat a =>
String
-> a
-> Int
-> (a -> [a] -> Either MathError [a])
-> Matrix a
-> [a]
-> Either MathError (Iteration a)
iterateEigen String
"rayleighIteration" a
tol Int
budget a -> [a] -> Either MathError [a]
action Matrix a
a
  where
    n :: Int
n = (Int, Int) -> Int
forall a b. (a, b) -> a
fst (Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a)
    action :: a -> [a] -> Either MathError [a]
action a
shift [a]
x = do
        shifted <- Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
forall a.
Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
build Int
n Int
n (\Int
i Int
j -> Matrix a -> Int -> Int -> a
forall a. Matrix a -> Int -> Int -> a
get Matrix a
a Int
i Int
j a -> a -> a
forall a. Num a => a -> a -> a
- if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then a
shift else a
0)
        case luDecomposeWithTolerance 0 shifted >>= (`solveLU` x) of
            Left (SingularMatrix Int
_) -> do
                -- The Rayleigh quotient can round to an exact eigenvalue
                -- before the vector residual reaches tolerance. Perturb the
                -- singular shift by one representable scale and continue.
                let delta :: a
delta = a
2 a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. RealFloat a => a -> a
epsilon a
shift a -> a -> a
forall a. Num a => a -> a -> a
* a -> a -> a
forall a. Ord a => a -> a -> a
max a
1 (a -> a
forall a. Num a => a -> a
abs a
shift)
                regular <- Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
forall a.
Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
build Int
n Int
n (\Int
i Int
j -> Matrix a -> Int -> Int -> a
forall a. Matrix a -> Int -> Int -> a
get Matrix a
a Int
i Int
j a -> a -> a
forall a. Num a => a -> a -> a
- if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then a
shift a -> a -> a
forall a. Num a => a -> a -> a
+ a
delta else a
0)
                luDecomposeWithTolerance 0 regular >>= (`solveLU` x)
            Either MathError [a]
result -> Either MathError [a]
result

{- | Krylov basis columns and upper-Hessenberg projection returned by arnoldi; exact
breakdown may shorten the basis.
-}
data Arnoldi a = Arnoldi
    { forall a. Arnoldi a -> Matrix a
arnoldiBasis :: Matrix a
    -- ^ Orthonormal Krylov basis vectors as matrix columns.
    , forall a. Arnoldi a -> Matrix a
arnoldiHessenberg :: Matrix a
    -- ^ Hessenberg projection relating the matrix action to the Arnoldi basis.
    }
    deriving (Arnoldi a -> Arnoldi a -> Bool
(Arnoldi a -> Arnoldi a -> Bool)
-> (Arnoldi a -> Arnoldi a -> Bool) -> Eq (Arnoldi a)
forall a. Eq a => Arnoldi a -> Arnoldi a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall a. Eq a => Arnoldi a -> Arnoldi a -> Bool
== :: Arnoldi a -> Arnoldi a -> Bool
$c/= :: forall a. Eq a => Arnoldi a -> Arnoldi a -> Bool
/= :: Arnoldi a -> Arnoldi a -> Bool
Eq, Int -> Arnoldi a -> ShowS
[Arnoldi a] -> ShowS
Arnoldi a -> String
(Int -> Arnoldi a -> ShowS)
-> (Arnoldi a -> String)
-> ([Arnoldi a] -> ShowS)
-> Show (Arnoldi a)
forall a. Show a => Int -> Arnoldi a -> ShowS
forall a. Show a => [Arnoldi a] -> ShowS
forall a. Show a => Arnoldi a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> Arnoldi a -> ShowS
showsPrec :: Int -> Arnoldi a -> ShowS
$cshow :: forall a. Show a => Arnoldi a -> String
show :: Arnoldi a -> String
$cshowList :: forall a. Show a => [Arnoldi a] -> ShowS
showList :: [Arnoldi a] -> ShowS
Show)

{- | Given positive finite tolerance,requested basis size,square matrix,then nonzero initial
vector, build Arnoldi basis/projection. Require 1\<=requested\<=matrix order. Invalid
input or numerical failures return Left MathError; exact Krylov breakdown terminates
the basis.
-}
arnoldi :: (RealFloat a) => a -> Int -> Matrix a -> [a] -> Either MathError (Arnoldi a)
arnoldi :: forall a.
RealFloat a =>
a -> Int -> Matrix a -> [a] -> Either MathError (Arnoldi a)
arnoldi a
tol Int
requested Matrix a
a [a]
initial
    | Int
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n = MathError -> Either MathError (Arnoldi a)
forall a b. a -> Either a b
Left ((Int, Int) -> MathError
NonSquareMatrix (Int
m, Int
n))
    | [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
initial Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n = MathError -> Either MathError (Arnoldi a)
forall a b. a -> Either a b
Left (Int -> Int -> MathError
VectorLengthMismatch Int
n ([a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
initial))
    | Int
requested Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
1 Bool -> Bool -> Bool
|| Int
requested Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
n =
        MathError -> Either MathError (Arnoldi a)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"arnoldi: requested basis must be between one and matrix order")
    | a
tol a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0 Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
tol Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
tol = MathError -> Either MathError (Arnoldi a)
forall a b. a -> Either a b
Left MathError
InvalidTolerance
    | Bool
otherwise = do
        Matrix a -> Either MathError ()
forall a. RealFloat a => Matrix a -> Either MathError ()
finiteMatrix Matrix a
a
        [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues [a]
initial
        if a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
beta Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
beta
            then MathError -> Either MathError (Arnoldi a)
forall a b. a -> Either a b
Left MathError
NonFiniteResult
            else
                if a
beta a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0
                    then MathError -> Either MathError (Arnoldi a)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"arnoldi: zero initial vector")
                    else [[a]] -> [[a]] -> Int -> Either MathError (Arnoldi a)
go [(a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
beta) [a]
initial] [] Int
0
  where
    (Int
m, Int
n) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a
    beta :: a
beta = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
initial
    go :: [[a]] -> [[a]] -> Int -> Either MathError (Arnoldi a)
go [[a]]
qs [[a]]
hs Int
k = do
        let w0 :: [a]
w0 = [[a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
dotV [a]
row ([[a]]
qs [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
k) | [a]
row <- Matrix a -> [[a]]
forall a. Matrix a -> [[a]]
entries Matrix a
a]
            orth :: ([a], [a]) -> [a] -> ([a], [a])
orth ([a]
w, [a]
coeffs) [a]
q = let h :: a
h = [a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
dotV [a]
q [a]
w in ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\a
x a
y -> a
x a -> a -> a
forall a. Num a => a -> a -> a
- a
h a -> a -> a
forall a. Num a => a -> a -> a
* a
y) [a]
w [a]
q, [a]
coeffs [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a
h])
            ([a]
w1, [a]
h1) = (([a], [a]) -> [a] -> ([a], [a]))
-> ([a], [a]) -> [[a]] -> ([a], [a])
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl ([a], [a]) -> [a] -> ([a], [a])
forall {a}. Num a => ([a], [a]) -> [a] -> ([a], [a])
orth ([a]
w0, []) [[a]]
qs
            ([a]
w2, [a]
h2) = (([a], [a]) -> [a] -> ([a], [a]))
-> ([a], [a]) -> [[a]] -> ([a], [a])
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl ([a], [a]) -> [a] -> ([a], [a])
forall {a}. Num a => ([a], [a]) -> [a] -> ([a], [a])
orth ([a]
w1, []) [[a]]
qs
            coefficients :: [a]
coefficients = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(+) [a]
h1 [a]
h2
            norm :: a
norm = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
w2
            column :: [a]
column = [a]
coefficients [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a
norm]
            done :: Bool
done = a
norm a -> a -> a
forall a. Fractional a => a -> a -> a
/ a -> a -> a
forall a. Ord a => a -> a -> a
max a
1 ([a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
w0) a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
tol Bool -> Bool -> Bool
|| Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
requested
            nextQs :: [[a]]
nextQs = if a
norm a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 then [[a]]
qs else [[a]]
qs [[a]] -> [[a]] -> [[a]]
forall a. [a] -> [a] -> [a]
++ [(a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
norm) [a]
w2]
        [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues ([a]
column [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a]
w2 [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [[a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
w0])
        if Bool
done
            then do
                let count :: Int
count = Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1; rows :: Int
rows = [[a]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[a]]
nextQs
                qmat <- Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
forall a.
Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
build Int
n Int
rows (\Int
i Int
j -> [[a]]
nextQs [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
i)
                hmat <-
                    build
                        rows
                        count
                        (\Int
i Int
j -> if Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (([[a]]
hs [[a]] -> [[a]] -> [[a]]
forall a. [a] -> [a] -> [a]
++ [[a]
column]) [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) then ([[a]]
hs [[a]] -> [[a]] -> [[a]]
forall a. [a] -> [a] -> [a]
++ [[a]
column]) [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
i else a
0)
                Right (Arnoldi qmat hmat)
            else [[a]] -> [[a]] -> Int -> Either MathError (Arnoldi a)
go [[a]]
nextQs ([[a]]
hs [[a]] -> [[a]] -> [[a]]
forall a. [a] -> [a] -> [a]
++ [[a]
column]) (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)

{- | Symmetric Lanczos with the three-term recurrence and two full
reorthogonalization passes. Small correction coefficients are retained in H
so the projection residual remains meaningful after floating-point rounding.
-}
lanczos :: (RealFloat a) => a -> Int -> Matrix a -> [a] -> Either MathError (Arnoldi a)
lanczos :: forall a.
RealFloat a =>
a -> Int -> Matrix a -> [a] -> Either MathError (Arnoldi a)
lanczos a
tol Int
requested Matrix a
a [a]
initial = do
    _ <- Matrix a -> Either MathError Int
forall a. RealFloat a => Matrix a -> Either MathError Int
checkSymmetric Matrix a
a
    if length initial /= n
        then Left (VectorLengthMismatch n (length initial))
        else
            if requested < 1 || requested > n
                then Left (InvalidArgument "lanczos: invalid basis size")
                else
                    if tol < 0 || isNaN tol || isInfinite tol
                        then Left InvalidTolerance
                        else do
                            finiteValues initial
                            let beta = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
initial
                            if isNaN beta || isInfinite beta
                                then Left NonFiniteResult
                                else
                                    if beta == 0
                                        then Left (DomainError "lanczos: zero initial vector")
                                        else go [map (/ beta) initial] [] 0 0
  where
    n :: Int
n = (Int, Int) -> Int
forall a b. (a, b) -> a
fst (Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a)
    go :: [[a]] -> [[a]] -> a -> Int -> Either MathError (Arnoldi a)
go [[a]]
qs [[a]]
columns a
previousBeta Int
k = do
        let q :: [a]
q = [[a]]
qs [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
k
            aq :: [a]
aq = [[a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
dotV [a]
row [a]
q | [a]
row <- Matrix a -> [[a]]
forall a. Matrix a -> [[a]]
entries Matrix a
a]
            alpha :: a
alpha = [a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
dotV [a]
q [a]
aq
            previous :: [a]
previous = if Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then Int -> a -> [a]
forall a. Int -> a -> [a]
replicate Int
n a
0 else [[a]]
qs [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
            residual :: [a]
residual = (a -> a -> a -> a) -> [a] -> [a] -> [a] -> [a]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 (\a
x a
y a
z -> a
x a -> a -> a
forall a. Num a => a -> a -> a
- a
alpha a -> a -> a
forall a. Num a => a -> a -> a
* a
y a -> a -> a
forall a. Num a => a -> a -> a
- a
previousBeta a -> a -> a
forall a. Num a => a -> a -> a
* a
z) [a]
aq [a]
q [a]
previous
            orth :: ([a], [a]) -> [a] -> ([a], [a])
orth ([a]
w, [a]
coeffs) [a]
basis = let h :: a
h = [a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
dotV [a]
basis [a]
w in ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\a
x a
y -> a
x a -> a -> a
forall a. Num a => a -> a -> a
- a
h a -> a -> a
forall a. Num a => a -> a -> a
* a
y) [a]
w [a]
basis, [a]
coeffs [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a
h])
            ([a]
w1, [a]
h1) = (([a], [a]) -> [a] -> ([a], [a]))
-> ([a], [a]) -> [[a]] -> ([a], [a])
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl ([a], [a]) -> [a] -> ([a], [a])
forall {a}. Num a => ([a], [a]) -> [a] -> ([a], [a])
orth ([a]
residual, []) [[a]]
qs
            ([a]
w2, [a]
h2) = (([a], [a]) -> [a] -> ([a], [a]))
-> ([a], [a]) -> [[a]] -> ([a], [a])
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl ([a], [a]) -> [a] -> ([a], [a])
forall {a}. Num a => ([a], [a]) -> [a] -> ([a], [a])
orth ([a]
w1, []) [[a]]
qs
            beta :: a
beta = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
w2
            base :: [a]
base = [if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
k then a
alpha else if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1 then a
previousBeta else a
0 | Int
i <- [Int
0 .. Int
k]]
            column :: [a]
column = (a -> a -> a -> a) -> [a] -> [a] -> [a] -> [a]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 (\a
x a
y a
z -> a
x a -> a -> a
forall a. Num a => a -> a -> a
+ a
y a -> a -> a
forall a. Num a => a -> a -> a
+ a
z) [a]
base [a]
h1 [a]
h2 [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a
beta]
            next :: [[a]]
next = if a
beta a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 then [[a]]
qs else [[a]]
qs [[a]] -> [[a]] -> [[a]]
forall a. [a] -> [a] -> [a]
++ [(a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
beta) [a]
w2]
        [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues ([a]
column [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a]
w2 [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [[a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
aq])
        if a
beta a -> a -> a
forall a. Fractional a => a -> a -> a
/ a -> a -> a
forall a. Ord a => a -> a -> a
max a
1 ([a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
aq) a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
tol Bool -> Bool -> Bool
|| Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
requested
            then do
                qmat <- Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
forall a.
Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
build Int
n ([[a]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[a]]
next) (\Int
i Int
j -> [[a]]
next [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
i)
                hmat <-
                    build
                        (length next)
                        (k + 1)
                        (\Int
i Int
j -> if Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (([[a]]
columns [[a]] -> [[a]] -> [[a]]
forall a. [a] -> [a] -> [a]
++ [[a]
column]) [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) then ([[a]]
columns [[a]] -> [[a]] -> [[a]]
forall a. [a] -> [a] -> [a]
++ [[a]
column]) [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
i else a
0)
                Right (Arnoldi qmat hmat)
            else [[a]] -> [[a]] -> a -> Int -> Either MathError (Arnoldi a)
go [[a]]
next ([[a]]
columns [[a]] -> [[a]] -> [[a]]
forall a. [a] -> [a] -> [a]
++ [[a]
column]) a
beta (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)

{- | Given positive finite tolerance,positive budget,then finite square matrix, iterate QR
to return @(Q,T)@: accumulated orthogonal vectors and the resulting iterate. Invalid
options/shape or failure to converge within budget returns Left MathError.
-}
qrIteration :: (RealFloat a) => a -> Int -> Matrix a -> Either MathError (Matrix a, Matrix a)
qrIteration :: forall a.
RealFloat a =>
a -> Int -> Matrix a -> Either MathError (Matrix a, Matrix a)
qrIteration a
tol Int
budget Matrix a
a
    | Int
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n = MathError -> Either MathError (Matrix a, Matrix a)
forall a b. a -> Either a b
Left ((Int, Int) -> MathError
NonSquareMatrix (Int
m, Int
n))
    | a
tol a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
0 Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
tol Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
tol = MathError -> Either MathError (Matrix a, Matrix a)
forall a b. a -> Either a b
Left MathError
InvalidTolerance
    | Int
budget Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
1 = MathError -> Either MathError (Matrix a, Matrix a)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"qrIteration: positive budget required")
    | Bool
otherwise = do
        Matrix a -> Either MathError ()
forall a. RealFloat a => Matrix a -> Either MathError ()
finiteMatrix Matrix a
a; q <- Int -> Either MathError (Matrix a)
forall a. Num a => Int -> Either MathError (Matrix a)
identityChecked Int
n; go 0 q a
  where
    (Int
m, Int
n) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a
    go :: Int
-> Matrix a -> Matrix a -> Either MathError (Matrix a, Matrix a)
go Int
k Matrix a
q Matrix a
t
        | (a -> Bool) -> [a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (\a
x -> a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
x Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
x) [a
offNorm, a
totalNorm] = MathError -> Either MathError (Matrix a, Matrix a)
forall a b. a -> Either a b
Left MathError
NonFiniteResult
        | a
offNorm a -> a -> a
forall a. Fractional a => a -> a -> a
/ a -> a -> a
forall a. Ord a => a -> a -> a
max a
1 a
totalNorm a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
tol = (Matrix a, Matrix a) -> Either MathError (Matrix a, Matrix a)
forall a b. b -> Either a b
Right (Matrix a
q, Matrix a
t)
        | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
budget = MathError -> Either MathError (Matrix a, Matrix a)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"qrIteration" Int
budget)
        | Bool
otherwise = do
            -- Unpivoted Givens QR is required for a similarity iteration.
            (u, r) <- Matrix a -> Either MathError (Matrix a, Matrix a)
forall a.
RealFloat a =>
Matrix a -> Either MathError (Matrix a, Matrix a)
givensQR Matrix a
t
            next <- multiplyChecked r u
            qs <- multiplyChecked q u
            finiteValues (toList qs ++ toList next)
            go (k + 1) qs next
      where
        offNorm :: a
offNorm = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [Matrix a -> Int -> Int -> a
forall a. Matrix a -> Int -> Int -> a
get Matrix a
t Int
i Int
j | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
j <- [Int
0 .. Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
        totalNorm :: a
totalNorm = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm (Matrix a -> [a]
forall a. Matrix a -> [a]
toList Matrix a
t)

{- | Verify supplied Jordan chains against A and lambda and construct J and P.
Computing a Jordan structure from floating input is ill posed; the caller
supplies chain vectors ordered eigenvector first. A P = P J is checked.
This API never pretends a Schur form is a Jordan form.
-}
jordanChains ::
    (RealFloat a) => a -> Matrix a -> [(a, [[a]])] -> Either MathError (Matrix a, Matrix a)
jordanChains :: forall a.
RealFloat a =>
a
-> Matrix a
-> [(a, [[a]])]
-> Either MathError (Matrix a, Matrix a)
jordanChains a
tol Matrix a
a [(a, [[a]])]
chains
    | Int
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n = MathError -> Either MathError (Matrix a, Matrix a)
forall a b. a -> Either a b
Left ((Int, Int) -> MathError
NonSquareMatrix (Int
m, Int
n))
    | a
tol a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0 Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
tol Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
tol = MathError -> Either MathError (Matrix a, Matrix a)
forall a b. a -> Either a b
Left MathError
InvalidTolerance
    | [[a]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[a]]
vectors Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n Bool -> Bool -> Bool
|| ([a] -> Bool) -> [[a]] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any ((Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n) (Int -> Bool) -> ([a] -> Int) -> [a] -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length) [[a]]
vectors Bool -> Bool -> Bool
|| ((a, [[a]]) -> Bool) -> [(a, [[a]])] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any ([[a]] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null ([[a]] -> Bool) -> ((a, [[a]]) -> [[a]]) -> (a, [[a]]) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a, [[a]]) -> [[a]]
forall a b. (a, b) -> b
snd) [(a, [[a]])]
chains =
        MathError -> Either MathError (Matrix a, Matrix a)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"jordanChains: chains must form a square basis")
    | Bool
otherwise = do
        Matrix a -> Either MathError ()
forall a. RealFloat a => Matrix a -> Either MathError ()
finiteMatrix Matrix a
a
        [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues ([a]
values [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [[a]] -> [a]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[a]]
vectors)
        p <- Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
forall a.
Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
build Int
n Int
n (\Int
i Int
j -> [[a]]
vectors [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
i)
        _ <- inverse p
        j <-
            build
                n
                n
                (\Int
i Int
col -> if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
col then [a]
values [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
i else if Int
col Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 Bool -> Bool -> Bool
&& Bool -> Bool
not (Int
col Int -> [Int] -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [Int]
starts) then a
1 else a
0)
        ap <- multiplyChecked a p
        pj <- multiplyChecked p j
        delta <- subtractMatrixChecked ap pj
        let deltaNorm = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm (Matrix a -> [a]
forall a. Matrix a -> [a]
toList Matrix a
delta); apNorm = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm (Matrix a -> [a]
forall a. Matrix a -> [a]
toList Matrix a
ap)
        finiteValues (toList p ++ toList j ++ toList ap ++ toList pj ++ toList delta ++ [deltaNorm, apNorm])
        if deltaNorm / max 1 apNorm > tol
            then Left (DomainError "jordanChains: chain residual exceeds tolerance")
            else Right (p, j)
  where
    (Int
m, Int
n) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a
    vectors :: [[a]]
vectors = ((a, [[a]]) -> [[a]]) -> [(a, [[a]])] -> [[a]]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (a, [[a]]) -> [[a]]
forall a b. (a, b) -> b
snd [(a, [[a]])]
chains
    values :: [a]
values = [[a]] -> [a]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [Int -> a -> [a]
forall a. Int -> a -> [a]
replicate ([[a]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[a]]
vs) a
value | (a
value, [[a]]
vs) <- [(a, [[a]])]
chains]
    starts :: [Int]
starts = (Int -> Int -> Int) -> Int -> [Int] -> [Int]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl Int -> Int -> Int
forall a. Num a => a -> a -> a
(+) Int
0 (((a, [[a]]) -> Int) -> [(a, [[a]])] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map ([[a]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([[a]] -> Int) -> ((a, [[a]]) -> [[a]]) -> (a, [[a]]) -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a, [[a]]) -> [[a]]
forall a b. (a, b) -> b
snd) [(a, [[a]])]
chains)

{- | Exact Jordan construction over Rational with a caller-provided complete
set of distinct rational eigenvalues (multiplicities are inferred). Unlike
floating Jordan algorithms this never uses a tolerance to decide structure.
Nullspaces of successive (A-lambda I) powers construct generalized chains.
An incomplete/nonrational spectrum is rejected rather than approximated.
-}
jordanForm :: Matrix Rational -> [Rational] -> Either MathError (Matrix Rational, Matrix Rational)
jordanForm :: Matrix Rational
-> [Rational]
-> Either MathError (Matrix Rational, Matrix Rational)
jordanForm Matrix Rational
a [Rational]
spectrum
    | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
m = MathError -> Either MathError (Matrix Rational, Matrix Rational)
forall a b. a -> Either a b
Left ((Int, Int) -> MathError
NonSquareMatrix (Int
n, Int
m))
    | [Rational] -> [Rational]
forall a. Eq a => [a] -> [a]
nub [Rational]
spectrum [Rational] -> [Rational] -> Bool
forall a. Eq a => a -> a -> Bool
/= [Rational]
spectrum =
        MathError -> Either MathError (Matrix Rational, Matrix Rational)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"jordanForm: provide distinct rational eigenvalues")
    | Bool
otherwise = do
        (vectors, blocks) <- (([[Rational]], [(Rational, Int)])
 -> Rational -> Either MathError ([[Rational]], [(Rational, Int)]))
-> ([[Rational]], [(Rational, Int)])
-> [Rational]
-> Either MathError ([[Rational]], [(Rational, Int)])
forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM ([[Rational]], [(Rational, Int)])
-> Rational -> Either MathError ([[Rational]], [(Rational, Int)])
eigenChains ([], []) [Rational]
spectrum
        if length vectors /= n
            then Left (DomainError "jordanForm: supplied rational spectrum is incomplete")
            else do
                p <- build n n (\Int
i Int
j -> [[Rational]]
vectors [[Rational]] -> Int -> [Rational]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j [Rational] -> Int -> Rational
forall a. HasCallStack => [a] -> Int -> a
!! Int
i)
                let values = [[Rational]] -> [Rational]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [Int -> Rational -> [Rational]
forall a. Int -> a -> [a]
replicate Int
size Rational
value | (Rational
value, Int
size) <- [(Rational, Int)]
blocks]
                    starts = (Int -> Int -> Int) -> Int -> [Int] -> [Int]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl Int -> Int -> Int
forall a. Num a => a -> a -> a
(+) Int
0 (((Rational, Int) -> Int) -> [(Rational, Int)] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map (Rational, Int) -> Int
forall a b. (a, b) -> b
snd [(Rational, Int)]
blocks)
                j <-
                    build
                        n
                        n
                        (\Int
i Int
col -> if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
col then [Rational]
values [Rational] -> Int -> Rational
forall a. HasCallStack => [a] -> Int -> a
!! Int
i else if Int
col Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 Bool -> Bool -> Bool
&& Int
col Int -> [Int] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [Int]
starts then Rational
1 else Rational
0)
                ap <- multiplyChecked a p
                pj <- multiplyChecked p j
                if ap /= pj || exactRank p /= n
                    then Left (DomainError "jordanForm: exact chain construction failed")
                    else Right (p, j)
  where
    (Int
n, Int
m) = Matrix Rational -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix Rational
a
    eigenChains :: ([[Rational]], [(Rational, Int)])
-> Rational -> Either MathError ([[Rational]], [(Rational, Int)])
eigenChains ([[Rational]]
basis, [(Rational, Int)]
blocks) Rational
lambda = do
        shifted <- Int
-> Int
-> (Int -> Int -> Rational)
-> Either MathError (Matrix Rational)
forall a.
Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
build Int
n Int
n (\Int
i Int
j -> Matrix Rational -> Int -> Int -> Rational
forall a. Matrix a -> Int -> Int -> a
get Matrix Rational
a Int
i Int
j Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
- if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then Rational
lambda else Rational
0)
        powers <- scanPowers shifted
        foldM (choose shifted powers lambda) (basis, blocks) (reverse [1 .. n])
    scanPowers :: Matrix a -> Either MathError [Matrix a]
scanPowers Matrix a
shifted = do
        i <- Int -> Either MathError (Matrix a)
forall a. Num a => Int -> Either MathError (Matrix a)
identityChecked Int
n
        foldM (\[Matrix a]
xs Int
_ -> do next <- Matrix a -> Matrix a -> Either MathError (Matrix a)
forall a.
Num a =>
Matrix a -> Matrix a -> Either MathError (Matrix a)
multiplyChecked ([Matrix a] -> Matrix a
forall a. HasCallStack => [a] -> a
last [Matrix a]
xs) Matrix a
shifted; Right (xs ++ [next])) [i] [1 .. n]
    choose :: Matrix Rational
-> [Matrix Rational]
-> a
-> ([[Rational]], [(a, Int)])
-> Int
-> Either MathError ([[Rational]], [(a, Int)])
choose Matrix Rational
shifted [Matrix Rational]
powers a
lambda ([[Rational]], [(a, Int)])
state Int
k = (([[Rational]], [(a, Int)])
 -> [Rational] -> Either MathError ([[Rational]], [(a, Int)]))
-> ([[Rational]], [(a, Int)])
-> [[Rational]]
-> Either MathError ([[Rational]], [(a, Int)])
forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM (Matrix Rational
-> a
-> Int
-> ([[Rational]], [(a, Int)])
-> [Rational]
-> Either MathError ([[Rational]], [(a, Int)])
forall {a}.
Matrix Rational
-> a
-> Int
-> ([[Rational]], [(a, Int)])
-> [Rational]
-> Either MathError ([[Rational]], [(a, Int)])
candidate Matrix Rational
shifted a
lambda Int
k) ([[Rational]], [(a, Int)])
state (Matrix Rational -> [[Rational]]
exactKernel ([Matrix Rational]
powers [Matrix Rational] -> Int -> Matrix Rational
forall a. HasCallStack => [a] -> Int -> a
!! Int
k))
    candidate :: Matrix Rational
-> a
-> Int
-> ([[Rational]], [(a, Int)])
-> [Rational]
-> Either MathError ([[Rational]], [(a, Int)])
candidate Matrix Rational
shifted a
lambda Int
k ([[Rational]]
basis, [(a, Int)]
blocks) [Rational]
top = do
        let descending :: [[Rational]]
descending = Int -> [[Rational]] -> [[Rational]]
forall a. Int -> [a] -> [a]
take Int
k (([Rational] -> [Rational]) -> [Rational] -> [[Rational]]
forall a. (a -> a) -> a -> [a]
iterate (\[Rational]
x -> [[Rational] -> [Rational] -> Rational
forall a. Num a => [a] -> [a] -> a
dotV [Rational]
row [Rational]
x | [Rational]
row <- Matrix Rational -> [[Rational]]
forall a. Matrix a -> [[a]]
entries Matrix Rational
shifted]) [Rational]
top)
            chain :: [[Rational]]
chain = [[Rational]] -> [[Rational]]
forall a. [a] -> [a]
reverse [[Rational]]
descending
        combined <- Int
-> Int
-> (Int -> Int -> Rational)
-> Either MathError (Matrix Rational)
forall a.
Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
build Int
n ([[Rational]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Rational]]
basis Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
k) (\Int
i Int
j -> ([[Rational]]
basis [[Rational]] -> [[Rational]] -> [[Rational]]
forall a. [a] -> [a] -> [a]
++ [[Rational]]
chain) [[Rational]] -> Int -> [Rational]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j [Rational] -> Int -> Rational
forall a. HasCallStack => [a] -> Int -> a
!! Int
i)
        if exactRank combined == length basis + k
            then Right (basis ++ chain, blocks ++ [(lambda, k)])
            else Right (basis, blocks)

exactRref :: Matrix Rational -> ([[Rational]], [Int])
exactRref :: Matrix Rational -> ([[Rational]], [Int])
exactRref Matrix Rational
a = Int -> Int -> [[Rational]] -> [Int] -> ([[Rational]], [Int])
forall {a}.
(Eq a, Fractional a) =>
Int -> Int -> [[a]] -> [Int] -> ([[a]], [Int])
go Int
0 Int
0 (Matrix Rational -> [[Rational]]
forall a. Matrix a -> [[a]]
entries Matrix Rational
a) []
  where
    (Int
r, Int
c) = Matrix Rational -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix Rational
a
    go :: Int -> Int -> [[a]] -> [Int] -> ([[a]], [Int])
go Int
i Int
j [[a]]
rows [Int]
pivots
        | Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
r Bool -> Bool -> Bool
|| Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
c = ([[a]]
rows, [Int]
pivots)
        | Bool
otherwise = case ([a] -> Bool) -> [[a]] -> Maybe Int
forall a. (a -> Bool) -> [a] -> Maybe Int
findIndex (\[a]
row -> [a]
row [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
j a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
0) (Int -> [[a]] -> [[a]]
forall a. Int -> [a] -> [a]
drop Int
i [[a]]
rows) of
            Maybe Int
Nothing -> Int -> Int -> [[a]] -> [Int] -> ([[a]], [Int])
go Int
i (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [[a]]
rows [Int]
pivots
            Just Int
offset ->
                let swapped :: [[a]]
swapped = Int -> Int -> [[a]] -> [[a]]
forall a. Int -> Int -> [a] -> [a]
swap Int
i (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
offset) [[a]]
rows
                    normalized :: [a]
normalized = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
forall a. Fractional a => a -> a -> a
/ ([[a]]
swapped [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
j)) ([[a]]
swapped [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i)
                    next :: [[a]]
next =
                        (Int -> [a] -> [a]) -> [Int] -> [[a]] -> [[a]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith
                            (\Int
k [a]
row -> if Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i then [a]
normalized else (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [a]
row ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (([a]
row [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) a -> a -> a
forall a. Num a => a -> a -> a
*) [a]
normalized))
                            [Int
0 ..]
                            [[a]]
swapped
                 in Int -> Int -> [[a]] -> [Int] -> ([[a]], [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) [[a]]
next ([Int]
pivots [Int] -> [Int] -> [Int]
forall a. [a] -> [a] -> [a]
++ [Int
j])
exactRank :: Matrix Rational -> Int
exactRank :: Matrix Rational -> Int
exactRank = [Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([Int] -> Int)
-> (Matrix Rational -> [Int]) -> Matrix Rational -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([[Rational]], [Int]) -> [Int]
forall a b. (a, b) -> b
snd (([[Rational]], [Int]) -> [Int])
-> (Matrix Rational -> ([[Rational]], [Int]))
-> Matrix Rational
-> [Int]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix Rational -> ([[Rational]], [Int])
exactRref
exactKernel :: Matrix Rational -> [[Rational]]
exactKernel :: Matrix Rational -> [[Rational]]
exactKernel Matrix Rational
a =
    [ [ if Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
free
            then Rational
1
            else case Int -> [(Int, Int)] -> Maybe Int
forall a b. Eq a => a -> [(a, b)] -> Maybe b
lookup Int
j ([Int] -> [Int] -> [(Int, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int]
pivots [Int
0 ..]) of Just Int
row -> Rational -> Rational
forall a. Num a => a -> a
negate ([[Rational]]
rows [[Rational]] -> Int -> [Rational]
forall a. HasCallStack => [a] -> Int -> a
!! Int
row [Rational] -> Int -> Rational
forall a. HasCallStack => [a] -> Int -> a
!! Int
free); Maybe Int
Nothing -> Rational
0
      | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
      ]
    | Int
free <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
    , Int
free Int -> [Int] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [Int]
pivots
    ]
  where
    (Int
_, Int
n) = Matrix Rational -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix Rational
a; ([[Rational]]
rows, [Int]
pivots) = Matrix Rational -> ([[Rational]], [Int])
exactRref Matrix Rational
a