C     PROGRAMA PRINCIPAL P/ CHAMADA DO ALGORITMO 
C               PENALIZACAO HIPERBOLICA
      PROGRAM PENHIPERB
      IMPLICIT REAL*8 (A-H,O-Z)                                                                              
      COMMON /VA13B/ IPRINT                                 
      DIMENSION XINIC(2) 
      
      character * 12 arqsaida        

      EXTERNAL BFGS,FO,GFO,FG,GFG                                           

C     X = PONTO INICIAL                                                     
      DATA XINIC /10.0, 10.0/    
C     N = NUMERO DE VARIAVEIS NO PROBLEMA                               
      DATA N /2/
C     M = NUMERO DE RESTRICOES                                          
      DATA M /3/
      U0=10.D0 
      TAL0=10.D0 

C  NAO VAI HAVER IMPRESSAO NA  BFGS !
      IPRINT=0
      LP = 4      
      arqsaida = 'hs116.out'   
                                                                                
C     CHAMADA DO ALGORITMO 'LAGRANGEANO HIPERBOLICO'                          
      CALL PHA1( N,M,XINIC,U0,TAL0,ARQSAIDA)                                                                                                                                                       
      
      STOP                                                              
                                                                        
      END PROGRAM                                                             
      
      SUBROUTINE PHA1(N,M,XINIC,U0,TAL0,ARQSAIDA)                     
C     SUBROTINA QUE EFETUA A MINIMINIZACAO SEGUNDO O METODO DO                  
C                      LAGRANGEANO HIPERBOLICO                                  
C   ENTRADAS:
C        INTEIRO : N - DIMENSOES DAS VARIAVEIS DE DECISAO
C        INTEIRO : M - QUANTIDADE DE RESTRICOES
C        REAL(N) : X - PONTO INICIAL
C		 CAHR(12): ARQSAIDA - NOME DO ARQUIVO DE SAIDA
                                                                             
																				
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
      COMMON /PENHIP/ U(200) , TAL(200)                                           
      DIMENSION XINIC(N)
      DIMENSION INVIAB(M)
      DIMENSION X(N), XOLD(N)                                   
      DIMENSION G(N)                             
      DIMENSION ALFAG(200)                                       
      DIMENSION GRESTR(200) ,  folga(200)                          
      DIMENSION XANT(100) , UANT(200) , folgaant(200)                                  
      DIMENSION Ucalc(200)
      DIMENSION W(6000) , SCALE(100)                                           
      LOGICAL INVIAVEL, CONVERGIU
      CHARACTER * 12 ARQSAIDA        


      EXTERNAL FUNMODA1 
c UNIDADE LOGICA DE SAIDA      
      DATA NULSAIDA /8/  
      DATA NULSAIDA2 /7/
      DATA PI/ 3.141592653589793 D0/                                            
                                                                                
      DATA   SCALE/ 100 * 1.D0 /                                                
                                                                                
      DATA EPS5/ 1.D-5/                                                        
     5     EPS6/ 0.1D-7/                                                       

C     NLHMAX  -  NUM. MAXIMO ITERACOES FASE 2 DA PEN. HIPERBOLICA               
      DATA NLHMAX/ 25  /                                                         
                                              
C  INICIALIZACOES
      X = XINIC
      XOLD = XINIC
      NLH = 0                                                                   
      REDUC_TAL = 10.D0   
      AUMENT_U = DSQRT(10.D0)
      CONVERGIU = .FALSE.
      
      OPEN(NULSAIDA,FILE='REL_CONV.TXT')
      OPEN(NULSAIDA2,FILE='REL_PARM.TXT')
      
      WRITE(NULSAIDA,2110) 
      WRITE(NULSAIDA,2112) 
      WRITE(NULSAIDA,2111) 
      WRITE(NULSAIDA,2112) 
      WRITE(*,2112) 
      WRITE(*,2111) 
      WRITE(*,2112) 
      WRITE(NULSAIDA2,660) 
      WRITE(NULSAIDA2,661) 
      WRITE(NULSAIDA2,662) 
      WRITE(NULSAIDA2,663) 
      WRITE(NULSAIDA2,664) 
    
      DO ND=1,200
        U(ND) = U0
        TAL(ND) = TAL0
      ENDDO

c PONTOS INICIAIS 
      CALL FO(X,FPURO)                                                        
      CALL FG(X,GRESTR) 
      CALL FOM ( X , N , M, FVALOR )	  

C VERIFICA SE O PONTO INICIAL E VIAVEL
        INVIAVEL = .FALSE.
         DO ND=1,M    
          INVIAB(ND) = 0
          IF ( GRESTR(ND).lt. -EPS5  ) THEN
            INVIAVEL = .TRUE.
            INVIAB(ND) = 1
          ENDIF
        ENDDO
   
      IF (INVIAVEL) THEN
        WRITE(*,2113) NLH,X(1),X(2),FPURO,FVALOR,'  TRUE  '  
        WRITE(NULSAIDA,2113) NLH,X(1),X(2),FPURO,FVALOR,'  TRUE  '  
      ELSE
        WRITE(*,2113) NLH,X(1),X(2),FPURO,FVALOR,'  FALSE '  
        WRITE(NULSAIDA,2113) NLH,X(1),X(2),FPURO,FVALOR,'  FALSE '  
      ENDIF
      WRITE(NULSAIDA2,665) NLH,(INVIAB(ND),TAL(ND),U(ND),ND=1,M)
      
                                                                            
