18.5 Specialized Solvers

x = bicg (A, b)

x = bicg (A, b, tol)

x = bicg (A, b, tol, maxit)

x = bicg (A, b, tol, maxit, M)

x = bicg (A, b, tol, maxit, M1, M2)

x = bicg (A, b, tol, maxit, M, [], x0)

x = bicg (A, b, tol, maxit, M1, M2, x0)

x = bicg (A, b, tol, maxit, M, [], x0, …)

x = bicg (A, b, tol, maxit, M1, M2, x0, …)

[x, flag, relres, iter, resvec] = bicg (A, b, …)

Solve the linear system of equations A * x = b by means of the Bi-Conjugate Gradient iterative method.

The input arguments are:

• A is the matrix of the linear system and it must be square. A can be passed as a matrix, function handle, or inline function Afun such that Afun (x, "notransp") = A * x and Afun (x, "transp") = A' * x. Additional parameters to Afun may be passed after x0.
• b is the right-hand side vector. It must be a column vector with the same number of rows as A.
• tol is the required relative tolerance for the residual error, b - A * x. The iteration stops if norm (b - A * x) ≤ tol * norm (b). If tol is omitted or empty, then a tolerance of 1e-6 is used.
• maxit is the maximum allowed number of iterations; if maxit is omitted or empty then a value of 20 is used.
• M1, M2 are the preconditioners. The preconditioner M is given as M = M1 * M2. Both M1 and M2 can be passed as a matrix or as a function handle or inline function g such that g (x, "notransp") = M1 \ x or g (x, "notransp") = M2 \ x and g (x, "transp") = M1' \ x or g (x, "transp") = M2' \ x. If M1 is omitted or empty, then preconditioning is not applied. The preconditioned system is theoretically equivalent to applying the bicg method to the linear system inv (M1) * A * inv (M2) * y = inv (M1) * b and inv (M2') * A' * inv (M1') * z = inv (M2') * b and then setting x = inv (M2) * y.
• x0 is the initial guess. If x0 is omitted or empty then the function sets x0 to a zero vector by default.

Any arguments which follow x0 are treated as parameters, and passed in an appropriate manner to any of the functions (Afun or Mfun) or that have been given to bicg.

The output parameters are:

• x is the computed approximation to the solution of A * x = b. If the algorithm did not converge, then x is the iteration which has the minimum residual.
• flag indicates the exit status:
  • 0: The algorithm converged to within the prescribed tolerance.
  • 1: The algorithm did not converge and it reached the maximum number of iterations.
  • 2: The preconditioner matrix is singular.
  • 3: The algorithm stagnated, i.e., the absolute value of the difference between the current iteration x and the previous is less than eps * norm (x,2).
  • 4: The algorithm could not continue because intermediate values became too small or too large for reliable computation.
• relres is the ratio of the final residual to its initial value, measured in the Euclidean norm.
• iter is the iteration which x is computed.
• resvec is a vector containing the residual at each iteration. The total number of iterations performed is given by length (resvec) - 1.

Consider a trivial problem with a tridiagonal matrix

n = 20;
A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ...
    toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
              sparse (1, 2, 1, 1, n) * n / 2);
