Prev: Selecting the first element of a array returned by pack()
Next: size of a derived type containing pointers...
From: glen herrmannsfeldt on 30 Mar 2010 15:14
Gordon Sande <Gordon.Sande(a)eastlink.ca> wrote:
> On 2010-03-30 10:09:18 -0300, Ingo Thies <ingo.thies(a)gmx.de> said:
(snip)
>> Well, DQAG is actually a wrapper for another routine that again is a
>> wrapper for the actuall quadrature routine. The function to be
>> integrated, f, is transferred down to this lowest level. Therefore I
>> would expect that I need to duplicate all routines to which f is passed
>> since they are again called inside f. All other routines, e.g. error
>> check, and everything that does not call or refer to f should indeed be
>> OK.

So, yes, all have to be duplicated if they don't have the RECURSIVE
attribute.

>> I think it is better to stick to the simple Gauss-Legendre quadrature
>> from Numerical Recipes (at least for one of the integration levels)
>> since there is only a single routine that refers to f, and the code is
>> quite simple (but requires to be initialized to get the coefficients).

But Gaussian quadrature is so easy to do in two (or more) dimensions
using loops.

> Numerical Recipes does no enjoy the highest reputation amounst those
> who operate at "professional" levels even if it is a good enjoyable
> read as an > introduction. Your choice!

True, but there aren't that many ways to do Gaussian quadrature.

As far as explaining how things work for the non-specialist, it
seems to be that NR does pretty well. If black-boxes are good
enough, and you have a nice expensive library available, then
don't use NR.

(snip)
> So DQAC is structured a couple levels deep to provide usable
> calling sequences as is sensible and common in good libraries.
> The concept is the same even if the techniicalities are slightly
> more elaborate. No need to duplicate all the utilities called
> atomicly which was you initial complaint.

I missed that one.

>> BTW why it is so complicated for the compiler to deal with
>> recursive calls? Technically, it should be a minor problem
>> to extend gfortran to handle such things.

It isn't on many systems, but it was on many that were popular
in the early days of Fortran. The calling sequence commonly
used on the PDP-8 stores the return address in the first word
of a subroutine, and then branches to the following word.
Not only does it fail in the case of recursion, it fails if the
program is stored in ROM.

> See RECURSIVE in the Fortran manual. It is just a bother on
> cooperative hardware if done from the beginnning and just a
> fair irritant otherwise if planned for. It was not planned
> for in 1960 and now many Fortran optimizations have gotten
> in the habit of not expecting it. So a new keyword warns the
> compiler now. It is often the case that the object code is
> almost recursive except for a few slipups. One is too many
> so the keyword matters as it tells the compiler to actively get
> it right. The less fortunate in hardware have to work harder.

Fortran I didn't even have user written SUBROUTINE and FUNCTION,
as those weren't added until Fortran II. The only choice for
subroutines was the ASSIGNed GOTO, along with appropriate
ASSIGN statements, with the extended range of the DO to allow
for such subroutines.

> The issue is that there is a fair amount of uncooperative
> hardware around or in the collective memory including its
> possible return for specialized purposes.

Also, it is often faster.

-- glen
From: glen herrmannsfeldt on 30 Mar 2010 15:21
Ingo Thies <ingo.thies(a)gmx.de> wrote:
> glen herrmannsfeldt wrote:

>> A better solution is a two dimensional integration routine.
>> In rectangular coordinates, Gaussian quadrature is separable,
>> so all you need is nested DO loops to loop over the two separately.

> If ever possible I would like not to modify the routines themselves but
> use them as a blackbox. But each routine does already contain loops over
> the integrated coordinate.

The first program I worked on that used Gaussian quadrature
did it inside loops. I don't believe that it was even recursive
in the sense here, but it ended up evaluating many of the same
integrals more than once, and was way too slow. With a little
moving code around, and saving the values of previous integrals,
it ran many times faster.

Personally, I don't like black boxes, but some do.

-- glen
First  |  Prev  | 
Pages: 1 2 3 4 5 6 7
Prev: Selecting the first element of a array returned by pack()
Next: size of a derived type containing pointers...