- :
`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

by means of the Bi-Conjugate Gradient iterative method.`A`*`x`=`b`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`Afcn`

such that`Afcn (x, "notransp") = A * x`

and`Afcn (x, "transp") = A' * x`

. Additional parameters to`Afcn`

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,

. The iteration stops if`b`-`A`*`x``norm (`

≤`b`-`A`*`x`)

. If`tol`* norm (`b`)`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

. Both`M`=`M1`*`M2``M1`and`M2`can be passed as a matrix or as a function handle or inline function`g`

such that`g (`

or`x`, "notransp") =`M1`\`x``g (`

and`x`, "notransp") =`M2`\`x``g (`

or`x`, "transp") =`M1`' \`x``g (`

. If`x`, "transp") =`M2`' \`x``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 (`

and`M1`) * A * inv (`M2`) *`y`= inv (`M1`) *`b``inv (`

and then setting`M2'`) * A' * inv (`M1'`) *`z`= inv (`M2'`) *`b`

.`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 (`Afcn`or`Mfcn`) or that have been given to`bicg`

.The output parameters are:

`x`is the computed approximation to the solution of

. If the algorithm did not converge, then`A`*`x`=`b``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; Afcn = @(x, string) strcmp (string, "notransp") * (A * x) + ... strcmp (string, "transp") * (A' * x); Mfcn = @(x, string) strcmp (string, "notransp") * (M \ x) + ... strcmp (string, "transp") * (M' \ x); M1fcn = @(x, string) strcmp (string, "notransp") * (M1 \ x) + ... strcmp (string, "transp") * (M1' \ x); M2fcn = @(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

and`A`*`x``A'`*`x`x = bicg (Afcn, 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 preconditionerx = bicg (Afcn, b, 1e-6, n, Mfcn)

EXAMPLE 5:

`bicg`

with preconditioner matrices`M1`and`M2`x = bicg (A, b, 1e-6, n, M1, M2)

EXAMPLE 6:

`bicg`

with functions as preconditionersx = bicg (Afcn, b, 1e-6, n, M1fcn, M2fcn)

EXAMPLE 7:

`bicg`

with as input a function requiring an argumentfunction 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 Apfcn = @(x, string, p) Ap (A, x, string, p); x = bicg (Apfcn, b, [], [], [], [], [], 2);

Reference:

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

- :
`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`Afcn`

such that`Afcn(x) = A * x`

. Additional parameters to`Afcn`

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,

. The iteration stops if`b`-`A`*`x``norm (`

≤`b`-`A`*`x`)

. If`tol`* norm (`b`)`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

. Both`M`=`M1`*`M2``M1`and`M2`can be passed as a matrix or as a function handle or inline function`g`

such that`g(`

or`x`) =`M1`\`x``g(`

. The technique used is the right preconditioning, i.e., it is solved`x`) =`M2`\`x`

and then`A`* inv (`M`) *`y`=`b`

.`x`= inv (`M`) *`y``x0`the initial guess, if not given or set to [] the default value`zeros (size (`

is used.`b`))

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(`

is possible to see the total number of (total) iterations performed.`resvec`) - 1) / 2

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; Afcn = @(x) A * x; Mfcn = @(x) M \ x; M1fcn = @(x) M1 \ x; M2fcn = @(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 (Afcn, 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 preconditionerx = bicgstab (Afcn, b, 1e-6, n, Mfcn)

EXAMPLE 5:

`bicgstab`

with preconditioner matrices`M1`and`M2`x = bicgstab (A, b, [], 1e-6, n, M1, M2)

EXAMPLE 6:

`bicgstab`

with functions as preconditionersx = bicgstab (Afcn, b, 1e-6, n, M1fcn, M2fcn)

EXAMPLE 7:

`bicgstab`

with as input a function requiring an argumentfunction y = Ap (A, x, z) # compute A^z * x y = x; for i = 1:z y = A * y; endfor endfunction Apfcn = @(x, string, p) Ap (A, x, string, p); x = bicgstab (Apfcn, 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

Reference:

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

- :
`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`Afcn`

such that`Afcn(x) = A * x`

. Additional parameters to`Afcn`

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

and then`A`*inv(`M`)*y = b

.`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(`

is possible to see the total number of iterations performed.`resvec`) - 1

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; Afcn = @(x) A * x; Mfcn = @(x) M \ x; M1fcn = @(x) M1 \ x; M2fcn = @(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 (Afcn, 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 preconditionerx = cgs (Afcn, b, 1e-6, n, Mfcn)

EXAMPLE 5:

`cgs`

with preconditioner matrices`M1`and`M2`x = cgs (A, b, [], 1e-6, n, M1, M2)

EXAMPLE 6:

`cgs`

with functions as preconditionersx = cgs (Afcn, b, 1e-6, n, M1fcn, M2fcn)

EXAMPLE 7:

`cgs`

with as input a function requiring an argumentfunction y = Ap (A, x, z) # compute A^z * x y = x; for i = 1:z y = A * y; endfor endfunction Apfcn = @(x, string, p) Ap (A, x, string, p); x = cgs (Apfcn, 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:

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

- :
`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`Afcn`

such that`Afcn(x) = A * x`

. Additional parameters to`Afcn`

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 (`

. The iteration stops if`M`) * (`b`-`a`*`x`)`norm (inv (`

. If`M`) * (`b`-`a`*`x`)) ≤`tol`* norm (inv (`M`) *`B`)`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,`

