Google
 

Trailing-Edge - PDP-10 Archives - decuslib10-12 - 43,50550/formal.doc
There are no other files named formal.doc in the archive.



                      GENERAL MATRIX ALGEBRA PROGRAMS
      SUBROUTINE ABS(X,Y,NR,NC)
      SUBROUTINE XABS(NRMX,NRMY,X,Y,NR,NC)
  TITLE: ABSolute value of a matrix
  PURPOSE: To make absolute value of matrix.
  SYMBOLS:
    NRMX,NRMY = number of rows of X and Y, respectively, in calling
      program
    X = input matrix
    Y = output matrix
    NR,NC = # of rows and columns, respectively, of X and Y
  PROCEDURE: Elements of X are made into absolute value
    and the result is stored in the corresponding elements of Y.
    X or Y may be the same matrix in the calling program.
      SUBROUTINE  ADD(X,Y,Z,NR,NC)
      SUBROUTINE XADD(NRMX,NRMY,NRMZ,X,Y,Z,NR,NC)
 TITLE: ADDition of 2 matrices
 PURPOSE: To make addition of 2 matrices.
 SYMBOLS:
   NRMX,NRMY,NRMZ = number of rows of X,Y,and Z, respectively,
     in calling program
   X = 1st input matrix
   Y = 2nd input matrix
   Z = resultant matrix
   NR,NC = # of rows and columns, respectively, of X, Y, and Z
 PROCEDURE: 
   Corresponding elements of X and Y are added and the 
   result is stored in the corresponding elements of Z.  Any of
   X, Y, or Z may be the same matrix in the calling program.
      SUBROUTINE  CHLSK(X,Y,N)
      SUBROUTINE XCHLSK(NRMX,NRMY,X,Y,N)
 TITLE: CHoLeSKy decomposition
 PURPOSE: To obtain a lower triangular matrix Y, such that YY' = X,
   a symmetric matrix.
 SYMBOLS:
   NRMX,NRMY = dimension of X and Y, respectively, in calling program.
     Must be .GE. 10.
   X = input symmetric matrix
   Y = resulting lower triangular matrix
   N = order of X and Y
 PROCEDURE:
   X is converted to a lower triangular matrix by transformation of
   successive rows of its lower triangle.  The result is stored in Y.
   The input matrix must be symmetric and positive definite.
   X and Y may be the same matrix in the calling program.
