Bessel function computation with SLATEC ( AMOS routines)

[iamjack]

Hi,
I am trying to compute the Bessel function(s) using the SLATEC Amos library available on Netlib. Here is the code I use. The arguments for routine CBESY by Amos are explained in comments header.

c****************************************************************************
PROGRAM BESSELTEST

IMPLICIT NONE
C***BEGIN PROLOGUE CBESY
C***PURPOSE Compute a sequence of the Bessel functions Y(a,z) for
C complex argument z and real nonnegative orders a=b,b+1,
C b+2,... where b>0. A scaling option is available to
C help avoid overflow.
C***LIBRARY SLATEC
C***CATEGORY C10A4
C***TYPE COMPLEX (CBESY-C, ZBESY-C)
C***KEYWORDS BESSEL FUNCTIONS OF COMPLEX ARGUMENT,
C BESSEL FUNCTIONS OF SECOND KIND, WEBER'S FUNCTION,
C Y BESSEL FUNCTIONS
C***AUTHOR Amos, D. E., (SNL)
C***DESCRIPTION
C
C On KODE=1, CBESY computes an N member sequence of complex
C Bessel functions CY(L)=Y(FNU+L-1,Z) for real nonnegative
C orders FNU+L-1, L=1,...,N and complex Z in the cut plane
C -pi<arg(Z)<=pi. On KODE=2, CBESY returns the scaled
C functions
C
C CY(L) = exp(-abs(Y))*Y(FNU+L-1,Z), L=1,...,N, Y=Im(Z)
C
C which remove the exponential growth in both the upper and
C lower half planes as Z goes to infinity. Definitions and
C notation are found in the NBS Handbook of Mathematical
C Functions (Ref. 1).
C
C Input
C Z - Nonzero argument of type COMPLEX
C FNU - Initial order of type REAL, FNU>=0
C KODE - A parameter to indicate the scaling option
C KODE=1 returns
C CY(L)=Y(FNU+L-1,Z), L=1,...,N
C =2 returns
C CY(L)=Y(FNU+L-1,Z)*exp(-abs(Y)), L=1,...,N
C where Y=Im(Z)
C N - Number of terms in the sequence, N>=1
C CWRK - A work vector of type COMPLEX and dimension N
C
C Output
C CY - Result vector of type COMPLEX
C NZ - Number of underflows set to zero
C NZ=0 Normal return
C NZ>0 CY(L)=0 for NZ values of L, usually on
C KODE=2 (the underflows may not be in an
C uninterrupted sequence)
C IERR - Error flag
C IERR=0 Normal return - COMPUTATION COMPLETED
C IERR=1 Input error - NO COMPUTATION
C IERR=2 Overflow - NO COMPUTATION
C (abs(Z) too small and/or FNU+N-1
C too large)
C IERR=3 Precision warning - COMPUTATION COMPLETED
C (Result has half precision or less
C because abs(Z) or FNU+N-1 is large)
C IERR=4 Precision error - NO COMPUTATION
C (Result has no precision because
C abs(Z) or FNU+N-1 is too large)
C IERR=5 Algorithmic error - NO COMPUTATION
C (Termination condition not met)

REAL:: FNU
INTEGER, PARAMETER :: N = 2
COMPLEX :: Z, CY(N), CWRK(N)

INTEGER :: IERR, NZ, KODE

FNU = 0.0
KODE = 2

Z = CMPLX (0.3, 0.4)

CALL CBESY(Z, FNU, KODE, N, CY, NZ, CWRK, IERR)

WRITE(*,*) 'NZ: ', NZ
WRITE(*,*) 'IERR: ', IERR
WRITE(*,*) 'CY: ', CY
STOP

END PROGRAM

c****************************************************************************


For the arguments used in the program above, I get following output:

*****************************************************************************
NZ: 0
IERR: 4
CY: ( 0.0000000 , 0.0000000 ) ( 0.0000000 , 0.0000000 )
*****************************************************************************


Ierr = 4 meaning there is some problem with the input itself. Although I can't figure out what's the problem. Any clues? Much thanks for the suggestions.

J.

[Dick Hendrickson]

Should CWRK come before CY in the argument list?

Dick Hendrickson

[iamjack]

> Dick Hendrickson

>
> Should CWRK come before CY in the argument list?
>

No, as per Netlib's site (http://www.netlib.org/slatec/src/cbesy.f), CWRK is always after CY.
J.

[Dick Hendrickson]

OK, I just suggested that based on the order they defined the arguments
in the comments.

1) I'd try writing out all of the arguments before and after the call.
Sometimes strange things happen that human beings just don't see.

2) Did you actually compile the the cbesy subroutine, or did you link
to it from some existing previously compiled library? Sometimes
libraries are compiled in a way that isn't compatible with the way you
compile the driver routine. Different compiler? Different stack (or
whatever) options?

3) Same question for the routine that cbesy calls. It looked to me
like cbesy was just a driver for some other routine (I've lost the link
to it's name, sigh).

4) Are you sure you (or the cbesy people) didn't compile with some
option like "auto-double" which magically turns ordinary real variables
into double precision?

