-- |
-- Copyright   : (c) Johannes Kropp
-- License     : BSD 3-Clause
-- Maintainer  : Johannes Kropp <jodak932@gmail.com>

module Math.Nuha.Algorithms where

import Prelude
import qualified Prelude as P
import Data.Vector.Unboxed (Unbox)
import Math.Nuha.Base
import Math.Nuha.Numeric
import Math.Nuha.Types

-- LITERATURE:
--  [1]: "Einführung in die Numerische Mathematik - Thomas Richter - 3.Auflage"

-------------
-- ** Solvers

{- | Solves the linear system A*x=b with A as a square matrix. The possible errors that can happen are 'NoSquareMatrixError', 'UnderdeterminedSystemError' and 'DimensionMismatchError'

> _A = matrix [[1,1,2],[2,-3,0],[2,4,-4]]
> b = vector [2,5,3]
> x = case solveLin _A b of
>     Left err -> error $ "solveLin : " ++ show err
>     Right x -> x
-}
solveLin :: (Unbox a, Real a, Floating a)
    => Holor a
    -> Holor a -> Either Error (Holor a)
solveLin :: Holor a -> Holor a -> Either Error (Holor a)
solveLin Holor a
_A Holor a
b = case Holor a -> Bool
forall a. Unbox a => Holor a -> Bool
isSquare Holor a
_A of
    Bool
False -> Error -> Either Error (Holor a)
forall a b. a -> Either a b
Left Error
NoSquareMatrixError
    Bool
True -> case Holor a -> Holor a -> Either Error (Holor a)
forall a.
(Unbox a, Real a, Floating a) =>
Holor a -> Holor a -> Either Error (Holor a)
solveLinLS Holor a
_A Holor a
b of
        Left Error
UnderdeterminedSystemError -> Error -> Either Error (Holor a)
forall a b. a -> Either a b
Left Error
UnderdeterminedSystemError
        Left Error
DimensionMismatchError -> Error -> Either Error (Holor a)
forall a b. a -> Either a b
Left Error
DimensionMismatchError
        Right Holor a
x -> Holor a -> Either Error (Holor a)
forall a b. b -> Either a b
Right Holor a
x

{-
With the QR-decomposition of A the solution can be found like this:
    A*x=b <=> Q*R*x=b <=> R*x=Q.T*b
The last equation can be solved easily by backward substitution.
We actually don't need the full QR-decomposition of A. It suffices if we have the householder reflection vectors from the incomplete QR decomposition. With that the right hand side of the equation can be expressed as:
    Q.T*b = S_(m-1)*S_(m-2)*...*S_1*b
For further explanation see [1], P.73
-}
{- | Solves the linear least squares problem, i.e. finds the least squares solution of the overdetermined linear equation system A*x=b. Possible errors are 'UnderdeterminedSystemError' and 'DimensionMismatchError'

> _A = matrix [[1,1,2],[2,-3,0],[2,4,-4],[2,4,-4.3]]
> b = vector [2,3,4,5]
> x = case solveLinLS _A b of
>     Left err -> error $ "solveLinLS : " ++ show err
>     Right x -> x
-}
solveLinLS :: (Unbox a, Real a, Floating a) => Holor a -> Holor a -> Either Error (Holor a)
solveLinLS :: Holor a -> Holor a -> Either Error (Holor a)
solveLinLS Holor a
_A Holor a
b = case Holor a -> Either Error ([Holor a], Holor a)
forall a.
(Unbox a, Real a, Floating a) =>
Holor a -> Either Error ([Holor a], Holor a)
facPreQR Holor a
_A of
    Left Error
TooFewRowsError -> Error -> Either Error (Holor a)
forall a b. a -> Either a b
Left Error
UnderdeterminedSystemError
    Left Error
NoMatrixError -> Error -> Either Error (Holor a)
forall a b. a -> Either a b
Left Error
DimensionMismatchError
    Right ([Holor a]
listV, Holor a
_R) -> Either Error (Holor a)
result where
        -- The product of a matrix S_j and a vector b can be expressed as a reflection with the corresponding housholder reflection vector v_j, see [1], S.71:
        reflection :: (Unbox a, Num a) => Holor a -> Holor a -> Holor a
        reflection :: Holor a -> Holor a -> Holor a
