1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief routines that optimize a functional using the limited memory bfgs 8!> quasi-newton method. 9!> The process set up so that a master runs the real optimizer and the 10!> others help then to calculate the objective function. 11!> The arguments for the objective function are physically present in 12!> every processor (nedeed in the actual implementation of pao). 13!> In the future tha arguments themselves could be distributed. 14!> \par History 15!> 09.2003 globenv->para_env, retain/release, better parallel behaviour 16!> \author Fawzi Mohamed 17!> @version 2.2002 18! ************************************************************************************************** 19MODULE cp_lbfgs_optimizer_gopt 20 USE cp_lbfgs, ONLY: setulb 21 USE cp_log_handling, ONLY: cp_get_default_logger,& 22 cp_logger_type,& 23 cp_to_string 24 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 25 cp_print_key_unit_nr 26 USE cp_para_env, ONLY: cp_para_env_release,& 27 cp_para_env_retain 28 USE cp_para_types, ONLY: cp_para_env_type 29 USE force_env_types, ONLY: force_env_type 30 USE gopt_f_methods, ONLY: gopt_f_io 31 USE gopt_f_types, ONLY: gopt_f_release,& 32 gopt_f_retain,& 33 gopt_f_type 34 USE gopt_param_types, ONLY: gopt_param_type 35 USE input_section_types, ONLY: section_vals_type 36 USE kinds, ONLY: dp 37 USE machine, ONLY: m_walltime 38 USE message_passing, ONLY: mp_bcast 39#include "../base/base_uses.f90" 40 41 IMPLICIT NONE 42 PRIVATE 43#include "gopt_f77_methods.h" 44 45 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE. 46 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_lbfgs_optimizer_gopt' 47 INTEGER, PRIVATE, SAVE :: last_lbfgs_optimizer_id = 0 48 49 ! types 50 PUBLIC :: cp_lbfgs_opt_gopt_type 51 52 ! core methods 53 54 ! special methos 55 56 ! underlying functions 57 PUBLIC :: cp_opt_gopt_create, cp_opt_gopt_release, & 58 cp_opt_gopt_next, & 59 cp_opt_gopt_stop 60 61 ! initialize the object 62 INTERFACE cp_create 63 MODULE PROCEDURE cp_opt_gopt_create 64 END INTERFACE 65 ! releases the given object 66 INTERFACE cp_release 67 MODULE PROCEDURE cp_opt_gopt_release 68 END INTERFACE 69 ! retains the given object 70 INTERFACE cp_retain 71 MODULE PROCEDURE cp_opt_gopt_retain 72 END INTERFACE 73 ! returns attributes about the object 74 INTERFACE cp_get 75 MODULE PROCEDURE cp_opt_gopt_get 76 END INTERFACE 77 ! goes to the next point, returns true if not converged or error 78 INTERFACE cp_next 79 MODULE PROCEDURE cp_opt_gopt_next 80 END INTERFACE 81 ! goes to the next point 82 INTERFACE cp_step 83 MODULE PROCEDURE cp_opt_gopt_step 84 END INTERFACE 85 ! stops the iteration 86 INTERFACE cp_stop 87 MODULE PROCEDURE cp_opt_gopt_stop 88 END INTERFACE 89 90! ************************************************************************************************** 91!> \brief info for the optimizer (see the description of this module) 92!> \param task the actual task of the optimizer (in the master it is up to 93!> date, in case of error also the slaves one get updated. 94!> \param csave internal character string used by the lbfgs optimizer, 95!> meaningful only in the master 96!> \param lsave logical array used by the lbfgs optimizer, updated only 97!> in the master 98!> On exit with task = 'NEW_X', the following information is 99!> available: 100!> lsave(1) = .true. the initial x did not satisfy the bounds; 101!> lsave(2) = .true. the problem contains bounds; 102!> lsave(3) = .true. each variable has upper and lower bounds. 103!> \param ref_count reference count (see doc/ReferenceCounting.html) 104!> \param id_nr identification number (unique) 105!> \param m the dimension of the subspace used to approximate the second 106!> derivative 107!> \param print_every every how many iterations output should be written. 108!> if 0 only at end, if print_every<0 never 109!> \param master the pid of the master processor 110!> \param max_f_per_iter the maximum number of function evaluations per 111!> iteration 112!> \param status 0: just initialized, 1: f g calculation, 113!> 2: begin new iteration, 3: ended iteration, 114!> 4: normal (converged) exit, 5: abnormal (error) exit, 115!> 6: daellocated 116!> \param n_iter the actual iteration number 117!> \param kind_of_bound an array with 0 (no bound), 1 (lower bound), 118!> 2 (both bounds), 3 (upper bound), to describe the bounds 119!> of every variable 120!> \param i_work_array an integer workarray of dimension 3*n, present only 121!> in the master 122!> \param isave is an INTEGER working array of dimension 44. 123!> On exit with task = 'NEW_X', it contains information that 124!> the user may want to access: 125!> \param isave (30) = the current iteration number; 126!> \param isave (34) = the total number of function and gradient 127!> evaluations; 128!> \param isave (36) = the number of function value or gradient 129!> evaluations in the current iteration; 130!> \param isave (38) = the number of free variables in the current 131!> iteration; 132!> \param isave (39) = the number of active constraints at the current 133!> iteration; 134!> \param f the actual best value of the object function 135!> \param wanted_relative_f_delta the wanted relative error on f 136!> (to be multiplied by epsilon), 0.0 -> no check 137!> \param wanted_projected_gradient the wanted error on the projected 138!> gradient (hessian times the gradient), 0.0 -> no check 139!> \param last_f the value of f in the last iteration 140!> \param projected_gradient the value of the sup norm of the projected 141!> gradient 142!> \param x the actual evaluation point (best one if converged or stopped) 143!> \param lower_bound the lower bounds 144!> \param upper_bound the upper bounds 145!> \param gradient the actual gradient 146!> \param dsave info date for lbfgs (master only) 147!> \param work_array a work array for lbfgs (master only) 148!> \param para_env the parallel environment for this optimizer 149!> \param obj_funct the objective function to be optimized 150!> \par History 151!> none 152!> \author Fawzi Mohamed 153!> @version 2.2002 154! ************************************************************************************************** 155 TYPE cp_lbfgs_opt_gopt_type 156 CHARACTER(len=60) :: task 157 CHARACTER(len=60) :: csave 158 LOGICAL :: lsave(4) 159 INTEGER :: m, print_every, master, max_f_per_iter, status, n_iter 160 INTEGER :: ref_count, id_nr 161 INTEGER, DIMENSION(:), POINTER :: kind_of_bound, i_work_array, isave 162 REAL(kind=dp) :: f, wanted_relative_f_delta, wanted_projected_gradient, & 163 last_f, projected_gradient, eold, emin, trust_radius 164 REAL(kind=dp), DIMENSION(:), POINTER :: x, lower_bound, upper_bound, & 165 gradient, dsave, work_array 166 TYPE(cp_para_env_type), POINTER :: para_env 167 TYPE(gopt_f_type), POINTER :: obj_funct 168 END TYPE cp_lbfgs_opt_gopt_type 169 170CONTAINS 171 172! ************************************************************************************************** 173!> \brief initializes the optimizer 174!> \param optimizer ... 175!> \param para_env ... 176!> \param obj_funct ... 177!> \param x0 ... 178!> \param m ... 179!> \param print_every ... 180!> \param wanted_relative_f_delta ... 181!> \param wanted_projected_gradient ... 182!> \param lower_bound ... 183!> \param upper_bound ... 184!> \param kind_of_bound ... 185!> \param master ... 186!> \param max_f_per_iter ... 187!> \param trust_radius ... 188!> \par History 189!> 02.2002 created [fawzi] 190!> 09.2003 refactored (retain/release,para_env) [fawzi] 191!> \author Fawzi Mohamed 192!> \note 193!> redirects the lbfgs output the the default unit 194! ************************************************************************************************** 195 SUBROUTINE cp_opt_gopt_create(optimizer, para_env, obj_funct, x0, m, print_every, & 196 wanted_relative_f_delta, wanted_projected_gradient, lower_bound, upper_bound, & 197 kind_of_bound, master, max_f_per_iter, trust_radius) 198 TYPE(cp_lbfgs_opt_gopt_type), POINTER :: optimizer 199 TYPE(cp_para_env_type), POINTER :: para_env 200 TYPE(gopt_f_type), POINTER :: obj_funct 201 REAL(kind=dp), DIMENSION(:), INTENT(in) :: x0 202 INTEGER, INTENT(in), OPTIONAL :: m, print_every 203 REAL(kind=dp), INTENT(in), OPTIONAL :: wanted_relative_f_delta, & 204 wanted_projected_gradient 205 REAL(kind=dp), DIMENSION(SIZE(x0)), INTENT(in), & 206 OPTIONAL :: lower_bound, upper_bound 207 INTEGER, DIMENSION(SIZE(x0)), INTENT(in), OPTIONAL :: kind_of_bound 208 INTEGER, INTENT(in), OPTIONAL :: master, max_f_per_iter 209 REAL(kind=dp), INTENT(in), OPTIONAL :: trust_radius 210 211 CHARACTER(len=*), PARAMETER :: routineN = 'cp_opt_gopt_create', & 212 routineP = moduleN//':'//routineN 213 214 INTEGER :: handle, lenwa, n 215 216 CALL timeset(routineN, handle) 217 218 ALLOCATE (optimizer) 219 NULLIFY (optimizer%kind_of_bound, & 220 optimizer%i_work_array, & 221 optimizer%isave, & 222 optimizer%x, & 223 optimizer%lower_bound, & 224 optimizer%upper_bound, & 225 optimizer%gradient, & 226 optimizer%dsave, & 227 optimizer%work_array, & 228 optimizer%para_env, & 229 optimizer%obj_funct) 230 optimizer%ref_count = 0 231 last_lbfgs_optimizer_id = last_lbfgs_optimizer_id + 1 232 optimizer%id_nr = last_lbfgs_optimizer_id 233 n = SIZE(x0) 234 optimizer%m = 4 235 IF (PRESENT(m)) optimizer%m = m 236 optimizer%master = para_env%source 237 optimizer%para_env => para_env 238 CALL cp_para_env_retain(para_env) 239 optimizer%obj_funct => obj_funct 240 CALL gopt_f_retain(obj_funct) 241 optimizer%max_f_per_iter = 20 242 IF (PRESENT(max_f_per_iter)) optimizer%max_f_per_iter = max_f_per_iter 243 optimizer%print_every = -1 244 optimizer%n_iter = 0 245 optimizer%f = -1.0_dp 246 optimizer%last_f = -1.0_dp 247 optimizer%projected_gradient = -1.0_dp 248 IF (PRESENT(print_every)) optimizer%print_every = print_every 249 IF (PRESENT(master)) optimizer%master = master 250 IF (optimizer%master == optimizer%para_env%mepos) THEN 251 !MK This has to be adapted for a new L-BFGS version possibly 252 lenwa = 2*optimizer%m*n + 5*n + 11*optimizer%m*optimizer%m + 8*optimizer%m 253 ALLOCATE (optimizer%kind_of_bound(n), optimizer%i_work_array(3*n), & 254 optimizer%isave(44)) 255 ALLOCATE (optimizer%x(n), optimizer%lower_bound(n), & 256 optimizer%upper_bound(n), optimizer%gradient(n), & 257 optimizer%dsave(29), optimizer%work_array(lenwa)) 258 optimizer%x = x0 259 optimizer%task = 'START' 260 optimizer%wanted_relative_f_delta = wanted_relative_f_delta 261 optimizer%wanted_projected_gradient = wanted_projected_gradient 262 optimizer%kind_of_bound = 0 263 IF (PRESENT(kind_of_bound)) optimizer%kind_of_bound = kind_of_bound 264 IF (PRESENT(lower_bound)) optimizer%lower_bound = lower_bound 265 IF (PRESENT(upper_bound)) optimizer%upper_bound = upper_bound 266 IF (PRESENT(trust_radius)) optimizer%trust_radius = trust_radius 267 268 CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, & 269 optimizer%lower_bound, optimizer%upper_bound, & 270 optimizer%kind_of_bound, optimizer%f, optimizer%gradient, & 271 optimizer%wanted_relative_f_delta, & 272 optimizer%wanted_projected_gradient, optimizer%work_array, & 273 optimizer%i_work_array, optimizer%task, optimizer%print_every, & 274 optimizer%csave, optimizer%lsave, optimizer%isave, & 275 optimizer%dsave, optimizer%trust_radius) 276 ELSE 277 NULLIFY ( & 278 optimizer%kind_of_bound, optimizer%i_work_array, optimizer%isave, & 279 optimizer%lower_bound, optimizer%upper_bound, optimizer%gradient, & 280 optimizer%dsave, optimizer%work_array) 281 ALLOCATE (optimizer%x(n)) 282 optimizer%x(:) = 0.0_dp 283 ALLOCATE (optimizer%gradient(n)) 284 optimizer%gradient(:) = 0.0_dp 285 END IF 286 CALL mp_bcast(optimizer%x, optimizer%master, optimizer%para_env%group) 287 optimizer%status = 0 288 optimizer%ref_count = 1 289 290 CALL timestop(handle) 291 292 END SUBROUTINE cp_opt_gopt_create 293 294! ************************************************************************************************** 295!> \brief retains the given optimizer (see doc/ReferenceCounting.html) 296!> \param optimizer the optimizer to retain 297!> \par History 298!> 08.2003 created [fawzi] 299!> \author Fawzi Mohamed 300! ************************************************************************************************** 301 SUBROUTINE cp_opt_gopt_retain(optimizer) 302 TYPE(cp_lbfgs_opt_gopt_type), POINTER :: optimizer 303 304 CHARACTER(len=*), PARAMETER :: routineN = 'cp_opt_gopt_retain', & 305 routineP = moduleN//':'//routineN 306 307 CPASSERT(ASSOCIATED(optimizer)) 308 CPASSERT(optimizer%ref_count > 0) 309 optimizer%ref_count = optimizer%ref_count + 1 310 END SUBROUTINE cp_opt_gopt_retain 311 312! ************************************************************************************************** 313!> \brief releases the optimizer (see doc/ReferenceCounting.html) 314!> \param optimizer the object that should be released 315!> \par History 316!> 02.2002 created [fawzi] 317!> 09.2003 dealloc_ref->release [fawzi] 318!> \author Fawzi Mohamed 319! ************************************************************************************************** 320 SUBROUTINE cp_opt_gopt_release(optimizer) 321 TYPE(cp_lbfgs_opt_gopt_type), POINTER :: optimizer 322 323 CHARACTER(len=*), PARAMETER :: routineN = 'cp_opt_gopt_release', & 324 routineP = moduleN//':'//routineN 325 326 INTEGER :: handle 327 328 CALL timeset(routineN, handle) 329 330 IF (ASSOCIATED(optimizer)) THEN 331 CPASSERT(optimizer%ref_count > 0) 332 optimizer%ref_count = optimizer%ref_count - 1 333 IF (optimizer%ref_count == 0) THEN 334 optimizer%status = 6 335 IF (ASSOCIATED(optimizer%kind_of_bound)) THEN 336 DEALLOCATE (optimizer%kind_of_bound) 337 END IF 338 IF (ASSOCIATED(optimizer%i_work_array)) THEN 339 DEALLOCATE (optimizer%i_work_array) 340 END IF 341 IF (ASSOCIATED(optimizer%isave)) THEN 342 DEALLOCATE (optimizer%isave) 343 END IF 344 IF (ASSOCIATED(optimizer%x)) THEN 345 DEALLOCATE (optimizer%x) 346 END IF 347 IF (ASSOCIATED(optimizer%lower_bound)) THEN 348 DEALLOCATE (optimizer%lower_bound) 349 END IF 350 IF (ASSOCIATED(optimizer%upper_bound)) THEN 351 DEALLOCATE (optimizer%upper_bound) 352 END IF 353 IF (ASSOCIATED(optimizer%gradient)) THEN 354 DEALLOCATE (optimizer%gradient) 355 END IF 356 IF (ASSOCIATED(optimizer%dsave)) THEN 357 DEALLOCATE (optimizer%dsave) 358 END IF 359 IF (ASSOCIATED(optimizer%work_array)) THEN 360 DEALLOCATE (optimizer%work_array) 361 END IF 362 CALL cp_para_env_release(optimizer%para_env) 363 CALL gopt_f_release(optimizer%obj_funct) 364 NULLIFY (optimizer%obj_funct) 365 DEALLOCATE (optimizer) 366 END IF 367 END IF 368 NULLIFY (optimizer) 369 CALL timestop(handle) 370 END SUBROUTINE cp_opt_gopt_release 371 372! ************************************************************************************************** 373!> \brief takes different valuse from the optimizer 374!> \param optimizer ... 375!> \param para_env ... 376!> \param obj_funct ... 377!> \param m ... 378!> \param print_every ... 379!> \param id_nr ... 380!> \param wanted_relative_f_delta ... 381!> \param wanted_projected_gradient ... 382!> \param x ... 383!> \param lower_bound ... 384!> \param upper_bound ... 385!> \param kind_of_bound ... 386!> \param master ... 387!> \param actual_projected_gradient ... 388!> \param n_var ... 389!> \param n_iter ... 390!> \param status ... 391!> \param max_f_per_iter ... 392!> \param at_end ... 393!> \param is_master ... 394!> \param last_f ... 395!> \param f ... 396!> \par History 397!> none 398!> \author Fawzi Mohamed 399!> @version 2.2002 400! ************************************************************************************************** 401 SUBROUTINE cp_opt_gopt_get(optimizer, para_env, & 402 obj_funct, m, print_every, id_nr, & 403 wanted_relative_f_delta, wanted_projected_gradient, & 404 x, lower_bound, upper_bound, kind_of_bound, master, & 405 actual_projected_gradient, & 406 n_var, n_iter, status, max_f_per_iter, at_end, & 407 is_master, last_f, f) 408 TYPE(cp_lbfgs_opt_gopt_type), POINTER :: optimizer 409 TYPE(cp_para_env_type), OPTIONAL, POINTER :: para_env 410 TYPE(gopt_f_type), OPTIONAL, POINTER :: obj_funct 411 INTEGER, INTENT(out), OPTIONAL :: m, print_every, id_nr 412 REAL(kind=dp), INTENT(out), OPTIONAL :: wanted_relative_f_delta, & 413 wanted_projected_gradient 414 REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: x, lower_bound, upper_bound 415 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: kind_of_bound 416 INTEGER, INTENT(out), OPTIONAL :: master 417 REAL(kind=dp), INTENT(out), OPTIONAL :: actual_projected_gradient 418 INTEGER, INTENT(out), OPTIONAL :: n_var, n_iter, status, max_f_per_iter 419 LOGICAL, INTENT(out), OPTIONAL :: at_end, is_master 420 REAL(kind=dp), INTENT(out), OPTIONAL :: last_f, f 421 422 CHARACTER(len=*), PARAMETER :: routineN = 'cp_opt_gopt_get', & 423 routineP = moduleN//':'//routineN 424 425! call timeset(routineN,handle) 426 427 CPASSERT(ASSOCIATED(optimizer)) 428 CPASSERT(optimizer%ref_count > 0) 429 430 IF (PRESENT(is_master)) is_master = optimizer%master == optimizer%para_env%mepos 431 IF (PRESENT(master)) master = optimizer%master 432 IF (PRESENT(status)) status = optimizer%status 433 IF (PRESENT(para_env)) para_env => optimizer%para_env 434 IF (PRESENT(id_nr)) id_nr = optimizer%id_nr 435 IF (PRESENT(obj_funct)) obj_funct = optimizer%obj_funct 436 IF (PRESENT(m)) m = optimizer%m 437 IF (PRESENT(max_f_per_iter)) max_f_per_iter = optimizer%max_f_per_iter 438 IF (PRESENT(wanted_projected_gradient)) & 439 wanted_projected_gradient = optimizer%wanted_projected_gradient 440 IF (PRESENT(wanted_relative_f_delta)) & 441 wanted_relative_f_delta = optimizer%wanted_relative_f_delta 442 IF (PRESENT(print_every)) print_every = optimizer%print_every 443 IF (PRESENT(x)) x => optimizer%x 444 IF (PRESENT(n_var)) n_var = SIZE(x) 445 IF (PRESENT(lower_bound)) lower_bound => optimizer%lower_bound 446 IF (PRESENT(upper_bound)) upper_bound => optimizer%upper_bound 447 IF (PRESENT(kind_of_bound)) kind_of_bound => optimizer%kind_of_bound 448 IF (PRESENT(n_iter)) n_iter = optimizer%n_iter 449 IF (PRESENT(last_f)) last_f = optimizer%last_f 450 IF (PRESENT(f)) f = optimizer%f 451 IF (PRESENT(at_end)) at_end = optimizer%status > 3 452 IF (PRESENT(actual_projected_gradient)) & 453 actual_projected_gradient = optimizer%projected_gradient 454 IF (optimizer%master == optimizer%para_env%mepos) THEN 455 IF (optimizer%isave(30) > 1 .AND. (optimizer%task(1:5) == "NEW_X" .OR. & 456 optimizer%task(1:4) == "STOP" .AND. optimizer%task(7:9) == "CPU")) THEN 457 ! nr iterations >1 .and. dsave contains the wanted data 458 IF (PRESENT(last_f)) last_f = optimizer%dsave(2) 459 IF (PRESENT(actual_projected_gradient)) & 460 actual_projected_gradient = optimizer%dsave(13) 461 ELSE 462 CPASSERT(.NOT. PRESENT(last_f)) 463 CPASSERT(.NOT. PRESENT(actual_projected_gradient)) 464 END IF 465 ELSE 466 IF (PRESENT(lower_bound) .OR. PRESENT(upper_bound) .OR. PRESENT(kind_of_bound)) & 467 CPWARN("asked undefined types") 468 END IF 469 470 ! call timestop(handle) 471 END SUBROUTINE cp_opt_gopt_get 472 473! ************************************************************************************************** 474!> \brief does one optimization step 475!> \param optimizer ... 476!> \param n_iter ... 477!> \param f ... 478!> \param last_f ... 479!> \param projected_gradient ... 480!> \param converged ... 481!> \param geo_section ... 482!> \param force_env ... 483!> \param gopt_param ... 484!> \par History 485!> none 486!> \author Fawzi Mohamed 487!> @version 2.2002 488!> \note 489!> use directly mainlb in place of setulb ?? 490! ************************************************************************************************** 491 RECURSIVE SUBROUTINE cp_opt_gopt_step(optimizer, n_iter, f, last_f, & 492 projected_gradient, converged, geo_section, force_env, & 493 gopt_param) 494 TYPE(cp_lbfgs_opt_gopt_type), POINTER :: optimizer 495 INTEGER, INTENT(out), OPTIONAL :: n_iter 496 REAL(kind=dp), INTENT(out), OPTIONAL :: f, last_f, projected_gradient 497 LOGICAL, INTENT(out), OPTIONAL :: converged 498 TYPE(section_vals_type), POINTER :: geo_section 499 TYPE(force_env_type), POINTER :: force_env 500 TYPE(gopt_param_type), POINTER :: gopt_param 501 502 CHARACTER(len=*), PARAMETER :: routineN = 'cp_opt_gopt_step', & 503 routineP = moduleN//':'//routineN 504 505 CHARACTER(LEN=5) :: wildcard 506 INTEGER :: dataunit, handle, its 507 LOGICAL :: conv, is_master, justEntred 508 REAL(KIND=dp) :: t_diff, t_now, t_old 509 REAL(KIND=dp), DIMENSION(:), POINTER :: xold 510 TYPE(cp_logger_type), POINTER :: logger 511 512 NULLIFY (logger, xold) 513 logger => cp_get_default_logger() 514 CALL timeset(routineN, handle) 515 justEntred = .TRUE. 516 is_master = optimizer%master == optimizer%para_env%mepos 517 IF (PRESENT(converged)) converged = optimizer%status == 4 518 ALLOCATE (xold(SIZE(optimizer%x))) 519 xold = optimizer%x 520 t_old = m_walltime() 521 522 CPASSERT(ASSOCIATED(optimizer)) 523 CPASSERT(optimizer%ref_count > 0) 524 525 IF (optimizer%status >= 4) THEN 526 CPWARN("status>=4, trying to restart") 527 optimizer%status = 0 528 IF (is_master) THEN 529 optimizer%task = 'START' 530 CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, & 531 optimizer%lower_bound, optimizer%upper_bound, & 532 optimizer%kind_of_bound, optimizer%f, optimizer%gradient, & 533 optimizer%wanted_relative_f_delta, & 534 optimizer%wanted_projected_gradient, optimizer%work_array, & 535 optimizer%i_work_array, optimizer%task, optimizer%print_every, & 536 optimizer%csave, optimizer%lsave, optimizer%isave, & 537 optimizer%dsave, optimizer%trust_radius) 538 END IF 539 END IF 540 541 DO 542 ifMaster: IF (is_master) THEN 543 IF (optimizer%task(1:7) == 'RESTART') THEN 544 ! restart the optimizer 545 optimizer%status = 0 546 optimizer%task = 'START' 547 CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, & 548 optimizer%lower_bound, optimizer%upper_bound, & 549 optimizer%kind_of_bound, optimizer%f, optimizer%gradient, & 550 optimizer%wanted_relative_f_delta, & 551 optimizer%wanted_projected_gradient, optimizer%work_array, & 552 optimizer%i_work_array, optimizer%task, optimizer%print_every, & 553 optimizer%csave, optimizer%lsave, optimizer%isave, & 554 optimizer%dsave, optimizer%trust_radius) 555 END IF 556 IF (optimizer%task(1:2) == 'FG') THEN 557 IF (optimizer%isave(36) > optimizer%max_f_per_iter) THEN 558 optimizer%task = 'STOP: CPU, hit max f eval in iter' 559 optimizer%status = 5 ! anormal exit 560 CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, & 561 optimizer%lower_bound, optimizer%upper_bound, & 562 optimizer%kind_of_bound, optimizer%f, optimizer%gradient, & 563 optimizer%wanted_relative_f_delta, & 564 optimizer%wanted_projected_gradient, optimizer%work_array, & 565 optimizer%i_work_array, optimizer%task, optimizer%print_every, & 566 optimizer%csave, optimizer%lsave, optimizer%isave, & 567 optimizer%dsave, optimizer%trust_radius) 568 ELSE 569 optimizer%status = 1 570 END IF 571 ELSE IF (optimizer%task(1:5) == 'NEW_X') THEN 572 IF (justEntred) THEN 573 optimizer%status = 2 574 CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, & 575 optimizer%lower_bound, optimizer%upper_bound, & 576 optimizer%kind_of_bound, optimizer%f, optimizer%gradient, & 577 optimizer%wanted_relative_f_delta, & 578 optimizer%wanted_projected_gradient, optimizer%work_array, & 579 optimizer%i_work_array, optimizer%task, optimizer%print_every, & 580 optimizer%csave, optimizer%lsave, optimizer%isave, & 581 optimizer%dsave, optimizer%trust_radius) 582 ELSE 583 optimizer%status = 3 584 END IF 585 ELSE IF (optimizer%task(1:4) == 'CONV') THEN 586 optimizer%status = 4 587 ELSE IF (optimizer%task(1:4) == 'STOP') THEN 588 optimizer%status = 5 589 CPWARN("task became stop in an unknown way") 590 ELSE IF (optimizer%task(1:5) == 'ERROR') THEN 591 optimizer%status = 5 592 ELSE 593 CPWARN("unknown task '"//optimizer%task//"'") 594 END IF 595 END IF ifMaster 596 CALL mp_bcast(optimizer%status, optimizer%master, optimizer%para_env%group) 597 ! Dump info 598 IF (optimizer%status == 3) THEN 599 its = 0 600 IF (is_master) THEN 601 ! Iteration level is taken into account in the optimizer external loop 602 its = optimizer%isave(30) 603 END IF 604 END IF 605 ! 606 SELECT CASE (optimizer%status) 607 CASE (1) 608 !op=1 evaluate f and g 609 CALL cp_eval_at(optimizer%obj_funct, x=optimizer%x, & 610 f=optimizer%f, & 611 gradient=optimizer%gradient, & 612 final_evaluation=.FALSE., & 613 master=optimizer%master, para_env=optimizer%para_env) 614 ! do not use keywords? 615 IF (is_master) THEN 616 CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, & 617 optimizer%lower_bound, optimizer%upper_bound, & 618 optimizer%kind_of_bound, optimizer%f, optimizer%gradient, & 619 optimizer%wanted_relative_f_delta, & 620 optimizer%wanted_projected_gradient, optimizer%work_array, & 621 optimizer%i_work_array, optimizer%task, optimizer%print_every, & 622 optimizer%csave, optimizer%lsave, optimizer%isave, & 623 optimizer%dsave, optimizer%trust_radius) 624 END IF 625 CALL mp_bcast(optimizer%x, optimizer%master, optimizer%para_env%group) 626 CASE (2) 627 !op=2 begin new iter 628 CALL mp_bcast(optimizer%x, optimizer%master, optimizer%para_env%group) 629 t_old = m_walltime() 630 CASE (3) 631 !op=3 ended iter 632 wildcard = "LBFGS" 633 dataunit = cp_print_key_unit_nr(logger, geo_section, & 634 "PRINT%PROGRAM_RUN_INFO", extension=".geoLog") 635 IF (is_master) its = optimizer%isave(30) 636 CALL mp_bcast(its, optimizer%master, optimizer%para_env%group) 637 638 ! Some IO and Convergence check 639 t_now = m_walltime() 640 t_diff = t_now - t_old 641 t_old = t_now 642 CALL gopt_f_io(optimizer%obj_funct, force_env, force_env%root_section, & 643 its, optimizer%f, dataunit, optimizer%eold, optimizer%emin, wildcard, gopt_param, & 644 SIZE(optimizer%x), optimizer%x - xold, optimizer%gradient, conv, used_time=t_diff) 645 CALL mp_bcast(conv, optimizer%master, optimizer%para_env%group) 646 CALL cp_print_key_finished_output(dataunit, logger, geo_section, & 647 "PRINT%PROGRAM_RUN_INFO") 648 optimizer%eold = optimizer%f 649 optimizer%emin = MIN(optimizer%emin, optimizer%eold) 650 xold = optimizer%x 651 IF (PRESENT(converged)) converged = conv 652 EXIT 653 CASE (4) 654 !op=4 (convergence - normal exit) 655 ! Specific L-BFGS convergence criteria.. overrides the convergence criteria on 656 ! stepsize and gradients 657 dataunit = cp_print_key_unit_nr(logger, geo_section, & 658 "PRINT%PROGRAM_RUN_INFO", extension=".geoLog") 659 IF (dataunit > 0) THEN 660 WRITE (dataunit, '(T2,A)') "" 661 WRITE (dataunit, '(T2,A)') "***********************************************" 662 WRITE (dataunit, '(T2,A)') "* Specific L-BFGS convergence criteria " 663 WRITE (dataunit, '(T2,A)') "* WANTED_PROJ_GRADIENT and WANTED_REL_F_ERROR " 664 WRITE (dataunit, '(T2,A)') "* satisfied .... run CONVERGED! " 665 WRITE (dataunit, '(T2,A)') "***********************************************" 666 WRITE (dataunit, '(T2,A)') "" 667 END IF 668 CALL cp_print_key_finished_output(dataunit, logger, geo_section, & 669 "PRINT%PROGRAM_RUN_INFO") 670 IF (PRESENT(converged)) converged = .TRUE. 671 EXIT 672 CASE (5) 673 ! op=5 abnormal exit () 674 CALL mp_bcast(optimizer%task, optimizer%master, & 675 optimizer%para_env%group) 676 CASE (6) 677 ! deallocated 678 CPABORT("step on a deallocated opt structure ") 679 CASE default 680 CALL cp_abort(__LOCATION__, & 681 "unknown status "//cp_to_string(optimizer%status)) 682 optimizer%status = 5 683 EXIT 684 END SELECT 685 IF (optimizer%status == 1 .AND. justEntred) THEN 686 optimizer%eold = optimizer%f 687 optimizer%emin = optimizer%eold 688 END IF 689 justEntred = .FALSE. 690 END DO 691 CALL mp_bcast(optimizer%x, optimizer%master, & 692 optimizer%para_env%group) 693 CALL cp_opt_gopt_bcast_res(optimizer, & 694 n_iter=optimizer%n_iter, & 695 f=optimizer%f, last_f=optimizer%last_f, & 696 projected_gradient=optimizer%projected_gradient) 697 698 DEALLOCATE (xold) 699 IF (PRESENT(f)) f = optimizer%f 700 IF (PRESENT(last_f)) last_f = optimizer%last_f 701 IF (PRESENT(projected_gradient)) & 702 projected_gradient = optimizer%projected_gradient 703 IF (PRESENT(n_iter)) n_iter = optimizer%n_iter 704 CALL timestop(handle) 705 706 END SUBROUTINE cp_opt_gopt_step 707 708! ************************************************************************************************** 709!> \brief returns the results (and broadcasts them) 710!> \param optimizer the optimizer object the info is taken from 711!> \param n_iter the number of iterations 712!> \param f the actual value of the objective function (f) 713!> \param last_f the last value of f 714!> \param projected_gradient the infinity norm of the projected gradient 715!> \par History 716!> none 717!> \author Fawzi Mohamed 718!> @version 2.2002 719!> \note 720!> private routine 721! ************************************************************************************************** 722 SUBROUTINE cp_opt_gopt_bcast_res(optimizer, n_iter, f, last_f, & 723 projected_gradient) 724 TYPE(cp_lbfgs_opt_gopt_type), POINTER :: optimizer 725 INTEGER, INTENT(out), OPTIONAL :: n_iter 726 REAL(kind=dp), INTENT(inout), OPTIONAL :: f, last_f, projected_gradient 727 728 CHARACTER(len=*), PARAMETER :: routineN = 'cp_opt_gopt_bcast_res', & 729 routineP = moduleN//':'//routineN 730 731 REAL(kind=dp), DIMENSION(4) :: results 732 733! call timeset(routineN,handle) 734 735 CPASSERT(ASSOCIATED(optimizer)) 736 CPASSERT(optimizer%ref_count > 0) 737 738 IF (optimizer%master == optimizer%para_env%mepos) THEN 739 results = (/REAL(optimizer%isave(30), kind=dp), & 740 optimizer%f, optimizer%dsave(2), optimizer%dsave(13)/) 741 END IF 742 CALL mp_bcast(results, optimizer%master, & 743 optimizer%para_env%group) 744 IF (PRESENT(n_iter)) n_iter = NINT(results(1)) 745 IF (PRESENT(f)) f = results(2) 746 IF (PRESENT(last_f)) last_f = results(3) 747 IF (PRESENT(projected_gradient)) & 748 projected_gradient = results(4) 749 750 ! call timestop(handle) 751 END SUBROUTINE cp_opt_gopt_bcast_res 752 753! ************************************************************************************************** 754!> \brief goes to the next optimal point (after an optimizer iteration) 755!> returns true if converged 756!> \param optimizer the optimizer that goes to the next point 757!> \param n_iter ... 758!> \param f ... 759!> \param last_f ... 760!> \param projected_gradient ... 761!> \param converged ... 762!> \param geo_section ... 763!> \param force_env ... 764!> \param gopt_param ... 765!> \return ... 766!> \par History 767!> none 768!> \author Fawzi Mohamed 769!> @version 2.2002 770!> \note 771!> if you deactivate convergence control it returns never false 772! ************************************************************************************************** 773 RECURSIVE FUNCTION cp_opt_gopt_next(optimizer, n_iter, f, last_f, & 774 projected_gradient, converged, geo_section, force_env, & 775 gopt_param) RESULT(res) 776 TYPE(cp_lbfgs_opt_gopt_type), POINTER :: optimizer 777 INTEGER, INTENT(out), OPTIONAL :: n_iter 778 REAL(kind=dp), INTENT(out), OPTIONAL :: f, last_f, projected_gradient 779 LOGICAL, INTENT(out) :: converged 780 TYPE(section_vals_type), POINTER :: geo_section 781 TYPE(force_env_type), POINTER :: force_env 782 TYPE(gopt_param_type), POINTER :: gopt_param 783 LOGICAL :: res 784 785 CHARACTER(len=*), PARAMETER :: routineN = 'cp_opt_gopt_next', & 786 routineP = moduleN//':'//routineN 787 788!call timeset(routineN,handle) 789 790 CPASSERT(ASSOCIATED(optimizer)) 791 CPASSERT(optimizer%ref_count > 0) 792 CALL cp_opt_gopt_step(optimizer, n_iter=n_iter, f=f, & 793 last_f=last_f, projected_gradient=projected_gradient, & 794 converged=converged, geo_section=geo_section, & 795 force_env=force_env, gopt_param=gopt_param) 796 res = (optimizer%status < 40) .AND. .NOT. converged 797 !call timestop(handle) 798 END FUNCTION cp_opt_gopt_next 799 800! ************************************************************************************************** 801!> \brief stops the optimization 802!> \param optimizer ... 803!> \par History 804!> none 805!> \author Fawzi Mohamed 806!> @version 2.2002 807!> \note 808!> necessary??? 809! ************************************************************************************************** 810 SUBROUTINE cp_opt_gopt_stop(optimizer) 811 TYPE(cp_lbfgs_opt_gopt_type), POINTER :: optimizer 812 813 CHARACTER(len=*), PARAMETER :: routineN = 'cp_opt_gopt_stop', & 814 routineP = moduleN//':'//routineN 815 816! call timeset(routineN,handle) 817 818 CPASSERT(ASSOCIATED(optimizer)) 819 CPASSERT(optimizer%ref_count > 0) 820 821 optimizer%task = 'STOPPED on user request' 822 optimizer%status = 4 ! normal exit 823 IF (optimizer%master == optimizer%para_env%mepos) THEN 824 CALL setulb(SIZE(optimizer%x), optimizer%m, optimizer%x, & 825 optimizer%lower_bound, optimizer%upper_bound, & 826 optimizer%kind_of_bound, optimizer%f, optimizer%gradient, & 827 optimizer%wanted_relative_f_delta, & 828 optimizer%wanted_projected_gradient, optimizer%work_array, & 829 optimizer%i_work_array, optimizer%task, optimizer%print_every, & 830 optimizer%csave, optimizer%lsave, optimizer%isave, & 831 optimizer%dsave, optimizer%trust_radius) 832 END IF 833 834 ! call timestop(handle) 835 END SUBROUTINE cp_opt_gopt_stop 836 837 ! template def put here so that line numbers in template and derived 838 ! files are almost the same (multi-line use change it a bit) 839 ! [template(type1,nametype1,USE,type1_retain,type1_release,type1_eval_at)] 840! ARGS: 841! INTERFACE = "#include "gopt_f77_methods.h"" 842! USE = "USE gopt_f_types, ONLY: gopt_f_type, gopt_f_release,gopt_f_retain" 843! nametype1 = "gopt" 844! type1 = "type(gopt_f_type)" 845! type1_eval_at = "cp_eval_at" 846! type1_release = "gopt_f_release" 847! type1_retain = "gopt_f_retain" 848 849END MODULE cp_lbfgs_optimizer_gopt 850