-- SPDX-License-Identifier: MPL-2.0

{- | Bounded reference iterative solvers for dense real matrices. Every returned
result is checked against the original residual. Square methods stop when
||b-Ax|| <= atol + rtol*||b||. Least-squares methods may also stop when
||A^T(b-Ax)|| <= atol + rtol*||A^T b||. Full Krylov bases are retained in
GMRES, MINRES and LSMR; these emphasize transparent mathematics over storage.
Nonfinite arithmetic, breakdown and exhausted budgets return structured errors.
-}
module HMath.Matrix.Iterative (
    IterativeOptions (..),
    defaultIterativeOptions,
    IterativeResult (..),
    Preconditioner,
    applyPreconditioner,
    identityPreconditioner,
    jacobiPreconditioner,
    blockJacobiPreconditioner,
    ilu0Preconditioner,
    ic0Preconditioner,
    conjugateGradient,
    preconditionedCG,
    gmres,
    minres,
    bicg,
    bicgstab,
    lsqr,
    lsmr,
    jacobiSolve,
    gaussSeidelSolve,
    sorSolve,
    iterativeRefinement,
) where

import Control.Monad (foldM)
import HMath.Error (MathError (..))
import HMath.Matrix (Matrix, fromListChecked, shape, transposeMatrix)
import HMath.Matrix.Internal (dotV, entries, finiteMatrix, finiteValues, stableNorm)
import HMath.Matrix.LU (luDecomposeWithTolerance, solveLU)
import HMath.Matrix.QR (leastSquares)

{- | Absolute and relative residual tolerances and maximum iterations. A zero
budget succeeds only when the supplied initial guess already meets tolerance.
-}
data IterativeOptions a = IterativeOptions
    { forall a. IterativeOptions a -> a
absoluteTolerance :: a
    -- ^ Finite nonnegative absolute residual tolerance.
    , forall a. IterativeOptions a -> a
relativeTolerance :: a
    -- ^ Finite nonnegative relative residual tolerance.
    , forall a. IterativeOptions a -> Int
maximumIterations :: Int
    -- ^ Nonnegative maximum iteration count.
    }
    deriving (IterativeOptions a -> IterativeOptions a -> Bool
(IterativeOptions a -> IterativeOptions a -> Bool)
-> (IterativeOptions a -> IterativeOptions a -> Bool)
-> Eq (IterativeOptions a)
forall a. Eq a => IterativeOptions a -> IterativeOptions a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall a. Eq a => IterativeOptions a -> IterativeOptions a -> Bool
== :: IterativeOptions a -> IterativeOptions a -> Bool
$c/= :: forall a. Eq a => IterativeOptions a -> IterativeOptions a -> Bool
/= :: IterativeOptions a -> IterativeOptions a -> Bool
Eq, Int -> IterativeOptions a -> ShowS
[IterativeOptions a] -> ShowS
IterativeOptions a -> String
(Int -> IterativeOptions a -> ShowS)
-> (IterativeOptions a -> String)
-> ([IterativeOptions a] -> ShowS)
-> Show (IterativeOptions a)
forall a. Show a => Int -> IterativeOptions a -> ShowS
forall a. Show a => [IterativeOptions a] -> ShowS
forall a. Show a => IterativeOptions a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> IterativeOptions a -> ShowS
showsPrec :: Int -> IterativeOptions a -> ShowS
$cshow :: forall a. Show a => IterativeOptions a -> String
show :: IterativeOptions a -> String
$cshowList :: forall a. Show a => [IterativeOptions a] -> ShowS
showList :: [IterativeOptions a] -> ShowS
Show)

-- | 1000 iterations, zero absolute tolerance and relative tolerance 1e-10.
defaultIterativeOptions :: (Fractional a) => IterativeOptions a
defaultIterativeOptions :: forall a. Fractional a => IterativeOptions a
defaultIterativeOptions = a -> a -> Int -> IterativeOptions a
forall a. a -> a -> Int -> IterativeOptions a
IterativeOptions a
0 a
1e-10 Int
1000

-- | Converged state with norms evaluated from the original matrix and rhs.
data IterativeResult a = IterativeResult
    { forall a. IterativeResult a -> [a]
iterativeSolution :: [a]
    -- ^ Solution in original variable order.
    , forall a. IterativeResult a -> Int
iterativeIterations :: Int
    -- ^ Number of completed updates.
    , forall a. IterativeResult a -> a
iterativeResidualNorm :: a
    -- ^ Euclidean norm of b-Ax.
    , forall a. IterativeResult a -> a
iterativeNormalResidualNorm :: a
    -- ^ Euclidean norm of A^T(b-Ax), useful for inconsistent least squares.
    }
    deriving (IterativeResult a -> IterativeResult a -> Bool
(IterativeResult a -> IterativeResult a -> Bool)
-> (IterativeResult a -> IterativeResult a -> Bool)
-> Eq (IterativeResult a)
forall a. Eq a => IterativeResult a -> IterativeResult a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall a. Eq a => IterativeResult a -> IterativeResult a -> Bool
== :: IterativeResult a -> IterativeResult a -> Bool
$c/= :: forall a. Eq a => IterativeResult a -> IterativeResult a -> Bool
/= :: IterativeResult a -> IterativeResult a -> Bool
Eq, Int -> IterativeResult a -> ShowS
[IterativeResult a] -> ShowS
IterativeResult a -> String
(Int -> IterativeResult a -> ShowS)
-> (IterativeResult a -> String)
-> ([IterativeResult a] -> ShowS)
-> Show (IterativeResult a)
forall a. Show a => Int -> IterativeResult a -> ShowS
forall a. Show a => [IterativeResult a] -> ShowS
forall a. Show a => IterativeResult a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> IterativeResult a -> ShowS
showsPrec :: Int -> IterativeResult a -> ShowS
$cshow :: forall a. Show a => IterativeResult a -> String
show :: IterativeResult a -> String
$cshowList :: forall a. Show a => [IterativeResult a] -> ShowS
showList :: [IterativeResult a] -> ShowS
Show)

-- | A checked inverse action M^-1; constructors retain their input dimension.
newtype Preconditioner a = Preconditioner ([a] -> Either MathError [a])

-- | Apply an inverse action; wrong vector lengths and nonfinite values fail.
applyPreconditioner :: Preconditioner a -> [a] -> Either MathError [a]
applyPreconditioner :: forall a. Preconditioner a -> [a] -> Either MathError [a]
applyPreconditioner (Preconditioner [a] -> Either MathError [a]
f) = [a] -> Either MathError [a]
f

-- | Identity inverse action in a nonnegative dimension.
identityPreconditioner :: (RealFloat a) => Int -> Either MathError (Preconditioner a)
identityPreconditioner :: forall a. RealFloat a => Int -> Either MathError (Preconditioner a)
identityPreconditioner Int
n
    | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = MathError -> Either MathError (Preconditioner a)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"negative preconditioner dimension")
    | Bool
otherwise = Preconditioner a -> Either MathError (Preconditioner a)
forall a b. b -> Either a b
Right (([a] -> Either MathError [a]) -> Preconditioner a
forall a. ([a] -> Either MathError [a]) -> Preconditioner a
Preconditioner (\[a]
xs -> Int -> [a] -> Either MathError ()
forall a. RealFloat a => Int -> [a] -> Either MathError ()
validateVector Int
n [a]
xs Either MathError () -> Either MathError [a] -> Either MathError [a]
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> [a] -> Either MathError [a]
forall a b. b -> Either a b
Right [a]
xs))

-- | Inverse diagonal preconditioner. Every diagonal entry must be nonzero.
jacobiPreconditioner :: (RealFloat a) => Matrix a -> Either MathError (Preconditioner a)
jacobiPreconditioner :: forall a.
RealFloat a =>
Matrix a -> Either MathError (Preconditioner a)
jacobiPreconditioner Matrix a
matrix = do
    n <- Matrix a -> Either MathError Int
forall a. RealFloat a => Matrix a -> Either MathError Int
squareFinite Matrix a
matrix
    let rows = Matrix a -> [[a]]
