Matlab code for Gauss-Seidel and Successive over relaxation iterative methods

代码matlab 维基原理

原始代码里面每迭代一次都求解线性方程组，不太合理。可以通过改写，只须求一次即可。
不过考虑到矩阵是上三角的，主要只涉及回代，也就未必不行。

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 splittingM01=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 convergencefunction [ 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

matrixSplit.m
sor.m