R16LIB
A Quadruple Precision Real Arithmetic Utility Library

R16LIB is a FORTRAN90 library which contains a number of utility routines using "R16" or "quadruple precision real" arithmetic.

At the moment, the status of quadruple precision real data in FORTRAN is somewhat vague. This library is very small, and right now is simply a bag containing the occasional results of forays into the world of quadruple real precision.

Issues include:

• compiler support: does the compiler allow quadruple real precision at all?
• declarations: how do you declare that a given value has type "quadruple real precision"? My choice has been "real ( kind = 16 )".
• conversion: how do you convert, for example, an integer value to quadruple real? My choice has been "x = real ( n, kind = 16 )".
• constants: how do you specify a literal decimal constant? If you are not careful, FORTRAN will treat a literal decimal constant as a single precision value, dropping all those extra digits you so carefully typed in. (Check this sometime!) My choice has been to follow the numeric value by the suffix "_16"?
• intrinsics: how do you call for the square root, cosine and so on. I assume that the generic intrinsic functions are available, so "y = sqrt ( x )" is presumed to carry out a quadruple precision real calculation if that is the type of x.
• non-intrinsics: Some FORTRAN non-intrinsic functions may be common enough that the user expects them to be available, even in oddball precisions. The GAMMA function, for instance, is not part of the standard, but is often supplied for real and double precision arguments. Is it available for quadruple precision?
• efficiency: R16 arithmetic is often done slowly, through software rather than in hardware. This can result in an enormous increase in computation time;
• complex values: many of the above questions must be repeated if the user is interested in quadruple precision complex arithmetic. It would seem natural to assume, for instance, that the corresponding declaration might be "complex ( kind = 16 )" but I have seen one or two cases in which it was asserted that the correct declaration would be "complex ( kind = 32 )" because one needs to count the storage for the real and imaginary parts.

Licensing:

The computer code and data files described and made available on this web page are distributed under the GNU LGPL license.

Languages:

R16LIB is available in a FORTRAN90 version.

Related Data and Programs:

C4LIB, a FORTRAN90 library which implements certain elementary functions for "C4" or single precision complex variables;

C8LIB, a FORTRAN90 library which implements certain elementary functions for "C8" or double precision complex variables;

I4LIB, a FORTRAN90 library which contains many utility routines, using "I4" or "single precision integer" arithmetic.

I8LIB, a FORTRAN90 library which contains many utility routines, using "I8" or "double precision integer" arithmetic.

R16_HERMITE_RULE, a FORTRAN90 program which can compute and print a Gauss-Hermite quadrature rule, using "quadruple precision real" arithmetic.

R16_INT_EXACTNESS_GEN_HERMITE, a FORTRAN90 program which tests the polynomial exactness of generalized Gauss-Hermite quadrature rules, using "quadruple precision real" arithmetic.

R4LIB, a FORTRAN90 library which contains many utility routines, using "R4" or "single precision real" arithmetic.

R8LIB, a FORTRAN90 library which contains many utility routines, using "R8" or "double precision real" arithmetic.

SUBPAK, a FORTRAN90 library which contains many utility routines;

Source Code:

• r16lib.f90, the source code;
• r16lib.sh, commands to compile the source code;

Examples and Tests:

• r16lib_test.f90, a sample calling program;
• r16lib_test.sh, commands to compile, link and run the sample calling program;
• r16lib_test.txt, the output file.

List of Routines:

