*----------------------------------------------------------------------- * tedium.f - very boring functions * by Andy Allinger, released to the public domain *----------------------------------------------------------------------- *----------------------------------------------------------------------- * test for infinity *----------------------------------------------------------------------- FUNCTION ISINF (X) IMPLICIT NONE REAL X LOGICAL ISINF REAL Y Y = X - X ISINF = Y .NE. Y RETURN END ! of ISINF *----------------------------------------------------------------------- * substitute for LEN_TRIM for old compilers *----------------------------------------------------------------------- FUNCTION LTRIM (STRING) IMPLICIT NONE CHARACTER*(*) STRING INTEGER LTRIM INTEGER J, L L = LEN(STRING) DO 50 J = L, 1, -1 IF (STRING(J:J) .EQ. ' ') GO TO 50 LTRIM = J RETURN 50 CONTINUE LTRIM = 0 RETURN END ! of LTRIM *----------------------------------------------------------------------- * first non-blank position in character variable *----------------------------------------------------------------------- FUNCTION KTRIM (STRING) IMPLICIT NONE CHARACTER*(*) STRING INTEGER KTRIM INTEGER J, L L = LEN(STRING) DO 50 J = 1, L, +1 IF (STRING(J:J) .EQ. ' ') GO TO 50 KTRIM = J RETURN 50 CONTINUE KTRIM = 0 RETURN END ! of KTRIM *----------------------------------------------------------------------- * BIT_SIZE replacement for old compilers *----------------------------------------------------------------------- FUNCTION BITSZ () IMPLICIT NONE INTEGER BITSZ INTEGER I, J J = 0 J = NOT(J) BITSZ = 0 DO I = 1, 256 J = ISHFT(J, -1) IF (.NOT. BTEST(J, 0)) THEN BITSZ = I RETURN END IF END DO END ! function BITSZ *----------------------------------------------------------------------- * Find an avaiable I/O unit *----------------------------------------------------------------------- FUNCTION IOUNIT () IMPLICIT NONE INTEGER IOUNIT INTEGER J, IFAULT LOGICAL ISOPEN, ISXIST IOUNIT = -1 DO 10 J = 1, 99 IF (J .EQ. 5 .OR. J .EQ. 6 .OR. J .EQ. 7) GO TO 10 INQUIRE (UNIT=J, EXIST=ISXIST, OPENED=ISOPEN, IOSTAT=IFAULT) IF (IFAULT .EQ. 0 .AND. ISXIST .AND. .NOT. ISOPEN) THEN IOUNIT = J RETURN ! success END IF 10 CONTINUE RETURN ! failure END ! of IOUNIT *======================================================================= * Convert to upper case *======================================================================= SUBROUTINE UCASE (STRING) IMPLICIT NONE CHARACTER*(*) STRING INTEGER I, J, LTRIM CHARACTER*26 LOWER, UPPER SAVE LOWER, UPPER DATA LOWER / 'abcdefghijklmnopqrstuvwxyz' / DATA UPPER / 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' / DO I = 1, LTRIM(STRING) J = INDEX(LOWER, STRING(I:I)) IF (J .NE. 0) STRING(I:I) = UPPER(J:J) END DO RETURN END ! of UCASE *----------------------------------------------------------------------- * Safe to divide? *----------------------------------------------------------------------- FUNCTION SAFDIV (NUMER, DENOM) IMPLICIT NONE REAL NUMER, DENOM LOGICAL SAFDIV REAL BIG ! local variable LOGICAL FIRST REAL RMACHINE ! external function SAVE BIG, FIRST DATA FIRST / .TRUE. / * Begin IF (FIRST) THEN FIRST = .FALSE. BIG = RMACHINE(3) END IF SAFDIV = ABS(DENOM) > MAX(0.0, ABS(NUMER / BIG)) RETURN END ! of SAFDIV *----------------------------------------------------------------------- * Safe to divide in double precision *----------------------------------------------------------------------- FUNCTION DSAFDIV (NUMER, DENOM) IMPLICIT NONE DOUBLE PRECISION NUMER, DENOM LOGICAL DSAFDIV DOUBLE PRECISION BIG ! local variable LOGICAL FIRST DOUBLE PRECISION DMACHINE ! external function SAVE BIG, FIRST DATA FIRST / .TRUE. / * Begin IF (FIRST) THEN FIRST = .FALSE. BIG = DMACHINE(3) END IF DSAFDIV = ABS(DENOM) > MAX(0.0D0, ABS(NUMER / BIG)) RETURN END ! of DSAFDIV *----------------------------------------------------------------------- * Almost equality test *----------------------------------------------------------------------- FUNCTION ALMEQ (R, S) IMPLICIT NONE REAL R, S LOGICAL ALMEQ REAL AR, AS, TWOEPS, LIT LOGICAL FIRST REAL RMACHINE ! external function SAVE TWOEPS, LIT, FIRST DATA FIRST / .TRUE. / IF (FIRST) THEN TWOEPS = 2. * RMACHINE(1) LIT = RMACHINE(2) FIRST = .FALSE. END IF AR = ABS(R) AS = ABS(S) IF (AR < LIT .AND. AS < LIT) THEN ! all denormal numbers are equal ALMEQ = .TRUE. ELSE ALMEQ = ABS(R - S) < TWOEPS * MAX(AR, AS) END IF RETURN END ! of ALMEQ *----------------------------------------------------------------------- * same thing in double precision *----------------------------------------------------------------------- FUNCTION DALMEQ (R, S) IMPLICIT NONE DOUBLE PRECISION R, S LOGICAL DALMEQ DOUBLE PRECISION AR, AS, TWOEPS, LIT LOGICAL FIRST DOUBLE PRECISION DMACHINE SAVE TWOEPS, LIT, FIRST DATA FIRST / .TRUE. / IF (FIRST) THEN TWOEPS = DMACHINE(1) LIT = DMACHINE(2) FIRST = .FALSE. END IF AR = ABS(R) AS = ABS(S) IF (AR < LIT .AND. AS < LIT) THEN DALMEQ = .TRUE. ELSE DALMEQ = ABS(R - S) < TWOEPS * MAX(AR, AS) END IF RETURN END ! of DALMEQ *----------------------------------------------------------------------- * Random permutation of 1...N * origin: Durstenfeld, 1964 *----------------------------------------------------------------------- SUBROUTINE SHUFFL (I, N) IMPLICIT NONE INTEGER I, N DIMENSION I(N) INTEGER J, K, L INTEGER HAT DO J = 1, N I(J) = J END DO DO J = N, 2, -1 K = HAT(J) L = I(K) I(K) = I(J) I(J) = L END DO RETURN END ! of SHUFFL *----------------------------------------------------------------------- * replacement for CEILING() *----------------------------------------------------------------------- FUNCTION CEIL (X) IMPLICIT NONE REAL X INTEGER CEIL REAL W ! whole parts LOGICAL ALMEQ ! almost-equal function * no support for negative numbers IF (X < 0.) THEN CEIL = 0 RETURN END IF W = AINT(X) IF (ALMEQ(X,W)) THEN CEIL = NINT(W) ELSE CEIL = NINT(W) + 1 END IF RETURN END ! of CEIL *??????????????????????????????????????????????????????????????????????? * succod - succint code: translate arbitrary integers to small integers * *_Name______Type_______In/Out____Description____________________________ * JNAM(L) Integer In Data array * TRAN(M) Integer Out Translation * IND(L) Integer Neither Workspace * L Integer In Length of JNAM * M Integer In Length of TRAN * N Integer Out Distinct entries in JNAM *??????????????????????????????????????????????????????????????????????? SUBROUTINE SUCCOD (JNAM, TRAN, IND, L, m, n) IMPLICIT NONE INTEGER JNAM, TRAN, IND, L, M, N DIMENSION JNAM(L), TRAN(M), IND(L) * local variables INTEGER I, J, K * Begin. Sort implicitly. CALL QSORTI (JNAM, IND, L) * Compare adjacent entries K = IND(1) N = 1 TRAN(1) = JNAM(K) DO I = 2, L J = IND(I) IF (JNAM(J) .NE. JNAM(K)) THEN K = J N = N + 1 TRAN(N) = JNAM(J) END IF END DO RETURN END ! of SUCCOD *??????????????????????????????????????????????????????????????????????? * nuniq - number of distinct integers in an array * *_Name______Type_______In/Out____Description____________________________ * JNAM(N) Integer In Data array * L Integer In length of JNAM * N Integer Out Distinct entries in JNAM *??????????????????????????????????????????????????????????????????????? SUBROUTINE NUNIQ (JNAM, L, N) IMPLICIT NONE INTEGER JNAM, L, N DIMENSION JNAM(L) * Local variables INTEGER I, J * Begin. Sort. CALL ISORT (JNAM, L) * compare adjacent entries N = 1 J = 1 DO I = 2, L IF (JNAM(I) .NE. JNAM(J)) THEN J = I N = N + 1 END IF END DO RETURN END ! of NUNIQ *----------------------------------------------------------------------- * Determinant of a 4 x 4 matrix, by expansion *----------------------------------------------------------------------- FUNCTION DET4(X) IMPLICIT NONE DOUBLE PRECISION X(4,4), DET4 DOUBLE PRECISION A, B, C, D, E, F, G, H, I, J A = X(3,3) * X(4,4) - X(3,4) * X(4,3) B = X(3,2) * X(4,4) - X(3,4) * X(4,2) C = X(3,2) * X(4,3) - X(3,3) * X(4,2) D = X(3,1) * X(4,4) - X(3,4) * X(4,1) E = X(3,1) * X(4,3) - X(3,3) * X(4,1) F = X(3,1) * X(4,2) - X(3,2) * X(4,1) G = X(1,1) * (X(2,2) * A - X(2,3) * B + X(2,4) * C) H = X(1,2) * (X(2,1) * A - X(2,3) * D + X(2,4) * E) I = X(1,3) * (X(2,1) * B - X(2,2) * D + X(2,4) * F) J = X(1,4) * (X(2,1) * C - X(2,2) * E + X(2,3) * F) DET4 = G - H + I - J RETURN END ! of DET4 *----------------------------------------------------------------------- * cfit - get a least-squares polynomial of order 3. * *___Name_____Type______In/Out___Description_____________________________ * W(N) Real In Importance weights * X(N) Real In Independent variable * Y(N) Real In Dependent variable * N Integer In # points in approximation * C(4) Real Out Coefficients * IFAULT Integer Out Error code *----------------------------------------------------------------------- SUBROUTINE CFIT (W, X, Y, N, C, IFAULT) IMPLICIT NONE INTEGER N, IFAULT REAL W(N), X(N), Y(N), C(4) * Local variables INTEGER I REAL EPS, XO, XF, DZDX ! rescale X DOUBLE PRECISION DEPS, ! very small number $ W1, Y1, X1, X2, X3, ! temporary storage $ SW, SX, SX2, SX3, SX4, SX5, SX6, ! sums in normal equations $ SY, SXY, SX2Y, SX3Y, $ SUB1, SUB2, SUB3, SUB4, SUB5, SUB6, DET, $ U(4,4), D(4), X0(4), Y0(4) * External functions REAL RMACHINE DOUBLE PRECISION DMACHINE, DET4 *----------------------------------------------------------------------- * Begin *----------------------------------------------------------------------- IFAULT = 0 EPS = RMACHINE(1) DEPS = 4 * DMACHINE(1) * Scrutinize input IF (N < 1) THEN PRINT *, 'cfit: no data' RETURN END IF *----------------------------------------------------------------------- * Rescale z = -1 + 2 @ (x - xo) / (xf - xo) * Thus dz/dx = 2 / (xf - xo) *----------------------------------------------------------------------- * Find limits of x XO = X(1) XF = X(1) DO I = 2, N XO = MIN(XO, X(I)) XF = MAX(XF, X(I)) END DO IF (XF - XO < EPS) THEN ! X all same, return average SY = 0. DO I = 1, N SY = SY + Y(I) END DO C(1) = SY / DBLE(N) C(2) = 0. C(3) = 0. C(4) = 0. RETURN END IF DZDX = 2. / (XF - XO) ! dz/dx *----------------------------------------------------------------------- * Normal equations for a cubic polynomial. Refer to * Murray Spiegel, Schaum's Outline of Statistics, 1988 [1961] *----------------------------------------------------------------------- SW = 0. SX = 0. SX2 = 0. SX3 = 0. SX4 = 0. SX5 = 0. SX6 = 0. SY = 0. SXY = 0. SX2Y = 0. SX3Y = 0. DO I = 1, N W1 = W(I) X1 = (X(I) - XO) * DZDX - 1. Y1 = Y(I) X2 = X1 * X1 X3 = X2 * X1 SW = SW + W1 ! SUM w SX = SX + X1 * W1 ! SUM x SX2 = SX2 + X2 * W1 ! SUM x^2 SX3 = SX3 + X3 * W1 ! SUM x^3 SX4 = SX4 + X2 * X2 * W1 ! SUM x^4 SX5 = SX5 + X3 * X2 * W1 ! SUM x^5 SX6 = SX6 + X3 * X3 * W1 ! SUM x^6 SY = SY + Y1 * W1 ! SUM y SXY = SXY + X1 * Y1 * W1 ! SUM x @ y SX2Y = SX2Y + X2 * Y1 * W1 ! SUM x^2 @ y SX3Y = SX3Y + X3 * Y1 * W1 ! SUM x^3 @ y END DO * Solve the system IF (N - 3) 20, 30, 40 *----------------------------------------------------------------------- 20 CONTINUE ! least-squares line *----------------------------------------------------------------------- DET = SW * SX2 - SX * SX IF (DET < DEPS) GO TO 50 D(1) = (SY * SX2 - SX * SXY) / DET D(2) = (SW * SXY - SX * SY) / DET C(4) = 0. C(3) = 0. C(2) = SNGL(D(2) * DZDX) C(1) = SNGL(D(1) - D(2) - C(2) * XO) RETURN *----------------------------------------------------------------------- 30 CONTINUE ! parabola *----------------------------------------------------------------------- SUB1 = SX2 * SX4 - SX3 * SX3 SUB2 = SX * SX4 - SX2 * SX3 SUB3 = SX * SX3 - SX2 * SX2 SUB4 = SXY * SX4 - SX2Y * SX3 SUB5 = SXY * SX3 - SX2 * SX2Y SUB6 = SX * SX2Y - SX2 * SXY DET = SW*SUB1 - SX*SUB2 + SX2*SUB3 IF (DET < DEPS) GO TO 20 D(1) = ( SY*SUB1 - SX*SUB4 + SX2*SUB5 ) / DET D(2) = ( SW*SUB4 - SY*SUB2 + SX2*SUB6 ) / DET D(3) = ( SW*(-SUB5) - SX*SUB6 + SY*SUB3 ) / DET * Interpolate the approximated polynomial X0(1) = XO X0(2) = (XO + XF) / 2. X0(3) = XF Y0(1) = D(1) - D(2) + D(3) Y0(2) = D(1) Y0(3) = D(1) + D(2) + D(3) CALL POLY(X0, Y0, D, 2) * Reduce precision DO I = 1, 3 C(I) = SNGL(D(I)) END DO C(4) = 0. RETURN *----------------------------------------------------------------------- 40 CONTINUE ! cubic *----------------------------------------------------------------------- U(1,1) = SW U(2,1) = SX U(3,1) = SX2 U(4,1) = SX3 U(1,2) = SX U(2,2) = SX2 U(3,2) = SX3 U(4,2) = SX4 U(1,3) = SX2 U(2,3) = SX3 U(3,3) = SX4 U(4,3) = SX5 U(1,4) = SX3 U(2,4) = SX4 U(3,4) = SX5 U(4,4) = SX6 DET = DET4(U) IF (DET < DEPS) GO TO 30 U(1,1) = SY U(2,1) = SXY U(3,1) = SX2Y U(4,1) = SX3Y D(1) = DET4(U) / DET U(1,1) = SW U(2,1) = SX U(3,1) = SX2 U(4,1) = SX3 U(1,2) = SY U(2,2) = SXY U(3,2) = SX2Y U(4,2) = SX3Y D(2) = DET4(U) / DET U(1,2) = SX U(2,2) = SX2 U(3,2) = SX3 U(4,2) = SX4 U(1,3) = SY U(2,3) = SXY U(3,3) = SX2Y U(4,3) = SX3Y D(3) = DET4(U) / DET U(1,3) = SX2 U(2,3) = SX3 U(3,3) = SX4 U(4,3) = SX5 U(1,4) = SY U(2,4) = SXY U(3,4) = SX2Y U(4,4) = SX3Y D(4) = DET4(U) / DET * Interpolate the approximated polynomial X0(1) = XO X0(2) = XO + (XF - XO) / 3. X0(3) = XO + 2. * (XF - XO) / 3. X0(4) = XF Y0(1) = D(1) - D(2) + D(3) - D(4) X1 = -1.D0 / 3.D0 Y0(2) = ((D(4) * X1 + D(3)) * X1 + D(2)) * X1 + D(1) X1 = 1.D0 / 3.D0 Y0(3) = ((D(4) * X1 + D(3)) * X1 + D(2)) * X1 + D(1) Y0(4) = D(1) + D(2) + D(3) + D(4) CALL POLY (X0, Y0, D, 3) * Reduce precision DO I = 1, 4 C(I) = SNGL(D(I)) END DO RETURN ! success *----------------------------------------------------------------------- 50 CONTINUE ! last chance. Endpoint interpolation *----------------------------------------------------------------------- PRINT *, 'cfit: using endpoints' C(3) = 0. C(4) = 0. C(2) = (Y(N) - Y(1)) / (XF - XO) C(1) = -C(2) * XO + Y(1) RETURN END ! of cfit *----------------------------------------------------------------------- * wlsl - draw a weighted least-squares line through a data sequence *----------------------------------------------------------------------- * I/O LIST *___Name_____Type______In/Out___Description_____________________________ * W(N) Real In Importance weights of data points * X(N) Real In Data independent variable * Y(N) Real In Data dependent variable * N Integer In Number of data points * CEPT Real Out Y-intercept * SLOPE Real Out Slope * IFAULT Integer Out Error code 0=no errors *----------------------------------------------------------------------- SUBROUTINE WLSL (W, X, Y, N, CEPT, SLOPE, IFAULT) IMPLICIT NONE INTEGER N, IFAULT REAL W, X, Y, SLOPE, CEPT DIMENSION W(N), X(N), Y(N) INTEGER J DOUBLE PRECISION $ DET, ! determinant of normal equations $ DEPS, DMACHINE, ! machine double precision $ WNEW, XNEW, YNEW, ! working values of w, x, y $ WX, SW, SWX, SWY, SWX2, SWXY ! least squares sums LOGICAL FIRST SAVE FIRST, DEPS DATA FIRST / .TRUE. / *----------------------------------------------------------------------- * Begin *----------------------------------------------------------------------- IFAULT = 0 IF (FIRST) THEN DEPS = DMACHINE(1) FIRST = .FALSE. END IF IF (N < 2) THEN ! not enough data IFAULT = 3 RETURN END IF * Weighted least squares sums SW = 0 SWX = 0 SWY = 0 SWX2 = 0 SWXY = 0 DO J = 1, N WNEW = DBLE(W(J)) XNEW = DBLE(X(J)) YNEW = DBLE(Y(J)) SW = SW + WNEW WX = WNEW * XNEW SWX = SWX + WX SWY = SWY + WNEW * YNEW SWX2 = SWX2 + WX * XNEW SWXY = SWXY + WX * YNEW END DO * Solve normal equations DET = SWX2 * SW - SWX * SWX IF (ABS(DET) > DEPS) THEN CEPT = SNGL((SWX2 * SWY - SWX * SWXY) / DET) SLOPE = SNGL((SW * SWXY - SWX * SWY) / DET) ELSE IFAULT = 1 PRINT *, 'wlsl: warning, zero determinant' SLOPE = (Y(N) - Y(1)) / (X(N) - X(1)) CEPT = Y(1) - SLOPE * X(1) END IF RETURN END ! of wlsl *----------------------------------------------------------------------- * drdi - Double Precision Rank, Determinant, and Inverse * * Intended for use inverting covariance matrices in Linear * Discriminant Analysis. These matrices are symmetric. The eigenvector * factorization is avoided because subroutines are known to fail in the * case of zero vectors, such as when there is a variable that doesn't vary, * which is expected to occur. It is required to find both U and V, * because the left and right singular vectors may have different signs. *----------------------------------------------------------------------- * TRANSFER LIST *___Name______Type_______________In/Out___Description___________________ * X(N,N) Double Precision Both Square matrix, returned inverted * N Integer In Size of matrix X * RANK Integer Out Pseudo-rank of X * LAD Double Precision Out Logarithm of the absolute value * of the determinant of X * U(N,N) Double Precision Work Left orthogonal factors of X * S(N) Double Precision Work Singular values of X * V(N,N) Double Precision Work Right orthogonal factors of X * E(N) Double Precision Work Used by DSVDC * WORK(N) Double Precision Work Workspace for DSVDC * IERR Integer Out Error code *----------------------------------------------------------------------- SUBROUTINE DRDI (X, N, RANK, LAD, U, S, V, E, WORK, IERR) IMPLICIT NONE INTEGER N, RANK, IERR DOUBLE PRECISION X, LAD, U, S, V, E, WORK DIMENSION X(N,N), U(N,N), S(N), V(N,N), E(N), WORK(N) * Local variables INTEGER I, J, K, JOB PARAMETER (JOB = 11) DOUBLE PRECISION $ PROD, ! dot product U . V $ TOL, ! a tolerance $ TOT ! a total * external function DOUBLE PRECISION DMACHINE *----------------------------------------------------------------------- * Begin. *----------------------------------------------------------------------- TOL = 0.D0 ! max of column sums of absolute values DO J = 1, N TOT = 0.D0 DO I = 1, N TOT = TOT + ABS(X(I,J)) END DO IF (TOT > TOL) TOL = TOT END DO TOL = TOL * DMACHINE(1) * Factor X = U @ S @ V^T CALL DSVDC (X, N, N, N, S, E, U, N, V, N, WORK, JOB, IERR) IF (IERR .NE. 0) THEN PRINT *, 'dsvdc: error #', IERR RETURN END IF !!!!!!!!!! debug: DO J = 1, N DO I = 1, N IF (U(I,J) .NE. U(I,J) .OR. V(I,J) .NE. V(I,J)) IERR = 1 END DO END DO 10 FORMAT (E12.4, $) 20 FORMAT () IF (IERR .EQ. 1) THEN PRINT *, 'X: ' DO J = 1, N DO I = 1, N PRINT 10, X(I,J) END DO PRINT 20 END DO PAUSE END IF !!!!!!!!!!!!!!!!!!!!! !! TOL = N * DMACHINE(1) ! n @ epsilon RANK = 0 LAD = 0.D0 DO K = 1, N PROD = 0. DO I = 1, N ! detect negative eigenvalues PROD = PROD + U(I,K) * V(I,K) END DO !!!!!!!!!!!!!!!! IF (PROD < 0.D0) PRINT *, 'drdi: dot-product is ', PROD !!!!!!!!!!!!!!!!!!!!!!! IF (S(K) > TOL .AND. PROD > 0.D0) THEN RANK = RANK + 1 ! count rank LAD = LAD + LOG(S(K)) ! multiply non-null entries of S S(K) = 1.D0 / S(K) ! invert S ELSE ! zero out small values S(K) = 0.D0 ! William Press: "It's not often you get to set 0=Inf." END IF END DO IF (RANK .EQ. 0) THEN IERR = -1 RETURN END IF * Multiply X^(-1) = U @ S^(-1) @ V^T DO J = 1, N ! clear DO I = 1, N X(I,J) = 0.D0 END DO END DO DO 50 K = 1, N IF (S(K) .EQ. 0.D0) GO TO 50 ! if comparison fails, no harm done DO J = 1, N DO I = 1, N X(I,J) = X(I,J) + U(I,K) * S(K) * V(J,K) END DO END DO 50 CONTINUE * Force symmetry DO J = 1, N DO I = J+1, N X(I,J) = 0.5D0 * (X(I,J) + X(J,I)) X(J,I) = X(I,J) END DO END DO * Force non-negative diagonal DO K = 1, N IF (X(K,K) < 0.D0) X(K,K) = 0.D0 END DO RETURN END ! of drdi