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

classic Classic list List threaded Threaded
27 messages Options
12
Reply | Threaded
Open this post in threaded view
|

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

Joza
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
Reply | Threaded
Open this post in threaded view
|

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

Joza
I should add, A only has elements on the main diagonal, and the next upper and lower diagonals. All other entries are zero.
Reply | Threaded
Open this post in threaded view
|

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

Sergei Steshenko
In reply to this post by Joza





----- Original Message -----

> From: Joza <[hidden email]>
> To: [hidden email]
> 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.
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

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

Sergei Steshenko




----- Original Message -----

> From: Sergei Steshenko <[hidden email]>
> To: Joza <[hidden email]>; "[hidden email]" <[hidden email]>
> Cc:
> Sent: Saturday, November 3, 2012 12:44 AM
> Subject: Re: Successive Over-Relaxation ... what is wrong? Improvements?
>
>
>
>
>
>
> ----- Original Message -----
>>  From: Joza <[hidden email]>
>>  To: [hidden email]
>>  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.


And, by the way, this is completely wrong to debug such problems on data having several thousand elements - because one will quickly get lost in this number of values.

Unless he/she is a genius making no mistakes.

So, if I understand correctly, you should first reduce your problem to, say, 4 x 4 sparse matrix and using debugger and/or print statements debug it.

Regards,
  Sergei.


>
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

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

Joza
In reply to this post by Joza
Sorry, I forgot I shouldn't have such numbers there ... I'll replace the likes of j=999:1000 with j = n-1:n.

The code works fine for a 10x10 matrix that I tried. But Octave won't churn out an answer with 1000x1000. I'm quite new to sparse matrices so I'm wondering if there are ways to speed up the code or optimize it, perhaps by taking advantage of the sparse matrix's properties?
Reply | Threaded
Open this post in threaded view
|

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

martin_helm
What about just adding a print statement which shows the value of iter
and norm_delta_x and

tol + 4.0*(10^-16)*norm_x

to see if you ever come near to your convergence criteria?

_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

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

Joza
Well, here's the read-out for the last iteration:

iters: 20374    norm: 2.0769230769230695

yet

iters: 100      norm: 2.0769230769230771
Reply | Threaded
Open this post in threaded view
|

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

Sergei Steshenko
In reply to this post by Joza




----- Original Message -----

> From: Joza <[hidden email]>
> To: [hidden email]
> Cc:
> Sent: Saturday, November 3, 2012 1:34 AM
> Subject: Re: Successive Over-Relaxation ... what is wrong? Improvements?
>
> Sorry, I forgot I shouldn't have such numbers there ... I'll replace the
> likes of j=999:1000 with j = n-1:n.
>
> The code works fine for a 10x10 matrix that I tried. But Octave won't churn
> out an answer with 1000x1000. I'm quite new to sparse matrices so I'm
> wondering if there are ways to speed up the code or optimize it, perhaps by
> taking advantage of the sparse matrix's properties?
>
>
>
> --
> View this message in context:
> http://octave.1599824.n4.nabble.com/Successive-Over-Relaxation-what-is-wrong-Improvements-tp4646041p4646046.html
> Sent from the Octave - General mailing list archive at Nabble.com.


Insert some diagnostic print statements allowing you to monitor progress.

Depending on how you count, you may say that 100 x 100 in the middle of 10 x 10 and 1000 x 1000 piece.

So, try 100 x 100 first.

Regards,
  Sergei.

_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

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

eem2314
In reply to this post by Joza


On Fri, Nov 2, 2012 at 5:13 PM, Joza <[hidden email]> wrote:
Well, here's the read-out for the last iteration:

iters: 20374    norm: 2.0769230769230695

yet

iters: 100      norm: 2.0769230769230771




--
View this message in context: http://octave.1599824.n4.nabble.com/Successive-Over-Relaxation-what-is-wrong-Improvements-tp4646041p4646048.html
Sent from the Octave - General mailing list archive at Nabble.com.
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave

It looks like it is oscillating about the solution, which probably means you need to fiddle with
the relaxation parameter

--
Ed Meyer


_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

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

Sergei Steshenko
In reply to this post by Joza

--- On Fri, 11/2/12, Ed Meyer <[hidden email]> wrote:

From: Ed Meyer <[hidden email]>
Subject: Re: Successive Over-Relaxation ... what is wrong? Improvements?
To: "Joza" <[hidden email]>
Cc: [hidden email]
Date: Friday, November 2, 2012, 8:37 PM