reflection Holor a
v Holor a
b = Holor a
b Holor a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => Holor a -> Holor a -> Holor a
|-| a
2a -> a -> a
forall a. Num a => a -> a -> a
*(Holor a -> Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> Holor a -> a
dot Holor a
v Holor a
b) a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => a -> Holor a -> Holor a
*| Holor a
v
        -- If the shapes of A and b match, the solution can be found by backward substitution in the System R*x=Q.T*b, see ([1], P.83):
        result :: Either Error (Holor a)
result = case Holor a -> [Int]
forall a. Holor a -> [Int]
hShape Holor a
b [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
== [Int
m,Int
1] of
            Bool
True -> case Holor a -> Holor a -> Either Error (Holor a)
forall a.
(Unbox a, Real a, Floating a) =>
Holor a -> Holor a -> Either Error (Holor a)
solveLinBack Holor a
_R Holor a
_QTb of
                Left Error
NoUpperTriError -> Error -> Either Error (Holor a)
forall a b. a -> Either a b
Left Error
UnderdeterminedSystemError
                Left Error
DimensionMismatchError -> Error -> Either Error (Holor a)
forall a b. a -> Either a b
Left Error
DimensionMismatchError
                Left Error
NoMatrixError -> Error -> Either Error (Holor a)
forall a b. a -> Either a b
Left Error
DimensionMismatchError
                Right Holor a
x -> Holor a -> Either Error (Holor a)
forall a b. b -> Either a b
Right Holor a
x
            Bool
False -> Error -> Either Error (Holor a)
forall a b. a -> Either a b
Left Error
DimensionMismatchError
            where
                m :: Int
m = [Int] -> Int
forall a. [a] -> a
head ([Int] -> Int) -> [Int] -> Int
forall a b. (a -> b) -> a -> b
$ Holor a -> [Int]
forall a. Holor a -> [Int]
hShape Holor a
_A
                _QTb :: Holor a
_QTb = (Holor a -> Holor a -> Holor a) -> Holor a -> [Holor a] -> Holor a
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr Holor a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => Holor a -> Holor a -> Holor a
reflection Holor a
b [Holor a]
listV


-- see [1], S.34
-- | Backward substitution for solving R*x=b with an upper triangular matrix R. Possible errors are 'DimensionMismatchError' and 'NoUpperTriError'
solveLinBack ::  (Unbox a, Real a, Floating a) => Holor a -> Holor a -> Either Error (Holor a)
solveLinBack :: Holor a -> Holor a -> Either Error (Holor a)
solveLinBack Holor a
_R Holor a
b = case Holor a -> [Int]
forall a. Holor a -> [Int]
hShape Holor a
_R of
    [Int
m,Int
n] -> case Holor a -> Bool
forall a. (Unbox a, Eq a, Num a) => Holor a -> Bool
isUpperTri Holor a
_R of
        Bool
False -> Error -> Either Error (Holor a)
forall a b. a -> Either a b
Left Error
NoUpperTriError
        Bool
True -> case Holor a -> [Int]
forall a. Holor a -> [Int]
hShape Holor a
b [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
== [Int
m,Int
1] of
            Bool
False -> Error -> Either Error (Holor a)
forall a b. a -> Either a b
Left Error
DimensionMismatchError
            Bool
True -> Holor a -> Either Error (Holor a)
forall a b. b -> Either a b
Right (Holor a -> Either Error (Holor a))
-> Holor a -> Either Error (Holor a)
forall a b. (a -> b) -> a -> b
$ Int -> [a] -> Holor a
recursion (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) []
            where
                recursion :: Int -> [a] -> Holor a
recursion Int
i [a]
xList = case Int
iInt -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
0 of
                    Bool
True -> [a] -> Holor a
forall a. Unbox a => [a] -> Holor a
vector [a]
xList'
                    Bool
False -> Int -> [a] -> Holor a
recursion (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) [a]
xList'
                    where
                        xList' :: [a]
xList' = a
xi a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a]
xList
                        xi :: a
xi = a
1a -> a -> a
forall a. Fractional a => a -> a -> a
/(Holor a
_RHolor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|![Int
i,Int
i]) a -> a -> a
forall a. Num a => a -> a -> a
*
                            (Holor a
bHolor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|![Int
i,Int
0] a -> a -> a
forall a. Num a => a -> a -> a
-
                                [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
P.sum ((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
(*) (Holor a -> [a]
forall a. Unbox a => Holor a -> [a]
toList (Holor a
_RHolor a -> [[Int]] -> Holor a
forall a. Unbox a => Holor a -> [[Int]] -> Holor a
|||![[Int
i],[Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1..Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]])) [a]
xList))
    [Int]
_ -> Error -> Either Error (Holor a)
forall a b. a -> Either a b
Left Error
DimensionMismatchError

--------------------
-- ** Factorizations

{- | Factorization of a rectangular matrix A with shape [m,n] where __m>=n__ into the matrices (Q,R) with A=Q*R where Q is orthogonal and R is of shape [m,n] with the first n rows in upper triangular form and the last m-n rows filled with zeros. Possible errors are 'TooFewRowsError' and 'NoMatrixError'

> _A = matrix [[1,1,2],[2,-3,0],[2,4,-4],[2,4,-5]]
> (_Q,_R) = case facQR _A of
>     Left err -> error $ "facQR : " ++ show err
>     Right (_Q,_R) -> (_Q,_R)
-}
facQR :: (Unbox a, Real a, Floating a) => Holor a -> Either Error (Holor a, Holor a)
facQR :: Holor a -> Either Error (Holor a, Holor a)
facQR Holor a
_A = case Holor a -> Either Error ([Holor a], Holor a)
forall a.
(Unbox a, Real a, Floating a) =>
Holor a -> Either Error ([Holor a], Holor a)
facPreQR Holor a
_A of
    Left Error
TooFewRowsError -> Error -> Either Error (Holor a, Holor a)
forall a b. a -> Either a b
Left Error
TooFewRowsError
    Left Error
NoMatrixError -> Error -> Either Error (Holor a, Holor a)
forall a b. a -> Either a b
Left Error
NoMatrixError
    Right ([Holor a]
listV, Holor a
_R) -> (Holor a, Holor a) -> Either Error (Holor a, Holor a)
forall a b. b -> Either a b
Right (Holor a
_Q,Holor a
_R) where
        _Q :: Holor a
_Q = Holor a -> Holor a
forall a. Unbox a => Holor a -> Holor a
transpose (Holor a -> Holor a) -> Holor a -> Holor a
forall a b. (a -> b) -> a -> b
$
            (Holor a -> Holor a -> Holor a) -> Holor a -> [Holor a] -> Holor a
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (\Holor a
v Holor a
_S -> (Int -> Holor a
forall a. (Unbox a, Num a) => Int -> Holor a
eye Int
m Holor a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => Holor a -> Holor a -> Holor a
|-| a
2a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => a -> Holor a -> Holor a
*| (Holor a
v Holor a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => Holor a -> Holor a -> Holor a
|.| Holor a -> Holor a
forall a. Unbox a => Holor a -> Holor a
transpose Holor a
v)) Holor a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => Holor a -> Holor a -> Holor a
|.| Holor a
_S) (Int -> Holor a
forall a. (Unbox a, Num a) => Int -> Holor a
eye Int
m) [Holor a]
listV
        [Int
m,Int
n] = Holor a -> [Int]
forall a. Holor a -> [Int]
hShape Holor a
_A


{- | Pre step of the QR factorization of a rectangular matrix A. The outputs of this function are a the list of the householder reflection vectors and the matrix R which is the same as in the full QR factorization. The matrix Q is implicitly stored in the householder reflection vectors by the rule: transpose(Q) = S_m * ... * S_1. For the matrices S_i applies: S_i = I - 2*(dot v_i v_i) where I is the Identity matrix and v_i the i-th householder vector. This function is sometimes useful where the Q Matrix isn't explicitly needed. Possible errors are 'TooFewRowsError' and 'NoMatrixError'

> _A = matrix [[1,1,2],[2,-3,0],[2,4,-4],[2,4,-5]]
> (reflectionVectors,_R) = case facPreQR _A of
>     Left err -> error $ "facPreQR : " ++ show err
>     Right (reflectionVectors,_R) -> (reflectionVectors,_R)
-}
facPreQR :: (Unbox a, Real a, Floating a) => Holor a -> Either Error ([Holor a], Holor a)
facPreQR :: Holor a -> Either Error ([Holor a], Holor a)
facPreQR Holor a
_A = case Holor a -> [Int]
forall a. Holor a -> [Int]
hShape Holor a
_A of
    [Int
m,Int
n] -> Either Error ([Holor a], Holor a)
guards where
        guards :: Either Error ([Holor a], Holor a)
guards
            | Int
mInt -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<Int
n = Error -> Either Error ([Holor a], Holor a)
forall a b. a -> Either a b
Left Error
TooFewRowsError
            | Int
mInt -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
n = ([Holor a], Holor a) -> Either Error ([Holor a], Holor a)
forall a b. b -> Either a b
Right ([Holor a]
listV, Holor a
_R)
            | Bool
otherwise = ([Holor a], Holor a) -> Either Error ([Holor a], Holor a)
forall a b. b -> Either a b
Right ([Holor a]
listV', Holor a
_R')
        -- For a square input matrix the recursion is called without an additional recursion step:
        ([Holor a]
listV, [[a]]
listR) = (Int, Int)
-> Bool
-> [Holor a]
-> [[a]]
-> Maybe (Holor a)
-> ([Holor a], [[a]])
forall a.
(Unbox a, Real a, Floating a) =>
(Int, Int)
-> Bool
-> [Holor a]
-> [[a]]
-> Maybe (Holor a)
-> ([Holor a], [[a]])
recursion (Int
m,Int
n) Bool
False [] [] (Holor a -> Maybe (Holor a)
forall a. a -> Maybe a
Just Holor a
_A)
        _R :: Holor a
_R = [[a]] -> Holor a
forall a. Unbox a => [[a]] -> Holor a
matrix ([[a]] -> Holor a) -> [[a]] -> Holor a
forall a b. (a -> b) -> a -> b
$ [[a]] -> [[a]]
forall a. [a] -> [a]
reverse [[a]]
listR
        -- For a rectangular input matrix there must me an additional recursion step. Notice that in this case there must be additional rows with zeros in the R-Matrix to go conform with the shape of Q:
        ([Holor a]
listV', [[a]]
listR') = (Int, Int)
-> Bool
-> [Holor a]
-> [[a]]
-> Maybe (Holor a)
-> ([Holor a], [[a]])
forall a.
(Unbox a, Real a, Floating a) =>
(Int, Int)
-> Bool
-> [Holor a]
-> [[a]]
-> Maybe (Holor a)
-> ([Holor a], [[a]])
recursion (Int
m,Int
n) Bool
True [] [] (Holor a -> Maybe (Holor a)
forall a. a -> Maybe a
Just Holor a
_A)
        _R' :: Holor a
_R' = [[a]] -> Holor a
forall a. Unbox a => [[a]] -> Holor a
matrix ([[a]] -> Holor a) -> [[a]] -> Holor a
forall a b. (a -> b) -> a -> b
$ [[a]] -> [[a]]
forall a. [a] -> [a]
reverse (Int -> [[a]] -> [[a]]
forall a. Int -> [a] -> [a]
take (Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
n) ([a] -> [[a]]
forall a. a -> [a]
repeat (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
n (a -> [a]
forall a. a -> [a]
repeat a
0))) [[a]] -> [[a]] -> [[a]]
forall a. [a] -> [a] -> [a]
++ [[a]]
listR')
    [Int]
_ -> Error -> Either Error ([Holor a], Holor a)
forall a b. a -> Either a b
Left Error
NoMatrixError
    where
        -- recursion to calculate the householder vectors and the elements of the R-Matrix, see [1], S.70ff. The parameter shapeIn is the shape of the input Matrix A. The parameter additionalStep determines if there should be an additional recursion step, which is the case for rectangular input matrices with more rows than columns:
        recursion :: (Unbox a, Real a, Floating a) =>
            (Int,Int) -> Bool -> [Holor a] -> [[a]] -> Maybe (Holor a) -> ([Holor a], [[a]])
        recursion :: (Int, Int)
-> Bool
-> [Holor a]
-> [[a]]
-> Maybe (Holor a)
-> ([Holor a], [[a]])
recursion (Int, Int)
shapeIn Bool
additionalStep [Holor a]
listV [[a]]
listR Maybe (Holor a)
maybeA = case Maybe (Holor a)
maybeA of
            Maybe (Holor a)
Nothing -> ([Holor a]
listV, [[a]]
listR)
            Just Holor a
_A -> (Int, Int)
-> Bool
-> [Holor a]
-> [[a]]
-> Maybe (Holor a)
-> ([Holor a], [[a]])
forall a.
(Unbox a, Real a, Floating a) =>
(Int, Int)
-> Bool
-> [Holor a]
-> [[a]]
-> Maybe (Holor a)
-> ([Holor a], [[a]])
recursion (Int, Int)
shapeIn Bool
additionalStep' [Holor a]
listV' [[a]]
listR' Maybe (Holor a)
maybeA' where
                -- shapes of the Input matrix and the submatrix of the current recursion step
                (Int
rows,Int
cols) = (Int, Int)
shapeIn
                [Int
m,Int
n] = Holor a -> [Int]
forall a. Holor a -> [Int]
hShape Holor a
_A
                -- top left element and first colomn of _A:
                a00 :: a
a00 = Holor a
_AHolor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
0,Int
0]
                a_0 :: Holor a
a_0 = Holor a
_AHolor a -> [[Int]] -> Holor a
forall a. Unbox a => Holor a -> [[Int]] -> Holor a
|||!![[Int
0..Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1],[Int
0]]
                -- next reflection vector, notice that the standard signum function can be zero, so a customized one must be used which can only bei -1 or 1:
                sign :: a -> p
sign a
x = if a
xa -> a -> Bool
forall a. Ord a => a -> a -> Bool
<a
0 then -p
1 else p
1
                v :: Holor a
v = Holor a -> Holor a
forall a. (Unbox a, Real a, Floating a) => Holor a -> Holor a
normalize2 (Holor a -> Holor a) -> Holor a -> Holor a
forall a b. (a -> b) -> a -> b
$ Holor a
a_0 Holor a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => Holor a -> Holor a -> Holor a
|+| [a] -> Holor a
forall a. Unbox a => [a] -> Holor a
vector ((a -> a
forall a p. (Ord a, Num a, Num p) => a -> p
sign a
a00)a -> a -> a
forall a. Num a => a -> a -> a
*(Holor a -> a
forall a b. (Unbox a, Real a, Unbox b, Floating b) => Holor a -> b
norm2 Holor a
a_0) a -> [a] -> [a]
forall a. a -> [a] -> [a]
: (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take (Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) (a -> [a]
forall a. a -> [a]
repeat a
0)))
                -- append leading zeros to the reflection vector to make sure all reflection vectors are of the same length and add it to the list of the existing ones:
                listV' :: [Holor a]
listV' = ([a] -> Holor a
forall a. Unbox a => [a] -> Holor a
vector (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take (Int
rows Int -> Int -> Int
forall a. Num a => a -> a -> a
- Holor a -> Int
forall a. Unbox a => Holor a -> Int
numElems Holor a
v) (a -> [a]
forall a. a -> [a]
repeat a
0) [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ Holor a -> [a]
forall a. Unbox a => Holor a -> [a]
toList Holor a
v)) Holor a -> [Holor a] -> [Holor a]
forall a. a -> [a] -> [a]
: [Holor a]
listV
                -- calculate the new list with the rows of the R-Matrix and maybe a new submatrix A'. If there is no next submatrix A', the recursion will stop at the begin of the next recursion:
                i :: Int
i = [Holor a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Holor a]
listV -- recursion step
                ([[a]]
listR', Maybe (Holor a)
maybeA', Bool
additionalStep') = case Int
iInt -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<Int
colsInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1 of
                    -- Do a standard recursion step. A new row of R is appended as the first row of S*A with leading zeros and the new submatrix A' is taken from S*A by skipping the first row and first column, see [1], P.71:
                    Bool
True ->
                        ( (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
i (a -> [a]
forall a. a -> [a]
repeat a
0) [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ (a
_SA00 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: ([[a]] -> [a]
forall a. [a] -> a
head [[a]]
_SA_withoutFirstColumn))) [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
: [[a]]
listR
                        , Holor a -> Maybe (Holor a)
forall a. a -> Maybe a
Just (Holor a -> Maybe (Holor a)) -> Holor a -> Maybe (Holor a)
forall a b. (a -> b) -> a -> b
$ [[a]] -> Holor a
forall a. Unbox a => [[a]] -> Holor a
matrix ([[a]] -> [[a]]
forall a. [a] -> [a]
tail [[a]]
_SA_withoutFirstColumn)
                        , Bool
additionalStep )
                    -- The last row of R is calculated different, depending if there should be an additional step at the end:
                    Bool
False -> case Bool
additionalStep of
                        -- If there should be no additional step the new row of R is calculated as the first row of S*A with leading zeros (which has only the last element as nonzero). There will be no next submatrix A':
                        Bool
False ->
                            ( (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take (Int
colsInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) (a -> [a]
forall a. a -> [a]
repeat a
0) [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a
_SA00]) [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
: [[a]]
listR
                            , Maybe (Holor a)
forall a. Maybe a
Nothing
                            , Bool
False )
                        -- If there should be an additional step, the matrix R doesn't change, but there is a next special submatrix A', see ([1], P.82-84). In addition the parameter additionalStep is set to False to finish at the next recursion:
                        Bool
True ->
                            ( [[a]]
listR
                            , Holor a -> Maybe (Holor a)
forall a. a -> Maybe a
Just (Holor a -> Maybe (Holor a)) -> Holor a -> Maybe (Holor a)
forall a b. (a -> b) -> a -> b
$ Holor a
_SA_firstColumnHolor
                            , Bool
False )
                    where
                        -- first element of the first column of S*A:
                        _SA00 :: a
_SA00 = - (a -> a
forall a p. (Ord a, Num a, Num p) => a -> p
sign a
a00) a -> a -> a
forall a. Num a => a -> a -> a
* (Holor a -> a
forall a b. (Unbox a, Real a, Unbox b, Floating b) => Holor a -> b
norm2 Holor a
a_0)
                        -- second to last column of S*A (from index j = 1 to n-1) as a nested list of rows:
                        _SA_withoutFirstColumn :: [[a]]
_SA_withoutFirstColumn = Holor a -> [[a]]
forall a. Unbox a => Holor a -> [[a]]
toList2 (Holor a -> [[a]]) -> Holor a -> [[a]]
forall a b. (a -> b) -> a -> b
$ Holor a -> Holor a
forall a. Unbox a => Holor a -> Holor a
transpose (Holor a -> Holor a) -> Holor a -> Holor a
forall a b. (a -> b) -> a -> b
$ [[a]] -> Holor a
forall a. Unbox a => [[a]] -> Holor a
matrix [
                            let
                                a :: Holor a
a = Holor a
_AHolor a -> [[Int]] -> Holor a
forall a. Unbox a => Holor a -> [[Int]] -> Holor a
|||!![[Int
0..Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1],[Int
j]]
                            in
                                Holor a -> [a]
forall a. Unbox a => Holor a -> [a]
toList (Holor a -> [a]) -> Holor a -> [a]
forall a b. (a -> b) -> a -> b
$ Holor a
a Holor a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => Holor a -> Holor a -> Holor a
|-| a
2a -> a -> a
forall a. Num a => a -> a -> a
*(Holor a -> Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> Holor a -> a
dot Holor a
v Holor a
a) a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => a -> Holor a -> Holor a
*| Holor a
v
                            | Int
j <- [Int
1..Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]]
                        -- new submatrix for the additional recursion step:
                        _SA_firstColumnHolor :: Holor a
_SA_firstColumnHolor = Holor a
a Holor a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => Holor a -> Holor a -> Holor a
|-| a
2a -> a -> a
forall a. Num a => a -> a -> a
*(Holor a -> Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> Holor a -> a
dot Holor a
v Holor a
a) a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => a -> Holor a -> Holor a
*| Holor a
v where
                            a :: Holor a
a = Holor a
_AHolor a -> [[Int]] -> Holor a
forall a. Unbox a => Holor a -> [[Int]] -> Holor a
|||!![[Int
0..Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1],[Int
0]]