Re: Successive Over-Relaxation ... what is wrong? Improvements?

[Sergei Steshenko]

----- Original Message -----
> From: Joza <address@hidden>
> To: address@hidden
> Cc: 
> Sent: Saturday, November 3, 2012 12:25 AM
> Subject: Successive Over-Relaxation ... what is wrong? Improvements?
> 
> Below is an implementation I have written of the Successive Over-Relaxation
> Method (for pedagogical purposes). I must use a starting guess of zeros and
> b in Ax=b is all ones. A is sparse. 
> 
> I've run Octave's CGS and it solved the problem in about 1 second ... my
> code has been running for well over 5 minutes and nothing ...
> 
> I'm sure there is something I'm missing here. Any ideas? Feel free to 
> run
> the code yourself ... w is the relaxation parameter ( try 1.35) and tol is
> the convergence tolerance. 
> 
> This is the second major algorithm I've had trouble with recently. Must be
> something up with mu approach! Anyway, I'd really like to hear all
> critiques!
> ************************************************************************************************************************************************
> 
> function [x, iters] = SOR(w, tol)
> 
> r = sparse([1 1], [1 2], [2 -1], 1, 1000);
> A = toeplitz(r);
> b = ones(1000,1);
> x = zeros(1000,1);
> 
> iters = 0;
> converged = 0;
> 
> while !converged
>     
>     norm_x = 0.0;
>     norm_delta_x = 0.0;
>     for i=1:1:1000;
>         sum = 0.0;
>         
>         if i==1
>             for j=1:2
>                 sum = sum + A(i,j)*x(j);
>             end
>         
>                 elseif i==1000
>             for j=999:1000
>                 sum = sum + A(i,j)*x(j);
>             end
>         
>                 else
>             for j=i-1:i+1
>                 sum = sum + A(i,j)*x(j);
>             end
>         end
> 
>         delta_x = w*(b(i) - sum)/A(i,i);
>         x(i) = x(i) + delta_x
>         iters = iters + 1;
> 
>         norm_x = max( norm_x, abs(x(i)) );
>         norm_delta_x = max(norm_delta_x, abs(delta_x) )
>     end
>     
>     if norm_delta_x <= tol + 4.0*(10^-16)*norm_x
>         converged = 1;
>     end
> end
> 


Without going into details I can definitely tell you shouldn't be writing code 
this way.

In any decent code review lines like:

    for i=1:1:1000;
             for j=999:1000

will be ruthlessly ruled out. This is because you are using constants, though 
you most likely mean to use something which is a function of matrices sizes.

I.e. your code is error prone because you need to change it whenever matrices 
sizes change.

And, by the way, I do not know ';' in

    for i=1:1:1000;

does.


Regards,
  Sergei.