C     FASE DO LAGRANGEANO HIPERBOLICO                                              
       DO WHILE (.NOT. CONVERGIU)
C NUM. ITER. DA PH MAIOR QUE MAXIMO
      IF (NLH.GT.NLHMAX) THEN                                         
         WRITE(NULSAIDA,1401)NLHMAX                                     
         WRITE(*,1401)NLHMAX                                     
1401     FORMAT(///,' NUMERO DE ITERACOES DO LAGRANG. HIPERB. ' ,         
     1      'ATINGIU O MAXIMO ESPECIFICADO = ',I3)                                    
         GO TO 9990                                                                
      ENDIF
C CONTADOR DE ITERACOES
      NLH = NLH + 1                                                             
C---- MINIMIZA O LAGRANGEANO HIPERBOLICO                                        
      ACC = TAL(1) / 10.D0**03                                                     
      ACC = DMIN1 ( ACC , 0.1 D-3 ) 
      CALL BFGS(N,X,F,G,SCALE,ACC,W)                           

      CALL FO(X,FPURO)                                                        
      CALL FG(X,GRESTR) 
      CALL FOM ( X , N , M, FVALOR )	  

C VERIFICA SE O PONTO ENCONTRATO E VIAVEL
        INVIAVEL = .FALSE.
         DO ND=1,M    
          INVIAB(ND) = 0
          IF ( GRESTR(ND).lt. -EPS5  ) THEN
            INVIAVEL = .TRUE.
            INVIAB(ND) = 1
          ENDIF
        ENDDO

         DO ND=1,M  
            IF (INVIAB(ND) == 0) THEN
c        reducao da penalidade (segunda fase do algoritmo penal. hiperb.)
              U(ND)=U(ND)
              TAL(ND)=TAL(ND)/REDUC_TAL
            ELSE
c         atualizacao "pen. hiperb." - aumenta angulo (primeira fase Pen. Hiperb.)
              U(ND)= AUMENT_U * U(ND)    
              TAL(ND)=TAL(ND)
            ENDIF
          ENDDO
 

      
        IF (INVIAVEL) THEN
          WRITE(*,2113) NLH,X(1),X(2),FPURO,FVALOR,'  TRUE  '  
          WRITE(NULSAIDA,2113) NLH,X(1),X(2),FPURO,FVALOR,'  TRUE  '  
        ELSE
          WRITE(*,2113) NLH,X(1),X(2),FPURO,FVALOR,'  FALSE '  
          WRITE(NULSAIDA,2113) NLH,X(1),X(2),FPURO,FVALOR,'  FALSE '  
        ENDIF
        WRITE(NULSAIDA2,665) NLH,(INVIAB(ND),TAL(ND),U(ND),ND=1,M)
                                          
C CRITERIO DE PARADA
        CONVERGIU=((SUM(DABS(X - XOLD)) .LT. EPS5).AND.(.NOT.INVIAVEL))
        XOLD = X
      ENDDO
                            
C     IMPRESSAO FINAL                                                           
9990  WRITE(NULSAIDA,2112)                                                            
      WRITE(*,2112)                                                            
      WRITE(NULSAIDA2,664)                                                            
      WRITE(*,'(A)')' FIM DO PROCESSAMENTO'
      READ(*,*)
      STOP

2110   FORMAT(///,'      RELATORIO DE CONVERGENDIA DO ALGORITMO DE PE'
     +  'ENALIDADE HIPERBOLICA  ',///)
2111   FORMAT(' IT       X1          X2           F           FOM    '
     +        '   INVIAB   ' )
2112   FORMAT('X---X------------X------------X------------X----------'
     +        '--X--------X ' )
2113   FORMAT(1X,I3,1X,F12.7,1X,F12.7,1X,F12.7,1X,F12.7,1X,A8)                

660    FORMAT(///,'      RELATORIO DO COMPORTAMENTO DOS PARAMETROS DE P'
     +  'ENALIDADE POR RESTRICAO   ',/// )
661    FORMAT('X---X-----------------------------X--------------------'
     +   '---------X-----------------------------X' )
662    FORMAT('           G1 = (X1 >= 0)                G2 = (X2 >= 0)'
     +   '             G3 = (X1^2 + X2^2 >= 1) ')
663    FORMAT(' IT  VIAV     TAU        LAMBDA    VIAV     TAU        '
     +   'LAMBDA    VIAV     TAU        LAMBDA    ' )
664    FORMAT('X---X---X------------X------------X---X------------X---'
     +   '---------X---X------------X------------X' )
665    FORMAT(1X,I3,6(1X,I3,1X,F12.7,1X,F12.7) )

      END                                                                       
