Dr. Dobb's is part of the Informa Tech Division of Informa PLC

This site is operated by a business or businesses owned by Informa PLC and all copyright resides with them. Informa PLC's registered office is 5 Howick Place, London SW1P 1WG. Registered in England and Wales. Number 8860726.


Channels ▼
RSS

Open Source

Undocumented Features of PC Fortran Libraries


JAN95: Undocumented Features of PC Fortran Libraries

Undocumented Features of PC Fortran Libraries

Multilanguage vendors sometimes give you more for your money

Kenneth G. Hamilton

Ken has a PhD in physics from the University of California, San Diego and has used numerous Fortran compilers in the pursuit of solutions to problems in solid-state theory, numerical hydrodynamics, signal processing, and random-number generation. He can be contacted on CompuServe at 72727,177.


Language wars are bogus. When you think about it, computers don't really run any high-level language, they run machine code. A compiler is actually just an interface to the programmer, and once things get beyond the first pass or two, all languages start looking the same.

Developers at DEC recognized this a number of years ago. As a result, the VAX came out with one common library (STARLET.OLB) supporting all of the languages. Some vendors of PC compilers are now starting down this path, so that language-support libraries often contain extra "goodies"--things for use with toolsets other than the one that you bought.

Cylindrical Bessel Functions

For example, Microsoft and Watcom both sell compilers for multiple languages, and at the moment, the Fortran customers are the lucky ones. At these companies, developers decided to support that language by taking the C library as a fundamental core, and then adding some modules to cover those requirements unique to Fortran.

Consequently, those companies' Fortran customers already have several useful features spinning around on their hard disks. Probably the most interesting, from the perspective of a numerical analyst, are the cylindrical Bessel functions. Explicit functions are available for J0(x), Y0(x), J1(x), Y1(x), as well as for Jn(x) and Yn(x), where n is a nonnegative integer order. All of these are provided only in DOUBLE PRECISION (REAL*8) form--there are no single-precision versions.

Listing One is a program that exercises Microsoft's Bessel functions. Fortran compilers normally convert the names of external references to upper case. Therefore, I have used a series of INTERFACE statements to make the connection to the library elements, which have lowercase entry-point names that start with underscore characters. The declarations also allow the program to use more descriptive names (containing the character string "BES") in the program itself.

When executed, this program displays several values from the Bessel function of the first kind (J), for orders n=0, 1, 2, and 3. The values for the first two sets of numbers come from specific J0 and J1 functions, while the others are produced by the general-purpose function for Jn.

A PAUSE statement is used to keep the data from scrolling off the screen, and after a key press, the Bessel function of the second kind (Y) is displayed. Again, Y0, Y1, Y2, and Y3 are output, the first two from individual functions, and the second two from the generalized Yn routine.

Watcom has the same features hidden inside its library, and the same program can be used, as long as the INTERFACE blocks are removed and replaced by the c$pragma declarations in Listing Two . In this case, the library members have entry-point names with trailing underscores, in addition to being in lower case, so the pragmas also perform a name conversion.

Watcom users have these features available in both the 16- and 32-bit compilers; Microsoft customers can find the capabilities in the Powerstation compilers for MS-DOS and Windows NT.

