maths_library Module

Routine to find y0 such that y0 = y(x0) given a set of values x(1:n), y(1:n) author: P J Knight, CCFE, Culham Science Centre x0 : input real : x value at which we want to find y x(1:n) : input real array : monotonically de/increasing x values y(1:n) : input real array : y values at each x n : input integer : size of array This routine performs a simple linear interpolation method to find the y value at x = x0. If x0 lies outside the range [x(1),x(n)], the y value at the nearest 'end' of the data is returned. None

Arguments

Return Value real(kind=dp)

Arguments

Return Value real(kind=dp)

Calculates the gamma function for arbitrary real x author: P J Knight, CCFE, Culham Science Centre x : input real : gamma function argument This routine evaluates the gamma function, using an asymptotic expansion based on Stirling's approximation. http://en.wikipedia.org/wiki/Gamma_function T&M/PKNIGHT/LOGBOOK24, p.5

Arguments

Return Value real(kind=dp)

Arguments

Return Value integer

Routine to compute XY where X and Y are vectors author: Jack Dongarra, Linpack author: P J Knight, CCFE, Culham Science Centre n : input integer : order of the matrices sx, sy sx(nincx) : input real array : first vector incx : input integer : interval in storage between sx array elements sy(nincy) : input real array : second vector incy : input integer : interval in storage between sy array elements This routine performs the dot product of two vectors, i.e. calculates the sum from i=1 to N, of X(i)Y(i). !

Arguments

Return Value real(kind=dp)

Routine to finds the index of the array element having the maximum absolute value author: Jack Dongarra, Linpack author: P J Knight, CCFE, Culham Science Centre n : input integer : order of the matrix sx sx(n*incx) : input real array : array being checked incx : input integer : interval in storage between sx array elements This routine finds the array element with the maximum absolute value, and returns the element index. !

Arguments

Return Value integer

Arguments

Return Value real(kind=dp)

Calculates machine (computing) parameters author: P J Knight, CCFE, Culham Science Centre i : input integer : Switch for return value: i=1 : B(1 - P), the machine precision i=2 : B(EMIN - 1), the smallest magnitude i=3 : BEMAX*(1 - B(-P)), the largest magnitude where the machine being used has P base B digits, and its smallest and largest exponents are EMIN and EMAX, respectively. This routine evaluates the numerical machine parameters of the computer being used to run the program, as defined above.

Note that the values of these parameters can be found for a given machine if the Mark 12 or later NAg library is installed on it.

SPMPAR(1) is equivalent to X02AJF();
SPMPAR(2) is equivalent to X02AKF();
SPMPAR(3) is equivalent to X02ALF(). Metcalf and Reid, Fortran 90/95 Explained, 2nd Edition (section 8.7.2)

Arguments

Return Value real(kind=dp)

Arguments

Return Value logical

Arguments

Return Value logical

Arguments

Return Value character(len=2)

Arguments

Return Value character(len=3)

Routine to integrate a 1-D array of y values, using a process similar to Simpson's Rule, and assuming equally-spaced x values. It returns the integral at all tabulated points. author: P J Knight, CCFE, Culham Science Centre dx : input real : (constant) spacing between adjacent x values y(1:n) : input real array : y values to be integrated inty(1:n) : input/output real array : calculated integral (see below) n : input integer : length of arrays y and inty This routine uses a process similar to (but not quite the same as) Simpson's Rule to integrate an array y, returning the integral up to point i in array element inty(i). Note that the first element of inty is not calculated, and must be set to the required value on entry. Usually, but not always, this value will be zero. The original source for this algorithm is not known...

Arguments

Routine to integrate a 1-D array of y values using the Extended Simpson's Rule, assuming equally-spaced x values author: P J Knight, CCFE, Culham Science Centre dx : input real : (constant) spacing between adjacent x values y(1:n) : input real array : y values to be integrated integral : output real : calculated integral n : input integer : length of array y This routine uses Simpson's Rule to integrate an array y. If n is even, routine sumup2 is called to perform the calculation.

Note: unlike sumup1 and sumup2, this routine returns only the complete integral, not the intermediate values as well. None

Arguments

Routine to extract the lower triangular part of a square matrix author: P J Knight, CCFE, Culham Science Centre a(n,n) : input real array : input matrix n : input integer : number of rows and columns in A a(n,n) : output real array : lower triangular part of A This routine extracts the lower triangular part of a square matrix, excluding the diagonal, into a new square matrix. The remainder of this matrix contains zeroes on exit. None