forall a. Matrix a -> [[a]]
entries Matrix a
matrix; diagonal = [[[a]] -> Int -> Int -> a
forall a. Num a => [[a]] -> Int -> Int -> a
value [[a]]
rows Int
i Int
i | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
    checkPivots diagonal
    pure
        ( Preconditioner
            ( \[a]
xs -> do
                Int -> [a] -> Either MathError ()
forall a. RealFloat a => Int -> [a] -> Either MathError ()
validateVector Int
n [a]
xs
                let scaled :: [a]
scaled = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Fractional a => a -> a -> a
(/) [a]
xs [a]
diagonal
                [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues [a]
scaled
                [a] -> Either MathError [a]
forall a. a -> Either MathError a
forall (f :: * -> *) a. Applicative f => a -> f a
pure [a]
scaled
            )
        )

{- | Block Jacobi from contiguous positive block sizes summing to the dimension.
Each diagonal block is factored with pivoted LU; off-block entries are ignored.
-}
blockJacobiPreconditioner ::
    (RealFloat a) => [Int] -> Matrix a -> Either MathError (Preconditioner a)
blockJacobiPreconditioner :: forall a.
RealFloat a =>
[Int] -> Matrix a -> Either MathError (Preconditioner a)
blockJacobiPreconditioner [Int]
sizes Matrix a
matrix = do
    n <- Matrix a -> Either MathError Int
forall a. RealFloat a => Matrix a -> Either MathError Int
squareFinite Matrix a
matrix
    if any (<= 0) sizes || sum (map toInteger sizes) /= toInteger n
        then Left (InvalidArgument "invalid block Jacobi partition")
        else do
            let offsets = (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 [Int]
sizes; rows = Matrix a -> [[a]]
forall a. Matrix a -> [[a]]
entries Matrix a
matrix
            factors <-
                mapM
                    ( \(Int
offset, Int
k) ->
                        Int -> Int -> [a] -> Either MathError (Matrix a)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
fromListChecked
                            Int
k
                            Int
k
                            [[[a]] -> Int -> Int -> a
forall a. Num a => [[a]] -> Int -> Int -> a
value [[a]]
rows (Int
offset Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
i) (Int
offset Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
j) | Int
i <- [Int
0 .. Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
j <- [Int
0 .. Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
                            Either MathError (Matrix a)
-> (Matrix a -> Either MathError (LU a)) -> Either MathError (LU a)
forall a b.
Either MathError a
-> (a -> Either MathError b) -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= a -> Matrix a -> Either MathError (LU a)
forall a. RealFloat a => a -> Matrix a -> Either MathError (LU a)
luDecomposeWithTolerance a
0
                    )
                    (zip offsets sizes)
            pure
                ( Preconditioner
                    ( \[a]
xs -> do
                        Int -> [a] -> Either MathError ()
forall a. RealFloat a => Int -> [a] -> Either MathError ()
validateVector Int
n [a]
xs
                        [[a]] -> [a]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat
                            ([[a]] -> [a]) -> Either MathError [[a]] -> Either MathError [a]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> [Either MathError [a]] -> Either MathError [[a]]
forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
forall (m :: * -> *) a. Monad m => [m a] -> m [a]
sequence
                                [ LU a -> [a] -> Either MathError [a]
forall a. RealFloat a => LU a -> [a] -> Either MathError [a]
solveLU LU a
factor (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
k (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
offset [a]
xs))
                                | (LU a
factor, Int
offset, Int
k) <- [LU a] -> [Int] -> [Int] -> [(LU a, Int, Int)]
forall a b c. [a] -> [b] -> [c] -> [(a, b, c)]
zip3 [LU a]
factors [Int]
offsets [Int]
sizes
                                ]
                    )
                )

{- | Incomplete LU with zero fill and no pivoting. The exact nonzero pattern of
A is retained, diagonal pivots are mandatory, and zero pivots fail explicitly.
-}
ilu0Preconditioner :: (RealFloat a) => Matrix a -> Either MathError (Preconditioner a)
ilu0Preconditioner :: forall a.
RealFloat a =>
Matrix a -> Either MathError (Preconditioner a)
ilu0Preconditioner Matrix a
matrix = do
    n <- Matrix a -> Either MathError Int
forall a. RealFloat a => Matrix a -> Either MathError Int
squareFinite Matrix a
matrix
    let original = Matrix a -> [[a]]
forall a. Matrix a -> [[a]]
entries Matrix a
matrix
        active Int
i Int
j = Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j Bool -> Bool -> Bool
|| [[a]] -> Int -> Int -> a
forall a. Num a => [[a]] -> Int -> Int -> a
value [[a]]
original Int
i Int
j a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
0
        update [[a]]
rows Int
k = do
            let pivot :: a
pivot = [[a]] -> Int -> Int -> a
forall a. Num a => [[a]] -> Int -> Int -> a
value [[a]]
rows Int
k Int
k
            if a
pivot a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0
                then MathError -> Either MathError [[a]]
forall a b. a -> Either a b
Left (Int -> MathError
SingularMatrix Int
k)
                else do
                    let lower :: [[a]]
lower =
                            [ [ if Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
k Bool -> Bool -> Bool
&& Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
k Bool -> Bool -> Bool
&& Int -> Int -> Bool
active Int
i Int
k then [[a]] -> Int -> Int -> a
forall a. Num a => [[a]] -> Int -> Int -> a
value [[a]]
rows Int
i Int
k a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
pivot else [[a]] -> Int -> Int -> a
forall a. Num a => [[a]] -> Int -> Int -> a
value [[a]]
rows Int
i Int
j | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                              ]
                            | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                            ]
                        next :: [[a]]
next =
                            [ [ if Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
k Bool -> Bool -> Bool
&& Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
k Bool -> Bool -> Bool
&& Int -> Int -> Bool
active Int
i Int
k Bool -> Bool -> Bool
&& Int -> Int -> Bool
active Int
i Int
j
                                    then [[a]] -> Int -> Int -> a
forall a. Num a => [[a]] -> Int -> Int -> a
value [[a]]
lower Int
i Int
j a -> a -> a
forall a. Num a => a -> a -> a
- [[a]] -> Int -> Int -> a
forall a. Num a => [[a]] -> Int -> Int -> a
value [[a]]
lower Int
i Int
k a -> a -> a
forall a. Num a => a -> a -> a
* [[a]] -> Int -> Int -> a
forall a. Num a => [[a]] -> Int -> Int -> a
value [[a]]
lower Int
k Int
j
                                    else [[a]] -> Int -> Int -> a
forall a. Num a => [[a]] -> Int -> Int -> a
value [[a]]
lower Int
i Int
j
                              | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                              ]
                            | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                            ]
                    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues ([[a]] -> [a]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[a]]
next)
                    [[a]] -> Either MathError [[a]]
forall a. a -> Either MathError a
forall (f :: * -> *) a. Applicative f => a -> f a
pure [[a]]
next
    packed <- foldM update original [0 .. n - 1]
    let lower =
            [ [if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then a
1 else if Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
j then [[a]] -> Int -> Int -> a
forall a. Num a => [[a]] -> Int -> Int -> a
value [[a]]
packed Int
i Int
j else a
0 | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
            | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
            ]
        upper = [[if Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
j then [[a]] -> Int -> Int -> a
forall a. Num a => [[a]] -> Int -> Int -> a
value [[a]]
packed Int
i Int
j else a
0 | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]] | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
    pure (triangularPreconditioner lower upper)

{- | Incomplete Cholesky with zero fill for a symmetric positive-definite input.
Dropping can itself cause a nonpositive pivot; that reports NotPositiveDefinite.
-}
ic0Preconditioner :: (RealFloat a) => Matrix a -> Either MathError (Preconditioner a)
ic0Preconditioner :: forall a.
RealFloat a =>
Matrix a -> Either MathError (Preconditioner a)
ic0Preconditioner Matrix a
matrix = do
    n <- Matrix a -> Either MathError Int
forall a. RealFloat a => Matrix a -> Either MathError Int
symmetricFinite Matrix a
matrix
    let a = Matrix a -> [[a]]
forall a. Matrix a -> [[a]]
entries Matrix a
matrix
        initial = Int -> [a] -> [[a]]
forall a. Int -> a -> [a]
replicate Int
n (Int -> a -> [a]
forall a. Int -> a -> [a]
replicate Int
n a
0)
        step [[a]]
rows (Int
i, Int
j)
            | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
j Bool -> Bool -> Bool
&& [[a]] -> Int -> Int -> a
forall a. Num a => [[a]] -> Int -> Int -> a
value [[a]]
a Int
i Int
j a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = [[a]] -> Either MathError [[a]]
forall a b. b -> Either a b
Right [[a]]
rows
            | Bool
otherwise = do
                let partial :: a
partial = [[a]] -> Int -> Int -> a
forall a. Num a => [[a]] -> Int -> Int -> a
value [[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]] -> Int -> Int -> a
forall a. Num a => [[a]] -> Int -> Int -> a
value [[a]]
rows Int
i Int
k a -> a -> a
forall a. Num a => a -> a -> a
* [[a]] -> Int -> Int -> a
forall a. Num a => [[a]] -> Int -> Int -> a
value [[a]]
rows Int
j Int
k | Int
k <- [Int
0 .. Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
                x <-
                    if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j
                        then if a
partial 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 a -> Either MathError a
forall a b. b -> Either a b
Right (a -> a
forall a. Floating a => a -> a
sqrt a
partial)
                        else
                            if [[a]] -> Int -> Int -> a
forall a. Num a => [[a]] -> Int -> Int -> a
value [[a]]
rows Int
j Int
j a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0
                                then MathError -> Either MathError a
forall a b. a -> Either a b
Left (Int -> MathError
SingularMatrix Int
j)
                                else a -> Either MathError a
forall a b. b -> Either a b
Right (a
partial a -> a -> a
forall a. Fractional a => a -> a -> a
/ [[a]] -> Int -> Int -> a
forall a. Num a => [[a]] -> Int -> Int -> a
value [[a]]
rows Int
j Int
j)
                finiteValues [x]
                pure (setValue i j x rows)
    lower <- foldM step initial [(i, j) | i <- [0 .. n - 1], j <- [0 .. i]]
    let upper = [[[[a]] -> Int -> Int -> a
forall a. Num a => [[a]] -> Int -> Int -> a
value [[a]]
lower Int
j Int
i | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]] | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
    pure (triangularPreconditioner lower upper)

triangularPreconditioner :: (RealFloat a) => [[a]] -> [[a]] -> Preconditioner a
triangularPreconditioner :: forall a. RealFloat a => [[a]] -> [[a]] -> Preconditioner a
triangularPreconditioner [[a]]
lower [[a]]
upper =
    ([a] -> Either MathError [a]) -> Preconditioner a
forall a. ([a] -> Either MathError [a]) -> Preconditioner a
Preconditioner
        ( \[a]
b -> do
            Int -> [a] -> Either MathError ()
forall a. RealFloat a => Int -> [a] -> Either MathError ()
validateVector ([[a]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[a]]
lower) [a]
b
            y <- Bool -> [[a]] -> [a] -> Either MathError [a]
forall a.
RealFloat a =>
Bool -> [[a]] -> [a] -> Either MathError [a]
triangularSolve Bool
True [[a]]
lower [a]
b
            triangularSolve False upper y
        )

triangularSolve :: (RealFloat a) => Bool -> [[a]] -> [a] -> Either MathError [a]
triangularSolve :: forall a.
RealFloat a =>
Bool -> [[a]] -> [a] -> Either MathError [a]
triangularSolve Bool
forward [[a]]
a [a]
b = ([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 -> Either MathError [a]
step (Int -> a -> [a]
forall a. Int -> a -> [a]
replicate Int
n a
0) [Int]
order
  where
    n :: Int
n = [[a]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[a]]
a
    order :: [Int]
order = if Bool
forward then [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1] else [Int] -> [Int]
forall a. [a] -> [a]
reverse [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
    step :: [a] -> Int -> Either MathError [a]
step [a]
x Int
i = do
        let pivot :: a
pivot = [[a]] -> Int -> Int -> a
forall a. Num a => [[a]] -> Int -> Int -> a
value [[a]]
a Int
i Int
i
            js :: [Int]
js = if Bool
forward then [Int
0 .. Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1] else [Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
        if a
pivot a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0
            then MathError -> Either MathError [a]
forall a b. a -> Either a b
Left (Int -> MathError
SingularMatrix Int
i)
            else do
                let y :: a
y = ([a] -> Int -> a
forall a. Num a => [a] -> Int -> a
vectorValue [a]
b 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]] -> Int -> Int -> a
forall a. Num a => [[a]] -> Int -> Int -> a
value [[a]]
a Int
i Int
j a -> a -> a
forall a. Num a => a -> a -> a
* [a] -> Int -> a
forall a. Num a => [a] -> Int -> a
vectorValue [a]
x Int
j | Int
j <- [Int]
js]) a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
pivot
                [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues [a
y]
                [a] -> Either MathError [a]
forall a. a -> Either MathError a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Int -> a -> [a] -> [a]
forall a. Int -> a -> [a] -> [a]
replaceAt Int
i a
y [a]
x)

{- | Conjugate gradient for symmetric positive-definite matrices. Symmetry is
checked; nonpositive Krylov curvature rejects indefinite matrices encountered.
-}
conjugateGradient ::
    (RealFloat a) =>
    IterativeOptions a ->
    Matrix a ->
    [a] ->
    [a] ->
    Either MathError (IterativeResult a)
conjugateGradient :: forall a.
RealFloat a =>
IterativeOptions a
-> Matrix a -> [a] -> [a] -> Either MathError (IterativeResult a)
conjugateGradient IterativeOptions a
options Matrix a
matrix [a]
b [a]
x = do
    preconditioner <- Int -> Either MathError (Preconditioner a)
forall a. RealFloat a => Int -> Either MathError (Preconditioner a)
identityPreconditioner ((Int, Int) -> Int
forall a b. (a, b) -> a
fst (Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
matrix))
    preconditionedCG options preconditioner matrix b x

{- | Preconditioned conjugate gradient. The inverse action must be SPD; a
nonpositive r^T M^-1 r or search curvature is reported as a domain error.
-}
preconditionedCG ::
    (RealFloat a) =>
    IterativeOptions a ->
    Preconditioner a ->
    Matrix a ->
    [a] ->
    [a] ->
    Either MathError (IterativeResult a)
preconditionedCG :: forall a.
RealFloat a =>
IterativeOptions a
-> Preconditioner a
-> Matrix a
-> [a]
-> [a]
-> Either MathError (IterativeResult a)
preconditionedCG IterativeOptions a
options Preconditioner a
preconditioner Matrix a
matrix [a]
b [a]
x0 = do
    _ <- Matrix a -> Either MathError Int
forall a. RealFloat a => Matrix a -> Either MathError Int
symmetricFinite Matrix a
matrix
    validateProblem True options matrix b x0
    let r0 = Matrix a -> [a] -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a] -> [a]
residual Matrix a
matrix [a]
b [a]
x0
    if converged options b r0
        then result matrix b 0 x0
        else do
            z0 <- applyPreconditioner preconditioner r0
            let rho0 = [a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
dotV [a]
r0 [a]
z0
            if rho0 <= 0
                then Left (DomainError "PCG preconditioner is not positive definite")
                else go 0 x0 z0 rho0
  where
    go :: Int -> [a] -> [a] -> a -> Either MathError (IterativeResult a)
go Int
k [a]
x [a]
p a
rho
        | Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== IterativeOptions a -> Int
forall a. IterativeOptions a -> Int
maximumIterations IterativeOptions a
options = MathError -> Either MathError (IterativeResult a)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"preconditioned CG" Int
k)
        | Bool
otherwise = do
            let ap :: [a]
ap = Matrix a -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a]
matvec Matrix a
matrix [a]
p; curvature :: a
curvature = [a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
dotV [a]
p [a]
ap
            if a
curvature a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
0
                then MathError -> Either MathError (IterativeResult a)
forall a b. a -> Either a b
Left (Int -> MathError
NotPositiveDefinite Int
k)
                else do
                    let alpha :: a
alpha = a
rho a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
curvature
                        next :: [a]
next = [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
addV [a]
x (a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scaleV a
alpha [a]
p)
                        rnext :: [a]
rnext = Matrix a -> [a] -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a] -> [a]
residual Matrix a
matrix [a]
b [a]
next
                    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues [a]
next
                    if IterativeOptions a -> [a] -> [a] -> Bool
forall a. RealFloat a => IterativeOptions a -> [a] -> [a] -> Bool
converged IterativeOptions a
options [a]
b [a]
rnext
                        then Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
forall a.
RealFloat a =>
Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
result Matrix a
matrix [a]
b (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [a]
next
                        else do
                            z <- Preconditioner a -> [a] -> Either MathError [a]
forall a. Preconditioner a -> [a] -> Either MathError [a]
applyPreconditioner Preconditioner a
preconditioner [a]
rnext
                            let rhonew = [a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
dotV [a]
rnext [a]
z
                            if rhonew <= 0
                                then Left (DomainError "PCG preconditioner is not positive definite")
                                else go (k + 1) next (addV z (scaleV (rhonew / rho) p)) rhonew

{- | Restarted GMRES with modified Gram--Schmidt Arnoldi and a second
orthogonalization pass. Restart length must be positive; total iterations are
bounded by options. Each projected least-squares problem is solved by QR.
-}
gmres ::
    (RealFloat a) =>
    IterativeOptions a ->
    Int ->
    Matrix a ->
    [a] ->
    [a] ->
    Either MathError (IterativeResult a)
gmres :: forall a.
RealFloat a =>
IterativeOptions a
-> Int
-> Matrix a
-> [a]
-> [a]
-> Either MathError (IterativeResult a)
gmres IterativeOptions a
options Int
restart Matrix a
matrix [a]
b [a]
x0
    | Int
restart Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0 = MathError -> Either MathError (IterativeResult a)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"GMRES restart must be positive")
    | Bool
otherwise = do
        Bool
-> IterativeOptions a
-> Matrix a
-> [a]
-> [a]
-> Either MathError ()
forall a.
RealFloat a =>
Bool
-> IterativeOptions a
-> Matrix a
-> [a]
-> [a]
-> Either MathError ()
validateProblem Bool
True IterativeOptions a
options Matrix a
matrix [a]
b [a]
x0
        Int -> [a] -> Either MathError (IterativeResult a)
outer Int
0 [a]
x0
  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
matrix)
    outer :: Int -> [a] -> Either MathError (IterativeResult a)
outer Int
done [a]
x
        | IterativeOptions a -> [a] -> [a] -> Bool
forall a. RealFloat a => IterativeOptions a -> [a] -> [a] -> Bool
converged IterativeOptions a
options [a]
b [a]
r = Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
forall a.
RealFloat a =>
Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
result Matrix a
matrix [a]
b Int
done [a]
x
        | Int
done Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== IterativeOptions a -> Int
forall a. IterativeOptions a -> Int
maximumIterations IterativeOptions a
options = MathError -> Either MathError (IterativeResult a)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"GMRES" Int
done)
        | Bool
otherwise = Int
-> [a]
-> a
-> [[a]]
-> [[a]]
-> Int
-> Either MathError (IterativeResult a)
cycleStep Int
done [a]
x ([a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
r) [[a] -> [a]
forall a. RealFloat a => [a] -> [a]
normalize [a]
r] [] Int
0
      where
        r :: [a]
r = Matrix a -> [a] -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a] -> [a]
residual Matrix a
matrix [a]
b [a]
x
    cycleStep :: Int
-> [a]
-> a
-> [[a]]
-> [[a]]
-> Int
-> Either MathError (IterativeResult a)
cycleStep Int
done [a]
base a
beta [[a]]
basis [[a]]
columns Int
k = do
        let q :: [a]
q = Int -> [[a]] -> [a]
forall a. Int -> [[a]] -> [a]
vectorAt Int
k [[a]]
basis
            ([a]
h, [a]
w) = [[a]] -> [a] -> ([a], [a])
forall a. RealFloat a => [[a]] -> [a] -> ([a], [a])
orthogonalize [[a]]
basis (Matrix a -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a]
matvec Matrix a
matrix [a]
q)
            nextNorm :: a
nextNorm = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
w
            column :: [a]
column = [a]
h [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a
nextNorm]
            cols :: [[a]]
cols = [[a]]
columns [[a]] -> [[a]] -> [[a]]
forall a. [a] -> [a] -> [a]
++ [[a]
column]
            width :: Int
width = Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1
        hmatrix <-
            Int -> Int -> [a] -> Either MathError (Matrix a)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
fromListChecked
                (Int
width Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
                Int
width
                [[a] -> Int -> a
forall a. Num a => [a] -> Int -> a
vectorValue [a]
col Int
i | Int
i <- [Int
0 .. Int
width], [a]
col <- [[a]]
cols]
        y <- leastSquares hmatrix (beta : replicate width 0)
        let x = [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
addV [a]
base ([[a]] -> [a] -> [a]
forall a. Num a => [[a]] -> [a] -> [a]
combine [[a]]
basis [a]
y)
            nextDone = Int
done Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1
        finiteValues x
        if converged options b (residual matrix b x)
            then result matrix b nextDone x
            else
                if nextDone == maximumIterations options
                    then Left (NonConvergence "GMRES" nextDone)
                    else
                        if nextNorm == 0
                            then Left (DomainError "GMRES invariant-subspace breakdown")
                            else
                                if width >= min restart n
                                    then outer nextDone x
                                    else cycleStep nextDone base beta (basis ++ [scaleV (1 / nextNorm) w]) cols width

{- | Minimum residual on the symmetric Lanczos Krylov space, including indefinite
systems. A full tridiagonal projected least-squares problem is solved by QR
each iteration (the mathematical MINRES objective, without short recurrences).
-}
minres ::
    (RealFloat a) =>
    IterativeOptions a ->
    Matrix a ->
    [a] ->
    [a] ->
    Either MathError (IterativeResult a)
minres :: forall a.
RealFloat a =>
IterativeOptions a
-> Matrix a -> [a] -> [a] -> Either MathError (IterativeResult a)
minres IterativeOptions a
options Matrix a
matrix [a]
b [a]
x0 = do
    _ <- Matrix a -> Either MathError Int
forall a. RealFloat a => Matrix a -> Either MathError Int
symmetricFinite Matrix a
matrix
    validateProblem True options matrix b x0
    let r = Matrix a -> [a] -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a] -> [a]
residual Matrix a
matrix [a]
b [a]
x0
    if converged options b r
        then result matrix b 0 x0
        else loop 0 (stableNorm r) [normalize r] [] []
  where
    loop :: Int
-> a -> [[a]] -> [a] -> [a] -> Either MathError (IterativeResult a)
loop Int
k a
beta0 [[a]]
basis [a]
alphas [a]
betas
        | Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== IterativeOptions a -> Int
forall a. IterativeOptions a -> Int
maximumIterations IterativeOptions a
options = MathError -> Either MathError (IterativeResult a)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"MINRES" Int
k)
        | Bool
otherwise = do
            let q :: [a]
q = Int -> [[a]] -> [a]
forall a. Int -> [[a]] -> [a]
vectorAt Int
k [[a]]
basis
                prior :: [a]
prior =
                    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 ([a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
q) a
0
                        else a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scaleV ([a] -> Int -> a
forall a. Num a => [a] -> Int -> a
vectorValue [a]
betas (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)) (Int -> [[a]] -> [a]
forall a. Int -> [[a]] -> [a]
vectorAt (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) [[a]]
basis)
                aq :: [a]
aq = [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
subV (Matrix a -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a]
matvec Matrix a
matrix [a]
q) [a]
prior
                alpha :: a
alpha = [a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
dotV [a]
q [a]
aq
                w :: [a]
w = [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
subV [a]
aq (a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scaleV a
alpha [a]
q)
                beta :: a
beta = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
w
                as :: [a]
as = [a]
alphas [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a
alpha]
                bs :: [a]
bs = [a]
betas [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a
beta]
                width :: Int
width = Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1
                entry :: Int -> Int -> a
entry Int
i Int
j
                    | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j = [a] -> Int -> a
forall a. Num a => [a] -> Int -> a
vectorValue [a]
as Int
j
                    | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 = [a] -> Int -> a
forall a. Num a => [a] -> Int -> a
vectorValue [a]
bs Int
j
                    | Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 = [a] -> Int -> a
forall a. Num a => [a] -> Int -> a
vectorValue [a]
bs Int
i
                    | Bool
otherwise = a
0
            projected <-
                Int -> Int -> [a] -> Either MathError (Matrix a)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
fromListChecked (Int
width Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int
width [Int -> Int -> a
entry Int
i Int
j | Int
i <- [Int
0 .. Int
width], Int
j <- [Int
0 .. Int
width Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
            y <- leastSquares projected (beta0 : replicate width 0)
            let x = [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
addV [a]
x0 ([[a]] -> [a] -> [a]
forall a. Num a => [[a]] -> [a] -> [a]
combine [[a]]
basis [a]
y)
            finiteValues x
            if converged options b (residual matrix b x)
                then result matrix b width x
                else
                    if beta == 0
                        then Left (DomainError "MINRES invariant-subspace breakdown")
                        else loop width beta0 (basis ++ [scaleV (1 / beta) w]) as bs

{- | Bi-conjugate gradient for nonsymmetric square matrices; the initial shadow
residual equals the initial residual. Biorthogonality breakdown fails explicitly.
-}
bicg ::
    (RealFloat a) =>
    IterativeOptions a ->
    Matrix a ->
    [a] ->
    [a] ->
    Either MathError (IterativeResult a)
bicg :: forall a.
RealFloat a =>
IterativeOptions a
-> Matrix a -> [a] -> [a] -> Either MathError (IterativeResult a)
bicg IterativeOptions a
options Matrix a
matrix [a]
b [a]
x0 = do
    Bool
-> IterativeOptions a
-> Matrix a
-> [a]
-> [a]
-> Either MathError ()
forall a.
RealFloat a =>
Bool
-> IterativeOptions a
-> Matrix a
-> [a]
-> [a]
-> Either MathError ()
validateProblem Bool
True IterativeOptions a
options Matrix a
matrix [a]
b [a]
x0
    let r :: [a]
r = Matrix a -> [a] -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a] -> [a]
residual Matrix a
matrix [a]
b [a]
x0
    if IterativeOptions a -> [a] -> [a] -> Bool
forall a. RealFloat a => IterativeOptions a -> [a] -> [a] -> Bool
converged IterativeOptions a
options [a]
b [a]
r
        then Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
forall a.
RealFloat a =>
Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
result Matrix a
matrix [a]
b Int
0 [a]
x0
        else Int
-> [a]
-> [a]
-> [a]
-> [a]
-> [a]
-> a
-> Either MathError (IterativeResult a)
go Int
0 [a]
x0 [a]
r [a]
r [a]
r [a]
r ([a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
dotV [a]
r [a]
r)
  where
    transposed :: Matrix a
transposed = Matrix a -> Matrix a
forall a. Matrix a -> Matrix a
transposeMatrix Matrix a
matrix
    go :: Int
-> [a]
-> [a]
-> [a]
-> [a]
-> [a]
-> a
-> Either MathError (IterativeResult a)
go Int
k [a]
x [a]
r [a]
shadow [a]
p [a]
phat a
rho
        | Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== IterativeOptions a -> Int
forall a. IterativeOptions a -> Int
maximumIterations IterativeOptions a
options = MathError -> Either MathError (IterativeResult a)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"BiCG" Int
k)
        | a
rho a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = MathError -> Either MathError (IterativeResult a)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"BiCG shadow breakdown")
        | Bool
otherwise = do
            let ap :: [a]
ap = Matrix a -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a]
matvec Matrix a
matrix [a]
p; atp :: [a]
atp = Matrix a -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a]
matvec Matrix a
transposed [a]
phat; denom :: a
denom = [a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
dotV [a]
phat [a]
ap
            if a
denom a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0
                then MathError -> Either MathError (IterativeResult a)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"BiCG projection breakdown")
                else do
                    let alpha :: a
alpha = a
rho a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
denom
                        xn :: [a]
xn = [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
addV [a]
x (a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scaleV a
alpha [a]
p)
                        rn :: [a]
rn = [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
subV [a]
r (a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scaleV a
alpha [a]
ap)
                        sn :: [a]
sn = [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
subV [a]
shadow (a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scaleV a
alpha [a]
atp)
                        rhon :: a
rhon = [a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
dotV [a]
sn [a]
rn
                    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues ([a]
xn [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a]
rn [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a]
sn)
                    if IterativeOptions a -> [a] -> [a] -> Bool
forall a. RealFloat a => IterativeOptions a -> [a] -> [a] -> Bool
converged IterativeOptions a
options [a]
b (Matrix a -> [a] -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a] -> [a]
residual Matrix a
matrix [a]
b [a]
xn)
                        then Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
forall a.
RealFloat a =>
Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
result Matrix a
matrix [a]
b (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [a]
xn
                        else
                            Int
-> [a]
-> [a]
-> [a]
-> [a]
-> [a]
-> a
-> Either MathError (IterativeResult a)
go
                                (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
                                [a]
xn
                                [a]
rn
                                [a]
sn
                                ([a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
addV [a]
rn (a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scaleV (a
rhon a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
rho) [a]
p))
                                ([a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
addV [a]
sn (a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scaleV (a
rhon a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
rho) [a]
phat))
                                a
rhon

{- | Stabilized BiCG with a fixed initial shadow residual and a minimal-residual
omega step. Zero projection/omega denominators return explicit breakdown errors.
-}
bicgstab ::
    (RealFloat a) =>
    IterativeOptions a ->
    Matrix a ->
    [a] ->
    [a] ->
    Either MathError (IterativeResult a)
bicgstab :: forall a.
RealFloat a =>
IterativeOptions a
-> Matrix a -> [a] -> [a] -> Either MathError (IterativeResult a)
bicgstab IterativeOptions a
options Matrix a
matrix [a]
b [a]
x0 = do
    Bool
-> IterativeOptions a
-> Matrix a
-> [a]
-> [a]
-> Either MathError ()
forall a.
RealFloat a =>
Bool
-> IterativeOptions a
-> Matrix a
-> [a]
-> [a]
-> Either MathError ()
validateProblem Bool
True IterativeOptions a
options Matrix a
matrix [a]
b [a]
x0
    let r :: [a]
r = Matrix a -> [a] -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a] -> [a]
residual Matrix a
matrix [a]
b [a]
x0; zeroV :: [a]
zeroV = 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]
x0) a
0
    if IterativeOptions a -> [a] -> [a] -> Bool
forall a. RealFloat a => IterativeOptions a -> [a] -> [a] -> Bool
converged IterativeOptions a
options [a]
b [a]
r then Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
forall a.
RealFloat a =>
Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
result Matrix a
matrix [a]
b Int
0 [a]
x0 else [a]
-> Int
-> [a]
-> [a]
-> [a]
-> [a]
-> a
-> a
-> a
-> Either MathError (IterativeResult a)
go [a]
r Int
0 [a]
x0 [a]
r [a]
zeroV [a]
zeroV a
1 a
1 a
1
  where
    go :: [a]
-> Int
-> [a]
-> [a]
-> [a]
-> [a]
-> a
-> a
-> a
-> Either MathError (IterativeResult a)
go [a]
shadow Int
k [a]
x [a]
r [a]
p [a]
v a
rhoOld a
alphaOld a
omegaOld
        | Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== IterativeOptions a -> Int
forall a. IterativeOptions a -> Int
maximumIterations IterativeOptions a
options = MathError -> Either MathError (IterativeResult a)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"BiCGSTAB" Int
k)
        | Bool
otherwise = do
            let rho :: a
rho = [a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
dotV [a]
shadow [a]
r
            if a
rho a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 Bool -> Bool -> Bool
|| a
omegaOld a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0
                then MathError -> Either MathError (IterativeResult a)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"BiCGSTAB shadow/omega breakdown")
                else do
                    let beta :: a
beta = (a
rho a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
rhoOld) a -> a -> a
forall a. Num a => a -> a -> a
* (a
alphaOld a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
omegaOld)
                        pn :: [a]
pn = [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
addV [a]
r (a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scaleV a
beta ([a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
subV [a]
p (a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scaleV a
omegaOld [a]
v)))
                        vn :: [a]
vn = Matrix a -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a]
matvec Matrix a
matrix [a]
pn
                        denom :: a
denom = [a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
dotV [a]
shadow [a]
vn
                    if a
denom a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0
                        then MathError -> Either MathError (IterativeResult a)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"BiCGSTAB projection breakdown")
                        else do
                            let alpha :: a
alpha = a
rho a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
denom
                                s :: [a]
s = [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
subV [a]
r (a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scaleV a
alpha [a]
vn)
                                xa :: [a]
xa = [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
addV [a]
x (a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scaleV a
alpha [a]
pn)
                            if IterativeOptions a -> [a] -> [a] -> Bool
forall a. RealFloat a => IterativeOptions a -> [a] -> [a] -> Bool
converged IterativeOptions a
options [a]
b (Matrix a -> [a] -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a] -> [a]
residual Matrix a
matrix [a]
b [a]
xa)
                                then Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
forall a.
RealFloat a =>
Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
result Matrix a
matrix [a]
b (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [a]
xa
                                else do
                                    let t :: [a]
t = Matrix a -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a]
matvec Matrix a
matrix [a]
s; tt :: a
tt = [a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
dotV [a]
t [a]
t
                                    if a
tt a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0
                                        then MathError -> Either MathError (IterativeResult a)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"BiCGSTAB residual breakdown")
                                        else do
                                            let omega :: a
omega = [a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
dotV [a]
t [a]
s a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
tt
                                                xn :: [a]
xn = [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
addV [a]
xa (a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scaleV a
omega [a]
s)
                                                rn :: [a]
rn = [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
subV [a]
s (a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scaleV a
omega [a]
t)
                                            [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues ([a]
xn [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a]
rn)
                                            if IterativeOptions a -> [a] -> [a] -> Bool
forall a. RealFloat a => IterativeOptions a -> [a] -> [a] -> Bool
converged IterativeOptions a
options [a]
b (Matrix a -> [a] -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a] -> [a]
residual Matrix a
matrix [a]
b [a]
xn)
                                                then Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
forall a.
RealFloat a =>
Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
result Matrix a
matrix [a]
b (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [a]
xn
                                                else [a]
-> Int
-> [a]
-> [a]
-> [a]
-> [a]
-> a
-> a
-> a
-> Either MathError (IterativeResult a)
go [a]
shadow (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [a]
xn [a]
rn [a]
pn [a]
vn a
rho a
alpha a
omega

{- | LSQR by Golub--Kahan bidiagonalization and plane-rotation recurrences.
Rectangular and rank-deficient systems are supported. The least-squares
stopping test uses the original normal residual as well as the true residual.
-}
lsqr ::
    (RealFloat a) =>
    IterativeOptions a ->
    Matrix a ->
    [a] ->
    [a] ->
    Either MathError (IterativeResult a)
lsqr :: forall a.
RealFloat a =>
IterativeOptions a
-> Matrix a -> [a] -> [a] -> Either MathError (IterativeResult a)
lsqr IterativeOptions a
options Matrix a
matrix [a]
b [a]
x0 = do
    Bool
-> IterativeOptions a
-> Matrix a
-> [a]
-> [a]
-> Either MathError ()
forall a.
RealFloat a =>
Bool
-> IterativeOptions a
-> Matrix a
-> [a]
-> [a]
-> Either MathError ()
validateProblem Bool
False IterativeOptions a
options Matrix a
matrix [a]
b [a]
x0
    let r :: [a]
r = Matrix a -> [a] -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a] -> [a]
residual Matrix a
matrix [a]
b [a]
x0
        beta :: a
beta = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
r
        u :: [a]
u = [a] -> [a]
forall a. RealFloat a => [a] -> [a]
normalize [a]
r
        atU :: [a]
atU = Matrix a -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a]
matvec Matrix a
transposed [a]
u
        alpha :: a
alpha = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
atU
        v :: [a]
v = [a] -> [a]
forall a. RealFloat a => [a] -> [a]
normalize [a]
atU
    if IterativeOptions a -> Matrix a -> [a] -> [a] -> Bool
forall a.
RealFloat a =>
IterativeOptions a -> Matrix a -> [a] -> [a] -> Bool
leastConverged IterativeOptions a
options Matrix a
matrix [a]
b [a]
r
        then Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
forall a.
RealFloat a =>
Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
result Matrix a
matrix [a]
b Int
0 [a]
x0
        else Int
-> [a]
-> [a]
-> [a]
-> [a]
-> a
-> a
-> a
-> Either MathError (IterativeResult a)
go Int
0 [a]
x0 [a]
u [a]
v [a]
v a
alpha a
beta a
alpha
  where
    transposed :: Matrix a
transposed = Matrix a -> Matrix a
forall a. Matrix a -> Matrix a
transposeMatrix Matrix a
matrix
    go :: Int
-> [a]
-> [a]
-> [a]
-> [a]
-> a
-> a
-> a
-> Either MathError (IterativeResult a)
go Int
k [a]
x [a]
u [a]
v [a]
w a
alpha a
phibar a
rhobar
        | Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== IterativeOptions a -> Int
forall a. IterativeOptions a -> Int
maximumIterations IterativeOptions a
options = MathError -> Either MathError (IterativeResult a)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"LSQR" Int
k)
        | Bool
otherwise = do
            let unscaled :: [a]
unscaled = [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
subV (Matrix a -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a]
matvec Matrix a
matrix [a]
v) (a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scaleV a
alpha [a]
u)
                beta :: a
beta = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
unscaled
                un :: [a]
un = [a] -> [a]
forall a. RealFloat a => [a] -> [a]
normalize [a]
unscaled
                vscaled :: [a]
vscaled = [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
subV (Matrix a -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a]
matvec Matrix a
transposed [a]
un) (a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scaleV a
beta [a]
v)
                alphan :: a
alphan = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
vscaled
                vn :: [a]
vn = [a] -> [a]
forall a. RealFloat a => [a] -> [a]
normalize [a]
vscaled
                rho :: a
rho = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a
rhobar, a
beta]
            if a
rho a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0
                then MathError -> Either MathError (IterativeResult a)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"LSQR bidiagonal breakdown")
                else do
                    let c :: a
c = a
rhobar a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
rho
                        s :: a
s = a
beta a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
rho
                        theta :: a
theta = a
s a -> a -> a
forall a. Num a => a -> a -> a
* a
alphan
                        phi :: a
phi = a
c a -> a -> a
forall a. Num a => a -> a -> a
* a
phibar
                        xn :: [a]
xn = [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
addV [a]
x (a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scaleV (a
phi a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
rho) [a]
w)
                        wn :: [a]
wn = [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
subV [a]
vn (a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scaleV (a
theta a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
rho) [a]
w)
                    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues ([a]
xn [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a]
wn)
                    if IterativeOptions a -> Matrix a -> [a] -> [a] -> Bool
forall a.
RealFloat a =>
IterativeOptions a -> Matrix a -> [a] -> [a] -> Bool
leastConverged IterativeOptions a
options Matrix a
matrix [a]
b (Matrix a -> [a] -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a] -> [a]
residual Matrix a
matrix [a]
b [a]
xn)
                        then Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
forall a.
RealFloat a =>
Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
result Matrix a
matrix [a]
b (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [a]
xn
                        else Int
-> [a]
-> [a]
-> [a]
-> [a]
-> a
-> a
-> a
-> Either MathError (IterativeResult a)
go (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [a]
xn [a]
un [a]
vn [a]
wn a
alphan (a
s a -> a -> a
forall a. Num a => a -> a -> a
* a
phibar) (-a
c a -> a -> a
forall a. Num a => a -> a -> a
* a
alphan)

{- | LSMR's minimum-normal-residual Krylov objective using explicit Golub--Kahan
bases and QR on the projected tridiagonal normal operator. A^T A is never
formed. This full-basis reference version uses more storage than short-recurrence
implementations; it minimizes ||A^T(b-Ax)|| over the generated Krylov space.
-}
lsmr ::
    (RealFloat a) =>
    IterativeOptions a ->
    Matrix a ->
    [a] ->
    [a] ->
    Either MathError (IterativeResult a)
lsmr :: forall a.
RealFloat a =>
IterativeOptions a
-> Matrix a -> [a] -> [a] -> Either MathError (IterativeResult a)
lsmr IterativeOptions a
options Matrix a
matrix [a]
b [a]
x0 = do
    Bool
-> IterativeOptions a
-> Matrix a
-> [a]
-> [a]
-> Either MathError ()
forall a.
RealFloat a =>
Bool
-> IterativeOptions a
-> Matrix a
-> [a]
-> [a]
-> Either MathError ()
validateProblem Bool
False IterativeOptions a
options Matrix a
matrix [a]
b [a]
x0
    let r :: [a]
r = Matrix a -> [a] -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a] -> [a]
residual Matrix a
matrix [a]
b [a]
x0
        beta0 :: a
beta0 = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
r
        u0 :: [a]
u0 = [a] -> [a]
forall a. RealFloat a => [a] -> [a]
normalize [a]
r
        vraw :: [a]
vraw = Matrix a -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a]
matvec Matrix a
transposed [a]
u0
        alpha0 :: a
alpha0 = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
vraw
        v0 :: [a]
v0 = [a] -> [a]
forall a. RealFloat a => [a] -> [a]
normalize [a]
vraw
    if IterativeOptions a -> Matrix a -> [a] -> [a] -> Bool
forall a.
RealFloat a =>
IterativeOptions a -> Matrix a -> [a] -> [a] -> Bool
leastConverged IterativeOptions a
options Matrix a
matrix [a]
b [a]
r
        then Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
forall a.
RealFloat a =>
Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
result Matrix a
matrix [a]
b Int
0 [a]
x0
        else Int
-> a
-> [a]
-> [[a]]
-> [a]
-> [a]
-> Either MathError (IterativeResult a)
loop Int
0 (a
alpha0 a -> a -> a
forall a. Num a => a -> a -> a
* a
beta0) [a]
u0 [[a]
v0] [a
alpha0] [a
beta0]
  where
    transposed :: Matrix a
transposed = Matrix a -> Matrix a
forall a. Matrix a -> Matrix a
transposeMatrix Matrix a
matrix
    loop :: Int
-> a
-> [a]
-> [[a]]
-> [a]
-> [a]
-> Either MathError (IterativeResult a)
loop Int
k a
gamma [a]
u [[a]]
basis [a]
alphas [a]
betas
        | Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== IterativeOptions a -> Int
forall a. IterativeOptions a -> Int
maximumIterations IterativeOptions a
options = MathError -> Either MathError (IterativeResult a)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"LSMR" Int
k)
        | Bool
otherwise = do
            let v :: [a]
v = Int -> [[a]] -> [a]
forall a. Int -> [[a]] -> [a]
vectorAt Int
k [[a]]
basis
                alpha :: a
alpha = [a] -> Int -> a
forall a. Num a => [a] -> Int -> a
vectorValue [a]
alphas Int
k
                rawU :: [a]
rawU = [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
subV (Matrix a -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a]
matvec Matrix a
matrix [a]
v) (a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scaleV a
alpha [a]
u)
                beta :: a
beta = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
rawU
                un :: [a]
un = [a] -> [a]
forall a. RealFloat a => [a] -> [a]
normalize [a]
rawU
                rawV :: [a]
rawV = [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
subV (Matrix a -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a]
matvec Matrix a
transposed [a]
un) (a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scaleV a
beta [a]
v)
                alphan :: a
alphan = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
rawV
                vn :: [a]
vn = [a] -> [a]
forall a. RealFloat a => [a] -> [a]
normalize [a]
rawV
                as :: [a]
as = [a]
alphas [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a
alphan]
                bs :: [a]
bs = [a]
betas [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a
beta]
                width :: Int
width = Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1
                entry :: Int -> Int -> a
entry Int
i Int
j
                    | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j = [a] -> Int -> a
forall a. Num a => [a] -> Int -> a
vectorValue [a]
as Int
j a -> a -> a
forall a. Num a => a -> a -> a
* [a] -> Int -> a
forall a. Num a => [a] -> Int -> a
vectorValue [a]
as Int
j a -> a -> a
forall a. Num a => a -> a -> a
+ [a] -> Int -> a
forall a. Num a => [a] -> Int -> a
vectorValue [a]
bs (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) a -> a -> a
forall a. Num a => a -> a -> a
* [a] -> Int -> a
forall a. Num a => [a] -> Int -> a
vectorValue [a]
bs (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
                    | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 = [a] -> Int -> a
forall a. Num a => [a] -> Int -> a
vectorValue [a]
bs (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) a -> a -> a
forall a. Num a => a -> a -> a
* [a] -> Int -> a
forall a. Num a => [a] -> Int -> a
vectorValue [a]
as (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
                    | Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 = [a] -> Int -> a
forall a. Num a => [a] -> Int -> a
vectorValue [a]
as Int
j a -> a -> a
forall a. Num a => a -> a -> a
* [a] -> Int -> a
forall a. Num a => [a] -> Int -> a
vectorValue [a]
bs Int
j
                    | Bool
otherwise = a
0
            projected <-
                Int -> Int -> [a] -> Either MathError (Matrix a)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
fromListChecked (Int
width Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int
width [Int -> Int -> a
entry Int
i Int
j | Int
i <- [Int
0 .. Int
width], Int
j <- [Int
0 .. Int
width Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
            y <- leastSquares projected (gamma : replicate width 0)
            let x = [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
addV [a]
x0 ([[a]] -> [a] -> [a]
forall a. Num a => [[a]] -> [a] -> [a]
combine [[a]]
basis [a]
y)
            finiteValues x
            if leastConverged options matrix b (residual matrix b x)
                then result matrix b width x
                else
                    if alphan == 0 || beta == 0
                        then Left (DomainError "LSMR invariant-subspace breakdown")
                        else loop width gamma un (basis ++ [vn]) as bs

{- | Simultaneous Jacobi stationary iteration. Convergence requires an appropriate
splitting (for example strict diagonal dominance); it is not assumed.
-}
jacobiSolve ::
    (RealFloat a) =>
    IterativeOptions a ->
    Matrix a ->
    [a] ->
    [a] ->
    Either MathError (IterativeResult a)
jacobiSolve :: forall a.
RealFloat a =>
IterativeOptions a
-> Matrix a -> [a] -> [a] -> Either MathError (IterativeResult a)
jacobiSolve = Bool
-> a
-> IterativeOptions a
-> Matrix a
-> [a]
-> [a]
-> Either MathError (IterativeResult a)
forall a.
RealFloat a =>
Bool
-> a
-> IterativeOptions a
-> Matrix a
-> [a]
-> [a]
-> Either MathError (IterativeResult a)
stationary Bool
False a
1

-- | Forward Gauss--Seidel stationary iteration.
gaussSeidelSolve ::
    (RealFloat a) =>
    IterativeOptions a ->
    Matrix a ->
    [a] ->
    [a] ->
    Either MathError (IterativeResult a)
gaussSeidelSolve :: forall a.
RealFloat a =>
IterativeOptions a
-> Matrix a -> [a] -> [a] -> Either MathError (IterativeResult a)
gaussSeidelSolve = Bool
-> a
-> IterativeOptions a
-> Matrix a
-> [a]
-> [a]
-> Either MathError (IterativeResult a)
forall a.
RealFloat a =>
Bool
-> a
-> IterativeOptions a
-> Matrix a
-> [a]
-> [a]
-> Either MathError (IterativeResult a)
stationary Bool
True a
1

-- | Successive over-relaxation, requiring a finite relaxation 0<omega<2.
sorSolve ::
    (RealFloat a) =>
    IterativeOptions a ->
    a ->
    Matrix a ->
    [a] ->
    [a] ->
    Either MathError (IterativeResult a)
sorSolve :: forall a.
RealFloat a =>
IterativeOptions a
-> a
-> Matrix a
-> [a]
-> [a]
-> Either MathError (IterativeResult a)
sorSolve IterativeOptions a
options a
omega = Bool
-> a
-> IterativeOptions a
-> Matrix a
-> [a]
-> [a]
-> Either MathError (IterativeResult a)
forall a.
RealFloat a =>
Bool
-> a
-> IterativeOptions a
-> Matrix a
-> [a]
-> [a]
-> Either MathError (IterativeResult a)
stationary Bool
True a
omega IterativeOptions a
options

stationary ::
    (RealFloat a) =>
    Bool ->
    a ->
    IterativeOptions a ->
    Matrix a ->
    [a] ->
    [a] ->
    Either MathError (IterativeResult a)
stationary :: forall a.
RealFloat a =>
Bool
-> a
-> IterativeOptions a
-> Matrix a
-> [a]
-> [a]
-> Either MathError (IterativeResult a)
stationary Bool
sequential a
omega IterativeOptions a
options Matrix a
matrix [a]
b [a]
x0
    | Bool -> Bool
not (a -> Bool
forall a. RealFloat a => a -> Bool
finite a
omega) Bool -> Bool -> Bool
|| a
omega a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
0 Bool -> Bool -> Bool
|| a
omega a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
2 =
        MathError -> Either MathError (IterativeResult a)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"SOR relaxation must lie in (0,2)")
    | Bool
otherwise = do
        Bool
-> IterativeOptions a
-> Matrix a
-> [a]
-> [a]
-> Either MathError ()
forall a.
RealFloat a =>
Bool
-> IterativeOptions a
-> Matrix a
-> [a]
-> [a]
-> Either MathError ()
validateProblem Bool
True IterativeOptions a
options Matrix a
matrix [a]
b [a]
x0
        [a] -> Either MathError ()
forall a. (Eq a, Num a) => [a] -> Either MathError ()
checkPivots [[[a]] -> Int -> Int -> a
forall a. Num a => [[a]] -> Int -> Int -> a
value [[a]]
a Int
i Int
i | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
        Int -> [a] -> Either MathError (IterativeResult a)
go Int
0 [a]
x0
  where
    a :: [[a]]
a = Matrix a -> [[a]]
forall a. Matrix a -> [[a]]
entries Matrix a
matrix
    n :: Int
n = [[a]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[a]]
a
    go :: Int -> [a] -> Either MathError (IterativeResult a)
go Int
k [a]
x
        | IterativeOptions a -> [a] -> [a] -> Bool
forall a. RealFloat a => IterativeOptions a -> [a] -> [a] -> Bool
converged IterativeOptions a
options [a]
b (Matrix a -> [a] -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a] -> [a]
residual Matrix a
matrix [a]
b [a]
x) = Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
forall a.
RealFloat a =>
Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
result Matrix a
matrix [a]
b Int
k [a]
x
        | Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== IterativeOptions a -> Int
forall a. IterativeOptions a -> Int
maximumIterations IterativeOptions a
options = MathError -> Either MathError (IterativeResult a)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"stationary iteration" Int
k)
        | Bool
otherwise = do
            let update :: [a] -> Int -> [a]
update [a]
acc Int
i =
                    let use :: [a]
use = if Bool
sequential then [a]
acc else [a]
x
                        rhs :: a
rhs = [a] -> Int -> a
forall a. Num a => [a] -> Int -> a
vectorValue [a]
b 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]] -> Int -> Int -> a
forall a. Num a => [[a]] -> Int -> Int -> a
value [[a]]
a Int
i Int
j a -> a -> a
forall a. Num a => a -> a -> a
* [a] -> Int -> a
forall a. Num a => [a] -> Int -> a
vectorValue [a]
use Int
j | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
i]
                        new :: a
new = (a
1 a -> a -> a
forall a. Num a => a -> a -> a
- a
omega) a -> a -> a
forall a. Num a => a -> a -> a
* [a] -> Int -> a
forall a. Num a => [a] -> Int -> a
vectorValue [a]
x Int
i a -> a -> a
forall a. Num a => a -> a -> a
+ a
omega a -> a -> a
forall a. Num a => a -> a -> a
* a
rhs a -> a -> a
forall a. Fractional a => a -> a -> a
/ [[a]] -> Int -> Int -> a
forall a. Num a => [[a]] -> Int -> Int -> a
value [[a]]
a Int
i Int
i
                     in Int -> a -> [a] -> [a]
forall a. Int -> a -> [a] -> [a]
replaceAt Int
i a
new [a]
acc
                xn :: [a]
xn = ([a] -> Int -> [a]) -> [a] -> [Int] -> [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] -> Int -> [a]
update [a]
x [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
            [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues [a]
xn
            Int -> [a] -> Either MathError (IterativeResult a)
go (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [a]
xn

{- | Residual-correction refinement using a reusable partial-pivoted LU.
Residuals and corrections use the scalar type's precision. This does not claim
mixed-precision accuracy; callers may provide an initial guess from any solver.
-}
iterativeRefinement ::
    (RealFloat a) =>
    IterativeOptions a ->
    Matrix a ->
    [a] ->
    [a] ->
    Either MathError (IterativeResult a)
iterativeRefinement :: forall a.
RealFloat a =>
IterativeOptions a
-> Matrix a -> [a] -> [a] -> Either MathError (IterativeResult a)
iterativeRefinement IterativeOptions a
options Matrix a
matrix [a]
b [a]
x0 = do
    Bool
-> IterativeOptions a
-> Matrix a
-> [a]
-> [a]
-> Either MathError ()
forall a.
RealFloat a =>
Bool
-> IterativeOptions a
-> Matrix a
-> [a]
-> [a]
-> Either MathError ()
validateProblem Bool
True IterativeOptions a
options Matrix a
matrix [a]
b [a]
x0
    factors <- a -> Matrix a -> Either MathError (LU a)
forall a. RealFloat a => a -> Matrix a -> Either MathError (LU a)
luDecomposeWithTolerance a
0 Matrix a
matrix
    let go Int
k [a]
x
            | IterativeOptions a -> [a] -> [a] -> Bool
forall a. RealFloat a => IterativeOptions a -> [a] -> [a] -> Bool
converged IterativeOptions a
options [a]
b (Matrix a -> [a] -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a] -> [a]
residual Matrix a
matrix [a]
b [a]
x) = Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
forall a.
RealFloat a =>
Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
result Matrix a
matrix [a]
b Int
k [a]
x
            | Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== IterativeOptions a -> Int
forall a. IterativeOptions a -> Int
maximumIterations IterativeOptions a
options = MathError -> Either MathError (IterativeResult a)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"iterative refinement" Int
k)
            | Bool
otherwise = do
                correction <- LU a -> [a] -> Either MathError [a]
forall a. RealFloat a => LU a -> [a] -> Either MathError [a]
solveLU LU a
factors (Matrix a -> [a] -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a] -> [a]
residual Matrix a
matrix [a]
b [a]
x)
                let xn = [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
addV [a]
x [a]
correction
                finiteValues xn
                go (k + 1) xn
    go 0 x0

validateProblem ::
    (RealFloat a) =>
    Bool ->
    IterativeOptions a ->
    Matrix a ->
    [a] ->
    [a] ->
    Either MathError ()
validateProblem :: forall a.
RealFloat a =>
Bool
-> IterativeOptions a
-> Matrix a
-> [a]
-> [a]
-> Either MathError ()
validateProblem Bool
square IterativeOptions a
options Matrix a
matrix [a]
b [a]
x = do
    if (a -> Bool) -> [a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (\a
v -> Bool -> Bool
not (a -> Bool
forall a. RealFloat a => a -> Bool
finite a
v) Bool -> Bool -> Bool
|| a
v a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0) [IterativeOptions a -> a
forall a. IterativeOptions a -> a
absoluteTolerance IterativeOptions a
options, IterativeOptions a -> a
forall a. IterativeOptions a -> a
relativeTolerance IterativeOptions a
options]
        then MathError -> Either MathError ()
forall a b. a -> Either a b
Left MathError
InvalidTolerance
        else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    if IterativeOptions a -> Int
forall a. IterativeOptions a -> Int
maximumIterations IterativeOptions a
options Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0
        then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"negative iteration budget")
        else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    let (Int
m, Int
n) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
matrix
    if Bool
square Bool -> Bool -> Bool
&& Int
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n then MathError -> Either MathError ()
forall a b. a -> Either a b
Left ((Int, Int) -> MathError
NonSquareMatrix (Int
m, Int
n)) else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    Matrix a -> Either MathError ()
forall a. RealFloat a => Matrix a -> Either MathError ()
finiteMatrix Matrix a
matrix
    Int -> [a] -> Either MathError ()
forall a. RealFloat a => Int -> [a] -> Either MathError ()
validateVector Int
m [a]
b
    Int -> [a] -> Either MathError ()
forall a. RealFloat a => Int -> [a] -> Either MathError ()
validateVector Int
n [a]
x

squareFinite :: (RealFloat a) => Matrix a -> Either MathError Int
squareFinite :: forall a. RealFloat a => Matrix a -> Either MathError Int
squareFinite Matrix a
matrix = do
    let (Int
m, Int
n) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
matrix
    if Int
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n then MathError -> Either MathError Int
forall a b. a -> Either a b
Left ((Int, Int) -> MathError
NonSquareMatrix (Int
m, Int
n)) else Matrix a -> Either MathError ()
forall a. RealFloat a => Matrix a -> Either MathError ()
finiteMatrix Matrix a
matrix 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

symmetricFinite :: (RealFloat a) => Matrix a -> Either MathError Int
symmetricFinite :: forall a. RealFloat a => Matrix a -> Either MathError Int
symmetricFinite Matrix a
matrix = do
    n <- Matrix a -> Either MathError Int
forall a. RealFloat a => Matrix a -> Either MathError Int
squareFinite Matrix a
matrix
    let a = Matrix a -> [[a]]
forall a. Matrix a -> [[a]]
entries Matrix a
matrix
    if and [value a i j == value a j i | i <- [0 .. n - 1], j <- [0 .. i - 1]]
        then Right n
        else Left (DomainError "symmetric solver requires an exactly symmetric matrix")

validateVector :: (RealFloat a) => Int -> [a] -> Either MathError ()
validateVector :: forall a. RealFloat a => Int -> [a] -> Either MathError ()
validateVector Int
n [a]
xs
    | [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
xs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n = MathError -> Either MathError ()
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]
xs))
    | Bool
otherwise = [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues [a]
xs

checkPivots :: (Eq a, Num a) => [a] -> Either MathError ()
checkPivots :: forall a. (Eq a, Num a) => [a] -> Either MathError ()
checkPivots [a]
xs = case [Int
i | (Int
i, a
x) <- [Int] -> [a] -> [(Int, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [a]
xs, a
x a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0] of
    Int
i : [Int]
_ -> MathError -> Either MathError ()
forall a b. a -> Either a b
Left (Int -> MathError
SingularMatrix Int
i)
    [] -> () -> Either MathError ()
forall a b. b -> Either a b
Right ()

result :: (RealFloat a) => Matrix a -> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
result :: forall a.
RealFloat a =>
Matrix a
-> [a] -> Int -> [a] -> Either MathError (IterativeResult a)
result Matrix a
matrix [a]
b Int
k [a]
x = do
    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues [a]
x
    let r :: [a]
r = Matrix a -> [a] -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a] -> [a]
residual Matrix a
matrix [a]
b [a]
x
        normal :: [a]
normal = Matrix a -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a]
matvec (Matrix a -> Matrix a
forall a. Matrix a -> Matrix a
transposeMatrix Matrix a
matrix) [a]
r
        rn :: a
rn = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
r
        an :: a
an = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
normal
    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues [a]
r
    [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues [a]
normal
    if (a -> Bool) -> [a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all a -> Bool
forall a. RealFloat a => a -> Bool
finite [a
rn, a
an] then IterativeResult a -> Either MathError (IterativeResult a)
forall a b. b -> Either a b
Right ([a] -> Int -> a -> a -> IterativeResult a
forall a. [a] -> Int -> a -> a -> IterativeResult a
IterativeResult [a]
x Int
k a
rn a
an) else MathError -> Either MathError (IterativeResult a)
forall a b. a -> Either a b
Left MathError
NonFiniteResult

converged :: (RealFloat a) => IterativeOptions a -> [a] -> [a] -> Bool
converged :: forall a. RealFloat a => IterativeOptions a -> [a] -> [a] -> Bool
converged IterativeOptions a
options [a]
b [a]
r =
    let norm :: a
norm = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
r
        threshold :: a
threshold = IterativeOptions a -> a
forall a. IterativeOptions a -> a
absoluteTolerance IterativeOptions a
options a -> a -> a
forall a. Num a => a -> a -> a
+ IterativeOptions a -> a
forall a. IterativeOptions a -> a
relativeTolerance IterativeOptions a
options a -> a -> a
forall a. Num a => a -> a -> a
* [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
b
     in (a -> Bool) -> [a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all a -> Bool
forall a. RealFloat a => a -> Bool
finite [a]
r Bool -> Bool -> Bool
&& (a -> Bool) -> [a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all a -> Bool
forall a. RealFloat a => a -> Bool
finite [a]
b Bool -> Bool -> Bool
&& a -> Bool
forall a. RealFloat a => a -> Bool
finite a
norm Bool -> Bool -> Bool
&& a -> Bool
forall a. RealFloat a => a -> Bool
finite a
threshold Bool -> Bool -> Bool
&& a
norm a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
threshold

leastConverged :: (RealFloat a) => IterativeOptions a -> Matrix a -> [a] -> [a] -> Bool
leastConverged :: forall a.
RealFloat a =>
IterativeOptions a -> Matrix a -> [a] -> [a] -> Bool
leastConverged IterativeOptions a
options Matrix a
matrix [a]
b [a]
r =
    (a -> Bool) -> [a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all a -> Bool
forall a. RealFloat a => a -> Bool
finite [a]
r
        Bool -> Bool -> Bool
&& ( IterativeOptions a -> [a] -> [a] -> Bool
forall a. RealFloat a => IterativeOptions a -> [a] -> [a] -> Bool
converged IterativeOptions a
options [a]
b [a]
r
                Bool -> Bool -> Bool
|| IterativeOptions a -> [a] -> [a] -> Bool
forall a. RealFloat a => IterativeOptions a -> [a] -> [a] -> Bool
converged IterativeOptions a
options (Matrix a -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a]
matvec (Matrix a -> Matrix a
forall a. Matrix a -> Matrix a
transposeMatrix Matrix a
matrix) [a]
b) (Matrix a -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a]
matvec (Matrix a -> Matrix a
forall a. Matrix a -> Matrix a
transposeMatrix Matrix a
matrix) [a]
r)
           )

finite :: (RealFloat a) => a -> Bool
finite :: forall a. RealFloat a => a -> Bool
finite a
x = Bool -> Bool
not (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)

matvec :: (Num a) => Matrix a -> [a] -> [a]
matvec :: forall a. Num a => Matrix a -> [a] -> [a]
matvec Matrix a
a [a]
x = ([a] -> a) -> [[a]] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map ([a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
`dotV` [a]
x) (Matrix a -> [[a]]
forall a. Matrix a -> [[a]]
entries Matrix a
a)

residual :: (Num a) => Matrix a -> [a] -> [a] -> [a]
residual :: forall a. Num a => Matrix a -> [a] -> [a] -> [a]
residual Matrix a
a [a]
b [a]
x = [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
subV [a]
b (Matrix a -> [a] -> [a]
forall a. Num a => Matrix a -> [a] -> [a]
matvec Matrix a
a [a]
x)

scaleV :: (Num a) => a -> [a] -> [a]
scaleV :: forall a. Num a => a -> [a] -> [a]
scaleV a
a = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a
a a -> a -> a
forall a. Num a => a -> a -> a
*)

addV :: (Num a) => [a] -> [a] -> [a]
addV :: forall a. Num a => [a] -> [a] -> [a]
addV = (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
(+)

subV :: (Num a) => [a] -> [a] -> [a]
subV :: forall a. Num a => [a] -> [a] -> [a]
subV = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-)

normalize :: (RealFloat a) => [a] -> [a]
normalize :: forall a. RealFloat a => [a] -> [a]
normalize [a]
xs = let n :: a
n = [a] -> a
forall a. RealFloat a => [a] -> a
stableNorm [a]
xs in if a
n a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 then (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
forall a b. a -> b -> a
const a
0) [a]
xs else a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scaleV (a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
n) [a]
xs

combine :: (Num a) => [[a]] -> [a] -> [a]
combine :: forall a. Num a => [[a]] -> [a] -> [a]
combine [[a]]
basis [a]
ys = ([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]
forall a. Num a => [a] -> [a] -> [a]
addV (Int -> a -> [a]
forall a. Int -> a -> [a]
replicate Int
dimension a
0) ((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]
scaleV [a]
ys [[a]]
basis)
  where
    dimension :: Int
dimension = case [[a]]
basis of [a]
row : [[a]]
_ -> [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
row; [] -> Int
0

orthogonalize :: (RealFloat a) => [[a]] -> [a] -> ([a], [a])
orthogonalize :: forall a. RealFloat a => [[a]] -> [a] -> ([a], [a])
orthogonalize [[a]]
basis [a]
w0 =
    let ([a]
hs, [a]
w) = [a] -> ([a], [a])
pass [a]
w0; ([a]
correction, [a]
wn) = [a] -> ([a], [a])
pass [a]
w
     in ((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]
hs [a]
correction, [a]
wn)
  where
    pass :: [a] -> ([a], [a])
pass [a]
w =
        (([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]
hs, [a]
current) [a]
q ->
                let h :: a
h = [a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
dotV [a]
q [a]
current
                 in ([a]
hs [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a
h], [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
subV [a]
current (a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scaleV a
h [a]
q))
            )
            ([], [a]
w)
            [[a]]
basis

vectorValue :: (Num a) => [a] -> Int -> a
vectorValue :: forall a. Num a => [a] -> Int -> a
vectorValue [a]
xs Int
i = [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
1 (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
i [a]
xs))

vectorAt :: Int -> [[a]] -> [a]
vectorAt :: forall a. Int -> [[a]] -> [a]
vectorAt Int
i [[a]]
rows = [[a]] -> [a]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat (Int -> [[a]] -> [[a]]
forall a. Int -> [a] -> [a]
take Int
1 (Int -> [[a]] -> [[a]]
forall a. Int -> [a] -> [a]
drop Int
i [[a]]
rows))

value :: (Num a) => [[a]] -> Int -> Int -> a
value :: forall a. Num a => [[a]] -> Int -> Int -> a
value [[a]]
rows Int
i Int
j = [a] -> Int -> a
forall a. Num a => [a] -> Int -> a
vectorValue (Int -> [[a]] -> [a]
forall a. Int -> [[a]] -> [a]
vectorAt Int
i [[a]]
rows) Int
j

replaceAt :: Int -> a -> [a] -> [a]
replaceAt :: forall a. Int -> a -> [a] -> [a]
replaceAt Int
i a
x [a]
xs = Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
i [a]
xs [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
: Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [a]
xs

setValue :: Int -> Int -> a -> [[a]] -> [[a]]
setValue :: forall a. Int -> Int -> a -> [[a]] -> [[a]]
setValue Int
i Int
j a
x [[a]]
rows = Int -> [a] -> [[a]] -> [[a]]
forall a. Int -> a -> [a] -> [a]
replaceAt Int
i (Int -> a -> [a] -> [a]
forall a. Int -> a -> [a] -> [a]
replaceAt Int
j a
x (Int -> [[a]] -> [a]
forall a. Int -> [[a]] -> [a]
vectorAt Int
i [[a]]
rows)) [[a]]
rows