445,819 Members | 1,185 Online Need help? Post your question and get tips & solutions from a community of 445,819 IT Pros & Developers. It's quick & easy.

# Error in Numerical Recipes lubksb routine

 P: n/a 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: Numerical Recipes sample pdf 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 Nov 13 '05 #1
5 Replies

 P: n/a 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: Numerical Recipes sample pdf 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 Nov 13 '05 #2

 P: n/a In article mma writes:I have been using the lubksb routine in Visual C++ 6.0 and noticedwhat looks like an error to me. The last section of the method lookslike 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 faultbecause j=i+1 = n+1 which is greater than n--the max index for both aand 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 email: forget about it http://web.torek.net/torek/index.html Reading email is like searching for food in the garbage, thanks to spammers. Nov 13 '05 #3

 P: n/a On 7 Dec 2003 13:57:55 -0800, mm****@foobox.net (mma) wrote: I have been using the lubksb routine in Visual C++ 6.0 and noticedwhat looks like an error to me. The last section of the method lookslike 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 page48 here:NumericalRecipes sample pdf for more information.It appears to me that the inner loop will cause a segmentation faultbecause j=i+1 = n+1 which is greater than n--the max index for both aand 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. <> Nov 13 '05 #4

 P: n/a 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 Nov 14 '05 #5

 P: n/a ! 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 .. Nov 14 '05 #6

### This discussion thread is closed

Replies have been disabled for this discussion. 