Function description

List of solvers in src.solvers

src.solvers.pd_gmres(A, b, mInitial, mMinMax, mStep, tol, maxit, xInitial, alphaPD, varargin)

PD-GMRES Proportional-Derivative GMRES(m)

Description:

pd_gmres is a modified implementation of the restarted Generalized Minimal Residual Error or GMRES(m) [1], performed by using a proportional-derivative control-inspired law to update adaptively the restarting parameter m before each restart. This implementation follows closely the one presented in [2].

Signature:

[x, flag, relresvec, kdvec, time] = pd_gmres(A, b, …

mInitial, mMinMax, mStep, tol, maxit, xInitial, alphaPD)

Input Parameters:

A: n-by-n matrix

Left-hand side of the linear system Ax = b.

b: n-by-1 vector

Right-hand side of the linear system Ax = b.

mInitial: int, optional

Initial restart parameter (similar to ‘restart’ in MATLAB). If empty, not given, or equal to n, then mInitial is not used and the full unrestarted gmres algorithm with maxit = min(n, 10) is employed. Note that we require 1 <= mInitial <= n.

mMinMax: 2-by-1 vector, optional

Minimum and maximum values of the restart paramter m. Default is [1; n]. We require 1 <= mMinMax(1) < mMinMax(2) <= n. Note that for large matrices (e.g., > 600), the default value mMinMax(2) = n might require large computational resources.

mStep: int, optional

Step size for increasing the restart parameter m when m < mMinMax(1). Default is 1 if n <= 10 and 3 otherwise.

tol: float, optional

Tolerance error threshold for the relative residual norm. Default is 1e-6.

maxit: int, optional

Maximum number of outer iterations.

xInitial: n-by-1 vector, optional

Vector of initial guess. Default is zeros(n, 1).

alphaPD: 2-by-1 vector, optional

Proportional and derivative coefficients from the proportional-derivative controller. Default is [-3, 5], following Ref. [2].

Output parameters:

x: n-by-1 vector

Approximate solution of the linear system.

flag: boolean

1 if the algorithm has converged, 0 otherwise.

relressvec: (1 up to maxit)-by-1 vector

Vector of relative residual norms of every outer iteration (cycles). The last relative residual norm is simply given by relresvec(end).

kdvec: (1 up to maxit)-by-1 vector

Vector of restart parameter values. In case the unrestarted algorithm is invoked, kdvec = NaN.

time: scalar

Computational time in seconds.

References:

[1] Saad, Y., & Schultz, M. H. (1986). GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM Journal on scientific and statistical computing, 7(3), 856-869.

[2] Nunez, R. C., Schaerer, C. E., & Bhaya, A. (2018). A proportional-derivative control strategy for restarting the GMRES(m) algorithm. Journal of Computational and Applied Mathematics, 337, 209-224.

src.solvers.lgmres(A, b, m, l, tol, maxit, xInitial, varargin)

LGMRES algorithm

Description:

LGMRES (“Loose GMRES”) is a modified implementation of the restarted Generalized Minimal Residual Error or GMRES(m) [1], performed by appending ‘l’ error approximation vectors to the restarting Krylov subspace, as a way to preserve information from previous discarted search subspaces from previous iterations of the method.

Augments the standard GMRES approximation space with approximations to the error from previous restart cycles as in [1].

Signature:

[x, flag, relresvec, time] = lgmres(A, b, m, l, tol, maxit, xInitial)

Input parameters:

A: n-by-n matrix

Left-hand side of the linear system Ax = b.

b: n-by-1 vector

Right-hand side of the linear system Ax = b.

m: int, optional

Restart parameter (similar to ‘restart’ in MATLAB). Default is min(n, 10). If m == n, built-in unrestarted gmres will be used.

l: int, optional

Number of error approximation vectors to be appended to the Krylov search subspace. Default is 3, but values between 1 and 5 are mostly used. If m < n AND l == 0, built-in gmres(m) will be used.

tol: float, optional

Tolerance error threshold for the relative residual norm. Default is 1e-6.

maxit: int, optional

Maximum number of outer iterations. Default is min(m, 10).

xInitial: n-by-1 vector, optional

Vector of initial guess. Default is zeros(n, 1).

Output parameters:

x: n-by-1 vector

Approximate solution to the linear system.

flag: boolean

1 if the algorithm has converged, 0 otherwise.

relresvec: (1 up to maxit)-by-1 vector

Vector of relative residual norms of every outer iteration (cycles). The last relative residual norm is simply given by relresvec(end).

kdvec: (1 up to maxit)-by-1 vector

For LGMRES, kdvec is a constant vector whose elements correspond to the size of the Krylov subspace, i.e., m + l. Note that in some cases, there could be a “happy breakdown” where the dimension of the search space < m + l.

time: scalar

Computational time in seconds.

References:

