*-----------------------------------------------------------------------
* 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