Those of little faith can compare the output from this program to the tables given in Chapter 9 of Handbook of Mathematical Functions (see "References"). It is a bit puzzling to me why these modules are standard in C libraries, but not in Fortran libraries: Fortran users--scientists and engineers--are the people most likely to have an application that involves Bessel functions. (For those readers who don't have an application handy right at this moment, Bessel functions are used in problems that involve oscillations with cylindrical symmetry: The normal modes of vibration of a drumhead, for example, are described in terms of the Jn functions.)

Floor, Ceiling, and Hypot

The definition of Fortran-90 includes a requirement for FLOOR and CEILING functions that, given a REAL value, return the next-lower and next-higher integers. Both Microsoft and Watcom have already provided FLOOR and CEILING functions to their customers--they just didn't say so in the manuals.

Listing Three shows a FLOOR and CEILING demonstrator. The only hitch is that, unlike the Fortran-90 definition, these two routines return DOUBLE PRECISION values, and require arguments of the same type. The demo program walks a variable from 4.2 to 6.1 in steps of 0.1, and shows what FLOOR and CEILING return for each of these values of the argument. Watcom users should replace the Microsoft INTERFACE declarations with the two c$pragma declarations in Listing Four .

C libraries also usually include a function called hypot that computes the hypotenuse of a triangle, presumably for those who have forgotten Pythagoras' rule. This is readily accessible from both Fortrans, using the demonstration program in Listing Five ; the Watcom replacement for the INTERFACE block is given in Listing Six .

Character Operations

Because Fortran has traditionally been a batch-oriented language, its default I/O handling is built to process whole records. This is the most efficient method of transferring large amounts of data, as the inevitable word-by-word loop is pushed down as far as possible in the software, as close to the hardware as it can be. When running interactively, however, this usually means that it is necessary to press an Enter key before input data is turned over to a running program.

This is a bit of a nuisance, and often programmers want an application to accept a single keypress without waiting for the Enter key. When running directly under DOS, it is always possible to execute an interrupt to perform such an operation, and higher-level environments (DOS extenders, Windows, Windows NT, OS/2) often claim to translate properly. (Then again, they may not--it's not the normal way of asking a big operating system for something.)

Programs written in C can read and write single characters and can even push a character back into the input buffer to be read again. Since we have a copy of the C library in these Fortran packages, we can do this, too. We'll simply connect to the same system service routines as the C programs_.

In Listing Seven , you can see a program for the Microsoft Powerstation compilers, in which a single character is read from the keyboard using the getche function, which echoes back to the screen. The ungetch function is then used to shove the character back into the keyboard input buffer, so that it can be read a second time by getch. Since getch does not echo the input to the screen, we then use putch to write it to the screen.

Listing Eight contains the c$pragma declarations that should replace the INTERFACE blocks, for use with one of the Watcom Fortran compilers.

C still handles characters as integers and arrays of integers, a bad habit picked up from old Fortran-66. Ever since the release of the ANSI-77 specification, however, Fortran programmers should have been using type CHARACTER variables. While it is still possible to use INTEGERs for character handling, I have inserted the CHAR() and ICHAR() explicit conversion functions in Listing Seven, whenever appropriate.

If you are really just dealing with single characters, then it may be somewhat simpler to stay with INTEGERs in the processing. On the other hand, if you intend to use any of the CHARACTER-oriented routines in the Fortran library (such as INDEX()) to concatenate strings or manipulate substrings, then it makes sense to get these values converted as soon as possible.

The C library also has a set of is functions, which can distinguish digits from letters, upper case from lower, and so on. These all have relatively straightforward names: isalpha reports whether or not a given character is a letter, for example.

In the second half of Listing Eight, I have used some of these routines and written out messages with "T" or "F," denoting which category the input character falls into. The is functions return zero if the condition is false and a nonzero value if the condition is true. I have coded in '.EQ.0' comparisons to produce officially correct LOGICAL values.

I have used the INTERFACE (or c$pragma) declarations to give the Fortran side slightly more readable names. Thus islower becomes IS_LOWER (rather than implying anything negative about program speed) and isupper is converted to IS_UPPER (instead of trying to make you hungry).

If a character is a lowercase letter, then there is a function that will convert it to upper; if it is already in upper case, then there is another that will shift it down. Depending upon the results of the upper- and lowercase tests, the sample program displays the opposite-case letter, using the tolower or toupper functions. These are also used in Listing Seven.

Microsoft has provided wrappers for a few of these keyboard routines, so that their library function GETCHARQQ really connects to getche, and PEEKCHARQQ allows you to check if there is anything there beforehand. They don't provide wrappers for the nonechoing function getch, the ungetch reversal routine, or any of the other character-manipulation functions.

It is theoretically possible to write through the C function printf but, since FORMAT is generally more powerful, I have not provided any examples of this. The adventuresome reader may find numerous other interesting items buried in the software.

Sorting

In its Powerstation compilers, Microsoft provides a wrapper for a quicksort routine. Watcom does not, however, even though its libraries (both 16- and 32-bit) include a sort module. Watcom users, wait no longer: The wrapper you need is in Listing Nine .

Again, I have used a c$pragma declaration to specify how each argument should be passed to QSORT. The interface differs slightly between real- and protected-mode worlds, and so the source code makes use of conditional compilation, activated by the symbol __386__ (which is automatically defined in Watcom's 32-bit compiler only).

The quicksort routine itself requires, as arguments, the name of the array to sort, the number of elements in the array, the size of a single array member (in bytes), and a pointer to a comparison function. To show the versatility of the library routine, my sample program first sorts an INTEGER*2 array, and then a REAL*4 array. For each type of data, you must set up a comparison function (such as ICMPI or ICMPR) and declare it to be EXTERNAL. (This causes its entry-point address to be passed.)

The actual calls to QSORT are in subroutines ISORT and RSORT, which insulate the call process from the main program. If you put ISORT and ICMPI themselves in a source file (with the c$pragma, of course), they can be treated as a black box: An application program could just call ISORT, oblivious to the whole ugly interface process. The same could be done with RSORT and ICMPR, giving a convenient wrapper for floating-point arrays. It should be quite straightforward for you to construct similar packages for other data types as needed using these examples.

Conclusion

I suspect that we will see more of the "common library" approach in PC software. Users of compilers from multilingual companies may wish to use a librarian utility to check for modules with interesting names. Here is a list of suggestions:

  • Salford Software's FTN77 is delivered with a large, well-documented library (including a sort). Apparently, the compiler is itself written in Fortran, so the items that developers needed for that task have been written up for the customers to use as well, leaving little to be discovered. Salford's new FTN90, however, is mingled with a C++ compiler, so that many of the routines I described might work with it.
  • Programs compiled with Silicon Valley Software's Fortran compiler must be linked with that company's Pascal library; this should provide some different possibilities.
  • Borland customers may well find C things in the Pascal box and vice versa.
You never know, there may be free software already on your disk. Good luck and happy spelunking!

References

Abramowitz, Milton and Irene A. Stegun. Handbook of Mathematical Functions. Washington, D.C.: National Bureau of Standards, 1972.

Listing One


      INTERFACE TO REAL*8 FUNCTION DBESJ0[C,ALIAS:"__j0"](X)
      real*8 x
      end
      INTERFACE TO REAL*8 FUNCTION DBESJ1[C,ALIAS:"__j1"](X)
      real*8 x
      end
      INTERFACE TO REAL*8 FUNCTION DBESY0[C,ALIAS:"__y0"](X)
      real*8 x
      end
      INTERFACE TO REAL*8 FUNCTION DBESY1[C,ALIAS:"__y1"](X)
      real*8 x
      end
      INTERFACE TO REAL*8 FUNCTION DBESJN[C,ALIAS:"__jn"](N,X)
      real*8 x
      integer*2 n
      end
      INTERFACE TO REAL*8 FUNCTION DBESYN[C,ALIAS:"__yn"](N,X)
      real*8 x
      integer*2 n
      end
      PROGRAM MBESSEL
      real*8 x,y,z,dbesj0,dbesj1,dbesy0,dbesy1,dbesjn,dbesyn
c
c     Bessel function demo
c     Kenneth G. Hamilton
c
      print 10
   10 format (1X/' Bessel functions of the first kind, orders 0 and 1')
      do i=0,10
        x=0.1D0*dfloat(i)
        y=dbesj0(x)
        z=dbesj1(x)
        print 20, x,y,z
   20   format (' X =',F5.2,', J0 =',1PD20.12,', J1 =',1PD20.12)
      enddo
c
      print 30
   30 format (1X/' Bessel functions of the first kind, orders 2 and 3')
      do i=0,10
        x=0.1D0*dfloat(i)
        y=dbesjn(2,x)
        z=dbesjn(3,x)
        print 40, x,y,z
   40   format (' X =',F5.2,', J2 =',1PD20.12,', J3 =',1PD20.12)
      enddo
      pause
c
      print 50
   50 format (1X/' Bessel functions of the second kind, orders 0 and 1')
      do i=1,10
        x=0.1D0*dfloat(i)
        y=dbesy0(x)
        z=dbesy1(x)
        print 60, x,y,z
   60   format (' X =',F5.2,', Y0 =',1PD20.12,', Y1 =',1PD20.12)
      enddo
c
      print 70
   70 format (1X/' Bessel functions of the second kind, orders 2 and 3')
      do i=1,10
        x=0.1D0*dfloat(i)
        y=dbesyn(2,x)
        z=dbesyn(3,x)
        print 80, x,y,z
   80   format (' X =',F5.2,', Y2 =',1PD20.12,', Y3 =',1PD20.12)
      enddo
c
      stop
      end



Listing Two


c$pragma aux dbesj0 "j0_" parm (value*8)
c$pragma aux dbesj1 "j1_" parm (value*8)
c$pragma aux dbesjn "jn_" parm (value*2, value*8)
c$pragma aux dbesy0 "y0_" parm (value*8)
c$pragma aux dbesy1 "y1_" parm (value*8)
c$pragma aux dbesyn "yn_" parm (value*2, value*8)




Listing Three


      INTERFACE TO REAL*8 FUNCTION FLOOR [C,ALIAS:"_floor"] (X)
      real*8 x
      end
      INTERFACE TO REAL*8 FUNCTION CEILING  [C,ALIAS:"_ceil"] (X)
      real*8 x
      end
      PROGRAM MFLOOR
      real*8 floor, ceiling
      real*8 t, tbelow, tabove
c
      do i=0,20
        t= 4.2D0 + 0.1D0*dfloat(i)
        tbelow = floor(t)
        tabove = ceiling(t)
        print 20, t, tbelow, tabove
   20   format (' T =',F5.2,', Below = ',F5.2,', Above = ',F5.2)
      enddo
c
      stop
      end
      


Listing Four


c$pragma aux floor    "floor_" parm (value*8)
c$pragma aux ceiling  "ceil_"  parm (value*8)



Listing Five


      INTERFACE TO REAL*8 FUNCTION HYPOT [C,ALIAS:"__hypot"] (X,Y)
      real*8 x, y
      end
      PROGRAM HYPE
      real*8 hypot
      real*8 a, b, c
c
      a = 3.D0
      b = 4.D0
      c = hypot(a,b)
      print *, 'a,b,c=',a,b,c
c
      stop
      end



Listing Six


c$pragma aux hypot "hypot_" parm (value*8, value*8)



Listing Seven


      INTERFACE TO INTEGER FUNCTION GETCH[C,ALIAS:"__getch"]()
      end
      INTERFACE TO INTEGER FUNCTION GETCHE[C,ALIAS:"__getche"]()
      end
      INTERFACE TO INTEGER FUNCTION PUTCH[C,ALIAS:"__putch"](IC)
      integer ic
      end
      INTERFACE TO INTEGER FUNCTION UNGETCH[C,ALIAS:"__ungetch"](IC)
      integer ic
      end
      INTERFACE TO INTEGER FUNCTION IS_ASCII[C,ALIAS:"___isascii"](IC)
      integer ic
      end
      INTERFACE TO INTEGER FUNCTION IS_ALNUM[C,ALIAS:"_isalnum"](IC)
      integer ic
      end
      INTERFACE TO INTEGER FUNCTION IS_ALPHA[C,ALIAS:"_isalpha"](IC)
      integer ic
      end
      INTERFACE TO INTEGER FUNCTION IS_CNTRL[C,ALIAS:"_iscntrl"](IC)
      integer ic
      end
      INTERFACE TO INTEGER FUNCTION IS_DIGIT[C,ALIAS:"_isdigit"](IC)
      integer ic
      end
      INTERFACE TO INTEGER FUNCTION IS_LOWER[C,ALIAS:"_islower"](IC)
      integer ic
      end
      INTERFACE TO INTEGER FUNCTION IS_PUNCT[C,ALIAS:"_ispunct"](IC)
      integer ic
      end
      INTERFACE TO INTEGER FUNCTION IS_SPACE[C,ALIAS:"_isspace"](IC)
      integer ic
      end
      INTERFACE TO INTEGER FUNCTION IS_UPPER[C,ALIAS:"_isupper"](IC)
      integer ic
      end
      INTERFACE TO INTEGER FUNCTION IS_XDIGIT[C,ALIAS:"_isxdigit"](IC)
      integer ic
      end
      INTERFACE TO INTEGER FUNCTION TO_LOWER[C,ALIAS:"_tolower"](IC)
      integer ic
      end
      INTERFACE TO INTEGER FUNCTION TO_UPPER[C,ALIAS:"_toupper"](IC)
      integer ic
      end
      PROGRAM CHARS
      implicit integer*4 (i-n)
      integer getch, getche, putch, ungetch
      integer is_ascii, is_alnum, is_alpha, is_cntrl, is_digit
      integer is_lower, is_punct, is_space, is_upper, is_xdigit
      integer to_lower, to_upper
      logical l_lower, l_upper
      character*1 c1, c2
c
c     Demonstration of C routines available in the Fortran library
c     Kenneth G. Hamilton
c
c     Perform console I/O using a single character
c
      c1 = char(getche())               ! Read one character, with echo
      print 10, ichar(c1)               ! Tell us what it is
   10 format (1X/' Character value is ',Z2,'(hex)')
      istat = ungetch(ichar(c1))        ! Put the character back
c
      c2 = char(getch())                ! Reread the "ungotten" character
      print 20                          ! And it it ...
   20 format (1X/' Reread character is:',$)
      istat = putch(ichar(c2))          ! this character!
c
c     What are the properties of the character?
c
      if (is_ascii(ichar(c2)).eq.0) then
        print 30                        ! This is not a good character
   30   format (1X/' The character is non-ASCII')
        stop
      endif
      print 40
   40 format (1X/' The character is in the ASCII set')
c

      print 50, 'Alphanumeric', (is_alnum(ichar(c2)).ne.0)
      print 50, 'Control', (is_cntrl(ichar(c2)).ne.0)
      print 50, 'Digit', (is_digit(ichar(c2)).ne.0)
      print 50, 'Hex Digit', (is_xdigit(ichar(c2)).ne.0)
      print 50, 'Punctuation', (is_punct(ichar(c2)).ne.0)
      print 50, 'White Space', (is_space(ichar(c2)).ne.0)
   50 format (1X,A15,'?',2X,L5)
c
      print 50, 'Alphabetic', (is_alpha(ichar(c2)).ne.0)
c
      l_lower = (is_lower(ichar(c2)).ne.0)      ! .TRUE. if lower case
      print 50, 'Lower Case', l_lower
      if (l_lower) print 60, char(to_upper(ichar(c2)))
   60 format (6X,'(Upper case equivalent is "',A1,'")')
c
      l_upper = (is_upper(ichar(c2)).ne.0)      ! .TRUE. if upper case
      print 50, 'Upper Case', l_upper
      if (l_upper) print 70, char(to_lower(ichar(c2)))
   70 format (6X,'(Lower case equivalent is "',A1,'")')
c
      stop
      end



Listing Eight


c$pragma aux getch     "getch_"    parm ()
c$pragma aux getche    "getche_"   parm ()
c$pragma aux putch     "putch_"    parm (value*4)
c$pragma aux ungetch   "ungetch_"  parm (value*4)
c$pragma aux is_ascii  "isascii_"  parm (value*4)
c$pragma aux is_alnum  "isalnum_"  parm (value*4)
c$pragma aux is_alpha  "isalpha_"  parm (value*4)
c$pragma aux is_cntrl  "iscntrl_"  parm (value*4)
c$pragma aux is_digit  "isdigit_"  parm (value*4)
c$pragma aux is_lower  "islower_"  parm (value*4)
c$pragma aux is_punct  "ispunct_"  parm (value*4)
c$pragma aux is_space  "isspace_"  parm (value*4)
c$pragma aux is_upper  "isupper_"  parm (value*4)
c$pragma aux is_xdigit "isxdigit_" parm (value*4)
c$pragma aux to_lower  "tolower_"  parm (value*4)
c$pragma aux to_upper  "toupper_"  parm (value*4)



Copyright © 1995, Dr. Dobb's Journal


Related Reading


More Insights






Currently we allow the following HTML tags in comments:

Single tags

These tags can be used alone and don't need an ending tag.

<br> Defines a single line break

<hr> Defines a horizontal line

Matching tags

These require an ending tag - e.g. <i>italic text</i>

<a> Defines an anchor

<b> Defines bold text

<big> Defines big text

<blockquote> Defines a long quotation

<caption> Defines a table caption

<cite> Defines a citation

<code> Defines computer code text

<em> Defines emphasized text

<fieldset> Defines a border around elements in a form

<h1> This is heading 1

<h2> This is heading 2

<h3> This is heading 3

<h4> This is heading 4

<h5> This is heading 5

<h6> This is heading 6

<i> Defines italic text

<p> Defines a paragraph

<pre> Defines preformatted text

<q> Defines a short quotation

<samp> Defines sample computer code text

<small> Defines small text

<span> Defines a section in a document

<s> Defines strikethrough text

<strike> Defines strikethrough text

<strong> Defines strong text

<sub> Defines subscripted text

<sup> Defines superscripted text

<u> Defines underlined text

Dr. Dobb's encourages readers to engage in spirited, healthy debate, including taking us to task. However, Dr. Dobb's moderates all comments posted to our site, and reserves the right to modify or remove any content that it determines to be derogatory, offensive, inflammatory, vulgar, irrelevant/off-topic, racist or obvious marketing or spam. Dr. Dobb's further reserves the right to disable the profile of any commenter participating in said activities.

 
Disqus Tips To upload an avatar photo, first complete your Disqus profile. | View the list of supported HTML tags you can use to style comments. | Please read our commenting policy.