is used. Note that, if`N`/`restart`)`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

. If both`restart`*`maxit``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(`

instead of`M`) *`A`*`x`= inv(`M`) *`b`

.`A`*`x`=`b``x0`is the initial guess, if not given or set to [], then the default value`zeros (size (`

is used.`b`))

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

. If`iter(2)`=`restart``restart`is empty or`N`, then 1 ≤`iter(2)`≤`maxit`.

To be more clear, the approximation

`x`is computed at the iteration`(`

. Since the output`iter(1)`- 1) *`restart`+`iter(2)``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; Afcn = @(x) A * x; Mfcn = @(x) M \ x; M1fcn = @(x) M1 \ x; M2fcn = @(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 (Afcn, b, [], [], n)

EXAMPLE 3: usage of

`gmres`

with the restartx = gmres (A, b, restart);

EXAMPLE 4:

`gmres`

with a preconditioner matrix`M`with and without restartx = gmres (A, b, [], 1e-06, n, M) x = gmres (A, b, restart, 1e-06, n, M)

EXAMPLE 5:

`gmres`

with a function as preconditionerx = gmres (Afcn, b, [], 1e-6, n, Mfcn)

EXAMPLE 6:

`gmres`

with preconditioner matrices`M1`and`M2`x = gmres (A, b, [], 1e-6, n, M1, M2)

EXAMPLE 7:

`gmres`

with functions as preconditionersx = gmres (Afcn, b, 1e-6, n, M1fcn, M2fcn)

EXAMPLE 8:

`gmres`

with as input a function requiring an argumentfunction y = Ap (A, x, p) # compute A^p * x y = x; for i = 1:p y = A * y; endfor endfunction Apfcn = @(x, p) Ap (A, x, p); x = gmres (Apfcn, 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

Reference:

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

- :
`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:

- R. Freund and N. Nachtigal, QMR: a quasi-minimal residual method for non-Hermitian linear systems, Numerische Mathematik, 1991, 60, pp. 315–339.
- 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.

- :
`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`Afcn`

such that`Afcn(x) = A * x`

. Additional parameters to`Afcn`

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 (`

is possible to see the total number of iterations performed.`resvec`) - 1

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; Afcn = @(x) A * x; Mfcn = @(x) M \ x; M1fcn = @(x) M1 \ x; M2fcn = @(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 (Afcn, 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 preconditionerx = tfqmr (Afcn, b, 1e-6, n, Mfcn)

EXAMPLE 5:

`tfqmr`

with preconditioner matrices`M1`and`M2`x = tfqmr (A, b, [], 1e-6, n, M1, M2)

EXAMPLE 6:

`tfmqr`

with functions as preconditionersx = tfqmr (Afcn, b, 1e-6, n, M1fcn, M2fcn)

EXAMPLE 7:

`tfqmr`

with as input a function requiring an argumentfunction y = Ap (A, x, z) # compute A^z * x y = x; for i = 1:z y = A * y; endfor endfunction Apfcn = @(x, string, p) Ap (A, x, string, p); x = tfqmr (Apfcn, 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

Reference:

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