# README
Gosl. la. Linear Algebra: vector, matrix, efficient sparse solvers, eigenvalues, decompositions, etc.
The la sub-package implements functions to perform linear algebra operations such as vector
addition or matrix-vector multiplications. It also wraps some efficient sparse linear systems
solvers such as Umfpack and
MUMPS (not the parotitis disease!).
Both Umfpack and MUMPS solvers are very efficient!
Sometimes, we call the lower level functions in la/oblas to improve performance.
Structures for sparse problems
In la, there are two types of structures to hold data for solving a sparse linear system:
triplet and column-compressed matrix. Both have a complex version and are named as follows:
TripletandTripletC(complex)CCMatrixandCCMatrixC(complex)
The Triplet is more convenient for data input; e.g. during the assemblage step in a finite element
code. The CCMatrix is better (faster) for computations and is the structure used by Umfpack.
Triplets are initialized by giving the size of the corresponding matrix and the number of non-zero
entries (components). Thus, only space for the non-zero entries is allocated. Afterwards, each
component can be set by calling the Put method after calling the Start method. The Start
method can be called several times. Therefore, this structure helps with the efficient
implementation of codes requiring multiple solutions of linear systems as in finite element
analyses.
To convert the Triplet into a Column-Compressed Matrix, call the ToMatrix method of Triplet. And
to convert the sparse matrix to a dense version (e.g. for reporting/printing), call the ToDense
method of CCMatrix.
The Triplet has also a very convenient method to copy the contents of another Triplet (i.e. sparse
matrix) into the positions starting with the maximum number of columns of this second matrix and the
positions starting with the maximum number of rows of this second matrix. This is done with the
method PutMatAndMatT. This operation is particularly common when solving mechanical problems with
Lagrange multipliers and can be illustrated as follows:
[... ... ... a00 a10 ...] 0
[... ... ... a01 a11 ...] 1
[... ... ... a02 a12 ...] 2 [. at .]
[a00 a01 a02 ... ... ...] 3 => [a . .]
[a10 a11 a12 ... ... ...] 4 [. . .]
[... ... ... ... ... ...] 5
The version in which the second matrix is a column-compressed matrix is named PutCCMatAndMatT.
Linear solvers for sparse problems
SparseSolver defines an interface for linear solvers in la. Two implementations satisfying this
interface are:
Umfpackwrapper to Umfpack; andMumpswrapper to MUMPS
There are also high level functions to solve linear systems with Umfpack:
SpSolve; andSpSolveCwith complex numbers
Note however that the high level functions shouldn't be used for repeated executions because memory would be constantly allocated and released.
Examples
Vectors and matrices
- source file Test Vector
- source file Test Matrix
- source file Test Matrix operations
BLAS1, 2 and 3 functions
- source file Test BLAS1 routines
- source file Test BLAS2 routines
- source file Test BLAS3 routines
General dense solver and Cholesky decomposition
- source file Test Dense Solver
Eigenvalues and eigenvectors of general matrix
- source file Test Eigenvalues/Eigenvectors
Eigenvalues of symmetric (3 x 3) matrix
- source file Test Jacobi iteration
Sparse BLAS functions
- source file Test sparse BLAS
Conversions related to sparse matrices
- source file Test conversions
Sparse Triplet and Matrix
- source file Test sparse Triplet and Matrix
Sparse linear solver using MUMPS
- source file Test sparse solver MUMPS
Sparse linear solver using UMFPACK
- source file Test sparse solver UMFPACK
Solutions using sparse solvers
- source file Test solutions of sparse linear systems
API
# Functions
CheckEigenVecL checks left eigenvector:
H H u [j] ⋅ A = λ[j] ⋅ u [j] LEFT eigenvectors
.
CheckEigenVecR checks right eigenvector:
A ⋅ v[j] = λ[j] ⋅ v[j] RIGHT eigenvectors
.
Cholesky returns the Cholesky decomposition of a symmetric positive-definite matrix
a = L * trans(L)
.
DenSolve solves dense linear system using LAPACK (OpenBLAS)
Given: A ⋅ x = b find x such that x = A⁻¹ ⋅ b
.
EigenVal computes eigenvalues of general matrix
A ⋅ v[j] = λ[j] ⋅ v[j]
INPUT: a -- general matrix
OUTPUT: w -- eigenvalues [pre-allocated]
.
EigenVecL computes eigenvalues and LEFT eigenvectors of general matrix
H H u [j] ⋅ A = λ[j] ⋅ u [j] LEFT eigenvectors
INPUT: a -- general matrix
OUTPUT: u -- matrix with the eigenvectors; each column contains one eigenvector [pre-allocated] w -- eigenvalues [pre-allocated]
.
EigenVecLR computes eigenvalues and LEFT and RIGHT eigenvectors of general matrix
A ⋅ v[j] = λ[j] ⋅ v[j] RIGHT eigenvectors
H H u [j] ⋅ A = λ[j] ⋅ u [j] LEFT eigenvectors
INPUT: a -- general matrix
OUTPUT: u -- matrix with the LEFT eigenvectors; each column contains one eigenvector [pre-allocated] v -- matrix with the RIGHT eigenvectors; each column contains one eigenvector [pre-allocated] w -- λ eigenvalues [pre-allocated]
.
EigenVecR computes eigenvalues and RIGHT eigenvectors of general matrix
A ⋅ v[j] = λ[j] ⋅ v[j]
INPUT: a -- general matrix
OUTPUT: v -- matrix with the eigenvectors; each column contains one eigenvector [pre-allocated] w -- eigenvalues [pre-allocated]
.
Jacobi performs the Jacobi transformation of a symmetric matrix to find its eigenvectors and eigenvalues.
MatAdd adds the scaled components of two matrices res := α⋅a + β⋅b ⇒ result[i][j] := α⋅a[i][j] + β⋅b[i][j].
MatCondNum returns the condition number of a square matrix using the inverse of this matrix; thus it is not as efficient as it could be, e.g.
MatInv computes the inverse of a general matrix (square or not).
MatInvSmall computes the inverse of small matrices of size 1x1, 2x2, or 3x3.
MatMatMul returns the matrix multiplication (scaled)
c := α⋅a⋅b ⇒ cij := α * aik * bkj
.
MatMatMulAdd returns the matrix multiplication (scaled)
c += α⋅a⋅b ⇒ cij += α * aik * bkj
.
MatMatTrMul returns the matrix multiplication (scaled) with transposed(b)
c := α⋅a⋅bᵀ ⇒ cij := α * aik * bjk
.
MatMatTrMulAdd returns the matrix multiplication (scaled) with transposed(b)
c += α⋅a⋅bᵀ ⇒ cij += α * aik * bjk
.
MatSvd performs the SVD decomposition Input: a -- matrix a copyA -- creates a copy of a; otherwise 'a' is modified Output: s -- diagonal terms [must be pre-allocated] len(s) = imin(a.M, a.N) u -- left matrix [must be pre-allocated] u is (a.M x a.M) vt -- transposed right matrix [must be pre-allocated] vt is (a.N x a.N).
MatTrMatMul returns the matrix multiplication (scaled) with transposed(a)
c := α⋅aᵀ⋅b ⇒ cij := α * aki * bkj
.
MatTrMatMulAdd returns the matrix multiplication (scaled) with transposed(a)
c += α⋅aᵀ⋅b ⇒ cij += α * aki * bkj
.
MatTrMatTrMul returns the matrix multiplication (scaled) with transposed(a) and transposed(b)
c := α⋅aᵀ⋅bᵀ ⇒ cij := α * aki * bjk
.
MatTrMatTrMulAdd returns the matrix multiplication (scaled) with transposed(a) and transposed(b)
c += α⋅aᵀ⋅bᵀ ⇒ cij += α * aki * bjk
.
MatTrVecMul returns the transpose(matrix)-vector multiplication
v = α⋅aᵀ⋅u ⇒ vi = α * aji * uj = α * uj * aji
.
MatVecMul returns the matrix-vector multiplication
v = α⋅a⋅u ⇒ vi = α * aij * uj
.
MatVecMulAdd returns the matrix-vector multiplication with addition
v += α⋅a⋅u ⇒ vi += α * aij * uj
.
MatVecMulAddC returns the matrix-vector multiplication with addition (complex version)
v += α⋅a⋅u ⇒ vi += α * aij * uj
.
MatVecMulC returns the matrix-vector multiplication (complex version)
v = α⋅a⋅u ⇒ vi = α * aij * uj
.
NewEquations creates a new Equations structure n -- total number of equations; i.e.
NewMatrix allocates a new (empty) Matrix with given (m,n) (row/col sizes).
NewMatrixC allocates a new (empty) MatrixC with given (m,n) (row/col sizes).
NewMatrixDeep2 allocates a new Matrix from given (Deep2) nested slice.
NewMatrixDeep2c allocates a new MatrixC from given (Deep2c) nested slice.
NewMatrixRaw creates a new Matrix using given raw data Input: rawdata -- data organized as column-major; e.g.
NewSparseConfig returns a new SparseConfig Input: comm -- may be nil.
NewSparseSolver finds a SparseSolver in database or panic kind -- "umfpack" or "mumps" NOTE: remember to call Free() to release allocated resources.
NewSparseSolverC finds a SparseSolver in database or panic NOTE: remember to call Free() to release allocated resources.
NewTriplet returns a new Triplet.
NewTripletC returns a new TripletC.
NewVector returns a new vector with size m.
NewVectorC returns a new vector with size m.
NewVectorMapped returns a new vector after applying a function over all of its components new: vi = f(i).
NewVectorMappedC returns a new vector after applying a function over all of its components new: vi = f(i).
NewVectorSlice returns a new vector from given Slice NOTE: This is equivalent to cast a slice to Vector as in: v := la.Vector([]float64{1,2,3}).
SolveRealLinSysSPD solves a linear system with real numbers and a Symmetric-Positive-Definite (SPD) matrix
x := inv(a) * b
NOTE: this function uses Cholesky decomposition and should be used for small systems.
SolveTwoRealLinSysSPD solves two linear systems with real numbers and Symmetric-Positive-Definite (SPD) matrices
x := inv(a) * b and X := inv(a) * B
NOTE: this function uses Cholesky decomposition and should be used for small systems.
SpAllocMatAddMat allocates a matrix 'c' to hold the result of the addition of 'a' and 'b'.
SpCheckDiag checks if all elements on the diagonal of "a" are present.
SpInitRc initializes two complex sparse matrices (residual correction) according to: Real: R := γ *I - J Complex: C := (α + βi)*I - J NOTE: "J" must include all diagonal elements.
SpInitSimilar initializes another matrix "b" with the same structure (Ap, Ai) of sparse matrix "a".
SpInitSimilarR2C initializes another matrix "b" (complex) with the same structure (Ap, Ai) of sparse matrix "a" (real).
SpMatAddI adds an identity matrix I to "a", scaled by α and β according to: r := α*a + β*I.
SpMatAddMat adds two sparse matrices.
SpMatAddMatC adds two real sparse matrices with two sets of coefficients in such a way that one real matrix (R) and another complex matrix (C) are obtained.
SpMatMatTrMul computes the multiplication of sparse matrix with itself transposed: b := α * a * aᵀ b_iq = α * a_ij * a_qj Note: b is symmetric.
SpMatTrVecMul returns the (sparse) matrix-vector multiplication with "a" transposed (scaled): v := α * transp(a) * u => vj = α * aij * ui NOTE: dense vector v will be first initialized with zeros.
SpMatTrVecMulAdd returns the (sparse) matrix-vector multiplication with addition and "a" transposed (scaled): v += α * transp(a) * u => vj += α * aij * ui.
SpMatTrVecMulAddC returns the (sparse/complex) matrix-vector multiplication with addition and "a" transposed (scaled): v += α * transp(a) * u => vj += α * aij * ui.
SpMatTrVecMulC returns the (sparse/complex) matrix-vector multiplication with "a" transposed (scaled): v := α * transp(a) * u => vj = α * aij * ui NOTE: dense vector v will be first initialized with zeros.
SpMatVecMul returns the (sparse) matrix-vector multiplication (scaled): v := α * a * u => vi = α * aij * uj NOTE: dense vector v will be first initialized with zeros.
SpMatVecMulAdd returns the (sparse) matrix-vector multiplication with addition (scaled): v += α * a * u => vi += α * aij * uj.
SpMatVecMulAddC returns the (sparse/complex) matrix-vector multiplication with addition (scaled): v += α * a * u => vi += α * aij * uj.
SpMatVecMulAddX returns the (sparse) matrix-vector multiplication with addition (scaled/extended): v += a * (α*u + β*w) => vi += aij * (α*uj + β*wj).
SpMatVecMulC returns the (sparse/complex) matrix-vector multiplication (scaled): v := α * a * u => vi = α * aij * uj NOTE: dense vector v will be first initialized with zeros.
SpSetRc sets the values within two complex sparse matrices (residual correction) according to: Real: R := γ *I - J Complex: C := (α + βi)*I - J NOTE: "J" must include all diagonal elements.
SpSolve solves a sparse linear system (using UMFPACK)
Given: A ⋅ x = b find x such that x = A⁻¹ ⋅ b
.
SpSolveC solves a sparse linear system (using UMFPACK) (complex version)
Given: A ⋅ x = b find x such that x = A⁻¹ ⋅ b
.
SpTriAdd adds two matrices in Triplet format: c := α*a + β*b NOTE: the output 'c' triplet must be able to hold all nonzeros of 'a' and 'b' actually the 'c' triplet is just expanded.
SpTriAddR2C adds two real matrices in Triplet format generating a complex triplet according to: c := (α+βi)*a + μ*b NOTE: the output 'c' triplet must be able to hold all nonzeros of 'a' and 'b'.
SpTriMatTrVecMul returns the matrix-vector multiplication with transposed matrix a in triplet format and two dense vectors x and y y := transpose(a) * x or y_I := a_JI * x_J or y_j := a_ij * x_i.
SpTriMatVecMul returns the matrix-vector multiplication with matrix a in triplet format and two dense vectors x and y y := a * x or y_i := a_ij * x_j.
SpTriSetDiag sets a (n x n) real triplet with diagonal values 'v'.
TestSolverResidual check the residual of a linear system solution.
TestSolverResidualC check the residual of a linear system solution (complex version).
TestSpSolver tests a sparse solver.
TestSpSolverC tests a sparse solver (complex version).
VecAdd adds the scaled components of two vectors res := α⋅u + β⋅v ⇒ result[i] := α⋅u[i] + β⋅v[i].
VecDot returns the dot product between two vectors: s := u・v.
VecMaxDiff returns the maximum absolute difference between two vectors maxdiff = max(|u - v|).
VecRmsError returns the scaled root-mean-square of the difference between two vectors with components normalized by a scaling factor __________________________ / ———— 2 / 1 \ / error[i] \ rms = \ / ——— / | —————————— | \/ N ———— \ scale[i] /
error[i] = |u[i] - v[i]|
scale[i] = a + m*|s[i]|
.
VecScaleAbs creates a "scale" vector using the absolute value of another vector scale := a + m ⋅ |x| ⇒ scale[i] := a + m ⋅ |x[i]|.
VecVecTrMul returns the matrix = vector-transpose(vector) multiplication (e.g.
# Structs
CCMatrix represents a sparse matrix using the so-called "column-compressed format".
CCMatrixC represents a sparse matrix using the so-called "column-compressed format".
Equations organises the identification numbers (IDs) of equations in a linear system: [A] ⋅ {x} = {b}.
Matrix implements a column-major representation of a matrix by using a linear array that can be passed to Fortran code
NOTE: the functions related to Matrix do not check for the limits of indices and dimensions.
MatrixC implements a column-major representation of a matrix of complex numbers by using a linear array that can be passed to Fortran code.
The SparseConfig structure holds configuration arguments for sparse solvers.
Triplet is a simple representation of a sparse matrix, where the indices and values of this matrix are stored directly.
TripletC is a simple representation of a sparse matrix, where the indices and values of this matrix are stored directly.
# Interfaces
SparseSolver solves sparse linear systems using UMFPACK or MUMPS
Given: A ⋅ x = b find x such that x = A⁻¹ ⋅ b
.
SparseSolverC solves sparse linear systems using UMFPACK or MUMPS (complex version)
Given: A ⋅ x = b find x such that x = A⁻¹ ⋅ b
.