NMS
Numerical Analysis Library
NMS
is a FORTRAN90 library which
is a good beginner's
mathematical library, with well tested routines for standard problems.
NMS accompanies the text
"Numerical Methods and Software". The book is a modern discussion
of current numerical algorithms and software. The software that comes
with the book has been extracted from standard software libraries,
particularly the SLATEC library. Thus the book is also a good
introduction to the use of a portion of the SLATEC library, which
does not have a widely available solid hardcopy reference.
In particular, NMS includes some or all of:

BLAS1, the Basis Linear Algebra Subprograms, Level 1;

DDRIV, for integrating systems of ODE's;

FFTPACK, for Fast Fourier Transforms,
(Paul Swarztrauber);

FNLIB, Wayne Fullerton's special function package,

FUNPACK, Daniel Amos's special function package,

LINPACK, for solving linear systems.

MACHINE, for reporting values of machine arithmetic constants.

PCHIP, Piecewise Cubic Hermite Interpolation Package,
(Fritsch and Carlson);

QUADPACK, for approximate integration;

UNCMIN, for the unconstrained minimization of a function
of several variables.

XERROR, for handling run time errors.
Languages:
NMS is available in
a FORTRAN90 version.
Related Data and Programs:
BERNSTEIN_POLYNOMIAL,
a FORTRAN90 library which
evaluates the Bernstein polynomials,
useful for uniform approximation of functions;
DISTANCE_TO_POSITION,
a FORTRAN90 program which
estimates the positions of cities based on a citytocity distance table.
It uses UNCMIN from NMS to solve this problem.
FFTPACK5,
a FORTRAN90 library which
contains version 5 of FFTPACK.
MACHINE,
a FORTRAN90 library which
reports the value of machine
arithmetic constants.
ODE,
a FORTRAN90 library which
implements the Shampine and Gordon ODE solver.
QR_SOLVE,
a FORTRAN90 library which
computes the least squares solution of a linear system A*x=b.
QUADPACK,
a FORTRAN90 library which
approximates integrals of functions.
REACTOR_SIMULATION,
a FORTRAN90 program which
is a simple Monte Carlo simulation of the shielding effect of a slab
of a certain thickness in front of a neutron source. This program was
provided as an example with the book "Numerical Methods and Software."
RKF45,
a FORTRAN90 library which
is a RungeKuttaFehlberg ODE solver.
SLATEC,
a FORTRAN90 library which
evaluates many special functions.
XERROR,
a FORTRAN90 library which
is designed to report and handle errors detected during program execution.
Reference:

Milton Abramowitz, Irene Stegun,
Handbook of Mathematical Functions,
National Bureau of Standards, 1964,
ISBN: 0486612724,
LC: QA47.A34.

Donald Amos, SL Daniel, MK Weston,
CDC 6600 subroutines IBESS and JBESS for Bessel functions
I(NU,X) and J(NU,X),
ACM Transactions on Mathematical Software,
Volume 3, pages 7692, 1977.

Paul Bratley, Bennett Fox, Linus Schrage,
A Guide to Simulation,
Second Edition,
Springer, 1987,
ISBN: 0387964673,
LC: QA76.9.C65.B73.

Richard Brent,
Algorithms for Minimization without Derivatives,
Dover, 2002,
ISBN: 0486419983,
LC: QA402.5.B74.

Roger Broucke,
Algorithm 446:
Ten Subroutines for the Manipulation of Chebyshev Series,
Communications of the ACM,
Volume 16, Number 4, April 1973, pages 254256.

Bill Buzbee,
The SLATEC Common Math Library,
in Sources and Development of Mathematical Software,
edited by Wayne Cowell,
PrenticeHall, 1984,
ISBN: 0138235015,
LC: QA76.95.S68.

Carl deBoor,
A Practical Guide to Splines,
Springer, 2001,
ISBN: 0387953663,
LC: QA1.A647.v27.

Jacob Dekker,
Finding a Zero by Means of Successive Linear Interpolation,
in Constructive Aspects of the Fundamental Theorem of Algebra,
edited by Bruno Dejon, Peter Henrici,
Wiley, 1969,
ISBN: 0471203009,
LC: QA212.C65.

John Dennis, Robert Schnabel,
Numerical Methods for Unconstrained Optimization
and Nonlinear Equations,
SIAM, 1996,
ISBN13: 9780898713640,
LC: QA402.5.D44.

Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
LINPACK User's Guide,
SIAM, 1979,
ISBN13: 9780898711721,
LC: QA214.L56.

