      MODULE subgdc
      
        !*..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..*
        !| .**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**. |
        !| |                                                                                  | |
        !| |      SUBGDC - THE AGGREGATE SUBGRADIENT METHOD FOR NONSMOOTH DC OPTIMIZATION     | | 
        !| |                                 (version 1)                                      | |
        !| |                                                                                  | |
        !| |                                                                                  | |
        !| |                                                                                  | |
        !| |                                                                                  | |
        !| .**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**. |
        !*..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..*     

 
        USE constants, ONLY : dp    ! Double precision (i.e. accuracy)
        USE functions               ! Contains INFORMATION from the USER
        
        IMPLICIT NONE    
                
        CONTAINS
                
        !*..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..*
        !| .**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**. |
        !| |                                                                                  | |
        !| |                            CONTAINS SUBROUTINES:                                 | | 
        !| |                                                                                  | |
        !| |    SUBGDC_algorithm   : The aggregate subgradient method for DC optimization.    | | 
        !| |                                                                                  | |       
        !| |                                                                                  | |
        !| .**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**. |
        !*..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..*   
        
    
    
        
        
        ! _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _  
        !| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | |
        ! START ** START ** START ** START ** START ** START ** START ** START ** START ** START  
        !|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|        
        ! <>  <>  <>  <>  <>  <>  <>  <>  <>  <>  <>  <>  <>  <>  <>  <>  <>  <>  <>  <>  <>  <>
        !***************************************************************************************
        !  ----------------------------------------------------------------------------------  |
        !  |                                                                                |  |
        !  |         THE AGGREGATE SUBGRDIENT METHOD FOR NONSMOOTH DC OPTIMIZATION          |  |
        !  |                                                                                |  |
        !  ----------------------------------------------------------------------------------  |
        !***************************************************************************************
        
           SUBROUTINE SUBGDC_algorithm(x_0, x_solution, f_solution, f0, mit, mrounds, termination, counter, &
                            &  CPUtime, time_limit, iprint,  consecutive) 
            ! 
            ! Solves the unconstrained nonsmooth DC minimization problem
            !
            ! INPUT:  * 'x_0'            : A starting point
            !         * 'mit'            : The maximum number of inner iterations      
            !         * 'mrouds'         : The maximum number of rounds during one 'inner iteration'
            !         * 'iprint'         : Specifies the print          
            !         * 'time_limit '    : the CPU time limit  
            !         * 'consecutive'    : .TRUE. if algorithm can be terminated if the difference between consecutive function values is small                     
            !
            ! OUTPUT: * 'x_solution' : The solution obtained to the minimizatin problem
            !         * 'f_solution' : The objective function value at the solution 'x_solution'
            !         * 'f0'         : The objective function value at the starting point 'x_0'
            !         * 'termination': The cause of the termination in the algorithm
            !         * 'counter'    : Gives the values of different counters
            !         * 'CPUtime'    : Gives the CPU time
            !
            ! NOTICE: * The dimension of vectors 'x_0' and 'x_solution' has to be 'user_n' ('user_n' is defined by USER in MODULE functions)
            !         * The dimension of the vector 'counter' has to be 6.
            !         * 'mit', 'mrounds' and 'mrounds_clarke' have to be integers.
            !         * IF ('mit' <= 0) THEN DEFAULT value 1000 is used
            !         * IF ('mrounds' <= 0 ) THEN DEFAULT value 5000 is used.
            !         * 'iprint' has to be -4, -3, -2, -1, 0, 1, 2, 3 or 4 . If it is NOT then DEFAULT value 1 is used. 
            !
            !***********************************************************************************
               IMPLICIT NONE
            !**************************** NEEDED FROM USER *************************************               
               REAL(KIND=dp), DIMENSION(:), INTENT(IN)  :: x_0         ! the starting point
               REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: x_solution  ! the solution obtained to the problem
               
               REAL(KIND=dp), INTENT(OUT) :: f_solution                ! the objective function value at the solution 'x_solution'
               REAL(KIND=dp), INTENT(OUT) :: f0                        ! the objective function value at the starting point 'x_0'
               REAL(KIND=dp), INTENT(OUT) :: CPUtime                   ! the CPU time
                              
               REAL(KIND=dp), INTENT(IN) :: time_limit          ! the CPU time limit
               
               INTEGER, INTENT(INOUT) :: mit                  ! the maximum number of 'main iterations'
               INTEGER, INTENT(INOUT) :: mrounds              ! the maximum number of rounds during one 'main iteration'
               
               INTEGER, INTENT(OUT) :: termination        ! 1 - the stopping condition satisfied
                                                          ! 2 - the biggest possible number of rounds executed in one 'inner iteration'
                                                          ! 3 - the biggest possible number of 'inner iterations' executed in the algorithm!
                                                          ! 4 - the CPU time limit is exceeded
                                                          ! 5 - the problem seem to be unbounded from below.
                                                          ! 6 - the difference in consecutive function values is small.                                                       
                                                          
               INTEGER, DIMENSION(6), INTENT(OUT) :: counter  ! contains the values of different counteres: 
                                                              !   counter(1) = serious_counter      the number of 'serious steps' executed
                                                              !   counter(2) = iter_counter         the number of 'inner iterations' executed
                                                              !   counter(3) = f_counter            the number of function values evaluated for DC component in the algorithm
                                                              !   counter(4) = subgrad1_counter     the number of subgradients calculated for f_1 in the algorithm
                                                              !   counter(5) = subgrad2_counter     the number of subgradients calculated for f_2 in the algortihm
                                                              !   counter(6) = stop_cond_counter    the number of times tau is decreased
                                                                                      
               INTEGER, INTENT(INOUT) :: iprint ! variable that specifies print option:
                                                !   iprint = 0 : print is suppressed
                                                !   iprint = 1 : basic print of final result 
                                                !   iprint = -1: basic print of final result (without the solution vector)
                                                !   iprint = 2 : extended print of final result 
                                                !   iprint = -2: extended print of final result (without the solution vector)
                                                !   iprint = 3 : basic print of intermediate results and extended print of final results
                                                !   iprint = -3: basic print of intermediate results and extended print of final results (without the solution vector)
                                                !   iprint = 4 : extended print of intermediate results and extended print of final results 
                                                !   iprint = -4: extended print of intermediate results and extended print of final results (without the solution vectors)
           
               LOGICAL, INTENT(IN) :: consecutive             ! .TRUE. if algorithm can be terminated if the difference between consecutive function values is small
           
           !***************************** LOCAL VARIABLES ************************************             
               REAL(KIND=dp), DIMENSION(user_n) :: x_current     ! the current iteration point (the dimension 'user_n' is the number of variables)
               REAL(KIND=dp), DIMENSION(user_n) :: x_new         ! the new iteration point used in line search to store an auxiliary result         
               REAL(KIND=dp), DIMENSION(user_n) :: x_tau         ! the point x_current + tau * dd             
               REAL(KIND=dp), DIMENSION(user_n) :: x_step        ! the point used in step-size detemination             
           
               REAL(KIND=dp), DIMENSION(user_n) :: grad          ! Combined subgradient (the dimenison 'user_n' is the length of subgradient)       
               REAL(KIND=dp), DIMENSION(user_n) :: grad_agg      ! Combined aggregated subgradient (the dimenison 'user_n' is the length of subgradient)       
               REAL(KIND=dp), DIMENSION(user_n) :: grad1, grad2  ! subgradients (the dimenison 'user_n' is the length of subgradient)       
               REAL(KIND=dp), DIMENSION(user_n) :: dd            ! the search direction obtained from the 'main iteration': dd=x_new -x_current
                                                                 ! (the dimension 'user_n' is the length of subgradient)                                                                         
                                                                 
               REAL(KIND=dp) :: tau            ! 'tau'=the search control parameter used in the method
               REAL(KIND=dp) :: eps1,eps2      ! 'eps1'=the optimality tolerance and 'eps2'=the optimality tolerance for the inner iteration
               REAL(KIND=dp) :: sigma          ! the descent parameter for tau 
               REAL(KIND=dp) :: sigma2         ! the descent parameter for tau 
               REAL(KIND=dp) :: c1, c2         ! 'c1'=the search control parameter and 'c2'=the line search parameter  
                         
               REAL(KIND=dp) :: lambda         ! the solution for the norm minimization problem    
               REAL(KIND=dp) :: norm           ! the value of the norm minimization problem    
               REAL(KIND=dp) :: n              ! the dimension of the problem 
               
               REAL(KIND=dp) :: step           ! the step-size determined by line search after inner iteration        
               REAL(KIND=dp) :: prod_diff      ! help variable        
                         
               REAL(KIND=dp) :: f_0                         ! the value of the objective function f=f_1-f_2 at the starting point x_0
               REAL(KIND=dp) :: f1_current, f2_current      ! the value of f_1 and f_2 at the current solution  
               REAL(KIND=dp) :: f1_new, f2_new              ! the value of f_1 and f_2 at the iteration point x_new (used in line search to store intermediate results)        
               REAL(KIND=dp) :: f1_tau, f2_tau              ! the value of f_1 and f_2 at the point x_tau (which is obtained in inner iteration)        
               REAL(KIND=dp) :: f1_step, f2_step            ! the values of f_1 and f_2 at the point x_step (used in line search determination)        
               REAL(KIND=dp) :: f_current, f_previous, f_c, f_p  ! the values of f at the point x_step (used in line search determination)        
               REAL(KIND=dp) :: norm_current, norm_previous ! the values of f at the point x_step (used in line search determination)        
               REAL(KIND=dp) :: change                      ! 'change' = f(x_new) - f(x_current)  (i.e. the change in the objective function value)
        
               REAL(KIND=dp) :: start_time, finish_time              ! start and finish CPU time
               REAL(KIND=dp) :: middle_time                          ! middle CPU time
               
               REAL(KIND=dp) :: elapsed_time                  ! elapsed 'clock' time in seconds
               INTEGER :: clock_start, clock_end, clock_rate  ! start and finish 'clock' time   
              
               INTEGER  :: line_counter        ! the number of line search rounds
               
               INTEGER  :: iter_counter        ! the number of inner iterations executed
               INTEGER  :: serious_counter     ! the number of serious steps executed
               INTEGER  :: inner_iter_counter  ! the number of inner iterations solved during the execution of the subgradient algorithm
               INTEGER  :: stop_cond_counter   ! the number of times 'Clarke stationary' algorithm is used during the algorithm
               INTEGER  :: f_counter           ! the number of function values evaluated for a DC component during 
                                               ! the execution of the bundle algorithm (same for the DC component f_1 and f_2) (without Clarke stationary algorithm)
               
               INTEGER  :: subgrad1_counter    ! the number of subgradients calculated for f_1 during the bundle algorithm (without Clarke stationary algorithm)
               INTEGER  :: subgrad2_counter    ! the number of subgradients calculated for f_2 during the bundle algorithm (without Clarke stationary algorithm)                 
               
               INTEGER :: reason_for_stop       ! the reason for stop during the 'inner iteration'
                                                ! 0 - a new descent direction found 
                                                ! 5 - the optimality tolerance reached
                                                ! 1 - the stopping condition satisfied
                                                ! 2 - the biggest possible number of rounds executed in the 'inner iteration'
                                                ! 3 - the biggest possible number of 'inner iterations' executed in the algorithm
               
               INTEGER :: i                     ! help variables 
                                                    
               LOGICAL :: extra                 ! .TRUE. if the subgradient method SUBGDC can be terminated                                             
               LOGICAL :: stop_alg              ! .TRUE. if the subgradient method SUBGDC can be terminated                                             
               LOGICAL :: stop_inner            ! .TRUE. if the current inner iteration can be terminated                                             
               LOGICAL :: stop_linesearch       ! .TRUE. if the current inner iteration can be terminated                                             


               CHARACTER*30 outfi/'results.txt'/
               OPEN(40,file=outfi)             
               
               CALL cpu_time(start_time)                 ! Start CPU timing     
               CALL cpu_time(middle_time)                ! Middle CPU timing     

               CALL SYSTEM_CLOCK(COUNT_RATE=clock_rate) ! Find the rate
               CALL SYSTEM_CLOCK(COUNT=clock_start)     ! Start timing 'clock' time            
                              
                 
           !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^     
           ! <>  <>  <>  <>  <>  <>  <>  MAIN ITERATION STARTS  <>  <>  <>  <>  <>  <>  <>  <>
           !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^      
           
               !*************** PARAMETER VALUES FROM USER ARE CHECKED **********************
               !
               !     IF VALUES ARE INCORRECT THEN DEFAULT VALUES ARE USED FOR PARAMETERS
               !          
               ! ***** descent parameter for TAU 'sigma' *****
               IF ( (user_sigma1<=0.0_dp) .OR. (user_sigma1>=1.0_dp) ) THEN
                   sigma = 0.2_dp
               ELSE       
                   sigma = user_sigma1 
               END IF
               
               ! ***** descent parameter for delta 'sigma2' *****
               IF ( (user_sigma2<=0.0_dp) .OR. (user_sigma2>1.0_dp+0.000000001_dp) ) THEN
                   sigma2 = 0.99_dp
               ELSE       
                   sigma2 = user_sigma2 
               END IF              
              
               ! ***** optimality tolerance 'eps1' *****
               IF (user_eps <= 0.0_dp) THEN
                   eps1 = (10.0_dp)**(-5)
               ELSE 
                   eps1 = user_eps
               END IF
               
               ! ***** optimality tolerance for inner iteration 'eps2' *****
               IF (user_delta <= 0.0_dp) THEN
                  eps2 = (10.0_dp)**(-5)
               ELSE 
                   eps2 = user_delta
               END IF              
              
              
               ! ***** the initial value for the search control parameter 'tau' *****
               IF (user_tau <= 0.0_dp) THEN
                  IF (user_n < 100) THEN
                     tau = 2.0_dp
                  ELSE IF (user_n<200) THEN
                     tau = 10.0_dp
                  ELSE 
                     tau = 50.0_dp        
                  END IF  
               ELSE
                   tau = user_tau
               END IF
 
               ! ***** search control parameter 'c1' *****               
               IF ((user_c1<=0.0_dp) .OR. (user_c1>1.0_dp)) THEN
                   c1 = 0.2_dp
               ELSE
                   c1 = user_c1
               END IF              
               
               ! *****  line search parameter 'c2' *****               
               IF ((user_c2>c1) .OR. (user_c2<=0.0_dp)) THEN
                   c2 = 0.25_dp*c1
               ELSE
                   c2 = user_c2
               END IF                  
               
               !----------------------------------------------------------------------------
               ! ***** maximum number of inner iterations 'mit' *****
               IF ( mit <= 0 ) THEN
                   mit = 1000
               END IF   
               
               ! ***** maximum number of rounds 'mrounds' in the 'inner iteration' *****
               IF ( mrounds <= 0 ) THEN
                   mrounds = 5000
               END IF                
               
               ! ***** print option 'iprint' *****
               IF ( (iprint < -4 ) .OR. (iprint > 4) ) THEN   !executed if print value is wrong
                   iprint = 1
               END IF              
              !----------------------------------------------------------------------------
              
           !_______________________________________________________________________________
           !************************ STEP 1: PARAMETER INITIALIZATION *********************                
           
               x_current = x_0                           ! the current iteration point is the starting point x_0
               f1_current = f1(x_0, problem1, user_n)    ! the value of the DC component f_1 at x_0
               f2_current = f2(x_0, problem2, user_n)    ! the value of the DC component f_2 at x_0
               f_counter = 1                             ! one function value evaluated for both f_1 and f_2
               
               f_0 = f1_current-f2_current  
               f0 = f_0            
               f_current = f_0
               f_c = f_0
               
               iter_counter = 0             ! the number of 'inner iterations' executed so far is zero
               serious_counter = 0          ! the number of 'serious steps' executed so far is zero
               stop_cond_counter = 0        ! the number of times the stopping condition is tested is zero
               stop_alg = .FALSE.           ! we cannot stop the subgradient algorithm
               
               subgrad1_counter = 0         ! the number of subgradients calculate for f1
               subgrad2_counter = 0         ! the number of subgradients calculate for f2
                               

               IF ( (ABS(iprint) >= 3) .AND. (ABS(iprint) <= 4) ) THEN        ! the basic/extended print of indermediate results     

                    WRITE(*,*) 'Start round:', 0, 'f(x):', f1_current - f2_current  

               END IF              
           !_______________________________________________________________________________
           !************************ STEP 1: END ******************************************    
           
           !----------------------------------------------------------------------------------   
           ! <>  <>  <>  <>  <>  <>  <> REPEATABLE PART BEGINS <>  <>  <>  <>  <>  <>  <>  <>
           !----------------------------------------------------------------------------------  
           ! is repeated until the aggregate subgradient method SUBGDC can be terminated 
            DO WHILE ( ( .NOT. stop_alg) .AND. (iter_counter < mit ) .AND. &
                        & ( (middle_time-start_time) < time_limit ) .AND.  &    ! The CPU time limit is not met
                        & ( (f1_current-f2_current) > -(10.0_dp)**8) )          ! The problem seems to be bounded below
               
              !------------------------------------------------------------------------------- 
              ! THE BEGINNING OF THE INNER ITERATION
              !-------------------------------------------------------------------------------            

               !_______________________________________________________________________________
               !************************ STEP 2: INNER ITERATION INITIALIZATION ***************
                inner_iter_counter = 0       ! the iteration counter for the current 'inner iterations' 
                stop_inner = .FALSE.         ! we cannot stop the inner iteration
                
                norm_current = (10.0_dp)**(9)
                extra = .FALSE.

                ! the initialization of the direction 'dd'
                DO i = 1, user_n
                   n = user_n
                   dd(i) = 1.0_dp/SQRT(n)
                END DO
                x_tau = x_current+tau*dd
               
                grad1 = subgradient_f1(x_tau, problem1, user_n)        ! the subgradient of f_1 at x_current+tau*dd
                subgrad1_counter = subgrad1_counter + 1                ! one more subgradient calculated for f_1
                grad2 = subgradient_f2(x_current, problem2, user_n)    ! the subgradient of f_2 at x_0
                subgrad2_counter = subgrad2_counter + 1                ! one more subgradient calculated for f_2
               
                grad = grad1 - grad2                           ! the combined subgradient 
                grad_agg = grad1 - grad2                       ! the initialization of the combined aggregated subgradient
                                   
               !_______________________________________________________________________________
               !************************ STEP 2: END ******************************************    

               DO WHILE ( (.NOT. stop_inner) .AND. (inner_iter_counter < mrounds )  )  ! the repeatable part of the inner iteration
               !_______________________________________________________________________________
               !************************ STEP 3 and 4: NEW AGGREGATED SUBGRADIENT *************
               
                inner_iter_counter = inner_iter_counter + 1
                
                prod_diff = DOT_PRODUCT(grad-grad_agg,grad-grad_agg)
                ! the solution to the norm minimization problem 
                IF (prod_diff>(10.0_dp)**(-4)) THEN 
                  lambda = -DOT_PRODUCT(grad_agg,grad-grad_agg)/prod_diff
                ELSE
                  lambda = 0.0_dp  
                END IF
                
                IF (lambda < 0.0_dp) THEN 
                   lambda = 0.0_dp
                ELSE 
                   IF (lambda > 1.0_dp) THEN
                       lambda = 1.0_dp
                   ELSE
                       lambda = lambda
                   END IF           
                END IF 
                
                ! the new aggregate subgradient
                grad_agg = lambda * grad + (1.0_dp - lambda)* grad_agg
                
                norm = SQRT(DOT_PRODUCT(grad_agg,grad_agg)) 
                
                IF (( MOD(inner_iter_counter,5)==0) ) THEN 
                   norm_previous = norm_current
                   norm_current = norm
                   IF (ABS(norm_previous-norm_current)<(10.0_dp)**(-2)) THEN 
                       norm = 0.0_dp
                       extra = .TRUE.
                   END IF
                END IF
 
                IF (norm <= eps1) THEN
                   stop_inner = .TRUE.  
                   reason_for_stop = 5               
               !_______________________________________________________________________________
               !************************ STEP 3 and 4: END ************************************    
                   
                ELSE
               !_______________________________________________________________________________
               !************************ STEP 5: SEARCH DIRECTION *****************************                   ´
                   dd = - grad_agg /norm
                   
               !_______________________________________________________________________________
               !************************ STEP 5: END ****************************************** 
               
               !_______________________________________________________________________________
               !************************ STEP 6: DESCENT TEST AND NEW SUBGRADIENT *************                

                   x_tau = x_current + tau*dd

                   f1_tau = f1(x_tau, problem1, user_n)    ! the value of the DC component f_1 at x_tau
                   f2_tau = f2(x_tau, problem2, user_n)    ! the value of the DC component f_2 at x_tau
                   f_counter = f_counter + 1               ! one more function value evaluated for both f_1 and f_2
                   
                   change = f1_tau - f2_tau -(f1_current - f2_current)
                   
                   IF (change <= -c1*tau*norm) THEN
                       stop_inner = .TRUE.
                       reason_for_stop = 0
                   ELSE
                       grad1 = subgradient_f1(x_tau, problem1, user_n)    ! the new subgradient of f_1 at x_current+tau*dd
                       subgrad1_counter = subgrad1_counter + 1            ! one subgradient calculated for f_1  
                       
                       grad = grad1 - grad2                               ! the new combined subgradient 
                       
                       IF (inner_iter_counter >= mrounds) THEN 
                          reason_for_stop = 2
                          stop_alg = .TRUE.
                       END IF 
                   END IF
                               
               !_______________________________________________________________________________
               !************************ STEP 6: END ******************************************         
 
                END IF
                
              END DO    

              !------------------------------------------------------------------------------- 
              ! THE END OF THE INNER ITERATION
              !-------------------------------------------------------------------------------  
              
               IF (reason_for_stop == 5) THEN   ! in this case, a stopping condition needs to be tested
               !_______________________________________________________________________________
               !************************ STEP 7: STOPPING CONDITION ***************************
               
                  stop_cond_counter = stop_cond_counter + 1
                    IF (tau < eps2) THEN
                       stop_alg = .TRUE.
                       reason_for_stop = 1
                    ELSE
                       tau = sigma * tau        ! parameter tau is decreased
                       eps2 = sigma2 * eps2     ! parameter delta is decreased
                    END IF
              
                   iter_counter = iter_counter + 1      ! one more inner iteration executed
                    
               !_______________________________________________________________________________
               !************************ STEP 7: END ******************************************       
               END IF 
               
               IF (reason_for_stop == 0) THEN     ! the descent direction is found             
               !_______________________________________________________________________________
               !************************ STEP 8: LINE SEARCH ********************************** 

                   x_new = x_tau                              ! update of the current iteration point
                   f1_new = f1_tau                            ! update of the function value f_1
                   f2_new = f2_tau                            ! update of the function value f_2
                   
                   stop_linesearch = .FALSE.
                   step = tau 

                   line_counter = 0                   
                   
                   
                   DO WHILE( (.NOT. stop_linesearch) .AND. (line_counter < 20) )                   
                      step = step*2.0_dp
                      x_step = x_current + step * dd
                    
                      f1_step = f1(x_step, problem1, user_n)    ! the value of the DC component f_1 at x_step
                      f2_step = f2(x_step, problem2, user_n)    ! the value of the DC component f_2 at x_step
                      f_counter = f_counter + 1                 ! one function value evaluated for both f_1 and f_2
                
                      change = f1_step - f2_step - (f1_current - f2_current)   

                    IF (change <= -c2*step*norm) THEN 
                      x_new = x_step                 ! a new better point found
                      f1_new = f1_step               ! a new function value for f1
                      f2_new = f2_step               ! a new function value for f2
                      line_counter = line_counter + 1
                    ELSE
                      stop_linesearch = .TRUE.       ! line search can be stopped since no improvement obtained
                    END IF
                   END DO
                   x_current = x_new                 ! update of the current iteration point
                   f1_current = f1_new               ! update of the function value f_1
                   f2_current = f2_new               ! update of the function value f_2

                   iter_counter = iter_counter + 1        ! one inner iteration done
                   serious_counter = serious_counter + 1  ! one serious step done
               
                   WRITE(*,*) 'Iteration:', iter_counter, 'f(x):', f1_current-f2_current, 'Inner_it:', inner_iter_counter 

                                                                         
                  IF ((MOD(iter_counter,5)==0) .AND. (tau > eps2)) THEN 
                     f_previous = f_current
                     f_current = f1_current-f2_current
                     IF (ABS(f_previous-f_current)<(10.0_dp)**(-2)) THEN 
                       tau = sigma * tau
                       !WRITE(*,*) 'extra tau update'
                     END IF
                  END IF
                   
                  IF ((MOD(serious_counter,10)==0)) THEN
                       f_p = f_c
                       f_c = f1_current-f2_current
                       IF (ABS(f_p-f_c)<10.0_dp**(-5)) THEN
                         IF (tau > eps2) THEN 
                           tau = sigma * tau
                         ELSE  
                          stop_alg = .TRUE.
                          reason_for_stop = 6
                         END IF  
                       END IF
                  END IF                   
               !_______________________________________________________________________________
               !************************ STEP 8: END ******************************************                      
               END IF

               CALL cpu_time(middle_time)         ! Middle CPU timing

               
           END DO
           !----------------------------------------------------------------------------------   
           ! <>  <>  <>  <>  <>  <>  <>  REPEATABLE PART ENDS  <>  <>  <>  <>  <>  <>  <>  <> 
           !----------------------------------------------------------------------------------         
           

           termination = reason_for_stop   ! the cause of the termination
           
           IF ( iter_counter >= mit ) THEN  ! the maximum number of 'main iterations' is executed and this causes the termination
               termination = 3 
           END IF
           
           IF ( middle_time-start_time > time_limit ) THEN  ! the CPU time limit is exceeded
               termination = 4 
           END IF

           IF ( (f1_current-f2_current) <= -(10.0_dp)**8 ) THEN  ! the problem seems to be unbounded below
               termination = 5 
           END IF          
           
           x_solution = x_current                ! the solution to the minimization problem
           f_solution = f1_current - f2_current  ! the objective function value at the solution
                   
           ! the values of different counters
           counter(1) = serious_counter 
           counter(2) = iter_counter
           counter(3) = f_counter
           counter(4) = subgrad1_counter
           counter(5) = subgrad2_counter 
           counter(6) = stop_cond_counter
            
           CALL cpu_time(finish_time)         ! Stop CPU timing
           CALL SYSTEM_CLOCK(COUNT=clock_end) ! Stop timing 'clock' time
           
          ! Calculate the elapsed 'clock' time in seconds:
          elapsed_time=(1.0_dp*clock_end-clock_start)/clock_rate     
          CPUtime = finish_time-start_time    

           IF ( (ABS(iprint) == 1) ) THEN ! basic print of the final result
               WRITE(*,*) '---------------------------------------------------------------------------'        
               WRITE(*,*) '           * * * * * BASIC PRINT OF THE FINAL SOLUTION * * * * * '
               WRITE(*,*) '---------------------------------------------------------------------------'            
               WRITE(*,*) 'The cause of termination: ', termination
               IF (iprint > 0) THEN
                   WRITE(*,*) '- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - '
                   DO i = 1, user_n
                       WRITE(*,*) 'x(', i ,')*=',  x_solution(i)
                   END DO
               END IF
               WRITE(*,*) '- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - '               
               WRITE(*,*) 'f*=', f_solution     
               WRITE(*,*) '---------------------------------------------------------------------------' 
               WRITE(*,*) 'CPU time used:', finish_time-start_time                              
               WRITE(*,*) '---------------------------------------------------------------------------'
               WRITE(*,*) 'Elapsed time:', elapsed_time
               WRITE(*,*) '---------------------------------------------------------------------------'            
           END IF
           
           IF ((ABS(iprint) /= 1) .AND. (iprint /= 0)) THEN  ! extended print of the final result
               WRITE(*,*) '---------------------------------------------------------------------------'        
               WRITE(*,*) '        * * * * * EXTENDED PRINT OF THE FINAL SOLUTION * * * * * '
               WRITE(*,*) '---------------------------------------------------------------------------'            
               WRITE(*,*) 'The cause of termination: ', termination
               IF (iprint > 0) THEN
                   WRITE(*,*) '- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - '
                   DO i = 1, user_n
                       WRITE(*,*) 'x(', i ,')*=',  x_solution(i)
                   END DO
               END IF
               WRITE(*,*) '- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - '               
               WRITE(*,*) 'f*=', f_solution                    
               WRITE(*,*) '- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - '   
               WRITE(*,*) 'The number of serious steps:', serious_counter
               WRITE(*,*) 'The number of inner iterations:', iter_counter             
               WRITE(*,*) 'The number of function values evaluated for DC component f_1:', f_counter 
               WRITE(*,*) 'The number of function values evaluated for DC component f_2:', f_counter 
               WRITE(*,*) 'The number of subgradients calculated for DC component f_1:', subgrad1_counter 
               WRITE(*,*) 'The number of subgradients calculated for DC component f_2:', subgrad2_counter 
               WRITE(*,*) '---------------------------------------------------------------------------'
               WRITE(*,*) 'The number of times tau is decreased:', stop_cond_counter  
               WRITE(*,*) '---------------------------------------------------------------------------'     
               WRITE(*,*) 'CPU time used:', finish_time-start_time                          
               WRITE(*,*) '--------------------------------------------------------------------------- '
               WRITE(*,*) 'Elapsed time:', elapsed_time
               WRITE(*,*) '---------------------------------------------------------------------------'                
           END IF 
           
       WRITE(40,*) 'The value of the objective at initial point:  ', f_0
       WRITE(40,*) 'The value of the objective at final point:    ',f_solution 
       WRITE(40,*) '-----------------------------------------------------------------------------&
         &---------------------------'            
      WRITE(40,42) iter_counter
  42  FORMAT(' The total number of times the inner iteration is executed solved:         ',i12)  
      WRITE(40,43) f_counter  
  43  FORMAT(' The total number of the objective evaluations in main iterations:  ',i12)
      WRITE(40,44) subgrad1_counter
  44  FORMAT(' The total number of subgradient evaluations for f1 in main iterations:',i12)
      WRITE(40,45) subgrad2_counter  
  45  FORMAT(' The total number of subgradient evaluations for f2 in main iterations:',i12)
        WRITE(40,*) '-----------------------------------------------------------------------------&
         &---------------------------'        
      WRITE(40,46) finish_time-start_time
  46  FORMAT(' The CPU time:                                 ',f10.4)
      WRITE(40,47) elapsed_time
  47  FORMAT(' The elapsed time:                             ',f10.4)

           
           CLOSE(40)   
           
           END SUBROUTINE SUBGDC_algorithm      
        !.......................................................................................           
        ! <>  <>  <>  <>  <>  <>  <>  <>  <>  <>  <>  <>  <>  <>  <>  <>  <>  <>  <>  <>  <>  <>   
        ! _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _  
        !| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | |        
        !** END ** END ** END ** END ** END ** END ** END ** END ** END ** END ** END ** END ** 
        !|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|
        



      END MODULE subgdc