Alan Miller's Fortran Software
[
Subset Selection 
Random Number Generation 
Quad Precision 
Applied Statistics Algorithms 
Logistic Regression 
TOMS Algorithms 
Naval Surface Warfare Center Code 
Miscellaneous 
10Byte Reals for NAS FortranPlus 
F code 
Linear Least Squares 
Links ]
This is an archived copy of the Fortran source code repository of
Alan Miller previously
located at http://users.bigpond.net.au/amiller/
. It is
hosted by Jason Blevins with
permission. The site has been slightly reformatted, but the source
code and descriptions below have not been modified.
All code written by Alan Miller is released into the public
domain. Code written by other authors or from other sources (e.g.,
academic journals) may be subject to other restrictions.
Subset Selection in Regression
The 2nd edition of my book on this subject was published by CRC
Press (Chapman & Hall) in April 2002 (ISBN 1584881712). Dr.
David Smith from the Medical College of Georgia, USA, has notified
me of certain missing references. Here they are in postscript and pdf
formats.
As I can no longer access my ozemail web site, I shall make the
latest versions of my subsets software available here.

lsq.f90 is the latest version of my
uncontrained weighted leastsquares module. It is an upgraded
version of Applied Statistics algorithm AS 274. It uses planar
rotations to produce an uppertriangular factorization. The routine
INCLUD is called once for each case in the data set. It is suitable
for situations in which the leastsquares calculations have to be
updated each time more observations become available. It has
routines for automatically setting up tolerances and testing for
singularities. It uses a version of rankrevealing QR decomposition
for this. This code is now compatible with ELF90.
If results for a subset of predictor variables are required,
those variables are moved to the first positions, and variables to
be excluded are ordered after them. Routine VMOVE moves one
variable; routine REORDR reorders the variables so that those
listed are in the first positions, though not necessarily in the
order specified.
For further details on how to use the module and on methods of
leastsquares calculation refer to the document lsq.txt.
There is a DEMO program which uses a
simple data set fuelcons.dat, and a
nasty test program test1.f90.
If you want to see more tests then download the zip file
tests.zip and the set of data files used by
these tests in testdata.

find_sub.f90 is a module for finding
subsets of variables using a variety of different algorithms.

subset.f90 is a driver program which
uses the above two modules. I have also made available the data set
pollute.dat of mortality rates against
socioeconomic, meteorological and pollution variables for 60
statistical areas in the USA.
For subset selection using the L1norm, that is minimizing the
sum of absolute residuals, here is toms615.f90 which is a translation of TOMS
algorithm 615 to make it ELF90 compatible. There is also a driver
program test615.f90, some test data
test615.dat, and the expected output
test615.out.
For linear regression, but when the regression coefficients must
be positive or zero, there is the Lawson & Hanson nonnegative
leastsquares routine nnls.f90. N.B. Two
call arguments have been removed from the Fortran 77 version. This
routine is ELF90compatible.
I have also added my own nonnegls.f90
routine which is called after a QRfactorization has been formed
using module LSQ. The file t_nnls.f90 uses
both of these routines.
Back to top
Uniform Random Number Generation
 luxury.f90 Another generator of
uniformly distributed random numbers. luxtst.f90 A program to test luxury.f90.
 taus88.f90 L'Ecuyer's 1996 Tausworthe
random number generator, and lfsr113.f90
L'Ecuyer's 1999 Tausworthe random number generator. The first has a
cycle of 2^88 while the second is a little slower but has a cycle
of 2^113. Both are translations from C. N.B. These both assume that
default integers are 32bit.
 lfsr258.f90 A 64bit random number
generator from Pierre L'Ecuyer with a cycle of about 2^258 or more
than 10^77.
 dprand.f90 Nick Maclaren's double
precision random number generator, translated into ELF90compatible
form.
 mt19937.f90 The 'Mersenne Twister'
random number generator from Japan with a cycle of length (2^19937
 1). mt19937a.f90 is a version for
