转贴 sor method， resource code
<P>function = 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>% = 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>% outputx 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 < 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 <= tol ), break, end % check convergence</P>
<P>end<BR>b = b / w; % restore rhs</P>
<P>if ( error > 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>% outputM DOUBLE PRECISION matrix<BR>% N DOUBLE PRECISION matrix such that A = M - N<BR>% b DOUBLE PRECISION rhs vector ( altered for SOR )</P>
<P> = 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>
回复：（xj2070）转贴 sor method， resource co...
这个最好转到算法版或者matlab版吧
页:
[1]