ALGLIB contains recursive matrix inversion algorithm. A, given by its triangular factorization (LU-decomposition or Cholesky decomposition) is partitioned into four parts and inverted according to the next formula:
Matrix-matrix products of general form and multiplications by U -1/L -1 are calculated with cache-efficient Level 3 ALGLIB BLAS subroutines. Expressions like (LU) -1 are calculated recursively. Base cases (matrix size is small enough to fit into L1 cache) are handled by another algorithm, which is based on explicit inversion of triangular factors of A. Recursive algorithm has following benefits:
These features allows to achieve the maximum performance possible with your compiler/language. However, different compilers and different programming languages have different limitations. As for linear algebra algorithms, C++ version of ALGLIB is the most efficient implementation. Implementations in other languages have significantly lower performance.
ALGLIB package contains following matrix inversion subroutines:
First four subroutines are used to invert real (R prefix), symmetric positive definite (SPD prefix), complex (C prefix) or Hermitian positive definite (HPD prefix) matrices. Next four subroutines solve same task, but for matrices given by corresponding triangular factorization (LU decomposition or Cholesky decomposition).
All subroutines check condition number before inversion begins. If matrix is degenerate (one of triangular factors have zero element on its diagonal) or is very ill-conditioned (condition number is so high that there is a risk of overflow during inversion) then subroutine returns -3 error code and zero matrix.
This article is intended for personal use only.
C++ source. MPFR/GMP is used.
GMP source is available from gmplib.org. MPFR source is available from www.mpfr.org.
Python version (CPython and IronPython are supported).
ALGLIB® - numerical analysis library, 1999-2013.