Bennett Fox,
Algorithm 647:
Implementation and Relative Efficiency of Quasirandom
Sequence Generators,
ACM Transactions on Mathematical Software,
Volume 12, Number 4, December 1986, pages 362376.

Leslie Fox, Ian Parker,
Chebyshev Polynomials in Numerical Analysis,
Oxford Press, 1968,
LC: QA297.F65.

Phyllis Fox, Andrew Hall, Norman Schryer,
Algorithm 528:
Framework for a Portable Library,
ACM Transactions on Mathematical Software,
Volume 4, Number 2, June 1978, page 176188.

Fred Fritsch, Judy Butland,
A Method for Constructing Local Monotone Piecewise
Cubic Interpolants,
SIAM Journal on Scientific and Statistical Computing,
Volume 5, Number 2, June 1984, pages 300304.

Fred Fritsch, Ralph Carlson,
Monotone Piecewise Cubic Interpolation,
SIAM Journal on Numerical Analysis,
Volume 17, Number 2, April 1980, pages 238246.

Charles Gear,
Numerical Initial Value Problems in Ordinary Differential
Equations,
PrenticeHall, 1971,
ISBN: 0136266061,
LC: QA372.G4.

Ron Jones, David Kahaner,
XERROR, The SLATEC Error Handling Package,
Software: Practice and Experience,
Volume 13, Number 3, 1983, pages 251257.

David Kahaner, Cleve Moler, Steven Nash,
Numerical Methods and Software,
Prentice Hall, 1989,
ISBN: 0136272584,
LC: TA345.K34.

Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
Algorithm 539:
Basic Linear Algebra Subprograms for Fortran Usage,
ACM Transactions on Mathematical Software,
Volume 5, Number 3, September 1979, pages 308323.

Pierre LEcuyer,
Random Number Generation,
in Handbook of Simulation,
edited by Jerry Banks,
Wiley, 1998,
ISBN: 0471134031,
LC: T57.62.H37.

Peter Lewis, Allen Goodman, James Miller,
A PseudoRandom Number Generator for the System/360,
IBM Systems Journal,
Volume 8, Number 2, 1969, pages 136143.

George Marsaglia, Wai Wan Tsang,
A fast, easily implemented method for sampling from decreasing or
symmetric unimodal density functions,
SIAM Journal of Scientific and Statistical Computing,
Volume 5, Number 2, June 1984, pages 349359.

Jorge More, Burton Garbow, Kenneth Hillstrom,
User Guide for MINPACK1,
Technical Report ANL8074,
Argonne National Laboratory, 1980.

Frank Olver,
Tables of Bessel Functions of Moderate or Large Orders,
NPL Mathematical Tables, Volume 6,
Her Majesty's Stationery Office, London, 1962.

Robert Piessens, Elise deDonckerKapenga,
Christian Ueberhuber, David Kahaner,
QUADPACK: A Subroutine Package for Automatic Integration,
Springer, 1983,
ISBN: 3540125531,
LC: QA299.3.Q36.

Michael Powell,
A Hybrid Method for Nonlinear Equations,
in Numerical Methods for Nonlinear Algebraic Equations,
edited by Philip Rabinowitz,
Gordon and Breach, 1970,
ISBN13: 9780677142302,
LC: QA218.N85.

Robert Schnabel, John Koontz, Barry Weiss,
A modular system of algorithms for unconstrained minimization,
Technical Report CUCS24082,
Computer Science Department,
University of Colorado at Boulder, 1982.

Lawrence Shampine, Herman Watts,
ZEROIN, a RootSolving Routine,
Technical Report: SCTM70631,
Sandia National Laboratories, September 1970.

Paul Swarztrauber,
Vectorizing the FFT's,
in Parallel Computations,
edited by Garry Rodrigue,
Academic Press, 1982,
ISBN: 0125921012,
LC: QA76.6.P348.

Stephen Wolfram,
The Mathematica Book,
Fourth Edition,
Cambridge University Press, 1999,
ISBN: 0521643147,
LC: QA76.95.W65.
Source Code:
Examples and Tests:
List of Routines:

AIRY_AI_VALUES returns some values of the Airy Ai(x) function.

AIRY_AI_PRIME_VALUES returns some values of the Airy function Ai'(x).

AIRY_BI_VALUES returns some values of the Airy Bi(x) function.

AIRY_BI_PRIME_VALUES returns some values of the Airy function Bi'(x).

