原始代码里面每迭代一次都求解线性方程组,不太合理。可以通过改写,只须求一次即可。
不过考虑到矩阵是上三角的,主要只涉及回代,也就未必不行。
function [x, error, iter, flag] = sor(A, x, b, w, max_it, tol)
% -- Iterative template routine --
% Univ. of Tennessee and Oak Ridge National Laboratory
% October 1, 1993
% Details of this algorithm are described in "Templates for the
% Solution of Linear Systems: Building Blocks for Iterative
% Methods", Barrett, Berry, Chan, Demmel, Donato, Dongarra,
% Eijkhout, Pozo, Romine, and van der Vorst, SIAM Publications,
% 1993. (ftp netlib2.cs.utk.edu; cd linalg; get templates.ps).
%
% [x, error, iter, flag] = sor(A, x, b, w, max_it, tol)
%
% sor.m solves the linear system Ax=b using the
% Successive Over-Relaxation Method (Gauss-Seidel method when omega = 1 ).
%
% input A REAL matrix
% x REAL initial guess vector
% b REAL right hand side vector
% w REAL relaxation scalar
% max_it INTEGER maximum number of iterations
% tol REAL error tolerance
%
% output x REAL solution vector
% error REAL error norm
% iter INTEGER number of iterations performed
% flag INTEGER: 0 = solution found to tolerance
% 1 = no convergence given max_it
flag = 0; % initialization
iter = 0;
bnrm2 = norm( b );
if ( bnrm2 == 0.0 ), bnrm2 = 1.0; end
r = b - A*x;
error = norm( r ) / bnrm2;
if ( error < tol ) return, end
[ M, N, b ] = MatriceSplit( A, b, w, 2 ); % matrix splitting
M01=M\N;%改写
M02=M\b;
for iter = 1:max_it % begin iteration
x_1 = x;
x = M01*x+M02; %M \ ( N*x + b ); % update approximation
error = norm( x - x_1 ) / norm( x ); % compute error
if ( error <= tol ), break, end % check convergence
end
b = b / w; % restore rhs
if ( error > tol ) flag = 1; end; % no convergence
function [ M, N, b ] = MatriceSplit( A, b, w, flag )
%
% function [ M, N, b ] = MatriceSplit( A, b, w, flag )
%
% MatriceSplit.m sets up the matrix splitting for the stationary
% iterative methods: jacobi and sor (gauss-seidel when w = 1.0 )
%
% input A DOUBLE PRECISION matrix
% b DOUBLE PRECISION right hand side vector (for SOR)
% w DOUBLE PRECISION relaxation scalar
% flag INTEGER flag for method: 1 = jacobi
% 2 = sor
%
% output M DOUBLE PRECISION matrix
% N DOUBLE PRECISION matrix such that A = M - N
% b DOUBLE PRECISION rhs vector ( altered for SOR )
[m,n] = size( A );
if ( flag == 1 ), % jacobi splitting
M = diag(diag(A));
N = diag(diag(A)) - A;
elseif ( flag == 2 ), % sor/gauss-seidel splitting
b = w * b;
M = w * tril( A, -1 ) + diag(diag( A ));
N = -w * triu( A, 1 ) + ( 1.0 - w ) * diag(diag( A ));
end;
% END MatriceSplit.m
% END sor.m