
12

Below is an implementation I have written of the Successive OverRelaxation 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=i1: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


I should add, A only has elements on the main diagonal, and the next upper and lower diagonals. All other entries are zero.


 Original Message 
> From: Joza < [hidden email]>
> To: [hidden email]
> Cc:
> Sent: Saturday, November 3, 2012 12:25 AM
> Subject: Successive OverRelaxation ... what is wrong? Improvements?
>
> Below is an implementation I have written of the Successive OverRelaxation
> 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=i1: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.
_______________________________________________
Helpoctave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/helpoctave


 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 OverRelaxation ... what is wrong? Improvements?
>
>
>
>
>
>
>  Original Message 
>> From: Joza < [hidden email]>
>> To: [hidden email]
>> Cc:
>> Sent: Saturday, November 3, 2012 12:25 AM
>> Subject: Successive OverRelaxation ... what is wrong? Improvements?
>>
>> Below is an implementation I have written of the Successive OverRelaxation
>> 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=i1: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.
>
_______________________________________________
Helpoctave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/helpoctave


Sorry, I forgot I shouldn't have such numbers there ... I'll replace the likes of j=999:1000 with j = n1: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?


Well, here's the readout for the last iteration:
iters: 20374 norm: 2.0769230769230695
yet
iters: 100 norm: 2.0769230769230771


 Original Message 
> From: Joza < [hidden email]>
> To: [hidden email]
> Cc:
> Sent: Saturday, November 3, 2012 1:34 AM
> Subject: Re: Successive OverRelaxation ... 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 = n1: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/SuccessiveOverRelaxationwhatiswrongImprovementstp4646041p4646046.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.
_______________________________________________
Helpoctave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/helpoctave


 On Fri, 11/2/12, Ed Meyer < [hidden email]> wrote:
From: Ed Meyer < [hidden email]>
Subject: Re: Successive OverRelaxation ... 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 readout 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/SuccessiveOverRelaxationwhatiswrongImprovementstp4646041p4646048.htmlSent from the Octave  General mailing list archive at Nabble.com.
_______________________________________________
Helpoctave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/helpoctaveIt looks like it is oscillating about the solution, which probably means you need to fiddle with
the relaxation parameter

Ed Meyer
Inline Attachment Follows
_______________________________________________
Helpoctave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/helpoctaveIIRC, the OP says Octave builtin algorithm converges, so the problem is most likely an error in the OP's implementation.
Regards,
Sergei.
_______________________________________________
Helpoctave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/helpoctave


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?


 Original Message 
> From: Joza < [hidden email]>
> To: [hidden email]
> Cc:
> Sent: Sunday, November 4, 2012 3:12 AM
> Subject: Re: Successive OverRelaxation ... 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.
_______________________________________________
Helpoctave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/helpoctave


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?


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 nonzero 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


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 1e15 in about 170 iterations.
c.
_______________________________________________
Helpoctave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/helpoctave


Would you mind posting the algorithm code?
I'd like to compare it to mine and see exactly what is going on.


 Original Message 
The code you've posted as an example in written in Pascallike 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 Octavespecific.
Regards,
Sergei.
_______________________________________________
Helpoctave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/helpoctave


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.


On Sun, Nov 4, 2012 at 11:31 AM, Ed Meyer <[hidden email]> wrote:
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
_______________________________________________
Helpoctave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/helpoctave

12
