Previous: Functions of a Matrix, Up: Linear Algebra [Contents][Index]
Solve the linear system of equations A * x = b
by means of the Bi-Conjugate Gradient iterative method.
The input arguments are:
Afun
such that Afun (x, "notransp") = A * x
and
Afun (x, "transp") = A' * x
. Additional parameters to
Afun
may be passed after x0.
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.
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
.
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:
A * x = b
. If the algorithm did not converge,
then x is the iteration which has the minimum residual.
eps * norm (x,2)
.
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:
Solve A x = b
using the stabilizied Bi-conjugate gradient iterative
method.
The input parameters are:
Afun
such that Afun(x) = A * x
. Additional
parameters to Afun
are passed after x0.
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.
min (20, numel (b))
is used.
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
.
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:
(A*x-b) / norm(b)
.
iter + 0.5
(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:
Solve A x = b
, where A is a square matrix, using the
Conjugate Gradients Squared method.
The input arguments are:
Afun
such that Afun(x) = A * x
. Additional
parameters to Afun
are passed after x0.
min (20, numel (b))
is used.
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
.
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:
(A*x-b) / norm(b)
.
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:
Solve A x = b
using the Preconditioned GMRES iterative method with
restart, a.k.a. PGMRES(restart).
The input arguments are:
Afun
such that Afun(x) = A * x
. Additional
parameters to Afun
are passed after x0.
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.
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)
.
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
.
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:
consecutive iterations is less than eps)
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
.
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:
Solve A x = b
using the Quasi-Minimal Residual iterative method
(without look-ahead).
min (20, numel (b))
is used.
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
(the value 2 is unused but skipped for compatibility).
References:
Solve A x = b
using the Transpose-Tree qmr method, based on the cgs.
The input parameters are:
Afun
such that Afun(x) = A * x
. Additional
parameters to Afun
are passed after x0.
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.
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
.
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:
(A*x-b) / norm (b)
.
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:
Previous: Functions of a Matrix, Up: Linear Algebra [Contents][Index]