NO NV Remarks: None.

- Number of observations - Number of variables for each observation
Subroutines and function subprograms required: None. Method: Each row (observation) of matrix A with corresponding non-zero element in S vector is tested. Observations are compared with specified lower and upper variable bounds and a count is kept in vectors under, between, and over.
SUBROUTINE BOUN01A.SOLO.BHI.UNDER.BETWOVER,NOOVI DIMENSION A111.51 1/0L0111.041111,UNDER1/1.8E7W11110VER111 CLEAR OUTPUT VECTORS. DO I Kwl.NV UNDER(10.0.0 BETWIK/A0.OVERIAI.D.0 TEST SUBSET VECTOR 00 B J.1,NO IJJ-410 IFISIJII 2,8.2 COMPARE OBSERVATIONS WITH BOUNDS 2 DO 7 1.1,88 IJAIJANO 1E141121-4%01111 5.3.IFTAIIJI-BH11111 4,4,6 COUNT 4 BETW111A8ETW11141.0 GO TO 7 S UNDER111UNDER117.1.0 GO TO OVERIII=OVER11141.CONTINUE B CONTINUE RETURN END
BOUND I BOUND 2 BOUND 3 BOUND 4 ROUND 5 ROUND S BOUND 7 ROUND B BOUND 9 ROUND 10 ROUND II BOUND 12 BOUND 13 BOUND 14 BOUND IS ROUND 16 BOUND 17 BO U ND BOUND 19 BO U ND 20

ROUND BOUND

22 BOUND 73 ROUND 24 ROUND 75 BOUND 76
SUBST Purpose: Derive a subset vector indicating which observations in a set have satisfied certain conditions on the variables. Usage: CALL SUBST (A, C, R, B, S, NO, NV, NC) Parameter B must be defined by an external statement in the calling program. Description of parameters: A - Observation matrix, NO by NV C - Input matrix, 3 by NC, of conditions to be considered. The first element of each column of C represents the number of the variable (column of the matrix A) to be tested, the second element of each column is a relational code as follows: 1. for less than 2. for less than or equal to 3. for equal to 4. for not equal to 5. for greater than or equal to 6. for greater than The third element of each column is a quantity to be used for comparison with the observation values. For example, the following column in C: 2. 5. 92.5 causes the second variable to be tested for greater than or equal to 92.5. R - Working vector used to store intermediate results of above tests on a single observation. If condition is satisfied, R(I) is set to 1. If it is not, R(I) is set to 0. Vector length is NC. B - Name of subroutine to be supplied by the user. It consists of a Boolean expression linking the intermediate values stored in vector R. The Boolean operators are ,4" for 'and', '+' for 'or ,. Example: SUBROUTINE BOOL (R, T) DIMENSION R(3) T= R(1)*(R,(2)+ R(3)) RETURN END The above expression is tested for R(1). AND. (R(2). OR. R(3))

where F. = frequency count in i-th interval. Then, the following are computed: First Moment (Mean): E F. [UBO 1 + (i-0. 5) UBO2] i=1
Subroutines and function subprograms required: None. Method: If S(I) contains a non-zero code, I-th observation is copied from the input matrix to the output matrix.
SUBROUTINE SuRmx 17,0,501,Nv,41 UIMENSION A1'10111,5111 1.4n LL41 OJ4I,Nv no 15 141,40 371.41 IFISIIII 14, 15, 10 IC LL=LI*1 DILL144ILI 15 CoNTINor 20 CONTINur COUNT NON-7271 CODES Is vrcT94 S N10 Ixi.NO 121511 I1 in, 30, 4741.CnNTINuF SUBMSUPMX SU pp X 50788

SUBmx SURmx SUBmx

I j-th Moment (Variance):

511280

ANS. - i= 1
(4) n \- F. [UBO + (i-0. 5) UBO2 - ANSI] 1 w-T

2, 3, 4

summx 13 SUPBX 14 sommx 15 SUpmx IS SUB m x /4 SuPmx 19 SUP 7 K 20
These moments are biased and not corrected for grouping
Subroutine. MOMEN Purpose: To find the first four moments for grouped data on equal class intervals. Usage: CALL MOMEN (F, UBO, NOP, ANS) Description of Parameters: - Grouped data (frequencies). Given as a F vector of length (UB0(3)-UB0(1))/ UBO(2) UBO - 3 cell vector, UBO(1) is lower bound and UBO(3) upper bound on data. UBO(2) is class interval. Note that UBO(3) must be greater than UB0(1). Option parameter. If NOP = 1, NOP ANS(1) = MEAN. If NOP = 2, ANS(2) = second moment. If NOP = 3, ANS(3) = third moment. If NOP = 4, ANS(4) = fourth moment. If NOP = 5, all four moments are filled in. ANS - Output vector of length 4 into which moments are put. Remarks: Note that the first moment is not central but the value of the mean itself. The mean is always calculated. Moments are biased and not corrected for grouping. Subroutines and function subprograms required: None. Method: Refer to M. G. Kendall, 'The Advanced Theory of Statistics', V.1, Hafner Publishing Company, 1958, Chapter 3.
SUBROUTINE MOMEN IF.0800401,4ANS/ DIMENSION F111.1J80131.ANS141 DO 100 1.1.4 ANSI 10.0 CALCULATE THE NUMBER OF CLASS INTERVALS NAIUB0131-UBOIIWUB0121+0.5 CALCULATE TOTAL FREQUENCY T=0.110 1=10 T=T4+1I1 IFINOR-51 lin. 12^, 115 NOR.5 JUMP.' GO ro 150 JU4P=2 FIRST 40mENT 00 Ino 1=1.5 RUMEN 1 MOMENM01 MORES 3 i4OMEN 4 MORES 5 MOMENMO2 MOMEN 7 HOMES R
TTSTT This subroutine computes certain t-statistics ,on the means of populations under various hypotheses. The sample means of A 1 , A2 ,. , ANA and B1, B2 ,. , BNB are normally found by the following formulas: NB NA E Bi E Ai i=1 NA

FIND NEAT LARIEST RANK 1.41.N IFIR111-41 30.30.10 in IFIRIIIRI 20.10.X=RIII INOAINCONTINUE IF ALL RANKS AAVE BEEN TESTES, RETURN IFIINO/ 90.90.40

40 Y=X

C C 80 90
CT=0.0 COUNT TIES Do 60 1=1,N IFIRIIIA) 60.50,60 CT=CT41.0 CONTINUE CALCULATE CORRECTION FACTOR IFICT1 70,500 IPIKT-1/ 75,00,75 TAT4CTAICT-1.1/2.0 GO TO 5 TT;,,, I CTACTACTCTI/12.0 GO RETURN ENO
TIE TIE TIE TIE TIE TIF TIE TIP TIE TIE TIE TIE TIE TIE TIE TIE TIE TIF TIF TIF TIE TIC TIE TIE TIE TIE TIE TIE TIF TIE TIE TIE