b = A * ones (n, 1);
restart = 5;
[M1, M2] = ilu (A);  # in this tridiag case, it corresponds to lu (A)
M = M1 * M2;
Afun = @(x, string) strcmp (string, "notransp") * (A * x) + ...
                     strcmp (string, "transp") * (A' * x);
Mfun = @(x, string) strcmp (string, "notransp") * (M \ x) + ...
                     strcmp (string, "transp") * (M' \ x);
M1fun = @(x, string) strcmp (string, "notransp") * (M1 \ x) + ...
                     strcmp (string, "transp") * (M1' \ x);
M2fun = @(x, string) strcmp (string, "notransp") * (M2 \ x) + ...
                     strcmp (string, "transp") * (M2' \ x);

EXAMPLE 1: simplest usage of bicg

x = bicg (A, b)

EXAMPLE 2: bicg with a function that computes A*x and A'*x

x = bicg (Afun, b, [], n)

EXAMPLE 3: bicg with a preconditioner matrix M

x = bicg (A, b, 1e-6, n, M)

EXAMPLE 4: bicg with a function as preconditioner

x = bicg (Afun, b, 1e-6, n, Mfun)

EXAMPLE 5: bicg with preconditioner matrices M1 and M2

x = bicg (A, b, 1e-6, n, M1, M2)

EXAMPLE 6: bicg with functions as preconditioners

x = bicg (Afun, b, 1e-6, n, M1fun, M2fun)

EXAMPLE 7: bicg with as input a function requiring an argument

function y = Ap (A, x, string, z)
  ## compute A^z * x or (A^z)' * x
  y = x;
  if (strcmp (string, "notransp"))
    for i = 1:z
      y = A * y;
    endfor
  elseif (strcmp (string, "transp"))
    for i = 1:z
      y = A' * y;
    endfor
  endif
endfunction

Apfun = @(x, string, p) Ap (A, x, string, p);
x = bicg (Apfun, b, [], [], [], [], [], 2);

References:

1. Y. Saad, Iterative Methods for Sparse Linear Systems, Second edition, 2003, SIAM.

See also: bicgstab, cgs, gmres, pcg, qmr, tfqmr.

x = bicgstab (A, b, tol, maxit, M1, M2, x0, …)

x = bicgstab (A, b, tol, maxit, M, [], x0, …)

[x, flag, relres, iter, resvec] = bicgstab (A, b, …)

Solve A x = b using the stabilized Bi-conjugate gradient iterative method.

The input parameters are:

• - A is the matrix of the linear system and it must be square. A can be passed as a matrix, function handle, or inline function Afun such that Afun(x) = A * x. Additional parameters to Afun are passed after x0.
• - b is the right hand side vector. It must be a column vector with the same number of rows as A.
• - tol is the required relative tolerance for the residual error, b - A * x. The iteration stops if norm (b - A * x) ≤ tol * norm (b). If tol is omitted or empty, then a tolerance of 1e-6 is used.
• - maxit the maximum number of outer iterations, if not given or set to [] the default value min (20, numel (b)) is used.
• - M1, M2 are the preconditioners. The preconditioner M is given as M = M1 * M2. Both M1 and M2 can be passed as a matrix or as a function handle or inline function g such that g(x) = M1 \ x or g(x) = M2 \ x. The technique used is the right preconditioning, i.e., it is solved A * inv (M) * y = b and then x = inv (M) * y.
• - x0 the initial guess, if not given or set to [] the default value zeros (size (b)) is used.

The arguments which follow x0 are treated as parameters, and passed in a proper way to any of the functions (A or M) which are passed to bicstab.

The output parameters are:

• - x is the approximation computed. If the method doesn’t converge then it is the iterated with the minimum residual.
• - flag indicates the exit status:
  • - 0: iteration converged to the within the chosen tolerance
  • - 1: the maximum number of iterations was reached before convergence
  • - 2: the preconditioner matrix is singular
  • - 3: the algorithm reached stagnation
  • - 4: the algorithm can’t continue due to a division by zero
• - relres is the relative residual obtained with as (A*x-b) / norm(b).
• - iter is the (possibly half) iteration which x is computed. If it is an half iteration then it is iter + 0.5
• - resvec is a vector containing the residual of each half and total iteration (There are also the half iterations since x is computed in two steps at each iteration). Doing (length(resvec) - 1) / 2 is possible to see the total number of (total) iterations performed.

Let us consider a trivial problem with a tridiagonal matrix

n = 20;
A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n))  + ...
    toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
    sparse (1, 2, 1, 1, n) * n / 2);
b = A * ones (n, 1);
restart = 5;
[M1, M2] = ilu (A); # in this tridiag case, it corresponds to lu (A)
M = M1 * M2;
Afun = @(x) A * x;
Mfun = @(x) M \ x;
M1fun = @(x) M1 \ x;
M2fun = @(x) M2 \ x;

EXAMPLE 1: simplest usage of bicgstab

x = bicgstab (A, b, [], n)

EXAMPLE 2: bicgstab with a function which computes A * x

x = bicgstab (Afun, b, [], n)

EXAMPLE 3: bicgstab with a preconditioner matrix M

x = bicgstab (A, b, [], 1e-06, n, M)

EXAMPLE 4: bicgstab with a function as preconditioner

x = bicgstab (Afun, b, 1e-6, n, Mfun)

EXAMPLE 5: bicgstab with preconditioner matrices M1 and M2

x = bicgstab (A, b, [], 1e-6, n, M1, M2)

EXAMPLE 6: bicgstab with functions as preconditioners

x = bicgstab (Afun, b, 1e-6, n, M1fun, M2fun)

EXAMPLE 7: bicgstab with as input a function requiring an argument

function y = Ap (A, x, z) # compute A^z * x
   y = x;
   for i = 1:z
     y = A * y;
   endfor
 endfunction
Apfun = @(x, string, p) Ap (A, x, string, p);
x = bicgstab (Apfun, b, [], [], [], [], [], 2);

EXAMPLE 8: explicit example to show that bicgstab uses a right preconditioner

[M1, M2] = ilu (A + 0.1 * eye (n)); # factorization of A perturbed
M = M1 * M2;

## reference solution computed by bicgstab after one iteration
[x_ref, fl] = bicgstab (A, b, [], 1, M)

## right preconditioning
[y, fl] = bicgstab (A / M, b, [], 1)
x = M \ y # compare x and x_ref

References:

1. Y. Saad, Iterative Methods for Sparse Linear Systems, Second edition, 2003, SIAM

See also: bicg, cgs, gmres, pcg, qmr, tfqmr.

x = cgs (A, b, tol, maxit, M1, M2, x0, …)

x = cgs (A, b, tol, maxit, M, [], x0, …)

[x, flag, relres, iter, resvec] = cgs (A, b, …)

Solve A x = b, where A is a square matrix, using the Conjugate Gradients Squared method.

The input arguments are:

• - A is the matrix of the linear system and it must be square. A can be passed as a matrix, function handle, or inline function Afun such that Afun(x) = A * x. Additional parameters to Afun are passed after x0.
• - b is the right hand side vector. It must be a column vector with same number of rows of A.
• - tol is the relative tolerance, if not given or set to [] the default value 1e-6 is used.
• - maxit the maximum number of outer iterations, if not given or set to [] the default value min (20, numel (b)) is used.
• - M1, M2 are the preconditioners. The preconditioner matrix is given as M = M1 * M2. Both M1 and M2 can be passed as a matrix or as a function handle or inline function g such that g(x) = M1 \ x or g(x) = M2 \ x. If M1 is empty or not passed then no preconditioners are applied. The technique used is the right preconditioning, i.e., it is solved A*inv(M)*y = b and then x = inv(M)*y.
• - x0 the initial guess, if not given or set to [] the default value zeros (size (b)) is used.

The arguments which follow x0 are treated as parameters, and passed in a proper way to any of the functions (A or P) which are passed to cgs.

The output parameters are:

• - x is the approximation computed. If the method doesn’t converge then it is the iterated with the minimum residual.
• - flag indicates the exit status:
  • - 0: iteration converged to the within the chosen tolerance
  • - 1: the maximum number of iterations was reached before convergence
  • - 2: the preconditioner matrix is singular
  • - 3: the algorithm reached stagnation
  • - 4: the algorithm can’t continue due to a division by zero
• - relres is the relative residual obtained with as (A*x-b) / norm(b).
• - iter is the iteration which x is computed.
• - resvec is a vector containing the residual at each iteration. Doing length(resvec) - 1 is possible to see the total number of iterations performed.

Let us consider a trivial problem with a tridiagonal matrix

n = 20;
A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n))  + ...
    toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
    sparse (1, 2, 1, 1, n) * n / 2);