compilers which stop when there are integer overflows, as some do
when compiler check options are enabled for debugging purposes.
Before using this module, you are advised to scan the Japanese web
site particularly under `Initialization'. mt19937.f90 was
revised on 5 February 2002; mt19937a.f90 has not been revised.
 rand3.f90 Yet another random number
generator; this one is based upon an algorithm by Donald Knuth
(1997).
 freq2d.f90 Pairs of uniform random
numbers should be uniformly distributed in the unit square. This is
a simple test program which failed using a new (1998) compiler.
ELF90compatible.
Back to top
Random numbers from other distributions
 random.f90 A module for generating
random numbers from a range of distributions. There is a test
program t_random.f90. Another library
for random number generation can be found at
randlib.
 ran_norm.f90 and rnorm.f90 Generate random numbers from the standard
normal distribution. ran_norm uses more random number calls, but
rnorm uses more logarithms and square roots. On those
machine/compiler combinations which have been tested recently,
rnorm has usually been faster (520%).
 r_gamma.f90 and rgamma.f90 Generate random numbers from the gamma
distribution. The first routine is adapted from Dagpunar's book,
while the second is by Marsaglia and Bray.
 ignpoi.f90 Generate random Poisson
deviates. Requires ran_exp.f90.
 ran_exp.f90 Generate random
exponential deviates, simple method.
 toms780.f90 Generate random
exponential deviates; two algorithms from TOMS algorithm 780.
 test780.f90 A program to compare
speeds of 3 different generators.
 t_kemp_b.f90 Generate random
binomial deviates, including a small driver program.
 zipf.f90 Generate discrete random
deviates having the Zipf (or zeta) distribution, with a small test
program.
 toms668.f90 The file for
r_hyperg.f90 contains TOMS algorithm 668 for generating random
values from the hypergeometric distribution.
 ziggurat.f90 George Marsaglia's
functions for generating random samples from the uniform, normal
and exponential distributions. Translated from C.
Back to top
Quadruple precision
There are multipleprecision (MP) software packages available
when double precision is not adequate, but these are extremely
slow. Quadruple precision gives about twice as many accurate digits
as double precision and is much faster than MP.
 quad.f90 Module for quadruple precision
arithmetic on PC's, using operator overloading. This is designed
for use with Lahey's LF90 compiler (you must use o0, i.e. no
optimization) or with Lahey's ELF90. For Compaq Visual Fortran 5.0
there is a separate version, quad_df.f90,
of the main module now for Compaq which gives the same accuracy as
Lahey's. This version works with SOME Unix f90 compilers and with
Absoft Pro Fortran (Windows version) and NAS FortranPlus. For F
(another subset of F90), Windows 95 version, there is quad.f. For the Lahey/Fujitsu LF95 compiler, use the
Compaq version. N.B. Quadruple precision is built into LF95. It is
easier to use than this version, and more accurate (about 1.8
decimals), but slower.
I have learnt that Visual Basic sets the FPU in the same way as
LF90 and F. If you are calling DLL's compiled in Fortran from VB,
it may be necessary to use quad.f90 even though quad_df.f90 is the
correct module for your Fortran compiler. This problem is likely to
occur with other mixed language programming. With the Portland
PGF90 compiler, I am told that you can use quad_df.f90, but you
must use the undocumented flags Kieee and PC 64, to `turn off
optimization'.
 There is a short test program t_quad.f90, for testing the basic arithmetic
operations (+*/), and other test programs t_cubert.f90, t_cst.f90
and t_logexp.f90 for testing logs,
exponential and trigonometric functions. To test the F version, use
t_quad.f.
 qcomplex.f90 Module for complex
quad. prec. arithmetic.
 From time to time I receive email messages telling me that the
quad. precision constants in my package are wrong. Before you do
the same, look at prtconst.f90.
 There are several example programs  hilbert.f90, which inverts a 10 x 10 Hilbert
matrix, long_gj.f90, for GaussJordan
elimination and sym_eig.f90 which finds
the eigenvalues and vectors of a symmetric matrix using Kaiser's
method.
 q_erf.f90 is a subroutine for
calculating the error function and its complement in quadruple
precision. Use q_erf.f if you are using
F.
 q_lngam.f90 is a function for the log
of the gamma function.
 hero.f90 A program for calculating
coefficients for halfHermite integration.
 Polynomial equations can be extremely difficult to solve, and