ALNGAM computes the log of the absolute value of the Gamma function.

ASYJY computes high order Bessel functions J and Y.

BAKSLV solves A'*x=b where A is a lower triangular matrix.

BERNSTEIN_POLY_VALUES returns some values of the Bernstein polynomials.

BESI0 computes the hyperbolic Bessel function of the first kind, order zero.

BESI0E computes the scaled hyperbolic Bessel function I0(X).

BESJ computes a sequence of J Bessel functions of increasing order.

BESSEL_I0_VALUES returns some values of the I0 Bessel function.

BESSEL_J0_VALUES returns some values of the J0 Bessel function.

BESSEL_J1_VALUES returns some values of the J1 Bessel function.

BESSEL_JN_VALUES returns some values of the Jn Bessel function.

BP01 evaluates the N+1 Bernstein basis functions of degree N on [0,1].

C8VEC_PRINT_SOME prints some of a C8VEC.

C8VEC_UNIFORM_01 returns a unit pseudorandom C8VEC.

CHFDV evaluates a cubic polynomial and its derivative given in Hermite form.

CHFEV evaluates a cubic polynomial given in Hermite form.

CHFIV evaluates the integral of a cubic polynomial in Hermite form.

CHFMC determines the monotonicity properties of a cubic polynomial.

CHKDER checks the gradients of M functions of N variables.

CHLHSN finds the L*L' decomposition of the perturbed model hessian matrix.

CHOLDC finds the perturbed L*L' decomposition of A+D.

COSQB computes the fast cosine transform of quarter wave data.

COSQB1 is a lower level routine used by COSQB.

COSQF computes the fast cosine transform of quarter wave data.

COSQF1 is a lower level routine used by COSQF.

COSQI initializes WSAVE, used in COSQF and COSQB.

COST computes the discrete Fourier cosine transform of an even sequence.

COSTI initializes WSAVE, used in COST.

CSEVL evaluates an N term Chebyshev series.

D1FCN is a dummy routine for evaluating the gradient vector.

D1MACH returns double precision machine constants.

D1MPYQ computes A*Q, where Q is the product of Householder transformations.

D2FCN is a dummy version of a routine that computes the second derivative.

D9LGMC computes the log gamma correction factor.

DAMAX returns the maximum absolute value of the entries in a vector.

DASUM takes the sum of the absolute values of a vector.

DAXPY computes constant times a vector plus a vector.

DDCOR computes corrections to the Y array of DDRIV3.

DDCST sets coefficients used by the core integrator DDSTP.

DDNTL sets parameters for DDSTP.

DDNTP interpolates the Kth derivative of the ODE solution Y at TOUT.

DDOT forms the dot product of two vectors.

DDPSC computes the predicted YH values.

DDPST is called to reevaluate the partial derivatives.

DDRIV1 solves a system of ordinary differential equations.

DDRIV2 solves a system of ordinary differential equations.

DDRIV3 solves a system of ordinary differential equations.

DDSCL rescales the YH array whenever the ODE step size is changed.

DDSTP performs one step of the integration of an ODE system.

DDZRO searches for a zero of a function in a given interval.

DFAULT sets default values for the optimization algorithm.

DFFTB computes a real periodic sequence from its Fourier coefficients.

DFFTB1 is a lower level routine used by DFFTB.

DFFTF computes the Fourier coefficients of a real periodic sequence.

DFFTF1 is a lower level routine used by DFFTF and SINT.

DFFTI initializes WSAVE, used in DFFTF and DFFTB.

DFFTI1 is a lower level routine used by DFFTI.

DGBFA factors a real band matrix by elimination.

DGBSL solves a real banded system factored by DGBCO or DGBFA.

DGECO factors a real matrix and estimates its condition number.

DGEFA factors a real matrix.

DGEFS solves a general N by N system of single precision linear equations.

DGESL solves a real general linear system A * X = B.

DNOR generates normal random numbers.

DSTART is an entry point used to initialize DNOR.

DNRM2 returns the euclidean norm of a vector.

DNSQ finds a zero of a system of N nonlinear functions in N variables.

DNSQE is the easytouse version of DNSQ.

DOGDRV finds the next Newton iterate by the double dogleg method.

DOGLEG finds the minimizing combination of GaussNewton and gradient steps.

DOGSTP finds a new step by the double dogleg algorithm.

DQRANK computes the QR factorization of a rectangular matrix.