Dick Hendrickson

[Louisa]

On Sep 17, 2:39 am, iamjack <gcdi...@gmail.com> wrote:

> Ierr = 4 meaning there is some problem with the input itself. Although I can't figure out what's the problem. Any clues? Much thanks for the suggestions.

Might be you have something that's not initialized.
Try initializing everything to zero.

[Andrew Mai]

On 09/16/2012 10:39 AM, iamjack wrote:

> NZ: 0
> IERR: 4
> CY: ( 0.0000000 , 0.0000000 ) ( 0.0000000 , 0.0000000 )
> *************************

I get some odd answers too:

1871 wsbk:/home/mai/amos > gfortran besseltest.f -L. -lamos
1872 wsbk:/home/mai/amos > a.out
NZ: 0
IERR: 4
CY: ( 2.80259693E-44, 0.0000000 ) (-1.25946278E-25, 4.59163468E-41)

I compiled the whole amos lib code with gfortran and put it into an archive. This is, to me, a
standard way to do things and always works.

Andy

[Paul Anton Letnes]

I thought Bessel functions were now a part of the standard?
http://fortranwiki.org/fortran/show/Intrinsic+procedures

Although, I realize that not many compilers support them.
http://fortranwiki.org/fortran/show/Fortran+2008+status
(See "Bessel functions", looks like gfortran and ifort only.)

Paul

[iamjack]

Hi Dick,

>OK, I just suggested that based on the order they defined the arguments
>in the comments.

I agree, it is always good to check the order of the arguments as there is always a possibility of one making an error.

>1) I'd try writing out all of the arguments before and after the call.
> Sometimes strange things happen that human beings just don't see.

Did that, everything checks out except the Bessel function result.

>2) Did you actually compile the the cbesy subroutine, or did you link
>to it from some existing previously compiled library? Sometimes
>libraries are compiled in a way that isn't compatible with the way you
>compile the driver routine. Different compiler? Different stack (or
>whatever) options?

Yes, I compiled the CBESY routine myself using gfortran. To be precise I compiled the whole SLATEC library which has over 1400 useful mathematical functions to create a library file. The SLATEC routines are in Fortran 77 but I believe Fortran 95 can compile these routines just fine. I don't think there is a compatibility issue as I followed similar process to create the LAPACK and BLAS libraries (although I used CMAKE GUI for the LAPACK and BLAS).

>3) Same question for the routine that cbesy calls. It looked to me
>like cbesy was just a driver for some other routine (I've lost the link
>to it's name, sigh).

I think one can call only the CBESY routine without any other routine being called. Of course, CBESY will call some other routines internally.

>4) Are you sure you (or the cbesy people) didn't compile with some
>option like "auto-double" which magically turns ordinary real variables
>into double precision?

No, I did not compiler it auto-double option. Honestly until now I did not know there exists such an option. Wouldn't it be nice setting that flag on so you don't have to go change the precision in the driver routines ( this however seems dangerous, not sure if that's a standard practice).

Many thanks for your speedy response, await further comments.
J.

[iamjack]

> I get some odd answers too:
>
>
>
> 1871 wsbk:/home/mai/amos > gfortran besseltest.f -L. -lamos
>
> 1872 wsbk:/home/mai/amos > a.out
>
> NZ: 0
>
> IERR: 4
>
> CY: ( 2.80259693E-44, 0.0000000 ) (-1.25946278E-25, 4.59163468E-41)
>
>
>
> I compiled the whole amos lib code with gfortran and put it into an archive. This is, to me, a
>
> standard way to do things and always works.
>
>
>
> Andy

It's relief to know that I am not alone :). Still looking.
J.

[iamjack]

>On Monday, September 17, 2012 7:10:22 AM UTC+1, Paul Anton Letnes wrote:
> I thought Bessel functions were now a part of the standard?
>
> http://fortranwiki.org/fortran/show/Intrinsic+procedures

> Although, I realize that not many compilers support them.
>
> http://fortranwiki.org/fortran/show/Fortran+2008+status
>
> (See "Bessel functions", looks like gfortran and ifort only.)
> Paul

Surprising that I did not come across this until now. Thanks a ton Paul for bringing it to notice. Will give this shot immediately and let you all know.
J.

[iamjack]

Problem solved! Fortran 2008 has bessel_jn(n, x) intrinsic.
http://fortranwiki.org/fortran/show/bessel_jn

I compiled my .f95 using this function and of course setting the standard as 2008, i.e.,

gfortran -c -std=f2008 .\testf2008.f95

and then

gfortran -o a .\testf2008.o

Result is bang on compared to Matlab/Mathematica. Thanks again Paul!
J.

[Tobias Burnus]

On 09/17/2012 09:04 AM, iamjack wrote:
>> On Monday, September 17, 2012 7:10:22 AM UTC+1, Paul Anton Letnes wrote:
>> I thought Bessel functions were now a part of the standard?
>> http://fortranwiki.org/fortran/show/Intrinsic+procedures

Yes, they were added in Fortran 2008.

>> Although, I realize that not many compilers support them.
>> http://fortranwiki.org/fortran/show/Fortran+2008+status
>> (See "Bessel functions", looks like gfortran and ifort only.)

I think many compiler support them as vendor extension since a long
time. For instance, g77 and gfortran support them under the name BESJ0,
BESY1 etc. Some other compilers supported also under the same name. Only
the name Bessel_J1 etc. is new since GCC 4.4.

There is one exception: The transformational three-argument version of
BESSEL_JN and BESSEL_YN weren't part of the vendor extension. (They are
useful if one needs all values between N1 and N2; at least in gfortran
they are obtained using a recurrence algorithm. The 3-arg variant is
supported in GCC/gfortran since version 4.6.)