NOTE:
   Algorithm reference: G.W.Stewart, Introduction to Matrix Computation
   Academic Press, 1973, pp.139-142.
      SUBROUTINE  CHSYM(A,N)
      SUBROUTINE XCHSYM(NRMA,A,N)
 TITLE: CHecks SYmmetry of Matrix
 PURPOSE: To check symmetry of matrix A
 SYMBOLS:
   NRMA   = number of rows of A in calling program
   A	  = matrix to be checked; input/output
   N	  = order of matrix; input/output
 PROCEDURE: The symmetry of matrix X is checked in the subroutine.
   The elements of lower half of trianglar matrix of X is checked
   with the elements of upper half of X.
      SUBROUTINE  COLSM(X,Y,NR,NC)
      SUBROUTINE XCOLSM(NRMX,NRMY,X,Y,NR,NC)
 TITLE: COLumn SuMmation
 PURPOSE: To add the elements down columns of the input matrix and
   store as a row vector.
 SYMBOLS:
   NRMX,NRMY = row dimensions of X and y, respectively, in calling program
   X = input matrix whose columns are to be summed
   Y = output matrix containing sums in first row
   NR = # of rows of X
   NC = # of cols of X & Y
 PROCEDURE: The columns of X are summed simultaneously and stored in Y.
   X & Y may be the same matrix in the calling program.
      SUBROUTINE  COPY(X,Y,NR,NC)
      SUBROUTINE XCOPY(NRMX,NRMY,X,Y,NR,NC)
 TITLE: COPY
 PURPOSE: To copy one matrix into another.
 SYMBOLS:
   NRMX,NRMY = row dimensions of X and Y, respectively, in calling 
     program
   X = matrix to be copied
   Y = matrix into which X is copied
   NR = # of rows of X and Y
   NC = # of cols of X and Y
      SUBROUTINE  DET(X,D,N,Y)
      SUBROUTINE XDET(NRMX,NCMX,NRMY,NCMY,X,D,N,Y)
 TITLE: DETerminant
 PURPOSE: To find the determinant of a square, real matrix.
 SYMBOLS:
   NRMX,NRMY = row dimsions of X and Y, respectively, in 
     the calling program
   NCMX,NCMY = column dimsions of X and Y, respectively, in the 
     calling program
   X = input square matrix
   D = determinant of X
   N = order of X
   Y = scratch matrix
   XMINVZ = matrix inversion routine which returns determinant
 PROCEDURE:
   The determinant is found by QMINVZ during the process of
   inversion.	Thus, if D is not 0, Y contains the inverse 
   of X on  output.  X and Y may be the same matrix in the 
   calling program.  Notice that X will be desreoyed if this 
   is so.  N must be strictly less than NRMX.
      SUBROUTINE  DIAG(X,Y,N)
      SUBROUTINE XDIAG(NRMX,NRMY,X,Y,N)
 TITLE: DIAGonal Matrix
 PURPOSE: To take the diagonal of a matrix and create a diagonal matrix.
 SYMBOLS:
   NRMX,NRMY = row dimensions of X and Y, respectively, in the calling 
     program
   X = input matrix
   Y = output diagonal matrix = Diag(X)
   N = # of rows and cols of X and Y
 PROCEDURE:
   For a given row and column of X, the corresponding off-diagnal
   are set to 0 in Y, and the corresponding diagnal element of Y is set 
   equal to the element of X.  X and Y can be the same matrix
   in the calling program.  Notice that X and Y are assumed to 
   be square and located in the upper left of the corresponding
   matrices in the calling program.
      SUBROUTINE  EIGEN(R,V,D,N,NVE)
      SUBROUTINE XEIGEN(NRMR,NRMV,NRMD,R,V,D,N,NVE)
 TITLE: EIGEN Solution
 PURPOSE: To find the eigen solution for a symmetric matrix,
      such that R = VDV'.
 SYMBOLS:
   NRMR,NRMV,NRMD = row dimensions of R, V, and D, respectively, in 
     calling program
   R = input symmetric matrix
   V = eigenvectors
   D = diagonal matrix of eigenvalues
   N = order of R
 PROCEDURE: Subroutine  QLRSM is invoked to do the real number crunching.
   After successful execution of that routine, all of the eigenvalues
   are stored in the diagonal elements of D and the off-diagonals are
   zeroed out.   If NVE is .GT. N, it is changed to N; if NVE is .LT. 1
   it is changed to 1.  V contains only NVE columns on out, one for each
   eigenvector.  The order of the input matrices (i.e., NRMR,etc) must 
   be .GE.N.  R, V, and D must be unique matrices in the calling program.
      SUBROUTINE ELEDI(X,Y,Z,NR,NC)
      SUBROUTINE XELEDI(NRMX,NRMY,NRMZ,X,Y,Z,NR,NC,ICHECK)
 TITLE: ELEmentwise DIvision
 PURPOSE: To find the elementwise division of 2 matrices.
 SYMBOLS:
    NRMX,NRMY,NRMZ = number of rows of X, Y, and Z, respectively, in
      calling program
    X = 1st input matrix
    Y = 2nd input matrix
    Z = result matrix
    NR = # of rows of X,Y, and Z
    NC = # of cols of X,Y, and Z
    ICHECK = flag to determine the effect of a zero entry in Y.
           = 0 (on input) stop program when zero is found. (default)
           = 1 (on input) continue but set Z(I,J) to 0.0D0.
  PROCEDURE: Divide an element of X by the corresponding element of Y,
    store in the corresponding element of Z. X, Y, and Z, or any pair
    combination, may be the same matrix in the calling program.
      SUBROUTINE  ELEML(X,Y,Z,NR,NC)
      SUBROUTINE XELEML(NRMX,NRMY,NRMZ,X,Y,Z,NR,NC)
 TITLE: ELEmentwise MuLtiplication
 PURPOSE: To find the elementwise product of 2 matrices.
 SYMBOLS:
   NRMX,NRMY,NRMZ = number of rows of X, Y, and Z, respectively, in
     calling program
   X = 1st input matrix
   Y = 2nd input matrix
   Z = result matrix
   NR = # of rows of X,Y, and Z
   NC = # of cols of X,Y, and Z
 PROCEDURE: Multiply an element of X by the corresponding element of Y
   store in the corresponding element of Z.  X, Y, and Z, or any 
   pair combination, may be the same matrix in the calling program.
      SUBROUTINE  GENR8(C,X,NR,NC)
      SUBROUTINE XGENR8(NRMX,C,X,NR,NC)
 TITLE: GENeRATE
 PURPOSE: To generate a matrix containing constants
 SYMBOLS:
   NRMX = number of rows of X in calling program
   C = the constant (real)
   X = the matrix to be generated containing all C values
   NR,NC = # of rows and cols of X
 PROCEDURE: Every element of X is set equal to C.  X may be any 
   sized matrix.  C may be a REAL*4 or double precision number.
      SUBROUTINE  GINVR(X,Y,NR,NC,A,D,C)
      SUBROUTINE XGINVR(NRMX,NRMY,NRMA,NRMD,NRMC,X,Y,NR,NC,A,D,C)
 TITLE: Generalized INVeRse
 PURPOSE: To compute the generalized inverse  (Y) of a matrix X,
   such that XYX = X.
 SYMBOLS:
   NRMX,NRMY,NRMA,NRMD,NRMC = number of rows of X, Y, A, D, and C,
     respectively, in calling program.  Must be .GE. 10.
   X = input matrix
   Y = generalized inverse of X
   NR,NC = # of rows & cols of X = # of cols & rows of Y
   A,D,C = temporary work matrices
   QLRSM = implicit QL eigen routine
 PROCEDURE: The matrix X'X is found and its eigenvalues (D) and
   eigenvectors (A) are computed.  Y is then simply A(D**-1)A'X'.
   Y may be the same as either A or D in the calling program,
   but no other combination (including A and D) may be the same.
      SUBROUTINE  HORIZ(X,Y,Z,NR,NCX,NCY)
      SUBROUTINE XHORIZ(NRMX,NRMY,NRMZ,X,Y,Z,NR,NCX,NCY)
 TITLE: HORIZontal concatenate
 PURPOSE: To concatenate the rows of 2 matrices to create a super matrix
 SYMBOLS:
   NRMX, NRMY, NRMZ = row dimensions of X, Y, and Z, respectively,  
     in calling program
   X = 1st input matrix
   Y = 2nd input matrix
   Z = resulting super matrix = (X:Y)
   NR = # of rows in both X and Y
   NCX = # of cols in X
   NCY = # of cols in Y
 PROCEDURE:
   A matrix Z = (X:Y) is created by storing X in cols 1 thru NCX,
   rows 1 thru NR of Z and then by storing Y in cols (NCX+1) thru NC 
   and row 1 thru NR of Z.  Thus Z will be NR X NC.  All three
   matrices, X, Y, and Z may be the same matrix in the calling 
   program; matrices X and Y may be the same but different from
   Z; matrices X and Z may be the same but different from Y; 
   however, matrices Y and Z may not be the same if X is different
   in the calling program.
      SUBROUTINE  IDENT(X,N)
      SUBROUTINE XIDENT(NRMX,X,N)
 TITLE: IDENTenty matrix generator
 PURPOSE: To create identity matrix of order N.
 SYMBOLS:
   NRMX = row dimension of X in calling program
   N = order of identity matrix
   X = identity matrix
 PROCEDURE: Creates a square identity matrix of order N
   N must be greater than 0.
      SUBROUTINE  INPUT(X)
      SUBROUTINE XINPUT(NRMX,X)
 TITLE: data INPUT from cards
 PURPOSE: To input data from cards.
 SYMBOLS:
   NRMX = row dimension of X in calling program
   X = matrix in which to store input data
   FMT = array for containing input format (read from the first card of
     data dack)
   C = labeled common areas containing NMAT, and DSRN
   NMAT = matrix counter
   DSRN = data set reference number
   QINCR = program number counter subroutine
   NR = number of rows input to X
   NC = number of cols input to X
 PROCEDURE: This subroutine  reads the header record from the data file,
   QINCR is the matrix counter, prints out NMAT,DSRN,NR,NC,& FMT for
   user, and, if all is OK, reads the actual data file.  If, for any 
   reason the end of the file is encountered before expected, process-
   ing terminate after printing an appropriate message.
      SUBROUTINE  INVRS(X,Y,N)
      SUBROUTINE XINVRS(NRMX,NCMX,NRMY,NCMY,X,Y,N)
 TITLE: INVeRSe
 PURPOSE: To invert a square, real matrix.
 SYMBOLS:
   NRMX,NRMY = row dimensions of X and Y, respectively, in calling 
     program
   NCMX,NCMY = column dimensions of X and Y, respectively, in
     calling program
   X = input square matrix
   Y = resulting inverse of X
   N = order of X
   QMINVZ = matrix inversion number cruncher
 PROCEDURE: The order of X is checked to make sure it is correct and
   QMINVZ is called.  X and Y may be the same matrix in the calling
   program.  N must be 2 less than NCMX & NCMY and less than or equal
   to NRMX & NRMY.  If X is singular, a message is printed and 
   execution ceases.
      SUBROUTINE  KRONK(X,Y,Z,NRX,NCX,NRY,NCY)
      SUBROUTINE XKRONK(NRMX,NRMY,NRMZ,X,Y,Z,NRX,NCX,NRY,NCY)
 TITLE: KRONecKer product
 PURPOSE: To find the Kronecker product of 2 matrices.
 SYMBOLS:
   NRMX, NRMY, NRMZ = row dimensions of X,Y, and Z, respectively,  
     in calling program
   X = 1st input matrix
   Y = 2nd input matrix
   Z = result matrix
   NRX,NCX = # of rows & cols of X
   NRY,NCY = # of rows & cols of Y
   NRZ,NCZ = # of rows & cols of Z
 PROCEDURE:
   A single element of X is scalar multiplied by Y and the resulting
   matrix is placed in a partition of Z corresponding to the relative 
   position of the single element of X.  This product is called a
   Kronecker product.  X and Y may be the same address, but Z must
   be different.
      SUBROUTINE  MATML(X,Y,Z,NRX,NXY,NCY)
      SUBROUTINE XMATML(NRMX,NRMY,NRMZ,X,Y,Z,NRX,NXY,NCY)
 TITLE: MATrix MuLtiply
 PURPOSE: To multiply 2 input matrices and return the result in a third
   matrix.
 SYMBOLS:
   NRMX, NRMY, NRMZ = row dimension of X, Y, and Z, respectively, 
     in calling program
   X,Y = the 2 input matrices
   Z = the result matrix = X*Y
   NRX = # of rows of X = # of rows of Z
   NXY = # of columns of X = # of rows of Y
   NCY = # of columns of Y = # of columns of Z
 PROCEDURE: Each row of X is vector multiplied by each column of Y to
   produce the successive elements of Z.  Note that X and Y may be the
   matrix in the calling program. However, Z must be a different matrix
   Z will have NRX rows and NCY columns.
      SUBROUTINE  PAGE
 TITLE: eject PAGE
 PURPOSE: To do a page eject, when called.
      SUBROUTINE  PARTT(X,Y,IBR,IER,IBC,IEC)
      SUBROUTINE XPARTT(NRMX,NRMY,X,Y,IBR,IER,IBC,IEC)
 TITLE: PARTiTion
 PURPOSE: To select off a partition of a matrix.
 SYMBOLS:
   NRMX, NRMY = row dimensions of X and Y, respectively,  
     in calling program
   X = input matrix to be partitioned.
   Y = partition selected out of X.
   IBR,IER = beginning & ending rows, respectively, of partition of X.
   IBC,IEC = beginning & ending cols, respectively, of partition of X.
 PROCEDURE: A partition of X is designated by IBR,IER,IBC,IEC.  The
   submatrix designated by this partition is copied into the upper left
   of Y.  A test is made to ascertain that IER(IEC).GE.IBR(IBC) and that
   IER(IEC) is not greater than NRMX.  If this test is not satisfied, a
   diagnostic error message is printed and execution is terminated.  X
   may be the same matrix in the calling program.
      SUBROUTINE  PRDDI(X,N,PD)
      SUBROUTINE XPRDDI(NRMX,X,N,PD)
 TITLE: PRoDuct of the DIagonal
 PURPOSE: To find the product of the elements in the major diagonal of
   input matrix.
 SYMBOLS:
   NRMX = row dimension of X in calling program
   X = input matrix
   N = # of elements in major diagonal
   PD = the resulting product
 PROCEDURE:
   The elements of the major diagonal of X are successively multiplied
   by the product of the preceding elements.  X may be non-square in 
   which case N, the length of the major diagonal, is equal to MIN(# of 
   rows of X, the number of cols. of X).
      SUBROUTINE  PRINT(X,NR,NC,IFLAG)
      SUBROUTINE XPRINT(NRMX,X,NR,NC,IFLAG)
 TITLE: matrix PRINT
 PURPOSE: To print a matrix in F or D format.
 SYMBOLS:
   NRMX = row dimension of X in calling program
   X = matrix to be printed
   NR = # of rows of X
   NC = # of cols of X
   IFLAG = flag specifying format type (either 'F' or 'D' on input to 
     print)
 PROCEDURE:
   Row and column headings are printed; a double space occurs 
   before returning.
      SUBROUTINE  PUNCH(X,NR,NC)
      SUBROUTINE XPUNCH(NRMX,X,NR,NC)
 TITLE: matrix PUNCH
 PURPOSE: To punch a matrix onto cards.
 SYMBOLS:
   NRMX = row dimension of X in calling program
   X = matrix to be punched
   NR = # of rows of X to be punched
   NC = # of cols of X to be punched
 PROCEDURE:
   A header card of the form recognized by the input subroutine
   formal is punched first.  The number of cards per row is calculated
   in the doubly nested do loop wherein each row of X is punched onto
   successive cards with sequence numbers.