this is one place where extra precision is often required. q_pzeros.f90 is a module for solving polynomial
equations, while q_poly20.f90 provides
an example program. Requires the qcomplex module. Click here for more details, and a link to a multiple
precision program.
I am indebted to Keith Briggs (previouly at University of
Cambridge) for access to his package in C++ for quadruple precision
which helped improve the algorithm for calculating
exponentials.
There is a newer package for doubledouble precision and
quaddouble precision (about 64 decimal digits accuracy) in C++,
which has a Fortran frontend, at Click here.
Back to top
Some Applied Statistics Algorithms
For many years, the Royal Statistical Society published
algorithms in its journal `Applied Statistics'. I have translated a
few of these to F90.
 as6.f90 Cholesky factorization of
symmetric positive definite matrix with only lower triangle stored
as a vector (2 versions).
 as7.f90 Inversion of symmetric positive
definite matrix with only lower triangle stored as a vector using
AS6 (2 versions). Includes short test program and interfaces.
 as27.f90 Upper tail area under Student's
tdistribution.
 as60.f90 Calculates the
eigenvalues/vectors of a real symmetric matrix.
 as63.f90 Computes the incomplete beta
function. The file includes a module for the log of the beta
function, which also has the log of the gamma function.
 as66.f90 Area under the normal
curve.
 as91.f90 Returns percentage points of
the chisquared distribution for input tail probabilities and
numbers of degrees of freedom.
 as110 Lpnorm fitting of straight
lines, particularly for 1 <= p <= 2, by extension of an
algorithm of Schlossmacher.
 chirp.f90 The CHIRPZ Fast Fourier
transform for series of arbitrary length. Refer to the journal
article for instructions on its use. Includes FFT for the usual
case when series length is a power of 2. this is AS120.
 as126 Calculates the distribution
function for the range for a sample of size n from the standard
normal distribution.
 as132.f90 Minimize the sum of absolute
deviations for a simple Y on X linear regression.
 as135.f90 MinMax (Linfinity)
estimates for linear multiple regression, includes a test
program.
 as136.f90 A Kmeans clustering
algorithm.
 as152.f90 Cumulative hypergeometric
probabilities.
 as154.f90 Combines AS 154 & 182,
for the exact maximum likelihood estimation of parameters of ARIMA
models using Kalman filtering.
 as155.f90 Distribution of a linear
combination of noncentral chisquared random variables. This code
is NOT compatible with Lahey's ELF90 compiler.
 as157.f90 The runsup and runsdown
tests.
 as177.f90 Calculates the expected
values of normal order statistics.
 as181.f90 Calculates the ShapiroWilks
statistic and its Pvalue. Actually uses ASR94.
 as190.f90 Distribution function and its
inverse, for the studentized range.
 as192.f90 Computes approximate
significance points of a Pearson curve with given first four
moments, or first three moments and left or right boundary.
 as205.f90 Enumeration of all R x C
contingency tables with given row and column totals, and
calculation of hypergeometric probability for each table.
 as207.f90 Fits a general loglinear
model. N.B. The F90 version has NOT been tested by me (Alan
Miller).
 as217.f90 Does the dip calculation for
an ordered vector X using the greatest convex minorant and the
least concave majorant, skipping through the data using the change
points of these distributions. It returns the dip statistic 'DIP'
and the modal interval (XL, XU).
 as227.f90 Generates all possible Nbit
binary codes.
 as239.f90 Integral of the gamma
distribution
 as241.f90 Calculate the normal deviate
(zscore) corresponding to a given area under the normal curve
 as245.f90 Log of the gamma
function
 as260.f90 Calculates the distribution
function of the square (R^2) of the multiple correlation
coefficient.
 as261.f90 Quantiles of the distribution
of R^2.
 as275.f90 Calculation of the
noncentral chisquared distribution function.
 as282.f90 Robust regression using the
least median of squares (LMS) criterion. There is also a test
program: t_as282.f90
 as285.f90 Multivariate normal
probabilities over regions defined in a usersupplied
function.
 as286.zip Parameter estimation for
nonlinear errorsinvariables problems. I have added a third
example, fitting a cylinder, to the two which were in the published
version. This is a ZIP file.
 as290.f90 Generates a rectangular grid
which can be used for contour plotting of confidence limits in
nonlinear regression. Also outputs the Fvalues for the confidence
levels selected. The test program: t_as290.f90 produces the grid for the beanroot
data used in the paper in Applied Statistics: beanroot.dat
 as295.f90 Algorithm for generating
designs of experiments which are close to Doptimal. The user
supplies the set of NCAND candidate design points and the algorithm
picks a subset of N of them. Tha main uses of this algorithm are
likely to be for choosing designs in situations in which there are
blocks of an awkward size, and for augmenting an experiment which
has already been carried out. Two driver programs are provided
here: fr_fac23.f90 for 2^p x 3^q
fractional factorial experiments with possibly unequal block sizes,
and ibquad.f90 for quadratic surface
designs with possibly unequal block sizes.
 as298.f90 Hybrid minimization routine
using simulated annealing and a usersupplied local minimizer. You
should also look at Dr. Tibor Csendes' global optimization
package.
 as304.f90 Fisher's nonparametric
randomization test for two small independent random samples.
 as310.f90 Incomplete beta
function.
 as319.f90 Variable metric unconstrained
function minimization without derivatives.
Back to top
Logistic Regression
I have received several requests for Fortran code to perform
logistic regression, that is to fit:
p = F/(1 + F)
where
p = the probability that a case is in one of two categories
F = exp(b0 + b1.X1 + b2.X2 + ... + bk.Xk)
X1, X2, ..., Xk is a set of k predictors, and
b0, b1, b2, ..., bk is a set of coefficients to be fitted.
 logistic.f90 Module for performing
the weighted least squares calculations. It requires my least
squares module lsq.f90
 t_lgstc1.f90 Driver program which
uses the data set birthwt.dat from
Appendix 1 of Hosmer & Lemeshow's book `Applied Logistic
Regression'.
 t_lgstc2.f90 Driver program which
uses the data set surgical.dat. This
illustrates the use of logistic regression with grouped data. There
is an error in the published data (deliberately NOT corrected here)
which has caused grief with several wellknown statistical
packages.
 t_lgstc3.f90 Driver program which
uses the data set clearcut.dat. This
illustrates the case in which there is a linear boundary such that
all cases on one side are in one category, and all cases on the
other side are in the other category.
 se_lgstc.f90 The standard errors
reported by logistic.f90 are often larger than those reported by
other packages for logistic regression. This simple simulation
program shows that those reported here are about right. It also
illustrates the bias in the slope parameters (they are always
biased towards being too large).
Back to top
Miscellaneous TOMS (and CACM) algorithms
N.B. I have been asked to provide a link to the copyright
policy of the ACM. Loosely paraphrased, this allows use, and
modification, of the TOMS algorithms for most noncommercial
purposes. It also emphasizes that the ACM accepts no responsibility
for the accuracy of the code.
I have updated some of the Transactions on Mathematical Software
(TOMS) algorithms to Fortran 90.
 cacm125.f90. Calculate ordinates and
weights for Gaussian integration.
 cacm395.f90. Calculate 2tailed
probabilities from Student's t distribution.
 cacm396.f90 calculates quantiles of
the tdistribution. These are translations of Geoff Hill's algol
code; Geoff was acting chief of CSIRO Division of Mathematical
Statistics at one time. There is also a test program which may be
useful test395.f90.
 cpoly.f90 This is the classic Jenkins