Arguments

Routine that calculates the complete elliptic integral of the first and second kinds author: P J Knight, CCFE, Culham Science Centre sqk : input real : square of the elliptic modulus kk : output real : complete elliptic integral of the first kind ek : output real : complete elliptic integral of the second kind This routine calculates the complete elliptic integral of the first and second kinds.

The method used is that described in the reference, and the code is taken from the Culham maglib library routine FN02A. Approximations for Digital Computers, C. Hastings, Princeton University Press, 1955

Arguments

Estimate the integral of fun(x) from a to b author: P J Knight, CCFE, Culham Science Centre fun : external function : integrand function subprogram fun(x) a : input real : lower limit of integration b : input real : upper limit of integration (b may be less than a) abserr : input real : absolute error tolerance (should be non-negative) relerr : input real : relative error tolerance (should be non-negative) result : output real : approximation to the integral hopefully satisfying the least stringent of the two error tolerances errest : output real : estimate of the magnitude of the actual error nofun : output integer : number of function values used in calculation flag : output real : Reliability indicator; if flag is zero, then result probably satisfies the error tolerance. If flag is xxx.yyy , then xxx = the number of intervals which have not converged and 0.yyy = the fraction of the interval left to do when the limit on nofun was approached. This routine estimates the integral of fun(x) from a to b to a user provided tolerance. An automatic adaptive routine based on the 8-panel Newton-Cotes rule. http://www.netlib.org/fmm/index.html : Computer Methods for Mathematical Computations, G E Forsythe, M A Malcolm, and C B Moler, Prentice-Hall, Englewood Cliffs, New Jersey 1977, ISBN 0-13-165332-6

Arguments

Routine to solve the linear equation system Ax = b author: P J Knight, CCFE, Culham Science Centre a(ndim,ndim) : in/out real array : array A ndim : input integer : dimension of a b(ndim) : in/out real array : RHS vector x(ndim) : output real array : solution for Ax = b This routine solves the linear equation Ax = b. It calls (local copies of) the linpack routines sgefa and sgesl. !

Arguments

Matrix inversion routine author: Roger L. Crane, Kenneth E. Hillstrom, Michael Minkoff; Linpack author: P J Knight, CCFE, Culham Science Centre h(ih,ih) : input/output real array : On input, matrix to be inverted On output, the calculated inverse ih : input integer : array size n : input integer : order of H; n <= ih ipvt(n) : output integer array : pivot vector info : output integer : info flag = 1 normal return = 2 H matrix is singular This routine inverts the matrix H by use of linpack software. !

Arguments

Calculates +/-C +/- (X.dot.Y) for arrays X, Y and scalar C author: Roger L. Crane, Kenneth E. Hillstrom, Michael Minkoff; Linpack author: P J Knight, CCFE, Culham Science Centre x(ixn) : input real array : X array ix : input integer : interval in storage between X array elements y(iyn) : input real array : Y array iy : input integer : interval in storage between Y array elements c : input real : C value total : output real : computed result n : input integer : number of terms in the dot product iflag : input integer : switch = 0 +c + (x dot y) is computed = 1 +c - (x dot y) is computed = 2 -c + (x dot y) is computed = 3 -c - (x dot y) is computed This subroutine computes total = (plus or minus c) plus or minus the dot product of x and y by invoking the basic linear algebra routine dot. !

Arguments

Routine to solve the the real system Ax = b or transp(A).x = b author: Cleve Moler, University of New Mexico, Argonne National Lab. author: P J Knight, CCFE, Culham Science Centre a(lda,n) : input real array : output from sgefa lda : input integer : leading dimension of the array A n : input integer : order of the matrix A ipvt(n) : input integer array : pivot vector from sgefa b(n) : input/output real array : RHS vector on input, solution vector x on output job : input integer : switch = 0 to solve Ax = b , = nonzero to solve transp(A)x = b where transp(A) is the transpose This routine solves the real system Ax = b or transp(A)x = b using the factors computed by sgefa.

A division by zero will occur if the input factor contains a zero on the diagonal. Technically this indicates singularity but it is often caused by improper arguments or improper setting of lda. !

Arguments

