[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Successive Over-Relaxation ... what is wrong? Improvements?
From: |
Sergei Steshenko |
Subject: |
Re: Successive Over-Relaxation ... what is wrong? Improvements? |
Date: |
Fri, 2 Nov 2012 15:44:53 -0700 (PDT) |
----- 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.