sparse-linear-algebra-0.2.2.0: Numerical computation in native Haskell

Safe Haskell None Haskell2010

Numeric.LinearAlgebra.Sparse

Synopsis

# Matrix factorizations

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

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) Source #

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 Source #

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 Source #

uses the R matrix from the QR factorization

# Householder reflection

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

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

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 Source #

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) Source #

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 ith iteration, it finds (i + 1) coefficients (the ith 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 Source #

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

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

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

Constructors

 GMRES_ CGNE_ TFQMR_ BCG_ CGS_ BICGSTAB_

Instances

 Source # Methods Source # MethodsshowList :: [LinSolveMethod] -> ShowS #

\ : linSolve using the BiCGSTAB method as default

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 Source #

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) Source #

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

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

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

# Matrix partitioning

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

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

# Random matrices and vectors

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

Dense SpMatrix

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

Dense SpVector

Sparse SpMatrix

Sparse SpVector

# Sparsify data

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

Sparsify an SpVector

# Iteration combinators

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

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] Source #

", NO convergence check

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

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

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