DQRDC computes the QR factorization of a real rectangular matrix.

DQRLS solves an linear system in the least squares sense.

DQRLSS solves a linear system in a least squares sense.

DQRSL computes transformations, projections, and least squares solutions.

DROT applies a plane rotation.

DROTG constructs a Givens plane rotation.

DSCAL scales a vector by a constant.

DSFTB computes a "slow" backward Fourier transform of real data.

DSFTF computes a "slow" forward Fourier transform of real data.

DSVDC computes the singular value decomposition of a real rectangular matrix.

DSWAP interchanges two vectors.

EA performs extrapolation to accelerate the convergence of a sequence.

ENORM computes the Euclidean norm of a vector.

ERF_VALUES returns some values of the ERF or "error" function.

ERROR_F computes the error function.

ERROR_FC computes the complementary error function.

EZFFTB computes a real periodic sequence from its Fourier coefficients.

EZFFTF computes the Fourier coefficients of a real periodic sequence.

EZFFTI initializes WSAVE, used in EZFFTF and EZFFTB.

EZFFTI1 is a lower level routine used by EZFFTI.

FDJAC1 estimates an N by N jacobian matrix using forward differences.

FMIN seeks a minimizer of a scalar function of a scalar variable.

FMIN_RC seeks a minimizer of a scalar function of a scalar variable.

FORSLV solves A*x=b where A is lower triangular matrix.

FSTOCD approximates the gradient of a function using central differences.

FSTOFD approximates a derivative by a first order approximation.

FZERO searches for a zero of a function F(X) in a given interval.

GAMLIM computes the minimum and maximum bounds for X in GAMMA(X).

GAMMA computes the gamma function.

GAMMA_VALUES returns some values of the Gamma function.

GL15T estimates the integral of a function over a finite interval.

GRDCHK checks an analytic gradient against an estimated gradient.

HESCHK checks an analytic hessian against a computed estimate.

HOOKDR finds the next Newton iterate by the MoreHebdon method.

HOOKST finds the new step by the MoreHebdon algorithm.

HSNINT provides initial hessian when using secant updates.

I1MACH returns integer machine constants.

I4_SWAP swaps two integer values.

I8_FACTOR factors an integer.

IDAMAX finds the index of the vector element of maximum absolute value.

INBIN takes a real value X and finds the correct bin for it.

INITS estimates the order of an orthogonal series for a given accuracy.

J4SAVE saves variables needed by the library error handling routines.

JAIRY computes the Airy function and its derivative.

LLTSLV solves A*x=b where A = L * L'.

LNSRCH finds a next Newton iterate by line search.

MVMLTL computes y = L * x where L is a lower triangular matrix stored in A.

MVMLTS computes y = A * x where A is a symmetric matrix.

MVMLTU computes y = L' * x where L is a lower triangular matrix.

NUMXER returns the most recent error number.

OPTCHK checks the input to the optimization routine.

OPTDRV is a driver for the nonlinear optimization package.

OPTIF0 provides a simple interface to the minimization package.

OPTSTP: unconstrained minimization stopping criteria

PASSB is a lower level routine used by CFFTB1.

PASSB2 is a lower level routine used by CFFTB1.

PASSB3 is a lower level routine used by CFFTB1.

PASSB4 is a lower level routine used by CFFTB1.

PASSB5 is a lower level routine used by CFFTB1.

PASSF is a lower level routine used by CFFTF1.

PASSF2 is a lower level routine used by CFFTF1.

PASSF3 is a lower level routine used by CFFTF1.

PASSF4 is a lower level routine used by CFFTF1.

PASSF5 is a lower level routine used by CFFTF1.

PCHCE is called by PCHIC to set end derivatives as requested by the user.

PCHCI sets derivatives for a monotone piecewise cubic Hermite interpolant.

PCHCS adjusts the curve produced by PCHIM so it is more "visually pleasing".

PCHDF approximates a derivative using divided differences of data.

PCHEV evaluates a piecewise cubic Hermite or spline function.

PCHEZ carries out easy to use spline or cubic Hermite interpolation.

PCHFD evaluates a piecewise cubic Hermite function.

PCHFE evaluates a piecewise cubic Hermite function at an array of points.

PCHIA evaluates the integral of a piecewise cubic Hermite function.

PCHIC sets derivatives for a piecewise monotone cubic Hermite interpolant.

PCHID evaluates the definite integral of a piecewise cubic Hermite function.

