forums.silverfrost.com

[Kenneth_Smith]

Mecej4,

Great work you are doing on this. I will certainly give the DLL a try when it becomes available. The functions for modified Bessel functions with complex arguments are of particular interest to myself (they are used extensively in Schelkunoff�s model of coaxial cables with cylindrical shields and/or armours).

[mecej4]

Kenneth_Smith,

Thanks for your kind words.

Yesterday, I posted a note about my progress:

http://forums.silverfrost.com/viewtopic.php?t=4745

I am very pleased with the quality of the SLATEC software, as well as the ability of my software tools to help me detect and fix the small number of bugs in the package ( a score or two, but that is minute in a package with 1900 source files).

The FTN95-compiled 64 bit DLL works very well, both in regard to accuracy as well as speed. There is one compiler bug (relating to complex variables, noted in this thread above) for Silverfrost to fix before the DLL can pass all the 54 SLATEC tests.

Do you need to evaluate Modified Bessel functions with double precision complex arguments? Or will single precision do? Note that there was no general support for double precision complex ("DCOMPLEX") in Fortran 77, so SLATEC may not provide the support that you may want. However, we can discuss the issues and consider whether it is worthwhile to supplement SLATEC by adding double precision complex functions using Fortran 200X intrinsics to the extent possible.

If you are willing to test the DLL and the few usage examples that I have prepared, notwithstanding the bug in FTN95, I can send you download links for the DLL and the example sources.

Yesterday I found this documentation page for SLATEC, and you may wish to browse through it.

http://www.cs.yorku.ca/~roumani/fortran/slatecAPI/keylist.index.html

[Kenneth_Smith]

mecej4,

This is interesting. The first Fortran 77 I used in the mid 1980s supported COMPLEX*16, and I was under the mistaken impression that this was part of the Fortran 77 standard for nearly 40 years!

Almost certainly the application I have in mind for the modified Bessel functions would run into floating point issues if the calculations were performed in single precision.

I have found the Fortran 77 routines from �FORTRAN routines for computation of Special Functions� by Shanjie Zhang, and Jianming Jin, to be useful (these do use COMPLEX*16 � perhaps reinforcing my erroneous understanding of what the 77 standard specified) . The original link for the source code files where I downloaded the routines does not work anymore. These routines are not in the public domain.

They do appear to have been collated into a single file by John Burkardt (I have not checked so see if they reflect the information I have in the errata document in the download I obtained about 8 years ago).

https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.f

Thanks for the link to the SLATEC documentation page, that is very useful � and easy to search.

I have some time at present to experiment with the DLL and I am more than happy to help with this where I can.

[mecej4]

Here is a link to a Zip file containing 32- and 64-bit SLATEC DLLs for use with FTN95. The Zip also contains a small number of test programs that exercise the DLLs.

https://drive.google.com/file/d/1-EGPTavu0RMBuu8ZHY2aTzeGn7vcBJj0/view?usp=sharing

If interested, you may download the various testNN.f files from https://netlib.org/slatec/chk/ and run them with the DLLs. You DO NOT need any of the other source files in that directory, since those routines are already included in the DLL.

To build and run, for example, Test25, the commands are

[Kenneth_Smith]

mecej4,

I have tested each of your examples and only have one comment:

xsdriv1.f90 Array subscripts out of bounds at line 22.

All the other examples run as expected.

I have also compared the SLATEC routines for J0, J1, Y0, and Y1 (real single and double precision arguments) with the corresponding intrinsic BESSEL_J0 etc and don�t see any unexpected results, with the 32 or 64 bit dlls.

It appears to be relatively straightforward to access the library functions.

I shall do some more testing of other functions later. I will also look at what happens when running from within a Plato project environment. I don�t use FTN95 from the command line that often.

[mecej4]

Thanks for the feedback.

Regarding the subscript overrun in XSDRIV1: it is my mistake. Please change the allocation statement to:

[Kenneth_Smith]

mecej4,

The integration subroutine DGAUS8 is also fine for both 32 and 64 bit dlls.

The following program (which you may wish to adapt as a further example) demonstrates DGAUS8 and was developed and run from within Plato.

[mecej4]

Excellent example of the power of integration by parts. I have added it to the examples collection.

Change the last line to "end" or "end function" or "end function fun".

[Kenneth_Smith]

mecej4,
I had some fun yesterday evening examining IFLAG from FZERO and came up with another example for you to consider.

[mecej4]

Nice example, added to collection.

The SLATEC FZERO routine does not require an initial bracket, and your trial values lead to success.

In general, the root found is not necessarily the one that is closest to the trial value provided, but in this case that is not a problem.

It so happens that the NAG library has a routine to do exactly this calculation:

http://monet.nag.co.uk/nagexample/?nagroutine=s17alf+%3A+nagf_specfun_bessel_zeros&loaddata=Load+default+data+file&inputdata=S17ALF+Example+Program+Data%0D%0A+0.0++5++1+%3A+Values+of+A%2C+N+and+MODE%0D%0A

Another implementation, in the CERN library:

https://cds.cern.ch/record/2050906/files/c345.pdf

I gave the online NAG program the input

[Kenneth_Smith]

mecej4,

My initial version was much longer and included a call to an initial routine to bracket the roots, rather than expanding the search space based on iflag returned from FZERO.

I was not ware of these other routines to do the same calculation. I chose BESJ0 since it was simply another function in the library which was well behaved for this purpose.

Here is another example which uses the routines CGECO and CGEDI. I did not expect this to run correctly due to the issues with complex arithmetic. What is interesting is that in this particular case, the results returned using the 64 bit compiler are correct while those with the 32 bit compiler are not correct. This is not what I was expecting (I thought it would be the other way round). Perhaps you can try this example?

Code corrected, as per next post. Now runs correctly.

[mecej4]

Kenneth,

The argument IPVT needs to be an array of dimension (3).

[Kenneth_Smith]

mecej4,

Thanks I should have spotted that. So no problems with this example.

The Scotland vs. Ireland rugby match as on the television when I was looking at this, which may have distracted me. The final score was very upsetting to a Scotsman!

[mecej4]

Kenneth,

If it is any consolation, 7/22 is a good approximation to 1 / Pi.

[PaulLaidler]

The bug relating to test04 above has now been fixed for the next release of FTN95.