GMRES sparse linear solver

The generalized minimum residual method (GMRES) is a popular iterative method for the numerical solution of a sparse nonsymmetric system of linear equations.

ALGLIB, a free and commercial open source numerical library, includes a high-quality GMRES implementation as well as many other large-scale sparse linear solvers available in C++, C#, Java and several other programming languages. The solver suite is dual-licensed, with free and commercial editions.

Programming languages supported

ALGLIB supports many programming languages, including C++, C#, Java, Python, and others:

• C++ version. Highly optimized, portable, self-contained library (x86/x64/ARM; works on Windows, Linux, and POSIX systems)
• C# version. Unified C# API for two backends: a fully managed C# implementation and a high-performance native computational core written in C
• Java, Python, and other languages. A wrapper for a high-performance native computational core written in C. Seamless integration with the target programming language

A distinctive feature of ALGLIB is that it provides the same API in all programming languages. This is achieved through our exclusive automatic code translation and wrapper generation technology.

GMRES algorithm

GMRES(k) is an iterative algorithm developed by Yousef Saad and Martin H. Schultz in 1986. The core idea of the algorithm is to approximate the exact solution of A·x=b by the vector x_k ∈ K_k(A,b) such that x_k minimizes the residual r = A·x-b. Here, K_k(A,b) is a k-th Krylov subspace, that is, an orthogonal basis spanning (r, A·r, A²·r, ..., A^k-1·r). One usually chooses k much less than the variables count n and iteratively repeats the solution process until the residual becomes small enough.

Although it is not limited by sparse linear equations, the GMRES(k) algorithm is most often used to solve sparse linear systems. It has modest and predictable memory requirements: O(k·N) in addition to the matrix itself. Additional running time overhead associated with each iteration is O(k²·N) (linear in N), which makes this method suitable for large-scale problems.

GMRES solver API

Convenience wrappers

The easiest way to use the GMRES solver is to call either sparsesolvegmres (general nonsymmetric matrix) or its symmetric version sparsesolvesymmetricgmres. A single function call captures all the necessary details: the system, its right part, restart frequency k, desired tolerance ε and iterations limit M.

Out-of-core mode

An out-of-core GMRES mode is intended for situations when the system matrix A is unavailable, but we can compute its product A·x for any x.

Some engineering problems allow one to perform a trial at any desired point, i.e. to compute A·x, but do not provide an easy way to generate A itself. Another situation when the out-of-core mode is necessary is when one solves a preconditioned problem A·P·x=b with A being available but the preconditioner P is given only by its product with x. Thus, although A is explicitly known, the entire system A·P is given only implicitly.

The out-of-core API is built around the sparsesolverstate structure. First, the structure is initialized with the sparsesolvercreate function and the algorithm is set to GMRES with the sparsesolversetalgogmres call. After that, a user may optionally set the stopping criteria with the sparsesolversetcond function and provide an initial point with the sparsesolversetstartingpoint function. Finally, the out-of-core mode of the solver is activated as demonstrated in the pseudocode below:

alglib.sparsesolveroocstart(state) while alglib.sparsesolverooccontinue(state) do alglib.sparsesolveroocgetrequestinfo(state, out RequestType) alglib.sparsesolveroocgetrequestdata(state, out X) if RequestType=0 then [calculate y=A*x] alglib.sparsesolveroocsendresult(state, in Y) alglib.sparsesolveroocstop(state, out X, out Report)

The pseudocode above can be easily adapted to any programming language supported by ALGLIB.

This article is licensed for personal use only.