🚀 go-pugleaf

iwelch@agsm.ucla.edu (Ivo Welch) writes:


*I  have  the 1989 version  of   Numerical Recipees  in  C, both   in book  and
*electronic form. I believe I have found a bug---or at least  a nuisance.  When
*a matrix  contains some  extreme values,  it   is  possible  to get  an access
*violation.

[ Code deleted. ]

*# note that l ends with a value of 0 if the loop exits. Typically,
*# rv1[l]=0, so the loop is exited early---no problem then. However,
*# if the matrix is poorly setup (and I have had such a case), rv1[1]
*# can become NaN (0/0), and the loop completes rather than breaks.
*# I believe rv1[1] turned into NaN, since my matrix contained an IEEE
*# -Inf.

I have had the same problem, but *only* when my input matrix contained
an Inf or NaN as appears to be the case in your example.  So, unless
you intend to work with infinities, it would appear that the code
works.  However, the fact that the program produces a segmentation
violation and core dumps as the result of a data error *is* annoying.

*# I presume this problem can be fixed if the condition reads l>=2,
*# and nm=l-1 is moved into both the start and iteration.

Perhaps you should just check for l == 0 at the end of the loop and
exit gracefully with an error message.  If you disallow l = 1 in the
body of the loop, the code may not work properly.  I think there is a
comment in the book (if not in the code) that rv1[1] is always zero so
that the loop will always exit before l = 0.  If rv1[1] is not zero,
there must be something wrong with the data.  Therefore, you could
check the value of rv1[1] to see if it is correct at various points in
the code.

Darrell

P.S. I would also be interested in seeing a bug list for Numerical
Recipes in C is anyone has it on line.

I have just found (perhaps) a more serious error in the same code. Look at nm:
it is  l-1 when the first  loop exits. So, the  loop exits with  l=1 and nm=0.
Look at the next use of nm. It is an index into the a[] array. ARGGGH!

Furthermore, perhaps someone can enlighten me as to what the statement

        if (fabs(rv1[l])+anorm == anorm)

tests. anorm is the max of abs(w[]) + abs(rv1[]). What ugly code.

/ivo

In article <1992Aug5.104832.27184@mic.ucla.edu> iwelch@agsm.ucla.edu (Ivo Welch) writes:
>
>I  have  the 1989 version  of   Numerical Recipees  in  C, both   in book  and
>electronic form. I believe I have found a bug---or at least  a nuisance.  When
>a matrix  contains some  extreme values,  it   is  possible  to get  an access
>violation.
>
>(Does anyone have a bug list of NRinC?)
>
>Note: The code below is copyrighted, but since it is only a small  fragment, I
>presume it is ok to post it to better explain the bug.
>
>----------------
>
>Description: On Chapter 2, page 70, the code looks like
>
>  for (k=n;k>=1;k--) {
>    for (its=1;its<=30;its++) {
>      flag=1;
>      for (l=k;l>=1;l--) {
>        nm=l-1;
>        if (fabs(rv1[l])+anorm == anorm) {
>          flag=0;
>          break;
>        }
>        if (fabs(w[nm])+anorm == anorm) break;
>      }
>
># note that l ends with a value of 0 if the loop exits. Typically,
># rv1[l]=0, so the loop is exited early---no problem then. However,
># if the matrix is poorly setup (and I have had such a case), rv1[1]
># can become NaN (0/0), and the loop completes rather than breaks.
># I believe rv1[1] turned into NaN, since my matrix contained an IEEE
># -Inf.
>
># I presume this problem can be fixed if the condition reads l>=2,
># and nm=l-1 is moved into both the start and iteration.
>
>      if (flag) {
>        c=0.0;
>        s=1.0;
>        for (i=l;i<=k;i++) {
>
># note that i starts with a value of l. If l can be 0
># this means that the subsequent indices refer to 0, not 1.
># Since NRinC allocates all matrices/vectors from 1, not 0,
># this results in an access violation error---if you are lucky.
>
>          f=s*rv1[i];
>          if (fabs(f)+anorm != anorm) {
>            g=w[i];
>            h=PYTHAG(f,g);
>            w[i]=h;
>            h=1.0/h;
>            c=g*h;
>            s=(-f*h);
>            for (j=1;j<=m;j++) {
>              y=a[j][nm];
>              z=a[j][i];
>              a[j][nm]=y*c+z*s;
>              a[j][i]=z*c-y*s;
>            }
>          }
>        }
>      }

In article <1992Aug6.085315.10239@mic.ucla.edu> iwelch@agsm.ucla.edu
(Ivo Welch) writes:
> [ referring to the "Numerical Recipes in C" SVDCMP function ]
>Furthermore, perhaps someone can enlighten me as to what the statement
>
>        if (fabs(rv1[l])+anorm == anorm)
>
>tests. anorm is the max of abs(w[]) + abs(rv1[]).

This is actually a nice way of testing whether or not  rv1[l]  is
"negligably small" relative to  w[] .  That is, in general the test

	if (x == (x + e))

(the  fabs()  actually isn't needed) tests whether or not  e  is
negligably small relative to  x .  More precisely, to within at
most a factor of two (which doesn't matter here) it tests for
 |e/x| < eps  where  eps  is the "machine epsilon", i.e. the
smallest positive floating point number such that  1.0+eps != 1.0
in floating point.  (eps = 5.97e-8 for IEEE single, 1.12e-16 for
IEEE double precision.)  This can be useful in testing whether or
not some iteration which updates  x  by adding  e  to it, has
converged.

This is a *nice* way of testing this because it *automatically*
adapts to the floating point format used.  That is, code written
using this style of test should work unchanged with any reasonable
floating point system (even Cray) *without* requiring changes to
"machine constants".

The same trick can also be used to test for things being "almost
negligable":

	if (x == (x + 0.01*e))

tests for  x + e  and  x  differing by at most 100 "units in the
last place" (100 times the value of the least significant bit in
 x ), i.e. for  |e/x| < 100*eps .  Here again the error tolerance
*automatically* adapts to the floating point precision.

(The one caveat in this is that the usual floating point precision
must be used for the computation/storage of  (x + e) , or at least
this quantity must be stored in a usual-precision register or memory
location.  If this computation is done/stored in (say) an extended-
precision register, all bets are off.  The usual hack to force
 (x + e)  to be computed in the "usual precision" is (using Fortran)

		double precision x, e
		double precision temp(2)
		common /extemp/ temp
		...
	c
	c	is  e  negligably small compared to  x ?
		temp(1) = x
		temp(2) = temp(1) + e
		if (temp(1) .eq. temp(2)) then
			... e  is negligably small compared to  x
		else
			... it isn't
		end if

(In C the array would be declared with "extern" storage class, i.e.
it would be declared at "top level", not inside a function, and would
*not* be declared "static".)

Here the use of an array, and its placement in a common block, are
intended to force the compiler to actually compute and store  x + e
in memory, hopefully rounding it to the "usual precision".  This is
rather ugly, but is safe against most optimizing compilers and
extended-precision-registers hardware.

- Jonathan Thornburg
  <jonathan@einstein.ph.utexas.edu> or <jonathan@hermes.chpc.utexas.edu>
  University of Texas at Austin / Physics Dept / Center for Relativity
  and (for a few more months) U of British Columbia / {Astronomy,Physics}