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

Hi,

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

Best wishes,
Saul Teukolsky

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.

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