• I4_LOG_10 returns the integer part of the logarithm base 10 of an I4.
• R16_ABS returns the absolute value of an R16.
• R16_ATAN computes the inverse tangent of the ratio Y / X.
• R16_CAS returns the "casine" of an R16.
• R16_CEILING rounds an R16 "up" (towards +oo) to the next I4.
• R16_CHOOSE computes the binomial coefficient C(N,K) as an R16.
• R16_CHOP chops an R16 to a given number of binary places.
• R16_CSQRT returns the complex square root of an R16.
• R16_CUBE_ROOT returns the cube root of an R16.
• R16_DIFF computes the difference of two R16's to a specified accuracy.
• R16_DIGIT returns a particular decimal digit of an R16.
• R16_EPSILON returns the R16 roundoff unit.
• R16_EXP computes the exponential of an R16, avoiding overflow and underflow.
• R16_FACTORIAL computes the factorial of N.
• R16_FACTORIAL2 computes the double factorial function.
• R16_FLOOR rounds an R16 "down" (towards -oo) to the next integer.
• R16_FRACTION returns the fraction part of an R16.
• R16_HUGE returns a very large R16.
• R16_IN_01 is TRUE if an R16 is in the range [0,1].
• R16_IS_INT determines if an R16 represents an integer value.
• R16_LOG_2 returns the logarithm base 2 of an R16.
• R16_LOG_10 returns the logarithm base 10 of an R16.
• R16_LOG_B returns the logarithm base B of an R8.
• R16_MANT computes the "mantissa" or "fraction part" of an R16.
• R16_MOD returns the remainder of R16 division.
• R16_MODP returns the nonnegative remainder of R16 division.
• R16_MOP returns the I-th power of -1 as an R16 value.
• R16_NINT returns the nearest integer to an R16.
• R16_NORMAL returns a scaled pseudonormal R16.
• R16_NORMAL_01 returns a unit pseudonormal R16.
• R16_PI returns the value of pi as an R16.
• R16_POWER computes the P-th power of an R16.
• R16_POWER_FAST computes an integer power of an R16.
• R16_PYTHAG computes sqrt ( A * A + B * B ), avoiding overflow and underflow.
• R16_ROUND2 rounds an R16 to a specified number of binary digits.
• R16_ROUNDB rounds an R16 to a given number of digits in a given base.
• R16_ROUNDX rounds an R16.
• R16_SIGN returns the sign of an R16.
• R16_SIGN_OPPOSITE is TRUE if two R16's are not of the same sign.
• R16_SIGN_OPPOSITE_STRICT is TRUE if two R16's are strictly of opposite sign.
• R16_SWAP swaps two R16's.
• R16_SWAP3 swaps three R16's.
• R16_TINY returns a very small but positive R16.
• R16_TO_R16_DISCRETE maps R to RD in [RMIN, RMAX] with NR possible values.
• R16_TO_DHMS converts decimal days into days, hours, minutes, seconds.
• R16_TO_I4 maps X in [XMIN, XMAX] to integer IX in [IXMIN, IXMAX].
• R16_UNIFORM returns a scaled pseudorandom R16.
• R16_UNIFORM_01 returns a unit pseudorandom R16.
• R16_UNSWAP3 unswaps three R16's.
• R16_WALSH_1D evaluates the Walsh function.
• R16MAT_BORDER_ADD adds a "border" to an R16MAT.
• R16MAT_BORDER_CUT cuts the "border" of an R16MAT.
• R16MAT_CHOLESKY_FACTOR computes the Cholesky factor of a symmetric matrix.
• R16MAT_CHOLESKY_SOLVE solves a Cholesky factored linear system A * x = b.
• R16MAT_CHORESKY_FACTOR computes the "Choresky" factor of a symmetric matrix.
• R16MAT_COPY copies an R16MAT.
• R16MAT_DET computes the determinant of an R16MAT.
• R16MAT_DET_2D computes the determinant of a 2 by 2 R16MAT.
• R16MAT_DET_3D computes the determinant of a 3 by 3 R16MAT.
• R16MAT_DET_4D computes the determinant of a 4 by 4 R16MAT.
• R16MAT_DET_5D computes the determinant of a 5 by 5 R16MAT.
• R16MAT_DIAG_ADD_SCALAR adds a scalar to the diagonal of an R16MAT.
• R16MAT_DIAG_ADD_VECTOR adds a vector to the diagonal of an R16MAT.
• R16MAT_DIAG_GET_VECTOR gets the value of the diagonal of an R16MAT.
• R16MAT_DIAG_SET_SCALAR sets the diagonal of an R16MAT to a scalar value.
• R16MAT_DIAG_SET_VECTOR sets the diagonal of an R16MAT to a vector.
• R16MAT_EXPAND_LINEAR linearly interpolates new data into an R16MAT.
• R16MAT_EXPAND_LINEAR2 expands an R16MAT by linear interpolation.
• R16MAT_GIVENS_POST computes the Givens postmultiplier rotation matrix.
• R16MAT_GIVENS_PRE computes the Givens premultiplier rotation matrix.
• R16MAT_HESS approximates a Hessian matrix via finite differences.
• R16MAT_HOUSE_AXH computes A*H where H is a compact Householder matrix.
• R16MAT_HOUSE_FORM constructs a Householder matrix from its compact form.
• R16MAT_HOUSE_HXA computes H*A where H is a compact Householder matrix.
• R16MAT_HOUSE_POST computes a Householder post-multiplier matrix.
• R16MAT_HOUSE_PRE computes a Householder pre-multiplier matrix.
• R16MAT_IDENTITY stores the identity matrix in an R16MAT.
• R16MAT_IN_01 is TRUE if the entries of an R16MAT are in the range [0,1].
• R16MAT_INDICATOR sets up an "indicator" R16MAT.
• R16MAT_INVERSE_2D inverts a 2 by 2 R16MAT using Cramer's rule.
• R16MAT_INVERSE_3D inverts a 3 by 3 R16MAT using Cramer's rule.
• R16MAT_INVERSE_4D inverts a 4 by 4 R16MAT using Cramer's rule.
• R16MAT_JAC estimates a dense jacobian matrix of the function FX.
• R16MAT_L_INVERSE inverts a lower triangular R16MAT.
• R16MAT_L_PRINT prints a lower triangular R16MAT.
• R16MAT_L_SOLVE solves a lower triangular linear system.
• R16MAT_L1_INVERSE inverts a unit lower triangular R16MAT.
• R16MAT_LT_SOLVE solves a transposed lower triangular linear system.
• R16MAT_LU computes the LU factorization of a rectangular R16MAT.
• R16MAT_MAX returns the maximum entry of an R16MAT.
• R16MAT_MAX_INDEX returns the location of the maximum entry of an R16MAT.
• R16MAT_MAXCOL_MINROW gets the maximum column minimum row of an M by N R16MAT.
• R16MAT_MAXROW_MINCOL gets the maximum row minimum column of an M by N R16MAT.
• R16MAT_MIN returns the minimum entry of an M by N R16MAT.
• R16MAT_MIN_INDEX returns the location of the minimum entry of an R16MAT.
• R16MAT_MINCOL_MAXROW gets the minimum column maximum row of an M by N R16MAT.
• R16MAT_MINROW_MAXCOL gets the minimum row maximum column of an M by N R16MAT.
• R16MAT_MM multiplies two R16MAT's.
• R16MAT_MTV multiplies a transposed matrix times a vector
• R16MAT_MV multiplies a matrix times a vector.
• R16MAT_NINT rounds the entries of an R16MAT.
• R16MAT_NORM_EIS returns the EISPACK norm of an R16MAT.
• R16MAT_NORM_FRO returns the Frobenius norm of an R16MAT.
• R16MAT_NORM_L1 returns the matrix L1 norm of an R16MAT.
• R16MAT_NORM_L2 returns the matrix L2 norm of an R16MAT.
• R16MAT_NORM_LI returns the matrix L-oo norm of an R16MAT.
• R16MAT_NULLSPACE computes the nullspace of a matrix.
• R16MAT_NULLSPACE_SIZE computes the size of the nullspace of a matrix.
• R16MAT_ORTH_UNIFORM returns a random orthogonal R16MAT.
• R16MAT_PLOT "plots" an R16MAT, with an optional title.
• R16MAT_PLOT_SYMBOL returns a symbol for an element of an R16MAT.
• R16MAT_POLY_CHAR computes the characteristic polynomial of an R16MAT.
• R16MAT_POWER computes a nonnegative power of an R16MAT.
• R16MAT_POWER_METHOD applies the power method to an R16MAT.
• R16MAT_PRINT prints an R16MAT.
• R16MAT_PRINT_SOME prints some of an R16MAT.
• R16MAT_PRINT2 prints an R16MAT.
• R16MAT_REF computes the row echelon form of a matrix.
• R16MAT_RREF computes the reduced row echelon form of a matrix.
• R16MAT_SOLVE uses Gauss-Jordan elimination to solve an N by N linear system.
• R16MAT_SOLVE_2D solves a 2 by 2 linear system using Cramer's rule.
• R16MAT_SOLVE_3D solves a 3 by 3 linear system using Cramer's rule.
• R16MAT_SOLVE2 computes the solution of an N by N linear system.
• R16MAT_SYMM_EIGEN returns a symmetric matrix with given eigensystem.
• R16MAT_SYMM_JACOBI applies Jacobi eigenvalue iteration to a symmetric matrix.
• R16MAT_TO_R16PLU factors a general R16MAT.
• R16MAT_TRACE computes the trace of an R16MAT.
• R16MAT_TRANSPOSE_IN_PLACE transposes a square matrix in place.
• R16MAT_TRANSPOSE_PRINT prints an R16MAT, transposed.
• R16MAT_TRANSPOSE_PRINT_SOME prints some of an R16MAT, transposed.
• R16MAT_U_INVERSE inverts an upper triangular R16MAT.
• R16MAT_U1_INVERSE inverts a unit upper triangular R16MAT.
• R16MAT_UNIFORM fills an R16MAT with scaled pseudorandom numbers.
• R16MAT_UNIFORM_01 fills an R16MAT with unit pseudorandom numbers.
• R16MAT_VAND2 returns the N by N row Vandermonde matrix A.
• R16MAT_ZERO zeroes an R16MAT.
• R16VEC_HEAP_A reorders an R16VEC into an ascending heap.
• R16VEC_HEAP_D reorders an R16VEC into an descending heap.
• R16VEC_PRINT prints an R16VEC.
• R16VEC_PRINT_PART prints "part" of an R16VEC.
• R16VEC_PRINT_SOME prints "some" of an R16VEC.
• R16VEC_UNIFORM_01 returns a unit pseudorandom R16VEC.
• TIMESTAMP prints the current YMDHMS date as a time stamp.

You can go up one level to the FORTRAN90 source codes.