b = A * ones (n, 1);
restart = 5;
[M1, M2] = ilu (A); # in this tridiag case it corresponds to chol (A)'
M = M1 * M2;
Afun = @(x) A * x;
Mfun = @(x) M \ x;
M1fun = @(x) M1 \ x;
M2fun = @(x) M2 \ x;

EXAMPLE 1: simplest usage of cgs

x = cgs (A, b, [], n)

EXAMPLE 2: cgs with a function which computes A * x

x = cgs (Afun, b, [], n)

EXAMPLE 3: cgs with a preconditioner matrix M

x = cgs (A, b, [], 1e-06, n, M)

EXAMPLE 4: cgs with a function as preconditioner

x = cgs (Afun, b, 1e-6, n, Mfun)

EXAMPLE 5: cgs with preconditioner matrices M1 and M2

x = cgs (A, b, [], 1e-6, n, M1, M2)

EXAMPLE 6: cgs with functions as preconditioners

x = cgs (Afun, b, 1e-6, n, M1fun, M2fun)

EXAMPLE 7: cgs with as input a function requiring an argument

function y = Ap (A, x, z) # compute A^z * x
   y = x;
   for i = 1:z
     y = A * y;
   endfor
 endfunction
Apfun = @(x, string, p) Ap (A, x, string, p);
x = cgs (Apfun, b, [], [], [], [], [], 2);