I 9 IO II IA IT 30 II 12

Statistics Nonparametric 59
Statistics - Random Number Generators RANDU Purpose: Computes uniformly distributed random floating point numbers between 0 and 1.0 and integers in the range 0 to 2**15. Usage: CALL RANDU (IX, IY, YFL) Description of parameters: IX - For the first entry this must contain any odd positive integer less than 32,768. After the first entry, IX should be the previous value of IY computed by this subroutine. IY - A resultant integer random number required for the next entry to this subroutine. The range of this number is from zero to 2**15. YFL - The resultant uniformly distributed, floating point, random number in the range 0 to 1.0. Remarks: This subroutine is specific to the IBM 1130. This subroutine should not repeat its cycle in less than 2 to the 13th entries. Note: If random bits are needed, the high order bits of IY should be chosen. Subroutines and function subprograms required: None. Method: Power residue method discussed in IBM manual Random Number Generation and Testing (C208011).
SUBROUTINE RANDUIIX,IY,YRL/ I7=IX,R9R [F11715.6,17=7B1.I7 YFL4YFL/32767. RETURN END RANDU I RANDU 2 RANDU 3 RANDII 4 RANDU RANDU 6 RANDU 7 RANDU
GAUSS This subroutine computes a normally distributed random number with a given mean and standard deviation. An approximation to normally distributed random numbers Y can be found from a sequence of uniform random numbers* using the formula:
(1) K/1 where X. is a uniformly distributed random number, 0< X. < 1 K is the number of values X. to be used Y approaches a true normal distribution asympototically as K approaches infinity. For this subroutine, K was chosen as 12 to reduce execution time. Equation (1) thus becomes: 12 Y = E - 6.0 i= 1 The adjustment for the required mean and standard deviation is then Y * S + AM (2)

C SUBROUTINE TPROIA,B,R,4,M,MSA,MSB,L1 DIMENSION 11111,81110111 SPECIAL CASE TOR DIAGONAL BY DIAGONAL MS=4SA10+mse ITIMS-221 30,10,30 on 70 11,N RIII.AIIIABIII RETURN MULTIPLY TRANSPOSE OF A BY B /R=I on 90 K.1,L J=1,M RIIR1.0.80 11,N IFIMSI 40,60,40 CALL LEICII,J,1401,M,MSAI CALL LOCII,K,I3,B,L,MS81 IFIIA1 50,80,50 1E1181 70,80,70 IA./S*13-11m1 18.6114-1161 RIIRI.R1/01.41/41.81101 CONTINUE IR.IRml RETURN END TPRO TPRD TPRD TPRD TPRO TPRO TPRO TURD TPRD TPRD TPRO TPRO TPRD TPRD TURD TURD TURD TROD 7/110 TURD TURD TURD TURD TURD TPRO TPRO 5 6

17 /17 IR 19 2n 25 26

MATA Purpose: Premultiply a matrix by its transpose to form a symmetric matrix. Usage: CALL MATA(A,R,N,M, MS) Description of parameters: A - Name of input matrix. R - Name of output matrix. N - Number of rows in A. M - Number of columns in A. Also number of rows and number of columns of R. MS - One digit number for storage mode of matrix A: 0 - General. 1 - Symmetric. 2 - Diagonal. Remarks: Matrix R cannot be in the same location as matrix A. Matrix R is always a symmetric matrix with a storage mode=1. Subroutines and function subprograms required: LOC Method: Calculation of (A transpose A) results in a symmetric matrix regardless of the storage mode of the input matrix. The elements of matrix A are not changed.
SUR4OUTINE MATAIA,R.N.M,MSI DIMENSION 4I11.411/ K.101 OXPIKPK-K//2 DJ=I.M IFIJKI 10.10,IRJJ+KX R114/40 DO 60 1.1.N 1E1451 20.40.CALL LOCII.J.IA.V.m.MSI CALL LOCII.K.IB.I.M.MS1 IFIIAI 30.60.IFIIB1 50.60.1445141J-1141 18.441K-RIIM/KRIIR/EAIIMPAIIBI 60 CONTINUE RETURN ENO KM 1 MATA 7 MATA 3 MATA 4 MATA 5 MATA 6 MATA 7 MATA MATO 9 MATA /0 MATO 11 MAYA 12 MATA 13 MATO /4 MATA /5 MATO 16 NATO 17 MATO 19 MATA 19 MATA 70
SADD Purpose: Add a scalar to each element of a matrix to form a resultant matrix. Usage: CALL SADD(A,C,R,N, M, MS) Description of parameters: A - Name of input matrix. C - Scalar. R - Name of output matrix. N - Number of rows in matrix A and R. M - Number of columns in matrix A and R. MS - One digit number for storage mode of matrix A (and R): 0 - General. 1 - Symmetric. 2 - Diagonal. Remarks: None. Subroutines and function subprograms required: LOC Method: Scalar is added to each element of matrix.

CCUT Purpose: Partition a matrix between specified columns to form two resultant matrices. Usage: CALL CCUT (A, L, R, S, N, M, MS) Description of parameters: A - Name of input matrix. L - Column of A to the left of which partitioning takes place. R - Name of matrix to be formed from left portion of A. S - Name of matrix to be formed from right portion of A. N - Number of rows in A. M - Number of columns in A. MS - One digit number for storage mode of matrix A: 0 - General. 1 - Symmetric. 2 - Diagonal. Remarks: Matrix R cannot be in same location as matrix A. Matrix S cannot be in same location as matrix A. Matrix R cannot be in same location as matrix S. Matrix R and matrix S are always general matrices. Subroutines and function subprograms required: LOC Method: Elements of matrix A to the left of column L are moved to form matrix R of N rows and L-1 columns. Elements of matrix A in column L and to the right of L are moved to form matrix S of N rows and M-L+ 1 columns.
SUBROUTINE CCUTIA.L.R.S.N.M.MSI DIMENSION A111011/tS111 IRO 10.0 DO 70 JI.M DO 70 191.N FINO LOCATION IN OUTPUT MATRIX AND SET Tr' ZERO IFIJ LI 20,10.STISI90.0 GO TO ,IR9IR.1 RIIRI90.0 LOCATE ELEMENT FOR ANY MATRIX STORAGE MODE 30 CALL LOCII.J.W.N.M.MSI TEST FOR ZERO ELEMENT IN DIAGONAL MATRIX IFIIJI 40.70,40 DETERMINE WHETHER RIGHT OR LEFT OF L 40 IFIJ L) 60.50.50 SO SIIS19411J/ GO TO RIIR11111JCONTINUE RETURN END CCUT CCUT CCUT CCUT CCITT CCUT CCUT CCUT CCUT CCUT CCUT CCUT CCUT CCUT CCUT CCUT CCUT CCUT CCUT CCUT CCUT CCUT CCUT CCUT CCUT S 10 II

C C C 30

Mathematics Matrices 79
RTIE Purpose: Adjoin two matrices with same column dimension to form one resultant matrix. (See Method.) Usage: CALL RTIE(A,B,R,N,M,MSA,MSB,L) Description of parameters: A - Name of first input matrix. B - Name of second input matrix. R - Name of output matrix. N - Number of rows in A. M - Number of columns in A, B, R. MSA - One digit number for storage mode of matrix A: 0 - General. 1 - Symmetric. 2 - Diagonal. MSB - Same as MSA except for matrix B. - Number of rows in B. L Remarks: Matrix R cannot be in the same location as matrices A or B. Matrix R is always a general matrix. Matrix A must have the same number of columns as matrix B. Subroutines and function subprograms required: LOC Method: Matrix B is attached to the bottom of matrix A. The resultant matrix R contains N+ L rows and M columns.
SUBROUTINE RTIEIAloBtRO,MoMSAvMSB,LI DIMENSION A111.81110111 NN.N
CTIE Purpose: Adjoin two matrices with same row dimension to form one resultant matrix. (See Method.) Usage: CALL CTIE(A, B, R, N, M, MSA, MSB, L) Description of parameters: A- Name of first input matrix. B- Name of second input matrix. R- Name of output matrix. N- Number of rows in A, B, R. M- Number of columns in A. MSA - One digit number for storage mode of matrix A: 0 - General. 1 - Symmetric. 2 - Diagonal. MSB Same as MSA except for matrix B. L Number of columns in B. Remarks: Matrix R cannot be in the same location as matrices A or B. Matrix R is always a general matrix. Matrix A must have the same number of rows as matrix B. Subroutines and function subprograms required: LOC Method: Matrix B is attached to the right of matrix A. The resultant matrix R contains N rows and M+ L columns.

SUBROUTINE CEL2IREStAK.A.B.IERI 1E71.0 TEST MODULUS GE0.14AK.AK IFIGE011.2.IER./ RETURN SET RESULT VALUE OVERFLOW 2 IFI813.5.RES.-1.E38 RETURN
1 CELCELCELCEL2 CELCEL7 CEL2 CELCELCELCELCEL2 12
CEL2 CEL2 CEL2 CEL2 CEL2 CEL2 CEL2 CEL2 CEL2 CEL2 CEL2 CEL2 CEL2 CEL2 CEL2 CEL2 CEL2 CEL2 CEL2 CEL2 CEL2 CEL2 CEL2
Mathematics Special Operations and Functions 107
EXPI This subroutine computes the exponential integral in the range from -4 to infinity. For positive x, the exponential integral is defined as:
The coefficients An are given in the article by Luke/Wimp. * Using only nine terms of the above infinite series results in a truncation error e (x) with:

E (x) = f e 1 x

dt, x> 0

-x <e

-x -8 < e 0.82 10
This function, E l (x), may be analytically continued throughout the complex plane, and defines a multivalued complex function. However, for any given real argument, this extended multivalued function has a unique real part. The subroutine EXPI computes this unique real number for x -4, x 0. For negative x, the real part of the extended exponential integral function is equal to -Ei (-x), where
Transformation of the shifted Chebyshev polynominals to ordinary polynomials finally leads to the approximation:

EXPI(x) = e-x

aV() for x ?. 4
The coefficients of this approximation given to eight signification digits are: a0 = 0.24999 999

Ei )(y = - f e

dt, y > 0
a = -0.l a2 = 0.a = -0.3 a = 0.4 a = -0.5 a = 0.6 a = -0.7 a = 0.2. Approximation in the range x I
(f denotes Cauchy principal value. ) For x = 0, a singularity of the function, the program returns 1. 0 x 10 38. No action is taken in case of an argument less than -4. Polynomial approximations which are close to Chebyshev approximations over their respective ranges are used for calculation. 1. Approximation in the range x 4. A polynomial approximation is obtained by means of truncation of the Expansion of E 1 (x) in terms of shifted Chebyshev Polynomials Tn* E l (x) e -x E AT *(-) , for 4 nn x x n=0
A polynomial approximation is obtained by means of telescoping of the Taylor series of the function: x f (e -t -1) dt = -lnx - C - E (x), 1
*Luke/Wimp, "Jacobi Polynomial expansion of a generalized hypergeometric function over a semiinfinite ray", Math. Comp. , Vol. 17, 1963, Iss. 84, p. 400.
0 where C = O. is Euler's constant. This results in the approximation: 14 v EXPI(x) = -ln Ixl + E bv x v=0
with a truncation error E absolutely less than 3 x 10- 8. The coefficients of this approximation given to eight significant digits are: b b b 0 l 4

= = = = = = =

Subroutines and function subprograms required: None. Method: Definition: RES=integral(EXP(-T)/T, , summed over T from integral(EXP(-T)/T X to infinity). Evaluation: Two different polynomial approximations are used for X greater than 4 and for ABS(X) equal or less than 4. Reference: Luke and Wimp, 'Jacobi Polynomial Expansions of a Generalized Hypergeometric Function over a Semi-Infinite Range', Mathematical Tables and Other Aids to Computation, Vol. 17, 1963, Issue 84, pp. 395-404.

PC LA Purpose: Move polynomial X to Y. Usage: CALL PC LA( Y, IDIMY, X, IDIMX) Description of parameters: Y - Vector of resultant coefficients, ordered from smallest to largest power. IDIMY - Dimension of Y.
- Vector of coefficients for polynomial, ordered from smallest to largest power. IDIMX - Dimension of X. Remarks: None. X Subroutines and function subprograms required: None. Method: IDIMY is replaced by IDIMX and vector X is moved to Y.
SUBROUTINE PCLA IY.IDIMY.11,1/11MX1 DIMENSION XIII.Y111 I0IMY.10IMX IFIIDINA/ 30.30.10 1.0 DO 20 (.1,1oimx 20 Y11/.XRETURN ENO PCLA PCLA PCLA PCLA PCLA PCLA PCLA PCLA I
SUBROUTINE PSUB12.101MZ.X.IDINX.Y.101MY/ DIMENSION 2111.X11/eY111 TEST DIMENSIONS OF SUMMANDS NDIN=IDIMX IF 1101101-1DINYI 10.20.NOIMIDINY 20 IF (NOM 90.90.DO 80 1.1001M IF 11-1DIMX1 40,40.IF 11-10INYI 50.50,111)=8(11Y111 GO TO 1111.Y11/ GO TO 2111.0(CONTINUE 90 101142NDIN RETURN END

POUR POUR

POUR 4 POUR 5 PSUB 6 PSUB 7 POUR 8 POUR 9 PSUB 10 POUR 11 POUR 12 POUR I3 POUR 14 PSUB 15 POUR 16 PSUB 17 POUR 18 PSUB 19
PMPY Purpose: Multiply two polynomials. Usage: CALL PMPY( Z, IDIMZ, X, IDIMX, Y, IDIMY) Description of parameters: - Vector of resultant coefficients, Z ordered from smallest to largest power. IDIMZ - Dimension of Z (calculated). X - Vector of coefficients for first polynomial, ordered from smallest to largest power. IDIMX - Dimension of X (degree is IDIMX-1). Y - Vector of coefficients for second polynomial, ordered from smallest to largest power. IDIMY - Dimension of Y (degree is IDIMY-1). Remarks: Z cannot be in the same location as X. Z cannot be in the same location as Y. Subroutines and function subprograms required: None. Method: Dimension of Z is calculated as IDIMX+IDIMY-1. The coefficients of Z are calculated as sum of products of coefficients of X and Y, whose exponents add up to the corresponding exponent of
Purpose: Subtract one polynomial from another. Usage: CALL PSUB(Z, IDIMZ, X, IDIMX, Y, IDIMY) Description of parameters: Z - Vector of resultant coefficients, ordered from smallest to largest power. IDIMZ - Dimension of Z (calculated). X - Vector of coefficients for first polynomial, ordered from smallest to largest power. IDIMX - Dimension of X (degree is IDIMX-1). Y - Vector of coefficients for second polynomial, ordered from smallest to largest power. IDIMY - Dimension of Y (degree is IDIMY-1). Remarks: Vector Z may be in same location as either vector X or vector Y only if the dimension of that vector is not less than the other input vector. The resultant polynomial may have trailing zero coefficients. Subroutines and function subprograms required: None. Method: Dimension of resultant vector IDIMZ is calculated as the larger of the two input vector dimensions. Coefficients in vector Y are then subtracted from corresponding coefficients in vector X.

PGCD Purpose: Determine greatest common divisor of two polynomials. Usage: CALL PGCD(X, IDIMX, Y, IDIMY, WORK, EPS, IER)
Description of parameters: - Vector of coefficients for first polyX nomial, ordered from smallest to largest power. IDIMX - Dimension of X. - Vector of coefficients for second Y polynomial, ordered from smallest to largest power. This is replaced by greatest common divisor. IDIMY - Dimension of Y. WORK - Working storage array. EPS - Tolerance value below which coefficient is eliminated during normalization. IER - Resultant error code where: No error. IER=O X or Y is zero polynoIER=1 mial. Remarks: IDIMX must be greater than IDIMY. IDIMY=1 on return means X and Y are prime, the GCD is a constant. Subroutines and function subprograms required: PDIV PNORM Method: Greatest common divisor of two polynomials X and Y is determined by means of Euclidean algorithm. Coefficient vectors X and Y are destroyed and greatest common divisor is generated in Y.
SUBROUTINE PGC0111.101MX,Y,IDIMY,WORK.EPS,IERI DIMENSION KI/1.Y11/sM00KIII DIMENSION REQJIRED FOR VECTOR NAMED WORK IS 1 CALL PDIVIWORK.NDIN.A.IDINX,TtIDINT.EPS.IFRI IFIIER1 5.2,IFIIDIMXI 5,593 INTERCHANGE X AND 1, 3 DO 4 J19101MY KORKI/14XIJI XIATIJI 4 11.114WORKIII NOIM4IDIMA
IDIMX - Dimension of X. It is replaced by final dimension. EPS - Tolerance below which coefficient is eliminated. Remarks: If all coefficients are less than EPS, result is a zero polynomial with IDEMX=0 but vector X remains intact. Subroutines and function subprograms required: None. Method: Dimension of vector X is reduced by one for each trailing coefficient with an absolute value less than or equal to EPS.
SUBROUTINE PNORMIX.IDINX.EPSI DIMENSION XIII 4 IEIIDIMXI 4,4,IPIABSIXIIOIMIM-EPSI 3,3,IDIMXIDIMX, GO To I 4 RETURN END
PNORM I PNORM 7 PNORM 3 PNORM 4 PNORM 5 PNORM 6 MORN 7 PNORM

101M/P-10 l

PGCD MCI) PGCD PGCD PGCD PGCD PG00 MD POLO /GC PGCD P GCD PGCD PGCD PGCD

PGCD EGCO

10101Y6NDIM GO TO 1 5' RETURN END
PNORM Purpose: Normalize coefficient vector of a polynomial. Usage: CALL PNORM(X, IDIMX, EPS) Description of parameters: X - Vector of original coefficients, ordered from smallest to largest power. It remains unchanged. Mathematic Polynomial Operation 127
APPENDIX A: ALPHABETIC GUIDE TO SUBROUTINES AND SAMPLE PROGRAMS, WITH STORAGE REQUIREMENTS The following alphabetic index lists the number of characters of storage required by each of the sub- routines in the Scientific Subroutine Package. The figures given were obtained by using 1130 Monitor FORTRAN, Version 2, Modification Level 1. Storage requirements are not given for the sample subroutines. Storage Required (Words) 94 Storage Required (Words) 166

SUB MX

MOMEN TTSTT ORDER AVDAT TRACE CHISQ UTEST
TWOAV QTEST SRANK KRA NK WTEST RANK TIE (111xe RANDU GAUSS
MSTR MFUN RE C P LOC CONVT ARRAY PADD PADDM PC LA PSUB PMPY PDIV PQSD PVAL PVSUB PC LD PILD PDER PINT PGCD PNORM
storage conversion matrix transformation by a function reciprocal function for MFUN location in compressed-stored matrix single precision, double precision conversion vector storage--double dimensioned conversion add two polynomials multiply polynomial by constant and add to another polynomial replace one polynomial by another subtract one polynomial from another multiply two polynomials divide one polynomial by another quadratic synthetic division of a polynomial value of a polynomial substitute variable of polynomial by another polynomial complete linear division evaluate polynomial and its first derivative derivative of a polynomial integral of a polynomial greatest common divisor of two polynomials normalize coefficient vector of polynomial
AVCAL MEANQ DMATX DISCR LOAD VARMX AUTO CROSS SMO EXSMO MINV EIGEN SIM Q QSF QATR RK1 RK2
FORIF FORIT RTWI RTMI RTNI
SUBROUTINES WHOSE ACCURACY IS DATA DEPENDENT The accuracy of the following subroutines cannot be predicted because it is dependent on the characteristics of the input data and on the size of the problem. The programmer using these subroutines must be aware of the limitations dictated by numerical analyses considerations. It cannot be assumed that the results are accurate simply because subroutine execution is completed. Subroutines in this category are: CORRE MULTR GDATA CANOR NROOT means, standard deviations, and correlations multiple regression and correlation data generation canonical correlation eigenvalues and eigenvectors of a special nonsymmetric matrix
and A operation mean square operation means and dispersion matrix discriminant functions factor loading varimax rotation autocovariances crosscovariances application of filter coefficients (weights) triple exponential smoothing matrix inversion eigenvalues and eigenvectors of a real, symmetric matrix solution of simultaneous linear, algebraic equations integral of tabulated function by Simpson's Rule integral of given function by trapezoidal rule integral of first-order differential equation by Runge-Kutta method tabulated integral of first-order differential equation by Runge-Kutta method solution of a system of first-order differential equations by Runge-Kutta method Fourier analysis of a given function Fourier analysis of a tabulated function refine estimate of root by Wegstein's iteration determine root within a range by Mueller's iteration refine estimate of root by Newton's iteration real and complex roots of polynomial

7.781777 <C77777.05

Figure 17. Output listing
the user with the program modification, the following general rules are supplied in terms of the sample problem: 1. Changes in the dimension statements of the main program, MCANO: a. The dimension of arrays XBAR, STD, CANR, CHISQ, and NDF must be greater than or equal to the total number of variables m (in = p + q, where p is the number of left-hand variables and q is the number of right-hand variables). Since there are seven variables, four on left and three on right, the value of m is 7. b. The dimension of array RX must be greater than or equal to the product of m x m. For the sample problem this product is 49 = 7 x 7. c. The dimension of array R must be greater than or equal to (rn + 1)m/2. For the sample problem this number is 28 = (7 + 1)7/2. d. The dimension of array COEFL must be greater than or equal to the product of p x q. For the sample problem this product is 12 = 4 x 3. e. The dimension of array COE FR must be greater than or equal to the product of q x q. For the sample problem this product is 9 = 3 x 3.
Sample Main Program for Canonical Correlation MCANO Purpose: (1) Read the problem parameter card for a canonical correlation, (2) Call two subroutines to calculate simple correlations, canonical correlations, chi-squares, degrees of freedom for chi-squares, and coefficients for left and right hand variables, namely canonical variates, and (3) Print the results. Remarks: I/O specifications transmitted to subroutines by COMMON. Input card: Column 2 MX - Logical unit number for output. Column 4 MY - Logical unit number for input. The number of left-hand variables must be greater than or equal to the number of righthand variables. Subroutines and function subprograms required: CORRE (which, in turn, calls the input subroutine named DATA. ) CANOR (which, in turn, calls the subroutines MINV and NROOT. NROOT, in turn, calls the subroutine EIGEN. ) Appendix D Sample Programs 159
Method: Refer to W. W. Cooley and P. R. Lohnes, Multivariate Procedures for the Behavioral Sciences', John Wiley and Sons, 1962, chapter 3.
// FOR *10CSICARD.TYPEWRITER.1132 PRINTER) * ONE WORD INTEGERS KANO 1 SAMPLE MAIN PROGRAM FUR CANONICAL CORRELATION - MCANU C THE FOLLOWING DIMENSIONS MUST BE GREATER THAN OR EQUAL TO THE MCANO. 2 TOTAL NUMBER OF VARIABLES M IM=NP0Q. WHERE MP IS THE NUMBER MCANJ 3 C MCANU 4 C OF LEFT HAND VARIABLES. ANU NO IS THE NUMBER OF RIGHT HAND MCANO 5 C VARIABLES). MCANO 6 LIMENSION XBAR19)01.0(9).CANA19).CHISQ19)00F(9) MCANO 7 C THE FOLLOWING DIMENSION MUSE BE GREATER THAN OR EQUAL TO THE PRODUCT OF M*M. MCANU 8 C MCANU 9 DIMENSION R%181) THE FULLOWING DIMENSION MUST BE GREATER THAN OR EQUAL TO MCANO 10 C NCANO LI IM*110/2. MCANO 12 DIMENSION 8145) MCANO 13 C THE FOLLOWING DIMENSION MUST dE GREATER THAN OR EQUAL TO THE PRODUCT OF MP*MQ MCANO 14 C MCANU L5 DIMENSION CCEFL181) THE FOLLOWING DIMENSION MUST BE GREATER THAN OR EQUAL TO THE MCANO 16 C SCANS 17 PRODUCT OF MQ*MJ. MCANO 18 DIMENSIEN COEFRI251 COMMON M200 MCANO FORMAT(A402.15.212) MCANU FORMATI////27H CANONICAL CORRELATION A402//22H NO. OF OBSEMCAND 21 INVATIONS0X.14/29H NO. OF LEFT HANG VARIABLES.15/30H NG. OF RIMCANU 22 2GHT HAND VARIABLESpI4/) MCANG FORMAT(//6H MEANS/(8F15.5)) MCANO FORMATI//20H STANDARD DEVIATIONS/18F15.5)) MCANU FORMAT)/125H CORRELATION COEFFICIENTS) MCANO FORMAT(//4H ROW.13/110F12.51) MCANO FORMAT)////12H NUMBER OF.74.7HLARGESTaX,13HCORRESPOND0G.31Ap7HMCANO 2B 10EGKEES/13H EIGENVALOES.50.10HEIGENVALUE.7Xp9HCANUNICAL.7X0HLAMBNCANO 24 20A.52.10HCHI-SQUARE.7%.2H0F/4X.7HREMOVED.70.9HREMAINING.UpLIHCORRMCANU 30 MCANO 51 3ELATION.32X.7HFREEDCM/) B FORMAT(/17019.5016.5.2F14.5.5X,15) NCANO FORMATI///22H CANONICAL CORRELAJION012.5) MCANU FORMATI//39H COEFFICIENTS FOR LEFT HANG VARIABLES/18F15.5)) MCANU 34 MCANG FORMAT(//40H COEFFICIENTS FOR RIGHT HAND VARIABLES/I8F15.51) 12 FORMATI2121 MCANJ 36 REA012,120X.MY MCANG 37 READ PROBLEM PARAMETER CARD MCANG 38 C MCANO REAP/ 07.1)PROR100P.MQ PR MEANS 40 PROBLEM MINDER (NAY BE ALPHAMERIC) MEANS 41 PRI PROBLEM NUMBER (CONTINUED) C NUMBER OF OBSERVATIONS SCANS 42 N MP NUMBER OF LEFT NANO VARIABLES MCANO 43 C MD NUMBER OF RIGHT HAND VARIABLES C SLANG 44 WRITE INX,20R,PRI.N.MP.MQ MCANJ 45 M=MP.MQ MCANU 46 MCANO 47 10=0 X=0.0 MCANO 48 MCANO 49 CALL CORRE IN0p1G00BAR.STO.RX0pCANR,CHISQ.COEFL) PRINT MEANS. STANDARD DEVIATIONS. AND CORRELATION MCANU 50 C COEFFICIENTS OF ALL VARIABLES MCANO 51 WRITE IMX.3)IXBARII).1=1-0) MCANO 52 WRITE IMR.411STD11),1=10) MCANU 53 WRITE (NX.5) MCANU 54 DO = 10 MCANO 55 DO 150 J=10 MCANO 56 IFII-3) 120. 130. 130 MCANO L=IfIJ4J-J)/2 MCANO 58 GO TO 140 SLANG 1.6.1*(1*1-I1/2 MCANO CANRIJI=RIL) MCANO CONTINUE MCANO WRITE (MA.6)1.(CANROI,J 6 10) MCANO 63 CALL CANOR (61,MP.MO,R,XBAR.STO,CANR,CHISO.NUF,GOEFR,GOEFL,RX1 MCANO 64 C PRINT EIGENVALUES. CANONICAL CORRELATIONS. LAMBDA. CHI-SQUARES MCANO 65 DEGREES OF FREEDOMS C MCANO 66 WRITE (MX,7) MCANG 67 DO 170 1=10Q MCANO 68 N161-1 MCANO 69TEST WHETHER EIGENVALUE IS GREATER THAN ZERO C MCANO 70 IFIXBARII)) 165, 165, 17G MCANO MM=N1 MCANO 72 GO TO 175 MCANO WRITE IMXIBINI.X8AR11).CANR(1)01. 011) 1 CHISQ(1).NDF11) MCANO 74 MM=MQ MCANO 75 PRINT CANONICAL COEFFICIENTS MCANO 76 C 175 N1=0 MCANO 77 N2=0 MCANO 78 DO 200 1=10M MCANO 79 WRITE (MX0)CANR(1) NCANO BO DO 180.)=101P MCANO 81 N1=5161 MCANO XBARI.1)=COEFL(NI) MCANO 83 WRITE IMX.101(XBARIJI.J=I I MP) NCANO 84 DO 190 J=100 MCANO 85 N20241 MCANO XBAR(J)=COEFR(N21 MCANO 87 WRITE IMX.11)1XBARIJ),J 6 100/ MLANO CONTINUE MCANO 89 GO TO 100 MCANO 90 END MCANJ 91 // DUP WS UA MCANO *STORE // %ER NCANO 01 *LOCALMCANO.CORRE.CANOR 12 SAMPLE/SO 151 188

Program

Description The analysis-of-variance sample program consists of a main routine, ANOVA, and three subroutines: AVDAT AVCAL MEANQ Capacity The capacity of the sample program and the format required for data input have been set up as follows: 1. Up to six-factor factorial experiment 2. Up to a total of 1600 data points. The total number of core locations for data points in a problem is calculated as follows: k T = H (LEVEL. + 1) i = 1 I where LEVELi = number of levels of i th factor k number of factors H = notation for repeated products 3. (12 F6.0) format for input data cards Therefore, if a problem satisfies the above conditions it is not necessary to modify the sample program. However, if there are more than six are from the Scientific Subroutine Package t17 - 22 - 27 - 32 - 67 - 70 7-- 15 116
If there are more than eleven factors, continue to the second card in the same manner.
Appendix D Sample Programs 161

Columns 1 2-5

Contents Label for the twelfth factor Number of levels for the twelfth factor
etc. Leading zeros are not required to be keypunched. Data Cards Data are keypunched in the following order: X1111, X2111, X3111, X4111, X1211, X2211, X3211, ,
X4332. In other words, the innermost subscript is changed first; namely, the first factor, and then second, third, and fourth subscripts. In the sample problem, the first subscript corresponds to factor A and the second, third, and fourth subscripts to factors B, C, and R. Since the number of data fields per cards is twelve, implied by the format (12F6. 0), each row in Table 5 is keypunched on a separate card. Deck Setup Deck setup is shown in Figure 18. Sample The listing of input cards for the sample problem is presented at the end of the sample main program.
Figure 18. Deck setup (analysis of variance)
Output Description The output of the sample analysis-of-variance program includes the numbers of levels of factors as input, the mean of all data, and the table of analysis of variance. In order to complete the analysis of variance properly, however, certain components in the table may need to be pooled. This is accomplished by means of summary instructions that specifically apply to the particular experiment as presented in Table 6.

// DUP

00020008001100 0.7507505 0.6194650 0.7086843 0.7241215 1.0000000 0.7986682 0.7446845 0.6144597 0.7029582 0.6244183 0.7986682 1.0000000 0.0.8152021 0.5838704 0.7468050 0.6439786 0.6963269 0.5745585 0.6573101 1.0000000 0.6644085 0.7601008 0.6751766 0.5571068 0.6373449 0.7284786 0.6010878 0.6876602 0.7446845 0.6144597 0.7029582 0.4299425 0.3547574.0.4058519 0.8635910 0.7125728 0.8152021 0.7601008 0.6271802 1.0000000 0.7601008 0.6271802 1.0000000 0.6963269 0.5745585 0.6573101
SUBROUTINE MATIN PURPOSE READS CONTROL CARO AND MATRIX DATA ELEMENTS. FROM LOGICAL UNIT S USAGE CALL MATINIICODE.A.ISIZEgIRCOvICOLIDIS.IER/ DESCRIPTION OF PARAMETERS ICODE-UPON RETURN. ICODE WILL CONTAIN FOUR DIGIT IDENTIFICATION CODE FROM MATRIX PARAMETER CARD -DATA AREA FOR INPUT MATRIX A ISIZE-NUMBER OF ELEMENTS DIMENSIONED BY USER FOR AREA A IRON -UPON RETURN. IRON WILL CONTAIN ROW DIMENSION FROM MATRIX PARAMETER CARD (COL -UPON RETURN, ICOL WILL CONTAIN COLUMN DIMENSION FROM MATRIX PARAMETER CARD IS -UPON RETURN. IS WILL CONTAIN STORAGE MODE CODE FROM MATRIX PARAMETER CARD WHERE IS.0 GENERAL'MATRIX IS./ SYMMETRIC MATRIX IS. 2 DIAGONAL MATRIX IER -UPON RETURN. IER WILL CONTAIN AN ERROR CODE WHERE IER.0 NO ERROR (SIZE IS LESS THAN NUMBER OF ELEMENTS IN IER.1 INPUT MATRIX 189-2 INCORRECT NUMBER OF DATA CARDS REMARKS NONE SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED LOC METHOD SUBROUTINE ASSUMES THAT INPUT MATRIX CONSISTS OF PARAMETER CARD FOLLOWED BY DATA CARDS PARAMETER CARD HAS THE FOLLOWING FORMAT COL. 1- 2 BLANK COL. 3- 6 UP TO FOUR DIGIT IDENTIFICATION CODE COL. 7-10 NUMBER OF ROWS IN MATRIX COL.11-04 NUMBER OF COLUMNS IN MATRIX COL.15-16 STORAGE MODE OF MATRIX WHERE 0 - GENERAL MATRIX 1- SYMMETRIC MATRIX 2 - DIAGONAL MATRIX DATA CARDS ARE ASSUMED TO HAVE SEVEN FIELDS OF TEN COLUMNS EACH. DECIMAL POINT MAY APPEAR ANYWHERE IN A FIELD. IF NO DECIMAL POINT IS INCLUDED. IT IS ASSUMED THAT THE DECIMAL POINT IS AT THE END OF THE 10 COLUMN FIELD. NUMBER IN EACH FIELD NAY BE PRECEDED BY BLANKS. DATA ELEMENTS MUST BE PUNCHED BY ROW. ROW MAY CONTINUE FROM CARD TO CARD. HOWEVER EACH NEW ROW MUST START IN THE FIRST FIELD OF THE NEXT CARD. ONLY THE UPPER TRIANGULAR PORTION OF A SYMMETRIC OR THE DIAGONAL ELEMENTS OF A DIAGONAL MATRIX ARE CONTAINED ON DATA CARDS. THE FIRST ELEMENT OF EACH NEW ROW WILL BE THE DIAGONAL ELEMENT FOR A MATRIX WITH SYMMETRIC OR DIAGONAL STORAGE MODE. COLUMNS 71-80 OF DATA CARDS MAY BE USED FOR IDENTIFICATION. SEQUENCE NUMBERING. ETC. THE LAST DATA CARD FOR ANY MATRIX MUST BE FOLLOWED BY A CARD WITH A 9 PUNCH IN COLUMN 1. A.1517E.IROW.ICOL.ISOER1 SUMROUTINE MATINIICODE, DIMENSION All) UIMENSION CARDI81 COMMON MX,MY I FORMATI7F10.01 FORMAT116.FORMAT1111 IDC=7 IER.0 READ) RY.211C00EOROM.ICOLoIS CALL LOCIIROW.101.1CNT.IROW.ICOL.ISI IFIISI2E-ICNT16.7,IER=IF 11CNT138.38.1COLT.ICOL IROCR=I COMPUTE NUMBER OF CARDS FOR THIS ROM 11 IRCDS.11COLT-11/1DC1 IFIES-1115.15.IRCOS.1 SET UP LOOP Elk NUMBER OF CARDS IN ROM 15 DO 31 K=1.IRCOS REAOIMY.I1ICARD111.11,10C1 SKIP THROUGH DATA CARDS IF INPUT AREA TOO SMALL IFIIER/16.16,L.0 COMPUTE COLUMN NUMBER FOR FIRST FIELD IN CURRENT CARD JS=1K-IIIDCICOL-ICOLT+1 MATIN 1 MATIN 2 MATIN 3 MATIN 4 MATIN 5 MATIN 6 MATIN 7 MATIN 9 MATIN 9 MATIN RATIN II MATIN 12 MATIN 13 MATIN /4 MAHN 15 MATIN 16 MATIN 17 RATIN 19 MATIN 19 MATIN 2^ MATIN 71 MATIN 22 MATIN 23 MATIN 24 MATIN 75 RATIN 76 MATIN 27 MATIN 7F

I/O Specification Card Each integration requires a parameter card with the following format: For Sample Problem 1.0

Columns 1 - 5

Contents Up to 5-digit numeric identification code
6 - 10 Number of tabulated values for this function 11 - 20 Interval between tabu- lated values

y out Anut 48

mxnUT 49 ABOUT 50

MxnuT 46

NUMERICAL QUADRATURE INTEGRATION Problem Description The tabulated values of a function for a given spacing are integrated. Multiple sets of tabulated values may be processed. Program Description The numerical quadrature integration program consists of a main routine QDINT, and one subroutine, QSF, from the Scientific Subroutine Package.
The first two parameters consist of up to five digits with no decimal point (FORMAT (215)). Note that the second parameter may not exceed 500. The third parameter consists of up to ten digits (FORMAT) (F10.0). Data Cards Data cards are assumed to be seven fields of ten columns each. The decimal point may appear anywhere in the field, or be omitted, but the number must be right-justified. The number in each field may be preceded by blanks. Columns 71 through 80 of the data cards may be used for identification, sequence numbering, etc. If there are more than seven tabulated values, the values should continue from card to card with seven values per card, until the number of values specified in the parameter card has been reached.
A blank card following the last set of data terminates the run. Deck Setup The deck setup is shown in Figure 28. Sample A listing of input cards for the sample problem is presented at the end of the sample main program. Output Description The identification code number, number of tabulated input values, the interval for the tabulated values, and the resultant integral values at each step are printed. Sample The output listing for the sample problem is shown in Figure 29. Program Modification Noting that storage problems may result, as previously discussed in "Sample Program Description" , the maximum number of tabulated values acceptable
to the sample program may be increased. Input data in a different format can also be handled by providing a specific format statement. 1. Modify the DIMENSION statement in QDINT so that the size of array Z is equal to the maximum number of tabulated values. 2. Changes to the format of the parameter cards and data cards may be made by modifying FORMAT statements 10 and 32, respectively, in QDINT.