Tobias

[Steven G. Kargl]

How can this solve the problem? The SLATEC routine CBESY takes a COMPLEX
argument while the F2008 intrinsic procedure take a REAL argument. The
two functions then compute a COMPLEX and REAL result, respectively. It
seems that you choose the wrong SLATEC routine. It seems that you wanted

http://www.netlib.org/slatec/src/besy.f

not cbesy.f.

--
steve

[Arjen Markus]

I was intrigued by the value 4 of IERR that CBESY returns.

After some printing and source diving, I found the cause:
I1MACH(9) returns 0. And that is because the function
I1MACH must adjusted for the machine you are working on. Uncomment the
right set of predefined values and it should work.

Regards,

Arjen

[Richard Maine]

Andrew Mai <m...@ucar.edu> wrote:

> I compiled the whole amos lib code with gfortran and put it into an
> archive. This is, to me, a standard way to do things and always works.

For some definitions of "always". While I do agree that this is the most
usual way of doing things (it is certainly how I tend to do it), one of
the many things I disliked about BOCK COMMON is that it often did not
work from an archive without special tricks. These days I'd recommend
avoiding COMMON at all, much less BLOCK COMMON. Many Fortran programmers
of today might not even be aware of BLOCK COMMON. If so, they should
feel lucky and should feel no need to learn about it's quirks unless
they have to deal with old code that uses it.

--
Richard Maine | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle | -- Mark Twain

[Steven G. Kargl]

> I was intrigued by the value 4 of IERR that CBESY returns.
>
> After some printing and source diving, I found the cause:
> I1MACH(9) returns 0. And that is because the function
> I1MACH must adjusted for the machine you are working on. Uncomment the
> right set of predefined values and it should work.

Yeah, that should help OP.

% gfc4x -o z test.f90 -L. -lslatec
% ./z
nz: 0
ierr: 0
cy: (-0.334009618 , 0.449142694 ) (-0.680322826 , 0.635805666 )


--
steve

[Gordon Sande]

On 2012-09-17 12:32:17 -0300, Richard Maine said:

> Andrew Mai <m...@ucar.edu> wrote:
>
>> I compiled the whole amos lib code with gfortran and put it into an
>> archive. This is, to me, a standard way to do things and always works.
>
> For some definitions of "always". While I do agree that this is the most
> usual way of doing things (it is certainly how I tend to do it), one of
> the many things I disliked about BOCK COMMON is that it often did not
> work from an archive without special tricks. These days I'd recommend
> avoiding COMMON at all, much less BLOCK COMMON. Many Fortran programmers
> of today might not even be aware of BLOCK COMMON. If so, they should
> feel lucky and should feel no need to learn about it's quirks unless
> they have to deal with old code that uses it.

Did you mean BLOCK DATA used to set values in a common block? Many
loaders/linkers would not find them in an object library. Several
different tricks were available but in the end they were tricks.
Common blocks have their own problems that seem to have been solved
by making them saved by default. They are a plausible but quite workable
way of structuring data. Much better to use a module with the same
data declared in it. Modules are always the same but common blocks
rely on good typing for all of them to be the same as is usually
intended. One can play some games with the contents which is one of
their problems.

By the way, bock beer is a seasonal delight but seems available year
round any more. (The internet at its extreme good/bad when things are
reduced to carping about spelling/typo errors! Long live spill chuckers!)

[Paul Anton Letnes]

Glad you figured it out!

Paul

[Paul Anton Letnes]

Oops. I didn't even think about that. Although, it seems like Jack
figured it out somehow.

Paul

[Richard Maine]

Gordon Sande <Gordon...@gmail.com> wrote:

> On 2012-09-17 12:32:17 -0300, Richard Maine said:

> > ...BLOCK COMMON....
>
> Did you mean BLOCK DATA...

Oops. Yes. Sorry. Brain fart. Or maybe just senility. And I agree with
all the comments you made about them.

[iamjack]

>
> How can this solve the problem? The SLATEC routine CBESY takes a COMPLEX
>
> argument while the F2008 intrinsic procedure take a REAL argument. The
>
> two functions then compute a COMPLEX and REAL result, respectively. It
>
> seems that you choose the wrong SLATEC routine. It seems that you wanted
>
> http://www.netlib.org/slatec/src/besy.f
>
> not cbesy.f.

> steve


I only wanted to test cbesy.f.

[iamjack]

Steven, Arjen -
Yup, this fixed the problem.