Error in Numerical Recipes lubksb routine

Discussion in 'C Programming' started by mma, Dec 7, 2003.

  1. mma

    mma Guest

    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;
    for(j=i+1;j<=n;j++)
    sum -= a[j]*b[j];
    b=sum/a;
    }

    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
    Recipes sample pdf</a> for more information.

    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
    mma, Dec 7, 2003
    #1
    1. Advertising

  2. mma

    Chris Torek Guest

    In article <news:>
    mma <> 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;
    > for(j=i+1;j<=n;j++)
    > sum -= a[j]*b[j];
    > b=sum/a;
    >}
    >
    >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 to sum / a. The variable sum is unmodified from the
    original b assigned to it, so this divides b by a.
    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
    email: forget about it http://web.torek.net/torek/index.html
    Reading email is like searching for food in the garbage, thanks to spammers.
    Chris Torek, Dec 8, 2003
    #2
    1. Advertising

  3. mma

    Tim Prince Guest

    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;
    > for(j=i+1;j<=n;j++)
    > sum -= a[j]*b[j];
    > b=sum/a;
    > }
    >
    > 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
    > Recipes sample pdf</a> for more information.
    >
    > 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
    the fault of your compiler.
    --
    Tim Prince
    Tim Prince, Dec 8, 2003
    #3
  4. On 7 Dec 2003 13:57:55 -0800, (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;
    > for(j=i+1;j<=n;j++)
    > sum -= a[j]*b[j];
    > b=sum/a;
    >}
    >
    >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
    >Recipes sample pdf</a> for more information.
    >
    >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>>
    Barry Schwarz, Dec 8, 2003
    #4
  5. 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;
    > for(j=i+1;j<=n;j++)
    > sum -= a[j]*b[j];
    > b=sum/a;


    (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
    glen herrmannsfeldt, Dec 9, 2003
    #5
  6. !
    On Tue, 09 Dec 2003 05:17:58 +0000, glen herrmannsfeldt wrote:

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

    >
    > (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
    ..
    Aleksander Nabaglo, Dec 9, 2003
    #6
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. Al Murphy

    Java numerical recipes online...

    Al Murphy, Feb 5, 2004, in forum: Java
    Replies:
    1
    Views:
    5,621
    Steve W. Jackson
    Feb 5, 2004
  2. Rouben Rostamian

    Re: numerical recipes book

    Rouben Rostamian, Jul 30, 2003, in forum: C Programming
    Replies:
    0
    Views:
    363
    Rouben Rostamian
    Jul 30, 2003
  3. E. Robert Tisdale

    Re: numerical recipes book

    E. Robert Tisdale, Jul 31, 2003, in forum: C Programming
    Replies:
    0
    Views:
    375
    E. Robert Tisdale
    Jul 31, 2003
  4. Edward Hua
    Replies:
    5
    Views:
    683
  5. Lew
    Replies:
    1
    Views:
    1,305
    Frank Cisco
    Feb 21, 2009
Loading...

Share This Page