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