Numeric.LinearAlgebra.Sparse

Matrix factorizations

qr :: (Epsilon a, Ord a, Floating a) => SpMatrix a -> (SpMatrix a, SpMatrix a)

Given a matrix A, returns a pair of matrices (Q, R) such that Q R = A, where Q is orthogonal and R is upper triangular. Applies Givens rotation iteratively to zero out sub-diagonal elements.

lu :: (Epsilon a, Fractional a, Real a) => SpMatrix a -> (SpMatrix a, SpMatrix a)

Given a matrix A, returns a pair of matrices (L, U) where L is lower triangular and U is upper triangular such that L U = A

chol :: (Epsilon a, Real a, Floating a) => SpMatrix a -> SpMatrix a

Given a positive semidefinite matrix A, returns a lower-triangular matrix L such that L L^T = A

Condition number

conditionNumberSM :: (Epsilon a, RealFloat a) => SpMatrix a -> a

uses the R matrix from the QR factorization

Householder reflection

hhMat :: Num a => a -> SpVector a -> SpMatrix a

hhRefl :: Num a => SpVector a -> SpMatrix a

a vector x uniquely defines an orthogonal plane; the Householder operator reflects any point v with respect to this plane: v' = (I - 2 x >< x) v

Givens' rotation

givens :: (Floating a, Epsilon a, Ord a) => SpMatrix a -> IxRow -> IxCol -> SpMatrix a

Givens method, row version: choose other row index i' s.t. i' is : * below the diagonal * corresponding element is nonzero

QR.C1 ) To zero out entry A(i, j) we must find row k such that A(k, j) is non-zero but A has zeros in row k for all columns less than j.

NB: the current version is quite inefficient in that: 1. the Givens' matrix G_i is different from Identity only in 4 entries 2. at each iteration i we multiply G_i by the previous partial result M. Since this corresponds to a rotation, and the givensCoef function already computes the value of the resulting non-zero component (output r), `G_i ## M` can be simplified by just changing two entries of M (i.e. zeroing one out and changing the other into r).

Arnoldi iteration

arnoldi :: (Epsilon a, Floating a, Eq a) => SpMatrix a -> SpVector a -> Int -> (SpMatrix a, SpMatrix a)

Given a matrix A, a vector b and a positive integer n, this procedure finds the basis of an order n Krylov subspace (as the columns of matrix Q), along with an upper Hessenberg matrix H, such that A = Q^T H Q. At the i`th iteration, it finds (i + 1) coefficients (the i`th column of the Hessenberg matrix H) and the (i + 1)`th Krylov vector.

Eigensolvers

eigsQR :: (Epsilon a, Real a, Floating a) => Int -> SpMatrix a -> SpVector a

`eigsQR n mm` performs n iterations of the QR algorithm on matrix mm, and returns a SpVector containing all eigenvalues

eigRayleigh :: Int -> SpMatrix Double -> (SpVector Double, Double) -> (SpVector Double, Double)

`eigsRayleigh n mm` performs n iterations of the Rayleigh algorithm on matrix mm and returns the eigenpair closest to the initialization. It displays cubic-order convergence, but it also requires an educated guess on the initial eigenpair eigRayleigh :: Int -- max # iterations -> SpMatrix Double -- matrix -> (SpVector Double, Double) -- initial guess of (eigenvector, eigenvalue) -> (SpVector Double, Double) -- final estimate of (eigenvector, eigenvalue)

Linear solvers

Iterative methods

linSolve :: LinSolveMethod -> SpMatrix Double -> SpVector Double -> SpVector Double

Linear solve with _deterministic_ starting vector (every component at 0.1)

data LinSolveMethod

Constructors

GMRES_

CGNE_

TFQMR_

BCG_

CGS_

BICGSTAB_

Instances

Eq LinSolveMethod

Methods

(==) :: LinSolveMethod -> LinSolveMethod -> Bool

(/=) :: LinSolveMethod -> LinSolveMethod -> Bool

Show LinSolveMethod

Methods

showsPrec :: Int -> LinSolveMethod -> ShowS

show :: LinSolveMethod -> String

showList :: [LinSolveMethod] -> ShowS

(<\>) :: SpMatrix Double -> SpVector Double -> SpVector Double

\ : linSolve using the BiCGSTAB method as default

pinv :: SpMatrix Double -> SpVector Double -> SpVector Double

Least-squares approximation of a rectangular system of equaitons. Uses \ for the linear solve

Direct methods

luSolve :: (Fractional a, Eq a, Epsilon a) => SpMatrix a -> SpMatrix a -> SpVector a -> SpVector a

Direct solver based on a triangular factorization of the system matrix.

Preconditioners

ilu0 :: (Epsilon a, Real a, Fractional a) => SpMatrix a -> (SpMatrix a, SpMatrix a)

Used for Incomplete LU : remove entries in m corresponding to zero entries in m2

mSsor :: Fractional a => SpMatrix a -> a -> (SpMatrix a, SpMatrix a)

`mSsor aa omega` : if `omega = 1` it returns the symmetric Gauss-Seidel preconditioner

Matrix partitioning

diagPartitions :: SpMatrix a -> (SpMatrix a, SpMatrix a, SpMatrix a)

Partition a matrix into strictly subdiagonal, diagonal and strictly superdiagonal parts

Random arrays

randArray :: PrimMonad m => Int -> Double -> Double -> m [Double]

Random matrices and vectors

randMat :: PrimMonad m => Int -> m (SpMatrix Double)

Dense SpMatrix

randVec :: PrimMonad m => Int -> m (SpVector Double)

Dense SpVector

Sparse "

randSpMat :: Int -> Int -> IO (SpMatrix Double)

Sparse SpMatrix

randSpVec :: Int -> Int -> IO (SpVector Double)

Sparse SpVector

Sparsify data

sparsifySV :: Epsilon a => SpVector a -> SpVector a

Sparsify an SpVector

Iteration combinators

modifyInspectN :: MonadState s m => Int -> ([s] -> Bool) -> (s -> s) -> m s

Keep a moving window buffer (length 2) of state x to assess convergence, stop when either a condition on that list is satisfied or when max # of iterations is reached (i.e. same thing as loopUntilAcc but this one runs in the State monad)

runAppendN' :: (t -> t) -> Int -> t -> [t]

", NO convergence check

untilConverged :: MonadState a m => (a -> SpVector Double) -> (a -> a) -> m a

iterate until convergence is verified or we run out of a fixed iteration budget

diffSqL :: Floating a => [a] -> a