Porting SLATEC Library -- Undocumented routines in SALFLIBC

About five years ago, I had noticed and reported ( https://forums.silverfrost.com/Forum/Topic/3293 ) the presence of the Netlib/PORT/SLATEC/BLAS/LAPACK routines D1MACH() and I1MACH() in the FTN95 RTL DLLs SALFLIBC.DLL and SALFLIBC64.DLL. These are undocumented, but appear to have the same functionality as the Netlib routines of the same name. They are not accurate, however, as this test program illustrates.

program tmach
integer i,i1mach
double precision d1mach
integer*8 i8
!
print *,'  i            D1MACH(i)           Expected'
print *
print 10,1,transfer(d1mach(1),i8), transfer(tiny(1.0d0),i8)
print 10,2,transfer(d1mach(2),i8), transfer(huge(1.0d0),i8)
print 10,3,transfer(d1mach(3),i8), transfer(epsilon(1.0d0)/2,i8)
print 10,4,transfer(d1mach(4),i8), transfer(epsilon(1.0d0),i8)
10 format(1x,i2,4x,2Z20.16)
end program

The output from the program (32-bit or 64-bit):

   i            D1MACH(i)           Expected

  1        0010000005017817    0010000000000000
  2        7FEFFFFFF596AC44    7FEFFFFFFFFFFFFF
  3        3C9FFFFFF4178DA3    3CA0000000000000
  4        3CAFFFFFF4178DA3    3CB0000000000000

Today, while attempting to build a SLATEC DLL using FTN95, I ran into a more serious problem. SLATEC makes heavy use of an error reporting routine called XERMSG, which has the following signature:

subroutine xermsg(Librar,Subrou,Messg,Nerr,Level)

The first three arguments are strings, and the last two are integers. It so happens, quite unexpectedly, that SALFLIBC/SALFLIC64 contain their own XERMSG, with a different signature. If the user is not aware of this, or is not able to control the order in which the linker searches for external routines, an EXE built with FTN95 may contain the FTN95 version instead of the Netlib version.

The following program demonstrates the problem.

program driver
integer ndeg,ierr
real coeff(2), root(1)
ndeg = 0
coeff(1) = 1.0
coeff(2) = 2.0

call rpqr79(ndeg,coeff,root,ierr)
end program

subroutine rpqr79(Ndeg,Coeff,Root,Ierr)
!
!
   integer Ierr, Ndeg
   real Coeff(*), Root(*)
!
!   real scale
   Ierr = 0
   if ( abs(Coeff(1)) == 0.0 ) then
      Ierr = 2
      call xermsg('SLATEC','RPQR79','LEADING COEFFICIENT IS ZERO.',2,1)
      return
   end if
!
   if ( Ndeg <= 0 ) then
      Ierr = 3
      call xermsg('SLATEC','RPQR79','DEGREE INVALID.',3,1)
      return
   end if
!
   if ( Ndeg == 1 ) then
      Root(1) = -Coeff(2)/Coeff(1)
      return
   end if
end subroutine rpqr79

If you forget to link in the dozens of SLATEC routines that report errors, the program uses the routine in the compiler RTL DLLs, and the result is a crash:

Invalid argument to numeric intrinsic function at address 1a0094dd

Within file rpqr79.exe
in RPQR79 in line 27, at address 195
in DRIVER in line 8, at address 69

Here is a test program to exhibit what can happen when a user's Fortran program uses one of the undocumented functions in SALFLIC(64).DLL.

program cfnck
implicit none
!
   COMPLEX ci, w23

#ifdef USEEXT
   complex, external :: cacosh          ! SLATEC external function
#else
   complex           :: cacosh          ! Undocumented FTN95 built-in
#endif
!
!     Constants to be used
!
   DATA ci/(0.E0,1.E0)/
!
   w23  = cacosh(1.E0-2.E0*ci)  ! which CACOSH will be linked in?
!
!     Check for possible errors
!
   print '(i3,2x,2es12.4,2x,2es12.4)',23,w23

END PROGRAM cfnck

#ifdef USEEXT
complex function cacosh(z)
   complex z, ci
   complex, external :: casin
   data ci/(0.,1.)/
   real :: piBy2 = Asin(1.0)

   cacosh = ci*(piBy2 - casin(z))
   print *,'SLATEC version of CACOSH used'

   return
end function
#endif

Compiling, running, and output:

ftn95 /cfp tacosh.f90 /lgo
 23    1.5286E+00 -1.1437E+00

Here we may note that the missing code for CACOSH did not cause a linker error, because the builtin function CACOSH# was used rather than the SLATEC routines CACOSH and CASIN. When building SLATEC with its 1900 source files, it would be easy to miss such events.

Running again, but using the SLATEC routine:

ftn95 /cfp /define USEEXT 1 tacosh.f90
ftn95 casin.f90
slink tacosh.obj casin.obj
tacosh
 SLATEC CASIN used
 User provided CACOSH used
 23   -1.5286E+00  1.1437E+00

mecej4

Thanks for the feedback. I think that both of these issues can be fixed. I just need to run the test suite to make sure that there are no side effects.

That sounds encouraging, Paul. When you wrote 'test suite', did you refer to the compiler test suite, or the SLATEC test suite?

Of the 54 test sets in SLATEC, one that gave me more trouble than the others is TEST04, which exercises the special function (FNLIB) part of SLATEC. For your convenience, I have prepared a Zip file of all the source files needed for just this test, and the Zip can be downloaded from:

 https://drive.google.com/file/d/1s8XzKeXZtNPWAHckMvocMyvVLZqwIsoe/view?usp=share_link

I have grouped the files into two directories:

(a) Xlib contains the error reporting package plus the machine parameter files.

(b) T04 contains the source files and a driver for TEST04.

mecej4

Thank you for this useful information.

In this context our test suite only includes functions such as erfc that are 200x standard intrinsics.

Depending on its size, one option would be to add the whole of SLATEC to salflibc.dll/salflibc64.dll.

In addition to the issues discussed earlier in this thread regarding building and running SLATEC tests using FTN95, there is a code generation error that occurs in programs that make heavy use of COMPLEX variables, such as TEST04 of the SLATEC test suite. I had noticed/suspected this behavior earlier, but it was elusive and not reportable.

Today, I was able to pin down the issue in a short reproducer. Whether the bug occurs or not is quite sensitive to the compiler options used. With /64, it occurs almost always regardless of the other options. Without /64, the code runs correctly almost always.

program test04
!
   complex c, ci, w
!
   data ci/(0.0, 1.0)/  ! sqrt(-1)
   data c/(5.98928,4.58112E-02)/   ! expected results
!
   w = cbeta(1.E0 + ci,1.E0 - 2*ci)
   print *,'CBETA returned w = ',w
   print *,'Expected       w = ',c

contains

   complex function cbeta(A,B)
!
      complex A, B
      real, parameter ::  xmax = 36.0
!
! The next two lines may cause the bug, and the bug may hide/disappear
! if the code is changed; for example, if the test clause is evaluated
! just once in an IF THEN ELSE block.
!
      if ( real(A)+real(B) < xmax ) cbeta = cgamma(A)*(cgamma(B)/cgamma(A+B))
      if ( real(A)+real(B) < xmax ) return
!
      cbeta = exp(clbeta(A,B))   ! This line is never reached, but its presence is
                                 ! needed for the bug to exist!
   end function cbeta

   complex function cgamma(z)
      complex z
      cgamma = 3.14*z+2.718
      return
   end function

   complex function clbeta(A,B)
      complex A,B
      clbeta = (0.0,0.0)*(A-B)  ! junk expression to fake function definition
      stop 'dummy function CLBETA entered'
   end function

end program

The output from ftn95 /64 /check:

 CBETA returned w =  (1.00000,0.0000)
 Expected       w =  (5.98928,4.581120E-02)

Without /64, and with other compilers, the output is similar to

 CBETA returned w =         (5.98927689,4.581117630E-02)
 Expected       w =         (5.98928022,4.581119865E-02)

I tried a REAL version of the same code, but that did not show any errors.

Quoted from PaulLaidler mecej4

Depending on its size, one option would be to add the whole of SLATEC to salflibc.dll/salflibc64.dll.

Paul, SLATEC contains about 1,800 files, amounting to 11 megabytes of Fortran 77 source code. The DLL that I built is over 6 megabytes, about two times the size of CLEARWIN64.DLL

mecej4

OK. I will look at the problem with COMPLEX variables.

Would it help if I send you new DLLs now with the fixes for the first two issues?

Would it help? Definitely, and thanks for the offer.

The last time I looked at SLATEC, I gave up after reading the following:

The huge size of the SLATEC library is both a plus (it has everything) and a minus (there's so much here I can't find what I'm looking for!) Since SLATEC is built, in large part, from a number of smaller, specialized libraries, I would strongly recommend that if your interests lie entirely within one of those libraries, you try to find a copy of that library! Another issue to keep in mind is the extraordinary complexity of some of the routines. It is not unusual for a single routine in the SLATEC library to call, directly or indirectly, thirty or forty routines. In part, this is a testimony to the modularity of the routines; however, it can make debugging a nightmare.

I ended up tracking down the source code for EISIPACK (eigenvalue computations), which got me out of a deep hole I had dug for myself. (I believe EISIPACK was superseded by LAPACK.)

Kenneth, your remarks are spot on (the quote is from John Burkardt's repository https://people.math.sc.edu/Burkardt/f_src/slatec/slatec.html ). Some of the tests (there are 54+1 of those in the suite) require up to 70 source files, and it takes a bit of work to find which subset of the 1900 source files are needed for a particular test. I developed scripts to do that task, which only needs to be done once. The distribution provides a couple of text files that display dependency 'trees'.

On the other hand, the source codes are of good quality (but old F77 fixed form), written to a standard pattern, and the number of remaining bugs is quite small. Once a library has been built for a particular compiler, using the library is quite simple. For instance, I tried out the NLS part of the library on the NIST problem set ( https://www.itl.nist.gov/div898/strd/nls/nls_main.shtml ), and found that the algorithms performed so well that I became willing to put some effort into making the software a bit easier to use.

Documentation is rather Spartan, but hypertext versions exist (for example, see http://sdpha2.ucsd.edu/slatec_Source/toc.htm ).

In a few days, we should be able to provide you with a DLL that you can try, and I hope that you will.

When you create the library and place it on main website don't forget to give it very good advertisement, so that it will stand out of crowd of many other libraries

A plus would be also to make as many as possible demo examples

mecej4

I have sent you a pm.

Thanks, I have obtained the new DLLs and run the SLATEC tests with them. The D1MACH results are now accurate and the built-in XERMSG appears to have the same interface as the SLATEC version. Once the complex variable issue discussed above gets fixed, all the SLATEC tests should report PASSED.

Here is a list of exported entry points in SALFLIB64/CLEARWIN64 that match entry names in SLATEC:

  ACOSH   ASINH   ATANH   CACOS   CACOSH  CASIN   CASINH  CATAN
  CATANH  CCOPY   CCOSH   CDIV    CSINH   CTAN    CTANH   D1MACH
  D9B0MP  D9B1MP  D9LGMC  DACOSH  DASINH  DASYJY  DATANH  DBESJ
  DBESY   DBSYNU  DCSEVL  DERF    DERFC   DGAMLM  DGAMMA  DJAIRY
  DSORT   DYAIRY  I1MACH  INITDS  ISORT   XERMSG  ZABS    ZDIV

With the exception of ZABS and ZDIV, the 32-bit SALFLIBC.DLL exports the same entry points.

These lists are shorter than I had expected. Some of these are subroutines rather than functions.

Perhaps I can fix the SLATEC sources so that these are all declared with the EXTERNAL attribute.

mecej4

Thanks for the list.

In a few cases the Silverfrost exports end with one or more $'s or @'s and these functions are unlikely to match similarly named SLATEC functions.

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).

Kenneth_Smith,

Thanks for your kind words.

Yesterday, I posted a note about my progress:

https://forums.silverfrost.com/Forum/Topic/4253

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

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.

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

ftn95 test25.f
slink test25.obj slatec32.dll
echo 3 | test25

The 'echo 3' is needed since each of the SLATEC test programs reads one integer before running. That integer specifies the verbosity of the output, and allowable values are 1, 2 or 3.

Feedback is welcome.

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.