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.

Contents

    1 Subroutines: general system
    2 Subroutines: SPD/HPD system
    3 Subroutines: non-square possibly singular system
    4 Performance
    5 Iterative refinement
    6 Downloads section

Subroutines: general system

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

These subroutines can be divided into following groups:

The latter requires some more explanations:

Subroutines: SPD/HPD system

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

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:

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 solutionx0 has been found, we calculate residual r0=b-Ax0 and use it to refine solution: Ad0=r0, x1=x0+d0 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-Ax0 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-Ax0 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.

Download ALGLIB for C++ / C# / Java / Python / ...

ALGLIB Project offers you two editions of ALGLIB:

ALGLIB Free Edition:
+delivered for free
+offers full set of numerical functionality
+extensive algorithmic optimizations
-no multithreading
-non-commercial license

ALGLIB Commercial Edition:
+flexible pricing
+offers full set of numerical functionality
+extensive algorithmic optimizations
+high performance (SMP, SIMD)
+commercial license with support plan

Links to download sections for Free and Commercial editions can be found below:

ALGLIB 4.01.0 for C++

C++ library.
Delivered with sources.
Monolithic design.
Extreme portability.
Editions:   FREE   COMMERCIAL

ALGLIB 4.01.0 for C#

C# library with native kernels.
Delivered with sources.
VB.NET and IronPython wrappers.
Extreme portability.
Editions:   FREE   COMMERCIAL

ALGLIB 4.01.0 for Java

Java wrapper around HPC core.
Delivered with sources.
Seamless integration with Java.
Editions:   FREE   COMMERCIAL

ALGLIB 4.01.0 for Delphi

Delphi wrapper around C core.
Delivered as precompiled binary.
Compatible with FreePascal.
Editions:   FREE   COMMERCIAL

ALGLIB 4.01.0 for CPython

CPython wrapper around C core.
Delivered as precompiled binary.
Editions:   FREE   COMMERCIAL