On Fri, Nov 2, 2012 at 5:13 PM, Joza <[hidden email]> wrote:


Well, here's the read-out for the last iteration:



iters: 20374    norm: 2.0769230769230695



yet



iters: 100      norm: 2.0769230769230771









--

View this message in context: http://octave.1599824.n4.nabble.com/Successive-Over-Relaxation-what-is-wrong-Improvements-tp4646041p4646048.html



Sent from the Octave - General mailing list archive at Nabble.com.

_______________________________________________

Help-octave mailing list

[hidden email]

https://mailman.cae.wisc.edu/listinfo/help-octave


It looks like it is oscillating about the solution, which probably means you need to fiddle with
the relaxation parameter

--
Ed Meyer



-----Inline Attachment Follows-----

_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave


IIRC, the OP says Octave built-in algorithm converges, so the problem is most likely an error in the OP's implementation.

Regards,
  Sergei.

_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

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

Joza
Yea....

I've tried a few values for the relaxation parameter and the matrix size ...

it seems that either the norm converges to a constant, or it keeps oscillating around a certain value.

As far as I am aware, I've followed the implementation as per textbook - I can't see anything obviously wrong with it. Any idea what could be causing this?
Reply | Threaded
Open this post in threaded view
|

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

Sergei Steshenko





----- Original Message -----

> From: Joza <[hidden email]>
> To: [hidden email]
> Cc:
> Sent: Sunday, November 4, 2012 3:12 AM
> Subject: Re: Successive Over-Relaxation ... what is wrong? Improvements?
>
> Yea....
>
> I've tried a few values for the relaxation parameter and the matrix size ...
>
> it seems that either the norm converges to a constant, or it keeps
> oscillating around a certain value.
>
> As far as I am aware, I've followed the implementation as per textbook - I
> can't see anything obviously wrong with it. Any idea what could be causing
> this?
>
>
>

Do you have any test case that converges ?

--Sergei.

_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

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

Joza
Yea, for A = 10x10, for example, here is the result, with relaxation parameter = 1.4:

...
 x: 4.999998     iters: 710      norm: 0.0000010676089660
 x: 4.999998     iters: 711      norm: 0.0000003683371855
 x: 4.999998     iters: 712      norm: 0.0000006345346101
 x: 4.999998     iters: 713      norm: 0.0000007962736616
 x: 4.999998     iters: 714      norm: 0.0000008603756717
 x: 4.999998     iters: 715      norm: 0.0000008603756717
 x: 4.999998     iters: 716      norm: 0.0000008603756717
 x: 4.999998     iters: 717      norm: 0.0000008603756717
 x: 4.999998     iters: 718      norm: 0.0000008603756717
 x: 4.999998     iters: 719      norm: 0.0000008603756717
 x: 4.999998     iters: 720      norm: 0.0000008603756717

ans =

    5.0000
    9.0000
   12.0000
   14.0000
   15.0000
   15.0000
   14.0000
   12.0000
    9.0000
    5.0000

So it seems to work here?

Reply | Threaded
Open this post in threaded view
|

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

Joza
Here is my textbook's implementation of the Algorithm which I followed:
The sparse matrix structure a contains val, col, row, d ie. 4 arrays containing non-zero values, column indices of those values, where each row starts, and the diagonal values. (I haven't done this part exactly as I figured it wouldn't make a difference since Octave stores the matrix I created anyway as a sparse matrix...right?...)

for k := 1 to maxits do
      normx := 0; dnormx := 0
      for i := 1 to n do
            sum : = 0
         for j :=rowstart[i ] to rowstart[i + 1] − 1 do
            sum := sum + val [ j ] × x [ col [ j ] ]
         endfor j
         delx := ω × (b[i ] − sum)/d[i ]
         x [i ] := x [i ] + delx
         normx := max(normx, abs( x [i ]) )
         dnormx := max(dnormx, abs(delx ) )
      endfor i
      if dnormx ≤ tol + 4 × m × normx then return (x, k)
endfor k

Reply | Threaded
Open this post in threaded view
|

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

c.-2
In reply to this post by Joza

On 4 Nov 2012, at 14:05, Joza wrote:

