Usage of SLATEC routine CBESY

[gcdiwan]

Hello all,

I am trying to compute the Bessel function of the second kind (Bessel_y) using the SLATEC's Amos library available on Netlib. Here is the SLATEC code I use Link. Below I have pasted my test program that calls SLATEC routine CBESY.

!************************************************************
PROGRAM BESSELTEST
IMPLICIT NONE
REAL:: FNU
INTEGER, PARAMETER :: N = 2, KODE = 1
COMPLEX,ALLOCATABLE :: CWRK ), CY )
COMPLEX:: Z, ci
INTEGER :: NZ, IERR

ALLOCATE(CWRK(N), CY(N))
ci = cmplx (0.0, 1.0)
FNU = 0.0e0
Z = CMPLX(0.3e0, 0.4e0)
CALL CBESY(Z, FNU, KODE, N, CY, NZ, CWRK, IERR)
WRITE(*,*) 'CY: ', CY
WRITE(*,*) 'IERR: ', IERR
STOP
END PROGRAM
!*****************************************************************

And here is the output of the above program:

!*****************************************************************
CY: ( 5.78591091E-39, 5.80327020E-39) ( 0.0000000 , 0.0000000 )
IERR: 4
!*****************************************************************

Ierr = 4 meaning there is some problem with the input itself. To be precise, the IERR = 4 means the following as per the header info in CBESY.f file:

!****************************************************************

! IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
! TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
! CANCE BY ARGUMENT REDUCTION

!****************************************************************

Clearly, CABS(Z) (which is 0.50) or FNU + N - 1 (which is 1.0) are not too large but still the routine CBESY throws the error message number 4 as above. The CY array should have following values for the argument given in above code:

******************************************************************
CY(1) = -0.4983 + 0.6700i
CY(2) = -1.0149 + 0.9485i
******************************************************************

These values are computed using Matlab.
I can't figure out what's the problem when I call CBESY from SLATEC library. Any clues? Much thanks for the suggestions/help.
PS: if it is of any help, I used gfortran to compile, link and then create the SLATEC library file ( the .a file ) which I keep in the same directory as my test program above.
shell command to execute above code:

******************************************************************
gfortran -c BesselTest.f95
gfortran -o a *.o libslatec.a
a
*******************************************************************
GD.

 

[gcdiwan]

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. Appears that the SLATEC library is not needed if one uses the intrinsic functions from Fortran 2008.
GD.

 

[mikrom]

I don't have SLATEC installed, so I couldn't try it, I only looked in the sources and think
that the values you provided doesn't pass the test for range in the subroutine CBESH which is called from CBESY.

http://www.netlib.org/amos/cbesh.f

look at the code:

SUBROUTINE CBESH(Z, FNU, KODE, M, N, CY, NZ, IERR)
...
      NN = N
...
      TOL = AMAX1(R1MACH(4),1.0E-18)
...
      FN = FNU + FLOAT(NN-1)
...
      AZ = CABS(Z)
C-----------------------------------------------------------------------
C     TEST FOR RANGE
C-----------------------------------------------------------------------
      AA = 0.5E0/TOL
      BB=FLOAT(I1MACH(9))*0.5E0
      AA=AMIN1(AA,BB)
      IF(AZ.GT.AA) GO TO 240
      IF(FN.GT.AA) GO TO 240
...

  240 CONTINUE
      NZ=0
      IERR=4
      RETURN
      END

In your case FN = 1.0, AZ = 0.5
IMHO, try to compute in the your program the value of AA and you will see what is the problem.
I1MACH and R1MACH should be contained in SLATEC

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

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

and AMIN1 a MAX1 are fortran intrinsic function (I didn't know it before

http://gcc.gnu.org/onlinedocs/gcc-3.4.6/g77/AMax1-Intrinsic.html

http://gcc.gnu.org/onlinedocs/gcc-3.4.2/g77/AMin1-Intrinsic.html

 

[mikrom]

Now I see you solved it.

 

[gcdiwan]

Yes, the problem is solved at the moment as I only needed Bessel functions for orders ranging from 0 to about 350 and Fortran 2008 does the job. I should still look into the problem you mentioned as at some point of time I might have to come back to SLATEC for other mathematical functions. I feel this because I am not sure if Fortran 2008 has all the mathematical functions from SLATEC implemented as intrinsics.

Thanks for your time and suggestions.