Dense linear solvers

Systems of linear equations Ax=b may be divided into two classes: those with square non-degenerate A and those with rectangular possibly rank deficient A. ALGLIB package have subroutines for both types of problems.

Subroutines: general system

The following subroutines can be used to solve real/complex systems with general form NxN non-singular matrix:

• rmatrixsolve
• rmatrixsolvem
• rmatrixlusolve
• rmatrixlusolvem
• rmatrixmixedsolve
• rmatrixmixedsolvem
• cmatrixsolve
• cmatrixsolvem
• cmatrixlusolve
• cmatrixlusolvem
• cmatrixmixedsolve
• cmatrixmixedsolvem

These subroutines can be divided into following groups:

• Those which accept real matrix (name starts with R) and those which accept complex matrix (name starts with C).
• Those which solve task with one right-hand part and those which solve task with multiple right-hand parts (name ends with M).
• Those which accept only A, those which acceps its LU decomposition (name contains LU) and those which accept both A and its LU decomposition (name contains Mixed).

The latter requires some more explanations:

• Subroutines which accept A build its LU decomposition themselves in O(N³) time. Solution found may be iteratively refined (if RFC parameter is True) using extra-precision arithmetic (see below).
• Subroutines which accept LU decomposition of A may solve system in O(N²) time. However they can't apply iterative refinement because it requires original matrix A, not its inprecise factorization.
• Subroutines which accept both A and its LU decomposition combine advantages of both variants described above: high precision and high speed.

Subroutines: SPD/HPD system

The following subroutines can be used to solve real/complex systems with symmetric/Hermitian positive definite matrix:

• spdmatrixsolve
• spdmatrixsolvem
• spdmatrixcholeskysolve
• spdmatrixcholeskysolvem
• hpdmatrixsolve
• hpdmatrixsolvem
• hpdmatrixcholeskysolve
• hpdmatrixcholeskysolvem

Naming conventions are the same as in previous section. It should be noted that subroutines from this section do not support iterative refinement.

Subroutines: non-square possibly singular system

rmatrixsolvels solves systems with rectangular possibly degenerate matrix. Is system matrix is rank deficient, solution in a least squares sense is searched for. Subroutine supports iterative refinement. Current version of ALGLIB doesn't include complex version of this subroutine or version which accept multiple right-hand parts.

Performance

Several factors influence linear solver performance:

• Performance of triangular factorization. Highly optimized cache-oblivious algorithms are used to calculate triangular factorization. However, performance of these algorithms depends on what version of ALGLIB you are using. C++ version is the most efficient, pure NET version and versions in other languages have significantly lower performance.
• Number of right-hand parts. ALGLIB solver is optimized for systems with one right-hand part or for systems with low number of right-hand parts. If right-hand part is too large then you'll experience a performance drop becauce there is no blocking for different levels of cache for the right-hand part. In such cases you can factorize matrix yourself using cache-efficient subroutines from trfac unit, and then you can use TRSM-subroutines from ablas unit for cache-efficient solution of your system.
• Iterative improvement. Its cost is negligible in comparison with the cost of the factorization - when you solve system with one right part. However, it can require a significant part of computational resources when matrix size is very small or when you have large number of right-hand sides. In such cases you can turn it off.
• Matrix copying during solution. Algorithms which accept LU decomposition don't allocate significant amounts of additional memory. But algorithms which accept original matrix require additional memory to create its copy. However, allocation cost is negligible when matrix size is significantly less than total amount of physical RAM.

Iterative refinement

Iterative refinement is a method proposed by James H. Wilkinson to improve the accuracy of numerical solutions to systems of linear equations. After system has been solved and solutionx₀ has been found, we calculate residual r₀=b-Ax₀ and use it to refine solution: Ad₀=r₀, x₁=x₀+d₀ Operation may be repeated several times. In ALGLIB package two-pass iterative refinement is used.

Note #1
Iterative refinement can decrease errors which arise due to floating point round-off. It can make solution as precise as condition number of A allows. However, refinement result depends on accuracy with which b-Ax₀ is calculated. For iterative refinement to be worth the effort this difference must be calculated with much higher precision than that was used to find a solution. ALGLIB solver uses portable extra-precise matrix multiplier. It allows to calculate b-Ax₀ with maximum precision possible given than A and b may have errors as large as one ULP. "Portable" in this context means that it doesn't use architecture-specific or compiler-specific extensions (like 80-bit Intel floating point type) nor it contains low level assembler code.

This article is licensed for personal use only.