PCHIM sets derivatives for a piecewise cubic Hermite interpolant.

PCHMC: piecewise cubic Hermite monotonicity checker.

PCHQA: easy to use cubic Hermite or spline integration.

PCHSP sets derivatives needed for Hermite cubic spline interpolant.

PCHST: PCHIP signtesting routine.

PCHSW: the PCHCS switch excursion limiter.

Q1DA approximates the definite integral of a function of one variable.

Q1DAX approximates the integral of a function of one variable.

QAGI approximates an integral over an infinite or semiinfinite interval.

QAGIE is called by QAGI to compute integrals on an infinite interval.

QFORM produces the explicit QR factorization of a matrix.

QK15 carries out a 15 point GaussKronrod quadrature rule.

QK15I applies a 15 point Gauss Kronrod quadrature rule.

QPSRT maintains the descending ordering in a list of integral error estimates.

QRAUX1 interchanges two rows of an upper Hessenberg matrix.

QRAUX2 premultiplies an upper Hessenberg matrix by a Jacobi rotation.

QRFAC computes a QR factorization using Householder transformations.

QRUPDT updates a QR factorization.

R1MACH returns single precision machine constants.

R1UPDT retriangularizes a matrix after a rank one update.

R8_SWAP swaps two R8's.

R8VEC_PRINT_SOME prints "some" of an R8VEC.

R8VEC_REVERSE reverses the elements of an R8VEC.

R8VEC_UNIFORM_01 returns a unit pseudorandom R8VEC.

RADB2 is a lower level routine used by RFFTB1.

RADB3 is a lower level routine used by RFFTB1.

RADB4 is a lower level routine used by RFFTB1.

RADB5 is a lower level routine used by RFFTB1.

RADBG is a lower level routine used by RFFTB1.

RADF2 is a lower level routine used by RFFTF1.

RADF3 is a lower level routine used by RFFTF1.

RADF4 is a lower level routine used by RFFTF1.

RADF5 is a lower level routine used by RFFTF1.

RADFG is a lower level routine used by RFFTF1.

RANDOM_INITIALIZE initializes the FORTRAN 90 random number seed.

RESULT prints information about the optimization process.

RUNGE evaluates Runge's function.

RUNGEP evaluates the derivative of Runge's function.

SECFAC updates the hessian by the BFGS factored method.

SECUNF updates a Hessian matrix by the BFGS unfactored method.

SINQB computes the fast sine transform of quarter wave data.

SINQF computes the fast sine transform of quarter wave data.

SINQI initializes WSAVE, used in SINQF and SINQB.

SINT computes the discrete Fourier sine transform of an odd sequence.

SINT1 is a lower level routine used by SINT.

SINTI initializes WSAVE, used in SINT.

SNDOFD approximates a Hessian with a second order finite difference.

TIMESTAMP prints the current YMDHMS date as a time stamp.

TREGUP decides whether to accept the next optimization iterate.

UNCMIN minimizes a smooth nonlinear function of N variables.

UNI generates real uniform random numbers on [0,1).

USTART is an entry into UNI that allows the user to specify the seed.

UNIB ?

XERABT aborts program execution and prints an error message.

XERCLR resets the current error number to zero.

XERCTL allows user control over handling of individual errors.

XERDMP prints the error tables and then clears them.

XERMAX sets the maximum number of times any error message is to be printed.

XERPRT prints a message on each file indicated by xgetua.

XERROR processes a diagnostic error message.

XERRWV processes an error message that includes numeric information.

XERSAV records that an error occurred.

XGETF returns current value of error control flag.

XGETUA returns the unit number(s) to which error messages are being sent.

XGETUN returns the (first) output file to which messages are being sent.

XSETF sets the error control flag.

XSETUA sets up to 5 unit numbers to which messages are to be sent.

XSETUN sets the output file to which error messages are to be sent.

ZFFTB computes the backward complex discrete Fourier transform.

ZFFTB1 is a lowerlevel routine used by ZFFTB.

ZFFTB_2D computes a backward two dimensional complex fast Fourier transform.

ZFFTF computes the forward complex discrete Fourier transform.

ZFFTF1 is a lower level routine used by ZFFTF.

ZFFTF_2D computes a two dimensional complex fast Fourier transform.

ZFFTI initializes WSAVE, used in ZFFTF and ZFFTB.

ZFFTI1 is a lower level routine used by ZFFTI.
You can go up one level to
the FORTRAN90 source codes.
Last revised on 02 April 2008.