& Traub algorithm (CACM 419) for the solution of complex
polynomials. q_cpoly.f90 is a version in
quadruple precision.
 rpoly.f90 TOMS algorithm 493 for the
solution of polynomial equations with real coefficients. q_rpoly.f90 is a version in quadruple
precision.
 toms513.f90 In situ transpose of an
MxN matrix.
 toms519.f90 KolmogorovSmirnov
probabilities.
 toms530.f90 Eigenvalues & vectors
of real skewsymmetric matrices (i.e. a(i,j) = a(j,i) ), and of
symmetric matrices with zero diagonal elements.
 toms533.zip A simple package for the
solution of sparse linear equations. It is even shorter and simpler
in Fortran 90 as the user does not supply work arrays. The zipped
file (use PKUNZIP or WinZip to access the contents) contains two
demo programs. The second of these contained an error in the F77
version, and the solution vector was completely wrong.
 toms615.f90 Fitting linear regression
models by minimizing the sum of absolute deviations. i.e. using the
L_1 norm. There are also, a short test program:
 test615.f90 A data file for the test
program:
 test615.dat and a file showing the
output:
 test615.out
 toms639.f90 Integration from zero to
infinity of oscillating functions. Includes driver program.
 toms642.zip Fit a cubic spline to
noisy data using generalized crossvalidation. Data may be unevenly
spaced, and have different weights. Bayesian point error estimates
are generated.
 fexact.zip (TOMS 643) Fisher's exact
test for twoway contingency tables. This is a zipped file
including a driver, data sets, and output. Use pkunzip or winzip to
unpack the files.
 toms644.zip Bessel and Airy functions
for complex arguments. This is a zipped file including test
programs. Use pkunzip or winzip to unpack the files.
 toms654 The incomplete gamma function
and its inverse. Test program: t_incgam.f90
 toms660.f90 Quadratic spline
interpolation of a bivariate function defined by irregularly spaced
data.
 toms661.f90 Quadratic spline
interpolation of a trivariate function defined by irregularly
spaced data.
 toms655.zip IQPACK (Interpolation /
Quadrature). Calculates nodes and weights for 8 types of Gaussian
quadrature (Legendre, Chebyshev, Hermite, rational polynomial,
etc.).
 toms667.f90 Global minimization using
a stochastic algorithm, with 37 test problems test667.f90. Many of the test problems require
hundreds of thousands of function evaluations. See also global.f90
below in the miscellaneous code section.
 toms668.f90 The file for
r_hyperg.f90 contains TOMS algorithm 668 for generating random
values from the hypergeometric distribution.
 toms683.f90 The exponential integral
En(x) for complex argument x.
 toms703.f90 A package for the
solution of systems of ordinary differential equations: DY/DT =
F(Y,T). There is also a test program: test703.f90
 toms707.f90 Confluent hypergeometric
function for complex arguments.
 toms715.zip Special functions by W.J.
Cody and others. Particularly Bessel functions, but also the gamma,
normal distribution, error, Dawson and psi functions.
 toms725.f90 Multivariate normal
integrals for dimension up to 20.
 toms726.zip Walter Gautschi's package
for generating orthogonal polynomials and Gausstype quadrature
rules. The zip file contains a double precision, and most of the
single precision code, plus the set of 11 test programs and
Gautschi's output.
 toms751.f90 Renka's TRIPACK package
for constrained 2D Delaunay triangulation, including test program
(with builtin data) and output files.
 toms757.zip Code for computing some
uncommon special functions, including Abramowitz functions, Airy
function variants, integrals of zeroorder Bessel functions, Debye
functions, Struve functions, Synchrotron radiation and transport
functions, Lobachevski and Stromgen integrals, etc. A test program
is included. This is a zip file.
 toms760.f90 Bicubic interpolation
from a rectangular grid of data. There is also a test program:
test760.f90
 toms766.zip Code for computing
Pade'Hermite and simultaneous Pade' approximants, plus two driver
programs.
 toms778.zip This subroutine solves
boundconstrained optimization problems by using the compact
formula of the limited memory BFGS updates.
 tensolve.zip This contains both TOMS
739 (UNCMIN) for unconstrained minimization, and TOMS 768
(TENSOLVE) for the solution of sets of nonlinear equations.
 toms790.zip Cubic interpolation of
scattered data in two dimensions with continuous second
derivatives.
 toms796.zip This package is for the
numerical inversion of Laplace transforms. The original was in very
poor Fortran 77 which used many nonstandard features, and it
contained many bugs. Each compiler I tried gave different errors,
and then different results when the errors were corrected. This
version is in standard Fortran and gives reproducible results.
 toms813.zip This package is for local
minimization, with or without convex constraints. It requires the
user to supply first derivatives. This version is in standard
Fortran. It uses the socalled spectral projected gradient (SPG)
method.
 toms819.zip This is for the
computation of the Airy functions Ai(z) and Bi(z) and their
derivatives for complex arguments.
Back to top
Code converted from the
Naval Surface Warfare Center Math. Library
 cgamma.f90 Complex gamma
function.
 erf.f90 The error function.
 dcerf.f90 Complex error function &
its complement.
 cexpli.f90 Exponential integral for
complex argument.
 cbsslj.f90 Complex Bessel function
J_{\nu}(z) where both the argument, z, and the order, \nu, are
complex.
 dple.f90 Solution of systems of linear
