C   BRENT.FOR                                                           
                                                                        
      SUBROUTINE BRENT(A,B,X,MH0I,MSEDI,DH0I,DSEDI,ADA)
C                                                                       
C MICHAEL B. PORTER  8/84                                               
C                                                                       
C FORTRAN CONVERSION OF ALGOL PROGRAM PUBLISHED IN                      
C    THE COMPUTER JOURNAL 14(4):422-425 (1971)                          
C    BY R. P. BRENT                                                     
C                                                                       
C RETURNS A ZERO X OF THE FUNCTION F IN THE GIVEN INTERVAL              
C    [A,B], TO WITHIN A TOLERANCE 6*MACHEP*ABS(X)+2*T, WHERE            
C    MACHEP IS THE RELATIVE MACHINE PRECISION AND T IS A POSITIVE       
C    TOLERANCE.  THE PROCEDURE ASSUMES THAT FUNCT(A) AND FUNCT(B) HAVE  
C    DIFFERENT SIGNS.                                                   
C                                                                       
                                                                        
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)                              
                                                                        
                                                                        
      DOUBLE PRECISION DH0I, DSEDI                                      
      DOUBLE PRECISION CON1, CON2, CON3, CON4, CON5, SEDK               
      DOUBLE PRECISION MACHEP, M, T, TREF                               
      DOUBLE PRECISION ADA( * )
                                                                        
      COMMON /LUNIT/ LUPLP, LUPLT, LUPRT
      COMMON /CONST/ CON1, CON2, CON3, CON4, CON5, SEDK                 
                                                                        
      DATA TREF/1.0D-13/                                                
                                                                        
      MACHEP = 1.0E-16                                                  
      TEN = 10.0                                                        
                                                                        
      T = TREF * MAX( ABS( A ), ABS( B ) )                              
                                                                        
C      CALL FUNCT(A,FA,IEXPA)                                           
C      CALL FUNCT(B,FB,IEXPB)                                           
      CALL CHARAC(A,MH0I,MSEDI,DH0I,DSEDI,FA,NOVFLA,ADA)
      CALL CHARAC(B,MH0I,MSEDI,DH0I,DSEDI,FB,NOVFLB,ADA)
      IEXPA=NOVFLA*20                                                   
      IEXPB=NOVFLB*20                                                   
      IF ( ( (FA .GT. 0.0) .AND. (FB .GT. 0.0) ) .OR.                   
     &     ( (FA .LT. 0.0) .AND. (FB .LT. 0.0) ) ) THEN
         WRITE(LUPRT,*)
     &   ' *** ZBRENT ERROR: FUNCT SGN SAME AT INTRVL ENDPTS'
         STOP                                                           
C         RETURN                                                        
      ENDIF                                                             
C                                                                       
C     INTERNAL ROOT                                                     
C                                                                       
 2000 C = A                                                             
      FC = FA                                                           
      IEXPC = IEXPA                                                     
      E = B-A                                                           
      D = E                                                             
C                                                                       
C     EXTERNAL ROOT                                                     
C                                                                       
      IF (IEXPA .LT. IEXPB) THEN                                        
         F1 = FC*TEN**(IEXPC-IEXPB)                                     
         F2 = FB                                                        
      ELSE                                                              
         F1 = FC                                                        
         F2 = FB*TEN**(IEXPB-IEXPC)                                     
      ENDIF                                                             
 3000 IF (ABS(F1) .LT. ABS(F2)) THEN                                    
         A = B                                                          
         B = C                                                          
         C = A                                                          
         FA = FB                                                        
         IEXPA = IEXPB                                                  
         FB = FC                                                        
         IEXPB = IEXPC                                                  
         FC = FA                                                        
         IEXPC = IEXPA                                                  
      ENDIF                                                             
      TOL = 2.0*MACHEP*ABS(B)+T                                         
      M = 0.5*(C-B)                                                     
      IF ((ABS(M) .GT. TOL) .AND. (FB .NE. 0.0)) THEN                   
C     ------ SEE IF A BISECTION IS FORCED                               
      IF (IEXPA .LT. IEXPB) THEN                                        
         F1 = FA*TEN**(IEXPA-IEXPB)                                     
         F2 = FB                                                        
      ELSE                                                              
         F1 = FA                                                        
         F2 = FB*TEN**(IEXPB-IEXPA)                                     
      ENDIF                                                             
         IF ( (ABS(E) .LT. TOL) .OR.                                    
     &       (ABS( F1 ) .LE. ABS(F2)) ) THEN                            
            E = M                                                       
            D = E                                                       
         ELSE                                                           
            S = FB/FA*TEN**(IEXPB-IEXPA)                                
            IF (A .EQ. C) THEN                                          
C              ------ LINEAR INTERPOLATION                              
               P = 2.0*M*S                                              
               Q = 1.0-S                                                
            ELSE                                                        
C              ------ INVERSE QUADRATIC INTERPOLATION                   
               Q = FA/FC*TEN**(IEXPA-IEXPC)                             
               R = FB/FC*TEN**(IEXPB-IEXPC)                             
               P = S*(2.0*M*Q*(Q-R)-(B-A)*(R-1.0))                      
               Q = (Q-1.0)*(R-1.0)*(S-1.0)                              
            ENDIF                                                       
            IF (P .GT. 0.0) THEN                                        
               Q = -Q                                                   
            ELSE                                                        
               P = -P                                                   
            ENDIF                                                       
            S = E                                                       
            E = D                                                       
            IF ((2.0*P .LT. 3.0*M*Q-ABS(TOL*Q)) .AND.                   
     &          (P .LT. ABS(0.5*S*Q))) THEN                             
               D = P/Q                                                  
            ELSE                                                        
               E = M                                                    
               D = E                                                    
            ENDIF                                                       
         ENDIF                                                          
         A = B                                                          
         FA = FB                                                        
         IEXPA = IEXPB                                                  
         IF (ABS(D) .GT. TOL) THEN                                      
            B = B+D                                                     
         ELSE                                                           
            IF (M .GT. 0.0) THEN                                        
               B = B+TOL                                                
            ELSE                                                        
               B = B-TOL                                                
            ENDIF                                                       
         ENDIF                                                          
      CALL CHARAC(B,MH0I,MSEDI,DH0I,DSEDI,FB,NOVFLB,ADA)
      IEXPB=NOVFLB*20                                                   
C         CALL FUNCT(B,FB,IEXPB)                                        
         IF ((FB .GT. 0.0) .EQV. (FC .GT. 0.0)) GOTO 2000               
         GOTO 3000                                                      
      ENDIF                                                             
      X = B                                                             
      RETURN                                                            
      END