EXAMPLE 8: explicit example to show that cgs uses a right preconditioner

[M1, M2] = ilu (A + 0.3 * eye (n)); # factorization of A perturbed
M = M1 * M2;

## reference solution computed by cgs after one iteration
[x_ref, fl] = cgs (A, b, [], 1, M)

## right preconditioning
[y, fl] = cgs (A / M, b, [], 1)
x = M \ y # compare x and x_ref

References:

1. Y. Saad, Iterative Methods for Sparse Linear Systems, Second edition, 2003, SIAM

See also: pcg, bicgstab, bicg, gmres, qmr, tfqmr.

x = gmres (A, b, restart, tol, maxit, M1, M2, x0, …)

x = gmres (A, b, restart, tol, maxit, M, [], x0, …)

[x, flag, relres, iter, resvec] = gmres (A, b, …)

Solve A x = b using the Preconditioned GMRES iterative method with restart, a.k.a. PGMRES(restart).

The input arguments are:

• - A is the matrix of the linear system and it must be square. A can be passed as a matrix, function handle, or inline function Afun such that Afun(x) = A * x. Additional parameters to Afun are passed after x0.
• - b is the right hand side vector. It must be a column vector with the same numbers of rows as A.
• - restart is the number of iterations before that the method restarts. If it is [] or N = numel (b), then the restart is not applied.
• - tol is the required relative tolerance for the preconditioned residual error, inv (M) * (b - a * x). The iteration stops if norm (inv (M) * (b - a * x)) ≤ tol * norm (inv (M) * B). If tol is omitted or empty, then a tolerance of 1e-6 is used.
• - maxit is the maximum number of outer iterations, if not given or set to [], then the default value min (10, N / restart) is used. Note that, if restart is empty, then maxit is the maximum number of iterations. If restart and maxit are not empty, then the maximum number of iterations is restart * maxit. If both restart and maxit are empty, then the maximum number of iterations is set to min (10, N).
• - M1, M2 are the preconditioners. The preconditioner M is given as M = M1 * M2. Both M1 and M2 can be passed as a matrix, function handle, or inline function g such that g(x) = M1 \ x or g(x) = M2 \ x. If M1 is [] or not given, then the preconditioner is not applied. The technique used is the left-preconditioning, i.e., it is solved inv(M) * A * x = inv(M) * b instead of A * x = b.
• - x0 is the initial guess, if not given or set to [], then the default value zeros (size (b)) is used.

The arguments which follow x0 are treated as parameters, and passed in a proper way to any of the functions (A or M or M1 or M2) which are passed to gmres.

The outputs are:

• - x the computed approximation. If the method does not converge, then it is the iterated with minimum residual.
• - flag indicates the exit status:

  0 : iteration converged to within the specified tolerance

  1 : maximum number of iterations exceeded

  2 : the preconditioner matrix is singular

  3 : algorithm reached stagnation (the relative difference between two

  consecutive iterations is less than eps)

• - relres is the value of the relative preconditioned residual of the approximation x.
• - iter is a vector containing the number of outer iterations and inner iterations performed to compute x. That is:
  • iter(1): number of outer iterations, i.e., how many times the method restarted. (if restart is empty or N, then it is 1, if not 1 ≤ iter(1) ≤ maxit).
  • iter(2): the number of iterations performed before the restart, i.e., the method restarts when iter(2) = restart. If restart is empty or N, then 1 ≤ iter(2) ≤ maxit.

  To be more clear, the approximation x is computed at the iteration (iter(1) - 1) * restart + iter(2). Since the output x corresponds to the minimal preconditioned residual solution, the total number of iterations that the method performed is given by length (resvec) - 1.

• - resvec is a vector containing the preconditioned relative residual at each iteration, including the 0-th iteration norm (A * x0 - b).

Let us consider a trivial problem with a tridiagonal matrix

