Velocity Reviews > Error in Numerical Recipes lubksb routine

# Error in Numerical Recipes lubksb routine

mma
Guest
Posts: n/a

 12-07-2003
I have been using the lubksb routine in Visual C++ 6.0 and noticed
what looks like an error to me. The last section of the method looks
like this:

for(i=n;i>=1;i--)
{
sum=b[i];
for(j=i+1;j<=n;j++)
sum -= a[i][j]*b[j];
b[i]=sum/a[i][i];
}

Where a is an n X n matrix (C array of arrays) which go from 1..n
(through the use of Numerical Recipes offset index method). See page
48 here:
<a href="http://www.library.cornell.edu/nr/bookcpdf/c2-3.pdf">Numerical

It appears to me that the inner loop will cause a segmentation fault
because j=i+1 = n+1 which is greater than n--the max index for both a
and b. Yet it seems to work--am I missing something?

--Mike

Chris Torek
Guest
Posts: n/a

 12-08-2003
In article <news:(E-Mail Removed). com>
mma <(E-Mail Removed)> writes:
>I have been using the lubksb routine in Visual C++ 6.0 and noticed
>what looks like an error to me. The last section of the method looks
>like this:
>
>for(i=n;i>=1;i--)
>{
> sum=b[i];
> for(j=i+1;j<=n;j++)
> sum -= a[i][j]*b[j];
> b[i]=sum/a[i][i];
>}
>
>Where a is an n X n matrix (C array of arrays) which go from 1..n
>(through the use of Numerical Recipes offset index method).

Be aware that this method is not strictly portable, although it
works in practice on most if not all systems. For more detail,
see the FAQ.

>It appears to me that the inner loop will cause a segmentation fault
>because j=i+1 = n+1 which is greater than n--the max index for both a
>and b. Yet it seems to work--am I missing something?

When i==n, the inner loop does this:

for (j = i + 1;

- This sets j to i + 1, which is also n + 1, as you note.

j <= n;

- Now compare j to n. Since j = n + 1, j > n, and j <= n is false
(produces the integer value 0 in an expression, and is treated
as false in a loop-controlling-expression). The inner loop
therefore does not run at all.

After the inner loop (immediately) terminates, the outer loop sets
b[i] to sum / a[i][i]. The variable sum is unmodified from the
original b[i] assigned to it, so this divides b[i] by a[i][i].
The outer loop then decrements i and resumes at its test expression
(which must have already been true exactly once).
--
In-Real-Life: Chris Torek, Wind River Systems
Salt Lake City, UT, USA (40°39.22'N, 111°50.29'W) +1 801 277 2603
Reading email is like searching for food in the garbage, thanks to spammers.

Tim Prince
Guest
Posts: n/a

 12-08-2003
mma wrote:

> I have been using the lubksb routine in Visual C++ 6.0 and noticed
> what looks like an error to me. The last section of the method looks
> like this:
>
> for(i=n;i>=1;i--)
> {
> sum=b[i];
> for(j=i+1;j<=n;j++)
> sum -= a[i][j]*b[j];
> b[i]=sum/a[i][i];
> }
>
> Where a is an n X n matrix (C array of arrays) which go from 1..n
> (through the use of Numerical Recipes offset index method). See page
> 48 here:
> <a href="http://www.library.cornell.edu/nr/bookcpdf/c2-3.pdf">Numerical
>
> It appears to me that the inner loop will cause a segmentation fault
> because j=i+1 = n+1 which is greater than n--the max index for both a
> and b. Yet it seems to work--am I missing something?
>

The condition j<=n causes the inner loop to end without doing anything for
the first value of i. While it might be slightly more efficient on some
platform to write a separate case for j=i+1, with no inner loop, the way it
is written is the most succinct. If your compiler fetches the operands
before testing whether they are needed, and this causes a segfault, that's
--
Tim Prince

Barry Schwarz
Guest
Posts: n/a

 12-08-2003
On 7 Dec 2003 13:57:55 -0800, http://www.velocityreviews.com/forums/(E-Mail Removed) (mma) wrote:

>I have been using the lubksb routine in Visual C++ 6.0 and noticed
>what looks like an error to me. The last section of the method looks
>like this:
>
>for(i=n;i>=1;i--)
>{
> sum=b[i];
> for(j=i+1;j<=n;j++)
> sum -= a[i][j]*b[j];
> b[i]=sum/a[i][i];
>}
>
>Where a is an n X n matrix (C array of arrays) which go from 1..n
>(through the use of Numerical Recipes offset index method). See page
>48 here:
><a href="http://www.library.cornell.edu/nr/bookcpdf/c2-3.pdf">Numerical
>
>It appears to me that the inner loop will cause a segmentation fault
>because j=i+1 = n+1 which is greater than n--the max index for both a
>and b. Yet it seems to work--am I missing something?

When i is set to n, then j=i+1 will indeed set j to n+1. However, the
terminating condition (j<=n) is examined at the top of the loop so the
loop will not execute at all.

If it did execute, it would cause undefined behavior. A segmentation
fault is one of the more user friendly manifestations of undefined
behavior but not guaranteed.

<<Remove the del for email>>

glen herrmannsfeldt
Guest
Posts: n/a

 12-09-2003
mma wrote:

> I have been using the lubksb routine in Visual C++ 6.0 and noticed
> what looks like an error to me. The last section of the method looks
> like this:
>
> for(i=n;i>=1;i--)
> {
> sum=b[i];
> for(j=i+1;j<=n;j++)
> sum -= a[i][j]*b[j];
> b[i]=sum/a[i][i];

(snip)

> It appears to me that the inner loop will cause a segmentation fault
> because j=i+1 = n+1 which is greater than n--the max index for both a
> and b. Yet it seems to work--am I missing something?

The j<=n; right after the j=i+1;

-- glen

Aleksander Nabaglo
Guest
Posts: n/a

 12-09-2003
!
On Tue, 09 Dec 2003 05:17:58 +0000, glen herrmannsfeldt wrote:

>> for(i=n;i>=1;i--)
>> {
>> sum=b[i];
>> for(j=i+1;j<=n;j++)
>> sum -= a[i][j]*b[j];
>> b[i]=sum/a[i][i];

>
> (snip)
>
>> It appears to me that the inner loop will cause a segmentation fault
>> because j=i+1 = n+1 which is greater than n--the max index for both a
>> and b. Yet it seems to work--am I missing something?

>
> The j<=n; right after the j=i+1;

for(i1; i2; i3) S;

is a shortcut for:

i1;
while(i2) { S; i3; }

--
A
..