🚀 go-pugleaf

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;
>            }
>          }
>        }
>      }