Routine to factor a real matrix by Gaussian elimination author: Cleve Moler, University of New Mexico, Argonne National Lab. author: P J Knight, CCFE, Culham Science Centre a(lda,n) : input/output real array : On entry, matrix to be factored. On exit, an upper triangular matrix and the multipliers which were used to obtain it. The factorization can be written A = L*U where L is a product of permutation and unit lower triangular matrices and U is upper triangular. lda : input integer : leading dimension of the array A n : input integer : order of the matrix A ipvt(n) : output integer array : pivot indices info : output integer : info flag = 0 normal completion = k if u(k,k) == 0.0 This routine factors a real matrix by Gaussian elimination. !

Arguments

Routine to compute the determinant and inverse of a matrix author: Cleve Moler, University of New Mexico, Argonne National Lab. author: P J Knight, CCFE, Culham Science Centre a(lda,n) : input/output real array : On entry, output from sgefa. On exit, the inverse if requested, otherwise unchanged lda : input integer : leading dimension of the array A n : input integer : order of the matrix A ipvt(n) : input integer array : pivot vector from sgefa det(2) : output real array : determinant of original matrix if requested, otherwise not referenced. Determinant = det(1) * 10.0**det(2) with 1.0 .le. abs(det(1)) .lt. 10.0 or det(1) .eq. 0.0 . job : input integer : switch for required outputs = 11 both determinant and inverse. = 01 inverse only. = 10 determinant only. This routine computes the determinant and inverse of a matrix using the factors computed by (SGECO or) SGEFA.

A division by zero will occur if the input factor contains a zero on the diagonal and the inverse is requested. !

Arguments

Routine to scale a vector by a constant author: Jack Dongarra, Linpack author: P J Knight, CCFE, Culham Science Centre n : input integer : order of the matrix sx sa : input real array : constant multiplier sx(n*incx) : input/output real array : On entry, matrix to be scaled; On exit, the scaled matrix incx : input integer : interval in storage between sx array elements This routine scales a vector by a constant, using unrolled loops for increments equal to 1. !

Arguments

Routine to scale a vector by a constant, then add another vector author: Jack Dongarra, Linpack author: P J Knight, CCFE, Culham Science Centre n : input integer : order of the matrices sx, sy sa : input real array : constant multiplier sx(nincx) : input real array : matrix to be scaled incx : input integer : interval in storage between sx array elements sy(nincy) : input/output real array : On entry, matrix being added; On exit, the final result incy : input integer : interval in storage between sy array elements This routine calculates sa*sx(:) + sy(:), using unrolled loops for increments equal to 1. !

Arguments

Routine to interchange two vectors author: Jack Dongarra, Linpack author: P J Knight, CCFE, Culham Science Centre n : input integer : order of the matrices sx, sy sx(nincx) : input/output real array : first vector incx : input integer : interval in storage between sx array elements sy(nincy) : input/output real array : second vector incy : input integer : interval in storage between sy array elements This routine swaps the contents of two vectors, using unrolled loops for increments equal to 1. !

Arguments

Singular Value Decomposition author: P J Knight, CCFE, Culham Science Centre author: B. S. Garbow, Applied Mathematics Division, Argonne National Laboratory nm : input integer : Max number of rows of arrays a, u, v; >= m,n m : input integer : Actual number of rows of arrays a, u n : input integer : Number of columns of arrays a, u, and the order of v a(nm,n) : input/output real array : On input matrix to be decomposed; on output, either unchanged or overwritten with u or v w(n) : output real array : The n (non-negative) singular values of a (the diagonal elements of s); unordered. If an error exit is made, the singular values should be correct for indices ierr+1,ierr+2,...,n. matu : input logical : Set to .true. if the u matrix in the decomposition is desired, and to .false. otherwise. u(nm,n) : output real array : The matrix u (orthogonal column vectors) of the decomposition if matu has been set to .true., otherwise u is used as a temporary array. u may coincide with a. If an error exit is made, the columns of u corresponding to indices of correct singular values should be correct. matv : input logical : Set to .true. if the v matrix in the decomposition is desired, and to .false. otherwise. v(nm,n) : output real array : The matrix v (orthogonal) of the decomposition if matv has been set to .true., otherwise v is not referenced. v may also coincide with a if u is not needed. If an error exit is made, the columns of v corresponding to indices of correct singular values should be correct. ierr : output integer : zero for normal return, or k if the k-th singular value has not been determined after 30 iterations. rv1(n) : output real array : work array This subroutine is a translation of the algol procedure SVD, Num. Math. 14, 403-420(1970) by Golub and Reinsch, Handbook for Auto. Comp., vol II - Linear Algebra, 134-151(1971).