*****************************************************************
*****************************************************************
**							       **
**   As Lederle does not have a card punch this program is     **
**   deactivated!!!					       **
**							       **
*****************************************************************
*****************************************************************
      SUBROUTINE  RANK(X,Y,NR,NC,RNK)
      SUBROUTINE XRANK(NRMX,NCMX,NRMY,X,Y,NR,NC,RNK)
 TITLE:RANK
 PURPOSE: To find the rank of a real matrix.
 SYMBOLS:
   NRMX,NRMY = row dimensions of X and Y, respectively, in calling 
     program
   NCMX = column dimension of X in calling program
   X = input matrix
   Y = scratch array
   NR = # of rows of X
   NC = # of cols of X
   RNK = smallest non-zero eigenvalue to be counted: input by user = 
     rank of X on output.
   XXTX = cross products subroutine
   QTRED1 = EISPACK tridiagonalization subroutine
   QIMTQL = EISPACK eigenvalues of a tridiagonal matrix subroutine
 PROCEDURE: The program calls XTX to find the cross-products of X. 
   The rank of X'X = rank of non-zero eigenvalues of X'X.
   Hence the EISPACK subroutines are issued for finding the
   eigenvalues of X'X.  X and Y may be the same for finding
   the eigenvalues of X'X.  X and Y may be the same 
   for finding the eigenvalues of X'X.  X and Y may be the same 
   address if the user doesn't mind that X will be destroyed.	
   The functional equivalent of 0 is used by the user for RNK on input.
   A number in the range 1.D-5 to 1.D-10 usually adequate.  (Mathematical
   note: the thing that makes this also work is the fact that X'X
   (or XX') is always positive definite when composed of real numbers.
      SUBROUTINE  RECIP(X,Y,NR,NC,IZERO)
      SUBROUTINE XRECIP(NRMX,NRMY,X,Y,NR,NC,IZERO)
 TITLE: RECIProcal
 PURPOSE: To find the reciprocal of every element of a matrix.
 SYMBOLS:
   NRMX,NRMY = row dimensions of X and Y, respectively, in calling
     program
   X = input matrix
   Y = output matrix containing the reciprocals of X
   NR,NC = # of rows and cols of X and Y
   IZERO = 0 to leave 0 elements as 0's;
           1 to terminate processing if a 0 element is found
   DABS = FORTRAN supplied double precision absolute value function
 PROCEDURE:
   Every element of X is tested against 10**-70 as the criteria of
   functional zero.  If the element is not functionally zero, the 
   reciprocal of the element is stored in the corresponding element of Y.   
   If the elementis functionally zero, it is set to 0 or processing is,
   terminated, depending on the input value of IZERO.
   X and Y may be the same matrix in the calling program.
      SUBROUTINE  RNKORD(V,S,IRP,IRO,N)
TITLE: RaNK ORDering of entries
PURPOSE: To put the rank order (low to high) of the entries of 
  a vector (V) in a vector Q.
SYMBOLS:
  V	   is a real vector; input/output
  S	   is a real vector containing on output the entries in V
           sorted from low to high
  IRP      is an integer vector containing on output the rank position
           of entries in  V
  IRO      is an integer vector containing on output entry identifications
           when entries are sorted into ascending order
  N	   is the number of entries in	V ; input/output
PROCEDURE: Sort entries in V into ascending order in S with
  original position in IRO.
      SUBROUTINE  RO2DI(X,Y,N)
      SUBROUTINE XRO2DI(NRMX,NRMY,X,Y,N)
 TITLE: ROw-to-DIagonal matrix
 PURPOSE: To make the first row of an inut matrix into a
   diagonal matrix.
 SYMBOLS:
   NRMX,NRMY = row dimensions of X and Y, respectively, in calling
               program.
   X = input matrix containing the input row in its first row.
   Y = resulting diagonal matrix.
   N = # of cols of X
 PROCEDURE:
   The elements of the first row of X are copied into the main diagonal
   of Y and the remaining elements of Y are zeroed out.  If X is 1 X 1,
   the copy is made and no zeroing occurs.  Note that X must be a 1 X 1,
   matrix, a scalar, or a 1 dimensional array in the calling program
   X cannot be a 2 dimensional array.  X and Y may be the same matrix
   in calling program.
      SUBROUTINE  RO2MT(X,Y,NR,NC)
      SUBROUTINE XRO2MT(NRMX,NRMY,X,Y,NR,NC)
 TITLE: ROw to MaTrix
 PURPOSE: To copy first row of an input matrix into a full matrix.
 SYMBOLS:
   NRMX,NRMY = row dimensions of X and Y, respectively, in calling program
   X = input matrix to be repeated
   Y = output matrix with each row equal to X
   NR = # of rows to be generated in Y
   NC = # of cols of X & Y
 PROCEDURE: The first row of X is copied into Y until Y has NR rows.
   X & Y may be the same matrix in the calling program.
      SUBROUTINE  SCAML(C,X,Y,NR,NC)
      SUBROUTINE XSCAML(NRMX,NRMY,C,X,Y,NR,NC)
 TITLE: SCAlar MuLtiplication
 PURPOSE: To multiply every element of a matrix by a scalar.
 SYMBOLS:
   NRMX,NRMY = row dimensions of X and Y, respectively, in calling program
   C = a real type scalar.  Note - C must be a real number, not an integer
   X = input matrix
   Y = output matrix containing C*X
   NR = # of rows of X
   NC = # of cols of X
 PROCEDURE:
   Every element of X is multiplied by C and the result is stored in Y.
   X and Y may be the same matrix in the calling program.  C must be a
   real number; if C has an integer value, a decimal point must be 
   present in order to indicate that C is a real type number.
   Example: use '1.' rather than just '1'. Also note that C may be a
   matrix, in which case the final element of the matrix is used as C.
      SUBROUTINE  SIMLN(X,Y,NR,NC)
      SUBROUTINE XSIMLN(NRMX,NRMY,X,Y,NR,NC)
 TITLE: solution of SIMultaneous LiNear equations
 PURPOSE: To find the solution set of several linear equations in several
   unknowns.
 SYMBOLS:
   NRMX,NRMY = row dimensions of X and Y, respectively, in calling program
   X = input matrix of the form (A:c) where A is NR X NR, C is NR X N1
   Y = resulting solution for Y in the equation A*Y = C
   NR = # of rows of X; # of rows & cols of A
   NC = # of cols of X
   QGELGZ = number cruncher routine
   REORD = a reordering routine for the result matrix
 PROCEDURE: The order of X is checked for appropriateness, namely the
   # of rows cannot be 1 or greater than or equal to the # of cols.
   The # of cols in the C portion of X is computed, the solution for Y
   in the equation A*Y = C is found and stored in the upper left corner
   of Y as a matrix the same order as C.   A and C are concatenated 
   (A:B), and input to this routine as X.  X and Y may be the same 
   matrix in the calling program.
      SUBROUTINE  SQRT(X,Y,NR,NC,EPS)
      SUBROUTINE XSQRT(NRMX,NRMY,X,Y,NR,NC,EPS)
 TITLE: SQuare RooT
 PURPOSE: To find the square root of all elements of a matrix
 SYMBOLS:
   NRMX,NRMY = row dimensions of X and Y, respectively, in calling program
   X = input matrix
   Y = output matrix containing square roots of X
   NR,NC = # of rows and cols of X and Y
   EPS = lower bound tolerance for 0. values (a negative number)
   DSQRT = FORTRAN supplied double precision square root function
 PROCEDURE:
   If any element of X is negative(less than EPS), an appropriate 
   error message will be written out and execution is terminated.
   Otherwise, the square root of every element of X is stored
   in the corresponding element of Y.  X and Y may be the same 
   matrix in the calling program.  Any element between 0 and
   EPS is set to 0.
      SUBROUTINE  SUB(X,Y,Z,NR,NC)
      SUBROUTINE XSUB(NRMX,NRMY,NRMZ,X,Y,Z,NR,NC)
 TITLE: SUBtraction of 2 matrices
 PURPOSE: To subtract 2 matrices.
 SYMBOLS:
   NRMX,NRMY,NRMZ = row dimensions of X, Y, and Z, respectively,
     in calling program
   X = 1st input matrix
   Y = 2nd input matrix
   Z = resultant matrix
   NR,NC = # of rows and columns, respectively, of X, Y, and Z
 PROCEDURE: Corresponding elements of X and Y are subtracted
   and the result is stored in the corresponding elements of Z.  Any of
   X, Y, or Z may be the same matrix in the calling program.
      SUBROUTINE  TINPT(XIN,X)
      SUBROUTINE XTINPT(NRMX,XIN,X)
TITLE: Temporary data INPuT from data unit 4
PURPOSE: To input data from temporary data unit 4.
SYMBOLS:
  XIN = Name of data file which contains the data to be READ.
  X = Matrix in which to store input data.
  FMT = Array for containing input format (read from the first card of
   data deck).
  C = Labeled common areas containing NMAT, and DSRN.
  NRMX = Dimensionality of X.
  NMAT = Matrix counter.
  QINCR = Program number counter subroutine.
  NR = Number of rows input to X.
  NC = Number of cols input to X.
  XDIMCH = Dimensionality check subroutine.
PROCEDURE: This subroutine opens up unit 4 a file which is CALLED
  'XIN', then reads the header record from the data file,
  increments the matrix counter, prints out NMAT,DSRN,NR,NC,& FMT for
  user, and, if all is ok, reads the actual data file.  If, for any 
  reason the end of the file is encountered before expected, 
  processing terminate after printing an appropriate message.
      SUBROUTINE  TITLE(CH,N)
 TITLE: TITLE printer
 PURPOSE: To print a character string specified by the user.
 SYMBOLS:
   CH = a character string of length N
   N = the length of CH including blanks and punctuation.  Note that 
     there is no limit on the length of CH.
      SUBROUTINE  TRACE(X,N,T)
      SUBROUTINE XTRACE(NRMX,X,N,T)
 TITLE: TRACE of a matrix
 PURPOSE: To find the trace (sum of the main diagonal elements) of a matrix
 SYMBOLS:
   NRMX = row dimension of X in calling program
   X = input matrix
   N = # of cols of X or # of rows of X, which ever is smaller
   T = trace of X, a scalar
 PROCEDURE: The main diagonal elements of X are summed and stored in T.
   Note that T may be a matrix in the calling program, in which case, T
   is returned in the first element of the matrix.
      SUBROUTINE  TRNSP(X,Y,NR,NC)
      SUBROUTINE XTRNSP(NRMX,NRMY,X,Y,NR,NC)
 TITLE: matrix TRaNSPose
 PURPOSE: To transpose a matrix.
 SYMBOLS:
   NRMX,NRMY = row dimensions of X and Y, respectively, in calling program
   X = matrix to be trnsposed
   Y = the transpose of X
   NR = # of rows of X = # of cols or Y
   NC = # of cols of X = # of rows or Y
 PROCEDURE:
   Row and column subscripts in X are interchanged while copying
   the data from X to Y.  X and Y must be different matrices in the 
   calling program.
      SUBROUTINE  VRTCL(X,Y,Z,NRX,NRY,NC)
      SUBROUTINE XVRTCL(NRMX,NRMY,NRMZ,X,Y,Z,NRX,NRY,NC)
 TITLE: VeRTiCaL concatenate
 PURPOSE: To concatenate the cols of 2 matrices to create a super matrix.
 SYMBOLS:
   NRMX,NRMY,NRMZ = row dimensions of X, Y, and Z, respectively,
     in calling program
   X = 1st input matrix
   Y = 2nd input matrix
   Z = resulting super matrix = (X..Y)
   NRY = # of rows in Y
   NRX = # of rows in X
   NC = # of cols in both X and Y
 PROCEDURE:
   A matrix Z' = (X':Y') is created by storing X in cols 1 thru NC,
   rows 1 thru NRX of Z and then by storing Y in cols 1 thru NC and row
   (NRX+1) thru NR of Z.  Thus Z will be NR X NC.  All 3 matrices, X,Y,
   and Z may be the same matrix in the calling program; matrices X 
   and Y may be same but different from Z; matrices X and Z may be the
   same but different from Y: however, matrices Y and Z may not be the
   same if X is different in the calling program.
      SUBROUTINE  XTX(X,Y,NR,NC)
      SUBROUTINE XXTX(NRMX,NRMY,X,Y,NR,NC)
 TITLE: matrix cross-products
 PURPOSE: To find a fast solution to the pre-multiplication of a matrix
   transpose.
 SYMBOLS:
   NRMX,NRMY = row dimensions of X and Y, respectively, in calling program
   X = input matrix
   Y = result matrix = X'X
   NR = # of rows of X
   NC = # of cols of X = # of rows and cols of Y, must be .LT. NRMX & NRMY
 PROCEDURE:
   The first col of X is cross multiplied by each successive column 
   and the results are stored in X.  X is copied into the lower
   triangular portion of a col. of Y.  Then the lower triangular
   portion of Y is copied into the upper triangular portion of Y
   up to the current column.  This process is repeated until the 
   product is completed.  
      SUBROUTINE XXT(X,Y,NR,NC)
      SUBROUTINE XXXT(NRMX,NRMY,X,Y,NR,NC)
 TITLE: MATRIX PRODUCT
 PURPOSE; To obtain matrix product Y = XX'.
 SYMBOLS:
     NRMX,NRMY = number of rows of X and Y, respectively, in 
       calling program
     X	  = data matrix; input/output.
     Y	  = product matrix; ouput.
     NR	  = number of rows of X and order of Y; input/output.
     NC	  = number of columns of X; input/output.
 PROCEDURE: The product matrix Y is to be zeroed.
   The row elements of X is multiplied by its elements of row
   up to the current column, and stored the value in Y.
   The process is continued until the product is completed.
   X and Y may not be the same matrix.




                      General Statistical - Monte Carlo Programs
      FUNCTION CHIVLT(P,NDF)
 TITLE: VaLue of CHI-square
 PURPOSE: To obtain the value of Chi square, with NDF degrees of
   freedom below which P percent of the area falls. 
 SYMBOLS:
   P   = proportion
   NDF = degree of freedom
 PROCEDURE:
   Uses  hyperbolic interpolation, and function CHIPRA used.
      SUBROUTINE  DAPPR(A,RNP,PJ,APP,G,VEC,S1,S2,HHH,NAT,NF,MCASE,
     1MPJ,MAPP,MPRIN,PCON,WCON,COSLT,NCYCLT)
      SUBROUTINE XDAPPR(NRMA,NRMRNP,NRMPJ,NCMPJ,NRMAPP,NRMG,NCMG,NRMVEC,
     2 NRMS1,NCMS1,NRMS2,NCMS2,A,RNP,PJ,APP,G,VEC,S1,S2,HHH,NAT,
     3 NF,MCASE,MPJ,MAPP,MPRIN,PCON,WCON,COSLT,NCYCLT) 
 TITLE: Direct method with APP weighting factor Rotation
 NOTE: APP is artificial personal probability of being in a hyperplane
   INPUT or:  APP = 1/(1 + ACON*(ABS(B) TO POWER PCON)
   Uses FORMAL subroutines: TITLE,INVRS,MATML,XTX,GENR8,SIMLN,COPY
   Other subroutines used: XMATML, XPRINT, XXTX, VARMX
 SYMBOLS:
   NRMA,NRMRNP,NRMPJ,NRMAPP,NRMG,NRMVEC,NRMS1,NRMS2 
          is the number of rows of A, RNP, PJ, APP, G, VEC,
          S1, and S2, respectively, in calling program
   NCMG,NCMS1,NCMS2  is the number of columns of G, S1,
          and S2, respectively, in calling program
   A	  is original factor matrix: input/output
   RNP    is FORMAL matrix N' of hyperplane normals (columns): output
   PJ	  is matrix P during solution, is converted to B for output
	      see note concerning MPJ
   APP    is matrix of APP: output
	      see note concerning MAPP
   G	  is formal matrix (NN') inverse during solution,
	  is converted to PHI for output
   VEC    is formal matrix having row vectors as follows:
	      1   mean absolute PJ for each rotated factor
	      2   sum APP*(PJ**2) over attributes for each rotated factor
	      3   shift vector in rotation of each factor in turn
	      4   square roots of diagonal G
	      5-I1  scratch
   S1 and S2 are scratch FORMAL matrices
          on output S1 contains correlation of items and factors,
          S2 contains T' 
   HHH    is a scratch vector of length I1
   NAT    is number of attributes: input/output
   NF	  is number of factors: input/output
   MCASE  is case of APP function: input/output
          must be specified when APP is to be iterated
	      = 1 for one sided function
	      = 2 for two sided function
   MPJ    is integer flag for initial weights hypothesized: input/output
	      = 0 for initial weights obtained by a normal VARIMAX rotation
	      = 1 for initial hypothesis factor weights input in PJ
	      default = 0
   MAPP   is integer flag concerning iteration of APP
	      = 0 for APP iteration
	      = 1 for fixed APP as input in APP
	      default = 0
   MPRIN  is integer flag concerning printed output
	      = 0 for no printed output
	      = 1 for printing matrices P and APP for each cycle
	      default = 0
   PCON   is power constant in APP function:input/output
	      default = 10.0
   WCON   is proportion of mean absolute factor weight for  APP = 0.5:
	      input/output; default = 0.5
   COSLT  is limit for COS between successive trial normals in test for
	      convergence: input/output; default = 0.99999
   NCYCLT is limit on number of iterative cycles: input/output
	      default = 20
      SUBROUTINE  DISCR(SH,SD,GPM,F,G,CV,TM1,TM2,NAT,NG,NDFH,NDFD,MOD,
     1NPS,NPE,ISWP)
      SUBROUTINE XDISCR(NRMSH,NRMSD,NRMGPM,NRMF,NRMG,NRMCV,NRMTM1,
     2NRMTM2,SH,SD,GPM,F,G,CV,TM1,TM2,NAT,NG,NDFH,NDFD,MOD,NPS,
     3ISWP) 
 TITLE: DISCRiminant analysis
 PURPOSE: To compute discriminant composite solution.
   solution space restricted to  NP  principal axes of estimated
   total covariance matrix, THETA, derived from components of
   covariance.
 NOTE: uses FORMAL; matrices must be at least (10,10), if NAT is
   larger than 9, the matrices must be at least (NAT+1,NAT+1)
 SYMBOLS:
   NRMSH,NRMSD,NRMGPM,NRMF,NRMG,NRMCV,NRMTM1,NRMTM2
           is the number of rows of SH, SD, GPM, F, G, CV, TM1,
           and TM2, respectively, in calling program
 MATRICES:
   SH    is hypothesis effect SSCP matrix; input/output
   SD    is denominator effect SSCP matrix; input/output
   GPM   is group means matrix (groups by hypothesis effect);
         input/output
         row for each group plus one row for total mean
   F, G, CV, TM1, TM2 are storage matrices
   NAT   is number of attributes
   NG    is number of groups (elements in hypothesis effect)
   NDFH  is number of degrees of freedom for hypothesis
   NDFD  is number of degrees of freedom for denominator
   MOD   is model identification number
         = 1 when total sample is divided into samples of
  	   subpopulations in proportion to relative frequencies
  	   of the subpopulations in the total population
         = 2 when the total sample is drawn at random from the total
  	   population
   NPS   is starting number of principal axes of estimated THETA
   NPE   is ending number of principal axes of estimated THETA
         separate solutions are obtained using  NP  principal axes
         of estimated THE from NPS  to  NPE  in steps of 1
         working values, NPSW and NPEW, of NPS and NPE are adjusted
         when necessary as follows:
         let NPAM be the maximum usable number of principal axes
         if NPS equals 0 (or less), NPSW is set equal to NPAM
         if NPS is greater than NPAM, NPSW is set equal to NPAM
         if NPE is less than NPSW, NPEW is set equal to NPSW
         if NPE is greater than NPAM, NPEW is set equal to NPAM
   PRINTING CONTROL:
     ISWP  =  0  for standard results
           =  1  for additional information:
   	      matrices  THETA-TILDA , L , F , G , A
     columns of CV are used as follows
	 1  =  1/D
	 2  =  1/SQRT(DELTA)
	 3  =  PHI
	 4  =  1/SQRT(1 - PHI)
	 5  =  LAMBDA
	 6  =  LN(1 + LAMBDA)
	 7  =  V STATISTIC
	 8  =  SQRT(SD(J,J)/NDFD)
      SUBROUTINE  FNFRT(A,NMP,TP,PHI,PJT,B,Q,NAT,NCF,ICN,ICT,ICP,
     1ICB)
      SUBROUTINE XFNFRT(NRMA,NRMNMP,NRMTP,NRMPHI,NCMPHI,NRMPJT,NRMB,
     2 NRMQ,NCMQ,A,NMP,TP,PHI,PJT,B,Q,NAT,NCF,ICN,ICT,ICP,ICB) 
 TITLE:  FiNal Factor RoTation
 PURPOSE: Performs final computations for factor analysis rotation of 
   matrix, results are printed and passed to main program.  
 NOTE: Number of common factors must be 2 or greater and at most 1 less
   than the matrix dimensions in FORMAL.  Uses FORMAL. 
 SYMBOLS:
   Input-output matrices
	A    = coordinate axes factor matrix
	NMP  = a double precision matrix for normals, N'
	     input form controlled by  ICN
	     ICN = 0	for N' input (normals are columns)
	   = 1	for N input (normals are rows)
	   = 2	for N' to be computed from T'
	     normals are assumed to be unit vectors when input
	     output  N'	(normals are columns)
	TP   = matrix for primary vectors
	     input form controlled by  ICT
	     ICT = 0 for T' input (primary vectors are columns)
	  = 1	for T input (primary vectors are rows)
	  = 2	for T' to be computed from N'
	     Primary vectors are assumed to be unit vectors when input
	     Output T' (primary vectors are columns)
	            Note: Either NMP or TP or both must be input
	PHI  = Factor intercorrelation matrix on output = TP'*TP
	     row NCF + 1 contains diagonal entries of matrix D
        D    = (Diag(PHI-INVERSE))**-.5
	PJT  = matrix for projections on normals = A*NMP
	     Input form controlled by  ICP
            ICP = 0 for PJT input
	  = 1 for PJT to be calculated
     PJT is output
	B    = factor weight matrix = PJT*D-INVERSE
     Input form controlled by  ICB
            ICB = 0 for bb input
	  = 1 for B to be computed
     B is output
	Q    = matrix of covariances of attributes with factors= B*PHI
	      Q is output
     Input scalars
         NRMA,NRMNMP,NRMTP,NRMPHI,NRMPJT,NRMB,NRMQ  = the number of rows
           of A,NMP,TP,PHI,PJT,B, and Q, respectively, in calling program
         NCMPHI,NCMQ = the number of columns of PHI and Q, respectively,
           in calling program
	 NAT  = number of attributes
	 NCF   = number of common factors
   NOTE: Control flags described above	 ICN, ICT, ICP, ICB
      SUBROUTINE  FAHIR(C,CMV,A,COM,NAT,NCF)
      SUBROUTINE XFAHIR(NRMC,NCMC,NRMCMV,NRMA,NCMA,NRMCOM,C,CMV,
     2 A,COM,NAT,NCF)
   TITLE:  Factor Analysis of C With highest covariance in diagonals
   PURPOSE: To perform factor analysis of C.  Runs on FORMAL.
   NOTE: Covariance matrix  C  is input/output in  C.
     Dimensionality must be less than  I1. Runs on FORMAL.
   SYMBOLS:
     NRMC,NRMCMV,NRMA,NRMCOM are the number of rows of C, CMV,
       A, and COM, respectively, in calling program
     NCMC,NCMA are the number of col. of C and A, respectively,
       in calling program
     CMV, A, COM are scratch on input
       CMV on output contains C with highest covariances on the diagonal
       Note: If  C  is not desired on output, C and CMV can be identical
       A	on output contains the factor matrix
       COM	on output contains communalities for each given number of factors
     NAT  is number of attributes; input/output
     NCF  is maximum number of factors to be obtained; if NCF is greater than
  	  the number of positive eigenvalues of C with highest covariances on the
  	  diagonal cells, NCF will be reduced to the number of positive
  	  eigenvalues.	NCF must be a symbolic representation of an integer.
      SUBROUTINE  FAWMV(C,CMV,A,COM,NAT,NCF)
      SUBROUTINE XFAWMV(NRMC,NCMC,NRMCMV,NRMA,NCMA,NRMCOM,C,CMV,
     2 A,COM,NAT,NCF)
   TITLE:  Factor Analysis of C With MV in diagonals
   PURPOSE: To perform factor analysis of C.
     runs on FORMAL
   NOTE: Covariance matrix  C  is input/output in  C.
     Dimensionality must be less than  I1. Runs on FORMAL.
     If matrix C is singular then the highest covariance is used
     as the communality estimate (See FAHIR).
   SYMBOLS:
     NRMC,NRMCMV,NRMA,NRMCOM are the number of rows of C, CMV,
       A, and COM, respectively, in calling program
     NCMC,NCMA are the number of col. of C and A, respectively,
       in calling program
     CMV, A, COM are scratch on input
       CMV on output contains C with multiple variances in diagonal
       Note: If  C  is not desired on output, C and CMV can be identical
       A	on output contains the factor matrix
       COM	on output contains communalities for each given number of factors
     NAT  is number of attributes; input/output
     NCF  is maximum number of factors to be obtained; if NCF is greater than
  	  the number of positive eigenvalues of C with multiple variances in the
  	  diagonal cells, NCF will be reduced to the number of positive
  	  eigenvalues.	NCF must be a symbolic representation of an integer.
      SUBROUTINE  ITPFA(C,A,HV,S1,S2,NAT,NF)
      SUBROUTINE XITPFA(NRMC,NRMA,NRMHV,NRMS1,NCMS1,NRMS2,NCMS2,
     2 C,A,HV,S1,S2,NAT,NF) 
   TITLE:  ITerative Principal axes Factor Analysis
   Note: This program does not have a protection against heywood cases
     Runs on FORMAL
 SYMBOLS:
   NRMC,NRMA,NRMHV,NRMS1,NRMS2 is the number of rows of C, A, HV, S1,
     and S2, respectively, in calling program
   NCMS1,NCMS2 is the number of col. of S1 and S2, respectively, in 
     calling program
   C	  is covariance matrix; input/output
   A	  is resulting factor matrix;output
   HV	  is a formal matrix with columns used as vectors:
    	  1   for communalities of previous trial,
 	  2   for communalities computed from factor matrix,
 	  3   for differences between input communalities and computed
               communalities for previous trial,
      	  4   for differences for present trial.
   S1	  contains the eigenvalues for the present trial in the diagonal; output
   S2	  contains the covariance matrix with the obtained communalities; output
   NAT    is the number of attributes; input/output
   NF	  is the number of factors; input/output
      SUBROUTINE  LSQHY(A,S,N,W,D,C,NR,NC,IFL)
      SUBROUTINE XLSQHY(NRMA,NRMS,NRMN,NRMW,NRMD,NRMC,A,S,
     2 N,W,D,C,NR,NC,IFL) 
 TITLE:  LeaSt SQuares HYperplane fit
 PURPOSE: To rotate a factor matrix to the best (least-squares)
   hyperplane fit using markers specified by the user.
 SYMBOLS:
   NRMA, NRMS, NRMN, NRMW, NRMD, NRMC = dimensions of input matrices
     in calling program,must be .GE. 10 
   A = input factor matrix
   S = selection matrix (all 1's and 0's)
   N = output factor normals matrix - transposed
   W,D,C = temporary work files
   NR = number of rows of A and S
   NC = number of cols of A and S, and dimensions of B
   IFL = o to skip intermediate printing,
         1 to print intermediate output
   QINCR = program counter subroutine
   XDIMCH = dimensionality check subroutine
   QLRSVM = implicit QL eigen routine
 PROCEDURE: According to Ledyard R Tucker,  none of the input matrices
   may be the same in the calling program.
      SUBROUTINE  MISNR (XDA,XBA,XVA,CXY,RXY,XMI,T,TV,XN,ZXY,
     2 NS,NV,IOPT,IONE)
      SUBROUTINE XMISNR (NRMXDA,NRMXBA,NRMXVA,NRMCXY,NRMRXY,NRMXM1,NRMT,
     2 NRMTV,NRMXN,NRMZXY,XDAT,XBAR,XVAR,CXY,RXY,XMISS,T,TV,XN,
     2 ZXY,NS,NV,IOPT,IONE)
 TITLE:  MISsiNg data coRrelation matrix
 PURPOSE:  To compute correlation matrix and vector of means
   and variances, when missing data is present.
 OPTIONS:

   1     compute vector of means and sample sizes only
         (output to diagonal of XBA). All other matrices zeroed.
   2  	 No missing data, will not check (most efficient
         if there is no missing data; will produce an
	 error if any data is missing).
   3	 Replace missing data by mean.
   4	 Deletion of all data for any subject who has
         any missing data present.
   5	 Compute covariances based on available pairs of
         data(default).  Means and covariances will be based on
         available pairwise data for that pair of variables
	 only.	Correlations will be based on these covariances
	 and variances from diagonal of cov. matrix.
   6	 Compute pairwise covariance as in option 5 above.
         Estimates for variances are based on available
	 pairwise data for that pair of variables only.
   7	 Covariances based on available pairs, however the
         full data means will be utilized.  Correlations
	 will be based on full data variances (diag of cov.
	 matrix).
   8	 Covariance matrix computed as in option 7 above.
         Variance for correlation matrix computed from available
	 pairs of data when full data mean is utilized.

VARIABLES USED
  NRMXDA,NRMXBA,NRMXVA,NRMCXY,NRMRXY,NRMXM1,NRMT,NRMTV,NRMXN,
        NRMZXY are the number of rows of XDA, XBA, XVA, CXY, RXY, XM1,
        T, TV, XN, and ZXY, respectively, in calling program
  XDA	is rectangular matrix containing raw data (dimensions input,
        NS - number of subjects  X,  NV - number of variables)
  XBA	is matrix of means of variables.  If options 5 or 6 are
        output, used this matrix will be mean for row variable for
        all cases when there was data for this variable and column variable.
  XVA	is matrix of variances of variables.  If options 6 or 8 are
        output, used this matrix will be variance for column variable for
        all cases when there was data for this variable and row variable.
  CXY	is a matrix of covariances.
        output
  RXY	is a matrix of correlations.
        output
  XMI	is a vector of length NV contains information as to how missing
        input data were coded (per variable). Missing data for any
	variable is identified by coding it as in TV(I).
	this subroutine  can not differentiate between -0.0 (or
	blanks) and 0.0 (or zero).
	If this procedure cannot describe how missing
	data is coded in the users study, write
	a new subroutine   SUBMIS  (see that subroutine
	for its parameters).
  T	a dummy double precision vector of length NV.
  TV 	A dummy double precision vector of length nv.
  XN 	A dummy matrix of size (NV X NV).
  ZXY	Matrix of size NV X NV containing sample size for
        output	 pairwise covariances.
  NS	number of subjects
        input
  NV	number of variables
        input
  IOPT	Option on how to handle missing data (see above)
        input
  IONE	Denominator of variances and covariances.
	input	if ione: -1 then denominator is N - 1 (unbiased),
			  0 then denominator is N (maximum likelihood), or
			  1 then denominator is N + 1 (minimum error).


      SUBROUTINE  PRDIF(V,N,MD)
 TITLE: PRints a column of values and first DIFferences
    these values may be either entries in a vector or diagonal entries of a
    matrix.	A common application is printing a series of eigenvalues.
 SYMBOLS:
    V	  is array containing values to be printed; input/output
    N	  is number of values; input/output
    MD  is a flag concerning location of values in V; input/output
	   = 0 for vector (includes a column of a two dimensional array)
	   = number of rows in main of array V when the diagonal entries
	     of V are to be printed.
      SUBROUTINE  PSTCR(P,Q,C,R,NIN,NAT,MPRN)
      SUBROUTINE XPSTCR(NRMP,NRMQ,NRMC,NRMR,P,Q,C,R,NIN,NAT,MPRN)
   PURPOSE: To compute means, standard deviations, covariance matrix, 
     correlation matrix from product matrix and attribute sums
     row dimension of arrays for P, Q, C, and R is I1 in common block B
     as in FORMAL, must be at least NAT+2
     column dimension of arrays for P, Q, C, and R must be at least 
     NAT+2.
   SYMBOLS:
      NRMP,NRMQ,NRMC,NRMR are the number of rows of P, Q, C, and R,
         respectively, in calling program
      P contains product matrix, input
      row  NAT + 1  contains attribute sums
      Q contains within sample SSCP matrix, row NAT+1 contains means; 
         output
      C to contain covariance matrix, output
         row NAT+1 contains means
         row NAT+2 contains standard deviations
         may use same array as P
      R to contain correlation matrix, output
         may use same array as C (C not output) or P
      NIN is number of individuals in sample
      NAT is number of attributes
      MPRN is an integer indicator for printing output
         = 0  for no printing;  = 1  for printed output
      SUBROUTINE  RANDQ(Q,CP,F,S,NIN,NAT,NCF,IRN,MF)
      SUBROUTINE XRANDQ(NRMQ,NRMCP,NRMF,NRMS,Q,CP,F,S,NIN,NAT,NCF,
     2 IRN,MF) 
   PURPOSE: To generate a sample SSCP matrix and mean vector
     The SSCP matrix contains sums of squares and cross products of
     measures deviated from sample means
     The population is multidimensional normal with a given mean vector
     and covariance matrix
   MATRICES;  must be at least NAT+1.
      Q	     contains sample SSCP matrix in first NAT rows and
   	     sample mean vector in row	NAT+1; output
      CP     contains population covariance matrix in first NAT rows and
   	     population mean vector in row  NAT+1; input/output
      F	     contanins factor matrix of population covariance matrix;
   	     may be input, otherwise is obtained by a cholesky
   	     decomposition of CP; output
      S	     is a scratch matrix
   SCALARS:
      NRMQ,NRMCP,NRMF,NRMS is the number of rows of Q, CP, 
             F, and S, respectively, in calling program
      NIN    is sample size; input/output
      NAT    is number of attributes; input/output
      NCF    is number of factors in  F ; input when  F  is input; output
      IRN    is a 15 digit integer ending in  1  used in generating random
   	     values; input, changed before output
      NCF    is number of factors in F ; input when F is input ; output
   	     = 0  when	F  is not input
   	     = 1  when	F  is input,
   	     on output	MF = 1
      SUBROUTINE  RNDCR(CS,RS,CP,F,NIN,NAT,NCF,IRN,MF)
      SUBROUTINE XRNDCR(NRMCS,NRMRS,NRMCP,NRMF,CS,RS,CP,
     2 F,NIN,NAT,NCF,IRN,MF) 
   PURPOSE:
     To generate a covariance matrix and correlation matrix among NAT 
     attributes for a sample size NIN drawn from a population 
     distributed multidimensional normal with given covariance matrix 
     and mean vector.
   MATRICES;  must be at least  NAT+2.
     CS  contains sample covariance matrix with sample means in row NAT+1
         and sample standard deviations in row NAT+2; output
     RS  contains sample correlation matrix; output
     CP  contains population covariance matrix in first  NAT  rows and
         population mean vector in row  NAT+1 ; input/output
     F   contanins factor matrix of population covariance matrix;
         may be input, otherwise is obtained by a cholesky
         decomposition of  CP ; output
   SCALARS:
     NRMCS,NRMRS,NRMCP,NRMF   is the number of rows of CS, RS,
           CP, and F, respectively, in calling program
     NIN   is sample size; input/output
     NAT   is number of attributes; input/output
     NCF   is number of factors in  F ; input when  F  is input; output
     IRN   is a 9 digit integer used in generating random
           values; input, changed before output
     MF    is a flag concerning matrix  F ; input:
             = 0  when F is not input
             = 1  when F is input,
           on output MF = 1
     RAN   is the subroutine  which furnishes the random numbers.
      SUBROUTINE  RNDQM(Q,NIN,NAT,IRN)
      SUBROUTINE XRNDQM(NRMQ,Q,NIN,NAT,IRN)
   PURPOSE: To generate a sample SSCP matrix and mean vector
     The SSCP matrix contains sums of squares and cross products of
     measures deviated from sample means
     population is multidimensional normal with mean vector = 0 and
     covariance matrix = I.
   SYMBOLS:
      NRMQ  is the number of row of Q in calling program, must be at
            least NAT + 1.
      Q     contains sample SSCP matrix in first NAT rows and
            sample mean vector in row NAT+1 ; output
      NIN   is sample size; input/output
      NAT   is number of attributes; input/output
      IRN   is a 15 digit integer ending in  1  used in generating random
            values; input, changed before output
      SUBROUTINE  RNDVC(RVEC,ISRT,N,IRN1,IRN2)
   PURPOSE: To generate a vector of random normal deviates, mean = 0,
     s.d. = 1 randomly sorted
   SYMBOLS:
        RVEC  contains the random normal deviates on output
        ISRT  is a vector of random integers used in sorting
        N     is the number of entries in  RVEC  and  ISRT
        IRN1  and  IRN2  are two	9  digit, odd integers used in
   	      generating the random deviates and integers
   	      their initial values are to be input, their values are
   	      changed randomly before output
      SUBROUTINE  RXTX(P,NIN,NAT)
      SUBROUTINE XRXTX(NRMP,P,MIN,NAT)
 TITLE: Read data from cards then compute XTX.
 PURPOSE: For a matrix X read from cards computes product matrix P = X'X
   column sums of X are in row NAT+1 of P
   row dimensions of array for P is NRMP and must be at least NAT + 1.
   column dimension of array for P must be at least NAT + 1
 SYMBOLS:
    NRMP    is number of row of P in calling program
    P	      is matrix containing results; output
    NIN     is number of objects (rows of X); read from header card for
	      data deck; output
    NAT     is number of attributes (columns of X); read from header card
	      for data deck; output
    X       is the matrix to be contained in a deck of data cards 
              with a record of one or more cards for each object 
              (row of X) a header card is to precede the data deck
              columns	1 -  4: NIN
              columns	5 -  8: NAT
              columns	9 - 80: format of data cards in parentheses
   	      data title may follow the format and will be printed
      SUBROUTINE  SUBMIS (NRMXN,T,ACC,XMISS,XN,NV)
TITLE:    SUBject's MISsing data (subroutine of MISNR)
PURPOSE:  To indicate which data points are missing.
  This subroutine  reads in a single subject's data - in a
  vector called T, sees whether it is equal to XMISS - a
  vector telling the program how missing data was coded -
  then returns a value of 1.0 in matrix XN(I,J) if the data
  for that pair of variables is good, or 0.0 if that pair
  had missing data.
SYMBOLS:
     NRMXN = number of rows dimensioned for XN in main
     T     = vector of data for a single subject
     ACC   = how small a difference between observed and missing
             value is  before considered missing
             (set in main to be [0.00001]). Related to 
             accuracy of computer.
     XMISS = vector containing elements indicating how missing data
             were coded per each variable.
     XN    = matrix containing elements of either zero (ignore this
             pair for this subject) or one (retain data for this
             pair for this subject).
     NV    = number of variables  correlated.
      SUBROUTINE VARMX(A,F,H,NV,NF)
      SUBROUTINE XVARMX(NRMA,NRMF,NCMF,A,F,NAT,NF,H,MRAW)
   PURPOSE: To perform VARIMAX rotation using subroutine QVARMX
   SYMBOLS:
      A	   is original loading matrix; input/output
      F	   is rotated loading matrix; output
      NAT  is number of attributes; input/output
      NF   is the number of factors; input/output
      H	   is vector containing communalities of attributes; output
   	      must have IZ elements
      MRAW is an integer concerning case of VARIMAX; input/output
   	      = 0  for unweighted (raw) solution
   	      = 1  for normalized solution
      SUBROUTINE  TSAMX (X,C,F,NIN,NAT,NCF,IRN,MF)
      SUBROUTINE XTSAMX (NRMX,NRMC,NRMF,X,C,F,NIN,NAT,NCF,IRN,MF)
   PURPOSE:
     To obtain sample data matrix from multidimensional normal population
     having given population mean vector and covariance matrix.
   Matrices:
      X     is sample data matrix on output
	     array to contain  X  must have at least NIN rows and NAT columns
      C     is population covariance matrix with mean vector in row NAT+1;
	     input/output
	     array to contain  C  must have at least NAT+1 rows and columns
      F     is factor matrix of population covariance matrix, may be input,
	     otherwise is obtained by a cholesky decomposition; OUTPUT
	     array to contain  F  must have at least NAT+1 rows and columns
   Scalars:
      NRMX, NRMC, NRMF  are number of rows in main of arrays to contain
	     matrices  X ,  C ,  F ; input/output
      NIN   is number of individuals in sample; input/output
      NAT   is number of attributes; input/output
      NCF   is number of columns of  F ; input if  F  is input; output
      IRN   is a 15 digit integer ending in 1; input, altered before output
      MF    is a flag concerning matrix  F ; input:
	     = 0  if  F  is not input,
	     = 1  if  F  is input
	     is set to 1 for output




                      GENERAL PURPOSE WORK-HORSE PROGRAMS
      SUBROUTINE QDIMCH(N,X,IFLAG)
      SUBROUTINE XDIMCH(ICHECK,N,X,IFLAG)
 TITLE: DIMensionality CHeck
 PURPOSE:
   This program tests a matrix dimension specified by the user and
   the maximum allowable dimension and outputs an appropriate error 
   message.




      FUNCTION QDVNRM (X)
 TITLE: NoRMal DeViate
 PURPOSE:
   This function will input A random uniform deviate and
   transform it into A normal deviate.
   it works via the inverse probability function.
   In essence, if X is an area under A normal curve, then
   QDVNRM will be the normal deviate which cuts the NORMAL
   distribution into this area.




      SUBROUTINE  QGELGZ(NRMB,NRMAA,B,AA,N,NR,X,A,EPS,IER)
 PURPOSE: To solve the simultaneous linear equations with >1 right 
   hand side by Gauss elimination with complete pivoting.




      SUBROUTINE  QIMTQL(N,D,E,IERR)
 TITLE: eigenvalues of tridiagonal matrix
 PURPOSE: 
   To determine eigenvalues of a symmetric tridiagonal matrix
   using the implicit QL method.




      SUBROUTINE  QLNEQT(NRMA,A,B,E,S,D,F,BP,N,MS,IERR,NCYC)
 PURPOSE: To solve for B in linear equation S*B = C
   by an iterative procedure may be used to improve an approximate
   solution.




      SUBROUTINE QLRSM(NRMA,N,A,D,E,IERR)
 TITLE: Latent Roots of a Symmetric Matrix
 PURPOSE: to find all latent roots (eigenvalues) of a symmetric matrix
   by the implicit QL method.




      SUBROUTINE  QLRSVM(NRMA,NRMZ,NRMRV,N,NRV,A,Z,W,D,E,E2,IND,RV,IERR)
TITLE: Latent Roots and Some Vectors of a symmetric Matrix.
PURPOSE: Find all latent roots (eigenvalues) and a smaller number of
  latent vectors (eigenvectors) of a symmetric matrix.




      SUBROUTINE QMILU(NRMA,NRMIDU,A,N,IDUM,V)
 PURPOSE: To find the inverse of the matrix whose LU decomposition 
   is furnished in A. 




      SUBROUTINE  QMINVZ(NRMAA,NRMA,NRMIDU,AA,N,DET,IDET,A,IDUM,IER,V)
 PURPOSE: To calculate the determinant of input matrix as a by-product
   of computation.  This program the inverse of the N by N real,
   general input matrix by Gauss Elimination.




      FUNCTION QRAND1(T)
 PURPOSE:
   To provide a pseudo random normal deviate (MEAN = 0, SD = 1)
   by inverse transformation of A pseudo random P provided by function
   RAN; function QDVNRM is used for the inverse transformation, QDVNRM
   gives a quick approximation to the inverse transformation
   a 'seed' integer is to be given function ran before the first use
   of RAN by subroutine  SETRAN(IRN) where IRN is a 9 digit integer.
   the final integer may be obtained by subroutine SAVRAN (IRN).




      SUBROUTINE  QTINVT(NRMZ,N,D,E,E2,M,W,IND,Z,IERR,
     1RV1,RV2,RV3,RV4,RV6)
 TITLE: eigenvectors of a Tridiagonal matrix by INVerse iTeration
 PURPOSE: To determine the eigenvectors of a tridiagonal matrix
   corresponding to a set of ordered eigenvalues, using inverse
   iteration.




      SUBROUTINE  QTRBAK(NRMA,NRMZ,N,A,E,M,Z)
 TITLE: back transformation of eigenvectors
 PURPOSE: To form the eigenvectors of a real symmetric matrix from the
   eigenvectors of that symmetric tridiagonal matrix determined by TRED




      SUBROUTINE  QTRED1 (NRMA,N,A,D,E,E2)
TITLE: reduce a symmetric matrix to tridiagonal form




      SUBROUTINE QVARMX(NRMAV,NCMA,M,K,A,NC,TV,H,KRAW)
 PURPOSE:
   To perform VARIMAX rotation revision of VARMX from IBM SSP.




      SUBROUTINE  QWRGDM(NR,NC)
 TITLE: error message routine for invrs,simln, and det




      SUBROUTINE XREORD(NRMX,X,NR,NC)
 PURPOSE:
   To accompany QGELGZ for purposes of rearranging G so it is 
   in the expected matrix form (i.e. such that AG=B).




      SUBROUTINE QINCR
Title: program counter
Purpose: This subroutine is called by all of the FORMAL
  subroutines in order to keep track of the number of subroutines 
  call by the user.  This number may be printed out when an error 
  condition has occured.  The first time QINCR is called, a FORMAL 
  heading is printed out, and data input and output files are opened
  to defaults of 'FRMLIN.DAT' and 'FRMLOT.DAT', respecively.
Symbols:
  A = labeled common area containing IPGM
  IPGM   = the current number of FORMAL call statements made from the
           main program.
  ISUB   = the number of subroutines used indirectly since
           the last call in main, will be zero at calls from
           main.
  IZINCR = will be 1 at all FORMAL calls made from the
           main program.  2 will be at all calls made one
           subroutine from the main, etc.
Procedure: QINCR adds 1 to IPGM, tests IPGM to see if this is the first
  QINCR and if so, it prints out a heading.




      SUBROUTINE IOPUT (XINPUT, OUTPUT, IFLAG)
TITLE: Input Output controller
Purpose: This subroutine will change the name of the input
  and/or output file to 'XINPUT' and 'OUTPUT', respectively.
Symbols:
  'XINPUT' = name of file containing input data (in the form
              of filename.ext - ext will default to DAT)
  'OUTPUT' = name of file which will contain the output data (in the
             form of filename.ext - ext will default to DAT)
  IFLAG    = flag denoting whether input and/or output will be
              done.
          1  input file - 'XINPUT' only will be opened; a close
              to the old input file will be done if necessary.
          2  output file - 'OUTPUT' only will be opened; a message
              to this newly opened file will be generated;  if any
              FORMAL calls were made previous to the first use of
              this program, then the previous output file (default
              is 'FRMLOT.DAT') will be closed and a message will
              be made where future output can be found; the old
              file will also be closed.
          3  both input and output file will be opened as by options
              1 and 2 above.
  A        = labeled common area containing IPGM.
  C        = labeled common area containing DSNR (name of input device) and NMAT.
  IPGM     = program counter to be incremented by this subroutine.
Procedure: Adds 1 to IPGM as in QINCR, tests to see if this is
  the the first usage, if so defaults will be
  'FRMLIN.DAT' and 'FRMLOT.DAT'.  Files will be opened with the
  appropriate header.