n = 20;
A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n))  + ...
    toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
    sparse (1, 2, 1, 1, n) * n / 2);
b = A * ones (n, 1);
restart = 5;
[M1, M2] = ilu (A); # in this tridiag case, it corresponds to lu (A)
M = M1 * M2;
Afun = @(x) A * x;
Mfun = @(x) M \ x;
M1fun = @(x) M1 \ x;
M2fun = @(x) M2 \ x;

EXAMPLE 1: simplest usage of gmres

x = gmres (A, b, [], [], n)

EXAMPLE 2: gmres with a function which computes A * x

x = gmres (Afun, b, [], [], n)

EXAMPLE 3: usage of gmres with the restart

x = gmres (A, b, restart);

EXAMPLE 4: gmres with a preconditioner matrix M with and without restart

x = gmres (A, b, [], 1e-06, n, M)
x = gmres (A, b, restart, 1e-06, n, M)

EXAMPLE 5: gmres with a function as preconditioner

x = gmres (Afun, b, [], 1e-6, n, Mfun)

EXAMPLE 6: gmres with preconditioner matrices M1 and M2

x = gmres (A, b, [], 1e-6, n, M1, M2)

EXAMPLE 7: gmres with functions as preconditioners

x = gmres (Afun, b, 1e-6, n, M1fun, M2fun)

EXAMPLE 8: gmres with as input a function requiring an argument

function y = Ap (A, x, p) # compute A^p * x
     y = x;
     for i = 1:p
       y = A * y;
     endfor
  endfunction
Apfun = @(x, p) Ap (A, x, p);
x = gmres (Apfun, b, [], [], [], [], [], [], 2);

EXAMPLE 9: explicit example to show that gmres uses a left preconditioner

[M1, M2] = ilu (A + 0.1 * eye (n)); # factorization of A perturbed
M = M1 * M2;

## reference solution computed by gmres after two iterations
[x_ref, fl] = gmres (A, b, [], [], 1, M)

## left preconditioning
[x, fl] = gmres (M \ A, M \ b, [], [], 1)
x # compare x and x_ref

References:

1. Y. Saad, Iterative Methods for Sparse Linear Systems, Second edition, 2003, SIAM

See also: bicg, bicgstab, cgs, pcg, pcr, qmr, tfqmr.

x = qmr (A, b, rtol, maxit, M1, M2, x0)

x = qmr (A, b, rtol, maxit, P)

[x, flag, relres, iter, resvec] = qmr (A, b, …)

Solve A x = b using the Quasi-Minimal Residual iterative method (without look-ahead).

• - rtol is the relative tolerance, if not given or set to [] the default value 1e-6 is used.
• - maxit the maximum number of outer iterations, if not given or set to [] the default value min (20, numel (b)) is used.
• - x0 the initial guess, if not given or set to [] the default value zeros (size (b)) is used.

A can be passed as a matrix or as a function handle or inline function f such that f(x, "notransp") = A*x and f(x, "transp") = A'*x.

The preconditioner P is given as P = M1 * M2. Both M1 and M2 can be passed as a matrix or as a function handle or inline function g such that g(x, "notransp") = M1 \ x or g(x, "notransp") = M2 \ x and g(x, "transp") = M1' \ x or g(x, "transp") = M2' \ x.

If called with more than one output parameter

• - flag indicates the exit status:
  • - 0: iteration converged to the within the chosen tolerance
  • - 1: the maximum number of iterations was reached before convergence
  • - 3: the algorithm reached stagnation

  (the value 2 is unused but skipped for compatibility).

• - relres is the final value of the relative residual.
• - iter is the number of iterations performed.
• - resvec is a vector containing the residual norms at each iteration.

References:

1. R. Freund and N. Nachtigal, QMR: a quasi-minimal residual method for non-Hermitian linear systems, Numerische Mathematik, 1991, 60, pp. 315–339.
2. R. Barrett, M. Berry, T. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhour, R. Pozo, C. Romine, and H. van der Vorst, Templates for the solution of linear systems: Building blocks for iterative methods, SIAM, 2nd ed., 1994.

See also: bicg, bicgstab, cgs, gmres, pcg.

x = tfqmr (A, b, tol, maxit, M1, M2, x0, …)

x = tfqmr (A, b, tol, maxit, M, [], x0, …)

