On Sep 25, 4:57 pm, "Dr Ivan D. Reid"
brunel.ac.uk> wrote:> [Oh, blast! I'd just finished typing a long reply when my connexion died
> just as I was about to hit
...]
> On Tue, 25 Sep 2007 11:30:00 -0700, monir mondenet.com>
> wrote in <1190745000.383043.172...@22g2000hsm.googlegroups.com>:
>>> Ivan;
>> Based on your recent reply, I'm absolutely certain that the code
>> you're quoting from and the code I'm using are two different things
>> sharing the same name!!
>
> Well, you're using the FORTRAN edition; the .pdfs on the Web are
> Fortran 77.
>
>> a) PZEXTR(): NR, p. 568
>> your ref statements (among others):
>> do 15 k1=1, iest-1
>> and in do 16
>> qcol(j,iest) = dy(j)
>> DO NOT EXIST in my NR code !!!
>
> Yes they do, but the variable name is different:
> ...
> M1=MIN(IEST,NUSE)
> ...
> DO 15 K1=1,M1-1
> ...
> DO 16 J=1,NV
> QCOL(J,M1)=DY(J)
> 16 CONTINUE
>
>> and THE ROUTINE DOES NOT HAVE ANY SAVE statements. Recently, however,
>> I added SAVE X to retain the values of X in the sequence calls.
>
> Remember:
> a) the FORTRAN code uses implicit naming rules; the Fortran 77 declares
> variable types (but didn't use IMPLICIT NONE as that was yet to be part of
> the standard)
> b) The way FORTRAN compilers and computers of the day worked, SAVE was not
> usually necessary -- storage was static. IIRC, strictly speaking SAVE _should_
> have been used, but since it wasn't necessary people tended to leave it out.
> With the new compilers/computers where allocation may be (e.g.) on the stack
> you _do_ need to specify what to save between invocations (as a general rule
> of thumb).
>
>> I've just quickly reviewed the link you kindly provided:
>>
http://www.nrbook.com/a/bookfpdf/f16-4.pdf
>> It is very clear that the posted codes for BSSTEP, RZEXTR and PZEXTR
>> are different from the ones I'm currently using!!
>> It's also evident that the statements I initially suspected have been
>> modified to retain certain values for subsequent calls.
>
> The only modification has been the elimination of the M1 variable
> and the addition of the SAVE statements.
>
>> (Assuming these are the corrected versions, are they from the latest
>> NR edition??? I've heard that the NR-Fortran 3rd Edition is
>> available!)
>
> No, as Mike said they are from the NR in F77 book, which has been
> obsoleted (in publishing terms) by the NR in F90, but you still need the
> (updated) F77 book for the text.
>
>> So my task now is to correct the code for the 3 routines (and maybe
>> more!!) according to link and try again.
>
> If you ignore the change to do...enddo style, the differences are
> minimal. I would strongly suggest you add IMPLICIT NONE statements to
> catch any possible typos. The code in the F90 book is more extensively
> modified, making extensive use of MODULEs and ALLOCATABLEs, for example.
>
> --
> Ivan Reid, School of Engineering & Design, _____________ CMS Collaboration,
> Brunel University. Ivan.Reid@[
brunel.ac.uk|cern.ch] Room 40-1-B12, CERN
> KotPT -- "for stupidity above and beyond the call of duty".
Hello;
1. Good news! The program (~ 250 routines and over 20,000 lines of
code) now works superbly with the CORRECTED versions of BSSTEP(),
PZEXTR() and RZEXTR() for the same test case!!
At least, reasonable results, no more extrapolated values of NaN, and
no more error estimates of 1.E+36!! What a relief!
2. The polynomial function extrapolation happens to work better for
the problem at hand. I've successfully tested both extrapolation
methods.
3. The alternative to calling BSSTEP() in ODEINT() is to call the
RKQC() adaptive quality-control stepsize routine. Although the RKQC
method is generally more suited for non-smooth functions, I'm tempted
to re-try it!! I had tried RKQC sometime ago, but it failed
miserably!
Why the renewed interest now? Well, the F77 version of RKQC posted
in:
http://www.nrbook.com/a/bookfpdf/f16-2.pdf
is considerably different from the one I had tried earlier (NR, p.
558). The new version includes, among other things, a call to a new
routine rkck() using the fifth-order Cash-Karp Runge-Kutta method to
advance the solution. Further, one can't help but to assume that the
RKQC version in the NR book might also have some coding errors!
4. I suppose it's a matter of style rather than substance. Although
F90 has eliminated the need for a label on the terminating statement
of a Do loop, together with the need to remember which statements may
not be terminating statements, it's my preference to use labelled Do/
Continue statements for clarity rather than using unlabeled do and
enddo.
5. The Numerical Recipes (Fortran) Book I have is (I guess!) the 1989
Edition (ISBN 0 521 38330 7), and all its routines are in F77 (to the
best of my knowledge!). I was simply enquiring whether the NR
(Fortran) 3rd Edition; scheduled for release late 2007/early 2008, is
available now.
Thanks again for your tremendous help.
Monir