> Yea, for A = 10x10, for example, here is the result, with relaxation
> parameter = 1.4:
>
> ...
> x: 4.999998     iters: 710      norm: 0.0000010676089660
> x: 4.999998     iters: 711      norm: 0.0000003683371855
> x: 4.999998     iters: 712      norm: 0.0000006345346101
> x: 4.999998     iters: 713      norm: 0.0000007962736616
> x: 4.999998     iters: 714      norm: 0.0000008603756717
> x: 4.999998     iters: 715      norm: 0.0000008603756717
> x: 4.999998     iters: 716      norm: 0.0000008603756717
> x: 4.999998     iters: 717      norm: 0.0000008603756717
> x: 4.999998     iters: 718      norm: 0.0000008603756717
> x: 4.999998     iters: 719      norm: 0.0000008603756717
> x: 4.999998     iters: 720      norm: 0.0000008603756717
>
> ans =
>
>    5.0000
>    9.0000
>   12.0000
>   14.0000
>   15.0000
>   15.0000
>   14.0000
>   12.0000
>    9.0000
>    5.0000
>
> So it seems to work here?
>

not really, it actually seems to be stagnating.

I just tried a very simple implementation of the algorithm wiyh your parameters
and it converges within a tolerance of 1e-15 in about 170 iterations.

c.







 

_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

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

Joza
Would you mind posting the algorithm code?

I'd like to compare it to mine and see exactly what is going on.
Reply | Threaded
Open this post in threaded view
|

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

Sergei Steshenko





----- Original Message -----

> From: Joza <[hidden email]>
> To: [hidden email]
> Cc:
> Sent: Sunday, November 4, 2012 8:29 PM
> Subject: Re: Successive Over-Relaxation ... what is wrong? Improvements?
>
> Would you mind posting the algorithm code?
>
> I'd like to compare it to mine and see exactly what is going on.
>
>
>
> --
> View this message in context:
> http://octave.1599824.n4.nabble.com/Successive-Over-Relaxation-what-is-wrong-Improvements-tp4646041p4646075.html
> Sent from the Octave - General mailing list archive at Nabble.com.
> _______________________________________________
> Help-octave mailing list
> [hidden email]
> https://mailman.cae.wisc.edu/listinfo/help-octave
>

The code you've posted as an example in written in Pascal-like language (if not in Pascal itself) and thus is trivially translatable into any other procedural language.

So, show some diligence and perseverance.

I.e. translate the code, take a simple case (4 x 4 matrix), do the calculations "manually" (say, using a calculator), compare the results to what's happening in your program.

This all is called debugging, and is not Octave-specific.

Regards,
  Sergei.

_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

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

Joza
I understand that much very well, but the code I posted in the OP was written Octave. I have already explained that the code works fine for a 10 by 10 sparse matrix, but not for a 1000 by 1000 sparse matrix.

Hence, I have been asking for help with my Octave code, from individuals here, whom, for the most part, are far more experience with the language than I am, and hopefully could shed some light on implementation problems with my code written in Octave.

Reply | Threaded
Open this post in threaded view
|

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

eem2314
In reply to this post by Joza


On Sat, Nov 3, 2012 at 6:12 PM, Joza <[hidden email]> wrote:
Yea....

I've tried a few values for the relaxation parameter and the matrix size ...

it seems that either the norm converges to a constant, or it keeps
oscillating around a certain value.

As far as I am aware, I've followed the implementation as per textbook - I
can't see anything obviously wrong with it. Any idea what could be causing
this?



--
View this message in context: http://octave.1599824.n4.nabble.com/Successive-Over-Relaxation-what-is-wrong-Improvements-tp4646041p4646059.html
Sent from the Octave - General mailing list archive at Nabble.com.
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave

It works for me if I use the formula for the overrelaxation factor
  w = 2/(1 + pi/n)
that is usually recommended

--
Ed Meyer


_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

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

eem2314


On Sun, Nov 4, 2012 at 11:31 AM, Ed Meyer <[hidden email]> wrote:


On Sat, Nov 3, 2012 at 6:12 PM, Joza <[hidden email]> wrote:
Yea....

I've tried a few values for the relaxation parameter and the matrix size ...

it seems that either the norm converges to a constant, or it keeps
oscillating around a certain value.

As far as I am aware, I've followed the implementation as per textbook - I
can't see anything obviously wrong with it. Any idea what could be causing
this?



--
View this message in context: http://octave.1599824.n4.nabble.com/Successive-Over-Relaxation-what-is-wrong-Improvements-tp4646041p4646059.html
Sent from the Octave - General mailing list archive at Nabble.com.
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave

It works for me if I use the formula for the overrelaxation factor
  w = 2/(1 + pi/n)
that is usually recommended

--
Ed Meyer

I would also point out the the example you are using is not diagonally dominant,
which is usually a requirement for convergence


--
Ed Meyer


_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
12