"Dan Nagle" <dannagle@[EMAIL PROTECTED]
> wrote in message
news:2008043018183343658-dannagle@[EMAIL PROTECTED]
> On 2008-04-30 13:17:01 -0400, "James Van Buskirk"
<not_valid@[EMAIL PROTECTED]
>
> said:
>> I further tried all functions that can accept a single REAL
>> variable as input. I was surprised to see in the docs that BESSEL_J0,
>> BESSEL_J1, BESSEL_Y0, and BESSEL_Y1, aren't elemental.
> The j0, j1, y0, y1 functions are elemental.
> The confusion come from the jn and yn functions.
Actually I was looking at the gfortran do***entation:
http://gcc.gnu.org/onlinedocs/gfortran/BESSEL_005fJ0.html#BESSEL_005fJ0
and you are referring to N1723.pdf. You can see there that gfortran
says for BESSEL_J0:
"Arguments:
X The type shall be REAL(*) and it shall be scalar."
and N1723.pdf says:
"Case (i): The result of BESSEL_JN(N,X) is scalar."
But that clearly isn't the right thing to say because
BESSEL_JN(N,X) is elemental, so maybe the verbiage should
be borrowed from ATAN2 which doesn't say anything like that,
so it's clear that ATAN2([1.0,2.0,3.0],1.0) and
ATAN2(1.0,[1.0,2.0,3.0]) have shape [3].
> There are two versions of these. One is elemental.
> That one is the straight y( :) = jn( n, x( :)) version.
> The transformational version is j/yn( n1, n2, x) version.
> This returns all j/y's between order n1 and order n2 inclusive.
> This version was added because many libraries compute
> the Bessel functions via recurrence relations, and you
> may as well get a bunch at once. The idea is that if
> you have an array of coefficients, coef( 1: m), then
> you can use the jn or yn with n2 - n1 + 1 = m, and pass
> the result and your coef array to a polynomial evaluation
> routine with high precision and efficiency.
Now, does transformational mean that BESSEL_JN(N1,N2,X) acts like
SPREAD if X is an array, or does it mean that N1, N2, and X must
all really be scalars in this form of the function? Here is an
example of the kind of confusion caused by the verbiage in N1723.pdf,
section 13.7.22:
C:\gfortran\clf\bessel_test>type bessel_test.f90
program bessel_test
implicit none
integer N(3)
! real(10) x(3)
real x(3)
N = [1,2,3]
x = [1,2,3]
write(*,*) bessel_j0(x(1))
write(*,*) bessel_j0(x)
write(*,*) bessel_jn(n(1),x(1))
write(*,*) bessel_jn(n,x(1))
write(*,*) bessel_jn(n(1),x)
write(*,*) bessel_jn(n,x)
! write(*,*) bessel_jn(n(1),n(3),x(1))
end program bessel_test
C:\gfortran\clf\bessel_test>gfortran bessel_test.f90 -obessel_test
C:\gfortran\clf\bessel_test>bessel_test
0.76519769
0.76519769 0.22389078 -0.26005197
0.44005057
0.44005057 0.11490349 1.95633546E-02
0.44005057 0.57672483 0.33905897
0.44005057 0.35283402 0.30906272
So gfortran treats the elemental case of BESSEL_J0 elementally
even though its do***entation says that X shall be scalar, and
it treats the elemental cases of BESSEL_JN elementally even though
its docs says N and X shall be scalar, and N1723.pdf says the result
of BESSEL_JN(N,X) is scalar.
gfortran doesn't do***ent nor does it implement the transformational
flavor of BESSEL_JN(N1,N2,X), partly because the N1723 do***entation
is vague (I think) and partly because transformational intrinsics
are a royal PITB. Look at how vendors are more than happy to
implement all sorts of elemental intrinsics of perhaps marginal
utility to keep perhaps isolated customers happy, but you don't
see any rush to create, say, a FIRST transformational intrinsic
no matter how useful it might be. In this case, N1723.pdf doesn't
restrict the KIND of N1 and N2 nor of X so that if by the time
N1723.pdf becomes law, gfortran may have 5 integer KINDs and 4
real KINDs, so they would have to have in place a testing regime
for 5X5X4 = 100 different flavors of this transformational
intrinsic, just assuming it doesn't start acting like SPREAD which
we can tell for sure by examining N1723.pdf.
As I was saying before I starting ranting, gfortran doesn't
implement the transformational form of BESSEL_JN, so that if
the line in the above program that invokes it like that is
uncommented, it doesn't compile. A real failure for gfortran
is the fact that the REAL(10) flavors of the Bessel functions
are implemented in the compiler but not in the library so that
if the real(10) declaration of x were in force, the program would
compile but fail at the linking stage -- something that should
be do***ented in the manual and with compiler error messages or
simply fixed.
Definitely some rough edges in the draft and gfortran's initial
implementation can be seen here. Par for the course at such an
early stage I should say.
--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end
|