[1] Baker, A. H., Jessup, E. R., & Manteuffel, T. (2005). A technique for accelerating the convergence of restarted GMRES. SIAM Journal on Matrix Analysis and Applications, 26(4), 962-984.

Copyright:

This file is part of the KrySBAS MATLAB Toolbox.

Copyright 2023 CC&MA - NIDTec - FP - UNA

KrySBAS is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

KrySBAS is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this file. If not, see <http://www.gnu.org/licenses/>.

src.solvers.gmres_e(A, b, m, d, tol, maxit, xInitial, eigstol, varargin)

GMRES-E algorithm

Description:

GMRES-E is a modified implementation of the restarted Generalized Minimal Residual Error or GMRES(m) [1], performed by appending ‘d’ eigenvectors corresponding to a few of the smallest eigenvalues in magnitude for each outer iteration. In practice, the approximate eigenvectors are the harmonic Ritz vectors associated to the harmonic Ritz values per outer iteration.

Signature:

[x, flag, relresvec, time] = …

gmres_e(A, b, m, k, tol, maxit, xInitial, eigstol)

Input Parameters:

A: n-by-n matrix

Left-hand side of the linear system Ax = b.

b: n-by-1 vector

Right-hand side of the linear system Ax = b.

m: int

Restart parameter (similar to ‘restart’ in MATLAB). If m == n, built-in unrestarted gmres will be used

d: int

Number of eigenvectors corresponding to a few of the smallest eigenvalues in magnitude for each outer iteration. Default is min(m, 3), but values between 1 and 5 are typical. According to [1], “even just a few eigenvectors can make a big difference if the matrix has both small and large eigenvalues”. If m < n AND d == 0, the built-in gmres(m) will be used. If d > m, an error is raised.

tol: float, optional

Tolerance error threshold for the relative residual norm. Default is 1e-6.

maxit: int, optional

Maximum number of outer iterations. Default is min(m, 10).

xInitial: n-by-1 vector, optional

Vector of initial guess. Default is zeros(n, 1).

eigstol: float, optional

Tolerance for computing eigenvectors using the built-in MATLAB function eigs. Default is 1e-6.

Output parameters:

x: n-by-1 vector

Approximate solution to the linear system.

flag: boolean

1 if the algorithm has converged, 0 otherwise.

relresvec: (1 up to maxit)-by-1 vector

Vector of relative residual norms of every outer iteration (cycles). The last relative residual norm is simply given by relresvec(end).

kdvec: (1 up to maxit)-by-1 vector

For LGMRES, kdvec is a constant vector whose elements correspond to the size of the Krylov subspace, i.e., m + d. Note that in some cases, there could be a “happy breakdown” where the dimension of the search space < m + d.

time: float

Computational time in seconds.

References:

[1] Morgan, R. B. (1995). A restarted GMRES method augmented with eigenvectors. SIAM Journal on Matrix Analysis and Applications, 16(4), 1154-1171.

Copyright:

This file is part of the KrySBAS MATLAB Toolbox.

Copyright 2023 CC&MA - NIDTec - FP - UNA

KrySBAS is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

KrySBAS is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this file. If not, see <http://www.gnu.org/licenses/>.

List of utility functions in src.utils

src.utils.augmented_gram_schmidt_arnoldi(A, v, m, appendV)

Modified Gram-Schmidt Arnoldi iteration, with Krylov subspace augmentation

Description:

Implementation of the augmented Gram-Schmidt Arnoldi iteration following [1], Algorithm 2.1.

Syntaxis:

[H, V, s] = augmented_gram_schmidt_arnoldi(A, v, m, appendV)

Input parameters:

A: n-by-n matrix

Matrix of coefficients.

v: n-by-1 vector

Normalized residual vector.

m: int

Restart parameter.

appendV: n-by-k matrix

Matrix of information vectors to append to the current Krylov subspace

Output parameters:

H: m+1-by-m matrix

Upper Hessenberg matrix.

V: n-by-m matrix

Orthonormal basis of Krylov subspace.

s: int

Size of the new search subspace

References:

[1] Saad, Y. (1997). Analysis of augmented Krylov subspace methods. SIAM Journal on Matrix Analysis and Applications, 18(2), 435-449.

Copyright:

This file is part of the KrySBAS MATLAB Toolbox.

Copyright 2023 CC&MA - NIDTec - FP - UNA

KrySBAS is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

KrySBAS is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this file. If not, see <http://www.gnu.org/licenses/>.

src.utils.harmonic_ritz_vectors(F, G, k, V, tol)

Harmonic Ritz Vectors function

Description:

This function is a modified implementation of the Rayleigh-Ritz method that finds good approximations to the smallest eigenvalues and its associated eigenvectors, and appends these ones to the next search subspaces. Refer to [1] for more information about this technique.

Signature:

dy = harmonic_ritz_vectors(F, G, k, V, tol)

Input Parameters:

F: s-by-s maxtrix

Matrix F from the generalized eigenvalue problem, as employed in step 5, p. 1161 of [1].

