itermeth.m

来自「sareli<matlab科学计算>配套程序」· M 代码 · 共 40 行

function [x, iter]= itermeth(A,b,x0,nmax,tol,P)
%ITERMETH    General iterative method
%   X = ITERMETH(A,B,X0,NMAX,TOL,P) attempts to solve the system of linear equations A*X=B
%   for X.  The N-by-N coefficient matrix A must be not singular
%   and the right hand side column vector B must have length N.
%   If P='J' the Jacobi method is used, if P='G' the Gauss-Seidel
%   method is selected. Otherwise, P is a N-by-N matrix that play
%   the role of a preconditioner. TOL specifies the tolerance of the method.
%   NMAX specifies the maximum number of iterations.  
[n,n]=size(A); 
if nargin == 6
  if ischar(P)==1
     if P=='J' 
        L = diag(diag(A)); 
        U = eye(n); beta = 1; alpha = 1;
     elseif P == 'G' 
        L = tril(A); 
        U = eye(n); beta = 1; alpha = 1;
     end
  else
     [L,U]=lu(P); beta = 0;
  end
else
  L = eye(n); U = L; beta = 0;
end
iter = 0; 
r = b - A * x0; 
r0 = norm(r); 
err = norm (r); x = x0;
while err {\ensuremath{>}} tol & iter {\ensuremath{<}} nmax
  iter = iter + 1;    
  z = L{\ensuremath{\backslash}}r; 
  z = U{\ensuremath{\backslash}}z;   
  if beta == 0
     alpha = z'*r/(z'*A*z);  
  end
  x = x + alpha*z;    
  r = b - A * x;      
  err = norm (r) / r0;
end