next up previous
Next: Graphics Up: Linking FORTRAN to SPLUS Previous: Linking FORTRAN to SPLUS

Is this really necessary?

One should be sure that it is really ``worth it'' before embarking on a FORTRAN project. I like to have an affirmative answer to each of the following questions:

Once you are convinced that this is a necessary evil you should muster all the available tools. Many would say that FORTRAN, like Latin, is a dead language and puts one at an immediate disadvantage. This is certainly true for some applications, particularly those involving a serious graphics component or character string manipulation. But for purely numerical applications FORTRAN continues to serve quite well. I prefer the FORTRAN dialect Ratfor, developed at Bell Labs by Kernighan and Plauger(1976). It provides much of the syntactical structure of C, but represents only a modest investment - the entire literature on learning the dialect is the classic 25 page tutorial written by Kernighan which is a model of clarity. We will not pretend to elaborate on how to write Ratfor; we simply illustrate the language through an example of Ratfor code.


#This is a ratfor implementation of the floyd-revest quantile algorithm-SELECT #Reference: CACM 1975, alg #489, p173, algol-68 version #As originally proposed: mmax=600, and cs=cd=.5 #Translation by Roger Koenker August, 1996. #Calls blas routine dswap subroutine select(n,x,l,r,k,mmax,cs,cd) integer n,m,l,r,k,ll,rr,i,j,mmax double precision x(n),z,s,d,t,cs,cd while(r>l){ if(r-l>mmax){ m=r-l+1 i=k-l+1 fm=dfloat(m) z=log(fm) s=cs*exp(2*z/3) d=cd*sqrt(z*s*(m-s)/fm)*sign(1.,i-m/2) ll=max(l,k-i*s/fm +d) rr=min(r,k+(m-i)*s/fm +d) call select(n,x,ll,rr,k,mmax,cs,cd) } t=x(k) i=l j=r call dswap(1,x(l),1,x(k),1) if(x(r)>t)call dswap(1,x(r),1,x(l),1) while(i<j){ call dswap(1,x(i),1,x(j),1) i=i+1 j=j-1 while(x(i)<t)i=i+1 while(x(j)>t)j=j-1 } if(x(l)==t) call dswap(1,x(l),1,x(j),1) else{ j=j+1 call dswap(1,x(j),1,x(r),1) } if(j<=k)l=j+1 if(k<=j)r=j-1 } return end

This is a pathbreaking algorithm development by Floyd and Rivest (1975) which computes the ordinary sample quantiles in an asymptotically linear number of comparisons in the sample size. The algorithm, as displayed, is a straightforward translation from the Algol given in the original paper. Note the recursive call.

Given the FORTRAN, how is it incorporated into SPLUS? This is inevitably system dependent; we will focus the discussion on the local environment on ragnar. To illustrate this we provide a simple function which calls select in order to compute a sample quantile.


"kuantile"<- function(x, p = 0.5, mmax = 600, cs = 0.5, cd = 0.5) { if(!is.loaded(symbol.For("kuantile"))) dyn.load("src/rqfn/fn.o") n <- length(x) if(p < 0 | p > 1) stop("p outside [0,1]") z <- .Fortran("kuantile", as.integer(n), as.double(x), p = as.double(p), q = double(1), as.integer(mmax), as.double(cs), as.double(cs)) return(z$q) }

This function calls the following ratfor function


#function to compute pth quantile of a sample of n observations subroutine kuantile(n,x,p,q,mmax,cs,cd) integer n,k,l,r,mmax double precision x(n),p,q,cs,cd if(p<0 | p>1) {call dblepr("sparsity bandwidth problem: p=",30,p,1);return} l=1 r=n k=nint(p*n) call select(n,x,l,r,k,mmax,cs,cd) q=x(k) return end

which in turns calls select. Note that the calling sequence requires the character name of the function to be the first argument, the remaining arguments are just as they appear in the FORTRAN. Debugging often requires us to print intermediate results from the FORTRAN, this is somewhat idiosyncratic in current versions of SPLUS, but can be accomplished with the calls

call intpr ("name",4,ivar,n)
call realpr ("rname",5,rvar,p)
call dblepr ("dname",5,dvar,p)
where the character string provides an identifying label, and the final integer argument specifies the number of elements of the variable we desire to print. See the ratfor code above for an example of the use of these calls.

In situations in which there are only one or two subroutines required for a FORTRAN function we may produce an object module corresponding to the source file f.r by invoking the command

f77 -c f.r
This produces a file f.o which can be, in most Unix systems at least, dynamically loaded into SPLUS with the command
dyn.load ("f.o")

In more complicated situations with many subroutines it is convenient to have a makefile to automate the compilation process. Now we illustrate the makefile for the function rqn which underlies the rqfn and RQFN functions mentioned above.


#This is a Makefile for the new frisch-newton RQ routine #The compile flags are intended to optimize for ragnar CFLAGS = -c -xarch=v8 -xchip=super2 -O4 LFLAGS = -r -dn /usr/local/SUNWspro/SC4.0/lib/v8/libsunperf.a

fn.o: fnc.o glob.o sparsity.o ld fnc.o glob.o sparsity.o (LFLAGS) -o fn.o fnc.o: fnc.r f77 (CFLAGS) fnc.r glob.o: glob.r f77 (CFLAGS) glob.r sparsity.o: sparsity.r f77 (CFLAGS) sparsity.r clean: rm fnc.o fn.o

Invoking make induces an evaluation of what elements of the code require recomplilation, and in the load step of the present example links several .o files together into one dynamically loadable module.

Note also in this makefile that the flags for the compile step are chosen to optimize performance for the Sparc 20 architecture of ragnar. In the load step the library libsunperf.a is referenced, enabling access to a broad array of linear algebra routines from LAPACK and elsewhere which are again tuned for good performance on the Sparc 20. See the files in the ragnar directory /usr/local/SUNWspro/READMEs for further details on this library.

The use of the library is an excuse to emphasize the obvious, but often overlooked, fact that it is always desirable to seek out and use well-established library routines rather than risk reinventing them yourself. This is particularly true of basic linear algebra subroutines which have reached an extremely refined state of development in the BLAS provided by LAPACK.


next up previous
Next: Graphics Up: Linking FORTRAN to SPLUS Previous: Linking FORTRAN to SPLUS

Roger Koenker
Sun Aug 31 15:47:37 CDT 1997