SECTION 1.12 MISCELLANEOUS
1.12.1 RANDOM NUMBER GENERATION
There are some differences between the Microsoft, Lahey, WATCOM, and RS/6000 versions of this module. It is to be noted that the Lahey compiler has built-in random number generator RND(), RRAND(), and RANDS(x); the RS/6000 has built-in random number generator RAND() Microsoft Power Station and Lahey LF90 both have RANDOM_NUMBER.
Microsoft and Digital Visual Fortran: A call to FUNCTION RAND as in X=RAND(J) produces a uniformly distributed random variable on the (0,1) interval, where J is an integer between 1 and 5 (inclusive) which determines which "seed" is used. A congruential method is used. A call to FUNCTION RRAND, as in X=RRAND(SD,L), produces a uniformly distributed random variable on the (0,1) interval. If SD=0.0 in the call, the "seed" is determined with reference to the system clock. If SD is nonzero (a 10-digit number should be used), it will be the seed. L should be set equal to zero before the first call. A call to FUNCTION SRAND, as in X=SRAND(SD,L), produces a uniformly distributed random variable on the (0,1) interval. Everything is the same as in the case of RRAND except on the first call SRAND generates a table of 100 random numbers; on subsequent calls it selects randomly from this table and then fills in the used array element with a new random number. This tends to reduce serial correlation in the numbers generated. SRAND calls RRAND.
Lahey version: In the Lahey version RRAND is replaced by RRND and SRAND is replaced by SRND. SRND uses RRND. This version also contains SRND1, which is the same as SRND, except that it uses the Lahey-provided random number generator. See also details above. The call are RRND(SEED,L), SRND(SEED,L), SRND1(DUMMY). Note that the Lahey LF90 compiler contains a system subroutine called RANDOM_NUMBER(x).
WATCOM version: WATCOM contains the random number generator URAND(ISEED) where ISEED is INTEGER. GQOPT provides RRND, to be called as X=RRND(SEED,L) (identical with usage under Lahey). It also contains RAND(J). It also contains SRND1(ISEED), ISEED an integer. SRND1 uses URAND.
RS/6000 version: This version contains RRND, SRND, and SRND1. SRND calls RRND and SRND1 calls the IBM-provided random number generator. Call RRND as RRND(SEED,L) and SRND1 as SRND1().
X=RANDOM() generates a uniformly distributed number between 0 and 1 from a fixed seed, using a congruential multiplicative generator. To reset the seed, call the SUBROUTINE RANSET as in CALL RANSET(SEED), where SEED is a 10-digit REAL*8 integer (e.g., 1234567899.0). If SEED is .LE. 0 seed is ignored and the fixed seed is used.
FUNCTION GAUS(DUMMY) generates a normally distributed variable with mean zero and standard deviation 1.0. The Box-Muller transformation is used. GAUS should be used with caution, since it can produce very nonnormal looking samples. The Lahey version uses the Lahey random number generator RND, the RS/6000 version uses the IBM provided random number generator RAND. The WATCOM version uses the WATCOM generator URAND. The Microsoft version uses RAND (provided in GQOPT).
FUNCTION ZNORM computes a normally distributed variable by adding together N uniformly distributed variables on the unit interval and normalizing for mean and variance. The call is X=ZNORM(N) for the Lahey, PowerStation, and RS/6000 versions, ZNORM(N,ISEED) for the WATCOM version (this uses URAND).
FUNCTION XMB1 produces a normally distributed variable with mean zero and variance 1 by using the Marsaglia-Bray trans- formation. The call is X=XMB1() in the RS/6000 version, XMB1(DUMMY) in the Lahey version, XMB1(SEED) in the Microsoft version (since this uses RRAND), and XMB1(ISEED) in the WATCOM version. In the Microsoft version the call to XMB1 is XMB1(ISEED), where ISEED is INTEGER*2.
SUBROUTINE MVARNORM(COV,XMU,Y,P,NDIM,NREP,IER) produces jointly normally distributed variables.
The MAIN program should include
COMMON/BPRINT/IPT,NFILE,NDIG,NPUNCH,JPT,MFILE
CALL DFLT
The arguments are:
COV = a DOUBLE PRECISION array dimensioned COV(NDIM,NDIM)
should contain the desired covariancve matrix
not checked for symmetry) (input)
XMU = a DOUBLE PRECISION array dimensioned XMU(NDIM) should contain
the desired mean vector (input)
Y = a DOUBLE PRECISION array dimensioned Y(NREP,NDIM) will
contain NREP rows, each of which is a drawing from the NDIM-
dimensional normal with require mean and covariance matrix
P = a DOUBLE PRECISION work array dimensioned P(NDIM,NDIM)
NDIM = the dimensionality of the normal (input); must be <= 100
NREP = the number of replicates (input)
IER = error return (output); = 0 normal return, =-3 NDIM exceeds
100, =-67 the eigenvalues of the covariance matrix are not
all positive
SUBROUTINE EXPON(A,X,IER). Each call produces an exponentially distributed variate x, with the density function
The MAIN program should include
COMMON/BPRINT/IPT,NFILE,NDIG,NPUNCH,JPT,MFILE
CALL DFLT
If the subroutine is called with a <= 0, IER is set to -3 (input error).
SUBROUTINE CAUCHY(S,X,IER). Each call produces a Cauchy-distributed variate x, with the density function
The MAIN program should include
COMMON/BPRINT/IPT,NFILE,NDIG,NPUNCH,JPT,MFILE
CALL DFLT
If the subroutine is called with s <= 0, IER is set to -3 (input error).
SUBROUTINE RAYLEI(S,X,IER). Each call produces a Rayleigh-distributed variate x, with the density function
The MAIN program should include
COMMON/BPRINT/IPT,NFILE,NDIG,NPUNCH,JPT,MFILE
CALL DFLT If the subroutine is called with s <= 0, IER is set to -3 (input error).
The MAIN program should include
COMMON/BPRINT/IPT,NFILE,NDIG,NPUNCH,JPT,MFILE
CALL DFLT If the subroutine is called with a <= 0 or b <= 0, IER is set to -3 (input error).
1.12.2 MATRIX INVERSION (VERS, GJR, MATVRT, GAJORE)
A. To invert a symmetric matrix A, issue the following call
CALL GJR(A,N,EPS,IER,DETER)
where A = N x N matrix to be inverted. It must be dimensioned 30 x 30 in calling program. If you wish to change this, change the dimensioning of A and arrays IP and JO in GJR as you wish. A is destroyed and will contain the inverse.
N = size of matrix to be inverted
EPS = a tolerance level (default 1.D-15)
IER = 0 inversion successful, otherwise singular.
DETER = value of the determinant
B.To invert a matrix that need not be symmetric, use SUBROUTINE VERS or MATVRT. The calling sequence for VERS is
CALL VERS(A,Y,N,NDIM,IER)
where
A = matrix to be inverted (destroyed in the process),
dimensioned NDIM by NDIM
Y = the inverse, dimensioned NDIM by NDIM
N = the order of A; must be less than or equal to NDIM and less
than or equal to 200
NDIM = the dimensions of A and Y (i.e., A(NDIM,NDIM))
IER = 0 normal return
IER = -1 input error
IER = -2 singularity .
The calling sequence for MATVRT is
CALL MATVRT(A,B,C,N,EPSV,ISING,DETER)
where
A = matrix to be inverted, dimensioned N by N (input and
output); input matrix A is destroyed in the process
B = a vector dimensioned N (input)
C = a vector dimensioned N (input)
EPSV = a pivot tolerance (set it at 1.D-15) (input)
ISING = 0 if matrix is successfully inverted,
ISING = 1 if singular (output)
DETER = the determinant (output)
C. To solve a system of simultaneous linear equations by Gauss-Jordan reduction, use GAJORE. Calling sequence is
CALL GAJORE(A,B,DETT,N,EPSSS,IER)
where
A = matrix of system of equations. dimensioned N by N (input)
B = vector of right hand sides, dimensioned N (input)
DETT = determinant (output)
N = size of matrix (input)
EPSSS = pivot tolerance (input), 1.D-15
IER = 0 if solution is successful,
IER = 1 if not (output).
See test program MATTST for a comparison of matrix inversion algorithms inverting the Hilbert matrix and some related matrices.
1.12.3 EIGENVALUE CALCULATION (MATEV2)
To find the eigenvalues and eigenvectors of a real symmetric matrix, issue the call
CALL MATEV2(A,B,N,IA,IB,ISW)
where
A = is a matrix dimensioned (IA, ) containing the matrix of
which the eigenvalues are to be found
B = is a matrix dimensioned (IB, )
N = is the order of the upper left N by N corner of A of which
the eigenvalues are to be found
ISW = 0 only the eigenvalues will be found,
ISW = 1 the eigenvectors will be found also
The eigenvalues will be placed in the diagonal positions of A in descend- ing order. Eigenvectors are placed in B by columns in corresponding order. For eigenvalues that are close, the eigenvectors are subject to error but will be orthogonal. Eigenvalues are found to an accuracy that depends on the parameter EPS1 (see listing of program).
1.12.4 COMPUTATION OF DETERMINANTS
To compute the determinant of a square array, the following call should be issued:
CALL DTRM(A,N,EPS,ISING,DETER)
where
A = name of array
N = number of rows (or columns) in the determinant array
EPS = a tolerance level (e.g.,1.D-15)
ISING = 0 if calculation is successful,
ISING =1 otherwise
DETER = the value of the determinant computed by DTRM
It should ne noted that A is a single dimensional array in DTRM. If in the user's calling program A is doubly dimensioned and if the dimensioned size of A is identical with the order N of the determinant array, nothing special needs to be done. EXAMPLE: User's calling program:
DIMENSION A(10,10)
. . .
EPS=1.D-15
N=10
CALL DTRM(A,N,EPS,ISING,DETER)
If the user's matrix containing the array elements is dimensioned larger than N, (i.e., if the relevant array is contained in the upper left hand corner of the matrix A, the matrix should first by transferred to a single dimensional array or to another array that is dimensioned N x N. EXAMPLE:
DIMENSION A(12,12),B100)
C statements here to fill the upper left 10x10 portion of A with numbers
EPS=1.D-15
K=1
DO J=1,10
DO I=1,10
B(K)=A(I,J)
K=K+1
END DO
END DO
CALL DTRM(B,10,EPS,ISING,DETER)
(See also Section 1.12.2)
Using the Microsoft Compilers or the WATCOM Compiler: The timing of computations (elapsed time) is accomplished by the SUBROUTINE ELTIM(TIME). ELTIM provides in the REAL*8 variable TIME the elapsed time, in seconds, since the last invocation of ELTIM. The first invocation of ELTIM in a program simply records the current state of the system clock. Typical usage is illustrated below:
CALL ELTIM(TIME1)
C . . . program steps to be timed
CALL ELTIM(TIME2)
WRITE (*,1000) TIME2
1000 FORMAT(' Time = ',f10.4)
Alternatively, for the WATCOM Compiler one may include the files TBEGWAT.FOR at the beginning and TENDWAT.FOR at the end of the program segment to be timed. These files are in the \MAIN subdirectory of the Installation diskettes. Using the Lahey Compilers: Two short Fortran files, TBEGIN.FOR and TEND.FOR are provided in the \MAIN subdirectory (where the test programs reside). Include the first of these just before the code segment to be timed and the second just after the code segment to be timed.
Five sort routines are provided. SORT1(A,N) sorts the first N items of the REAL*8 array A into ascending order. SORT1A does the same for a REAL*4 array A. SORT2(A,B,N) sorts the first N items of the REAL*4 array A into ascending order and performs the corresponding permutation on the REAL*4 array B. SORT3A(A,IP,N) sorts the first N items of a REAL*4 array A and performs the corresponding permutation on the INTEGER*4 array IP. For details, see Section 15.2.
The same program is called MOPTL for the Lahey, Microsoft, and WATCOM versions of GQOPT and is called MAKEOPT for the Microsoft Power Station and RS/6000 versions. In what follows, we refer only to MAKEOPT, but all remarks pertain to MOPTL as well. The MAKEOPT utility makes it easier to write a GQOPT optimization program. Optimization programs characteristically share anywhere from 40 to 100 program steps, including DIMENSION and COMMON statements, other statements that initialize certain variables required by GQOPT, and so on. MAKEOPT is a utility that asks the user a certain number of questions about the program to be written and then generates the appropriate FORTRAN code. MAKEOPT is straigthforward to use and is started by typing MAKEOPT (assuming that you are in the appropriate directory). The only nonobvious question that is addressed by MAKEOPT to the user is "Create standard input file?" The objective of this is to have the optimization program read certain parameters out of a file rather than have these parameters be set in the program itself, which would necessitate recompiling every time that these parameters are changed. Obviously, the standard input file created by MAKEOPT can then also be used by the user for other parameters that he or she finds useful. MAKEOPT is written in FORTRAN (except for the RS/6000 version, in which it is written in Kornshell) and the appropriate files are MAKEOPT.EXE and MOPTL.EXE.
1.12.8 COMPILING AND LIBRARY ROUTINES
For each PC version of GQOPT we provide batch programs to (1) compile all of GQOPT, (2) make a library file out of the compiled code.
Lahey LF95: To
compile GQOPT: CMPGQ710.BAT To make a library: MKLIB710.BAT
WATCOM: To compile GQOPT:
CMWAT710.BAT To make a library: MKLIB710.BAT
DVF: To compile GQOPT:
CMDVF710.BAT To make a library: MLIB710.BAT
All versions are easily modified with respect to the way they work and with respect to what compiler switches are turned on.
1.12.9 ENCRYPTION AND DECRYPTION (ENC, DEC. CRYPT)
ENC, DEC and CRYPT are programs that encrypt, respectively decrypt, ASCII files containing no more than 80 columns and arbitrary number of rows. For added security, no source program is provided for these routines. Some versions of GQOPT contain ENC and DEC to encrypt and decrypt, while some others use CRYPT for both purposes. In the latter case the user is prompted for which operation is desired. To encrypt a file, type ENC. ENC will prompt you for the name of the file to be encrypted, and then will prompt you (twice) for a keyword to be used in encryption. The keyword may not exceed 8 characters and any printable characters may be used. It is strongly recommended that at least 7 characters be used in the keyword and that no dictionary words or names be employed. The output of encryption is placed in a file called ENCRP.TXT and the original file remains unaltered. You need to erase the input file yourself. To decrypt a file, type DEC and you will be prompted for the name of the file to be decrypted and the keyword to be used. The output will be placed in a file called DECRP.TXT and the original encrypted file is left unaltered. Alternately, if your version contains CRYPT, just type CRYPT for either encryption or decryption. For added security, no source program is provided for ENC, DEC, and CRYPT. Note: These programs are not available for purchasers outside the United States.
1.12.10 ASCII TO INTEGER AND INTEGER TO ASCII CONVERSION (ASCINT, INTASC)
It is frequently desirable to convert integers (Type INTEGER) to their ASCII equivalent formats or to convert numbers entered in ASCII to their INTEGER form. To do the former, use
CALL INTASC(N,ANS)
where N is INTEGER*4 and is the input quantity, and where ANS should be CHARACTER*1(10) and will contain the ASCII equivalent. To convert from ASCII to INTEGER, use
CALL ASCINT(REPLY,N,ICV)
where REPLY should be CHARACTER*1(10) and is the input, N is INTEGER*4 and is the output. ICV is an INTEGER error code with the following meaning:
ICV = 0 Normal return and conversion
ICV = 1 Leading or embedded blank found in REPLY
ICV = 2 A character other than a digit was found
ICV = 3 Number is too large to be converted (>2147483647)
1.12.11 BIVARIATE STANDARD NORMAL DENSITY (BINORM)
To calculate the standard normal bivariate density (unit variances, zero means, and zero correlation) use the FUNCTION BINORM(X,Y). Thus,
DENS = BINORM(X,Y)
sets DENS equal to the density at X,Y.
1.12.12 COMPUTATION OF THE ROOTS OF POLYNOMIALS (PROOTS)
The calling sequence is
CALL PROOTS(XCOF,COF,M,ROOTR,ROOTI,ITERL,IER)
where
XCOF = The coefficients of the polynomial ordered from
smallest to largest power (for x3
+ 2x2 + 3x + 4=0, XCOF(1)=4,
XCOF(2)=3, XCOF(3)=2, XCOF(4)=1). To be dimensioned M+1. (input)
COF = Working array to be dimensioned M+1.
M = degree of the
polynomial (3 in the previous example)(input)
ROOTR = Real parts
of the roots, dimensioned M (output)
ROOTI = Imaginary parts of the roots, dimensioned M (output)
ITERL = Number of permitted iterations (input)
IER = 0 successful solution,
IER = 1 no solution (output)