G: s-by-1 vector

Matrix G from the generalized eigenvalue problem, as employed in step 5, p. 1161 of [1].

k: int

Number of eigenvectors corresponding to a few of the smallest eigenvalues in magnitude. According to [1], “even just a few eigenvectors can make a big difference if the matrix has both small and large eigenvalues”.

V: n-by-s matrix

A matrix from the modified Gram-Schmidt Arnoldi method, whose column vectors span the search subspace.

tol: float

Convergence tolerance for the built-in ‘eigs’ algorithm.

Output parameters:

dy: n-by-k matrix

Matrix whose column vectors are the approximate eigenvectors of the matrix A.

References:

[1] Morgan, R. B. (1995). A restarted GMRES method augmented with eigenvectors. SIAM Journal on Matrix Analysis and Applications, 16(4), 1154-1171.

Copyright:

This file is part of the KrySBAS MATLAB Toolbox.

Copyright 2023 CC&MA - NIDTec - FP - UNA

KrySBAS is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

KrySBAS is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this file. If not, see <http://www.gnu.org/licenses/>.

src.utils.modified_gram_schmidt_arnoldi(A, v, m)

Modified Gram-Schmidt Arnoldi iteration

Description:

Implementation of the modified Gram-Schmidt Arnoldi iteration following [1], Algorithm 6.2.

Syntaxis:

[H, V, mUpdated] = modified_gram_schmidt_arnoldi(A, v, m)

Input parameters:

A: n-by-n matrix

Matrix of coefficients.

v: n-by-1 vector

Normalized residual vector.

m: int

Restart parameter.

Output parameters:

H: m+1-by-m matrix

Upper Hessenberg matrix.

V: n-by-m matrix

Orthonormal basis of Krylov subspace.

mUpdated: int

Updated value of the restart parameter ‘m’. This parameter changes only if the last element of H is exactly 0 and the matrix is trunctated. See, in particular, the if-else statement inside the Modified Gram-Schmidt block.

References:

[1] Saad, Y. (2003). Iterative methods for sparse linear systems. Society for Industrial and Applied Mathematics.

Copyright:

This file is part of the KrySBAS MATLAB Toolbox.

Copyright 2023 CC&MA - NIDTec - FP - UNA

KrySBAS is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

KrySBAS is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this file. If not, see <http://www.gnu.org/licenses/>.

src.utils.pd_rule(m, n, mInitial, mMin, mMax, mStep, res, iter, alphaP, alphaD)

Rule for the Proportional-Derivative law. Algorithm 1 of [1].

Signature:

miter = (m, n, mInitial, mMin, mMax, mStep, …

res, iter, alphaP, alphaD)

Input Parameters:

m: int

Restart parameter at the last cycle

n: int

size of the square matrix A

mInitial: int

Initial restart parameter at the last cycle.

mMin: int

Minimum value of the restart parameter m.

res number of cycles-by-1 vector

Vector of residual norms of every outer iteration (cycles).

iter number of cycles-by-1 vector

Number of restart cycles previous.

mStep: int

Step size for increasing the mInitial when m < mMinMax(1).

mMax: int

Maximum value of the restart paramter m.

maxit: int

Maximum number of outer iterations.

alphaP: int

Proportional coefficient from PD controller.

alphaD: int

Derivative coefficient from PD controller.

Output parameters:

miter: 1-by-2 vector

Restart parameter and Initial restart parameter at the new cycle.

% References:

[1] Nunez, R. C., Schaerer, C. E., & Bhaya, A. (2018). A proportional-derivative control strategy for restarting the GMRES(m) algorithm. Journal of Computational and Applied Mathematics, 337, 209-224.

Copyright:

This file is part of the KrySBAS MATLAB Toolbox.

Copyright 2023 CC&MA - NIDTec - FP - UNA

KrySBAS is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

KrySBAS is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this file. If not, see <http://www.gnu.org/licenses/>.

src.utils.plane_rotations(H, beta)

Performs plane rotations.

Description:

Implementation of plane rotations (a.k.a Givens rotation) following [1].

Syntaxis:

[HUpTri, g] = plane_rotations(H, beta)

Input parameters:

H: m+1-by-m matrix

Upper Hessenberg matrix.

beta: int

Norm of the last cycle residual vector.

Output parameters:

HUpTri: m+1-by-m matrix

Upper-triangular Hessenberg matrix.

g: m+1-by-1 vector

Resulting right-hand-side vector.

Notes:

HUpTri and g (with the last rows deleted) will enter into the least-squares problem.

References:

[1] Saad, Y. (2003). Iterative methods for sparse linear systems. Society for Industrial and Applied Mathematics.

Copyright:

This file is part of the KrySBAS MATLAB Toolbox.

Copyright 2023 CC&MA - NIDTec - FP - UNA

KrySBAS is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

KrySBAS is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this file. If not, see <http://www.gnu.org/licenses/>.