Hello everybody,


It appears to me that Numerical recipes has a mistake in the ALGORITHM, and in the code, of course in the Toeplitz linear equation solver (chapter 2.8.2). I have noticed it in the second edition. Of course, I do not have the third edition, but it happens that google-books has the right pages :)

The mistake (if confirmed) is still there.


The problem is that vector indeceis run out of vector bounds in eq 2.8.23 and 2.824, for example.
In the second edtion, the corresponding terms with R_{m+1} and R_{-m-1} are out of range of the vector size.
If m reaches its maximum of N-1 than the addressed
components are R_{N} and R_{-N}, but the vector index
is between -(N-1) and (N-1)!

In the third edition, the corresonding terms in 2.8.23 and 2.8.24 look like: R_{m+2} and R_{-m-2}. The difference
is due to shift in the indexing of the vectors, that is all.
Look also at the code of the function. The line which has
comments "compute numerator and denominator for G and H eqs 2.8.24 and 2.8.23". Trace what the maximum value
of n1-m1-1 and n1+m1+1 and you will see that it is out
of range of the vector size.

Well, these are my thoughts. Either I am wrong, and I would like to understand why or since for 15 years nobody noticed/complained/fixed it, the algorithm presented is not
of interest. In the latter case, could someone suggest a different algorithm which is described in details similar to
that presented in the NR book. Or maybe some code which
can be used? I just need to solve many (100) large(1000-2000 size) Toeplitz equations...

Thanks a lot for your help.

ZF




Thanks a lot,

ZF

Dear Z,

The algorithm and the code are both fine. Equations (2.8.23) and (2.8.24) are recurrence relations to determine G and H for M+1 in terms of G and H for M. Since the maximum value of M is N-1, the equations as written are valid up to M+1 = N-1, or M = N-2. Thus there are no out of range indices.

Similarly, in the code the statement
if (m1 == n1) return;
right before the statements you were concerned about ensures that the indices are never out of range.

Best wishes,
Saul Teukolsky

please forgive me, but you say that:

"Since the maximum value of M is N-1, the equations as written are valid up to M+1 = N-1, or M = N-2. Thus there are no out of range indices."

Now, if M=N-1 than how is it possible thatn M+1=N-1???
The last step is when M+1=N thus M=N-1 is the last iteration.
Or if the iterations go only upto M=N-2 but they have to go
upto M=N-1 what is done in the last step?

Sorry, I do not understand something.
Could you comment on that.

Thanks a lot

ZF

Equations (2.8.23) and (2.8.24) are to calculate quantities G and H with an index M+1 on the left-hand side. Therefore the maximum value of M that must be used in these equations is N-2, which will then give G and H with an index N-1 as required.

Saul Teukolsky

Dear Saul Teukolsky,

Thank you very much for your clarifications. I think I know
what caused the problem.

After equation 2.8.26 there is a description of the iteration process.

"At each stage M we use equations (2.8.23) and (2.8.24) to find H_{M+1}^{M+1} and G_{M+1}^{M+1}, and then equation (2.8.25) to find the other components of H^{M+1}, G^{M+1}. From there the vectors x^{M+1}
and/or z^{M+1} are easily calculated."

The problem I had stemmed the fact that this description is
incorrect, in my opinion. The order of operations at each
iteration is reversed.

First the x^{M+1} or z^{M+1} vectors are computed
using current vectors x^{M}, G^{M} and H^{M} from 2.8.19 and 2.8.16 (in that order). After that, the new values of G^{M+1} and H^{M+1} are computed using 2.8.23, 2.8.24 and 2.8.25 from G^{M} and H^{M}.

So, it appears to me that the code is correct, the equations are correct, but the description of the usage of these equations needs to be clarified.

Thanks a lot for your help,

ZF