(j3.2006) BESSEL_JN in draft of Fortran 2008

[Van Snyder]

On Mon, 16 Jun 2008 20:21:13 +0200 Tobias Burnus wrote:

Fortran 2008 adds the BESSEL_YN and BESSEL_JN intrinsic procedures. 
These are special in so far that they belong into two classes:

Class.
Case (i): BESSEL_JN (N,X) is an elemental function.
Case (ii): BESSEL_JN (N1,N2,X) is a transformational function.

Having a procedure of dual nature is confusing and conceptually 
inconsistent. Furthermore, it increases the implementation burden for 
the compiler vendor.

Additionally, I do not not see the need for the transformational 
function; the same functionality can easily gained by using the 
elemental function.


The reason for the transformational function is that Bessel functions
are evaluated by recurrence, and in particular by a Miller algorith.
When Bessel_JN is needed for n1, n1+1, ..., n2, it is far more efficient
to ask for that, than to ask for Bessel_JN(n1,x), Bessel_JN(n1
+1,x), ..., Bessel_JN(n2) either separately or as an elemental function
reference.  It is probably too much to ask a processor to notice this
special case when presented with an elemental reference in which N is a
rank-1 array having values that are consecutive integers in ascending
order, and X is a scalar.

Processors frequently implement elemental procedure references by
building loops around scalar procedure references, not by passing a
descriptor to a procedure (or a different procedure for each rank) that
has the loops within it.  This is especially true for elemental
procedures having more than one argument.  With the maximum rank being
increased to 15, the latter strategy would require the processor to have
256 procedures on hand.  Thus, a requirement for Bessel_JN(n1,x),
Bessel_JN(n1+1,x), ..., Bessel_JN(n2) that is expressed as an elemental
reference would almost surely be satisfied by references to a scalar
procedure n2-n1+1 times.  This multiplies the amount of work by roughly
(n2-n1+1).

Best regards,
Van Snyder