[x, flag, relres, iter, resvec] = tfqmr (A, b, …)

Solve A x = b using the Transpose-Tree qmr method, based on the cgs.

The input parameters are:

• - A is the matrix of the linear system and it must be square. A can be passed as a matrix, function handle, or inline function Afun such that Afun(x) = A * x. Additional parameters to Afun are passed after x0.
• - b is the right hand side vector. It must be a column vector with the same number of rows as A.
• - tol is the relative tolerance, if not given or set to [] the default value 1e-6 is used.
• - maxit the maximum number of outer iterations, if not given or set to [] the default value min (20, numel (b)) is used. To be compatible, since the method as different behaviors in the iteration number is odd or even, is considered as iteration in tfqmr the entire odd-even cycle. That is, to make an entire iteration, the algorithm performs two sub-iterations: the odd one and the even one.
• - M1, M2 are the preconditioners. The preconditioner M is given as M = M1 * M2. Both M1 and M2 can be passed as a matrix or as a function handle or inline function g such that g(x) = M1 \ x or g(x) = M2 \ x. The technique used is the right-preconditioning, i.e., it is solved A*inv(M)*y = b and then x = inv(M)*y, instead of A x = b.
• - x0 the initial guess, if not given or set to [] the default value zeros (size (b)) is used.

The arguments which follow x0 are treated as parameters, and passed in a proper way to any of the functions (A or M) which are passed to tfqmr.

The output parameters are:

• - x is the approximation computed. If the method doesn’t converge then it is the iterated with the minimum residual.
• - flag indicates the exit status:
  • - 0: iteration converged to the within the chosen tolerance
  • - 1: the maximum number of iterations was reached before convergence
  • - 2: the preconditioner matrix is singular
  • - 3: the algorithm reached stagnation
  • - 4: the algorithm can’t continue due to a division by zero
• - relres is the relative residual obtained as (A*x-b) / norm (b).
• - iter is the iteration which x is computed.
• - resvec is a vector containing the residual at each iteration (including norm (b - A x0)). Doing length (resvec) - 1 is possible to see the total number of iterations performed.

Let us consider a trivial problem with a tridiagonal matrix

n = 20;
A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n))  + ...
    toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
    sparse (1, 2, 1, 1, n) * n / 2);
b = A * ones (n, 1);
restart = 5;
[M1, M2] = ilu (A); # in this tridiag case it corresponds to chol (A)'
M = M1 * M2;
Afun = @(x) A * x;
Mfun = @(x) M \ x;
M1fun = @(x) M1 \ x;
M2fun = @(x) M2 \ x;

EXAMPLE 1: simplest usage of tfqmr

x = tfqmr (A, b, [], n)

EXAMPLE 2: tfqmr with a function which computes A * x

x = tfqmr (Afun, b, [], n)

EXAMPLE 3: tfqmr with a preconditioner matrix M

x = tfqmr (A, b, [], 1e-06, n, M)

EXAMPLE 4: tfqmr with a function as preconditioner

x = tfqmr (Afun, b, 1e-6, n, Mfun)

EXAMPLE 5: tfqmr with preconditioner matrices M1 and M2

x = tfqmr (A, b, [], 1e-6, n, M1, M2)

EXAMPLE 6: tfmqr with functions as preconditioners

x = tfqmr (Afun, b, 1e-6, n, M1fun, M2fun)

EXAMPLE 7: tfqmr with as input a function requiring an argument

function y = Ap (A, x, z) # compute A^z * x
   y = x;
   for i = 1:z
     y = A * y;
   endfor
 endfunction
Apfun = @(x, string, p) Ap (A, x, string, p);
x = tfqmr (Apfun, b, [], [], [], [], [], 2);

EXAMPLE 8: explicit example to show that tfqmr uses a right preconditioner

[M1, M2] = ilu (A + 0.3 * eye (n)); # factorization of A perturbed
M = M1 * M2;

## reference solution computed by tfqmr after one iteration
[x_ref, fl] = tfqmr (A, b, [], 1, M)

## right preconditioning
[y, fl] = tfqmr (A / M, b, [], 1)
x = M \ y # compare x and x_ref

References:

1. Y. Saad, Iterative Methods for Sparse Linear Systems, Second edition, 2003, SIAM

See also: bicg, bicgstab, cgs, gmres, pcg, qmr, pcr.