转贴 sor method， resource code

xj2070

<P>function [x, error, iter, flag]  = sor(A, x, b, w, max_it, tol)</P>
<P>%  -- Iterative template routine --<BR>%     Univ. of Tennessee and Oak Ridge National Laboratory<BR>%     October 1, 1993<BR>%     Details of this algorithm are described in "Templates for the<BR>%     Solution of Linear Systems: Building Blocks for Iterative<BR>%     Methods", Barrett, Berry, Chan, Demmel, Donato, Dongarra,<BR>%     Eijkhout, Pozo, Romine, and van der Vorst, SIAM Publications,<BR>%     1993. (ftp netlib2.cs.utk.edu; cd linalg; get templates.ps).<BR>%<BR>% [x, error, iter, flag]  = sor(A, x, b, w, max_it, tol)<BR>%<BR>% sor.m solves the linear system Ax=b using the <BR>% Successive Over-Relaxation Method (Gauss-Seidel method when omega = 1 ).<BR>%<BR>% input   A        REAL matrix<BR>%         x        REAL initial guess vector<BR>%         b        REAL right hand side vector<BR>%         w        REAL relaxation scalar<BR>%         max_it   INTEGER maximum number of iterations<BR>%         tol      REAL error tolerance<BR>%<BR>% output  x        REAL solution vector<BR>%         error    REAL error norm<BR>%         iter     INTEGER number of iterations performed<BR>%         flag     INTEGER: 0 = solution found to tolerance<BR>%                           1 = no convergence given max_it</P>
<P>  flag = 0;                                   % initialization<BR>  iter = 0;</P>
<P>  bnrm2 = norm( b );<BR>  if  ( bnrm2 == 0.0 ), bnrm2 = 1.0; end</P>
<P>  r = b - A*x;<BR>  error = norm( r ) / bnrm2;<BR>  if ( error &lt; tol ) return, end</P>
<P>  [ M, N, b ] = split( A, b, w, 2 );          % matrix splitting</P>
<P>  for iter = 1:max_it                         % begin iteration</P>
<P>     x_1 = x;<BR>     x   = M \ ( N*x + b );                   % update approximation</P>
<P>     error = norm( x - x_1 ) / norm( x );     % compute error<BR>     if ( error &lt;= tol ), break, end          % check convergence</P>
<P>  end<BR>  b = b / w;                                  % restore rhs</P>
<P>  if ( error &gt; tol ) flag = 1; end;           % no convergence</P>
<P>% END sor.m</P>
<P>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<BR>function [ M, N, b ] = split( A, b, w, flag )<BR>%<BR>% function [ M, N, b ] = split( A, b, w, flag )<BR>%<BR>% split.m sets up the matrix splitting for the stationary<BR>% iterative methods: jacobi and sor (gauss-seidel when w = 1.0 )<BR>%<BR>% input   A        DOUBLE PRECISION matrix<BR>%         b        DOUBLE PRECISION right hand side vector (for SOR)<BR>%         w        DOUBLE PRECISION relaxation scalar<BR>%         flag     INTEGER flag for method: 1 = jacobi<BR>%                                           2 = sor<BR>%<BR>% output  M        DOUBLE PRECISION matrix<BR>%         N        DOUBLE PRECISION matrix such that A = M - N<BR>%         b        DOUBLE PRECISION rhs vector ( altered for SOR )</P>
<P>  [m,n] = size( A );<BR>       <BR>  if ( flag == 1 ),                   % jacobi splitting</P>
<P>     M = diag(diag(A));<BR>     N = diag(diag(A)) - A;</P>
<P>  elseif ( flag == 2 ),               % sor/gauss-seidel splitting</P>
<P>     b = w * b;<BR>     M =  w * tril( A, -1 ) + diag(diag( A ));<BR>     N = -w * triu( A,  1 ) + ( 1.0 - w ) * diag(diag( A ));</P>
<P>  end;</P>
<P>% END split.m</P>

gghhjj

这个最好转到算法版或者matlab版吧