MODULE mdbdc !*..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..* !| .**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**. | !| | | | !| | MDBDC - THE MULTIOBJECTIVE DOUBLE BUNDLE METHOD | | !| | FOR NONSMOOTH DC OPTIMIZATION | | !| | (version 1) | | !| | | | !| | | | !| | Features : | | !| | | | !| | * Possibility to use simple stepsize determination after | | !| | each 'main iteration'. | | !| | | | | | !| | * During each round of 'main iteration' OpenMP can be utilized | | !| | to calculate subproblems in parallel. However if you DO NOT | | !| | WANT to use this feature then | | !| | 1) in Makefile DELETE '-fopenmp' from line 5 | | !| | 2) in mdbdc.f95 COMMENT lines 615-618 | | !| | | | !| | | | !| | | | !| .**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**. | !*..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..* !| | !| Utilizes the new version of PLQDF1 by Ladislav Luksan as a quadratic solver. | !| | !| | !| .**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**. | !*..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..* USE omp_lib USE constants, ONLY : dp ! Double precision (i.e. accuracy) USE bundle1 ! The BUNDLE of the DC component f_1 USE bundle2 ! The BUNDLE of the DC component f_2 USE functions ! Contains INFORMATION from the USER IMPLICIT NONE EXTERNAL PLQDF1 ! The QUADRATIC SOLVER by Ladislav Luksan CONTAINS !*..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..* !| .**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**. | !| | | | !| | CONTAINS SUBROUTINES: | | !| | | | !| | bundle_algorithm : The double bundle method for nonsmooth | | !| | multiobjective DC optimization. | | !| | main_iteration : The main iteration algorithm needed in the double | | !| | bundle method for DC optimization. | | !| | guaranteeing_clarke : The algorithm guaranteeing Clarke stationarity | | !| | for a solution obtained. | | !| | | | !| | quadratic_solver : The solver for the quadratic norm minimization | | !| | problem. Needed in the Clarke stationary algorithm | | !| | subproblem_solver : The solver for subproblems in the search direction | | !| | problem. Needed in the main iteration algorithm. | | !| | | | !| | | | !| | | | !| .**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**. | !*..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..**..* ! _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ !| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | ! START ** START ** START ** START ** START ** START ** START ** START ** START ** START !|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_| ! <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> !*************************************************************************************** ! ---------------------------------------------------------------------------------- | ! | | | ! | THE DOUBLE BUNDLE ALGORITHM FOR NONSMOOTH DC OPTIMIZATION | | ! | | | ! ---------------------------------------------------------------------------------- | !*************************************************************************************** SUBROUTINE bundle_algorithm( x_0, x_solution, f_solution, mit,& & mrounds, mrounds_clarke, termination, counter, & & stepsize_used, iprint, scale_func) ! ! Solves the unconstrained nonsmooth DC minimization problem ! ! INPUT: * 'x_0' : A starting point ! * 'mit' : The maximum number of 'main iterations' ! * 'mrouds' : The maximum number of rounds during one 'main iteration' ! * 'mrouds_clarke' : The maximum number of rounds during one 'Clarke stationary' algorithm ! * 'stepsize_used' : If .TRUE. then simple stepsize determination is used in the algorithm ! * 'iprint' : Specifies the print ! ! OUTPUT: * 'x_solution' : The solution obtained to the minimizatin problem ! * 'f_solution' : The objective function value at the solution 'x_solution' ! * 'termination': The cause of the termination in the algorithm ! * 'counter' : Gives the values of different counters ! ! NOTICE: * The dimension of vectors 'x_0' and 'x_solution' has to be 'user_n' defined by USER in MODULE functions. ! * The dimension of the vector 'counter' has to be 5. ! * '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. ! * IF ('mrounds_clarke' <= 0 ) THEN DEFAULT value 5000 is used. ! * 'iprint' has to be -4, -3, -2, -1, 0, 1, 2, 3 or 4 (5). If it is NOT then DEFAULT value 1 is used. ! * If 'iprint = 5' then everything is printed step by step (this is NOT recommended!!!) ! !*********************************************************************************** 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' 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(INOUT) :: mrounds_clarke ! the maximum number of rounds during one 'Clarke stationary' algorithm INTEGER, INTENT(OUT) :: termination ! 1 - the stopping condition is satisfied (i.e. approximate Clarke stationarity) ! 2 - the approximate stopping condition is satisfied (i.e. the step-length beta* < beta_0) ! 3 - the maximum number 'mrounds' of rounds executed in one main iteration ! 4 - the maximum number of 'main iterations' is executed ! 5 - the maximum number 'mrounds_clarke' of rounds executed in one 'Clarke stationary' alqorithm INTEGER, DIMENSION(8), INTENT(OUT) :: counter ! contains the values of different counteres: ! counter(1) = iter_counter the number of 'main iterations' executed ! counter(2) = subprob_counter the number of subproblems solved ! counter(3) = f_counter the number of function values evaluated for DC component ! counter(4) = subgrad1_counter the number of subgradients calculated for f_1 ! counter(5) = subgrad2_counter the number of subgradients calculated for f_2 !-------------------------------------------------------------------------------------------------------------------------- ! counter(6) = stop_cond_counter the number of times 'Clarke stationary algorithm' is executed ! counter(7) = clarke_f_counter the number of function values evaluated for f (in 'Clarke stationary algorithms') ! counter(8) = clarke_sub_counter the number of subgradients caluculated for f (in 'Clarke stationary algorithms') LOGICAL, INTENT(IN) :: stepsize_used ! .TRUE. if simple stepsize determination is used in the algorithm. Otherwise .FALSE. LOGICAL, INTENT(IN) :: scale_func ! .TRUE. if scaling procedure is used in the algorithm. Otherwise .FALSE. 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) ! iprint = 5 : prints each step of the bundle algorithm (i.e. everything is printed step by step) !***************************** LOCAL VARIABLES ************************************ TYPE(kimppu1) :: B1 ! The bundle B_1 for the DC component f_1 TYPE(kimppu2) :: B2 ! The bundle B_2 for the DC component f_2 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 (obtained from the previous 'main iteration') REAL(KIND=dp), DIMENSION(user_n) :: vect ! a 'help' vector REAL(KIND=dp), DIMENSION(number_of_func, user_n) :: grad1_k_all ! the subgradients of the DC components f1 and g1 at x_k REAL(KIND=dp), DIMENSION(number_of_func, user_n) :: grad2_k_all ! the subgradients of the DC components f2 and g2 at x_k REAL(KIND=dp), DIMENSION(number_of_func, user_n) :: grad1_new_all ! the subgradients of the DC components f1 and g1 at x_new REAL(KIND=dp), DIMENSION(number_of_func, user_n) :: grad2_new_all ! the subgradients of the DC components f2 and g2 at x_new REAL(KIND=dp), DIMENSION(number_of_func) :: fg1_k_all ! the values of the DC components f1 and g1 at x_k REAL(KIND=dp), DIMENSION(number_of_func) :: fg2_k_all ! the values of the DC components f2 and g2 at x_k REAL(KIND=dp), DIMENSION(number_of_func) :: fg_k_all ! the values of the DC components f and g at x_k REAL(KIND=dp), DIMENSION(number_of_func) :: fg1_new_all ! the values of the DC components f1 and g1 at x_new REAL(KIND=dp), DIMENSION(number_of_func) :: fg2_new_all ! the values of the DC components f2 and g2 at x_new REAL(KIND=dp), DIMENSION(number_of_func) :: factor ! the scaling factors for f_i and g_l REAL(KIND=dp), DIMENSION(number_of_func) :: factor_old ! the old scaling factors for f_i and g_l INTEGER, DIMENSION(number_of_func) :: fac_ind ! the scaling factor indices for f_i and g_l INTEGER, DIMENSION(number_of_func) :: ero_ind ! the difference between scaling factor indices for f_i and g_l REAL(KIND=dp), DIMENSION(number_of_obj) :: f_0_all ! the values of f_i at x_0 REAL(KIND=dp), DIMENSION(number_of_func) :: AandB_k_all ! the values of the components A_i and B_l at x_k REAL(KIND=dp), DIMENSION(number_of_func) :: AandB_new_all ! the values of the components A_i and B_l at x_new 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) :: c1, c2, c3 ! Decrease parameter used in 'main iteration' REAL(KIND=dp) :: crit_tol, eps ! 'crit_tol'=the stopping tolerance and 'eps'=the proximity measure REAL(KIND=dp) :: m ! the descent parameter used in 'main iteration' REAL(KIND=dp) :: m_extra ! the extra descent parameter used in 'main iteration' REAL(KIND=dp) :: r_dec, r_inc ! 'r_dec'=the decrease parameter and 'R_inc'=the increase parameter REAL(KIND=dp) :: alpha ! help variable when linearization error is calculated REAL(KIND=dp) :: decrease ! approximation of decrease REAL(KIND=dp) :: dec_tol ! Decraese tolerance used in main iteration REAL(KIND=dp) :: step_tol ! Step-length tolerance used in 'Clarke stationary' algorithm REAL(KIND=dp) :: m_clarke ! the descent parameter used in 'Clarke stationary' algorithm REAL(KIND=dp) :: t ! the value of the proximity parameter REAL(KIND=dp) :: t_min ! the lower bound for the proximity parameter REAL(KIND=dp) :: t_max ! the upper bound for the proximity parameter REAL(KIND=dp) :: norm1 ! the norm of the subgradient for H_1 at the current iteration point REAL(KIND=dp) :: max_norm ! the maximum norm of subgradient in B_2 INTEGER :: i_t ! the index used in the update process of the proximity parameter REAL(KIND=dp) :: H1_current, H2_current ! the value of H_1 and H_2 at the current solution REAL(KIND=dp) :: H1_new, H2_new ! the value of H_1 and H_2 at the new iteration point x_new (obtained from the previous main iteration) REAL(KIND=dp) :: change ! 'change' = H(x_new) - H(x_current) (i.e. the change in the objective function value) REAL(KIND=dp) :: change1 ! 'change1' = H1(x_new) - H1(x_current) (i.e. the change in the DC component H1) REAL(KIND=dp) :: change2 ! 'change2' = H2(x_new) - H2(x_current) (i.e. the change in the DC component H2) REAL(KIND=dp) :: norm ! the distance between two consecutive iteration points REAL(KIND=dp) :: start_time, finish_time ! start and finish CPU time REAL(KIND=dp) :: start_time_main_it, finish_time_main_it ! start and finish CPU time in one 'main iteration' REAL(KIND=dp) :: elapsed_time ! elapsed 'clock' time in seconds INTEGER :: clock_start, clock_end, clock_rate ! start and finish 'clock' time INTEGER :: size_b1 , size_b2 ! The biggest possible size of the bundles B_1 and B_2 ! Parameter used in proximity parameter update procedure REAL(KIND=dp) :: u REAL(KIND=dp) :: u_k REAL(KIND=dp) :: t_next REAL(KIND=dp) :: u_next REAL(KIND=dp) :: u_min LOGICAL :: extra_dec INTEGER :: iter_counter ! the number of main iterations executed INTEGER :: subprob_counter ! the number of subproblems solved during the execution of the bundle 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 :: clarke_f_counter ! the number of function values evaluated for f in 'Clarke stationary' algorithm during the execution of the bundle algorithm INTEGER :: clarke_sub_counter ! the number of subgradients evaluated for f in 'Clarke stationary' algorithm during the execution of the bundle algorithm INTEGER :: subgrad1_counter ! the number of subgradients calculated for f_1 during the bundle algorithm INTEGER :: subgrad2_counter ! the number of subgradients calculated for f_2 during the bundle algorithm INTEGER :: help_mit_iter_counter ! the number of iteration rounds executed in one 'main iteration' INTEGER :: help_subprob_counter ! the number of subproblems solved in one 'main iteration' INTEGER :: help_f_counter ! the number of function values evaluated for a DC component in one 'main iteration' (same for f_1 and f_2) INTEGER :: help_subgrad1_counter ! the number of subgradients calculated for f_1 in one 'main iteration' INTEGER :: help_subgrad2_counter ! the number of subgradients calculated for f_2 in one 'main iteration' INTEGER :: help_stop_cond_counter ! the number of times 'Clarke stationarity' was tested during the 'main iteration' INTEGER :: help_clarke_f_counter ! the number of function values evaluated for f_i in 'Clarke stationary' algorithm INTEGER :: help_clarke_sub_counter ! the number of subgradients evaluated for f_i in 'Clarke stationary' algorithm INTEGER :: reason_for_stop ! the reason for stop during the 'main iteration' ! 0 - a new iteration point found ! 1 - the stopping condition satisfied (i.e. approximate Clarke stationarity) ! 2 - approximate stopping condition satisfied (i.e. the step-length beta* < beta_0) ! 3 - the biggest possible number of rounds executed in the 'main iteration' ! 5 - the biggest possible number of rounds executed in the 'Clarke stationary' algorithm INTEGER :: max_threads ! the maximum number of threads that can be used in parallellization INTEGER :: threads ! the number of threads used in parallellization INTEGER :: max_sub_prob ! the maximum number of subproblems solved at the same time INTEGER :: i, j INTEGER :: ind INTEGER :: small_ind LOGICAL :: stop_alg ! .TRUE. if the proximal bundle algorithm can be terminated CHARACTER*30 outfi/'results.txt'/ OPEN(40,file=outfi) CALL cpu_time(start_time) ! Start 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 LOOKED THROUGH ******************** ! ! IF VALUES ARE SET WRONG THEN DEFAULT PARAMETER VALUES ARE USED ! ! ***** descent parameter 'm' ***** IF ( (user_m<=0.0_dp) .OR. (user_m>=1.0_dp) ) THEN m = 0.2_dp ELSE m = user_m END IF ! ***** extra descent parameter 'm' ***** IF ( (user_m_extra<=m) .OR. (user_m_extra>=1.0_dp) ) THEN m_extra = 2.0_dp*m IF (m_extra > 1.0_dp) THEN m_extra = 0.5_dp*(m+1.0_dp) END IF ELSE m_extra = user_m_extra END IF ! ***** decrease parameter 'c1' ***** IF ( (user_c1<=0.0_dp) .OR. (user_c1>=1.0_dp) ) THEN c1 = 0.5_dp ELSE c1 = user_c1 END IF ! ***** decrease parameter 'c2' ***** IF ( (user_c2<-0.0000000000001_dp) .OR. (user_c2>=1.0_dp) ) THEN IF (number_of_obj==2) THEN IF (user_n <= 10) THEN c2 = 0.4_dp ELSE IF (user_n > 10 .AND. user_n < 100) THEN c2 = 0.25_dp ELSE IF (user_n>99 .AND. user_n < 200) THEN c2 = 0.1_dp ELSE IF (user_n > 199 .AND. user_n < 500) THEN c2 = 0.01_dp ELSE IF (user_n > 499 ) THEN c2 = 0.001_dp END IF ELSE IF (user_n < 100) THEN c2 = 0.5_dp ELSE IF (user_n>99 .AND. user_n < 200) THEN c2 = MIN(0.1_dp*(number_of_obj-1),0.5_dp) ELSE IF (user_n > 199 .AND. user_n < 500) THEN c2 = MIN(0.01_dp*(number_of_obj-1),0.5_dp) ELSE IF (user_n > 499 ) THEN c2 = MIN(0.001_dp*(number_of_obj-1),0.5_dp) END IF END IF ELSE c2 = user_c2 END IF ! ***** decrease parameter 'c3' ***** IF ( (user_c3<-0.0000000000000001_dp) .OR. (user_c3>=1.0_dp) ) THEN c3 = 0.1_dp ELSE c3 = user_c3 END IF ! ***** stopping tolerance 'crit_tol' ***** IF (user_crit_tol <= 0.0_dp .OR. user_crit_tol >= 1.0_dp ) THEN crit_tol = (10.0_dp)**(-5) ELSE crit_tol = user_crit_tol END IF ! ***** proximity parameter 'eps' ***** IF (user_eps <= 0.0_dp) THEN eps = 0.00005_dp ELSE eps = user_eps END IF ! ***** decrease parameter 'r_dec' ***** IF ((user_r_dec <= 0.0_dp) .OR. (user_r_dec >= 1.0_dp) ) THEN IF( user_n < 10) THEN r_dec = 0.75_dp ELSE IF ( user_n >= 10 .AND. user_n <300 ) THEN r_dec = user_n/(user_n+5.0_dp) !WRITE(*,*) r_dec ELSE IF ( user_n >= 300 ) THEN r_dec = 0.99 END IF ELSE r_dec = user_r_dec END IF ! ***** increase parameter 'r_inc' ***** IF ( user_r_inc <= 1.0_dp ) THEN r_inc = (10.0_dp)**10 ELSE r_inc = user_r_inc END IF ! ***** bundle B1 size ***** IF ( user_size_b1 <= 0 ) THEN size_b1 = MIN((user_n+5)*number_of_func,1000) ELSE size_b1 = user_size_b1 END IF ! ***** bundle B2 size ***** IF ( user_size_b2 <= 0 ) THEN size_b2 = 3 ELSE size_b2 = user_size_b2 END IF !-------------------------------------------------------------------------- ! ***** Descent parameter 'm_clarke' ***** IF ((user_m_clarke<=0.0_dp) .OR. (user_m_clarke>=1.0_dp)) THEN m_clarke = 0.01_dp ELSE m_clarke = user_m_clarke END IF ! ***** Step-length tolerance 'step_tol' ***** IF (user_step_tol <= 0.0_dp) THEN step_tol = 0.0001_dp ELSE step_tol = user_step_tol END IF ! ***** Decrease tolerance 'dec_tol' ***** IF (user_dec_tol > 0.0_dp) THEN IF (user_n <=100) THEN dec_tol = 0.0_dp ELSE IF (user_n>100 .AND. user_n <250) THEN dec_tol = -(10.0_dp)**(-5) ELSE IF (user_n>=250 ) THEN dec_tol = -(10.0_dp)**(-4) END IF ELSE dec_tol = user_dec_tol END IF !---------------------------------------------------------------------------- ! ***** maximum number of main iterations 'mit' ***** IF ( mit <= 0 ) THEN mit = 1000 END IF ! ***** maximum number of rounds 'mrounds' in the 'main iteration' ***** IF ( mrounds <= 0 ) THEN mrounds = 5000 END IF ! ***** maximum number of rounds 'mrounds' in the 'main iteration' ***** IF ( mrounds_clarke <= 0 ) THEN mrounds_clarke = 5000 END IF ! ***** print option 'iprint' ***** IF ( (iprint < -4 ) .OR. (iprint > 5) ) THEN !executed if print value is wrong iprint = 1 END IF !---------------------------------------------------------------------------- !_______________________________________________________________________________ !************************ STEP 0: PARAMETER INITIALIZATION ********************* factor = 1.0_dp ! the initial scaling factor x_current = x_0 ! the current iteration point is the starting point x_0 DO i = 1, number_of_obj ! the values of DC components f1_i and f2_i at x_0 fg1_k_all(i) = f1(x_0,problem1(i),factor(i)) fg2_k_all(i) = f2(x_0,problem2(i),factor(i)) END DO DO i = number_of_obj+1, number_of_func ! the values of DC components g1_l and g2_l at x_0 fg1_k_all(i) = g1(x_0,problem1(i),factor(i)) fg2_k_all(i) = g2(x_0,problem2(i),factor(i)) END DO fg_k_all = fg1_k_all - fg2_k_all small_ind = 1 DO i = 2, number_of_obj IF (ABS(fg_k_all(small_ind))> ABS(fg_k_all(i))) THEN small_ind = i END IF END DO IF (scale_func) THEN DO i = 1, number_of_obj IF ( ABS(fg_k_all(i)) < (10.0_dp)**(0)) THEN fac_ind(i) = 0 ELSE IF ( ABS(fg_k_all(i)) < (10.0_dp)**(1)) THEN fac_ind(i) = 1 ELSE IF ( ABS(fg_k_all(i)) < (10.0_dp)**(2)) THEN fac_ind(i) = 2 ELSE IF ( ABS(fg_k_all(i)) < (10.0_dp)**(3)) THEN fac_ind(i) = 3 ELSE IF ( ABS(fg_k_all(i)) < (10.0_dp)**(4)) THEN fac_ind(i) = 4 ELSE IF ( ABS(fg_k_all(i)) < (10.0_dp)**(5)) THEN fac_ind(i) = 5 ELSE IF ( ABS(fg_k_all(i)) < (10.0_dp)**(6)) THEN fac_ind(i) = 6 ELSE IF ( ABS(fg_k_all(i)) < (10.0_dp)**(7)) THEN fac_ind(i) = 7 ELSE IF ( ABS(fg_k_all(i)) < (10.0_dp)**(8)) THEN fac_ind(i) = 8 ELSE IF ( ABS(fg_k_all(i)) < (10.0_dp)**(9)) THEN fac_ind(i) = 9 ELSE IF ( ABS(fg_k_all(i)) < (10.0_dp)**(10)) THEN fac_ind(i) = 10 ELSE IF ( ABS(fg_k_all(i)) < (10.0_dp)**(11)) THEN fac_ind(i) = 11 ELSE IF ( ABS(fg_k_all(i)) < (10.0_dp)**(12)) THEN fac_ind(i) = 12 ELSE IF ( ABS(fg_k_all(i)) < (10.0_dp)**(13)) THEN fac_ind(i) = 13 ELSE IF ( ABS(fg_k_all(i)) < (10.0_dp)**(14)) THEN fac_ind(i) = 14 END IF END DO DO i = 1, number_of_obj ero_ind(i) = fac_ind(i)-fac_ind(small_ind) IF(ero_ind(i) > 1) THEN ero_ind(i) = ero_ind(i)-1 END IF factor(i) = (10.0_dp)**(-ero_ind(i)) END DO END IF DO i = 1, number_of_obj ! the values of DC components f1_i and f2_i at x_0 fg1_k_all(i) = f1(x_0,problem1(i),factor(i)) fg2_k_all(i) = f2(x_0,problem2(i),factor(i)) f_0_all(i) = fg1_k_all(i) - fg2_k_all(i) ! the values of the objectives f_i at x_0 END DO DO i = number_of_obj+1, number_of_func ! the values of DC components g1_l and g2_l at x_0 fg1_k_all(i) = g1(x_0,problem1(i),factor(i)) fg2_k_all(i) = g2(x_0,problem2(i),factor(i)) END DO f_counter = 1 ! one function value evaluated for each DC component f1_i, f2_i, g1_l and g2_l DO i = 1, number_of_obj ! the values of subgradients for DC components f1_i and f2_i at x_0 grad1 = subgradient_f1(x_0, problem1(i),factor(i)) grad2 = subgradient_f2(x_0, problem2(i),factor(i)) DO j = 1, user_n grad1_k_all(i,j) = grad1(j) grad2_k_all(i,j) = grad2(j) END DO END DO DO i = number_of_obj+1, number_of_func ! the values of subgradients for DC components g1_l and g2_l at x_0 grad1 = subgradient_g1(x_0, problem1(i),factor(i)) grad2 = subgradient_g2(x_0, problem2(i),factor(i)) DO j = 1, user_n grad1_k_all(i,j) = grad1(j) grad2_k_all(i,j) = grad2(j) END DO END DO subgrad1_counter = 1 ! one subgradient calculated for each DC component f1_i and g1_l subgrad2_counter = 1 ! one subgradient calculated for each DC component f2_i and g2_l DO i = 1, number_of_func ! the values of components A_i and B_l for H_1 at x_0 AandB_k_all(i) = AorB_i(fg1_k_all, fg2_k_all, fg1_k_all, fg2_k_all ,i) END DO H1_current = H1(AandB_k_all) ! the value of the DC component H_1 at x_0 H2_current = H2(fg2_k_all) ! the value of the DC component H_2 at x_0 ind = index_H1(AandB_k_all) ! the index of the component A_i or B_l yielding the subgradient of H_1 grad1 = subgradient_AorB_i(grad1_k_all, grad2_k_all, ind) ! the subgradient of H_1 at x_0 grad2 = subgradient_H2(grad2_k_all) ! the subgradient of H_2 at x_0 ! The bundles B_1 and B_2 are initialized CALL init_bundle_b1(B1, size_b1, user_n) CALL init_bundle_b2(B2, size_b2, user_n) ! The first bundle element is added into the bundle B_1 and B_2 (i.e. the one corresponding to the starting point) CALL add_first_element_b1(B1, grad1, ind) CALL add_first_element_b2(B2, grad2) DO j = 1, number_of_func IF ( j/= ind) THEN grad1 = subgradient_AorB_i(grad1_k_all, grad2_k_all, j) ! the subgradients of A_i and B_l at x_0 alpha = 0 CALL add_element_b1(B1, grad1, alpha, j) END IF END DO iter_counter = 0 ! the number of 'main iterations' executed so far is zero subprob_counter = 0 ! the number of 'subproblems' solved so far is also zero stop_cond_counter = 0 ! the number of times 'Clarke stationary' algorithm is used is zero stop_alg = .FALSE. ! we cannot stop the proximal bundle algorithm clarke_f_counter = 0 ! the initialization of counter for 'Clarke stationary' algorithm clarke_sub_counter = 0 ! the initialization of counter for 'Clarke stationary' algorithm ! -- Proximity parameter initialization vect = give_subgrad_b1(B1,0) ! the vector \bxi_1(x_k) norm1 = SQRT(DOT_PRODUCT(vect,vect)) ! the value of the norm ||\bxi_1(x_k)|| max_norm = max_norm_value(B2) ! the value of the maximum subgradient norm in the bundle B_2 t_min = (0.5_dp * r_dec * eps) / ( norm1 + max_norm ) ! the lower bound for t t_max = r_inc * t_min ! the upper bound for t i_t = 0 ! index initialization IF (user_n < 250) THEN t = 1.0_dp ! the parameter t is selected ELSE IF (number_of_obj == 2) THEN t = 2.0_dp ELSE t = 10.0_dp END IF END IF ! --- --- --- Needed in OpenMP when we use PARALLELLIZATION --- --- --- max_threads = omp_get_max_threads() max_sub_prob = give_max_size_b2(B2)+1 threads = MIN(max_threads, max_sub_prob) CALL omp_set_num_threads(threads) ! ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- IF ( (ABS(iprint) >= 3) .AND. (ABS(iprint) <= 5) ) THEN ! the basic/extended print of indermediate results WRITE(*,*) 'Main Iter:', 0, 'f(x):', fg1_k_all - fg2_k_all ! WRITE(*,*) 'x:', x_current ! WRITE(*,*) '------------------------------------------------------------------------' END IF !_______________________________________________________________________________ !************************ STEP 0: END ****************************************** !---------------------------------------------------------------------------------- ! <> <> <> <> <> <> <> REPEATABLE PART BEGINS <> <> <> <> <> <> <> <> !---------------------------------------------------------------------------------- DO WHILE ( ( .NOT. stop_alg) .AND. (iter_counter < mit ) ) ! is repeated until the bundle algorithm can be terminated !_______________________________________________________________________________ !************************ STEP 1: MAIN ITERATION ******************************* iter_counter = iter_counter + 1 ! a new 'main iteration' is executed IF(iprint==5) THEN !prints everything WRITE(*,*) '-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*' WRITE(*,*) ' MAIN ITERATION', iter_counter, ' STARTS' WRITE(*,*) '-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*' WRITE(*,*) ' ' END IF CALL cpu_time(start_time_main_it) CALL main_iteration( x_current, H1_current, H2_current,& & fg1_k_all, fg2_k_all, x_new, H1_new, H2_new, & & fg1_new_all, fg2_new_all, f_0_all, AandB_k_all, m_extra, & & B1, B2, c1, c2, c3, t, t_min, t_max, i_t, decrease, & & crit_tol, eps, m, r_dec, r_inc, m_clarke, step_tol, mrounds,& & mrounds_clarke, iprint, reason_for_stop, help_mit_iter_counter,& & help_subprob_counter, help_f_counter, help_subgrad1_counter, & & help_subgrad2_counter, help_stop_cond_counter, & & help_clarke_f_counter, help_clarke_sub_counter, & & stepsize_used, extra_dec, factor, dec_tol) IF(iprint==5) THEN !prints everything WRITE(*,*) '-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*' WRITE(*,*) ' MAIN ITERATION', iter_counter, ' ENDS' WRITE(*,*) '-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*' WRITE(*,*) ' ' END IF CALL cpu_time(finish_time_main_it) ! Different counters are updated subprob_counter = subprob_counter + help_subprob_counter f_counter = f_counter + help_f_counter stop_cond_counter = stop_cond_counter + help_stop_cond_counter subgrad1_counter = subgrad1_counter + help_subgrad1_counter subgrad2_counter = subgrad2_counter + help_subgrad2_counter clarke_f_counter = clarke_f_counter + help_clarke_f_counter clarke_sub_counter = clarke_sub_counter + help_clarke_sub_counter !_______________________________________________________________________________ !************************ STEP 1: END ****************************************** IF (reason_for_stop == 0) THEN ! in this case, a new iteration point is found in the 'main iteration' !_______________________________________________________________________________ !************************ STEP 2: BUNDLE UPDATE******************************* !---------------------------------------------------------- ! Update of the proximity parameter BEGINS !---------------------------------------------------------- u_k = 1.0_dp / t u = u_k u_min = 1.0_dp/t_max IF (extra_dec .AND. i_t>0 ) THEN u = (H1_new-H2_new-(H1_current-H2_current)) u = u / decrease u = 2.0_dp * u_k * (1-u) ELSE IF ( i_t > 3 ) THEN u = u_k * 0.5_dp END IF u_next = MAX(u, 0.1_dp * u_k ) u_next = MAX(u_next, u_min) i_t = MAX(i+1,1) IF (u_k /= u_next) THEN i_t = 1 END IF t_next = 1.0_dp / u_next IF (t_next >= t_max) THEN t = t_max ELSE IF (t_next < t_max .AND. t_next >= t_min) THEN t = t_next ELSE IF (t_next < t_min) THEN t = t_min END IF !---------------------------------------------------------- ! Update of the proximity parameter ENDS !---------------------------------------------------------- DO i = 1, number_of_func ! the values of components A_i and B_l for H_1 at x_new AandB_new_all(i) = AorB_i(fg1_new_all, fg2_new_all, fg1_new_all, fg2_new_all, i) END DO H1_new = H1(AandB_new_all) ! a new function value for H_1 H2_new = H2(fg2_new_all) ! a new function value for H_2 DO i = 1, number_of_obj ! the values of subgradients for DC components f1_i and f2_i at x_0 grad1 = subgradient_f1(x_new, problem1(i),factor(i)) grad2 = subgradient_f2(x_new, problem2(i),factor(i)) DO j = 1, user_n grad1_new_all(i,j) = grad1(j) grad2_new_all(i,j) = grad2(j) END DO END DO DO i = number_of_obj+1, number_of_func ! the values of subgradients for DC components g1_l and g2_l at x_0 grad1 = subgradient_g1(x_new, problem1(i),factor(i)) grad2 = subgradient_g2(x_new, problem2(i),factor(i)) DO j = 1, user_n grad1_new_all(i,j) = grad1(j) grad2_new_all(i,j) = grad2(j) END DO END DO subgrad1_counter = subgrad1_counter + 1 ! one subgradient calculated for each DC component f1_i and g1_l subgrad2_counter = subgrad2_counter + 1 ! one subgradient calculated for each DC component f2_i and g2_l ! Subgradients of H_1 and H_2 at the new point x_new ind = index_H1(AandB_new_all) grad1 = subgradient_AorB_i(grad1_new_all, grad2_new_all, ind) change1 = H1_new - H1_current grad2 = subgradient_H2(grad2_new_all) change2 = H2_new - H2_current dd = x_new - x_current ! the search direction change = H1_new - H2_new - (H1_current -H2_current) ! the change in the objective function value CALL update_b1(B1, grad1, ind, dd, AandB_new_all, AandB_k_all, & ! bundle update for B_1 & fg1_k_all, fg2_k_all, fg1_new_all, fg2_new_all) CALL update_b2(B2, grad2, dd, change2) ! bundle update for B_2 DO j = 1, number_of_func ! The rest of the new elements are added into the bundle B_1 IF ( j/= ind) THEN grad1 = subgradient_AorB_i(grad1_new_all, grad2_new_all, j) ! the subgradient of A_i and B_l at x_new alpha = 0.0_dp CALL add_element_b1(B1, grad1, alpha, j) END IF END DO x_current = x_new ! update of the current iteration point H1_current = H1_new ! update of the function value H_1 H2_current = H2_new ! update of the function value H_2 fg1_k_all = fg1_new_all ! update of the values f1_i and g1_l fg2_k_all = fg2_new_all ! update of the values f2_i and g2_l fg_k_all = fg1_k_all -fg2_k_all ! update of the values f_i and g_l AandB_k_all = AandB_new_all ! update of the component values A_i and B_l norm = SQRT(DOT_PRODUCT(dd,dd)) ! Distance between two consecutive iteration points !---------------------------------- ! Alteranative stopping condition utilizing the distance between consecutive iteration points ! IF ( change > - (10.0_dp)**(-6) ) THEN ! The difference between objective function values in consecutive steps is small enough ! IF ( norm < (10.0_dp)**(-6) ) THEN ! The difference between iteration points in consecutive steps is small enough ! stop_alg = .TRUE. ! reason_for_stop = 6 ! END IF !---------------------------------- IF ( (ABS(iprint) >= 3) .AND. (ABS(iprint) <= 5) ) THEN ! the basic/extended print of indermediate results IF (number_of_const == 0) THEN WRITE(*,*) 'Main Iter:', iter_counter, 'f(x):', fg_k_all, 't:', t ! WRITE(*,*) 'x:', x_current !WRITE(*,*) '------------------------------------------------------------------------' ELSE WRITE(*,*) 'Main Iter:', iter_counter, 'f(x) and g(x):', fg1_k_all - fg2_k_all , 't:', t ! WRITE(*,*) 'x:', x_current ! WRITE(*,*) '------------------------------------------------------------------------' END IF IF (iprint > 3) THEN ! IF iprint = 4 or 5 then the iteration point is printes WRITE(*,*) 'The new iteration point:' DO i = 1, user_n WRITE(*,*) 'x(', i ,')*=', x_current(i) END DO WRITE(*,*) ' ' END IF IF ( (ABS(iprint) == 4 ) .OR. (ABS(iprint) == 5 ) ) THEN WRITE(*,*) '--------------------------------------------------------------------------- ' WRITE(*,*) ' * * * * * DETAILIED INFORMATION ABOUT THE MAIN ITERATION', iter_counter ,' * * * * * ' WRITE(*,*) 'the number of rounds needed:', help_mit_iter_counter WRITE(*,*) 'the number of subproblems solved:', help_subprob_counter WRITE(*,*) 'the number of function values calculated for DC components:', help_f_counter WRITE(*,*) 'the number of subgradients calculated for DC components:', help_subgrad1_counter WRITE(*,*) 'the number of times Clarke stationary algorithm was tested:',help_stop_cond_counter IF (help_stop_cond_counter > 0) THEN WRITE(*,*) '--------------------------------------------------------------------------- ' WRITE(*,*) ' * * * * * * * * DETAILS ABOUT CLARKE STATIONARY ALGORITHM * * * * * * * * ' WRITE(*,*) '--------------------------------------------------------------------------- ' WRITE(*,*) 'The number of function values calculated for f', help_clarke_f_counter WRITE(*,*) 'The number of subgradients calculated for f', help_clarke_sub_counter WRITE(*,*) '--------------------------------------------------------------------------- ' END IF WRITE(*,*) '--------------------------------------------------------------------------- ' WRITE(*,*) 'CPU time used:', finish_time_main_it-start_time_main_it WRITE(*,*) '--------------------------------------------------------------------------- ' WRITE(*,*) ' ' END IF END IF !_______________________________________________________________________________ !************************ STEP 2: END **************************************** ELSE ! one of the stopping conditions is fulfilled stop_alg = .TRUE. ! the minimization algorithm can be STOPPED IF((ABS(iprint) >= 4) .AND. (ABS(iprint) <= 5)) THEN ! the extended print of indermediate results WRITE(*,*) ' ' IF ( (ABS(iprint) == 4 ) .OR. (ABS(iprint) == 5 ) ) THEN WRITE(*,*) '--------------------------------------------------------------------------- ' WRITE(*,*) ' * * * * * DETAILIED INFORMATION ABOUT THE MAIN ITERATION', iter_counter ,' * * * * * ' WRITE(*,*) 'the number of rounds needed:', help_mit_iter_counter WRITE(*,*) 'the number of subproblems solved:', help_subprob_counter WRITE(*,*) 'the number of function values calculated for DC components:', help_f_counter WRITE(*,*) 'the number of subgradients calculated for DC components:', help_subgrad1_counter WRITE(*,*) 'the number of times Clarke stationary algorithm was tested:',help_stop_cond_counter IF (help_stop_cond_counter > 0) THEN WRITE(*,*) '--------------------------------------------------------------------------- ' WRITE(*,*) ' * * * * * * * * DETAILS ABOUT CLARKE STATIONARY ALGORITHM * * * * * * * * ' WRITE(*,*) '--------------------------------------------------------------------------- ' WRITE(*,*) 'The number of function values calculated for f', help_clarke_f_counter WRITE(*,*) 'The number of subgradients calculated for f', help_clarke_sub_counter WRITE(*,*) '--------------------------------------------------------------------------- ' END IF WRITE(*,*) '--------------------------------------------------------------------------- ' WRITE(*,*) 'CPU time used:', finish_time_main_it-start_time_main_it WRITE(*,*) '--------------------------------------------------------------------------- ' END IF WRITE(*,*) ' ' WRITE(*,*) 'During Main iteration round:', iter_counter WRITE(*,*) 'Some of the stopping conditions is fulfilled and the algorithm is stopped.' WRITE(*,*) ' ' END IF END IF 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 = 4 END IF x_solution = x_current ! the solution to the minimization problem f_solution = H1_current - H2_current ! the objective function value at the solution factor_old = factor factor = 1.0_dp ! the values of different counters counter(1) = iter_counter counter(2) = subprob_counter counter(3) = f_counter counter(4) = subgrad1_counter counter(5) = subgrad2_counter counter(6) = stop_cond_counter counter(7) = clarke_f_counter counter(8) = clarke_sub_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 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(*,*) 'Objective functions:' DO i = 1, number_of_obj WRITE(*,*) 'f(x)', problem1(i), ':', f1(x_solution,problem1(i),factor(i))-f2(x_solution,problem2(i),factor(i)) END DO WRITE(*,*) '- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ' IF (number_of_const > 0) THEN WRITE(*,*) 'Constraint functions:' DO i = number_of_obj+1, number_of_func WRITE(*,*) 'g(x)', problem1(i), ':', g1(x_solution,problem1(i),factor(i))-g2(x_solution,problem2(i),factor(i)) END DO END IF 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(*,*) 'Objective functions:' DO i = 1, number_of_obj WRITE(*,*) 'f(x)', problem1(i), ':', f1(x_solution,problem1(i),factor(i))-f2(x_solution,problem2(i),factor(i)) END DO WRITE(*,*) '- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ' IF (number_of_const > 0) THEN WRITE(*,*) 'Constraint functions:' DO i = number_of_obj+1, number_of_func WRITE(*,*) 'g(x)', problem1(i), ':', g1(x_solution,problem1(i),factor(i))-g2(x_solution,problem2(i),factor(i)) END DO END IF WRITE(*,*) '---------------------------------------------------------------------------' WRITE(*,*) 'scaling:', factor_old WRITE(*,*) '---------------------------------------------------------------------------' WRITE(*,*) 'The number of main iterations:', iter_counter WRITE(*,*) 'The number of subproblems solved:', subprob_counter WRITE(*,*) 'The number of function values evaluated for DC components:', f_counter WRITE(*,*) 'The number of subgradients calculated for DC components:', subgrad1_counter WRITE(*,*) '---------------------------------------------------------------------------' WRITE(*,*) 'The number of times Clarke stationary algorithm (CSA) was used:', stop_cond_counter WRITE(*,*) 'The number of function values computed for f in CSA:', clarke_f_counter WRITE(*,*) 'The number of subgradients calculated for f in CSA: ', clarke_sub_counter WRITE(*,*) '---------------------------------------------------------------------------' WRITE(*,*) 'CPU time used:', finish_time-start_time WRITE(*,*) '--------------------------------------------------------------------------- ' WRITE(*,*) 'Elapsed time:', elapsed_time WRITE(*,*) '---------------------------------------------------------------------------' !WRITE(*,*) 'total n_f: ', f_counter+clarke_f_counter !WRITE(*,*) 'total n_xi ', subgrad1_counter+clarke_sub_counter !WRITE(*,*) '---------------------------------------------------------------------------' END IF WRITE(40,*) 'The values of the objectives at final point: ' DO i = 1, number_of_obj WRITE(40,*) 'f(x)', problem1(i), ':', f1(x_solution,problem1(i),factor(i))-f2(x_solution,problem2(i),factor(i)) END DO WRITE(40,*) '-----------------------------------------------------------------------------& &---------------------------' WRITE(40,142) subprob_counter 142 FORMAT(' The total number of subproblems solved: ',i12) WRITE(40,143) f_counter + clarke_f_counter 143 FORMAT(' The total number of the objective evaluations: ',i12) WRITE(40,144) subgrad1_counter + clarke_sub_counter 144 FORMAT(' The total number of gradient evaluations for DC components:',i12) WRITE(40,*) '-----------------------------------------------------------------------------& &---------------------------' WRITE(40,145) f_counter 145 FORMAT(' The number of the objective evaluations in main iterations: ',i12) WRITE(40,146) subgrad1_counter 146 FORMAT(' The number of gradient evaluations for DC components in main iterations:',i12) WRITE(40,*) '-----------------------------------------------------------------------------& &---------------------------' WRITE(40,147) clarke_f_counter 147 FORMAT(' The number of the objective evaluations in Clarke stationary algorithms: ',i12) WRITE(40,148) clarke_sub_counter 148 FORMAT(' The number of gradient evaluations for DC components in Clarke stationary algortihms:',i12) WRITE(40,*) '-----------------------------------------------------------------------------& &---------------------------' WRITE(40,149) finish_time-start_time 149 FORMAT(' The CPU time: ',f10.4) WRITE(40,150) elapsed_time 150 FORMAT(' The elapsed time: ',f10.4) CLOSE(40) END SUBROUTINE bundle_algorithm !....................................................................................... ! <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> ! _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ !| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | !** END ** END ** END ** END ** END ** END ** END ** END ** END ** END ** END ** END ** !|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_| ! _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ !| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | ! START ** START ** START ** START ** START ** START ** START ** START ** START ** START !|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_| ! <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> !*************************************************************************************** ! ---------------------------------------------------------------------------------- | ! | | | ! | THE MAIN ITERATION | | ! | | | ! ---------------------------------------------------------------------------------- | !*************************************************************************************** SUBROUTINE main_iteration( x_k, H1_k, H2_k, fg1_k_all, fg2_k_all,& & x_new, H1_new, H2_new, fg1_new_all, fg2_new_all, f_0_all,& & AB_k_all, m_extra, B1, B2, c1, c2, c3, t, t_min, t_max, i_t, & & decrease, crit_tol, eps, m, r_dec, r_inc, m_clarke, step_tol, & & mrounds, mrounds_clarke, iprint, reason_for_stop, iter_counter, & & subprob_counter, fi_counter, subgrad1_counter, subgrad2_counter,& & stop_cond_counter, clarke_f_counter, clarke_sub_counter, & & stepsize_used, extra_dec, factor, dec_tol) ! ! Executes the 'main iteration' algorithm. It is needed to find the search direction at the current iteration point. ! ! INPUT: * 'x_k' : The current iteration point ! * 'H1_k' and 'H2_k' : The values of the DC components H_1 and H_2 at the current iteration point ! * 'fg1_k_all' : The values of the DC components f1_i and g1_l at the current iteration point ! * 'fg2_k_all' : The values of the DC components f2_i and g2_l at the current iteration point ! * 'AB_k_all' : The values of the components A_i and B_l at the current iteration point ! * 'f_val' : The value used in descent condition ! * 'f_0_all' : The values of the objective functions f_i at the starting point ! * 'factor' : The values of sclaing factors of the objective functions f_i and constraints g_l ! ! * 'c1', 'c2' and 'c3' : The decrease parameters used in 'main iteration' ! * 'crit_tol' and 'eps': The stopping tolerance and the proximity measure ! * 'm' : The descent parameter ! * 'm_etxra' : The extra descent parameter ! * 'r_dec' and 'r_inc' : The decrease and increase parameters ! * 'm_clarke' : The descent parameter used in 'Clarke stationary' algorithm ! * 'step_tol' : The step-length parameter used in 'Clarke stationary' algorithm ! * 'dec_tol' : The tolerace for the approximate decrease obtained from the model ! ! * 'mrounds' : The maximum number of possible rounds in the 'main iteration' ! * 'mrounds_clarke' : The maximum number of possible rounds in the 'Clarke stationary' algorithm ! * 'iprint' : Specifies print ! ! * 'stepsize_used' : If .TRUE. then simple stepsize determination is used in the algorithm after each main iteration ! ! ! OUTPUT: * 'x_new' : the new iteration point ! * 'H1_new' and 'H2_new': the new value of DC components H_1 and H_2 ! * 'fg1_new_all' : The values of the DC components f1_i and g1_l at the new iteration point ! * 'fg2_new_all' : The values of the DC components f2_i and g2_l at the new iteration point ! * 'reason_for_stop' : Indicates the cause of the termination in the 'main iteration' ! * 'decrease ' : approximation of descent in objective H ! * 'extra_dec' : .TRUE. if the decrease in the objective satisfies the extra condtion ! ! * 'iter_counter' : the number of rounds needed in 'main iteration' ! * 'subprob_counter' : the number of subproblems solved in 'main iteration' ! * 'fi_counter' : the number of function values evaluated for a DC component in 'main iteration' (same for f_1 and f_2) ! * 'subgrad1_counter' : the number of subgradients calculated for f_1 in 'main iteration' ! * 'subgrad2_counter' : the number of subgradients calculated for f_2 in 'main iteration' ! * 'stop_cond_counter' : the number of times approximate stopping condition is tested during the 'main iteration' ! * 'clarke_f_counter' : the number of function values evaluated for f in 'Clarke stationary' algorithm ! * 'clarke_sub_counter' : the number of subgradients evaluated for f in 'Clarke stationary' algorithm ! ! INOUT: * 'B_1' and B_2' : the bundles of the DC components f_1 and f_2 ! * 't' : the value of the proximity parameter ! * 't_min' and 't_max' : the lower and upper bounds of the proximity parameter 't' ! * 'i_t' : the index used in the update procedure of the proximity parameter 't' ! ! NOTICE: * The dimensions of vectors 'x_k' and 'x_new' has to be same. (i.e. the dimensio of 'x_k' and 'x_new' has to be ! 'user_n' when SUBROUTINE main_iteration() is used in SUBROUTINE bundle_method. ! * The dimension of vectors 'fg1_k_all', 'fg2_k_all', 'fg1_new_all', 'fg2_new_all' and AB_k_all' has to be 'number_of_func' ! * The dimention of vector 'f_0_all' has to be 'number_of_obj' ! !*********************************************************************************** IMPLICIT NONE !**************************** NEEDED FROM THE USER ********************************* TYPE(kimppu1), INTENT(INOUT) :: B1 ! the bundle B_1 for the DC component f_1 TYPE(kimppu2), INTENT(INOUT) :: B2 ! the bundle B_2 for the DC component f_2 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: x_k ! the current iteration point REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: fg1_k_all ! the values of DC components f1_i at the current iteration point REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: fg2_k_all ! the values of DC components f2_i at the current iteration point REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: x_new ! the new iteration point if 'reason_for_stop'=0 in the 'main iteration' REAL(KIND=dp), INTENT(IN) :: H1_k, H2_k ! the value of H_1 and H_2 at x_k REAL(KIND=dp), INTENT(OUT) :: H1_new, H2_new ! the value of H_1 and H_2 at x_new if 'reason_for_stop=0' in the 'main iteration' REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: fg1_new_all ! the values of DC components f2_i at the new iteration point REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: fg2_new_all ! the values of DC components f2_i at the new iteration point REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: f_0_all ! the values of DC components f_i at the starting point x_0 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: AB_k_all ! the values of components A_i at the starting point x_k REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: factor ! the values of scaling factors of f_i and g_l REAL(KIND=dp), INTENT(IN) :: c1, c2, c3 ! the decrease parameters used in main iteration REAL(KIND=dp), INTENT(IN) :: crit_tol, eps ! crit_tol=stopping tolerance and eps=proximity measure (needed in stopping conditions) REAL(KIND=dp), INTENT(IN) :: m ! descent parameter REAL(KIND=dp), INTENT(IN) :: m_extra ! extra descent parameter REAL(KIND=dp), INTENT(IN) :: r_dec, r_inc ! r_dec=decrease parameter and r_inc=increase parameter REAL(KIND=dp), INTENT(IN) :: m_clarke ! The descent parameter used in 'Clarke stationary' algorithm REAL(KIND=dp), INTENT(IN) :: step_tol ! The step-length parameter used in 'Clarke stationary' algorithm REAL(KIND=dp), INTENT(IN) :: dec_tol ! The tolerace for the approximate decrease obtained from the model REAL(KIND=dp), INTENT(OUT) :: decrease ! approximation of descent in objective H REAL(KIND=dp), INTENT(INOUT) :: t ! the proximity parameter t REAL(KIND=dp), INTENT(INOUT) :: t_min, t_max ! the bounds for proximity parameter t INTEGER, INTENT(INOUT) :: i_t ! the index used in update proxcedure of proximity parameter t INTEGER :: mrounds ! the maximum number of possible rounds in the 'main iteration' ! If mrounds<=0, then DEFAULT value 500 is used (This is also done in the SUBROUTINE bundle_method() INTEGER :: mrounds_clarke ! the maximum number of possible rounds in the 'Clarke statinonary' algorithm INTEGER, INTENT(OUT) :: reason_for_stop ! 0 - a new iteration point found ! 1 - stopping condition satisfied (i.e. approximate Clarke stationarity) ! 2 - approximate stopping condition satisfied (i.e. the step-length beta* < beta_0) ! 3 - maximum number of rounds executed in 'main iteration' ! 5 - maximum number of rounds executed in 'Clarke stationary' algorithm INTEGER, INTENT(OUT) :: iter_counter ! the number of rounds used in 'main iteration' INTEGER, INTENT(OUT) :: subprob_counter ! the number of subproblems solved in 'main iteration' INTEGER, INTENT(OUT) :: fi_counter ! the number of function values evaluated for a DC component in 'main iteration' (same for f_1 and f_2) INTEGER, INTENT(OUT) :: subgrad1_counter ! the number of subgradients calculated for f_1 in 'main iteration' INTEGER, INTENT(OUT) :: subgrad2_counter ! the number of subgradients calculated for f_2 in 'main iteration' INTEGER, INTENT(OUT) :: stop_cond_counter ! the number of times approximate stopping condition is tested during the 'main iteration' INTEGER, INTENT(OUT) :: clarke_f_counter ! the number of function values evaluated for f in 'Clarke stationary' algorithm INTEGER, INTENT(OUT) :: clarke_sub_counter ! the number of subgradients evaluated for f in 'Clarke stationary' algorithm INTEGER, INTENT(IN) :: iprint ! the variable that specifies print: ! iprint = 0 : print is suppressed ! iprint = 1 : basic print of final result (nothing is printed in main iteration) ! iprint = -1: basic print of final result without the solution vector (nothing is printed in main iteration) ! iprint = 2 : extended print of final result (nothing is printed in main iteration) ! iprint = -2: extended print of final result without the solution vector (nothing is printed in main iteration) ! iprint = 3 : basic print of intermediate results and extended print of final results (nothing is printed in main iteration) ! iprint = -3: basic print of intermediate results and extended print of final results without the solution vector (nothing is printed in main iteration) ! iprint = 4 : extended print of intermediate and final results (something is printed about the 'main iteration' (this print is done in the bundle algorithm)) ! iprint = -4: extended print of intermediate and final results without the solution vector (something is printed about the 'main iteration' (this print is done in the bundle algorithm)) ! iprint = 5 : everything is printed step by step (this is the only option which causes print in the 'main iteration') LOGICAL, INTENT(OUT) :: extra_dec ! If .TRUE. then the descent satisfies an extra condition LOGICAL, INTENT(IN) :: stepsize_used ! If .TRUE. then simple stepsize determination is used in the algorithm after each main iteration ! !***************************** LOCAL VARIABLES ************************************ REAL(KIND=dp), DIMENSION(number_of_func, user_n) :: grad1_y_all ! subgradients for DC components f1_i and g1_l at y REAL(KIND=dp), DIMENSION(number_of_func, user_n) :: grad2_y_all ! subgradients for DC components f2_i and g2_l at y REAL(KIND=dp), DIMENSION(number_of_func) :: fg1_y_all ! values of DC components f1_i and g1_l at y REAL(KIND=dp), DIMENSION(number_of_func) :: fg2_y_all ! values of DC components f2_i and g2_l at y REAL(KIND=dp), DIMENSION(number_of_func) :: AB_y_all ! values of components A_i at y REAL(KIND=dp), DIMENSION(number_of_func) :: fg1_koe_all ! values of DC components f1_i and g1_l at x_koe REAL(KIND=dp), DIMENSION(number_of_func) :: fg2_koe_all ! values of DC components f2_i and g2_l at x_koe REAL(KIND=dp), DIMENSION(number_of_func) :: AB_koe_all ! values of components A_i and B_l at x_koe REAL(KIND=dp), DIMENSION(give_n_b1(B1)) :: d_t ! the search direction REAL(KIND=dp), DIMENSION(give_n_b1(B1)) :: y ! the new auxilary point REAL(KIND=dp), DIMENSION(give_n_b1(B1)) :: new_grad1, new_grad2 ! the subgradients of A_i or B_l and H_2 at y REAL(KIND=dp), DIMENSION(give_n_b1(B1)) :: vect ! 'help' vector REAL(KIND=dp) :: norm1 ! ||\bxi_1(x_k)|| REAL(KIND=dp) :: max_norm ! the value of the maximum subgradient norm ||\bxi_{2,max}|| REAL(KIND=dp) :: d_norm ! the norm of the direction vector d_t REAL(KIND=dp) :: delta1, delta2 ! the predicted changes of H_1 and H_2, respectively REAL(KIND=dp) :: H1_y, H2_y ! the function values of H_1 and H_2 at y REAL(KIND=dp) :: H_k ! the function value of H at x_k REAL(KIND=dp) :: real_decrease ! the real decrease of the objective function H REAL(KIND=dp) :: lin_err1, lin_err2 ! the linearization errors of A_i or B_l and H_2 (calculated using y) REAL(KIND=dp) :: erotus ! the biggest difference obtained when function values at x_0 and at the new iteration point are compared (each f_i and g_l are taken into account) REAL(KIND=dp) :: apu_erotus ! the help variable used to obtain 'erotus' REAL(KIND=dp) :: help ! 'help' variable INTEGER :: f_counter ! help counter for the number of function values INTEGER :: sub_counter ! help counter for the number of subgradients INTEGER :: clarke_iter_counter ! help counter for the number of subgradients INTEGER :: s_counter ! the number of subproblems solved during the current iteration round of the 'main iteration' algorithm INTEGER :: i, j ! 'help' variables INTEGER :: i_dec ! the number of bundle update steps made successively LOGICAL :: was_b1_full ! .TRUE. if during the previous round something was added to B_1 and B_1 was full before this insertion LOGICAL :: stop_main_it ! .TRUE. if the current 'main iteration' can be stopped LOGICAL :: stop_askelpituus ! .TRUE. if the current 'line_search' can be stopped REAL(KIND=dp) :: askelpituus ! stepsize determined during 'line search' into descent direction REAL(KIND=dp) :: koe_H1 ! the value of f1 at the point tested in 'line search' REAL(KIND=dp) :: koe_H2 ! the value of f2 at the point tested in 'line search' REAL(KIND=dp), DIMENSION(give_n_b1(B1)) :: koe_y ! a point tested during 'line search' !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! <> <> <> <> <> <> <> MAIN ITERATION STARTS <> <> <> <> <> <> <> <> !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ !__________________________________________________________________________________ !************************ STEP 0: STOPPING CONDITION ****************************** iter_counter = 1 ! the first round is executed subprob_counter = 0 ! the number of subproblems solved fi_counter = 0 ! the number of function values evaluated subgrad1_counter = 0 ! the number of subgradients calculated for f_1 subgrad2_counter = 0 ! the number of subgradients calculated for f_2 stop_cond_counter = 0 ! the number of times approximate stopping condition is tested clarke_f_counter = 0 ! the number of function values evaluated for f in 'Clarke stationary' algorithm clarke_sub_counter = 0 ! the number of subgradients evaluated for f in 'Clarke stationary' algorithm i_dec = 0 ! the number of bundle update steps made successively stop_main_it = .FALSE. ! We are NOT rady to stop vect = give_subgrad_b1(B1,0) - give_subgrad_b2(B2,0) ! vector: ' \bxi_1(x_k) - \bxi_2(x_k) ' start: IF ( SQRT(DOT_PRODUCT(vect,vect)) < (crit_tol)) THEN ! ||\bxi_1(x_k) - \bxi_2(x_k)|| < crit_tol !**>>**>>**>> APPROXIMATE CRITICALITY OBTAINED <<**<<**<<** d_t = 1.0_dp CALL guaranteeing_clarke( x_k, H1_k, H2_k, fg1_k_all, fg2_k_all, & & reason_for_stop, x_new , H1_new, H2_new, fg1_new_all, fg2_new_all,& & d_t, crit_tol, m_clarke, step_tol, iprint, mrounds_clarke, & & clarke_iter_counter, f_counter, sub_counter, factor ) clarke_f_counter = clarke_f_counter + f_counter clarke_sub_counter = clarke_sub_counter + sub_counter stop_cond_counter = stop_cond_counter + 1 ! The times that 'Clarke stationary' algorithm has been executed is increased by one stop_main_it = .TRUE. ! the 'main iteration' can be stopped extra_dec = .FALSE. i_t = 0 IF(iprint == 5) THEN ! print everything WRITE(*,*) 'At STEP 0: the Clarke stationary algorithm is executed' SELECT CASE(reason_for_Stop) CASE(0) WRITE(*,*) ' ' WRITE(*,*) 'A BETTER iteration point is found.' WRITE(*,*) ' ' CASE(1) WRITE(*,*) ' ' WRITE(*,*) 'The solution satisfies the CLARKE STATIONARITY CONDITION.' WRITE(*,*) ' ' CASE(2) WRITE(*,*) ' ' WRITE(*,*) 'The solution satisfies the APPROXIMATE STOPPING CONDITION.' WRITE(*,*) ' ' CASE(5) WRITE(*,*) ' ' WRITE(*,*) 'The maximum number of rounds is executed in Clarke stationary algorithm.' WRITE(*,*) ' ' END SELECT END IF !__________________________________________________________________________________ !************************ STEP 0: END ********************************************* ELSE start !_______________________________________________________________________________ !************************ STEP 1: PARAMETER INITIALIZATION ********************* H_k = H1_k - H2_k ! the function value of H at x_k vect = give_subgrad_b1(B1,0) ! the vector \bxi_1(x_k) norm1 = SQRT(DOT_PRODUCT(vect,vect)) ! the value of the norm ||\bxi_1(x_k)|| max_norm = max_norm_value(B2) ! the value of the maximum subgradient norm in the bundle B_2 t_min = (0.5_dp * r_dec * eps) / ( norm1 + max_norm ) ! the lower bound for t t_max = r_inc * t_min ! the upper bound for t ! Checking that the parameter t is from the right interval IF (t < t_min) THEN t = t_min END IF IF (t > t_max) THEN t = t_max END IF ! print everything IF(iprint == 5) THEN WRITE(*,*) ' ' WRITE(*,*) 'STEP 1: INITIALIZATION OF PARAMETERS' WRITE(*,*) ' ' WRITE(*,*) 't=', t WRITE(*,*) 't is from the interval [t_min, t_max]=[', t_min, ',', t_max, ']' WRITE(*,*) 'eps=', eps WRITE(*,*) 'The maximum subgradient norm in B_2: ', max_norm, ' The norm |||\bxi_1(x_k)||: ', norm1 WRITE(*,*) ' ' WRITE(*,*) '-------------------------------------------------------------------------------' END IF was_b1_full = .FALSE. ! B_1 was not full during the previous round since this is the first round of 'main iteration' IF ( mrounds <= 0 ) THEN mrounds = 5000 ! the DEFAULT value for 'mrounds' is used END IF !_______________________________________________________________________________ !************************ STEP 1: END ****************************************** END IF start !---------------------------------------------------------------------------------- ! <> <> <> <> <> <> <> REPEATABLE PART BEGINS <> <> <> <> <> <> <> <> !---------------------------------------------------------------------------------- DO WHILE ( (.NOT. stop_main_it ) .AND. (iter_counter <= mrounds )) ! is repeated until the 'main iteration' can be terminated !______________________________________________________________________________ !************************ STEP 2: SEARCH DIRECTION **************************** CALL subproblem_solver(x_k, AB_k_all, give_n_b1(B1),& ! subproblems are solved & give_size_b1(B1)+1, B1, B2, t, s_counter) subprob_counter = subprob_counter + s_counter ! the number of subproblems solved so far CALL add_glob_index(B2) ! calculates the index of the subproblem yielding the global solution d_t = give_solution(B2) ! the global solution d_t is selected d_norm = SQRT(DOT_PRODUCT(d_t,d_t)) ! the norm ||d_t|| of the solution d_t is calculated ! print everything IF(iprint == 5) THEN WRITE(*,*) ' ' WRITE(*,*) 'STEP 2: SEARCH DIRECTION' WRITE(*,*) ' ' WRITE(*,*) 'Round:' , iter_counter , 'paramter t = ', t WRITE(*,*) 'Search direction: [', d_t , ']' WRITE(*,*) 'The norm of search direction: ', d_norm WRITE(*,*) 'The number of subproblems solved during this round: ', s_counter WRITE(*,*) 'The number of subproblems solved till now: ', subprob_counter WRITE(*,*) ' ' WRITE(*,*) '-------------------------------------------------------------------------------' END IF !______________________________________________________________________________ !************************ STEP 2: END ***************************************** !->->->->->-> EXECUTION OF BRANCH BEGINS (2 POSSIBLE BRANCHES) <-<-<-<-<-<-<-<- branches: IF (d_norm < crit_tol .OR. (give_decrease(B2)-H1_k)>dec_tol ) THEN !->->->->->->->->->->->->->-> BRANCH 1 BEGIN <-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<- !__________________________________________________________________________ !******************** STEP 3: CLARKE STATIONARITY CHECK ******************* CALL guaranteeing_clarke( x_k, H1_k, H2_k, fg1_k_all, fg2_k_all,& & reason_for_stop, x_new, H1_new, H2_new, fg1_new_all, fg2_new_all,& & d_t, crit_tol, m_clarke, step_tol, iprint, mrounds_clarke, & & clarke_iter_counter, f_counter, sub_counter, factor ) !Counters are updated clarke_f_counter = clarke_f_counter + f_counter clarke_sub_counter = clarke_sub_counter + sub_counter stop_cond_counter = stop_cond_counter + 1 ! the update of the stopping condition counter stop_main_it = .TRUE. extra_dec = .FALSE. i_t = 0 ! print everything IF(iprint == 5) THEN WRITE(*,*) ' ' WRITE(*,*) 'STEP 3: Clarke stationarity test' WRITE(*,*) ' ' WRITE(*,*) 'The number of the Clarke stationary algorithm' , stop_cond_counter SELECT CASE(reason_for_Stop) CASE(0) WRITE(*,*) ' ' WRITE(*,*) 'A BETTER iteration point is found.' WRITE(*,*) ' ' CASE(1) WRITE(*,*) ' ' WRITE(*,*) 'The solution satisfies the CLARKE STATIONARITY CONDITION.' WRITE(*,*) ' ' CASE(2) WRITE(*,*) ' ' WRITE(*,*) 'The solution satisfies the APPROXIMATE STOPPING CONDITION.' WRITE(*,*) ' ' CASE(5) WRITE(*,*) ' ' WRITE(*,*) 'The maximum number of rounds is executed in algorithm.' WRITE(*,*) ' ' END SELECT END IF !__________________________________________________________________________ !******************** STEP 3: APPROXIMATE STOPPING CONDITION ************** !->->->->->->->->->->->->->-> BRANCH 1 END <-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<- ELSE branches !->->->->->->->->->->->->->-> BRANCH 2 BEGIN <-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<- !__________________________________________________________________________ !******************** STEP 4: DESCENT TEST ******************************** y = x_k + d_t ! a new auxilary point y DO i = 1, number_of_obj ! the values of DC components f1_i and f2_i at y fg1_y_all(i) = f1(y,problem1(i),factor(i)) fg2_y_all(i) = f2(y,problem2(i),factor(i)) END DO DO i = number_of_obj+1, number_of_func ! the values of DC components g1_l and g2_l at y fg1_y_all(i) = g1(y,problem1(i),factor(i)) fg2_y_all(i) = g2(y,problem2(i),factor(i)) END DO fi_counter = fi_counter + 1 ! one more objective function value calculated for a DC component DO i = 1, number_of_func ! the values of components A_i and B_l for H_1 at y AB_y_all(i) = AorB_i(fg1_y_all, fg2_y_all, fg1_k_all, fg2_k_all ,i) END DO H1_y = H1(AB_y_all) ! the value of H_1 at y H2_y = H2(fg2_y_all) ! the value of H_2 at y real_decrease = ( H1_y - H2_y ) - H_k ! the real decrease in the value of H ! print everything IF(iprint == 5) THEN WRITE(*,*) ' ' WRITE(*,*) 'STEP 4: DESCENT TEST' WRITE(*,*) ' ' WRITE(*,*) 'A new AUXILIARY point: [', y, ']' WRITE(*,*) ' ' WRITE(*,*) 'Real decrease:' , real_decrease, ' Predicted decrease: ', give_decrease(B2) WRITE(*,*) ' Is DECREASE enough: ', real_decrease, '<=', (m *( give_decrease(B2)-H1_k)), '?' WRITE(*,*) ' ' END IF branch2: IF ( real_decrease <= (m * (give_decrease(B2)-H1_k)))THEN ! is the descent in the function H enough? !**>>**>>**>> NEW ITERATION POINT FOUND <<**<<**<<** IF(stepsize_used) THEN !-**--**- STEPSIZE DETERMINATION -**--**- IF (real_decrease < -0.1_dp) THEN askelpituus = 1.0_dp stop_askelpituus = .FALSE. DO WHILE((.NOT. stop_askelpituus)) koe_y = y + d_t DO i = 1, number_of_obj ! the values of DC components f1_i and f2_i at koe_y fg1_koe_all(i) = f1(koe_y,problem1(i),factor(i)) fg2_koe_all(i) = f2(koe_y,problem2(i),factor(i)) END DO DO i = number_of_obj+1, number_of_func ! the values of DC components g1_l and g2_l at koe_y fg1_koe_all(i) = g1(koe_y,problem1(i),factor(i)) fg2_koe_all(i) = g2(koe_y,problem2(i),factor(i)) END DO fi_counter = fi_counter + 1 ! one more objective function value calculated for a DC component DO i = 1, number_of_func ! the values of components A_i and B_l for H_1 at y AB_koe_all(i) = AorB_i(fg1_koe_all, fg2_koe_all, fg1_k_all, fg2_k_all ,i) END DO koe_H1 = H1(AB_koe_all) ! the value of H_1 at y koe_H2 = H2(fg2_koe_all) ! the value of H_2 at y IF ( (koe_H1 - koe_H2 - H_k) <= ( 0.7_dp * m * (give_decrease(B2)-H1_k))) THEN !the case where the stepsize can be increased askelpituus = askelpituus + 1.0_dp y = koe_y H1_y = koe_H1 H2_y = koe_H2 fg1_y_all = fg1_koe_all fg2_y_all = fg2_koe_all ELSE stop_askelpituus = .TRUE. END IF END DO END IF !-**--**- STEPSIZE DETERMINATION END -**--**- END IF x_new = y ! the new iteration point is the current auxilary point H1_new = H1_y ! the value of H_1 at x_new H2_new = H2_y ! the value of H_2 at x_new fg1_new_all = fg1_y_all ! the values of DC components f1_i and g1_l at x_new fg2_new_all = fg2_y_all ! the values of DC components f2_i and g2_l at x_new stop_main_it = .TRUE. ! the 'main iteration' can be stopped reason_for_stop = 0 ! the reason for stopping is the new iteration point decrease = give_decrease(B2)-H1_k IF ( real_decrease <= ( m_extra * decrease )) THEN extra_dec = .TRUE. ELSE extra_dec = .FALSE. END IF ! print everything IF(iprint == 5) THEN WRITE(*,*) 'A NEW ITERATION POINT is found: x_new= [', x_new, ']' WRITE(*,*) 'A NEW OBJECTIVE VALUE: f_new=', H1_new - H2_new WRITE(*,*) 'The old point: [', x_k,']' WRITE(*,*) 'The old objective value:', H_k WRITE(*,*) ' ' END IF !__________________________________________________________________________ !******************** STEP 4: END ***************************************** ELSE branch2 !----------- SOME UPDATES IN VARIABLES BEGINS ------------------------- ! print everything IF (iprint == 5) THEN WRITE(*,*) '-----------------------------------------------------------------------------------' END IF !__________________________________________________________________________ !******************** STEP 5: BUNDLE UPDATE ******************************* erotus = fg1_y_all(1)-fg2_y_all(1)-f_0_all(1) !the biggest difference between the new auxiliary point y and x_0 when we compare function values of f_i and g_l DO i =2, number_of_obj apu_erotus = fg1_y_all(i)-fg2_y_all(i)-f_0_all(i) IF (apu_erotus > erotus) THEN erotus = apu_erotus END IF END DO update: IF (( erotus>0 ) .AND. & & (d_norm > eps)) THEN !----------- PARAMETER t REDUCTION -------------------------------- t = t - c1 * ( t - t_min ) ! the parameter t is updated i_dec = 0 ! the numebr of bundle update steps made successively i_t = -1 ! null step is done in main iteration ! print everything IF(iprint == 5) THEN WRITE(*,*) ' ' WRITE(*,*) 'STEP 5a: PARAMETER t UPDATE' WRITE(*,*) ' ' WRITE(*,*) 'Updated parameter t=', t WRITE(*,*) '-----------------------------------------------------------------------------------' END IF ELSE update !----------- BUNDLE AND PARAMETER UPDATE -------------------------- DO i = 1, number_of_obj ! the values of subgradients for DC components f1_i and f2_i at y new_grad1 = subgradient_f1(y, problem1(i),factor(i)) new_grad2 = subgradient_f2(y, problem2(i),factor(i)) DO j = 1, user_n grad1_y_all(i,j) = new_grad1(j) grad2_y_all(i,j) = new_grad2(j) END DO END DO DO i = number_of_obj+1, number_of_func ! the values of subgradients for DC components g1_l and g2_l at y new_grad1 = subgradient_g1(y, problem1(i),factor(i)) new_grad2 = subgradient_g2(y, problem2(i),factor(i)) DO j = 1, user_n grad1_y_all(i,j) = new_grad1(j) grad2_y_all(i,j) = new_grad2(j) END DO END DO subgrad1_counter = subgrad1_counter + 1 ! the subgradient counter is updated subgrad2_counter = subgrad2_counter + 1 ! the subgradient counter is updated DO j = 1, number_of_func new_grad1 = subgradient_AorB_i(grad1_y_all, grad2_y_all, j) ! a subgradient of A_i or B_l at y lin_err1 = AB_k_all(j) - AB_y_all(j) + & & DOT_PRODUCT(d_t,new_grad1) ! a new linearization error for A_i or B_l !->>>>>>>>>>>> BUNDLE B_1 UPDATE BEGINS <<<<<<<<<<<<<<<<<<<<<<<<<<- bundle1: IF (is_full_b1(B1)) THEN ! bundle B_1 is full and a new element is added to the bundle B_1 CALL add_element_b1(B1, new_grad1, lin_err1, j) was_b1_full = .TRUE. ELSE bundle1 ! bundle B_1 is NOT full and a new element is added to the bundle B_1 CALL add_element_b1(B1, new_grad1, lin_err1, j) was_b1_full = .FALSE. END IF bundle1 !->>>>>>>>>>>> BUNDLE B_1 UPDATE ENDS <<<<<<<<<<<<<<<<<<<<<<<<<<<- END DO ! print everything IF(iprint == 5) THEN WRITE(*,*) ' ' WRITE(*,*) 'STEP 5b: BUNDLE UPDATE' WRITE(*,*) ' ' WRITE(*,*) 'delta2: ', delta2, ' delta1: ', delta1, ' glob solution index in B2: ', i WRITE(*,*) ' ' WRITE(*,*) 'A new bundle element is inserted into the bundle B_1. ' END IF !->>>>>>>>>>>> BUNDLE B_2 AND PARAMETER UPDATE BEGINS <<<<<<<<<<<- !bundle2: IF( delta2 >= 0 .AND. (give_max_size_b2(B2)>0) ) THEN ! bundle B_2 is updated (happens only when delta2 >= 0 and bundle size is larger than 1) (SOME OTHER INSERTION RULES ALSO POSSIBLE!) bundle2: IF( (give_max_size_b2(B2)>0) ) THEN ! bundle B_2 is updated (happens only when delta2 >= 0 and bundle size is larger than 1) (SOME OTHER INSERTION RULES ALSO POSSIBLE!) new_grad2 = subgradient_H2(grad2_y_all) ! a subgradient of f_2 at y lin_err2 = H2_k - H2_y + DOT_PRODUCT(d_t,new_grad2) ! a new linearization error for f_2 CALL add_element_b2(B2, new_grad2, lin_err2) ! a new element is inserted into the bundle B_2 !In the algorithm we never overwrite the 'bundle element' yielding the previous global solution of the search direction problem. ! print everything IF(iprint == 5) THEN WRITE(*,*) 'A new bundle element is inserted into the bundle B_2' WRITE(*,*) ' ' WRITE(*,*) ' subgrad2 = [', new_grad2,']' WRITE(*,*) ' ' WRITE(*,*) '--------------------& &------------------------------& &--------------------------------' END IF !__________________________________________________________________________ !******************** STEP 5: END ***************************************** !__________________________________________________________________________ !******************** STEP 6: PARAMETER UPDATE **************************** help = SQRT(DOT_PRODUCT(new_grad2,new_grad2)) ! the norm of the new subgradient of f_2 (subgradient is calculated ! at the new auxilary point y) IF(help > max_norm) THEN max_norm = help ! the updated value of the maximum subgradient norm in the bundle B_2 t_min = (0.5_dp * r_dec * eps) / ( norm1 + max_norm ) ! the updated value of the lower bound t_min END IF ! print everything IF(iprint == 5) THEN WRITE(*,*) '---------------------& &----------------------------------------------------& &----------' WRITE(*,*) ' ' WRITE(*,*) 'STEP 6: PARAMETER UPDATE' WRITE(*,*) ' ' END IF !__________________________________________________________________________ !******************** STEP 6: END ***************************************** END IF bundle2 ! --> the update of the proximity parameter BEGINS <-- i_t = -1 ! Null step is done in main iteration IF (i_dec <= 50) THEN t = t - c2 * ( t - t_min ) ! the parameter t is decreased ELSE IF (i_dec >= 51 ) THEN t = t - c3 * ( t - t_min ) ! the parameter t is decreased !WRITE(*,*) 'taalla' END IF i_dec = i_dec + 1 ! one more bundle update step executed ! --> the update of the proximity parameter ENDS <-- ! print everything IF(iprint == 5) THEN WRITE(*,*) '-----------------------------------------------------------------------------------' END IF !->>>>>>>>>>>> BUNDLE B_2 AND PARAMETER UPDATE ENDS <<<<<<<<<<<<<- END IF update !------ PARAMETER t REDUCTION & BUNDLE AND PARAMETER UPDATE ENDS ----- iter_counter = iter_counter + 1 ! update of the iteration counter END IF branch2 !--------- NEW ITERATION POINT & SOME UPDATES IN VARIABLES ENDS --------- !->->->->->->->->->->->->->-> BRANCH 2 END <-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<- END IF branches !->->->->->->-> EXECUTION OF BRANCH ENDS (2 BRANCHES) <-<-<-<-<-<-<-<-<-<-<-<- END DO !---------------------------------------------------------------------------------- ! <> <> <> <> <> <> <> REPEATABLE PART ENDS <> <> <> <> <> <> <> <> !---------------------------------------------------------------------------------- IF ( (.NOT. stop_main_it ) ) THEN reason_for_stop = 3 ! the maximum number of rounds have been executed iter_counter = iter_counter -1 END IF END SUBROUTINE main_iteration !....................................................................................... ! <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> ! _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ !| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | !** END ** END ** END ** END ** END ** END ** END ** END ** END ** END ** END ** END ** !|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_| ! _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ !| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | ! START ** START ** START ** START ** START ** START ** START ** START ** START ** START !|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_| ! <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> !*************************************************************************************** ! ---------------------------------------------------------------------------------- | ! | | | ! | Guaranteeing Clarke stationarity | | ! | | | ! ---------------------------------------------------------------------------------- | !*************************************************************************************** SUBROUTINE guaranteeing_clarke( x_k, H1_k, H2_k, fg1_k_all, fg2_k_all, & & reason_for_stop, x_new , H1_new, H2_new, fg1_new_all, fg2_new_all,& & d_cause, crit_tol, m, step_tol, iprint, mrounds, & & iter_counter, f_counter, subgrad_counter, factor ) ! ! Executes the 'Clarke stationary' algorithm. It is needed to guarantee Clarke stationarity at the current iteration point. ! If the current point is not Clarke stationary then the algorithm generates a descent direction. ! ! INPUT: * 'x_k' : The current iteratioin point ! * 'H1_k' and 'H2_k' : The values of DC components H_1 and H_2 at the current iteration point ! ! * 'fg1_k_all' and 'fg2_k_all' : The values of DC components f1_i, f2_i, g1_l and g2_l at the current iteration point ! ! * 'factor' : The values of scaling factors for f_i and g_l ! * 'd_cause' : The search direction which caused the execution of the clarke stationary algorithm ! * 'crit_tol' : The stopping tolerance ! * 'm' : The descent parameter ! * 'step_tol' : The step-length tolerance ! ! * 'mrounds' : The maximum number of possible rounds in the 'Clarke stationary' algorithm ! * 'iprint' : Specifies print ! ! ! OUTPUT: * 'x_new' : Tthe new iteration point (if it is found, i.e., 'reason_for_stop' = 0) ! * 'H1_new' and 'H2_new' : The values of DC components H_1 and H_2 at the new iteration point ! * 'fg1_new_all' and 'fg2_new_all': the new value of DC components f1_i, f2_i, g1_l and g2_l (if new iteration point is found, i.e., 'reason_for_stop' = 0) ! ! * 'reason_for_stop' : Indicates the cause of the termination in the 'Clarke stationarity' algorithm ! ! * 'iter_counter' : the number of rounds needed in 'Clarke stationary' algorithm ! * 'f_counter' : the number of function values evaluated for f in 'Clarke stationary' algorithm (same for f_1 and f_2) ! * 'subgrad_counter' : the number of subgradients calculated for f in 'Clarke stationary' algorithm (same for f_1 and f_2) ! ! NOTICE: * The dimensions of vectors 'x_k', 'x_new' and 'd_cause' has to be same. (i.e. the dimension of 'x_k', 'x_new' and 'd_cause' has to be ! 'user_n' when SUBROUTINE guaranteeing_clarke() is used in SUBROUTINE main_iteration(). ! * The dimension of vectors 'fg1_k_all', 'fg2_k_all', 'fg1_new_all' and 'fg2_new_all' has to be 'number_of_func' ! !*********************************************************************************** IMPLICIT NONE !**************************** NEEDED FROM THE USER ********************************* REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: x_k ! the current iteration point REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: d_cause ! the search direction REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: x_new ! the new iteration point if 'reason_for_stop'=0 in the 'Clarke stationary' algorithm REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: fg1_k_all ! the values of DC components f1_i and g1_l at x_k REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: fg2_k_all ! the values of DC components f2_i and g2_l at x_k REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: fg1_new_all ! the values of DC components f1_i and g1_l at x_new REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: fg2_new_all ! the values of DC components f2_i and g2_l at x_new REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: factor ! the values of scaling factors for f_i and g_l REAL(KIND=dp), INTENT(IN) :: H1_k, H2_k ! the value of H_1 and H_2 at x_k REAL(KIND=dp), INTENT(OUT) :: H1_new, H2_new ! the value of H_1 and H_2 at x_new if 'reason_for_stop=0' in the 'Clarke stationary' algorithm REAL(KIND=dp), INTENT(IN) :: crit_tol ! 'crit_tol'=stopping tolerance (needed in stopping condition) REAL(KIND=dp), INTENT(IN) :: m ! a descent parameter REAL(KIND=dp), INTENT(IN) :: step_tol ! a step-length tolerance INTEGER :: mrounds ! the maximum number of possible rounds in the 'Clarke statinonary' algorithm INTEGER, INTENT(OUT) :: reason_for_stop ! 0 - a new iteration point found ! 1 - stopping condition satisfied (approximate Clarke stationarity) ! 2 - approximate stopping condition satisfied (i.e. the step-length beta* < beta_0) ! 5 - the maximum number of rounds executed in 'Clarke stationary' algorithm INTEGER, INTENT(OUT) :: iter_counter ! the number of rounds used in 'Clarke stationary' algorithm INTEGER, INTENT(OUT) :: f_counter ! the number of function values evaluated for H in 'Clarke stationary' algorithm (same for f_1 and f_2) INTEGER, INTENT(OUT) :: subgrad_counter ! the number of subgradients calculated for H in 'Clarke stationary' algorithm INTEGER, INTENT(IN) :: iprint ! the variable that specifies print: ! iprint = 0 : print is suppressed ! iprint = 1 : basic print of final result (nothing is printed in main iteration) ! iprint = -1: basic print of final result without the solution vector (nothing is printed in main iteration) ! iprint = 2 : extended print of final result (nothing is printed in main iteration) ! iprint = -2: extended print of final result without the solution vector (nothing is printed in main iteration) ! iprint = 3 : basic print of intermediate results and extended print of final results (nothing is printed in main iteration) ! iprint = -3: basic print of intermediate results and extended print of final results without the solution vector (nothing is printed in main iteration) ! iprint = 4 : extended print of intermediate and final results (something is printed about the 'main iteration' (this print is done in the bundle algorithm)) ! iprint = -4: extended print of intermediate and final results without the solution vector (something is printed about the 'main iteration' (this print is done in the bundle algorithm)) ! iprint = 5 : everything is printed step by step (this is the only option which causes print in the 'main iteration') !***************************** LOCAL VARIABLES ************************************ TYPE(kimppu1) :: B ! the bundle B of H containing subgradients from \partial H(x) REAL(KIND=dp), DIMENSION(user_n) :: d ! the search direction REAL(KIND=dp), DIMENSION(user_n) :: new_grad ! the subgradient of H at x REAL(KIND=dp), DIMENSION(user_n) :: new_grad1, new_grad2 ! the subgradients of H_1 and H_2 at x_d REAL(KIND=dp), DIMENSION(user_n) :: u_k ! the vector yielding minimum norm for quadratic norm minimization problem REAL(KIND=dp), DIMENSION(user_n) :: y ! the new auxiliary point REAL(KIND=dp), DIMENSION(user_n) :: test_y ! the test auxiliary point REAL(KIND=dp), DIMENSION(user_n) :: lower_y ! the 'hepl' auxiliary point REAL(KIND=dp), DIMENSION(number_of_func) :: fg1_y_all ! the values of DC components f1_i and g1_l at y REAL(KIND=dp), DIMENSION(number_of_func) :: fg2_y_all ! the values of DC components f2_i and g2_l at y REAL(KIND=dp), DIMENSION(number_of_func) :: AB_y_all ! the values of components A_i and B_l at y REAL(KIND=dp), DIMENSION(number_of_func) :: fg1_test_all ! the values of DC components f1_i and g1_l at test_y REAL(KIND=dp), DIMENSION(number_of_func) :: fg2_test_all ! the values of DC components f2_i and g2_l at test_y REAL(KIND=dp), DIMENSION(number_of_func) :: AB_test_all ! the values of components A_i and B_l at test_y REAL(KIND=dp), DIMENSION(number_of_func) :: fg1_lower_all ! the values of DC components f1_i and g1_l at lower_y REAL(KIND=dp), DIMENSION(number_of_func) :: fg2_lower_all ! the values of DC components f2_i and g2_l at lower_y REAL(KIND=dp), DIMENSION(number_of_func) :: AB_lower_all ! the values of components A_i and B_l at lower_y REAL(KIND=dp), DIMENSION(number_of_func,user_n) :: fg1_grad_all ! the subgradients of DC components f1_i and g1_l REAL(KIND=dp), DIMENSION(number_of_func,user_n) :: fg2_grad_all ! the subgradients of DC components f2_i and g2_l REAL(KIND=dp) :: H_k ! the value of the objective function H at the current iteration point REAL(KIND=dp) :: obj_u ! the optimal value of the quadratin norm minimization problem REAL(KIND=dp) :: norm_u ! the norm for the vector u_k REAL(KIND=dp) :: eps ! stepsize used in subgradient calculation REAL(KIND=dp) :: div_eps ! 1.0_dp / eps REAL(KIND=dp) :: real_decrease ! the real decrease in the objective function f REAL(KIND=dp) :: H1_y, H2_y ! the function values of H_1 and H_2 at y REAL(KIND=dp) :: test_H1, test_H2 ! the function values of H_1 and H_2 at test_y REAL(KIND=dp) :: lower_H1, lower_H2 ! the function values of H_1 and H_2 at lower_y REAL(KIND=dp) :: direc_der ! directional derivative of H REAL(KIND=dp) :: stepsize ! stepsize REAL(KIND=dp) :: upper, lower ! bounds for stepsize REAL(KIND=dp) :: descent_app ! approximation of descent in the objective function value INTEGER :: size_b ! the biggest possible bundle size for B INTEGER :: N ! the current bundle size for B with the current element INTEGER :: i, j, l LOGICAL :: ylos ! .TRUE. if 'Clarke stationary' algorithm can be stopped LOGICAL :: stop_alg ! .TRUE. if 'Clarke stationary' algorithm can be stopped LOGICAL :: stop_step_det ! .TRUE. if the stepsize determination can be stopped !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! <> <> <> <> <> <> <> ALGORITHM STARTS <> <> <> <> <> <> <> <> !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ IF(iprint==5) THEN !prints everything WRITE(*,*) '-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*' WRITE(*,*) ' CLARKE STATIONARY ALGORITHM STARTS' WRITE(*,*) '-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*' WRITE(*,*) ' ' END IF !_______________________________________________________________________________ !************************ STEP 0: INITIALIZATION ******************************* ! the size of the bundle B IF (user_size < 2) THEN size_b = (user_n + 5)*2 ELSE size_b = user_size END IF eps = (10.0_dp)**(-5) ! Stepsize used when subgradient is calculated div_eps = 1.0_dp / eps ! Initialization of iteration counters iter_counter = 1 ! The first iteration round is executed f_counter = 0 ! The number of function values evaluated for f subgrad_counter = 0 ! The numeber of subgradients evaluated for f l = 0 H_k = H1_k - H2_k d = -d_cause / SQRT(DOT_PRODUCT(d_cause, d_cause)) ! The bundles B is initialized CALL init_bundle_b1(B, user_size, user_n) ! Nothing is stored into bundle B N = 0 stop_alg = .FALSE. IF(iprint==5) THEN !prints everything WRITE(*,*)'' WRITE(*,*)' STEP 0: INITIALIATION ' WRITE(*,*)'' END IF !_______________________________________________________________________________ !******************** STEP 0: END ********************************************** y = x_k + (eps) * d DO i = 1, number_of_obj ! the values of DC components f1_i and f2_i at y fg1_y_all(i) = f1(y,problem1(i),factor(i)) fg2_y_all(i) = f2(y,problem2(i),factor(i)) END DO DO i = number_of_obj+1, number_of_func ! the values of DC components g1_l and g2_l at y fg1_y_all(i) = g1(y,problem1(i),factor(i)) fg2_y_all(i) = g2(y,problem2(i),factor(i)) END DO f_counter = f_counter + 1 ! one more objective function value calculated for DC components DO i = 1, number_of_func ! the values of components A_i and B_l for H_1 at y AB_y_all(i) = AorB_i(fg1_y_all, fg2_y_all, fg1_k_all, fg2_k_all ,i) END DO !---------------------------------------------------------------------------------- ! <> <> <> <> <> <> <> REPEATABLE PART BEGINS <> <> <> <> <> <> <> <> !---------------------------------------------------------------------------------- DO WHILE ( (.NOT. stop_alg) .AND. (iter_counter <= mrounds)) !_______________________________________________________________________________ !************************ STEP 1: NEW SUBGRADIENT ****************************** IF (iter_counter > 100) THEN eps = (10.0_dp)**(-6) END IF DO i = 1, number_of_obj ! the values of subgradients for DC components f1_i and f2_i at x_0 new_grad1 = subgradient_f1(y, problem1(i),factor(i)) new_grad2 = subgradient_f2(y, problem2(i),factor(i)) DO j = 1, user_n fg1_grad_all(i,j) = new_grad1(j) fg2_grad_all(i,j) = new_grad2(j) END DO END DO DO i = number_of_obj+1,number_of_func ! the values of subgradients for DC components g1_l and g2_l at x_0 new_grad1 = subgradient_g1(y, problem1(i),factor(i)) new_grad2 = subgradient_g2(y, problem2(i),factor(i)) DO j = 1, user_n fg1_grad_all(i,j) = new_grad1(j) fg2_grad_all(i,j) = new_grad2(j) END DO END DO subgrad_counter = subgrad_counter+1 ! one more subgradient calculated for DC components new_grad1 = subgradient_H1(fg1_grad_all, fg2_grad_all, AB_y_all) new_grad2 = subgradient_H2(fg2_grad_all) ! new subgradient for H at x new_grad = new_grad1 - new_grad2 IF (N == 0) THEN CALL add_first_element_b1(B, new_grad,0) N = N + 1 ELSE CALL add_element_b1(B, new_grad, 0.0_dp,0) END IF IF(iprint==5) THEN !prints everything WRITE(*,*)'-------------------------------------------------------------------------------' WRITE(*,*)'' WRITE(*,*)' STEP 1: NEW SUBGRADIENT ' WRITE(*,*)'' WRITE(*,*)' Iteration round: ', iter_counter WRITE(*,*)' New subgardient = [', new_grad, ']' END IF !_______________________________________________________________________________ !******************** STEP 1: END ********************************************** !_______________________________________________________________________________ !************************ STEP 2: CALRKE STATIONARITY ************************** N = give_size_b1(B) + 1 ! the current bundle size CALL quadratic_solver(x_k, user_n, N, B, 1.0_dp, u_k, obj_u) ! The solution u_k is has already negative sign. Thus d = u_k/norm_u norm_u = SQRT( DOT_PRODUCT(u_k,u_k)) IF(iprint==5) THEN !prints everything WRITE(*,*)'-------------------------------------------------------------------------------' WRITE(*,*)'' WRITE(*,*)' STEP 2: CLARKE STATIONARITY ' WRITE(*,*)'' WRITE(*,*)' The value of minimum norm = ', norm_u WRITE(*,*)' The solution to quadratic norm minimization problem = [',u_k, ']' END IF IF (norm_u < crit_tol) THEN reason_for_stop = 1 stop_alg = .TRUE. IF(iprint==5) THEN !prints everything WRITE(*,*)'' WRITE(*,*)' CLARKE STATIONARITY ACHIEVED' WRITE(*,*)'' END IF !_______________________________________________________________________________ !******************** STEP 2: END ********************************************** ELSE !_______________________________________________________________________________ !************************ STEP 3: SERACH DIRECTION ***************************** d = u_k / norm_u y = x_k + eps * d DO i = 1, number_of_obj ! the values of DC components f1_i anf f2_i at y fg1_y_all(i) = f1(y,problem1(i),factor(i)) fg2_y_all(i) = f2(y,problem2(i),factor(i)) END DO DO i = number_of_obj+1,number_of_func ! the values of DC components g1_l and g2_l at y fg1_y_all(i) = g1(y,problem1(i),factor(i)) fg2_y_all(i) = g2(y,problem2(i),factor(i)) END DO f_counter = f_counter + 1 ! one more objective function value calculated for DC components DO i = 1, number_of_func ! the values of components A_i and B_l for H_1 at y AB_y_all(i) = AorB_i(fg1_y_all, fg2_y_all, fg1_k_all, fg2_k_all ,i) END DO H1_y = H1(AB_y_all) H2_y = H2(fg2_y_all) real_decrease = H1_y - H2_y - H_k direc_der = div_eps * real_decrease IF(iprint==5) THEN !prints everything WRITE(*,*)'-------------------------------------------------------------------------------' WRITE(*,*)'' WRITE(*,*)' STEP 3: SEARCH DIRECTION ' WRITE(*,*)'' WRITE(*,*)' The search direction d = [', d, ']' WRITE(*,*)'' END IF IF (direc_der <= (-m * norm_u)) THEN reason_for_stop = 10 ! we will execute step-length determination at Step 4 stop_alg = .TRUE. IF(iprint==5) THEN !prints everything WRITE(*,*)'' WRITE(*,*)' SUFFICIENT DESCENT IS OBTAINED.' WRITE(*,*)'' END IF END IF iter_counter = iter_counter +1 END IF !_______________________________________________________________________________ !******************** STEP 3: END ********************************************** END DO !---------------------------------------------------------------------------------- ! <> <> <> <> <> <> <> REPEATABLE PART ENDS <> <> <> <> <> <> <> <> !---------------------------------------------------------------------------------- !_______________________________________________________________________________ !************************ STEP 4: STEP-LENGTH ********************************** IF ( (reason_for_stop > 1) .AND. stop_alg ) THEN !-**--**- STEPSIZE DETERMINATION START -**--**- upper = 1.0_dp lower = user_step_tol stepsize = user_step_tol stop_step_det = .FALSE. ylos = .FALSE. lower_y = x_k + stepsize * d DO i = 1, number_of_obj ! the values of DC components f1_i and f2_i at lower_y fg1_lower_all(i) = f1(lower_y,problem1(i),factor(i)) fg2_lower_all(i) = f2(lower_y,problem2(i),factor(i)) END DO DO i = number_of_obj+1,number_of_func ! the values of DC components g1_l and g2_l at lower_y fg1_lower_all(i) = g1(lower_y,problem1(i),factor(i)) fg2_lower_all(i) = g2(lower_y,problem2(i),factor(i)) END DO f_counter = f_counter+1 DO i = 1, number_of_func ! the values of components A_i and B_l for H_1 at lower_y AB_lower_all(i) = AorB_i(fg1_lower_all, fg2_lower_all, & & fg1_k_all, fg2_k_all ,i) END DO lower_H1 = H1(AB_lower_all) ! the value of H_1 at test_y lower_H2 = H2(fg2_lower_all) ! the value of H_2 at test_y descent_app = -m * stepsize * norm_u IF ((lower_H1 - lower_H2 - H_k) > descent_app ) THEN stop_step_det = .TRUE. stepsize = -1.0_dp END IF l = 1 DO WHILE((.NOT. stop_step_det)) IF (l == 1) THEN stepsize = 100.0_dp * stepsize ELSE IF (l == 2) THEN IF (upper < 0.02_dp) THEN stepsize = 0.1_dp * stepsize !stepsize = 0.5_dp*(upper +lower) ELSE stepsize = 10.0_dp*stepsize END IF ELSE IF (l > 2) THEN stepsize = 0.5_dp*(upper +lower) END IF test_y = x_k + stepsize * d DO i = 1, number_of_obj ! the values of DC components f1_i and f2_i at test_y fg1_test_all(i) = f1(test_y,problem1(i),factor(i)) fg2_test_all(i) = f2(test_y,problem2(i),factor(i)) END DO DO i = number_of_obj+1,number_of_func ! the values of DC components g1_l and g2_l at test_y fg1_test_all(i) = g1(test_y,problem1(i),factor(i)) fg2_test_all(i) = g2(test_y,problem2(i),factor(i)) END DO f_counter = f_counter + 1 ! one more objective function value calculated for a DC component DO i = 1, number_of_func ! the values of components A_i and B_l for H_1 at test_y AB_test_all(i) = AorB_i(fg1_test_all, fg2_test_all, & & fg1_k_all, fg2_k_all ,i) END DO test_H1 = H1(AB_test_all) ! the value of H_1 at test_y test_H2 = H2(fg2_test_all) ! the value of H_2 at test_y IF (l < 7) THEN descent_app = -m * stepsize * norm_u IF ( (test_H1 - test_H2 - H_k) < descent_app) THEN fg1_lower_all = fg1_test_all fg2_lower_all = fg2_test_all lower_H1 = test_H1 lower_H2 = test_H2 lower_y = test_y l = l + 1 lower = stepsize ELSE upper = stepsize l = l+1 END IF ELSE stop_step_det = .TRUE. y = lower_y H1_y = lower_H1 H2_y = lower_H2 fg1_y_all = fg1_lower_all fg2_y_all = fg2_lower_all END IF END DO !-**--**- STEPSIZE DETERMINATION END -**--**- IF(iprint==5) THEN !prints everything WRITE(*,*)'-------------------------------------------------------------------------------' WRITE(*,*)'' WRITE(*,*)' STEP 4: STEP-LENGTH ' WRITE(*,*)'' WRITE(*,*)' The step-lenght t =', stepsize END IF IF (stepsize >= step_tol) THEN ! the step-length is big enough x_new = y ! the new iteration point is the auxilary point y H1_new = H1_y ! the value of H_1 at x_new H2_new = H2_y ! the value of H_2 at x_new fg1_new_all = fg1_y_all ! the values of f1_i and g1_l at x_new fg2_new_all = fg2_y_all ! the values of f2_i and g2_l at x_new reason_for_stop = 0 ! the reason for stopping is the new iteration point IF(iprint==5) THEN !prints everything WRITE(*,*)' New itaration point is found and x_new = [', x_new, ']' WRITE(*,*)' New objective function value f_new = ', H1_new -H2_new WRITE(*,*)'' END IF ELSE reason_for_stop = 2 IF(iprint==5) THEN !prints everything WRITE(*,*)' APPROXIMATE STOPPING CONDITION HOLD at the current iteration point.' WRITE(*,*)'' END IF END IF iter_counter = iter_counter -1 END IF !_______________________________________________________________________________ !******************** STEP 4: END ********************************************** IF ( (.NOT. stop_alg ) ) THEN reason_for_stop = 5 ! the maximum number of rounds have been executed IF(iprint==5) THEN !prints everything WRITE(*,*)'--------------------------------------& &-----------------------------------------' WRITE(*,*)'' WRITE(*,*)' THE MAXIMUM NUMBER OF ROUNDS HAS BEEN EXECUTED.' WRITE(*,*)'' END IF END IF !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! <> <> <> <> <> <> <> ALGORITHM ENDS <> <> <> <> <> <> <> <> !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ IF(iprint==5) THEN !prints everything WRITE(*,*) '-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*' WRITE(*,*) ' CLARKE STATIONARY ALGORITHM ENDS' WRITE(*,*) '-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*' WRITE(*,*) ' ' END IF END SUBROUTINE guaranteeing_clarke !....................................................................................... ! <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> ! _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ !| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | !** END ** END ** END ** END ** END ** END ** END ** END ** END ** END ** END ** END ** !|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_| ! _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ !| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | ! START ** START ** START ** START ** START ** START ** START ** START ** START ** START !|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_| ! <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> !*************************************************************************************** ! ---------------------------------------------------------------------------------- | ! | | | ! | QUADRATIC NORM SOLVER | | ! | | | ! ---------------------------------------------------------------------------------- | !*************************************************************************************** SUBROUTINE quadratic_solver(x, NF, NA, B, t, d, obj) ! Solves the quadratic norm problem in the 'Clarke stationary' algorithm. ! Calls PLQDF1 by Ladislav Luksan IMPLICIT NONE !**************************** NEEDED FROM USER *********************************** TYPE(kimppu1), INTENT(IN) :: B ! the bundle B REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: x ! the current iteration point REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: d ! the vector giving minimum norm REAL(KIND=dp), INTENT(IN) :: t ! the proximity parameter t REAL(KIND=dp), INTENT(OUT) :: obj ! the optimal value of the quadratic norm problem INTEGER, INTENT(IN) :: NF ! the number of variables INTEGER, INTENT(IN) :: NA ! the bundle size of B is 'give_size_b1(B) + 1' !*************************** LOCAL VARIABLES ************************************ ! .. Parameters used in PLQDF1 .. REAL(KIND=dp), PARAMETER :: ETA0 = 1.0E-15_dp REAL(KIND=dp), PARAMETER :: ETA2 = 1.0E-12_dp REAL(KIND=dp), PARAMETER :: ETA9 = 1.0E+60_dp REAL(KIND=dp), PARAMETER :: EPS7 = 1.0E-14_dp REAL(KIND=dp), PARAMETER :: EPS9 = 1.0E-12_dp ! .. Scalar Arguments in PLQDF1 .. REAL(KIND=dp) :: GMAX ! output of PLQDF1: maximum absolute value of a partial derivative REAL(KIND=dp) :: UMAX ! output of PLQDF1: maximum absolute value of a negative lagrange multiplier REAL(KIND=dp) :: XNORM ! output of PLQDF1: value of linearized minimax function = delta_1 + delta_2 REAL(KIND=dp) :: u ! 'help' variable INTEGER :: j, k, l INTEGER, PARAMETER :: NC = 0 ! number of constraints is zero !INTEGER :: IDECF=10 ! IDECF=10 diagonal matrix !INTEGER :: MFP=2 ! MFP=2 optimum feasible point INTEGER :: KBC=0, KBF=0 ! KBC=0 no linear constraints; KBF=0 no simple bounds INTEGER :: ITERQ, N ! output values of PLQDF1: ITERQ=type of feasible point ! N=dimension of manifold defined by active constraints ! .. Array Arguments in PLQDF1.. INTEGER, DIMENSION(NF) :: IX ! vector containing types of bounds INTEGER, DIMENSION(NA) :: IA ! vector containing types of deviations INTEGER, DIMENSION(NC) :: IC ! vector containing types of constraints. NOT significant because NC=0 INTEGER, DIMENSION(NF+1) :: IAA ! Output of PLQDF1: vector containing indicies of active functions REAL(KIND=dp), DIMENSION(NF) :: XL, XU ! lower and upper bounds for x. NOT significant because variables are unbounded REAL(KIND=dp), DIMENSION(NF) :: H ! diagonal matrix (1/t)*I REAL(KIND=dp), DIMENSION(NC) :: CF, CL, CU ! NOT significant since NC=0 (are related to constraints) REAL(KIND=dp), DIMENSION(NA) :: AF ! vector of bundle function values (-alpha) REAL(KIND=dp), DIMENSION(NA) :: AFD ! Output of PLQDF1: vector containing increments of the approximated functions REAL(KIND=dp), DIMENSION(NF+1) :: AZ ! Output of PLQDF1: vector of Lagrange multipliers REAL(KIND=dp), DIMENSION(NF+1) :: S ! Output of PLQDF1: direction vector REAL(KIND=dp), DIMENSION(NF+1) :: G ! Output of PLQDF1: gradient of the Lagrangian function REAL(KIND=dp), DIMENSION(NF*NA) :: AG ! matrix whose columns are bundle subgradients REAL(KIND=dp), DIMENSION((NF+1)*(NF+2)/2) :: AR ! Output of PLQDF1: triangular decomposition of kernel of the orthogonal projection REAL(KIND=dp), DIMENSION(NF*NC) :: CG ! NOT significant since NC=0. matrix whose columns are normals of the linear constraints ! .. Some other varibles .. REAL(KIND=dp), DIMENSION(NF*NA) :: grad_m_b ! subgradient matrix of B !************************** SUBPROBLEM SOLVER STARTS ********************************* !****************************** INITIALIZATIONS ************************************** IX = 0 ! types of bounds: 0 - unbounded variables IA = 2 ! types of deviations u = 1.0_dp / t H = u ! diagonal matrix grad_m_b = grad_matrix(B) ! subgradient matrix of B ! NO linearization errors AF = 0.0_dp ! matrix containing subgradients DO j = 1, NA k = (j-1)*NF DO l = 1, NF AG(k+l) = grad_m_b(k+l) END DO END DO !Calls PLQDF1 by Ladislav Luksan CALL PLQDF1(NF,NA,NC,X,IX,XL,XU,AF,AFD,IA,IAA, & & AG,AR,AZ,CF,IC,CL,CU,CG,G,H,S,2,KBF, & & KBC,10,ETA0,ETA2,ETA9,EPS7,EPS9,XNORM, & & UMAX,GMAX,N,ITERQ) ! the vector giving minimum norm DO j = 1, NF d(j) = S(j) END DO ! the value of minimum norm obj = 0.5_dp * DOT_PRODUCT(d, d) END SUBROUTINE quadratic_solver !....................................................................................... ! <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> ! _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ !| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | !** END ** END ** END ** END ** END ** END ** END ** END ** END ** END ** END ** END ** !|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_| ! _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ !| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | ! START ** START ** START ** START ** START ** START ** START ** START ** START ** START !|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_| ! <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> !*************************************************************************************** ! ---------------------------------------------------------------------------------- | ! | | | ! | SUBPROBLEM SOLVER | | ! | | | ! ---------------------------------------------------------------------------------- | !*************************************************************************************** SUBROUTINE subproblem_solver(x, AB_k_all, NF, NA, B1, B2, t, subprob_counter) ! Solves the subproblems in the original search direction problem. ! Calls PLQDF1 by Ladislav Luksan IMPLICIT NONE !**************************** NEEDED FROM USER *********************************** TYPE(kimppu1), INTENT(IN) :: B1 ! bundle B_1 TYPE(kimppu2), INTENT(INOUT) :: B2 ! bundle B_2 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: x ! current iteration point REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: AB_k_all ! the values of components A_i and B_l at the current iteration point REAL(KIND=dp), INTENT(IN) :: t ! proximity parameter t INTEGER, INTENT(IN) :: NF ! number of variables INTEGER, INTENT(IN) :: NA ! bundle size of B1 INTEGER, INTENT(OUT) :: subprob_counter ! number of subproblems solved !*************************** LOCAL VARIABLES ************************************ ! .. Parameters used in PLQDF1 .. REAL(KIND=dp), PARAMETER :: ETA0 = 1.0E-15_dp REAL(KIND=dp), PARAMETER :: ETA2 = 1.0E-12_dp REAL(KIND=dp), PARAMETER :: ETA9 = 1.0E+60_dp REAL(KIND=dp), PARAMETER :: EPS7 = 1.0E-14_dp REAL(KIND=dp), PARAMETER :: EPS9 = 1.0E-12_dp ! .. Scalar Arguments in PLQDF1 .. REAL(KIND=dp) :: GMAX ! output of PLQDF1: maximum absolute value of a partial derivative REAL(KIND=dp) :: UMAX ! output of PLQDF1: maximum absolute value of a negative lagrange multiplier REAL(KIND=dp) :: XNORM ! output of PLQDF1: value of linearized minimax function = delta_1 + delta_2 REAL(KIND=dp) :: alpha_b2 ! a linearization error of B_2 REAL(KIND=dp) :: u, obj ! 'help' variable REAL(KIND=dp) :: a INTEGER :: i, j, k, l INTEGER, PARAMETER :: NC = 0 ! number of constraints is zero !INTEGER :: IDECF=10 ! IDECF=10 diagonal matrix; !INTEGER :: MFP=2 ! MFP=2 optimum feasible point INTEGER :: KBC=0, KBF=0 ! KBC=0 no linear constraints; KBF=0 no simple bounds INTEGER :: ITERQ, N ! output values of PLQDF1: ITERQ=type of feasible point ! N=dimension of manifold defined by active constraints ! .. Array Arguments in PLQDF1.. INTEGER, DIMENSION(NF) :: IX ! vector containing types of bounds INTEGER, DIMENSION(NA) :: IA ! vector containing types of deviations INTEGER, DIMENSION(NC) :: IC ! vector containing types of constraints. NOT significant because NC=0 INTEGER, DIMENSION(NF+1) :: IAA ! Output of PLQDF1: vector containing indicies of active functions REAL(KIND=dp), DIMENSION(NF) :: XL, XU ! lower and upper bounds for x. NOT significant because variables are unbounded REAL(KIND=dp), DIMENSION(NF) :: H ! diagonal matrix (1/t)*I REAL(KIND=dp), DIMENSION(NC) :: CF, CL, CU ! NOT significant since NC=0 (are related to constraints) REAL(KIND=dp), DIMENSION(NA) :: AF ! vector of bundle function values (-alpha) REAL(KIND=dp), DIMENSION(NA) :: AFD ! Output of PLQDF1: vector containing increments of the approximated functions REAL(KIND=dp), DIMENSION(NF+1) :: AZ ! Output of PLQDF1: vector of Lagrange multipliers REAL(KIND=dp), DIMENSION(NF+1) :: S ! Output of PLQDF1: direction vector REAL(KIND=dp), DIMENSION(NF) :: direction ! Actual direction vector used (notice dimension) REAL(KIND=dp), DIMENSION(NF+1) :: G ! Output of PLQDF1: gradient of the Lagrangian function REAL(KIND=dp), DIMENSION(NF*NA) :: AG ! matrix whose columns are bundle subgradients REAL(KIND=dp), DIMENSION((NF+1)*(NF+2)/2) :: AR ! Output of PLQDF1: triangular decomposition of kernel of the orthogonal projection REAL(KIND=dp), DIMENSION(NF*NC) :: CG ! NOT significant since NC=0. matrix whose columns are normals of the linear constraints ! .. Some other varibles .. REAL(KIND=dp), DIMENSION(NF*NA) :: grad_m_b1 ! subgradient matrix of B_1 REAL(KIND=dp), DIMENSION(NA) :: alpha_m_b1 ! linearization error matrix of B_1 REAL(KIND=dp), DIMENSION(NF) :: grad_b2 ! a subgradient of B_2 !************************** SUBPROBLEM SOLVER STARTS ********************************* !****************************** INITIALIZATIONS ************************************** IX = 0 ! types of bounds: 0 - unbounded variables IA = 2 ! types of deviations u = 1.0_dp / t H = u ! diagonal matrix grad_m_b1 = grad_matrix(B1) ! subgradient matrix of B_1 alpha_m_b1 = lin_err_and_f1_matrix(B1,AB_k_all) ! linearization error matrix of B_1 with component values A_i(x_h) and B_l(x_h) subprob_counter = give_size_b2(B2) !$OMP PARALLEL DO PRIVATE(grad_b2,alpha_b2,direction,a,obj) & !$OMP FIRSTPRIVATE(NA,X,IX,XL,XU,AFD,IA,IAA) & !$OMP FIRSTPRIVATE(AR,AZ,CF,IC,CL,CU,CG,G,H,S,KBF) & !$OMP FIRSTPRIVATE(KBC,XNORM,UMAX,GMAX,N,ITERQ) & !$OMP PRIVATE(i, AG, AF ) & !$OMP SHARED(grad_m_b1,alpha_m_b1,B2) !->->->->->->->->->->->->->-> EACH SUBPROBLEM SOLVED <-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<- subproblems1: DO i = 0, subprob_counter ! each subproblem is looked through grad_b2 = give_subgrad_b2(B2, i) ! subgradient of B_2 in the subproblem i alpha_b2 = give_linerr_b2(B2, i) ! linearization error of B_2 in the subproblem i !$OMP CRITICAL DO j = 1, NA k = (j-1)*NF DO l = 1, NF AG(k+l) = grad_m_b1(k+l) - grad_b2(l) END DO AF(j) = - alpha_m_b1(j) + alpha_b2 END DO !$OMP END CRITICAL !Calls PLQDF1 by Ladislav Luksan CALL PLQDF1(NF,NA,NC,X,IX,XL,XU,AF,AFD,IA,IAA, & & AG,AR,AZ,CF,IC,CL,CU,CG,G,H,S,2,KBF, & & KBC,10,ETA0,ETA2,ETA9,EPS7,EPS9,XNORM, & & UMAX,GMAX,N,ITERQ) DO j = 1, NF direction(j) = S(j) END DO a = DOT_PRODUCT(direction, direction) obj = XNORM + (a * u) / 2 !$OMP CRITICAL CALL add_solution(B2, i , direction, XNORM, obj ) !$OMP END CRITICAL END DO subproblems1 !->->->->->->->->->->->->->-> EACH SUBPROBLEM SOLVED END <-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<-<- !$OPM END PARALLEL DO subprob_counter = subprob_counter + 1 END SUBROUTINE subproblem_solver !....................................................................................... ! <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> ! _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ !| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | !** END ** END ** END ** END ** END ** END ** END ** END ** END ** END ** END ** END ** !|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_| END MODULE mdbdc