It determines the singular value decomposition a=usv^t of a real m by n rectangular matrix. Householder bidiagonalization and a variant of the QR algorithm are used. None

Arguments

Arguments

Arguments

Arguments

Arguments

Arguments

Arguments

Routine to calculate the inboard, outboard and total volumes of a toroidal shell comprising two elliptical sections author: P J Knight, CCFE, Culham Science Centre rshell : input real : major radius of centre of both ellipses (m) rmini : input real : horizontal distance from rshell to outer edge of inboard elliptical shell (m) rmino : input real : horizontal distance from rshell to inner edge of outboard elliptical shell (m) zminor : input real : vertical internal half-height of shell (m) drin : input real : horiz. thickness of inboard shell at midplane (m) drout : input real : horiz. thickness of outboard shell at midplane (m) dz : input real : vertical thickness of shell at top/bottom (m) vin : output real : volume of inboard section (m3) vout : output real : volume of outboard section (m3) vtot : output real : total volume of shell (m3) This routine calculates the volume of the inboard and outboard sections of a toroidal shell defined by two co-centred semi-ellipses. Each section's internal and external surfaces are in turn defined by two semi-ellipses. The volumes of each section are calculated as the difference in those of the volumes of revolution enclosed by their inner and outer surfaces.

See also eshellarea Internal CCFE note T&M/PKNIGHT/PROCESS/009, P J Knight: Surface Area and Volume Calculations for Toroidal Shells

Arguments

Routine to calculate the inboard, outboard and total volumes of a D-shaped toroidal shell author: P J Knight, CCFE, Culham Science Centre rmajor : input real : major radius to outer point of inboard straight section of shell (m) rminor : input real : horizontal internal width of shell (m) zminor : input real : vertical internal half-height of shell (m) drin : input real : horiz. thickness of inboard shell at midplane (m) drout : input real : horiz. thickness of outboard shell at midplane (m) dz : input real : vertical thickness of shell at top/bottom (m) vin : output real : volume of inboard straight section (m3) vout : output real : volume of outboard curved section (m3) vtot : output real : total volume of shell (m3) This routine calculates the volume of the inboard and outboard sections of a D-shaped toroidal shell defined by the above input parameters. The inboard section is assumed to be a cylinder of uniform thickness. The outboard section's internal and external surfaces are defined by two semi-ellipses, centred on the outer edge of the inboard section; its volume is calculated as the difference in those of the volumes of revolution enclosed by the two surfaces.

See also dshellarea Internal CCFE note T&M/PKNIGHT/PROCESS/009, P J Knight: Surface Area and Volume Calculations for Toroidal Shells

Arguments

Routine to calculate the inboard, outboard and total surface areas of a D-shaped toroidal shell author: P J Knight, CCFE, Culham Science Centre rmajor : input real : major radius of inboard straight section (m) rminor : input real : horizontal width of shell (m) zminor : input real : vertical half-height of shell (m) ain : output real : surface area of inboard straight section (m3) aout : output real : surface area of outboard curved section (m3) atot : output real : total surface area of shell (m3) This routine calculates the surface area of the inboard and outboard sections of a D-shaped toroidal shell defined by the above input parameters. The inboard section is assumed to be a cylinder. The outboard section is defined by a semi-ellipse, centred on the major radius of the inboard section.

See also dshellvol Internal CCFE note T&M/PKNIGHT/PROCESS/009, P J Knight: Surface Area and Volume Calculations for Toroidal Shells

Arguments

Routine to calculate the inboard, outboard and total surface areas of a toroidal shell comprising two elliptical sections author: P J Knight, CCFE, Culham Science Centre rshell : input real : major radius of centre of both ellipses (m) rmini : input real : horizontal distance from rshell to inboard elliptical shell (m) rmino : input real : horizontal distance from rshell to outboard elliptical shell (m) zminor : input real : vertical internal half-height of shell (m) ain : output real : surface area of inboard section (m3) aout : output real : surface area of outboard section (m3) atot : output real : total surface area of shell (m3) This routine calculates the surface area of the inboard and outboard sections of a toroidal shell defined by two co-centred semi-ellipses.

See also eshellvol Internal CCFE note T&M/PKNIGHT/PROCESS/009, P J Knight: Surface Area and Volume Calculations for Toroidal Shells

Arguments

Arguments

Arguments