load watt1
load west0479
A = watt1;
As = (A + A')/2;
Au = (A - A')/2;
norm(As,'fro')
ans =
11.3137
norm(As,'fro') / norm(A,'fro')
ans =
1.0000
norm(Au, 'fro')
ans =
5.4619e-007
norm(Au, 'fro') / norm(A,'fro')
ans =
4.8277e-008
[E,V] = eig(full(A));
nnz(E)
ans =
3207284
nnz(V)
ans =
1856
v = diag(V);
plot(v)
whos
Name Size Bytes Class
A 1856x1856 143748 double array (sparse)
As 1856x1856 280068 double array (sparse)
Au 1856x1856 280068 double array (sparse)
E 1856x1856 27557888 double array
V 1856x1856 27557888 double array
ans 1x1 8 double array
v 1856x1 14848 double array
watt1 1856x1856 143748 double array (sparse)
west0479 479x479 24840 double array (sparse)
Grand total is 6961399 elements using 56003104 bytes
max(v)
ans =
1
min(v)
ans =
-1.3855e-006
hist(v)
hist(log10(v))
??? Edges vector must be monotonically non-decreasing.
Error in ==> hist at 79
nn = histc(y,[-inf bins],1);
hist(log10(abs(v)))
nnz( v > 0)
ans =
128
nnz( v < 0)
ans =
1728
n = max(size(A))
n =
1856
b = ones(n,1);
[x,flag,relres,iter,resvec] = bicgstab(A,b,tol,maxit);
??? Undefined function or variable 'tol'.
tol = 10^-12
tol =
1.0000e-012
maxit = 1000
maxit =
1000
[x,flag,relres,iter,resvec] = bicgstab(A,b,tol,maxit);
flag
flag =
3
help bicgstab
BICGSTAB BiConjugate Gradients Stabilized Method.
X = BICGSTAB(A,B) attempts to solve the system of linear equations
A*X=B for X. The N-by-N coefficient matrix A must be square and the
right hand side column vector B must have length N.%
X = BICGSTAB(AFUN,B) accepts a function handle AFUN instead of the
matrix A. AFUN(X) accepts a vector input X and returns the
matrix-vector product A*X. In all of the following syntaxes, you can
replace A by AFUN.
X = BICGSTAB(A,B,TOL) specifies the tolerance of the method. If TOL is
[] then BICGSTAB uses the default, 1e-6.
X = BICGSTAB(A,B,TOL,MAXIT) specifies the maximum number of iterations.
If MAXIT is [] then BICGSTAB uses the default, min(N,20).
X = BICGSTAB(A,B,TOL,MAXIT,M) and X = BICGSTAB(A,B,TOL,MAXIT,M1,M2) use
preconditioner M or M=M1*M2 and effectively solve the system
inv(M)*A*X = inv(M)*B for X. If M is [] then a preconditioner is not
applied. M may be a function handle returning M\X.
X = BICGSTAB(A,B,TOL,MAXIT,M1,M2,X0) specifies the initial guess. If
X0 is [] then BICGSTAB uses the default, an all zero vector.
[X,FLAG] = BICGSTAB(A,B,...) also returns a convergence FLAG:
0 BICGSTAB converged to the desired tolerance TOL within MAXIT iterations.
1 BICGSTAB iterated MAXIT times but did not converge.
2 preconditioner M was ill-conditioned.
3 BICGSTAB stagnated (two consecutive iterates were the same).
4 one of the scalar quantities calculated during BICGSTAB became
too small or too large to continue computing.
[X,FLAG,RELRES] = BICGSTAB(A,B,...) also returns the relative residual
NORM(B-A*X)/NORM(B). If FLAG is 0, then RELRES <= TOL.
[X,FLAG,RELRES,ITER] = BICGSTAB(A,B,...) also returns the iteration
number at which X was computed: 0 <= ITER <= MAXIT. ITER may be an
integer + 0.5, indicating convergence half way through an iteration.
[X,FLAG,RELRES,ITER,RESVEC] = BICGSTAB(A,B,...) also returns a vector
of the residual norms at each half iteration, including NORM(B-A*X0).
Example:
n = 21; A = gallery('wilk',n); b = sum(A,2);
tol = 1e-12; maxit = 15; M = diag([10:-1:1 1 1:10]);
x = bicgstab(A,b,tol,maxit,M);
Or, use this matrix-vector product function
%-----------------------------------------------------------------%
function y = afun(x,n)
y = [0; x(1:n-1)] + [((n-1)/2:-1:0)'; (1:(n-1)/2)'].*x+[x(2:n); 0];
%-----------------------------------------------------------------%
and this preconditioner backsolve function
%------------------------------------------%
function y = mfun(r,n)
y = r ./ [((n-1)/2:-1:1)'; 1; (1:(n-1)/2)'];
%------------------------------------------%
as inputs to BICGSTAB:
x1 = bicgstab(@(x)afun(x,n),b,tol,maxit,@(x)mfun(x,n));
Class support for inputs A,B,M1,M2,X0 and the output of AFUN:
float: double
See also bicg, cgs, gmres, lsqr, minres, pcg, qmr, symmlq, luinc,
function_handle.
Reference page in Help browser
doc bicgstab
semilogy(resvec / resvec(0))
??? Subscript indices must either be real positive integers or logicals.
semilogy(resvec / resvec(1))
condest(A)
ans =
5.3847e+009
tol = 10^-10
tol =
1.0000e-010
[x,flag,relres,iter,resvec] = bicgstab(A,b,tol,maxit);
flag
flag =
0
norm(A*x-b)
ans =
2.4401e-009
semilogy(resvec / resvec(1))
help gmres
GMRES Generalized Minimum Residual Method.
X = GMRES(A,B) attempts to solve the system of linear equations A*X = B
for X. The N-by-N coefficient matrix A must be square and the right
hand side column vector B must have length N. This uses the unrestarted
method with MIN(N,10) total iterations.
X = GMRES(AFUN,B) accepts a function handle AFUN instead of the matrix
A. AFUN(X) accepts a vector input X and returns the matrix-vector
product A*X. In all of the following syntaxes, you can replace A by
AFUN.
X = GMRES(A,B,RESTART) restarts the method every RESTART iterations.
If RESTART is N or [] then GMRES uses the unrestarted method as above.
X = GMRES(A,B,RESTART,TOL) specifies the tolerance of the method. If
TOL is [] then GMRES uses the default, 1e-6.
X = GMRES(A,B,RESTART,TOL,MAXIT) specifies the maximum number of outer
iterations. Note: the total number of iterations is RESTART*MAXIT. If
MAXIT is [] then GMRES uses the default, MIN(N/RESTART,10). If RESTART
is N or [] then the total number of iterations is MAXIT.
X = GMRES(A,B,RESTART,TOL,MAXIT,M) and
X = GMRES(A,B,RESTART,TOL,MAXIT,M1,M2) use preconditioner M or M=M1*M2
and effectively solve the system inv(M)*A*X = inv(M)*B for X. If M is
[] then a preconditioner is not applied. M may be a function handle
returning M\X.
X = GMRES(A,B,RESTART,TOL,MAXIT,M1,M2,X0) specifies the first initial
guess. If X0 is [] then GMRES uses the default, an all zero vector.
[X,FLAG] = GMRES(A,B,...) also returns a convergence FLAG:
0 GMRES converged to the desired tolerance TOL within MAXIT iterations.
1 GMRES iterated MAXIT times but did not converge.
2 preconditioner M was ill-conditioned.
3 GMRES stagnated (two consecutive iterates were the same).
[X,FLAG,RELRES] = GMRES(A,B,...) also returns the relative residual
NORM(B-A*X)/NORM(B). If FLAG is 0, then RELRES <= TOL. Note with
preconditioners M1,M2, the residual is NORM(M2\(M1\(B-A*X))).
[X,FLAG,RELRES,ITER] = GMRES(A,B,...) also returns both the outer and
inner iteration numbers at which X was computed: 0 <= ITER(1) <= MAXIT
and 0 <= ITER(2) <= RESTART.
[X,FLAG,RELRES,ITER,RESVEC] = GMRES(A,B,...) also returns a vector of
the residual norms at each inner iteration, including NORM(B-A*X0).
Note with preconditioners M1,M2, the residual is NORM(M2\(M1\(B-A*X))).
Example:
n = 21; A = gallery('wilk',n); b = sum(A,2);
tol = 1e-12; maxit = 15; M = diag([10:-1:1 1 1:10]);
x = gmres(A,b,10,tol,maxit,M);
Or, use this matrix-vector product function
%-----------------------------------------------------------------%
function y = afun(x,n)
y = [0; x(1:n-1)] + [((n-1)/2:-1:0)'; (1:(n-1)/2)'].*x+[x(2:n); 0];
%-----------------------------------------------------------------%
and this preconditioner backsolve function
%------------------------------------------%
function y = mfun(r,n)
y = r ./ [((n-1)/2:-1:1)'; 1; (1:(n-1)/2)'];
%------------------------------------------%
as inputs to GMRES:
x1 = gmres(@(x)afun(x,n),b,10,tol,maxit,@(x)mfun(x,n));
Class support for inputs A,B,M1,M2,X0 and the output of AFUN:
float: double
See also bicg, bicgstab, cgs, lsqr, minres, pcg, qmr, symmlq, luinc,
function_handle.
Reference page in Help browser
doc gmres
[x,flag,relres,iter,resvec] = gmres(A,b,[],tol,2000);
flag
flag =
0
iter
iter =
1 320
hold on
semilogy(resvec / resvec(1),'r')
resvec(end-2 : end)
ans =
1.0e-004 *
0.0001
0.0000
0.3476
relres
relres =
8.0691e-007
resvec(end-1)
ans =
4.5306e-009
[x,flag,relres,iter,resvec] = gmres(A,b,20,tol,2000);
iter
iter =
2000 3
semilogy(resvec / resvec(1),'g')
[x,flag,relres,iter,resvec] = gmres(A,b,50,tol,2000);
Error in ==> gmres at 316
Jtemp = speye(n);
[x,flag,relres,iter,resvec] = gmres(A,b,50,tol,200);
flag
flag =
0
iter
iter =
200 2
semilogy(resvec / resvec(1),'c')
close all
[x,flag,relres,iter,resvec] = bicgstab(A,b,tol,maxit);
semilogy(resvec / resvec(1))
hold on
[L,U,P] = luinc(A,'0');
figure, spy(A)
spy(L)
spy(U)
spy(P)
nnz(diag(P))
ans =
1856
n
n =
1856
AD = diag(diag(A));
size(AD)
ans =
1856 1856
AOD = A - AD;
sum(abs(AOD),2);
rs = sum(abs(AOD),2);
de = abs(diag(AD));
hist( de / rs)
size(de)
ans =
1856 1
size(rs)
ans =
1856 1
hist( de ./ rs)
Warning: Divide by zero.
hist( de ./ max(rs,1))
hist( rs ./ de)
figure
hist( rs ./ de)
e = diag(E);
size(e)
ans =
1856 1
plot(e,'.')
plot(log10(abs(e)),'.')
Warning: Log of zero.
> In log10 at 20
[L,U,P] = luinc(A,'0');
nnz(L+U)
ans =
11360
nnz(A)
ans =
11360
[x,flag,relres,iter,resvec] = bicgstab(A,b,tol,maxit,L,U);
iter
iter =
33
figure(1)
hold on
semilogy(resvec / resvec(1),'r')
help luinc
LUINC Sparse Incomplete LU factorization.
LUINC produces two different kinds of incomplete LU factorizations -- the
drop tolerance and the 0 level of fill-in factorizations. These factors
may be useful as preconditioners for a system of linear equations being
solved by an iterative method such as BICG (BiConjugate Gradients).
LUINC(X,DROPTOL) performs the incomplete LU factorization of
X with drop tolerance DROPTOL.
LUINC(X,OPTS) allows additional options to the incomplete LU
factorization. OPTS is a structure with up to four fields:
droptol --- the drop tolerance of incomplete LU
milu --- modified incomplete LU
udiag --- replace zeros on the diagonal of U
thresh --- the pivot threshold (see also LU)
Only the fields of interest need to be set.
droptol is a non-negative scalar used as the drop
tolerance for the incomplete LU factorization. This factorization
is computed in the same (column-oriented) manner as the
LU factorization except after each column of L and U has
been calculated, all entries in that column which are smaller
in magnitude than the local drop tolerance, which is
droptol * NORM of the column of X, are "dropped" from L or U.
The only exception to this dropping rule is the diagonal of the
upper triangular factor U which remains even if it is too small.
Note that entries of the lower triangular factor L are tested
before being scaled by the pivot. Setting droptol = 0
produces the complete LU factorization, which is the default.
milu stands for modified incomplete LU factorization. Its
value is either 0 (unmodified, the default) or 1 (modified).
Instead of discarding those entries from the newly-formed
column of the factors, they are subtracted from the diagonal
of the upper triangular factor, U.
udiag is either 0 or 1. If it is 1, any zero diagonal entries
of the upper triangular factor U are replaced by the local drop
tolerance in an attempt to avoid a singular factor. The default
is 0.
thresh is a pivot threshold in [0,1]. Pivoting occurs
when the diagonal entry in a column has magnitude less
than thresh times the magnitude of any sub-diagonal entry in
that column. thresh = 0 forces diagonal pivoting. thresh = 1 is
the default.
Example:
load west0479
A = west0479;
nnz(A)
nnz(lu(A))
nnz(luinc(A,1e-6))
This shows that A has 1887 nonzeros, its complete LU factorization
has 16777 nonzeros, and its incomplete LU factorization with a
drop tolerance of 1e-6 has 10311 nonzeros.
[L,U,P] = LUINC(X,'0') produces the incomplete LU factors of a sparse
matrix with 0 level of fill-in (i.e. no fill-in). L is unit lower
trianglar, U is upper triangular and P is a permutation matrix. U has the
same sparsity pattern as triu(P*X). L has the same sparsity pattern as
tril(P*X), except for 1's on the diagonal of L where P*X may be zero. Both
L and U may have a zero because of cancellation where P*X is nonzero. L*U
differs from P*X only outside of the sparsity pattern of P*X.
[L,U] = LUINC(X,'0') produces upper triangular U and L is a permutation of
unit lower triangular matrix. Thus no comparison can be made between the
sparsity patterns of L,U and X, although nnz(L) + nnz(U) = nnz(X) + n. L*U
differs from X only outside of its sparsity pattern.
LU = LUINC(X,'0') returns "L+U-I", where L is unit lower triangular, U is
upper triangular and the permutation information is lost.
Example:
load west0479
A = west0479;
[L,U,P] = luinc(A,'0');
isequal(spones(U),spones(triu(P*A)))
spones(L) ~= spones(tril(P*A))
D = (L*U) .* spones(P*A) - P*A
spones(L) differs from spones(tril(P*A)) at some positions on the
diagonal and at one position in L where cancellation zeroed out a
nonzero element of P*A. The entries of D are of the order of eps.
LUINC works only for sparse matrices.
See also lu, cholinc, bicg.
Reference page in Help browser
doc luinc
opt.milu = 0
opt =
milu: 0
opt.udiag = 0
opt =
milu: 0
udiag: 0
opt.thresh = 0
opt =
milu: 0
udiag: 0
thresh: 0
opt.droptol = .1
opt =
milu: 0
udiag: 0
thresh: 0
droptol: 0.1000
[L,U,P] = luinc(A,opt);
nnz(diag(P))
ans =
1856
nnz(L+)
??? nnz(L+)
|
Error: Unbalanced or misused parentheses or brackets.
nnz(L+U)
ans =
7699
nnz(A)
ans =
11360
figure,spy(A)
spy(L+U)
opt.droptol = .05
opt =
milu: 0
udiag: 0
thresh: 0
droptol: 0.0500
[L,U,P] = luinc(A,opt);
nnz(A)
ans =
11360
nnz(L+U)
ans =
10556
opt.droptol = .04
opt =
milu: 0
udiag: 0
thresh: 0
droptol: 0.0400
[L,U,P] = luinc(A,opt);
nnz(L+U)
ans =
11618
nnz(A)
ans =
11360
[x,flag,relres,iter,resvec] = bicgstab(A,b,tol,maxit,L,U);
figure 1
??? Error using ==> figure
Single input must be an existing figure handle or a scalar integer from 1 to 2030832096
figure(1)
hold on
semilogy(resvec / resvec(1),'g')
A = west0497
??? Undefined function or variable 'west0497'.
A = west0479;
n = length(A)
n =
479
close all
spy(A)
condest(A)
ans =
1.9326e+013
b = ones(n,1);
x = A \ b;
norm(A*x - b)
ans =
1.9912e-008
tol = 10^-6
tol =
1.0000e-006
[E,V] = eig(full(A));
nnz(V)
ans =
479
v = diag(V);
whos
Name Size Bytes Class
A 479x479 24840 double array (sparse)
AD 1856x1856 29700 double array (sparse)
AOD 1856x1856 166020 double array (sparse)
As 1856x1856 280068 double array (sparse)
Au 1856x1856 280068 double array (sparse)
E 479x479 3671056 double array (complex)
L 1856x1856 143748 double array (sparse)
P 1856x1856 29700 double array (sparse)
U 1856x1856 143748 double array (sparse)
V 479x479 3671056 double array (complex)
ans 1x1 8 double array
b 479x1 3832 double array
de 1856x1 22280 double array (sparse)
e 1856x1 14848 double array
flag 1x1 8 double array
iter 1x1 8 double array
maxit 1x1 8 double array
n 1x1 8 double array
opt 1x1 528 struct array
relres 1x1 8 double array
resvec 49x1 392 double array
rs 1856x1 20744 double array (sparse)
tol 1x1 8 double array
v 479x1 7664 double array (complex)
watt1 1856x1856 143748 double array (sparse)
west0479 479x479 24840 double array (sparse)
x 479x1 3832 double array
Grand total is 566091 elements using 8682768 bytes
figure, plot(v)
figure, plot(v,'.')
size(E)
ans =
479 479
condest(E)
ans =
4.0523e+010
[x,flag,relres,iter,resvec] = bicgstab(A,b,tol,maxit);
flag
flag =
1
iter
iter =
0
size(b)
ans =
479 1
tol
tol =
1.0000e-006
maxit
maxit =
1000
help bicgstab
BICGSTAB BiConjugate Gradients Stabilized Method.
X = BICGSTAB(A,B) attempts to solve the system of linear equations
A*X=B for X. The N-by-N coefficient matrix A must be square and the
right hand side column vector B must have length N.%
X = BICGSTAB(AFUN,B) accepts a function handle AFUN instead of the
matrix A. AFUN(X) accepts a vector input X and returns the
matrix-vector product A*X. In all of the following syntaxes, you can
replace A by AFUN.
X = BICGSTAB(A,B,TOL) specifies the tolerance of the method. If TOL is
[] then BICGSTAB uses the default, 1e-6.
X = BICGSTAB(A,B,TOL,MAXIT) specifies the maximum number of iterations.
If MAXIT is [] then BICGSTAB uses the default, min(N,20).
X = BICGSTAB(A,B,TOL,MAXIT,M) and X = BICGSTAB(A,B,TOL,MAXIT,M1,M2) use
preconditioner M or M=M1*M2 and effectively solve the system
inv(M)*A*X = inv(M)*B for X. If M is [] then a preconditioner is not
applied. M may be a function handle returning M\X.
X = BICGSTAB(A,B,TOL,MAXIT,M1,M2,X0) specifies the initial guess. If
X0 is [] then BICGSTAB uses the default, an all zero vector.
[X,FLAG] = BICGSTAB(A,B,...) also returns a convergence FLAG:
0 BICGSTAB converged to the desired tolerance TOL within MAXIT iterations.
1 BICGSTAB iterated MAXIT times but did not converge.
2 preconditioner M was ill-conditioned.
3 BICGSTAB stagnated (two consecutive iterates were the same).
4 one of the scalar quantities calculated during BICGSTAB became
too small or too large to continue computing.
[X,FLAG,RELRES] = BICGSTAB(A,B,...) also returns the relative residual
NORM(B-A*X)/NORM(B). If FLAG is 0, then RELRES <= TOL.
[X,FLAG,RELRES,ITER] = BICGSTAB(A,B,...) also returns the iteration
number at which X was computed: 0 <= ITER <= MAXIT. ITER may be an
integer + 0.5, indicating convergence half way through an iteration.
[X,FLAG,RELRES,ITER,RESVEC] = BICGSTAB(A,B,...) also returns a vector
of the residual norms at each half iteration, including NORM(B-A*X0).
Example:
n = 21; A = gallery('wilk',n); b = sum(A,2);
tol = 1e-12; maxit = 15; M = diag([10:-1:1 1 1:10]);
x = bicgstab(A,b,tol,maxit,M);
Or, use this matrix-vector product function
%-----------------------------------------------------------------%
function y = afun(x,n)
y = [0; x(1:n-1)] + [((n-1)/2:-1:0)'; (1:(n-1)/2)'].*x+[x(2:n); 0];
%-----------------------------------------------------------------%
and this preconditioner backsolve function
%------------------------------------------%
function y = mfun(r,n)
y = r ./ [((n-1)/2:-1:1)'; 1; (1:(n-1)/2)'];
%------------------------------------------%
as inputs to BICGSTAB:
x1 = bicgstab(@(x)afun(x,n),b,tol,maxit,@(x)mfun(x,n));
Class support for inputs A,B,M1,M2,X0 and the output of AFUN:
float: double
See also bicg, cgs, gmres, lsqr, minres, pcg, qmr, symmlq, luinc,
function_handle.
Reference page in Help browser
doc bicgstab
maxit
maxit =
1000
iter
iter =
0
resvec
resvec =
1.0e+018 *
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0000
0.0001
0.0000
0.0001
0.0001
0.0001
0.0000
0.0002
0.0002
0.0000
0.0000
0.0001
0.0001
0.0000
0.0000
0.0000
0.0000
0.0001
0.0001
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0002
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0004
0.0004
0.0001
0.0001
0.0021
0.0018
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0003
0.0002
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0001
0.0001
0.0001
0.0001
0.0000
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0020
0.0016
0.0002
0.0002
0.0002
0.0001
0.0000
0.0000
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0003
0.0003
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0009
0.0005
0.0001
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0000
0.0000
0.0000
0.0012
0.0008
0.0003
0.0003
0.0002
0.0001
0.0002
0.0002
0.0002
0.0002
0.0002
0.0002
0.0002
0.0002
0.0002
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0010
0.0004
0.0004
0.0004
0.0006
0.0006
0.0009
0.0009
0.0040
0.0026
0.0103
0.0102
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0002
0.0001
0.0002
0.0001
0.0002
0.0001
0.0005
0.0003
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0003
0.0003
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0001
0.0001
0.0000
0.0002
0.0002
0.0001
0.0000
0.0000
0.0000
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0001
0.0000
0.0000
0.0000
0.0000
0.0003
0.0003
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0001
0.0001
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0010
0.0005
0.0011
0.0010
0.0007
0.0007
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0000
0.0001
0.0001
0.0001
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0001
0.0001
0.0001
0.0055
0.0031
0.0002
0.0002
0.0002
0.0002
0.0006
0.0006
0.0003
0.0003
0.0003
0.0003
0.0005
0.0003
0.0004
0.0003
0.0009
0.0008
0.0001
0.0001
0.0010
0.0007
0.0007
0.0006
0.0013
0.0008
0.0008
0.0007
0.0007
0.0007
0.0008
0.0007
0.0007
0.0007
0.0012
0.0008
0.0071
0.0043
0.0046
0.0026
0.0026
0.0024
0.0027
0.0025
0.0013
0.0011
0.0160
0.0143
0.0102
0.0100
0.0008
0.0008
0.0019
0.0008
0.0008
0.0008
0.0112
0.0068
0.0067
0.0065
0.0075
0.0045
0.0050
0.0046
0.0047
0.0045
0.0093
0.0070
0.2162
0.1768
0.1662
0.1408
0.0224
0.0224
0.0174
0.0173
0.0184
0.0086
0.0102
0.0080
0.0070
0.0069
0.0072
0.0070
0.0135
0.0129
0.0463
0.0311
0.1095
0.1060
0.0099
0.0083
0.1190
0.1047
0.0214
0.0197
0.0636
0.0358
0.0514
0.0249
0.0697
0.0611
0.0115
0.0096
0.0102
0.0102
0.0157
0.0099
0.0120
0.0111
0.0585
0.0291
0.0047
0.0027
0.0027
0.0027
0.0049
0.0039
0.0036
0.0009
0.0021
0.0013
0.0014
0.0013
0.0014
0.0014
0.0089
0.0087
0.0033
0.0022
0.0018
0.0017
0.0019
0.0016
0.0016
0.0015
0.0041
0.0031
0.0022
0.0019
0.0034
0.0032
0.0015
0.0015
0.0013
0.0012
0.0019
0.0012
0.0096
0.0095
0.0106
0.0088
0.0071
0.0055
0.0294
0.0288
0.0053
0.0051
0.0490
0.0314
0.0223
0.0120
0.0157
0.0155
0.0113
0.0079
0.0080
0.0079
0.0012
0.0011
0.0002
0.0002
0.0013
0.0008
0.0007
0.0004
0.0023
0.0023
0.0005
0.0004
0.0005
0.0004
0.0004
0.0004
0.0004
0.0004
0.0375
0.0374
0.0004
0.0004
0.0005
0.0004
0.0013
0.0007
0.0007
0.0007
0.0006
0.0006
0.0006
0.0005
0.0005
0.0005
0.0004
0.0003
0.0003
0.0003
0.0003
0.0003
0.0006
0.0006
0.0016
0.0013
0.0004
0.0004
0.0004
0.0004
0.0004
0.0004
0.0003
0.0003
0.0003
0.0003
0.0003
0.0003
0.0003
0.0003
0.0005
0.0004
0.0011
0.0008
0.0056
0.0045
0.0052
0.0051
0.0056
0.0056
0.0393
0.0391
0.0009
0.0008
0.0009
0.0008
0.0018
0.0017
0.0017
0.0015
0.0019
0.0017
0.0017
0.0017
0.0017
0.0017
0.0026
0.0026
0.0018
0.0018
0.0004
0.0004
0.0004
0.0004
0.0004
0.0004
0.0001
0.0001
0.0002
0.0001
0.0006
0.0005
0.0001
0.0001
0.0002
0.0002
0.0002
0.0002
0.0002
0.0002
0.0002
0.0002
0.0002
0.0002
0.0003
0.0003
0.0011
0.0010
0.0001
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0002
0.0001
0.0001
0.0000
0.0000
0.0000
0.0001
0.0001
0.0001
0.0001
0.0001
0.0000
0.0001
0.0001
0.0000
0.0000
0.0001
0.0000
0.0000
0.0000
0.0001
0.0000
0.0001
0.0000
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0003
0.0003
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0003
0.0002
0.0002
0.0002
0.0003
0.0003
0.0002
0.0002
0.0002
0.0002
0.0000
0.0000
0.0000
0.0000
0.0002
0.0001
0.0001
0.0000
0.0000
0.0000
0.0001
0.0000
0.0000
0.0000
0.0001
0.0000
0.0001
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0006
0.0005
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0016
0.0008
0.0001
0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0063
0.0052
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0002
0.0001
4.3044
1.9697
0.0387
0.0381
0.0003
0.0003
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0003
0.0002
0.0004
0.0003
0.0003
0.0003
0.0000
0.0000
0.0000
0.0000
0.0002
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
0.0001
figure
semilogy(resvec / resvec(1))
resvec(1)
ans =
21.8861
norm(b)
ans =
21.8861
sqrt(n)
ans =
21.8861
[x,flag,relres,iter,resvec] = gmres(A,b,[],tol,maxit);
??? Attempted to access w(480); index out of bounds because numel(w)=479.
Error in ==> gmres at 327
normr = abs(w(initer+1));
[x,flag,relres,iter,resvec] = gmres(A,b,[],tol,479);
flag
flag =
1
iter
iter =
1 479
figure(1)
figure(3)
hold on
semilogy(resvec / resvec(1),'r')
[L,U,P] = luinc(A,'0');
nnz(L+U)
ans =
1893
nnz(A)
ans =
1910
figure(1)
spy(P)
[x,flag,relres,iter,resvec] = gmres(P*A,P*b,[],tol,479,L,U);
Warning: Matrix is singular to working precision.
> In sparfun\private\iterapp at 35
In gmres at 226
n
n =
479
spy(L)
spy(U)
opt
opt =
milu: 0
udiag: 0
thresh: 0
droptol: 0.0400
opt.udiag = 1
opt =
milu: 0
udiag: 1
thresh: 0
droptol: 0.0400
opt.thresh = 1
opt =
milu: 0
udiag: 1
thresh: 1
droptol: 0.0400
[L,U,P] = luinc(A,opt);
Warning: Incomplete upper triangular factor had 35 zero diagonals replaced
by local drop tolerance.
It may not make a good preconditioner for an iterative method.
[x,flag,relres,iter,resvec] = gmres(P*A,P*b,[],tol,479,L,U);
flag
flag =
0
iter
iter =
1 9
figure(3)
semilogy(resvec / resvec(1),'g')
resvec
resvec =
1.0e+006 *
4.6781
4.6780
4.6669
3.5960
3.5640
3.5631
3.4687
0.7603
0.6700
0.6728
resvec / resvec(1)
ans =
1.0000
1.0000
0.9976
0.7687
0.7618
0.7617
0.7415
0.1625
0.1432
0.1438
tol
tol =
1.0000e-006
nnz(A)
ans =
1910
nnz(L + U)
ans =
2279
norm(A*x - b)
ans =
1.2390e+005
norm(b)
ans =
21.8861
norm (P*A*x - P*b)
ans =
1.2390e+005
[x,flag,relres,iter,resvec] = bicgstab(P*A,P*b,tol,479,L,U);
flag
flag =
1
iter
iter =
0
tol
tol =
1.0000e-006
opt
opt =
milu: 0
udiag: 1
thresh: 1
droptol: 0.0400
opt.droptol = .01
opt =
milu: 0
udiag: 1
thresh: 1
droptol: 0.0100
[L,U,P] = luinc(A,opt);
Warning: Incomplete upper triangular factor had 32 zero diagonals replaced
by local drop tolerance.
It may not make a good preconditioner for an iterative method.
nnz(L+U)
ans =
3185
[x,flag,relres,iter,resvec] = gmres(P*A,P*b,[],tol,479,L,U);
iter
iter =
1 70
relres
relres =
3.9947e-005
figure(3)
semilogy(resvec / resvec(1),'c')
nnz(A)
ans =
1910
[l,u,p] = lu(A);
nnz(l+u)
ans =
16186
help gmres
GMRES Generalized Minimum Residual Method.
X = GMRES(A,B) attempts to solve the system of linear equations A*X = B
for X. The N-by-N coefficient matrix A must be square and the right
hand side column vector B must have length N. This uses the unrestarted
method with MIN(N,10) total iterations.
X = GMRES(AFUN,B) accepts a function handle AFUN instead of the matrix
A. AFUN(X) accepts a vector input X and returns the matrix-vector
product A*X. In all of the following syntaxes, you can replace A by
AFUN.
X = GMRES(A,B,RESTART) restarts the method every RESTART iterations.
If RESTART is N or [] then GMRES uses the unrestarted method as above.
X = GMRES(A,B,RESTART,TOL) specifies the tolerance of the method. If
TOL is [] then GMRES uses the default, 1e-6.
X = GMRES(A,B,RESTART,TOL,MAXIT) specifies the maximum number of outer
iterations. Note: the total number of iterations is RESTART*MAXIT. If
MAXIT is [] then GMRES uses the default, MIN(N/RESTART,10). If RESTART
is N or [] then the total number of iterations is MAXIT.
X = GMRES(A,B,RESTART,TOL,MAXIT,M) and
X = GMRES(A,B,RESTART,TOL,MAXIT,M1,M2) use preconditioner M or M=M1*M2
and effectively solve the system inv(M)*A*X = inv(M)*B for X. If M is
[] then a preconditioner is not applied. M may be a function handle
returning M\X.
X = GMRES(A,B,RESTART,TOL,MAXIT,M1,M2,X0) specifies the first initial
guess. If X0 is [] then GMRES uses the default, an all zero vector.
[X,FLAG] = GMRES(A,B,...) also returns a convergence FLAG:
0 GMRES converged to the desired tolerance TOL within MAXIT iterations.
1 GMRES iterated MAXIT times but did not converge.
2 preconditioner M was ill-conditioned.
3 GMRES stagnated (two consecutive iterates were the same).
[X,FLAG,RELRES] = GMRES(A,B,...) also returns the relative residual
NORM(B-A*X)/NORM(B). If FLAG is 0, then RELRES <= TOL. Note with
preconditioners M1,M2, the residual is NORM(M2\(M1\(B-A*X))).
[X,FLAG,RELRES,ITER] = GMRES(A,B,...) also returns both the outer and
inner iteration numbers at which X was computed: 0 <= ITER(1) <= MAXIT
and 0 <= ITER(2) <= RESTART.
[X,FLAG,RELRES,ITER,RESVEC] = GMRES(A,B,...) also returns a vector of
the residual norms at each inner iteration, including NORM(B-A*X0).
Note with preconditioners M1,M2, the residual is NORM(M2\(M1\(B-A*X))).
Example:
n = 21; A = gallery('wilk',n); b = sum(A,2);
tol = 1e-12; maxit = 15; M = diag([10:-1:1 1 1:10]);
x = gmres(A,b,10,tol,maxit,M);
Or, use this matrix-vector product function
%-----------------------------------------------------------------%
function y = afun(x,n)
y = [0; x(1:n-1)] + [((n-1)/2:-1:0)'; (1:(n-1)/2)'].*x+[x(2:n); 0];
%-----------------------------------------------------------------%
and this preconditioner backsolve function
%------------------------------------------%
function y = mfun(r,n)
y = r ./ [((n-1)/2:-1:1)'; 1; (1:(n-1)/2)'];
%------------------------------------------%
as inputs to GMRES:
x1 = gmres(@(x)afun(x,n),b,10,tol,maxit,@(x)mfun(x,n));
Class support for inputs A,B,M1,M2,X0 and the output of AFUN:
float: double
See also bicg, bicgstab, cgs, lsqr, minres, pcg, qmr, symmlq, luinc,
function_handle.
Reference page in Help browser
doc gmres
close all
quit