Numerical Recipes Forum  

Go Back   Numerical Recipes Forum > Obsolete Editions Forum > Methods: Chapters 7, 8, and 20

Reply
 
Thread Tools Display Modes
  #1  
Old 10-13-2005, 07:04 AM
p.minguzzi p.minguzzi is offline
Member
 
Join Date: Oct 2005
Location: Italy
Posts: 2
Question Shell sorting

I found a problem with the Shell subroutine of Chapter 8.1
(Fortran 77): in the rare case when the dimension N of the array A is exactly (3^k-1)/2 the algorithm accesses the location A(N+1), which in principle is undefined and can contain any value. If it happens that this alien value is smaller than all the elements of the array it will be inserted at the first place in the sorted array.
These statements can be checked by the following simple example. Let us assume that we have a 5 element set A={1,5,2,7,0} and we want to sort the subset {1,5,2,7} of the first 4 elements of A, then a call Shell(4,A) will produce this output: {0,2,5,7}.
For just this case a correct answer is otherwise obtained if in the first IF statement .LE. is changed to .LT. ; however it is not clear to me whether this is the correct fix for every situation.

P. Minguzzi
Reply With Quote
  #2  
Old 10-13-2005, 01:18 PM
Saul Teukolsky Saul Teukolsky is offline
Numerical Recipes Author
 
Join Date: Jan 2002
Posts: 212
Hi,

I am unable to reproduce your result. The array is correctly sorted by shell.

Best wishes,
Saul Teukolsky
Reply With Quote
  #3  
Old 10-17-2005, 10:49 AM
p.minguzzi p.minguzzi is offline
Member
 
Join Date: Oct 2005
Location: Italy
Posts: 2
The problem that I reported was not caused by a trivial misuse of the SHELL routine, so after the short answer of the Author, I tried to get a deeper understanding of what was happening.
There can be little doubt that, given N = 4, after the

INC = 1
1 INC = 3*INC+1
if (INC.le.N) goto 1
INC = INC/3

sequence of statements, INC has the value of 4, so at the start of the DO loop the counter I = INC+1 = 5.
Unfortunately I am not a beginner and in the old days of Fortran DO loops were always executed at least once, since the test of completion was made at the end of the loop (incidentally this is also the low-level operation of most modern processors, e.g. the LOOP instruction of x86). So I was not particularly surprised of getting a bad result, since A(5)=0. However many recent Fortran compilers will not enter the DO loop if the iteration count, as established from the parameter list, is zero. In some cases the compiler has an option to select whether at least one execution is made or not: usually the default is NOT.
In my case the default was disabled and I did not realize it, because the behaviour was quite logical for me.
What will the general user learn from my post ? In some rare cases the correct operation of the SHELL routine depends critically on the default option of the compiler: so be careful.
Reply With Quote
  #4  
Old 10-17-2005, 11:00 AM
Saul Teukolsky Saul Teukolsky is offline
Numerical Recipes Author
 
Join Date: Jan 2002
Posts: 212
Thanks for tracking your problem down. In our defense, the "always execute a do loop once" option was the standard in f66. Since f77, the default has been the opposite. 28 years is a long time!

Best wishes,
Saul Teukolsky
Reply With Quote
Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is Off
HTML code is Off

Forum Jump


All times are GMT -5. The time now is 03:16 AM.


Powered by vBulletin® Version 3.8.6
Copyright ©2000 - 2014, Jelsoft Enterprises Ltd.