equations using the HendersonWassyng partial pivot algorithm.
Includes a test program to solve 100 simultaneous equations. The
user must provide a subroutine to supply any requested row of the
matrix.
 dspslv.f90 Solution of sparse systems
of linear equations using Gaussian elimination.
 dsvdc.f90 Calculates the singularvalue
decomposition (SVD) of a real matrix. There is also a test program:

 t_svd.f90
 fft.f90 Fast Fourier Transform (FFT) for
any length of series which has no prime factor greater than 23.
Also the inverse and multivariate FFT.
 hbrd.f90 Solve sets of nonlinear
equations using Powell's Hybrid algorithm.
 specfunc.zip A zip file containing
all of the special functions from the NSWC library.
 polyarea.f90 Calculates the area of
a polygon.
 p_intcpt.f90 Finds the crossing
points of a finite line and a polygon.
 bnd_solv.f90 Solve banded linear
equations using compact storage of the banded matrix. There is a
complex version of this  cbnd_slv.f90
 big_solv.f90 Solves a large set of n
general linear equations using outofcore methods, requiring
storage for about n^2/4 values.
 qagi.f90 Adaptive onedimensional
integration over infinite or semiinfinite ranges (adapted from
QUADPACK).
 qxgs.f90 Adaptive onedimensional
integration over finite ranges (adapted from TOMS algorithm
691).
 inc_gam.f90 The incomplete gamma
function and its inverse. Test program: t_incgam.f90 based on TOMS algorithm 654.
 constant.f90 This is a module of
constants used by some of the functions in the NSWC collection of
routines.
 qsortd.f90 A subroutine which
implements a quicksort algorithm without changing the input array.
It returns an integer array containing the order.
 smplx.f90 Linear programming using the
simplex algorithm. This is a translation of the Fortran 66 program
from the NSWC (Naval Surface Warfare Center) library written by
Alfred Morris. There is also a simple test program t_smplx.f90. Needs the module constant.f90 which defines the precision and
returns certain machine constants.
 dmexp.f90 Calculates the exponential of
a matrix.
 fprob.f90 Cumulative F distribution.
Requires bratio.f90 for the incomplete
gamma function, and constant.f90. N.B.
bratio.f90 contains code for a number of special functions
including the error function, the logarithm of the gamma function,
the logarithm of the beta function, and the digamma function.
bratio was translated from the NSWC library.
 dceigv.f90 Calculate the eigenvalues
and vectors of a general complex matrix. Another routine from the
NSWC library. Needs constant.f90. Based
upon EISPACK routines. There is a test program: t_dceigv.f90
 locpt.f90 Is a point inside a
polygon?
 qtcrt.f90 Solve quadratic, cubic and
quartic equations. Includes a short driver program, and hence
includes the interfaces needed for your program. Adapted from the
NSWC library.
 toeplitz.f90 Solution of Toeplitz
systems of linear equations.
 zeroin.f90 Finds a zero of a
usersupplied function in a specified range (a, b).
Back to top
Miscellaneous code
 qsort.f90 A version of the quicksort
algorithm adapted from Walt Brainerd's code.
 cobyla.f90 Mike Powell's routine for
minimization of nonlinear functions with smooth nonlinear
constraints, using local linear approximations to the
constraints.
 tron.zip Newton's method for large
boundconstrained optimization problems by ChihJen Lin & Jorge
More', a MINPACK2 project. Use PKUNZIP or WINZIP to unpack the
file.
 lm.zip LevenbergMarquardt algorithm for
nonlinear least squares (unconstrained). This is a translation of
the MINPACK routines, LMDER & LMDIF. Use LMDER for functions
which can be differentiated, and LMDIF when it is necessary to use
differences. The ZIPped file includes the MINPACK test programs,
and a simple example fitting a 4parameter logistic.
 conmin.zip The classic CONMIN package
for constrained minimization updated to Fortran 90. Test examples
and the manual are included in the ZIP file. N.B. CONMIN is
included as just one of the algorithms in TOMS algorithm 734 which
can be downloaded from netlib (http://www.netlib.org).
 minim.f90 The NelderMead simplex
algorithm for unconstrained minimization. It does NOT require or
use derivatives. N.B. This is NOT for linear programming! t_minim.f90 is a very simple test program for
minim.f90 which may help users. This is now ELF90compatible.
 primefac.f90 A program to find prime
factors.
 uobyqa.f90 Mike Powell's package for
unconstrained minimization when derivatives are not available.
There is a test/driver program:
t_uobyqa.f90 and a file of results from
the test program:
test.out The documentation, which is in
gzipped postscript, can be downloaded from the Optimization
Decision Tree.
 tn.zip Stephen Nash's truncatedNewton
code for the minimization of continuous functions. It can use
differences instead of derivatives, and bounds may be imposed on
the parameters.
 xdlegf.f90 Legendre functions and
polynomials, from the CMLIB library.
 datesub.f90 Some date manipulation
routines collected together by H.D. Knoble.
 global.f90 At Arnold Neumaier's web
site, this is recommended as the most successful of the global
optimization packages. There is a sample program fit.f90 and the original documentation global.txt for the f77 version. I have included
testfunc.f90 which will eventually
contain all of Neumaier's 30 test functions. N.B. Users of local
optimization packages usually obtain satisfactory convergence after
10s or sometimes 100s of function evaluations. Global optimization
routines usually require many 1000s of function evaluations.
 tensolve.zip A package for solving
sets of nonlinear equations using Robert Schnabel's tensor method.
This is a translation of TOMS algorithm 768. The file was
compressed using PKZIP. It is now compatible with version 4.0 of
ELF90, but misbehaves on example 2.
 ks2.f90 Calculates 1 and 2tail
probabilities for the singlesample KolmogorovSmirnov statistic.
For 2tail probabilities, it uses a combination of the first
algorithm from CACM 487, double the singletail probability, and
the asymptotic distribution. See also TOMS algorithm 519.
 ncr.f90 Calculates number of combinations
of r out of n.
 update.f90 Three very short
subroutines to update the sample mean and sum of squares of
deviations about the mean (and hence update variances or std.
deviations), when one observation is added, dropped or replaced
with another. Designed for fast, repeated use with no checks.
 mvnpack.f90 Alan Genz's package of 4
methods of evaluating multivariate normal integrals.
 bivnorm.f90 A function for
calculating bivariate normal probabilities, extracted from Alan
Genz's package for multivariate normal integration.
 genz2d3d.f90 This is code for
bivariate and trivariate normal integrals which I discovered on
Alan Genz's web site. I have made it compatible with ELF90. It is
more recent (January 2001) than bivnorm above.
 dcuhre.f90 Alan Genz's program for
general multivariate integration, not just of one function but
simultaneously for a vector of functions over the same multivariate
region. There is also a test program dtest1.f90 and a text file dcuhre.txt. which contains the results from the
test program (run in single precision) N.B. While the Fortran 77
version of this code is still at Alan Genz's web site, he is
referring users to the CUBPACK project at: CUBPACK
 tfunc.f90 A module for calculating
bivariate normal probabilities using code by Baughman and
Patefield.
 t_bivnor.f90 Comparison of 3
functions for the bivariate normal  those of Donnelly (CACM 462),
Genz, and Baughman/Patefield.
 hash.f90 and hashord.f90 Hashing routines by Richard Brent and
Donald Knuth, taken from Herman Noble's web site
(http://ftp.cac.psu.edu/pub/ger/fortran/).
 elsunc.f90 Per Lindstroem's package
for nonlinear least squares with upper & lower bounds on
parameter values. d_elsunc.f90
demonstrates this package.
 solvopt.zip The SolvOpt package
minimizes or maximizes nonlinear functions, which may be nonsmooth
and may have constraints, using the socalled method of exact
penalization. This is a zip file containing several driver
programs.
 ga.zip A package for global minimization
using genetic algorithms. Please note the copyright conditions if
you want to use this for commercial purposes. This is a ZIP
file.
 pikaia.zip Another genetic algorithm
for global optimization without derivatives. This one comes from
the High Altitude Observatory of NCAR. This is a ZIP file,
containing a number of interesting examples. There is an excellent
manual, but you will have to download it from their
web site. N.B. The manual is 120 pages long, and in postscript.
Where does the name Pikaia come from? Pikaia gracilens is a
flattened wormlike beast about 5cm long which crawled in the mud
on the sea floor about 530 million years ago!
 cdfcor.zip A package for rational
(Pade) approximation in one and two dimensions. Includes a driver
for 9 test problems, input data and output. This is a ZIP
file.
 sym_band.f90 Find the
eigenvalues/vectors of a symmetric banded matrix stored in compact
form. Based upon EISPACK code.
 strassen.f90 The Strassen fast
matrix multiply algorithm for large matrices. Code downloaded from
the comp.lang.fortran newsgroup.
 varpro.f90 The VARPRO package for
separable nonlinear least squares is for fitting models of the
kind
Y = a1.F1(x, b) + a2.F2(x, b) + ...
in which the a's are linear parameters, the vector b is a vector of
nonlinear parameters, and the F's are nonlinear functions. An
important application is fitting sums of exponential decay terms. A
driver program twoexp.f90 is provided, as
well as the data.
 bvls.f90 Fits a linear model using least
squares but with upper and lower bounds as constraints on each
regression coefficient.
 foldat73.f90 There is a compiler
which accepts Fortran code in fixed format extending beyond column
72, other than sequence numbers which occupy these columns in fixed
form. This simple program takes such code as input data and breaks
long lines after column 72 starting with a continuation character
in column 6.
 DIEHARD A version of George
Marsaglia's random number tests in standard Fortran. You will also
need the files tests.txt and operm5d.ata, and you will need to generate a
binary file containing just over 11 million random 32bit integers
using the random number generator which you want to test. I have
provided a short note diehard.txt and an
example t_taus88.f90 showing how to
generate the binary file using Pierre L'Ecuyer's TAUS88 random
number generator.
 sortchar.f90 Code for sorting
character strings. This was made available on the comp.lang.fortran
newsgroup. Author unknown.
 total_ls.f90 Van Huffel's Total
Least Squares. It can be used for leastsquares regressions in
which there are errors in both the X and Y variables, but the user
must have first scaled the variables so that the errors in all
variables are all the same. Can be used for ODR (orthogonal
distance regression) after the same scaling, otherwise use ODRPACK
(a much larger package) from netlib, or Applied Statistics
algorithm AS 286. test_tls.f90 A test
program for total_ls requiring test_tls.dat. ptlsdoc.txt is Van Huffel's docfile.
 k_smooth.zip Eva Hermann's software
for kernel smoothing, both local (lokern) and global (glkern).
 ewma.f90 Code for updating
exponentiallyweighted moving averages, including updating a
residual sum of squares.
 adventur.zip This is a Fortran 90
version of the classic Adventure game.
 zhangjin.zip The source code from
`Computation of Special Functions' by Zhang & Jin, published by
Wiley, 1996. As well as the `usual' functions such as gamma, error
function, Bessel and Airy functions, there are the confluent
hypergeometric, parabolic cylinder, Mathieu, spheroidal wave, and
various exponential integrals, etc.
 dli.f90 Code from the Slatec library for
the logarithmic integral, Li(x), and the exponential integrals,
Ei(x) and E1(x).
 dcosint.f90 & dsinint.f90 are for evaluating the cosine, Ci(x),
and sine, Si(x), integrals respectively. They are translations of
Numerical Algorithm NA 12. There are corresponding test programs
dcitest.f90 & dsitest.f90
 r_zeta.f90 Riemann's zeta function for
real argument. Adapted from DRIZET in the MATHLIB library from
CERN.
 cincgam.f90 The incomplete gamma
function for complex arguments.
 assndx.f90 The Munkres algorithm for
solution of the assignment problem. Adapted from a routine in the
MATHLIB library from CERN.
 easter.f90 When is Easter in year
XXXX?
 nnes.zip Code for the solution of
simultaneous nonlinear equations using a variety of algorithms. No
documentation, but 3 example programs, one of which (MDR) contains
10 different problems. Warning: This code contains the statement
`Copyright R.S.Bain (1991)', but attempts to contact Rod Bain have
been unsuccesful. This too is a zipped file.
 lanczos.f90 A simple algorithm for
the logarithm of the gamma function.
 to_f90.f90 This program takes Fortran
77 code and converts it to make it look more like Fortran 90.
 kaiser.f90 A simple routine to
calculate the eigenvalues and eigenvectors of a symmetric positive
definite, e.g. a covariance matrix.
 singlton.f90 The classic Singleton
multidimensional FFT algorithm for series whose length is not
necessarily a power of 2.
 fft_simple.f90 A simple Fast
Fourier routine for the case in which the series length is a power
of 2.
 chirp.f90 The ChirpZ algorithm for the
FFT of a series of any length.
 fft235.f90 Fast Fourier Transform for
the case in which the series length is a multiple of some or all of
the integers 2, 3 and 5 (e.g. length = 60 or 240).
 hartly2d.f90 Hartley 2D Fast Fourier
Transform.
 pzeros.f90 Solve polynomial equations using
Aberth's method. Translated from a Fortran 77 algorithm by Dario
Bini published in Numerical Algorithms, vol.13 (1996). A short test
program poly20.f90 is also available.
Dario Bini has a package called MPSolve which does it in multiple
precision. It can be downloaded by ftp from:
http://fibonacci.dm.unipi.it/~bini/ric.html.
 envelope.f90 A simple but efficient
routine for finding 2D convex hulls, i.e. in finding the minimum
polygon to enclose a set of points.
 median.f90 Finds the median of a set
of numbers using a truncated quicksort algorithm.
 hermite.f90 Hermite integration, that
is integration of f(x).p(x) from minus infinity to plus infinity,
where f(x) is the user's function and p(x) = exp(x^2).
 hh.f90 HalfHermite integration, that is
from zero to plus infinity. hh_test.f90
is a test program which also serves as an illustration of how to
use hh. See below for a program to calculate in quadruple precision
the weights and abscissae.
 chi_sq.f90 Chisquared distribution
function. Requires as239.f90 and lanczos.f90.
 hyperg.f90 Calculate hypergeometric
probabilities.
 rhohat.f90 Maximum likelihood
estimation of the shape parameter of the gamma (Erlang)
distribution, including calculation of the derivatives of the
loggamma function.
 fnprod.f90 Computes the distribution
function for the product of two correlated normal variates using
the algorithm of Meeker & Escobar. Needs constant.f90 and qxgs.f90.
 twodqd.f90 A translation of Kahaner
& Rechard's 1984 program for bivariate integration over
triangular regions, including a test program.
 simann.f90 is a module and test
program for simulated annealing based upon the algorithm of
Corana.
 ives.f90 is a routine for generating all
combinations of n objects.
 mace.f90 is an F90 version of Jerry
Friedman's program to estimate multiple optimal transformations for
regression and correlation by alternating conditional expectation
estimates.
Back to top
10byte REAL code for NAS FortranPlus
This exploits the 10byte REAL data type which is supported by
this compiler. At the moment, it only contains a special version of
my quadruple precision package which gives about 38 decimal digit
representation of quadruple precision numbers. N.B. It is extremely
unlikely that this will give correct answers with any other
compiler (e.g. Salford) which supports 10byte REALs.
 quad_nas.f90 Module for quadruple
precision arithmetic using operator overloading.
 There is a short test program t_quad_nas.f90, for testing the basic
arithmetic operations (+*/), and other test programs t_cubert_nas.f90, t_cst_nas.f90, t_logexp_nas.f90 and t_arc_nas.f90,for testing logs, exponential and
trigonometric functions.
 qcomplex_nas.f90 Module for
complex arithmetic in quadruple precision.
 q_cpoly_nas.f90 The classic
Jenkins/Traub algorithm for finding the zeroes of a polynomial
equation with complex coefficients. Originally published as CACM
algorithm 419.
 q_rpoly_nas.f90 The classic
Jenkins/Traub algorithm for finding the zeroes of a polynomial
equation with real coefficients. Originally published as CACM
algorithm 493.
 hilbert_nas.f90 Simple program to
invert a Hilbert matrix.
Back to top
Code for Imagine1's F compiler
Quadruple Precision Code
 quad.f The main module for quadruple
precision.
N.B. This version includes functions for complex arithmetic. I have
only tested it with the Windows version of F. I am told that it
also works with Linux on PCs. On machines which do not use the
Intel (or similar) CPU with an 80bit FPU, try deleting that `11'
from the expression for the real (DP), parameter :: const
The IEEE double precision format has a 53bit mantissa, compared
with 64 bits for the Intel FPU, hence that 11. Problems with other
computers are most likely to occur with the exponential and trig.
functions, so make sure that you test these.
 t_quad.f A simple test program. Various
test quantities, e.g. (a^2  b^2) = (ab)*(a+b) are calculated in
two different ways and the results compared. All of the differences
should be of the order of 10^(30) or zero.
 t_logexp.f and t_cst.f Simple tests of the logarithmic and
trigonometric functions.
 q_erf.f The error function and its
complement. The file includes a test program which differences the
calculated values and compares it with the derivative.
 q_lngam.f The loggamma function,
including a test program.
 q_pzeros.f A module for calculating
the zeroes of polynomials. There is also a demo program q_poly20.f.
Other code for F
 krout.f A module for solving linear
equations and inverting matrices.
 lsq.f A module for multiple linear
regression using least squares, based upon Applied Statistics
algorithm AS 274. There is a demo program
which requires a data set
Back to top
Some linear leastsquares examples and tests
 fit_poly.f90 Fit a polynomial to a
set of (x,y) data.
 quadsurf.f90 Fit a quadratic surface
to a set of (x,y,z) data, i.e. fit:
Z = b0 + b10.X + b01.Y + b20.X^2 + b11.X.Y + b02.Y^2
There is a module ridge.f90 for ridge
regression / regularization. It uses the output from the module lsq
to form the singular value decomposition (SVD) and offers a choice
of 4 different methods of regularization, or the user can input
their own vector to add to the diagonal of the X'Xmatrix.
 wtd_quin.f90 This is a program to
fit a quintic polynomial with exponential weighting of past values.
It was developed for someone who had reason to believe that his
process was well approximated locally by a quintic, but wanted to
give progressively less weight to old observations. The slope and
2nd derivative (acceleration) of the smoothing polynomial are
output after each new case is read. This is suitable for real time
applications. It demonstrates the use of the leastsquares
package's updating algorithm.
 spline5.f90 This program fits quintic
splines with userchosen but evenlyspaced knots. The slope and 2nd
derivative (acceleration) of the smoothing polynomial are output,
but only after all of the data have been read. The smoothed fit and
derivatives at any point are based upon the data from both before
and after the point.
Back to top
Some other useful web sites:
Try the Fortran
Market for general information on Fortran compilers, tutorials,
books and access to some sources of Fortran code, Gary Scott's
Fortran Library web
site.
The Fortran90 FAQ (frequently asked questions) can be obtained
from: F90FAQ.
For an extensive set of routines for sorting and ranking real
numbers see OrderPack.
NAG F90
Software Repository is a source of useful Fortran 90
code.
The ACM collection of TOMS algorithms is a
source of refereed code, mainly in Fortran, for a wide range of
numerical calculations.
A collection of functions and subroutines covering a wide area of
mathematical
John Monahan's site contains the software from his book `Numerical
Methods of Statistics'. It covers a very similar are to this web
site, and is at: Monahan's web
site. operations can be found at JeanPierre Moreau's web site.
For a guide to available mathematical software, refer to GAMS, and for the numerical analysis
FAQ numafaq.
If you are interested in objectoriented programming in Fortran,
you should see oof90.html.
For statistical software in a variety of languages try statlib. For distribution functions,
random number generation and other statistical programs, try
cdflib/ranlib.
For optimization, particularly constrained optimization, see
the Optimization
Decision Tree, or Arnold Neumaier's Global
Optimization web site. The latter also contains many general
links to mathematical and statistical software.
For multivariate normal integrals, and for multiple integration in
general, look at Alan
Genz's home page.
Michel Olagnon's ORDERPACK is for sorting and ranking. It can be
found at: ORDERPACK.
Australian Fortran users will find the web page of Computer Transition
Systems useful. (This company distributes Lahey, Salford,
Edinburgh Portable Compilers, Digital Visual Fortran and other
compilers in Australia.)
To download Lahey's cheap Fortran 90 compiler click on ELF90.
There is an interpreter for a subset of Fortran 90 available from:
Georg Petrich's interpreter. I
have not tried it myself, so comments would be welcome. It can be
downloaded freely, but I understand it is a shareware
product.
For a wide range of code for the FFT, in several computer
languages, see: The FFT
Home Page. and: Another FFT Home Page.
email: amiller@bigpond.net.au
Back to top