From 0325e59ebef7fa8c0c6757d126215ee28db3021d Mon Sep 17 00:00:00 2001 From: ydamour Date: Tue, 18 Apr 2023 11:22:04 +0200 Subject: [PATCH] remove old utils_trust_region --- src/utils_trust_region/EZFIO.cfg | 89 - src/utils_trust_region/NEED | 1 - src/utils_trust_region/README.rst | 5 - src/utils_trust_region/TANGLE_org_mode.sh | 7 - src/utils_trust_region/algo_trust.irp.f | 248 --- src/utils_trust_region/algo_trust.org | 593 ------ .../apply_mo_rotation.irp.f | 85 - src/utils_trust_region/apply_mo_rotation.org | 86 - src/utils_trust_region/mat_to_vec_index.irp.f | 61 - src/utils_trust_region/mat_to_vec_index.org | 63 - src/utils_trust_region/pi.h | 2 - src/utils_trust_region/rotation_matrix.irp.f | 443 ----- src/utils_trust_region/rotation_matrix.org | 454 ----- .../sub_to_full_rotation_matrix.irp.f | 64 - .../sub_to_full_rotation_matrix.org | 65 - .../trust_region_expected_e.irp.f | 119 -- .../trust_region_expected_e.org | 121 -- .../trust_region_optimal_lambda.irp.f | 1655 ---------------- .../trust_region_optimal_lambda.org | 1665 ----------------- src/utils_trust_region/trust_region_rho.irp.f | 121 -- src/utils_trust_region/trust_region_rho.org | 123 -- .../trust_region_step.irp.f | 716 ------- src/utils_trust_region/trust_region_step.org | 726 ------- src/utils_trust_region/vec_to_mat_index.irp.f | 71 - src/utils_trust_region/vec_to_mat_index.org | 72 - src/utils_trust_region/vec_to_mat_v2.irp.f | 39 - src/utils_trust_region/vec_to_mat_v2.org | 40 - 27 files changed, 7734 deletions(-) delete mode 100644 src/utils_trust_region/EZFIO.cfg delete mode 100644 src/utils_trust_region/NEED delete mode 100644 src/utils_trust_region/README.rst delete mode 100755 src/utils_trust_region/TANGLE_org_mode.sh delete mode 100644 src/utils_trust_region/algo_trust.irp.f delete mode 100644 src/utils_trust_region/algo_trust.org delete mode 100644 src/utils_trust_region/apply_mo_rotation.irp.f delete mode 100644 src/utils_trust_region/apply_mo_rotation.org delete mode 100644 src/utils_trust_region/mat_to_vec_index.irp.f delete mode 100644 src/utils_trust_region/mat_to_vec_index.org delete mode 100644 src/utils_trust_region/pi.h delete mode 100644 src/utils_trust_region/rotation_matrix.irp.f delete mode 100644 src/utils_trust_region/rotation_matrix.org delete mode 100644 src/utils_trust_region/sub_to_full_rotation_matrix.irp.f delete mode 100644 src/utils_trust_region/sub_to_full_rotation_matrix.org delete mode 100644 src/utils_trust_region/trust_region_expected_e.irp.f delete mode 100644 src/utils_trust_region/trust_region_expected_e.org delete mode 100644 src/utils_trust_region/trust_region_optimal_lambda.irp.f delete mode 100644 src/utils_trust_region/trust_region_optimal_lambda.org delete mode 100644 src/utils_trust_region/trust_region_rho.irp.f delete mode 100644 src/utils_trust_region/trust_region_rho.org delete mode 100644 src/utils_trust_region/trust_region_step.irp.f delete mode 100644 src/utils_trust_region/trust_region_step.org delete mode 100644 src/utils_trust_region/vec_to_mat_index.irp.f delete mode 100644 src/utils_trust_region/vec_to_mat_index.org delete mode 100644 src/utils_trust_region/vec_to_mat_v2.irp.f delete mode 100644 src/utils_trust_region/vec_to_mat_v2.org diff --git a/src/utils_trust_region/EZFIO.cfg b/src/utils_trust_region/EZFIO.cfg deleted file mode 100644 index 9c9f6248..00000000 --- a/src/utils_trust_region/EZFIO.cfg +++ /dev/null @@ -1,89 +0,0 @@ -[thresh_delta] -type: double precision -doc: Threshold to stop the optimization if the radius of the trust region delta < thresh_delta -interface: ezfio,provider,ocaml -default: 1.e-10 - -[thresh_rho] -type: double precision -doc: Threshold for the step acceptance in the trust region algorithm, if (rho .geq. thresh_rho) the step is accepted, else the step is cancelled and a smaller step is tried until (rho .geq. thresh_rho) -interface: ezfio,provider,ocaml -default: 0.1 - -[thresh_eig] -type: double precision -doc: Threshold to consider when an eigenvalue is 0 in the trust region algorithm -interface: ezfio,provider,ocaml -default: 1.e-12 - -[thresh_model] -type: double precision -doc: If if ABS(criterion - criterion_model) < thresh_model, the program exit the trust region algorithm -interface: ezfio,provider,ocaml -default: 1.e-12 - -[absolute_eig] -type: logical -doc: If True, the algorithm replace the eigenvalues of the hessian by their absolute value to compute the step (in the trust region) -interface: ezfio,provider,ocaml -default: false - -[thresh_wtg] -type: double precision -doc: Threshold in the trust region algorithm to considere when the dot product of the eigenvector W by the gradient v_grad is equal to 0. Must be smaller than thresh_eig by several order of magnitude to avoid numerical problem. If the research of the optimal lambda cannot reach the condition (||x|| .eq. delta) because (||x|| .lt. delta), the reason might be that thresh_wtg is too big or/and thresh_eig is too small -interface: ezfio,provider,ocaml -default: 1.e-6 - -[thresh_wtg2] -type: double precision -doc: Threshold in the trust region algorithm to considere when the dot product of the eigenvector W by the gradient v_grad is 0 in the case of avoid_saddle .eq. true. There is no particular reason to put a different value that thresh_wtg, but it can be useful one day -interface: ezfio,provider,ocaml -default: 1.e-6 - -[avoid_saddle] -type: logical -doc: Test to avoid saddle point, active if true -interface: ezfio,provider,ocaml -default: false - -[version_avoid_saddle] -type: integer -doc: cf. trust region, not stable -interface: ezfio,provider,ocaml -default: 3 - -[thresh_rho_2] -type: double precision -doc: Threshold for the step acceptance for the research of lambda in the trust region algorithm, if (rho_2 .geq. thresh_rho_2) the step is accepted, else the step is rejected -interface: ezfio,provider,ocaml -default: 0.1 - -[thresh_cc] -type: double precision -doc: Threshold to stop the research of the optimal lambda in the trust region algorithm when (dabs(1d0-||x||^2/delta^2) < thresh_cc) -interface: ezfio,provider,ocaml -default: 1.e-6 - -[thresh_model_2] -type: double precision -doc: if (ABS(criterion - criterion_model) < thresh_model_2), i.e., the difference between the actual criterion and the predicted next criterion, during the research of the optimal lambda in the trust region algorithm it prints a warning -interface: ezfio,provider,ocaml -default: 1.e-12 - -[version_lambda_search] -type: integer -doc: Research of the optimal lambda in the trust region algorithm to constrain the norm of the step by solving: 1 -> ||x||^2 - delta^2 .eq. 0, 2 -> 1/||x||^2 - 1/delta^2 .eq. 0 -interface: ezfio,provider,ocaml -default: 2 - -[nb_it_max_lambda] -type: integer -doc: Maximal number of iterations for the research of the optimal lambda in the trust region algorithm -interface: ezfio,provider,ocaml -default: 100 - -[nb_it_max_pre_search] -type: integer -doc: Maximal number of iterations for the pre-research of the optimal lambda in the trust region algorithm -interface: ezfio,provider,ocaml -default: 40 diff --git a/src/utils_trust_region/NEED b/src/utils_trust_region/NEED deleted file mode 100644 index 1a65ce38..00000000 --- a/src/utils_trust_region/NEED +++ /dev/null @@ -1 +0,0 @@ -hartree_fock diff --git a/src/utils_trust_region/README.rst b/src/utils_trust_region/README.rst deleted file mode 100644 index 6a0689b6..00000000 --- a/src/utils_trust_region/README.rst +++ /dev/null @@ -1,5 +0,0 @@ -============ -trust_region -============ - -The documentation can be found in the org files. diff --git a/src/utils_trust_region/TANGLE_org_mode.sh b/src/utils_trust_region/TANGLE_org_mode.sh deleted file mode 100755 index 059cbe7d..00000000 --- a/src/utils_trust_region/TANGLE_org_mode.sh +++ /dev/null @@ -1,7 +0,0 @@ -#!/bin/sh - -list='ls *.org' -for element in $list -do - emacs --batch $element -f org-babel-tangle -done diff --git a/src/utils_trust_region/algo_trust.irp.f b/src/utils_trust_region/algo_trust.irp.f deleted file mode 100644 index eac17275..00000000 --- a/src/utils_trust_region/algo_trust.irp.f +++ /dev/null @@ -1,248 +0,0 @@ -! Algorithm for the trust region - -! step_in_trust_region: -! Computes the step in the trust region (delta) -! (automatically sets at the iteration 0 and which evolves during the -! process in function of the evolution of rho). The step is computing by -! constraining its norm with a lagrange multiplier. -! Since the calculation of the step is based on the Newton method, an -! estimation of the gain in energy is given using the Taylors series -! truncated at the second order (criterion_model). -! If (DABS(criterion-criterion_model) < 1d-12) then -! must_exit = .True. -! else -! must_exit = .False. - -! This estimation of the gain in energy is used by -! is_step_cancel_trust_region to say if the step is accepted or cancelled. - -! If the step must be cancelled, the calculation restart from the same -! hessian and gradient and recomputes the step but in a smaller trust -! region and so on until the step is accepted. If the step is accepted -! the hessian and the gradient are recomputed to produce a new step. - -! Example: - - -! !### Initialization ### -! delta = 0d0 -! nb_iter = 0 ! Must start at 0 !!! -! rho = 0.5d0 -! not_converged = .True. -! -! ! ### TODO ### -! ! Compute the criterion before the loop -! call #your_criterion(prev_criterion) -! -! do while (not_converged) -! ! ### TODO ## -! ! Call your gradient -! ! Call you hessian -! call #your_gradient(v_grad) (1D array) -! call #your_hessian(H) (2D array) -! -! ! ### TODO ### -! ! Diagonalization of the hessian -! call diagonalization_hessian(n,H,e_val,w) -! -! cancel_step = .True. ! To enter in the loop just after -! ! Loop to Reduce the trust radius until the criterion decreases and rho >= thresh_rho -! do while (cancel_step) -! -! ! Hessian,gradient,Criterion -> x -! call trust_region_step_w_expected_e(tmp_n,W,e_val,v_grad,prev_criterion,rho,nb_iter,delta,criterion_model,tmp_x,must_exit) -! -! if (must_exit) then -! ! ### Message ### -! ! if step_in_trust_region sets must_exit on true for numerical reasons -! print*,'algo_trust1 sends the message : Exit' -! !### exit ### -! endif -! -! !### TODO ### -! ! Compute x -> m_x -! ! Compute m_x -> R -! ! Apply R and keep the previous MOs... -! ! Update/touch -! ! Compute the new criterion/energy -> criterion -! -! call #your_routine_1D_to_2D_antisymmetric_array(x,m_x) -! call #your_routine_2D_antisymmetric_array_to_rotation_matrix(m_x,R) -! call #your_routine_to_apply_the_rotation_matrix(R,prev_mos) -! -! TOUCH #your_variables -! -! call #your_criterion(criterion) -! -! ! Criterion -> step accepted or rejected -! call trust_region_is_step_cancelled(nb_iter,prev_criterion, criterion, criterion_model,rho,cancel_step) -! -! ! ### TODO ### -! !if (cancel_step) then -! ! Cancel the previous step (mo_coef = prev_mos if you keep them...) -! !endif -! #if (cancel_step) then -! #mo_coef = prev_mos -! #endif -! -! enddo -! -! !call save_mos() !### depend of the time for 1 iteration -! -! ! To exit the external loop if must_exit = .True. -! if (must_exit) then -! !### exit ### -! endif -! -! ! Step accepted, nb iteration + 1 -! nb_iter = nb_iter + 1 -! -! ! ### TODO ### -! !if (###Conditions###) then -! ! no_converged = .False. -! !endif -! #if (#your_conditions) then -! # not_converged = .False. -! #endif -! -! enddo - - - -! Variables: - -! Input: -! | n | integer | m*(m-1)/2 | -! | m | integer | number of mo in the mo_class | -! | H(n,n) | double precision | Hessian | -! | v_grad(n) | double precision | Gradient | -! | W(n,n) | double precision | Eigenvectors of the hessian | -! | e_val(n) | double precision | Eigenvalues of the hessian | -! | criterion | double precision | Actual criterion | -! | prev_criterion | double precision | Value of the criterion before the first iteration/after the previous iteration | -! | rho | double precision | Given by is_step_cancel_trus_region | -! | | | Agreement between the real function and the Taylor series (2nd order) | -! | nb_iter | integer | Actual number of iterations | - -! Input/output: -! | delta | double precision | Radius of the trust region | - -! Output: -! | criterion_model | double precision | Predicted criterion after the rotation | -! | x(n) | double precision | Step | -! | must_exit | logical | If the program must exit the loop | - - -subroutine trust_region_step_w_expected_e(n,H,W,e_val,v_grad,prev_criterion,rho,nb_iter,delta,criterion_model,x,must_exit) - - include 'pi.h' - - BEGIN_DOC - ! Compute the step and the expected criterion/energy after the step - END_DOC - - implicit none - - ! in - integer, intent(in) :: n, nb_iter - double precision, intent(in) :: H(n,n), W(n,n), v_grad(n) - double precision, intent(in) :: rho, prev_criterion - - ! inout - double precision, intent(inout) :: delta, e_val(n) - - ! out - double precision, intent(out) :: criterion_model, x(n) - logical, intent(out) :: must_exit - - ! internal - integer :: info - - must_exit = .False. - - call trust_region_step(n,nb_iter,v_grad,rho,e_val,W,x,delta) - - call trust_region_expected_e(n,v_grad,H,x,prev_criterion,criterion_model) - - ! exit if DABS(prev_criterion - criterion_model) < 1d-12 - if (DABS(prev_criterion - criterion_model) < thresh_model) then - print*,'' - print*,'###############################################################################' - print*,'DABS(prev_criterion - criterion_model) <', thresh_model, 'stop the trust region' - print*,'###############################################################################' - print*,'' - must_exit = .True. - endif - - if (delta < thresh_delta) then - print*,'' - print*,'##############################################' - print*,'Delta <', thresh_delta, 'stop the trust region' - print*,'##############################################' - print*,'' - must_exit = .True. - endif - - ! Add after the call to this subroutine, a statement: - ! "if (must_exit) then - ! exit - ! endif" - ! in order to exit the optimization loop - -end subroutine - - - -! Variables: - -! Input: -! | nb_iter | integer | actual number of iterations | -! | prev_criterion | double precision | criterion before the application of the step x | -! | criterion | double precision | criterion after the application of the step x | -! | criterion_model | double precision | predicted criterion after the application of x | - -! Output: -! | rho | double precision | Agreement between the predicted criterion and the real new criterion | -! | cancel_step | logical | If the step must be cancelled | - - -subroutine trust_region_is_step_cancelled(nb_iter,prev_criterion, criterion, criterion_model,rho,cancel_step) - - include 'pi.h' - - BEGIN_DOC - ! Compute if the step should be cancelled - END_DOC - - implicit none - - ! in - double precision, intent(in) :: prev_criterion, criterion, criterion_model - - ! inout - integer, intent(inout) :: nb_iter - - ! out - logical, intent(out) :: cancel_step - double precision, intent(out) :: rho - - ! Computes rho - call trust_region_rho(prev_criterion,criterion,criterion_model,rho) - - if (nb_iter == 0) then - nb_iter = 1 ! in order to enable the change of delta if the first iteration is cancelled - endif - - ! If rho < thresh_rho -> give something in output to cancel the step - if (rho >= thresh_rho) then !0.1d0) then - ! The step is accepted - cancel_step = .False. - else - ! The step is rejected - cancel_step = .True. - print*, '***********************' - print*, 'Step cancel : rho <', thresh_rho - print*, '***********************' - endif - -end subroutine diff --git a/src/utils_trust_region/algo_trust.org b/src/utils_trust_region/algo_trust.org deleted file mode 100644 index aa836f98..00000000 --- a/src/utils_trust_region/algo_trust.org +++ /dev/null @@ -1,593 +0,0 @@ -* Algorithm for the trust region - -step_in_trust_region: -Computes the step in the trust region (delta) -(automatically sets at the iteration 0 and which evolves during the -process in function of the evolution of rho). The step is computing by -constraining its norm with a lagrange multiplier. -Since the calculation of the step is based on the Newton method, an -estimation of the gain in energy is given using the Taylors series -truncated at the second order (criterion_model). -If (DABS(criterion-criterion_model) < 1d-12) then - must_exit = .True. -else - must_exit = .False. - -This estimation of the gain in energy is used by -is_step_cancel_trust_region to say if the step is accepted or cancelled. - -If the step must be cancelled, the calculation restart from the same -hessian and gradient and recomputes the step but in a smaller trust -region and so on until the step is accepted. If the step is accepted -the hessian and the gradient are recomputed to produce a new step. - -Example: - -#+BEGIN_SRC f90 :comments org :tangle algo_trust.irp.f -! !### Initialization ### -! delta = 0d0 -! nb_iter = 0 ! Must start at 0 !!! -! rho = 0.5d0 -! not_converged = .True. -! -! ! ### TODO ### -! ! Compute the criterion before the loop -! call #your_criterion(prev_criterion) -! -! do while (not_converged) -! ! ### TODO ## -! ! Call your gradient -! ! Call you hessian -! call #your_gradient(v_grad) (1D array) -! call #your_hessian(H) (2D array) -! -! ! ### TODO ### -! ! Diagonalization of the hessian -! call diagonalization_hessian(n,H,e_val,w) -! -! cancel_step = .True. ! To enter in the loop just after -! ! Loop to Reduce the trust radius until the criterion decreases and rho >= thresh_rho -! do while (cancel_step) -! -! ! Hessian,gradient,Criterion -> x -! call trust_region_step_w_expected_e(tmp_n,W,e_val,v_grad,prev_criterion,rho,nb_iter,delta,criterion_model,tmp_x,must_exit) -! -! if (must_exit) then -! ! ### Message ### -! ! if step_in_trust_region sets must_exit on true for numerical reasons -! print*,'algo_trust1 sends the message : Exit' -! !### exit ### -! endif -! -! !### TODO ### -! ! Compute x -> m_x -! ! Compute m_x -> R -! ! Apply R and keep the previous MOs... -! ! Update/touch -! ! Compute the new criterion/energy -> criterion -! -! call #your_routine_1D_to_2D_antisymmetric_array(x,m_x) -! call #your_routine_2D_antisymmetric_array_to_rotation_matrix(m_x,R) -! call #your_routine_to_apply_the_rotation_matrix(R,prev_mos) -! -! TOUCH #your_variables -! -! call #your_criterion(criterion) -! -! ! Criterion -> step accepted or rejected -! call trust_region_is_step_cancelled(nb_iter,prev_criterion, criterion, criterion_model,rho,cancel_step) -! -! ! ### TODO ### -! !if (cancel_step) then -! ! Cancel the previous step (mo_coef = prev_mos if you keep them...) -! !endif -! #if (cancel_step) then -! #mo_coef = prev_mos -! #endif -! -! enddo -! -! !call save_mos() !### depend of the time for 1 iteration -! -! ! To exit the external loop if must_exit = .True. -! if (must_exit) then -! !### exit ### -! endif -! -! ! Step accepted, nb iteration + 1 -! nb_iter = nb_iter + 1 -! -! ! ### TODO ### -! !if (###Conditions###) then -! ! no_converged = .False. -! !endif -! #if (#your_conditions) then -! # not_converged = .False. -! #endif -! -! enddo -#+END_SRC - -Variables: - -Input: -| n | integer | m*(m-1)/2 | -| m | integer | number of mo in the mo_class | -| H(n,n) | double precision | Hessian | -| v_grad(n) | double precision | Gradient | -| W(n,n) | double precision | Eigenvectors of the hessian | -| e_val(n) | double precision | Eigenvalues of the hessian | -| criterion | double precision | Actual criterion | -| prev_criterion | double precision | Value of the criterion before the first iteration/after the previous iteration | -| rho | double precision | Given by is_step_cancel_trus_region | -| | | Agreement between the real function and the Taylor series (2nd order) | -| nb_iter | integer | Actual number of iterations | - -Input/output: -| delta | double precision | Radius of the trust region | - -Output: -| criterion_model | double precision | Predicted criterion after the rotation | -| x(n) | double precision | Step | -| must_exit | logical | If the program must exit the loop | - -#+BEGIN_SRC f90 :comments org :tangle algo_trust.irp.f -subroutine trust_region_step_w_expected_e(n,H,W,e_val,v_grad,prev_criterion,rho,nb_iter,delta,criterion_model,x,must_exit) - - include 'pi.h' - - BEGIN_DOC - ! Compute the step and the expected criterion/energy after the step - END_DOC - - implicit none - - ! in - integer, intent(in) :: n, nb_iter - double precision, intent(in) :: H(n,n), W(n,n), v_grad(n) - double precision, intent(in) :: rho, prev_criterion - - ! inout - double precision, intent(inout) :: delta, e_val(n) - - ! out - double precision, intent(out) :: criterion_model, x(n) - logical, intent(out) :: must_exit - - ! internal - integer :: info - - must_exit = .False. - - call trust_region_step(n,nb_iter,v_grad,rho,e_val,W,x,delta) - - call trust_region_expected_e(n,v_grad,H,x,prev_criterion,criterion_model) - - ! exit if DABS(prev_criterion - criterion_model) < 1d-12 - if (DABS(prev_criterion - criterion_model) < thresh_model) then - print*,'' - print*,'###############################################################################' - print*,'DABS(prev_criterion - criterion_model) <', thresh_model, 'stop the trust region' - print*,'###############################################################################' - print*,'' - must_exit = .True. - endif - - if (delta < thresh_delta) then - print*,'' - print*,'##############################################' - print*,'Delta <', thresh_delta, 'stop the trust region' - print*,'##############################################' - print*,'' - must_exit = .True. - endif - - ! Add after the call to this subroutine, a statement: - ! "if (must_exit) then - ! exit - ! endif" - ! in order to exit the optimization loop - -end subroutine -#+END_SRC - -Variables: - -Input: -| nb_iter | integer | actual number of iterations | -| prev_criterion | double precision | criterion before the application of the step x | -| criterion | double precision | criterion after the application of the step x | -| criterion_model | double precision | predicted criterion after the application of x | - -Output: -| rho | double precision | Agreement between the predicted criterion and the real new criterion | -| cancel_step | logical | If the step must be cancelled | - -#+BEGIN_SRC f90 :comments org :tangle algo_trust.irp.f -subroutine trust_region_is_step_cancelled(nb_iter,prev_criterion, criterion, criterion_model,rho,cancel_step) - - include 'pi.h' - - BEGIN_DOC - ! Compute if the step should be cancelled - END_DOC - - implicit none - - ! in - double precision, intent(in) :: prev_criterion, criterion, criterion_model - - ! inout - integer, intent(inout) :: nb_iter - - ! out - logical, intent(out) :: cancel_step - double precision, intent(out) :: rho - - ! Computes rho - call trust_region_rho(prev_criterion,criterion,criterion_model,rho) - - if (nb_iter == 0) then - nb_iter = 1 ! in order to enable the change of delta if the first iteration is cancelled - endif - - ! If rho < thresh_rho -> give something in output to cancel the step - if (rho >= thresh_rho) then !0.1d0) then - ! The step is accepted - cancel_step = .False. - else - ! The step is rejected - cancel_step = .True. - print*, '***********************' - print*, 'Step cancel : rho <', thresh_rho - print*, '***********************' - endif - -end subroutine -#+END_SRC - -** Template for MOs -#+BEGIN_SRC f90 :comments org :tangle trust_region_template_mos.txt -subroutine algo_trust_template(tmp_n, tmp_list_size, tmp_list) - - implicit none - - ! Variables - - ! In - integer, intent(in) :: tmp_n, tmp_list_size, tmp_list(tmp_list_size) - - ! Out - ! Rien ou un truc pour savoir si ça c'est bien passé - - ! Internal - double precision, allocatable :: e_val(:), W(:,:), tmp_R(:,:), R(:,:), tmp_x(:), tmp_m_x(:,:) - double precision, allocatable :: prev_mos(:,:) - double precision :: criterion, prev_criterion, criterion_model - double precision :: delta, rho - logical :: not_converged, cancel_step, must_exit, enforce_step_cancellation - integer :: nb_iter, info, nb_sub_iter - integer :: i,j,tmp_i,tmp_j - - allocate(W(tmp_n, tmp_n),e_val(tmp_n),tmp_x(tmp_n),tmp_m_x(tmp_list_size, tmp_list_size)) - allocate(tmp_R(tmp_list_size, tmp_list_size), R(mo_num, mo_num)) - allocate(prev_mos(ao_num, mo_num)) - - ! Provide the criterion, but unnecessary because it's done - ! automatically - PROVIDE C_PROVIDER H_PROVIDER g_PROVIDER cc_PROVIDER - - ! Initialization - delta = 0d0 - nb_iter = 0 ! Must start at 0 !!! - rho = 0.5d0 ! Must start at 0.5 - not_converged = .True. ! Must be true - - ! Compute the criterion before the loop - prev_criterion = C_PROVIDER - - do while (not_converged) - - print*,'' - print*,'******************' - print*,'Iteration', nb_iter - print*,'******************' - print*,'' - - ! The new hessian and gradient are computed at the end of the previous iteration - ! Diagonalization of the hessian - call diagonalization_hessian(tmp_n, H_PROVIDER, e_val, W) - - cancel_step = .True. ! To enter in the loop just after - nb_sub_iter = 0 - - ! Loop to Reduce the trust radius until the criterion decreases and rho >= thresh_rho - do while (cancel_step) - - print*,'-----------------------------' - print*,'Iteration:', nb_iter - print*,'Sub iteration:', nb_sub_iter - print*,'-----------------------------' - - ! Hessian,gradient,Criterion -> x - call trust_region_step_w_expected_e(tmp_n, H_PROVIDER, W, e_val, g_PROVIDER, & - prev_criterion, rho, nb_iter, delta, criterion_model, tmp_x, must_exit) - - if (must_exit) then - ! if step_in_trust_region sets must_exit on true for numerical reasons - print*,'trust_region_step_w_expected_e sent the message : Exit' - exit - endif - - ! 1D tmp -> 2D tmp - call vec_to_mat_v2(tmp_n, tmp_list_size, tmp_x, tmp_m_x) - - ! Rotation submatrix (square matrix tmp_list_size by tmp_list_size) - call rotation_matrix(tmp_m_x, tmp_list_size, tmp_R, tmp_list_size, tmp_list_size, info, enforce_step_cancellation) - - if (enforce_step_cancellation) then - print*, 'Forces the step cancellation, too large error in the rotation matrix' - rho = 0d0 - cycle - endif - - ! tmp_R to R, subspace to full space - call sub_to_full_rotation_matrix(tmp_list_size, tmp_list, tmp_R, R) - - ! Rotation of the MOs - call apply_mo_rotation(R, prev_mos) - - ! touch mo_coef - call clear_mo_map ! Only if you are using the bi-electronic integrals - ! mo_coef becomes valid - ! And avoid the recomputation of the providers which depend of mo_coef - TOUCH mo_coef C_PROVIDER H_PROVIDER g_PROVIDER cc_PROVIDER - - ! To update the other parameters if needed - call #update_parameters() - - ! To enforce the program to provide new criterion after the update - ! of the parameters - FREE C_PROVIDER - PROVIDE C_PROVIDER - criterion = C_PROVIDER - - ! Criterion -> step accepted or rejected - call trust_region_is_step_cancelled(nb_iter, prev_criterion, criterion, criterion_model, rho, cancel_step) - - ! Cancellation of the step ? - if (cancel_step) then - ! Replacement by the previous MOs - mo_coef = prev_mos - ! call save_mos() ! depends of the time for 1 iteration - - ! No need to clear_mo_map since we don't recompute the gradient and the hessian - ! mo_coef becomes valid - ! Avoid the recomputation of the providers which depend of mo_coef - TOUCH mo_coef H_PROVIDER g_PROVIDER C_PROVIDER cc_PROVIDER - else - ! The step is accepted: - ! criterion -> prev criterion - - ! The replacement "criterion -> prev criterion" is already done - ! in trust_region_rho, so if the criterion does not have a reason - ! to change, it will change nothing for the criterion and will - ! force the program to provide the new hessian, gradient and - ! convergence criterion for the next iteration. - ! But in the case of orbital optimization we diagonalize the CI - ! matrix after the "FREE" statement, so the criterion will change - - FREE C_PROVIDER H_PROVIDER g_PROVIDER cc_PROVIDER - PROVIDE C_PROVIDER H_PROVIDER g_PROVIDER cc_PROVIDER - prev_criterion = C_PROVIDER - - endif - - nb_sub_iter = nb_sub_iter + 1 - enddo - - ! call save_mos() ! depends of the time for 1 iteration - - ! To exit the external loop if must_exit = .True. - if (must_exit) then - exit - endif - - ! Step accepted, nb iteration + 1 - nb_iter = nb_iter + 1 - - ! Provide the convergence criterion - ! Provide the gradient and the hessian for the next iteration - PROVIDE cc_PROVIDER - - ! To exit - if (dabs(cc_PROVIDER) < thresh_opt_max_elem_grad) then - not_converged = .False. - endif - - if (nb_iter > optimization_max_nb_iter) then - not_converged = .False. - endif - - if (delta < thresh_delta) then - not_converged = .False. - endif - - enddo - - ! Save the final MOs - call save_mos() - - ! Diagonalization of the hessian - ! (To see the eigenvalues at the end of the optimization) - call diagonalization_hessian(tmp_n, H_PROVIDER, e_val, W) - - deallocate(e_val, W, tmp_R, R, tmp_x, prev_mos) - -end -#+END_SRC - -** Cartesian version -#+BEGIN_SRC f90 :comments org :tangle trust_region_template_xyz.txt -subroutine algo_trust_cartesian_template(tmp_n) - - implicit none - - ! Variables - - ! In - integer, intent(in) :: tmp_n - - ! Out - ! Rien ou un truc pour savoir si ça c'est bien passé - - ! Internal - double precision, allocatable :: e_val(:), W(:,:), tmp_x(:) - double precision :: criterion, prev_criterion, criterion_model - double precision :: delta, rho - logical :: not_converged, cancel_step, must_exit - integer :: nb_iter, nb_sub_iter - integer :: i,j - - allocate(W(tmp_n, tmp_n),e_val(tmp_n),tmp_x(tmp_n)) - - PROVIDE C_PROVIDER X_PROVIDER H_PROVIDER g_PROVIDER - - ! Initialization - delta = 0d0 - nb_iter = 0 ! Must start at 0 !!! - rho = 0.5d0 ! Must start at 0.5 - not_converged = .True. ! Must be true - - ! Compute the criterion before the loop - prev_criterion = C_PROVIDER - - do while (not_converged) - - print*,'' - print*,'******************' - print*,'Iteration', nb_iter - print*,'******************' - print*,'' - - if (nb_iter > 0) then - PROVIDE H_PROVIDER g_PROVIDER - endif - - ! Diagonalization of the hessian - call diagonalization_hessian(tmp_n, H_PROVIDER, e_val, W) - - cancel_step = .True. ! To enter in the loop just after - nb_sub_iter = 0 - - ! Loop to Reduce the trust radius until the criterion decreases and rho >= thresh_rho - do while (cancel_step) - - print*,'-----------------------------' - print*,'Iteration:', nb_iter - print*,'Sub iteration:', nb_sub_iter - print*,'-----------------------------' - - ! Hessian,gradient,Criterion -> x - call trust_region_step_w_expected_e(tmp_n, H_PROVIDER, W, e_val, g_PROVIDER, & - prev_criterion, rho, nb_iter, delta, criterion_model, tmp_x, must_exit) - - if (must_exit) then - ! if step_in_trust_region sets must_exit on true for numerical reasons - print*,'trust_region_step_w_expected_e sent the message : Exit' - exit - endif - - ! New coordinates, check the sign - X_PROVIDER = X_PROVIDER - tmp_x - - ! touch X_PROVIDER - TOUCH X_PROVIDER H_PROVIDER g_PROVIDER cc_PROVIDER - - ! To update the other parameters if needed - call #update_parameters() - - ! New criterion - PROVIDE C_PROVIDER ! Unnecessary - criterion = C_PROVIDER - - ! Criterion -> step accepted or rejected - call trust_region_is_step_cancelled(nb_iter, prev_criterion, criterion, criterion_model, rho, cancel_step) - - ! Cancel the previous step - if (cancel_step) then - ! Replacement by the previous coordinates, check the sign - X_PROVIDER = X_PROVIDER + tmp_x - - ! Avoid the recomputation of the hessian and the gradient - TOUCH X_PROVIDER H_PROVIDER g_PROVIDER C_PROVIDER cc_PROVIDER - endif - - nb_sub_iter = nb_sub_iter + 1 - enddo - - ! To exit the external loop if must_exit = .True. - if (must_exit) then - exit - endif - - ! Step accepted, nb iteration + 1 - nb_iter = nb_iter + 1 - - PROVIDE cc_PROVIDER - - ! To exit - if (dabs(cc_PROVIDER) < thresh_opt_max_elem_grad) then - not_converged = .False. - endif - - if (nb_iter > optimization_max_nb_iter) then - not_converged = .False. - endif - - if (delta < thresh_delta) then - not_converged = .False. - endif - - enddo - - deallocate(e_val, W, tmp_x) - -end -#+END_SRC - -** Script template -#+BEGIN_SRC bash :tangle script_template_mos.sh -#!/bin/bash - -your_file= - -your_C_PROVIDER= -your_H_PROVIDER= -your_g_PROVIDER= -your_cc_PROVIDER= - -sed "s/C_PROVIDER/$your_C_PROVIDER/g" trust_region_template_mos.txt > $your_file -sed -i "s/H_PROVIDER/$your_H_PROVIDER/g" $your_file -sed -i "s/g_PROVIDER/$your_g_PROVIDER/g" $your_file -sed -i "s/cc_PROVIDER/$your_cc_PROVIDER/g" $your_file -#+END_SRC - -#+BEGIN_SRC bash :tangle script_template_xyz.sh -#!/bin/bash - -your_file= - -your_C_PROVIDER= -your_X_PROVIDER= -your_H_PROVIDER= -your_g_PROVIDER= -your_cc_PROVIDER= - -sed "s/C_PROVIDER/$your_C_PROVIDER/g" trust_region_template_xyz.txt > $your_file -sed -i "s/X_PROVIDER/$your_X_PROVIDER/g" $your_file -sed -i "s/H_PROVIDER/$your_H_PROVIDER/g" $your_file -sed -i "s/g_PROVIDER/$your_g_PROVIDER/g" $your_file -sed -i "s/cc_PROVIDER/$your_cc_PROVIDER/g" $your_file -#+END_SRC - diff --git a/src/utils_trust_region/apply_mo_rotation.irp.f b/src/utils_trust_region/apply_mo_rotation.irp.f deleted file mode 100644 index e274ec11..00000000 --- a/src/utils_trust_region/apply_mo_rotation.irp.f +++ /dev/null @@ -1,85 +0,0 @@ -! Apply MO rotation -! Subroutine to apply the rotation matrix to the coefficients of the -! MOs. - -! New MOs = Old MOs . Rotation matrix - -! *Compute the new MOs with the previous MOs and a rotation matrix* - -! Provided: -! | mo_num | integer | number of MOs | -! | ao_num | integer | number of AOs | -! | mo_coef(ao_num,mo_num) | double precision | coefficients of the MOs | - -! Intent in: -! | R(mo_num,mo_num) | double precision | rotation matrix | - -! Intent out: -! | prev_mos(ao_num,mo_num) | double precision | MOs before the rotation | - -! Internal: -! | new_mos(ao_num,mo_num) | double precision | MOs after the rotation | -! | i,j | integer | indexes | - -subroutine apply_mo_rotation(R,prev_mos) - - include 'pi.h' - - BEGIN_DOC - ! Compute the new MOs knowing the rotation matrix - END_DOC - - implicit none - - ! Variables - - ! in - double precision, intent(in) :: R(mo_num,mo_num) - - ! out - double precision, intent(out) :: prev_mos(ao_num,mo_num) - - ! internal - double precision, allocatable :: new_mos(:,:) - integer :: i,j - double precision :: t1,t2,t3 - - print*,'' - print*,'---apply_mo_rotation---' - - call wall_time(t1) - - ! Allocation - allocate(new_mos(ao_num,mo_num)) - - ! Calculation - - ! Product of old MOs (mo_coef) by Rotation matrix (R) - call dgemm('N','N',ao_num,mo_num,mo_num,1d0,mo_coef,size(mo_coef,1),R,size(R,1),0d0,new_mos,size(new_mos,1)) - - prev_mos = mo_coef - mo_coef = new_mos - - !if (debug) then - ! print*,'New mo_coef : ' - ! do i = 1, mo_num - ! write(*,'(100(F10.5))') mo_coef(i,:) - ! enddo - !endif - - ! Save the new MOs and change the label - mo_label = 'MCSCF' - !call save_mos - call ezfio_set_determinants_mo_label(mo_label) - - !print*,'Done, MOs saved' - - ! Deallocation, end - deallocate(new_mos) - - call wall_time(t2) - t3 = t2 - t1 - print*,'Time in apply mo rotation:', t3 - print*,'---End apply_mo_rotation---' - -end subroutine diff --git a/src/utils_trust_region/apply_mo_rotation.org b/src/utils_trust_region/apply_mo_rotation.org deleted file mode 100644 index 918581b7..00000000 --- a/src/utils_trust_region/apply_mo_rotation.org +++ /dev/null @@ -1,86 +0,0 @@ -* Apply MO rotation -Subroutine to apply the rotation matrix to the coefficients of the -MOs. - -New MOs = Old MOs . Rotation matrix - -*Compute the new MOs with the previous MOs and a rotation matrix* - -Provided: -| mo_num | integer | number of MOs | -| ao_num | integer | number of AOs | -| mo_coef(ao_num,mo_num) | double precision | coefficients of the MOs | - -Intent in: -| R(mo_num,mo_num) | double precision | rotation matrix | - -Intent out: -| prev_mos(ao_num,mo_num) | double precision | MOs before the rotation | - -Internal: -| new_mos(ao_num,mo_num) | double precision | MOs after the rotation | -| i,j | integer | indexes | -#+BEGIN_SRC f90 :comments org :tangle apply_mo_rotation.irp.f -subroutine apply_mo_rotation(R,prev_mos) - - include 'pi.h' - - BEGIN_DOC - ! Compute the new MOs knowing the rotation matrix - END_DOC - - implicit none - - ! Variables - - ! in - double precision, intent(in) :: R(mo_num,mo_num) - - ! out - double precision, intent(out) :: prev_mos(ao_num,mo_num) - - ! internal - double precision, allocatable :: new_mos(:,:) - integer :: i,j - double precision :: t1,t2,t3 - - print*,'' - print*,'---apply_mo_rotation---' - - call wall_time(t1) - - ! Allocation - allocate(new_mos(ao_num,mo_num)) - - ! Calculation - - ! Product of old MOs (mo_coef) by Rotation matrix (R) - call dgemm('N','N',ao_num,mo_num,mo_num,1d0,mo_coef,size(mo_coef,1),R,size(R,1),0d0,new_mos,size(new_mos,1)) - - prev_mos = mo_coef - mo_coef = new_mos - - !if (debug) then - ! print*,'New mo_coef : ' - ! do i = 1, mo_num - ! write(*,'(100(F10.5))') mo_coef(i,:) - ! enddo - !endif - - ! Save the new MOs and change the label - mo_label = 'MCSCF' - !call save_mos - call ezfio_set_determinants_mo_label(mo_label) - - !print*,'Done, MOs saved' - - ! Deallocation, end - deallocate(new_mos) - - call wall_time(t2) - t3 = t2 - t1 - print*,'Time in apply mo rotation:', t3 - print*,'---End apply_mo_rotation---' - -end subroutine -#+END_SRC diff --git a/src/utils_trust_region/mat_to_vec_index.irp.f b/src/utils_trust_region/mat_to_vec_index.irp.f deleted file mode 100644 index 35e12232..00000000 --- a/src/utils_trust_region/mat_to_vec_index.irp.f +++ /dev/null @@ -1,61 +0,0 @@ -! Matrix to vector index - -! *Compute the index i of a vector element from the indexes p,q of a -! matrix element* - -! Lower diagonal matrix (p,q), p > q -> vector (i) - -! If a matrix is antisymmetric it can be reshaped as a vector. And the -! vector can be reshaped as an antisymmetric matrix - -! \begin{align*} -! \begin{pmatrix} -! 0 & -1 & -2 & -4 \\ -! 1 & 0 & -3 & -5 \\ -! 2 & 3 & 0 & -6 \\ -! 4 & 5 & 6 & 0 -! \end{pmatrix} -! \Leftrightarrow -! \begin{pmatrix} -! 1 & 2 & 3 & 4 & 5 & 6 -! \end{pmatrix} -! \end{align*} - -! !!! Here the algorithm only work for the lower diagonal !!! - -! Input: -! | p,q | integer | indexes of a matrix element in the lower diagonal | -! | | | p > q, q -> column | -! | | | p -> row, | -! | | | q -> column | - -! Input: -! | i | integer | corresponding index in the vector | - - -subroutine mat_to_vec_index(p,q,i) - - include 'pi.h' - - implicit none - - ! Variables - - ! in - integer, intent(in) :: p,q - - ! out - integer, intent(out) :: i - - ! internal - integer :: a,b - double precision :: da - - ! Calculation - - a = p-1 - b = a*(a-1)/2 - - i = q+b - -end subroutine diff --git a/src/utils_trust_region/mat_to_vec_index.org b/src/utils_trust_region/mat_to_vec_index.org deleted file mode 100644 index 50840584..00000000 --- a/src/utils_trust_region/mat_to_vec_index.org +++ /dev/null @@ -1,63 +0,0 @@ -* Matrix to vector index - -*Compute the index i of a vector element from the indexes p,q of a -matrix element* - -Lower diagonal matrix (p,q), p > q -> vector (i) - -If a matrix is antisymmetric it can be reshaped as a vector. And the -vector can be reshaped as an antisymmetric matrix - -\begin{align*} -\begin{pmatrix} -0 & -1 & -2 & -4 \\ -1 & 0 & -3 & -5 \\ -2 & 3 & 0 & -6 \\ -4 & 5 & 6 & 0 -\end{pmatrix} -\Leftrightarrow -\begin{pmatrix} -1 & 2 & 3 & 4 & 5 & 6 -\end{pmatrix} -\end{align*} - -!!! Here the algorithm only work for the lower diagonal !!! - -Input: -| p,q | integer | indexes of a matrix element in the lower diagonal | -| | | p > q, q -> column | -| | | p -> row, | -| | | q -> column | - -Input: -| i | integer | corresponding index in the vector | - -#+BEGIN_SRC f90 :comments org :tangle mat_to_vec_index.irp.f -subroutine mat_to_vec_index(p,q,i) - - include 'pi.h' - - implicit none - - ! Variables - - ! in - integer, intent(in) :: p,q - - ! out - integer, intent(out) :: i - - ! internal - integer :: a,b - double precision :: da - - ! Calculation - - a = p-1 - b = a*(a-1)/2 - - i = q+b - -end subroutine -#+END_SRC - diff --git a/src/utils_trust_region/pi.h b/src/utils_trust_region/pi.h deleted file mode 100644 index bbfabfec..00000000 --- a/src/utils_trust_region/pi.h +++ /dev/null @@ -1,2 +0,0 @@ - !logical, parameter :: debug=.False. - double precision, parameter :: pi = 3.1415926535897932d0 diff --git a/src/utils_trust_region/rotation_matrix.irp.f b/src/utils_trust_region/rotation_matrix.irp.f deleted file mode 100644 index 4738fd67..00000000 --- a/src/utils_trust_region/rotation_matrix.irp.f +++ /dev/null @@ -1,443 +0,0 @@ -! Rotation matrix - -! *Build a rotation matrix from an antisymmetric matrix* - -! Compute a rotation matrix $\textbf{R}$ from an antisymmetric matrix $$\textbf{A}$$ such as : -! $$ -! \textbf{R}=\exp(\textbf{A}) -! $$ - -! So : -! \begin{align*} -! \textbf{R}=& \exp(\textbf{A}) \\ -! =& \sum_k^{\infty} \frac{1}{k!}\textbf{A}^k \\ -! =& \textbf{W} \cdot \cos(\tau) \cdot \textbf{W}^{\dagger} + \textbf{W} \cdot \tau^{-1} \cdot \sin(\tau) \cdot \textbf{W}^{\dagger} \cdot \textbf{A} -! \end{align*} - -! With : -! $\textbf{W}$ : eigenvectors of $\textbf{A}^2$ -! $\tau$ : $\sqrt{-x}$ -! $x$ : eigenvalues of $\textbf{A}^2$ - -! Input: -! | A(n,n) | double precision | antisymmetric matrix | -! | n | integer | number of columns of the A matrix | -! | LDA | integer | specifies the leading dimension of A, must be at least max(1,n) | -! | LDR | integer | specifies the leading dimension of R, must be at least max(1,n) | - -! Output: -! | R(n,n) | double precision | Rotation matrix | -! | info | integer | if info = 0, the execution is successful | -! | | | if info = k, the k-th parameter has an illegal value | -! | | | if info = -k, the algorithm failed | - -! Internal: -! | B(n,n) | double precision | B = A.A | -! | work(lwork,n) | double precision | work matrix for dysev, dimension max(1,lwork) | -! | lwork | integer | dimension of the syev work array >= max(1, 3n-1) | -! | W(n,n) | double precision | eigenvectors of B | -! | e_val(n) | double precision | eigenvalues of B | -! | m_diag(n,n) | double precision | diagonal matrix with the eigenvalues of B | -! | cos_tau(n,n) | double precision | diagonal matrix with cos(tau) values | -! | sin_tau(n,n) | double precision | diagonal matrix with sin cos(tau) values | -! | tau_m1(n,n) | double precision | diagonal matrix with (tau)^-1 values | -! | part_1(n,n) | double precision | matrix W.cos_tau.W^t | -! | part_1a(n,n) | double precision | matrix cos_tau.W^t | -! | part_2(n,n) | double precision | matrix W.tau_m1.sin_tau.W^t.A | -! | part_2a(n,n) | double precision | matrix W^t.A | -! | part_2b(n,n) | double precision | matrix sin_tau.W^t.A | -! | part_2c(n,n) | double precision | matrix tau_m1.sin_tau.W^t.A | -! | RR_t(n,n) | double precision | R.R^t must be equal to the identity<=> R.R^t-1=0 <=> norm = 0 | -! | norm | integer | norm of R.R^t-1, must be equal to 0 | -! | i,j | integer | indexes | - -! Functions: -! | dnrm2 | double precision | Lapack function, compute the norm of a matrix | -! | disnan | logical | Lapack function, check if an element is NaN | - - - -subroutine rotation_matrix(A,LDA,R,LDR,n,info,enforce_step_cancellation) - - implicit none - - BEGIN_DOC - ! Rotation matrix to rotate the molecular orbitals. - ! If the rotation is too large the transformation is not unitary and must be cancelled. - END_DOC - - include 'pi.h' - - ! Variables - - ! in - integer, intent(in) :: n,LDA,LDR - double precision, intent(inout) :: A(LDA,n) - - ! out - double precision, intent(out) :: R(LDR,n) - integer, intent(out) :: info - logical, intent(out) :: enforce_step_cancellation - - ! internal - double precision, allocatable :: B(:,:) - double precision, allocatable :: work(:,:) - double precision, allocatable :: W(:,:), e_val(:) - double precision, allocatable :: m_diag(:,:),cos_tau(:,:),sin_tau(:,:),tau_m1(:,:) - double precision, allocatable :: part_1(:,:),part_1a(:,:) - double precision, allocatable :: part_2(:,:),part_2a(:,:),part_2b(:,:),part_2c(:,:) - double precision, allocatable :: RR_t(:,:) - integer :: i,j - integer :: info2, lwork ! for dsyev - double precision :: norm, max_elem, max_elem_A, t1,t2,t3 - - ! function - double precision :: dnrm2 - logical :: disnan - - print*,'' - print*,'---rotation_matrix---' - - call wall_time(t1) - - ! Allocation - allocate(B(n,n)) - allocate(m_diag(n,n),cos_tau(n,n),sin_tau(n,n),tau_m1(n,n)) - allocate(W(n,n),e_val(n)) - allocate(part_1(n,n),part_1a(n,n)) - allocate(part_2(n,n),part_2a(n,n),part_2b(n,n),part_2c(n,n)) - allocate(RR_t(n,n)) - -! Pre-conditions - -! Initialization -info=0 -enforce_step_cancellation = .False. - -! Size of matrix A must be at least 1 by 1 -if (n<1) then - info = 3 - print*, 'WARNING: invalid parameter 5' - print*, 'n<1' - return -endif - -! Leading dimension of A must be >= n -if (LDA < n) then - info = 25 - print*, 'WARNING: invalid parameter 2 or 5' - print*, 'LDA < n' - return -endif - -! Leading dimension of A must be >= n -if (LDR < n) then - info = 4 - print*, 'WARNING: invalid parameter 4' - print*, 'LDR < n' - return -endif - -! Matrix elements of A must by non-NaN -do j = 1, n - do i = 1, n - if (disnan(A(i,j))) then - info=1 - print*, 'WARNING: invalid parameter 1' - print*, 'NaN element in A matrix' - return - endif - enddo -enddo - -do i = 1, n - if (A(i,i) /= 0d0) then - print*, 'WARNING: matrix A is not antisymmetric' - print*, 'Non 0 element on the diagonal', i, A(i,i) - call ABORT - endif -enddo - -do j = 1, n - do i = 1, n - if (A(i,j)+A(j,i)>1d-16) then - print*, 'WANRING: matrix A is not antisymmetric' - print*, 'A(i,j) /= - A(j,i):', i,j,A(i,j), A(j,i) - print*, 'diff:', A(i,j)+A(j,i) - call ABORT - endif - enddo -enddo - -! Fix for too big elements ! bad idea better to cancel if the error is too big -!do j = 1, n -! do i = 1, n -! A(i,j) = mod(A(i,j),2d0*pi) -! if (dabs(A(i,j)) > pi) then -! A(i,j) = 0d0 -! endif -! enddo -!enddo - -max_elem_A = 0d0 -do j = 1, n - do i = 1, n - if (ABS(A(i,j)) > ABS(max_elem_A)) then - max_elem_A = A(i,j) - endif - enddo -enddo -print*,'max element in A', max_elem_A - -if (ABS(max_elem_A) > 2 * pi) then - print*,'' - print*,'WARNING: ABS(max_elem_A) > 2 pi ' - print*,'' -endif - -! B=A.A -! - Calculation of the matrix $\textbf{B} = \textbf{A}^2$ -! - Diagonalization of $\textbf{B}$ -! W, the eigenvectors -! e_val, the eigenvalues - - -! Compute B=A.A - -call dgemm('N','N',n,n,n,1d0,A,size(A,1),A,size(A,1),0d0,B,size(B,1)) - -! Copy B in W, diagonalization will put the eigenvectors in W -W=B - -! Diagonalization of B -! Eigenvalues -> e_val -! Eigenvectors -> W -lwork = 3*n-1 -allocate(work(lwork,n)) - -print*,'Starting diagonalization ...' - -call dsyev('V','U',n,W,size(W,1),e_val,work,lwork,info2) - -deallocate(work) - -if (info2 == 0) then - print*, 'Diagonalization : Done' -elseif (info2 < 0) then - print*, 'WARNING: error in the diagonalization' - print*, 'Illegal value of the ', info2,'-th parameter' -else - print*, "WARNING: Diagonalization failed to converge" -endif - -! Tau^-1, cos(tau), sin(tau) -! $$\tau = \sqrt{-x}$$ -! - Calculation of $\cos(\tau)$ $\Leftrightarrow$ $\cos(\sqrt{-x})$ -! - Calculation of $\sin(\tau)$ $\Leftrightarrow$ $\sin(\sqrt{-x})$ -! - Calculation of $\tau^{-1}$ $\Leftrightarrow$ $(\sqrt{-x})^{-1}$ -! These matrices are diagonals - -! Diagonal matrix m_diag -do j = 1, n - if (e_val(j) >= -1d-12) then !0.d0) then !!! e_avl(i) must be < -1d-12 to avoid numerical problems - e_val(j) = 0.d0 - else - e_val(j) = - e_val(j) - endif -enddo - -m_diag = 0.d0 -do i = 1, n - m_diag(i,i) = e_val(i) -enddo - -! cos_tau -do j = 1, n - do i = 1, n - if (i==j) then - cos_tau(i,j) = dcos(dsqrt(e_val(i))) - else - cos_tau(i,j) = 0d0 - endif - enddo -enddo - -! sin_tau -do j = 1, n - do i = 1, n - if (i==j) then - sin_tau(i,j) = dsin(dsqrt(e_val(i))) - else - sin_tau(i,j) = 0d0 - endif - enddo -enddo - -! Debug, display the cos_tau and sin_tau matrix -!if (debug) then -! print*, 'cos_tau' -! do i = 1, n -! print*, cos_tau(i,:) -! enddo -! print*, 'sin_tau' -! do i = 1, n -! print*, sin_tau(i,:) -! enddo -!endif - -! tau^-1 -do j = 1, n - do i = 1, n - if ((i==j) .and. (e_val(i) > 1d-16)) then!0d0)) then !!! Convergence problem can come from here if the threshold is too big/small - tau_m1(i,j) = 1d0/(dsqrt(e_val(i))) - else - tau_m1(i,j) = 0d0 - endif - enddo -enddo - -max_elem = 0d0 -do i = 1, n - if (ABS(tau_m1(i,i)) > ABS(max_elem)) then - max_elem = tau_m1(i,i) - endif -enddo -print*,'max elem tau^-1:', max_elem - -! Debug -!print*,'eigenvalues:' -!do i = 1, n -! print*, e_val(i) -!enddo - -!Debug, display tau^-1 -!if (debug) then -! print*, 'tau^-1' -! do i = 1, n -! print*,tau_m1(i,:) -! enddo -!endif - -! Rotation matrix -! \begin{align*} -! \textbf{R} = \textbf{W} \cos(\tau) \textbf{W}^{\dagger} + \textbf{W} \tau^{-1} \sin(\tau) \textbf{W}^{\dagger} \textbf{A} -! \end{align*} -! \begin{align*} -! \textbf{Part1} = \textbf{W} \cos(\tau) \textbf{W}^{\dagger} -! \end{align*} -! \begin{align*} -! \textbf{Part2} = \textbf{W} \tau^{-1} \sin(\tau) \textbf{W}^{\dagger} \textbf{A} -! \end{align*} - -! First: -! part_1 = dgemm(W, dgemm(cos_tau, W^t)) -! part_1a = dgemm(cos_tau, W^t) -! part_1 = dgemm(W, part_1a) -! And: -! part_2= dgemm(W, dgemm(tau_m1, dgemm(sin_tau, dgemm(W^t, A)))) -! part_2a = dgemm(W^t, A) -! part_2b = dgemm(sin_tau, part_2a) -! part_2c = dgemm(tau_m1, part_2b) -! part_2 = dgemm(W, part_2c) -! Finally: -! Rotation matrix, R = part_1+part_2 - -! If $R$ is a rotation matrix: -! $R.R^T=R^T.R=\textbf{1}$ - -! part_1 -call dgemm('N','T',n,n,n,1d0,cos_tau,size(cos_tau,1),W,size(W,1),0d0,part_1a,size(part_1a,1)) -call dgemm('N','N',n,n,n,1d0,W,size(W,1),part_1a,size(part_1a,1),0d0,part_1,size(part_1,1)) - -! part_2 -call dgemm('T','N',n,n,n,1d0,W,size(W,1),A,size(A,1),0d0,part_2a,size(part_2a,1)) -call dgemm('N','N',n,n,n,1d0,sin_tau,size(sin_tau,1),part_2a,size(part_2a,1),0d0,part_2b,size(part_2b,1)) -call dgemm('N','N',n,n,n,1d0,tau_m1,size(tau_m1,1),part_2b,size(part_2b,1),0d0,part_2c,size(part_2c,1)) -call dgemm('N','N',n,n,n,1d0,W,size(W,1),part_2c,size(part_2c,1),0d0,part_2,size(part_2,1)) - -! Rotation matrix R -R = part_1 + part_2 - -! Matrix check -! R.R^t and R^t.R must be equal to identity matrix -do j = 1, n - do i=1,n - if (i==j) then - RR_t(i,j) = 1d0 - else - RR_t(i,j) = 0d0 - endif - enddo -enddo - -call dgemm('N','T',n,n,n,1d0,R,size(R,1),R,size(R,1),-1d0,RR_t,size(RR_t,1)) - -norm = dnrm2(n*n,RR_t,1) -print*, 'Rotation matrix check, norm R.R^T = ', norm - -! Debug -!if (debug) then -! print*, 'RR_t' -! do i = 1, n -! print*, RR_t(i,:) -! enddo -!endif - -! Post conditions - -! Check if R.R^T=1 -max_elem = 0d0 -do j = 1, n - do i = 1, n - if (ABS(RR_t(i,j)) > ABS(max_elem)) then - max_elem = RR_t(i,j) - endif - enddo -enddo - -print*, 'Max error in R.R^T:', max_elem -print*, 'e_val(1):', e_val(1) -print*, 'e_val(n):', e_val(n) -print*, 'max elem in A:', max_elem_A - -if (ABS(max_elem) > 1d-12) then - print*, 'WARNING: max error in R.R^T > 1d-12' - print*, 'Enforce the step cancellation' - enforce_step_cancellation = .True. -endif - -! Matrix elements of R must by non-NaN -do j = 1,n - do i = 1,LDR - if (disnan(R(i,j))) then - info = 666 - print*, 'NaN in rotation matrix' - call ABORT - endif - enddo -enddo - -! Display -!if (debug) then -! print*,'Rotation matrix :' -! do i = 1, n -! write(*,'(100(F10.5))') R(i,:) -! enddo -!endif - -! Deallocation, end - -deallocate(B) - deallocate(m_diag,cos_tau,sin_tau,tau_m1) - deallocate(W,e_val) - deallocate(part_1,part_1a) - deallocate(part_2,part_2a,part_2b,part_2c) - deallocate(RR_t) - - call wall_time(t2) - t3 = t2-t1 - print*,'Time in rotation matrix:', t3 - - print*,'---End rotation_matrix---' - -end subroutine diff --git a/src/utils_trust_region/rotation_matrix.org b/src/utils_trust_region/rotation_matrix.org deleted file mode 100644 index 73ba0298..00000000 --- a/src/utils_trust_region/rotation_matrix.org +++ /dev/null @@ -1,454 +0,0 @@ -* Rotation matrix - -*Build a rotation matrix from an antisymmetric matrix* - -Compute a rotation matrix $\textbf{R}$ from an antisymmetric matrix $$\textbf{A}$$ such as : -$$ -\textbf{R}=\exp(\textbf{A}) -$$ - -So : -\begin{align*} -\textbf{R}=& \exp(\textbf{A}) \\ -=& \sum_k^{\infty} \frac{1}{k!}\textbf{A}^k \\ -=& \textbf{W} \cdot \cos(\tau) \cdot \textbf{W}^{\dagger} + \textbf{W} \cdot \tau^{-1} \cdot \sin(\tau) \cdot \textbf{W}^{\dagger} \cdot \textbf{A} -\end{align*} - -With : -$\textbf{W}$ : eigenvectors of $\textbf{A}^2$ -$\tau$ : $\sqrt{-x}$ -$x$ : eigenvalues of $\textbf{A}^2$ - -Input: -| A(n,n) | double precision | antisymmetric matrix | -| n | integer | number of columns of the A matrix | -| LDA | integer | specifies the leading dimension of A, must be at least max(1,n) | -| LDR | integer | specifies the leading dimension of R, must be at least max(1,n) | - -Output: -| R(n,n) | double precision | Rotation matrix | -| info | integer | if info = 0, the execution is successful | -| | | if info = k, the k-th parameter has an illegal value | -| | | if info = -k, the algorithm failed | - -Internal: -| B(n,n) | double precision | B = A.A | -| work(lwork,n) | double precision | work matrix for dysev, dimension max(1,lwork) | -| lwork | integer | dimension of the syev work array >= max(1, 3n-1) | -| W(n,n) | double precision | eigenvectors of B | -| e_val(n) | double precision | eigenvalues of B | -| m_diag(n,n) | double precision | diagonal matrix with the eigenvalues of B | -| cos_tau(n,n) | double precision | diagonal matrix with cos(tau) values | -| sin_tau(n,n) | double precision | diagonal matrix with sin cos(tau) values | -| tau_m1(n,n) | double precision | diagonal matrix with (tau)^-1 values | -| part_1(n,n) | double precision | matrix W.cos_tau.W^t | -| part_1a(n,n) | double precision | matrix cos_tau.W^t | -| part_2(n,n) | double precision | matrix W.tau_m1.sin_tau.W^t.A | -| part_2a(n,n) | double precision | matrix W^t.A | -| part_2b(n,n) | double precision | matrix sin_tau.W^t.A | -| part_2c(n,n) | double precision | matrix tau_m1.sin_tau.W^t.A | -| RR_t(n,n) | double precision | R.R^t must be equal to the identity<=> R.R^t-1=0 <=> norm = 0 | -| norm | integer | norm of R.R^t-1, must be equal to 0 | -| i,j | integer | indexes | - -Functions: -| dnrm2 | double precision | Lapack function, compute the norm of a matrix | -| disnan | logical | Lapack function, check if an element is NaN | - - -#+BEGIN_SRC f90 :comments org :tangle rotation_matrix.irp.f -subroutine rotation_matrix(A,LDA,R,LDR,n,info,enforce_step_cancellation) - - implicit none - - BEGIN_DOC - ! Rotation matrix to rotate the molecular orbitals. - ! If the rotation is too large the transformation is not unitary and must be cancelled. - END_DOC - - include 'pi.h' - - ! Variables - - ! in - integer, intent(in) :: n,LDA,LDR - double precision, intent(inout) :: A(LDA,n) - - ! out - double precision, intent(out) :: R(LDR,n) - integer, intent(out) :: info - logical, intent(out) :: enforce_step_cancellation - - ! internal - double precision, allocatable :: B(:,:) - double precision, allocatable :: work(:,:) - double precision, allocatable :: W(:,:), e_val(:) - double precision, allocatable :: m_diag(:,:),cos_tau(:,:),sin_tau(:,:),tau_m1(:,:) - double precision, allocatable :: part_1(:,:),part_1a(:,:) - double precision, allocatable :: part_2(:,:),part_2a(:,:),part_2b(:,:),part_2c(:,:) - double precision, allocatable :: RR_t(:,:) - integer :: i,j - integer :: info2, lwork ! for dsyev - double precision :: norm, max_elem, max_elem_A, t1,t2,t3 - - ! function - double precision :: dnrm2 - logical :: disnan - - print*,'' - print*,'---rotation_matrix---' - - call wall_time(t1) - - ! Allocation - allocate(B(n,n)) - allocate(m_diag(n,n),cos_tau(n,n),sin_tau(n,n),tau_m1(n,n)) - allocate(W(n,n),e_val(n)) - allocate(part_1(n,n),part_1a(n,n)) - allocate(part_2(n,n),part_2a(n,n),part_2b(n,n),part_2c(n,n)) - allocate(RR_t(n,n)) -#+END_SRC - -** Pre-conditions -#+BEGIN_SRC f90 :comments org :tangle rotation_matrix.irp.f - ! Initialization - info=0 - enforce_step_cancellation = .False. - - ! Size of matrix A must be at least 1 by 1 - if (n<1) then - info = 3 - print*, 'WARNING: invalid parameter 5' - print*, 'n<1' - return - endif - - ! Leading dimension of A must be >= n - if (LDA < n) then - info = 25 - print*, 'WARNING: invalid parameter 2 or 5' - print*, 'LDA < n' - return - endif - - ! Leading dimension of A must be >= n - if (LDR < n) then - info = 4 - print*, 'WARNING: invalid parameter 4' - print*, 'LDR < n' - return - endif - - ! Matrix elements of A must by non-NaN - do j = 1, n - do i = 1, n - if (disnan(A(i,j))) then - info=1 - print*, 'WARNING: invalid parameter 1' - print*, 'NaN element in A matrix' - return - endif - enddo - enddo - - do i = 1, n - if (A(i,i) /= 0d0) then - print*, 'WARNING: matrix A is not antisymmetric' - print*, 'Non 0 element on the diagonal', i, A(i,i) - call ABORT - endif - enddo - - do j = 1, n - do i = 1, n - if (A(i,j)+A(j,i)>1d-16) then - print*, 'WANRING: matrix A is not antisymmetric' - print*, 'A(i,j) /= - A(j,i):', i,j,A(i,j), A(j,i) - print*, 'diff:', A(i,j)+A(j,i) - call ABORT - endif - enddo - enddo - - ! Fix for too big elements ! bad idea better to cancel if the error is too big - !do j = 1, n - ! do i = 1, n - ! A(i,j) = mod(A(i,j),2d0*pi) - ! if (dabs(A(i,j)) > pi) then - ! A(i,j) = 0d0 - ! endif - ! enddo - !enddo - - max_elem_A = 0d0 - do j = 1, n - do i = 1, n - if (ABS(A(i,j)) > ABS(max_elem_A)) then - max_elem_A = A(i,j) - endif - enddo - enddo - print*,'max element in A', max_elem_A - - if (ABS(max_elem_A) > 2 * pi) then - print*,'' - print*,'WARNING: ABS(max_elem_A) > 2 pi ' - print*,'' - endif - -#+END_SRC - -** Calculations - -*** B=A.A - - Calculation of the matrix $\textbf{B} = \textbf{A}^2$ - - Diagonalization of $\textbf{B}$ - W, the eigenvectors - e_val, the eigenvalues - - #+BEGIN_SRC f90 :comments org :tangle rotation_matrix.irp.f - ! Compute B=A.A - - call dgemm('N','N',n,n,n,1d0,A,size(A,1),A,size(A,1),0d0,B,size(B,1)) - - ! Copy B in W, diagonalization will put the eigenvectors in W - W=B - - ! Diagonalization of B - ! Eigenvalues -> e_val - ! Eigenvectors -> W - lwork = 3*n-1 - allocate(work(lwork,n)) - - print*,'Starting diagonalization ...' - - call dsyev('V','U',n,W,size(W,1),e_val,work,lwork,info2) - - deallocate(work) - - if (info2 == 0) then - print*, 'Diagonalization : Done' - elseif (info2 < 0) then - print*, 'WARNING: error in the diagonalization' - print*, 'Illegal value of the ', info2,'-th parameter' - else - print*, "WARNING: Diagonalization failed to converge" - endif - #+END_SRC - -*** Tau^-1, cos(tau), sin(tau) - $$\tau = \sqrt{-x}$$ - - Calculation of $\cos(\tau)$ $\Leftrightarrow$ $\cos(\sqrt{-x})$ - - Calculation of $\sin(\tau)$ $\Leftrightarrow$ $\sin(\sqrt{-x})$ - - Calculation of $\tau^{-1}$ $\Leftrightarrow$ $(\sqrt{-x})^{-1}$ - These matrices are diagonals - #+BEGIN_SRC f90 :comments org :tangle rotation_matrix.irp.f - ! Diagonal matrix m_diag - do j = 1, n - if (e_val(j) >= -1d-12) then !0.d0) then !!! e_avl(i) must be < -1d-12 to avoid numerical problems - e_val(j) = 0.d0 - else - e_val(j) = - e_val(j) - endif - enddo - - m_diag = 0.d0 - do i = 1, n - m_diag(i,i) = e_val(i) - enddo - - ! cos_tau - do j = 1, n - do i = 1, n - if (i==j) then - cos_tau(i,j) = dcos(dsqrt(e_val(i))) - else - cos_tau(i,j) = 0d0 - endif - enddo - enddo - - ! sin_tau - do j = 1, n - do i = 1, n - if (i==j) then - sin_tau(i,j) = dsin(dsqrt(e_val(i))) - else - sin_tau(i,j) = 0d0 - endif - enddo - enddo - - ! Debug, display the cos_tau and sin_tau matrix - !if (debug) then - ! print*, 'cos_tau' - ! do i = 1, n - ! print*, cos_tau(i,:) - ! enddo - ! print*, 'sin_tau' - ! do i = 1, n - ! print*, sin_tau(i,:) - ! enddo - !endif - - ! tau^-1 - do j = 1, n - do i = 1, n - if ((i==j) .and. (e_val(i) > 1d-16)) then!0d0)) then !!! Convergence problem can come from here if the threshold is too big/small - tau_m1(i,j) = 1d0/(dsqrt(e_val(i))) - else - tau_m1(i,j) = 0d0 - endif - enddo - enddo - - max_elem = 0d0 - do i = 1, n - if (ABS(tau_m1(i,i)) > ABS(max_elem)) then - max_elem = tau_m1(i,i) - endif - enddo - print*,'max elem tau^-1:', max_elem - - ! Debug - !print*,'eigenvalues:' - !do i = 1, n - ! print*, e_val(i) - !enddo - - !Debug, display tau^-1 - !if (debug) then - ! print*, 'tau^-1' - ! do i = 1, n - ! print*,tau_m1(i,:) - ! enddo - !endif - #+END_SRC - -*** Rotation matrix - \begin{align*} - \textbf{R} = \textbf{W} \cos(\tau) \textbf{W}^{\dagger} + \textbf{W} \tau^{-1} \sin(\tau) \textbf{W}^{\dagger} \textbf{A} - \end{align*} - \begin{align*} - \textbf{Part1} = \textbf{W} \cos(\tau) \textbf{W}^{\dagger} - \end{align*} - \begin{align*} - \textbf{Part2} = \textbf{W} \tau^{-1} \sin(\tau) \textbf{W}^{\dagger} \textbf{A} - \end{align*} - - First: - part_1 = dgemm(W, dgemm(cos_tau, W^t)) - part_1a = dgemm(cos_tau, W^t) - part_1 = dgemm(W, part_1a) - And: - part_2= dgemm(W, dgemm(tau_m1, dgemm(sin_tau, dgemm(W^t, A)))) - part_2a = dgemm(W^t, A) - part_2b = dgemm(sin_tau, part_2a) - part_2c = dgemm(tau_m1, part_2b) - part_2 = dgemm(W, part_2c) - Finally: - Rotation matrix, R = part_1+part_2 - - If $R$ is a rotation matrix: - $R.R^T=R^T.R=\textbf{1}$ - #+BEGIN_SRC f90 :comments org :tangle rotation_matrix.irp.f - ! part_1 - call dgemm('N','T',n,n,n,1d0,cos_tau,size(cos_tau,1),W,size(W,1),0d0,part_1a,size(part_1a,1)) - call dgemm('N','N',n,n,n,1d0,W,size(W,1),part_1a,size(part_1a,1),0d0,part_1,size(part_1,1)) - - ! part_2 - call dgemm('T','N',n,n,n,1d0,W,size(W,1),A,size(A,1),0d0,part_2a,size(part_2a,1)) - call dgemm('N','N',n,n,n,1d0,sin_tau,size(sin_tau,1),part_2a,size(part_2a,1),0d0,part_2b,size(part_2b,1)) - call dgemm('N','N',n,n,n,1d0,tau_m1,size(tau_m1,1),part_2b,size(part_2b,1),0d0,part_2c,size(part_2c,1)) - call dgemm('N','N',n,n,n,1d0,W,size(W,1),part_2c,size(part_2c,1),0d0,part_2,size(part_2,1)) - - ! Rotation matrix R - R = part_1 + part_2 - - ! Matrix check - ! R.R^t and R^t.R must be equal to identity matrix - do j = 1, n - do i=1,n - if (i==j) then - RR_t(i,j) = 1d0 - else - RR_t(i,j) = 0d0 - endif - enddo - enddo - - call dgemm('N','T',n,n,n,1d0,R,size(R,1),R,size(R,1),-1d0,RR_t,size(RR_t,1)) - - norm = dnrm2(n*n,RR_t,1) - print*, 'Rotation matrix check, norm R.R^T = ', norm - - ! Debug - !if (debug) then - ! print*, 'RR_t' - ! do i = 1, n - ! print*, RR_t(i,:) - ! enddo - !endif - #+END_SRC - -*** Post conditions - #+BEGIN_SRC f90 :comments org :tangle rotation_matrix.irp.f - ! Check if R.R^T=1 - max_elem = 0d0 - do j = 1, n - do i = 1, n - if (ABS(RR_t(i,j)) > ABS(max_elem)) then - max_elem = RR_t(i,j) - endif - enddo - enddo - - print*, 'Max error in R.R^T:', max_elem - print*, 'e_val(1):', e_val(1) - print*, 'e_val(n):', e_val(n) - print*, 'max elem in A:', max_elem_A - - if (ABS(max_elem) > 1d-12) then - print*, 'WARNING: max error in R.R^T > 1d-12' - print*, 'Enforce the step cancellation' - enforce_step_cancellation = .True. - endif - - ! Matrix elements of R must by non-NaN - do j = 1,n - do i = 1,LDR - if (disnan(R(i,j))) then - info = 666 - print*, 'NaN in rotation matrix' - call ABORT - endif - enddo - enddo - - ! Display - !if (debug) then - ! print*,'Rotation matrix :' - ! do i = 1, n - ! write(*,'(100(F10.5))') R(i,:) - ! enddo - !endif - #+END_SRC - -** Deallocation, end - #+BEGIN_SRC f90 :comments org :tangle rotation_matrix.irp.f - deallocate(B) - deallocate(m_diag,cos_tau,sin_tau,tau_m1) - deallocate(W,e_val) - deallocate(part_1,part_1a) - deallocate(part_2,part_2a,part_2b,part_2c) - deallocate(RR_t) - - call wall_time(t2) - t3 = t2-t1 - print*,'Time in rotation matrix:', t3 - - print*,'---End rotation_matrix---' - -end subroutine - #+END_SRC - diff --git a/src/utils_trust_region/sub_to_full_rotation_matrix.irp.f b/src/utils_trust_region/sub_to_full_rotation_matrix.irp.f deleted file mode 100644 index bdd1f6ba..00000000 --- a/src/utils_trust_region/sub_to_full_rotation_matrix.irp.f +++ /dev/null @@ -1,64 +0,0 @@ -! Rotation matrix in a subspace to rotation matrix in the full space - -! Usually, we are using a list of MOs, for exemple the active ones. When -! we compute a rotation matrix to rotate the MOs, we just compute a -! rotation matrix for these MOs in order to reduce the size of the -! matrix which has to be computed. Since the computation of a rotation -! matrix scale in $O(N^3)$ with $N$ the number of MOs, it's better to -! reuce the number of MOs involved. -! After that we replace the rotation matrix in the full space by -! building the elements of the rotation matrix in the full space from -! the elements of the rotation matrix in the subspace and adding some 0 -! on the extradiagonal elements and some 1 on the diagonal elements, -! for the MOs that are not involved in the rotation. - -! Provided: -! | mo_num | integer | Number of MOs | - -! Input: -! | m | integer | Size of tmp_list, m <= mo_num | -! | tmp_list(m) | integer | List of MOs | -! | tmp_R(m,m) | double precision | Rotation matrix in the space of | -! | | | the MOs containing by tmp_list | - -! Output: -! | R(mo_num,mo_num | double precision | Rotation matrix in the space | -! | | | of all the MOs | - -! Internal: -! | i,j | integer | indexes in the full space | -! | tmp_i,tmp_j | integer | indexes in the subspace | - - -subroutine sub_to_full_rotation_matrix(m,tmp_list,tmp_R,R) - - BEGIN_DOC - ! Compute the full rotation matrix from a smaller one - END_DOC - - implicit none - - ! in - integer, intent(in) :: m, tmp_list(m) - double precision, intent(in) :: tmp_R(m,m) - - ! out - double precision, intent(out) :: R(mo_num,mo_num) - - ! internal - integer :: i,j,tmp_i,tmp_j - - ! tmp_R to R, subspace to full space - R = 0d0 - do i = 1, mo_num - R(i,i) = 1d0 ! 1 on the diagonal because it is a rotation matrix, 1 = nothing change for the corresponding orbital - enddo - do tmp_j = 1, m - j = tmp_list(tmp_j) - do tmp_i = 1, m - i = tmp_list(tmp_i) - R(i,j) = tmp_R(tmp_i,tmp_j) - enddo - enddo - -end diff --git a/src/utils_trust_region/sub_to_full_rotation_matrix.org b/src/utils_trust_region/sub_to_full_rotation_matrix.org deleted file mode 100644 index 16434dc8..00000000 --- a/src/utils_trust_region/sub_to_full_rotation_matrix.org +++ /dev/null @@ -1,65 +0,0 @@ -* Rotation matrix in a subspace to rotation matrix in the full space - -Usually, we are using a list of MOs, for exemple the active ones. When -we compute a rotation matrix to rotate the MOs, we just compute a -rotation matrix for these MOs in order to reduce the size of the -matrix which has to be computed. Since the computation of a rotation -matrix scale in $O(N^3)$ with $N$ the number of MOs, it's better to -reuce the number of MOs involved. -After that we replace the rotation matrix in the full space by -building the elements of the rotation matrix in the full space from -the elements of the rotation matrix in the subspace and adding some 0 -on the extradiagonal elements and some 1 on the diagonal elements, -for the MOs that are not involved in the rotation. - -Provided: -| mo_num | integer | Number of MOs | - -Input: -| m | integer | Size of tmp_list, m <= mo_num | -| tmp_list(m) | integer | List of MOs | -| tmp_R(m,m) | double precision | Rotation matrix in the space of | -| | | the MOs containing by tmp_list | - -Output: -| R(mo_num,mo_num | double precision | Rotation matrix in the space | -| | | of all the MOs | - -Internal: -| i,j | integer | indexes in the full space | -| tmp_i,tmp_j | integer | indexes in the subspace | - -#+BEGIN_SRC f90 :comments org :tangle sub_to_full_rotation_matrix.irp.f -subroutine sub_to_full_rotation_matrix(m,tmp_list,tmp_R,R) - - BEGIN_DOC - ! Compute the full rotation matrix from a smaller one - END_DOC - - implicit none - - ! in - integer, intent(in) :: m, tmp_list(m) - double precision, intent(in) :: tmp_R(m,m) - - ! out - double precision, intent(out) :: R(mo_num,mo_num) - - ! internal - integer :: i,j,tmp_i,tmp_j - - ! tmp_R to R, subspace to full space - R = 0d0 - do i = 1, mo_num - R(i,i) = 1d0 ! 1 on the diagonal because it is a rotation matrix, 1 = nothing change for the corresponding orbital - enddo - do tmp_j = 1, m - j = tmp_list(tmp_j) - do tmp_i = 1, m - i = tmp_list(tmp_i) - R(i,j) = tmp_R(tmp_i,tmp_j) - enddo - enddo - -end -#+END_SRC diff --git a/src/utils_trust_region/trust_region_expected_e.irp.f b/src/utils_trust_region/trust_region_expected_e.irp.f deleted file mode 100644 index b7d849d1..00000000 --- a/src/utils_trust_region/trust_region_expected_e.irp.f +++ /dev/null @@ -1,119 +0,0 @@ -! Predicted energy : e_model - -! *Compute the energy predicted by the Taylor series* - -! The energy is predicted using a Taylor expansion truncated at te 2nd -! order : - -! \begin{align*} -! E_{k+1} = E_{k} + \textbf{g}_k^{T} \cdot \textbf{x}_{k+1} + \frac{1}{2} \cdot \textbf{x}_{k+1}^T \cdot \textbf{H}_{k} \cdot \textbf{x}_{k+1} + \mathcal{O}(\textbf{x}_{k+1}^2) -! \end{align*} - -! Input: -! | n | integer | m*(m-1)/2 | -! | v_grad(n) | double precision | gradient | -! | H(n,n) | double precision | hessian | -! | x(n) | double precision | Step in the trust region | -! | prev_energy | double precision | previous energy | - -! Output: -! | e_model | double precision | predicted energy after the rotation of the MOs | - -! Internal: -! | part_1 | double precision | v_grad^T.x | -! | part_2 | double precision | 1/2 . x^T.H.x | -! | part_2a | double precision | H.x | -! | i,j | integer | indexes | - -! Function: -! | ddot | double precision | dot product (Lapack) | - - -subroutine trust_region_expected_e(n,v_grad,H,x,prev_energy,e_model) - - include 'pi.h' - - BEGIN_DOC - ! Compute the expected criterion/energy after the application of the step x - END_DOC - - implicit none - - ! Variables - - ! in - integer, intent(in) :: n - double precision, intent(in) :: v_grad(n),H(n,n),x(n) - double precision, intent(in) :: prev_energy - - ! out - double precision, intent(out) :: e_model - - ! internal - double precision :: part_1, part_2, t1,t2,t3 - double precision, allocatable :: part_2a(:) - - integer :: i,j - - !Function - double precision :: ddot - - print*,'' - print*,'---Trust_e_model---' - - call wall_time(t1) - - ! Allocation - allocate(part_2a(n)) - -! Calculations - -! part_1 corresponds to the product g.x -! part_2a corresponds to the product H.x -! part_2 corresponds to the product 0.5*(x^T.H.x) - -! TODO: remove the dot products - - -! Product v_grad.x - part_1 = ddot(n,v_grad,1,x,1) - - !if (debug) then - print*,'g.x : ', part_1 - !endif - - ! Product H.x - call dgemv('N',n,n,1d0,H,size(H,1),x,1,0d0,part_2a,1) - - ! Product 1/2 . x^T.H.x - part_2 = 0.5d0 * ddot(n,x,1,part_2a,1) - - !if (debug) then - print*,'1/2*x^T.H.x : ', part_2 - !endif - - print*,'prev_energy', prev_energy - - ! Sum - e_model = prev_energy + part_1 + part_2 - - ! Writing the predicted energy - print*, 'Predicted energy after the rotation : ', e_model - print*, 'Previous energy - predicted energy:', prev_energy - e_model - - ! Can be deleted, already in another subroutine - if (DABS(prev_energy - e_model) < 1d-12 ) then - print*,'WARNING: ABS(prev_energy - e_model) < 1d-12' - endif - - ! Deallocation - deallocate(part_2a) - - call wall_time(t2) - t3 = t2 - t1 - print*,'Time in trust e model:', t3 - - print*,'---End trust_e_model---' - print*,'' - -end subroutine diff --git a/src/utils_trust_region/trust_region_expected_e.org b/src/utils_trust_region/trust_region_expected_e.org deleted file mode 100644 index 58c8f804..00000000 --- a/src/utils_trust_region/trust_region_expected_e.org +++ /dev/null @@ -1,121 +0,0 @@ -* Predicted energy : e_model - -*Compute the energy predicted by the Taylor series* - -The energy is predicted using a Taylor expansion truncated at te 2nd -order : - -\begin{align*} -E_{k+1} = E_{k} + \textbf{g}_k^{T} \cdot \textbf{x}_{k+1} + \frac{1}{2} \cdot \textbf{x}_{k+1}^T \cdot \textbf{H}_{k} \cdot \textbf{x}_{k+1} + \mathcal{O}(\textbf{x}_{k+1}^2) -\end{align*} - -Input: -| n | integer | m*(m-1)/2 | -| v_grad(n) | double precision | gradient | -| H(n,n) | double precision | hessian | -| x(n) | double precision | Step in the trust region | -| prev_energy | double precision | previous energy | - -Output: -| e_model | double precision | predicted energy after the rotation of the MOs | - -Internal: -| part_1 | double precision | v_grad^T.x | -| part_2 | double precision | 1/2 . x^T.H.x | -| part_2a | double precision | H.x | -| i,j | integer | indexes | - -Function: -| ddot | double precision | dot product (Lapack) | - -#+BEGIN_SRC f90 :comments org :tangle trust_region_expected_e.irp.f -subroutine trust_region_expected_e(n,v_grad,H,x,prev_energy,e_model) - - include 'pi.h' - - BEGIN_DOC - ! Compute the expected criterion/energy after the application of the step x - END_DOC - - implicit none - - ! Variables - - ! in - integer, intent(in) :: n - double precision, intent(in) :: v_grad(n),H(n,n),x(n) - double precision, intent(in) :: prev_energy - - ! out - double precision, intent(out) :: e_model - - ! internal - double precision :: part_1, part_2, t1,t2,t3 - double precision, allocatable :: part_2a(:) - - integer :: i,j - - !Function - double precision :: ddot - - print*,'' - print*,'---Trust_e_model---' - - call wall_time(t1) - - ! Allocation - allocate(part_2a(n)) -#+END_SRC - -** Calculations - -part_1 corresponds to the product g.x -part_2a corresponds to the product H.x -part_2 corresponds to the product 0.5*(x^T.H.x) - -TODO: remove the dot products - -#+BEGIN_SRC f90 :comments org :tangle trust_region_expected_e.irp.f - ! Product v_grad.x - part_1 = ddot(n,v_grad,1,x,1) - - !if (debug) then - print*,'g.x : ', part_1 - !endif - - ! Product H.x - call dgemv('N',n,n,1d0,H,size(H,1),x,1,0d0,part_2a,1) - - ! Product 1/2 . x^T.H.x - part_2 = 0.5d0 * ddot(n,x,1,part_2a,1) - - !if (debug) then - print*,'1/2*x^T.H.x : ', part_2 - !endif - - print*,'prev_energy', prev_energy - - ! Sum - e_model = prev_energy + part_1 + part_2 - - ! Writing the predicted energy - print*, 'Predicted energy after the rotation : ', e_model - print*, 'Previous energy - predicted energy:', prev_energy - e_model - - ! Can be deleted, already in another subroutine - if (DABS(prev_energy - e_model) < 1d-12 ) then - print*,'WARNING: ABS(prev_energy - e_model) < 1d-12' - endif - - ! Deallocation - deallocate(part_2a) - - call wall_time(t2) - t3 = t2 - t1 - print*,'Time in trust e model:', t3 - - print*,'---End trust_e_model---' - print*,'' - -end subroutine -#+END_SRC diff --git a/src/utils_trust_region/trust_region_optimal_lambda.irp.f b/src/utils_trust_region/trust_region_optimal_lambda.irp.f deleted file mode 100644 index f71bb405..00000000 --- a/src/utils_trust_region/trust_region_optimal_lambda.irp.f +++ /dev/null @@ -1,1655 +0,0 @@ -! Newton's method to find the optimal lambda - -! *Compute the lambda value for the trust region* - -! This subroutine uses the Newton method in order to find the optimal -! lambda. This constant is added on the diagonal of the hessian to shift -! the eiganvalues. It has a double role: -! - ensure that the resulting hessian is positive definite for the -! Newton method -! - constrain the step in the trust region, i.e., -! $||\textbf{x}(\lambda)|| \leq \Delta$, where $\Delta$ is the radius -! of the trust region. -! We search $\lambda$ which minimizes -! \begin{align*} -! f(\lambda) = (||\textbf{x}_{(k+1)}(\lambda)||^2 -\Delta^2)^2 -! \end{align*} -! or -! \begin{align*} -! \tilde{f}(\lambda) = (\frac{1}{||\textbf{x}_{(k+1)}(\lambda)||^2}-\frac{1}{\Delta^2})^2 -! \end{align*} -! and gives obviously 0 in both cases. \newline - -! There are several cases: -! - If $\textbf{H}$ is positive definite the interval containing the -! solution is $\lambda \in (0, \infty)$ (and $-h_1 < 0$). -! - If $\textbf{H}$ is indefinite ($h_1 < 0$) and $\textbf{w}_1^T \cdot -! \textbf{g} \neq 0$ then the interval containing -! the solution is $\lambda \in (-h_1, \infty)$. -! - If $\textbf{H}$ is indefinite ($h_1 < 0$) and $\textbf{w}_1^T \cdot -! \textbf{g} = 0$ then the interval containing the solution is -! $\lambda \in (-h_1, \infty)$. The terms where $|h_i - \lambda| < -! 10^{-12}$ are not computed, so the term where $i = 1$ is -! automatically removed and this case becomes similar to the previous one. - -! So to avoid numerical problems (cf. trust_region) we start the -! algorithm at $\lambda=\max(0 + \epsilon,-h_1 + \epsilon)$, -! with $\epsilon$ a little constant. -! The research must be restricted to the interval containing the -! solution. For that reason a little trust region in 1D is used. - -! The Newton method to find the optimal $\lambda$ is : -! \begin{align*} -! \lambda_{(l+1)} &= \lambda_{(l)} - f^{''}(\lambda)_{(l)}^{-1} f^{'}(\lambda)_{(l)}^{} \\ -! \end{align*} -! $f^{'}(\lambda)_{(l)}$: the first derivative of $f$ with respect to -! $\lambda$ at the l-th iteration, -! $f^{''}(\lambda)_{(l)}$: the second derivative of $f$ with respect to -! $\lambda$ at the l-th iteration.\newline - -! Noting the Newton step $y = - f^{''}(\lambda)_{(l)}^{-1} -! f^{'}(\lambda)_{(l)}^{}$ we constrain $y$ such as -! \begin{align*} -! y \leq \alpha -! \end{align*} -! with $\alpha$ a scalar representing the trust length (trust region in -! 1D) where the function $f$ or $\tilde{f}$ is correctly describe by the -! Taylor series truncated at the second order. Thus, if $y > \alpha$, -! the constraint is applied as -! \begin{align*} -! y^* = \alpha \frac{y}{|y|} -! \end{align*} -! with $y^*$ the solution in the trust region. - -! The size of the trust region evolves in function of $\rho$ as for the -! trust region seen previously cf. trust_region, rho_model. -! The prediction of the value of $f$ or $\tilde{f}$ is done using the -! Taylor series truncated at the second order cf. "trust_region", -! "trust_e_model". - -! The first and second derivatives of $f(\lambda) = (||\textbf{x}(\lambda)||^2 - -! \Delta^2)^2$ with respect to $\lambda$ are: -! \begin{align*} -! \frac{\partial }{\partial \lambda} (||\textbf{x}(\lambda)||^2 - \Delta^2)^2 -! = 2 \left(\sum_{i=1}^n \frac{-2(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^3} \right) -! \left( - \Delta^2 + \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i+ \lambda)^2} \right) -! \end{align*} -! \begin{align*} -! \frac{\partial^2}{\partial \lambda^2} (||\textbf{x}(\lambda)||^2 - \Delta^2)^2 -! = 2 \left[ \left( \sum_{i=1}^n 6 \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^4} \right) \left( - \Delta^2 + \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^2} \right) + \left( \sum_{i=1}^n -2 \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^3} \right)^2 \right] -! \end{align*} - -! The first and second derivatives of $\tilde{f}(\lambda) = (1/||\textbf{x}(\lambda)||^2 - -! 1/\Delta^2)^2$ with respect to $\lambda$ are: -! \begin{align*} -! \frac{\partial}{\partial \lambda} (1/||\textbf{x}(\lambda)||^2 - 1/\Delta^2)^2 -! &= 4 \frac{\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3}} -! {(\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^3} -! - \frac{4}{\Delta^2} \frac{\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3)}} -! {(\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^2} \\ -! &= 4 \sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3} -! \left( \frac{1}{(\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^3} -! - \frac{1}{\Delta^2 (\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^2} \right) -! \end{align*} - -! \begin{align*} -! \frac{\partial^2}{\partial \lambda^2} (1/||\textbf{x}(\lambda)||^2 - 1/\Delta^2)^2 -! &= 4 \left[ \frac{(\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3)})^2} -! {(\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^4} -! - 3 \frac{\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^4}} -! {(\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^3} \right] \\ -! &- \frac{4}{\Delta^2} \left[ \frac{(\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2} -! {(h_i + \lambda)^3)})^2}{(\sum_ {i=1}^n\frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^3} -! - 3 \frac{\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^4}} -! {(\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^2} \right] -! \end{align*} - -! Provided in qp_edit: -! | thresh_rho_2 | -! | thresh_cc | -! | nb_it_max_lambda | -! | version_lambda_search | -! | nb_it_max_pre_search | -! see qp_edit for more details - -! Input: -! | n | integer | m*(m-1)/2 | -! | e_val(n) | double precision | eigenvalues of the hessian | -! | tmp_wtg(n) | double precision | w_i^T.v_grad(i) | -! | delta | double precision | delta for the trust region | - -! Output: -! | lambda | double precision | Lagrange multiplier to constrain the norm of the size of the Newton step | -! | | | lambda > 0 | - -! Internal: -! | d1_N | double precision | value of d1_norm_trust_region | -! | d2_N | double precision | value of d2_norm_trust_region | -! | f_N | double precision | value of f_norm_trust_region | -! | prev_f_N | double precision | previous value of f_norm_trust_region | -! | f_R | double precision | (norm(x)^2 - delta^2)^2 or (1/norm(x)^2 - 1/delta^2)^2 | -! | prev_f_R | double precision | previous value of f_R | -! | model | double precision | predicted value of f_R from prev_f_R and y | -! | d_1 | double precision | value of the first derivative | -! | d_2 | double precision | value of the second derivative | -! | y | double precision | Newton's step, y = -f''^-1 . f' = lambda - prev_lambda | -! | prev_lambda | double precision | previous value of lambda | -! | t1,t2,t3 | double precision | wall time | -! | i | integer | index | -! | epsilon | double precision | little constant to avoid numerical problem | -! | rho_2 | double precision | (prev_f_R - f_R)/(prev_f_R - model), agreement between model and f_R | -! | version | integer | version of the root finding method | - -! Function: -! | d1_norm_trust_region | double precision | first derivative with respect to lambda of (norm(x)^2 - Delta^2)^2 | -! | d2_norm_trust_region | double precision | first derivative with respect to lambda of (norm(x)^2 - Delta^2)^2 | -! | d1_norm_inverse_trust_region | double precision | first derivative with respect to lambda of (1/norm(x)^2 - 1/Delta^2)^2 | -! | d2_norm_inverse_trust_region | double precision | second derivative with respect to lambda of (1/norm(x)^2 - 1/Delta^2)^2 | -! | f_norm_trust_region | double precision | value of norm(x)^2 | - - - -subroutine trust_region_optimal_lambda(n,e_val,tmp_wtg,delta,lambda) - - include 'pi.h' - - BEGIN_DOC - ! Research the optimal lambda to constrain the step size in the trust region - END_DOC - - implicit none - - ! Variables - - ! in - integer, intent(in) :: n - double precision, intent(inout) :: e_val(n) - double precision, intent(in) :: delta - double precision, intent(in) :: tmp_wtg(n) - - ! out - double precision, intent(out) :: lambda - - ! Internal - double precision :: d1_N, d2_N, f_N, prev_f_N - double precision :: prev_f_R, f_R - double precision :: model - double precision :: d_1, d_2 - double precision :: t1,t2,t3 - integer :: i - double precision :: epsilon - double precision :: y - double precision :: prev_lambda - double precision :: rho_2 - double precision :: alpha - integer :: version - - ! Functions - double precision :: d1_norm_trust_region,d1_norm_trust_region_omp - double precision :: d2_norm_trust_region, d2_norm_trust_region_omp - double precision :: f_norm_trust_region, f_norm_trust_region_omp - double precision :: d1_norm_inverse_trust_region - double precision :: d2_norm_inverse_trust_region - double precision :: d1_norm_inverse_trust_region_omp - double precision :: d2_norm_inverse_trust_region_omp - - print*,'' - print*,'---Trust_newton---' - print*,'' - - call wall_time(t1) - - ! version_lambda_search - ! 1 -> ||x||^2 - delta^2 = 0, - ! 2 -> 1/||x||^2 - 1/delta^2 = 0 (better) - if (version_lambda_search == 1) then - print*, 'Research of the optimal lambda by solving ||x||^2 - delta^2 = 0' - else - print*, 'Research of the optimal lambda by solving 1/||x||^2 - 1/delta^2 = 0' - endif - ! Version 2 is normally better - - - -! Resolution with the Newton method: - - -! Initialization - epsilon = 1d-4 - lambda =MAX(0d0, -e_val(1)) - - ! Pre research of lambda to start near the optimal lambda - ! by adding a constant epsilon and changing the constant to - ! have ||x(lambda + epsilon)|| ~ delta, before setting - ! lambda = lambda + epsilon - print*, 'Pre research of lambda:' - print*,'Initial lambda =', lambda - f_N = f_norm_trust_region_omp(n,e_val,tmp_wtg,lambda + epsilon) - print*,'||x(lambda)||=', dsqrt(f_N),'delta=',delta - i = 1 - - ! To increase lambda - if (f_N > delta**2) then - print*,'Increasing lambda...' - do while (f_N > delta**2 .and. i <= nb_it_max_pre_search) - - ! Update the previous norm - prev_f_N = f_N - ! New epsilon - epsilon = epsilon * 2d0 - ! New norm - f_N = f_norm_trust_region_omp(n,e_val,tmp_wtg,lambda + epsilon) - - print*, 'lambda', lambda + epsilon, '||x||', dsqrt(f_N), 'delta', delta - - ! Security - if (prev_f_N < f_N) then - print*,'WARNING, error: prev_f_N < f_N, exit' - epsilon = epsilon * 0.5d0 - i = nb_it_max_pre_search + 1 - endif - - i = i + 1 - enddo - - ! To reduce lambda - else - print*,'Reducing lambda...' - do while (f_N < delta**2 .and. i <= nb_it_max_pre_search) - - ! Update the previous norm - prev_f_N = f_N - ! New epsilon - epsilon = epsilon * 0.5d0 - ! New norm - f_N = f_norm_trust_region_omp(n,e_val,tmp_wtg,lambda + epsilon) - - print*, 'lambda', lambda + epsilon, '||x||', dsqrt(f_N), 'delta', delta - - ! Security - if (prev_f_N > f_N) then - print*,'WARNING, error: prev_f_N > f_N, exit' - epsilon = epsilon * 2d0 - i = nb_it_max_pre_search + 1 - endif - - i = i + 1 - enddo - endif - - print*,'End of the pre research of lambda' - - ! New value of lambda - lambda = lambda + epsilon - - print*, 'e_val(1):', e_val(1) - print*, 'Staring point, lambda =', lambda - - ! thresh_cc, threshold for the research of the optimal lambda - ! Leaves the loop when ABS(1d0-||x||^2/delta^2) > thresh_cc - ! thresh_rho_2, threshold to cancel the step in the research - ! of the optimal lambda, the step is cancelled if rho_2 < thresh_rho_2 - print*,'Threshold for the CC:', thresh_cc - print*,'Threshold for rho_2:', thresh_rho_2 - - print*, 'w_1^T . g =', tmp_wtg(1) - - ! Debug - !if (debug) then - ! print*, 'Iteration rho_2 lambda delta ||x|| |1-(||x||^2/delta^2)|' - !endif - - ! Initialization - i = 1 - f_N = f_norm_trust_region_omp(n,e_val,tmp_wtg,lambda) ! Value of the ||x(lambda)||^2 - model = 0d0 ! predicted value of (||x||^2 - delta^2)^2 - prev_f_N = 0d0 ! previous value of ||x||^2 - prev_f_R = 0d0 ! previous value of (||x||^2 - delta^2)^2 - f_R = 0d0 ! value of (||x||^2 - delta^2)^2 - rho_2 = 0d0 ! (prev_f_R - f_R)/(prev_f_R - m) - y = 0d0 ! step size - prev_lambda = 0d0 ! previous lambda - - ! Derivatives - if (version_lambda_search == 1) then - d_1 = d1_norm_trust_region_omp(n,e_val,tmp_wtg,lambda,delta) ! first derivative of (||x(lambda)||^2 - delta^2)^2 - d_2 = d2_norm_trust_region_omp(n,e_val,tmp_wtg,lambda,delta) ! second derivative of (||x(lambda)||^2 - delta^2)^2 - else - d_1 = d1_norm_inverse_trust_region_omp(n,e_val,tmp_wtg,lambda,delta) ! first derivative of (1/||x(lambda)||^2 - 1/delta^2)^2 - d_2 = d2_norm_inverse_trust_region_omp(n,e_val,tmp_wtg,lambda,delta) ! second derivative of (1/||x(lambda)||^2 - 1/delta^2)^2 - endif - - ! Trust length - alpha = DABS((1d0/d_2)*d_1) - - ! Newton's method - do while (i <= 100 .and. DABS(1d0-f_N/delta**2) > thresh_cc) - print*,'--------------------------------------' - print*,'Research of lambda, iteration:', i - print*,'--------------------------------------' - - ! Update of f_N, f_R and the derivatives - prev_f_N = f_N - if (version_lambda_search == 1) then - prev_f_R = (prev_f_N - delta**2)**2 - d_1 = d1_norm_trust_region_omp(n,e_val,tmp_wtg,lambda,delta) ! first derivative of (||x(lambda)||^2 - delta^2)^2 - d_2 = d2_norm_trust_region_omp(n,e_val,tmp_wtg,lambda,delta) ! second derivative of (||x(lambda)||^2 - delta^2)^2 - else - prev_f_R = (1d0/prev_f_N - 1d0/delta**2)**2 - d_1 = d1_norm_inverse_trust_region_omp(n,e_val,tmp_wtg,lambda,delta) ! first derivative of (1/||x(lambda)||^2 - 1/delta^2)^2 - d_2 = d2_norm_inverse_trust_region_omp(n,e_val,tmp_wtg,lambda,delta) ! second derivative of (1/||x(lambda)||^2 - 1/delta^2)^2 - endif - write(*,'(a,E12.5,a,E12.5)') ' 1st and 2nd derivative: ', d_1,', ', d_2 - - ! Newton's step - y = -(1d0/DABS(d_2))*d_1 - - ! Constraint on y (the newton step) - if (DABS(y) > alpha) then - y = alpha * (y/DABS(y)) ! preservation of the sign of y - endif - write(*,'(a,E12.5)') ' Step length: ', y - - ! Predicted value of (||x(lambda)||^2 - delta^2)^2, Taylor series - model = prev_f_R + d_1 * y + 0.5d0 * d_2 * y**2 - - ! Updates lambda - prev_lambda = lambda - lambda = prev_lambda + y - print*,'prev lambda:', prev_lambda - print*,'new lambda:', lambda - - ! Checks if lambda is in (-h_1, \infty) - if (lambda > MAX(0d0, -e_val(1))) then - ! New value of ||x(lambda)||^2 - f_N = f_norm_trust_region_omp(n,e_val,tmp_wtg,lambda) - - ! New f_R - if (version_lambda_search == 1) then - f_R = (f_N - delta**2)**2 ! new value of (||x(lambda)||^2 - delta^2)^2 - else - f_R = (1d0/f_N - 1d0/delta**2)**2 ! new value of (1/||x(lambda)||^2 -1/delta^2)^2 - endif - - if (version_lambda_search == 1) then - print*,'Previous value of (||x(lambda)||^2 - delta^2)^2:', prev_f_R - print*,'Actual value of (||x(lambda)||^2 - delta^2)^2:', f_R - print*,'Predicted value of (||x(lambda)||^2 - delta^2)^2:', model - else - print*,'Previous value of (1/||x(lambda)||^2 - 1/delta^2)^2:', prev_f_R - print*,'Actual value of (1/||x(lambda)||^2 - 1/delta^2)^2:', f_R - print*,'Predicted value of (1/||x(lambda)||^2 - 1/delta^2)^2:', model - endif - - print*,'previous - actual:', prev_f_R - f_R - print*,'previous - model:', prev_f_R - model - - ! Check the gain - if (DABS(prev_f_R - model) < thresh_model_2) then - print*,'' - print*,'WARNING: ABS(previous - model) <', thresh_model_2, 'rho_2 will tend toward infinity' - print*,'' - endif - - ! Will be deleted - !if (prev_f_R - f_R <= 1d-16 .or. prev_f_R - model <= 1d-16) then - ! print*,'' - ! print*,'WARNING: ABS(previous - model) <= 1d-16, exit' - ! print*,'' - ! exit - !endif - - ! Computes rho_2 - rho_2 = (prev_f_R - f_R)/(prev_f_R - model) - print*,'rho_2:', rho_2 - else - rho_2 = 0d0 ! in order to reduce the size of the trust region, alpha, until lambda is in (-h_1, \infty) - print*,'lambda < -e_val(1) ===> rho_2 = 0' - endif - - ! Evolution of the trust length, alpha - if (rho_2 >= 0.75d0) then - alpha = 2d0 * alpha - elseif (rho_2 >= 0.5d0) then - alpha = alpha - elseif (rho_2 >= 0.25d0) then - alpha = 0.5d0 * alpha - else - alpha = 0.25d0 * alpha - endif - write(*,'(a,E12.5)') ' New trust length alpha: ', alpha - - ! cancellaion of the step if rho < 0.1 - if (rho_2 < thresh_rho_2) then !0.1d0) then - lambda = prev_lambda - f_N = prev_f_N - print*,'Rho_2 <', thresh_rho_2,', cancellation of the step: lambda = prev_lambda' - endif - - print*,'' - print*,'lambda, ||x||, delta:' - print*, lambda, dsqrt(f_N), delta - print*,'CC:', DABS(1d0 - f_N/delta**2) - print*,'' - - i = i + 1 - enddo - - ! if trust newton failed - if (i > nb_it_max_lambda) then - print*,'' - print*,'######################################################' - print*,'WARNING: i >', nb_it_max_lambda,'for the trust Newton' - print*,'The research of the optimal lambda has failed' - print*,'######################################################' - print*,'' - endif - - print*,'Number of iterations :', i - print*,'Value of lambda :', lambda - print*,'Error on the trust region (1d0-f_N/delta**2) (Convergence criterion) :', 1d0-f_N/delta**2 - print*,'Error on the trust region (||x||^2 - delta^2)^2) :', (f_N - delta**2)**2 - print*,'Error on the trust region (1/||x||^2 - 1/delta^2)^2)', (1d0/f_N - 1d0/delta**2)**2 - - ! Time - call wall_time(t2) - t3 = t2 - t1 - print*,'Time in trust_newton:', t3 - - print*,'' - print*,'---End trust_newton---' - print*,'' - -end subroutine - -! OMP: First derivative of (||x||^2 - Delta^2)^2 - -! *Function to compute the first derivative of (||x||^2 - Delta^2)^2* - -! This function computes the first derivative of (||x||^2 - Delta^2)^2 -! with respect to lambda. - -! \begin{align*} -! \frac{\partial }{\partial \lambda} (||\textbf{x}(\lambda)||^2 - \Delta^2)^2 -! = -4 \left(\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3} \right) -! \left( - \Delta^2 + \sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i+ \lambda)^2} \right) -! \end{align*} - -! \begin{align*} -! \text{accu1} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2} \\ -! \text{accu2} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3} -! \end{align*} - -! Provided: -! | mo_num | integer | number of MOs | - -! Input: -! | n | integer | mo_num*(mo_num-1)/2 | -! | e_val(n) | double precision | eigenvalues of the hessian | -! | W(n,n) | double precision | eigenvectors of the hessian | -! | v_grad(n) | double precision | gradient | -! | lambda | double precision | Lagrange multiplier | -! | delta | double precision | Delta of the trust region | - -! Internal: -! | accu1 | double precision | first sum of the formula | -! | accu2 | double precision | second sum of the formula | -! | tmp_accu1 | double precision | temporary array for the first sum | -! | tmp_accu2 | double precision | temporary array for the second sum | -! | tmp_wtg(n) | double precision | temporary array for W^t.v_grad | -! | i,j | integer | indexes | - -! Function: -! | d1_norm_trust_region | double precision | first derivative with respect to lambda of (norm(x)^2 - Delta^2)^2 | - - -function d1_norm_trust_region_omp(n,e_val,tmp_wtg,lambda,delta) - - use omp_lib - include 'pi.h' - - BEGIN_DOC - ! Compute the first derivative with respect to lambda of (||x(lambda)||^2 - Delta^2)^2 - END_DOC - - implicit none - - ! in - integer, intent(in) :: n - double precision, intent(in) :: e_val(n) - double precision, intent(in) :: tmp_wtg(n) - double precision, intent(in) :: lambda - double precision, intent(in) :: delta - - ! Internal - double precision :: wtg,accu1,accu2 - integer :: i,j - double precision, allocatable :: tmp_accu1(:), tmp_accu2(:) - - ! Functions - double precision :: d1_norm_trust_region_omp - - ! Allocation - allocate(tmp_accu1(n), tmp_accu2(n)) - - ! OMP - call omp_set_max_active_levels(1) - - ! OMP - !$OMP PARALLEL & - !$OMP PRIVATE(i,j) & - !$OMP SHARED(n,lambda, e_val, thresh_eig,& - !$OMP tmp_accu1, tmp_accu2, tmp_wtg, accu1,accu2) & - !$OMP DEFAULT(NONE) - - !$OMP MASTER - accu1 = 0d0 - accu2 = 0d0 - !$OMP END MASTER - - !$OMP DO - do i = 1, n - tmp_accu1(i) = 0d0 - enddo - !$OMP END DO - - !$OMP DO - do i = 1, n - tmp_accu2(i) = 0d0 - enddo - !$OMP END DO - - !$OMP DO - do i = 1, n - if (ABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - tmp_accu1(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**2 - endif - enddo - !$OMP END DO - - !$OMP MASTER - do i = 1, n - accu1 = accu1 + tmp_accu1(i) - enddo - !$OMP END MASTER - - !$OMP DO - do i = 1, n - if (ABS(e_val(i)) > thresh_eig) then - tmp_accu2(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**3 - endif - enddo - !$OMP END DO - - !$OMP MASTER - do i = 1, n - accu2 = accu2 + tmp_accu2(i) - enddo - !$OMP END MASTER - - !$OMP END PARALLEL - - call omp_set_max_active_levels(4) - - d1_norm_trust_region_omp = -4d0 * accu2 * (accu1 - delta**2) - - deallocate(tmp_accu1, tmp_accu2) - -end function - -! OMP: Second derivative of (||x||^2 - Delta^2)^2 - -! *Function to compute the second derivative of (||x||^2 - Delta^2)^2* - -! This function computes the second derivative of (||x||^2 - Delta^2)^2 -! with respect to lambda. -! \begin{align*} -! \frac{\partial^2 }{\partial \lambda^2} (||\textbf{x}(\lambda)||^2 - \Delta^2)^2 -! = 2 \left[ \left( \sum_{i=1}^n 6 \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^4} \right) \left( - \Delta^2 + \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^2} \right) + \left( \sum_{i=1}^n -2 \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^3} \right)^2 \right] -! \end{align*} - -! \begin{align*} -! \text{accu1} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^2} \\ -! \text{accu2} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^3} \\ -! \text{accu3} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^4} -! \end{align*} - -! Provided: -! | m_num | integer | number of MOs | - -! Input: -! | n | integer | mo_num*(mo_num-1)/2 | -! | e_val(n) | double precision | eigenvalues of the hessian | -! | W(n,n) | double precision | eigenvectors of the hessian | -! | v_grad(n) | double precision | gradient | -! | lambda | double precision | Lagrange multiplier | -! | delta | double precision | Delta of the trust region | - -! Internal: -! | accu1 | double precision | first sum of the formula | -! | accu2 | double precision | second sum of the formula | -! | accu3 | double precision | third sum of the formula | -! | tmp_accu1 | double precision | temporary array for the first sum | -! | tmp_accu2 | double precision | temporary array for the second sum | -! | tmp_accu2 | double precision | temporary array for the third sum | -! | tmp_wtg(n) | double precision | temporary array for W^t.v_grad | -! | i,j | integer | indexes | - -! Function: -! | d2_norm_trust_region | double precision | second derivative with respect to lambda of (norm(x)^2 - Delta^2)^2 | - - -function d2_norm_trust_region_omp(n,e_val,tmp_wtg,lambda,delta) - - use omp_lib - include 'pi.h' - - BEGIN_DOC - ! Compute the second derivative with respect to lambda of (||x(lambda)||^2 - Delta^2)^2 - END_DOC - - implicit none - - ! Variables - - ! in - integer, intent(in) :: n - double precision, intent(in) :: e_val(n) - double precision, intent(in) :: tmp_wtg(n) - double precision, intent(in) :: lambda - double precision, intent(in) :: delta - - ! Functions - double precision :: d2_norm_trust_region_omp - double precision :: ddot - - ! Internal - double precision :: accu1,accu2,accu3 - double precision, allocatable :: tmp_accu1(:), tmp_accu2(:), tmp_accu3(:) - integer :: i, j - - ! Allocation - allocate(tmp_accu1(n), tmp_accu2(n), tmp_accu3(n)) - - call omp_set_max_active_levels(1) - - ! OMP - !$OMP PARALLEL & - !$OMP PRIVATE(i,j) & - !$OMP SHARED(n,lambda, e_val, thresh_eig,& - !$OMP tmp_accu1, tmp_accu2, tmp_accu3, tmp_wtg, & - !$OMP accu1, accu2, accu3) & - !$OMP DEFAULT(NONE) - - ! Initialization - - !$OMP MASTER - accu1 = 0d0 - accu2 = 0d0 - accu3 = 0d0 - !$OMP END MASTER - - !$OMP DO - do i = 1, n - tmp_accu1(i) = 0d0 - enddo - !$OMP END DO - !$OMP DO - do i = 1, n - tmp_accu2(i) = 0d0 - enddo - !$OMP END DO - !$OMP DO - do i = 1, n - tmp_accu3(i) = 0d0 - enddo - !$OMP END DO - - ! Calculations - - ! accu1 - !$OMP DO - do i = 1, n - if (ABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - tmp_accu1(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**2 - endif - enddo - !$OMP END DO - - !$OMP MASTER - do i = 1, n - accu1 = accu1 + tmp_accu1(i) - enddo - !$OMP END MASTER - - ! accu2 - !$OMP DO - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - tmp_accu2(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**3 - endif - enddo - !$OMP END DO - - ! accu3 - !$OMP MASTER - do i = 1, n - accu2 = accu2 + tmp_accu2(i) - enddo - !$OMP END MASTER - - !$OMP DO - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - tmp_accu3(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**4 - endif - enddo - !$OMP END DO - - !$OMP MASTER - do i = 1, n - accu3 = accu3 + tmp_accu3(i) - enddo - !$OMP END MASTER - - !$OMP END PARALLEL - - d2_norm_trust_region_omp = 2d0 * (6d0 * accu3 * (- delta**2 + accu1) + (-2d0 * accu2)**2) - - deallocate(tmp_accu1, tmp_accu2, tmp_accu3) - -end function - -! OMP: Function value of ||x||^2 - -! *Compute the value of ||x||^2* - -! This function computes the value of ||x(lambda)||^2 - -! \begin{align*} -! ||\textbf{x}(\lambda)||^2 = \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^2} -! \end{align*} - -! Provided: -! | m_num | integer | number of MOs | - -! Input: -! | n | integer | mo_num*(mo_num-1)/2 | -! | e_val(n) | double precision | eigenvalues of the hessian | -! | W(n,n) | double precision | eigenvectors of the hessian | -! | v_grad(n) | double precision | gradient | -! | lambda | double precision | Lagrange multiplier | - -! Internal: -! | tmp_wtg(n) | double precision | temporary array for W^T.v_grad | -! | tmp_fN | double precision | temporary array for the function | -! | i,j | integer | indexes | - - -function f_norm_trust_region_omp(n,e_val,tmp_wtg,lambda) - - use omp_lib - - include 'pi.h' - - BEGIN_DOC - ! Compute ||x(lambda)||^2 - END_DOC - - implicit none - - ! Variables - - ! in - integer, intent(in) :: n - double precision, intent(in) :: e_val(n) - double precision, intent(in) :: tmp_wtg(n) - double precision, intent(in) :: lambda - - ! functions - double precision :: f_norm_trust_region_omp - - ! internal - double precision, allocatable :: tmp_fN(:) - integer :: i,j - - ! Allocation - allocate(tmp_fN(n)) - - call omp_set_max_active_levels(1) - - ! OMP - !$OMP PARALLEL & - !$OMP PRIVATE(i,j) & - !$OMP SHARED(n,lambda, e_val, thresh_eig,& - !$OMP tmp_fN, tmp_wtg, f_norm_trust_region_omp) & - !$OMP DEFAULT(NONE) - - ! Initialization - - !$OMP MASTER - f_norm_trust_region_omp = 0d0 - !$OMP END MASTER - - !$OMP DO - do i = 1, n - tmp_fN(i) = 0d0 - enddo - !$OMP END DO - - ! Calculations - !$OMP DO - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - tmp_fN(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**2 - endif - enddo - !$OMP END DO - - !$OMP MASTER - do i = 1, n - f_norm_trust_region_omp = f_norm_trust_region_omp + tmp_fN(i) - enddo - !$OMP END MASTER - - !$OMP END PARALLEL - - deallocate(tmp_fN) - -end function - -! First derivative of (||x||^2 - Delta^2)^2 -! Version without omp - -! *Function to compute the first derivative of ||x||^2 - Delta* - -! This function computes the first derivative of (||x||^2 - Delta^2)^2 -! with respect to lambda. - -! \begin{align*} -! \frac{\partial }{\partial \lambda} (||\textbf{x}(\lambda)||^2 - \Delta^2)^2 -! = 2 \left(-2\sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^3} \right) -! \left( - \Delta^2 + \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i+ \lambda)^2} \right) -! \end{align*} - -! \begin{align*} -! \text{accu1} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^2} \\ -! \text{accu2} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^3} -! \end{align*} - -! Provided: -! | m_num | integer | number of MOs | - -! Input: -! | n | integer | mo_num*(mo_num-1)/2 | -! | e_val(n) | double precision | eigenvalues of the hessian | -! | W(n,n) | double precision | eigenvectors of the hessian | -! | v_grad(n) | double precision | gradient | -! | lambda | double precision | Lagrange multiplier | -! | delta | double precision | Delta of the trust region | - -! Internal: -! | accu1 | double precision | first sum of the formula | -! | accu2 | double precision | second sum of the formula | -! | wtg | double precision | temporary variable to store W^T.v_grad | -! | i,j | integer | indexes | - -! Function: -! | d1_norm_trust_region | double precision | first derivative with respect to lambda of (norm(x)^2 - Delta^2)^2 | -! | ddot | double precision | blas dot product | - - -function d1_norm_trust_region(n,e_val,w,v_grad,lambda,delta) - - include 'pi.h' - - BEGIN_DOC - ! Compute the first derivative with respect to lambda of (||x(lambda)||^2 - Delta^2)^2 - END_DOC - - implicit none - - ! Variables - - ! in - integer, intent(in) :: n - double precision, intent(in) :: e_val(n) - double precision, intent(in) :: w(n,n) - double precision, intent(in) :: v_grad(n) - double precision, intent(in) :: lambda - double precision, intent(in) :: delta - - ! Internal - double precision :: wtg, accu1, accu2 - integer :: i, j - - ! Functions - double precision :: d1_norm_trust_region - double precision :: ddot - - ! Initialization - accu1 = 0d0 - accu2 = 0d0 - - do i = 1, n - wtg = 0d0 - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - do j = 1, n - wtg = wtg + w(j,i) * v_grad(j) - enddo - !wtg = ddot(n,w(:,i),1,v_grad,1) - accu1 = accu1 + wtg**2 / (e_val(i) + lambda)**2 - endif - enddo - - do i = 1, n - wtg = 0d0 - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - do j = 1, n - wtg = wtg + w(j,i) * v_grad(j) - enddo - !wtg = ddot(n,w(:,i),1,v_grad,1) - accu2 = accu2 - 2d0 * wtg**2 / (e_val(i) + lambda)**3 - endif - enddo - - d1_norm_trust_region = 2d0 * accu2 * (accu1 - delta**2) - -end function - -! Second derivative of (||x||^2 - Delta^2)^2 -! Version without OMP - -! *Function to compute the second derivative of ||x||^2 - Delta* - - -! \begin{equation} -! \frac{\partial^2 }{\partial \lambda^2} (||\textbf{x}(\lambda)||^2 - \Delta^2)^2 -! = 2 \left[ \left( \sum_{i=1}^n 6 \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^4} \right) \left( - \Delta^2 + \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^2} \right) + \left( \sum_{i=1}^n -2 \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^3} \right)^2 \right] -! \end{equation} - -! \begin{align*} -! \text{accu1} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^2} \\ -! \text{accu2} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^3} \\ -! \text{accu3} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^4} -! \end{align*} -! Provided: -! | m_num | integer | number of MOs | - -! Input: -! | n | integer | mo_num*(mo_num-1)/2 | -! | e_val(n) | double precision | eigenvalues of the hessian | -! | W(n,n) | double precision | eigenvectors of the hessian | -! | v_grad(n) | double precision | gradient | -! | lambda | double precision | Lagrange multiplier | -! | delta | double precision | Delta of the trust region | - -! Internal: -! | accu1 | double precision | first sum of the formula | -! | accu2 | double precision | second sum of the formula | -! | accu3 | double precision | third sum of the formula | -! | wtg | double precision | temporary variable to store W^T.v_grad | -! | i,j | integer | indexes | - -! Function: -! | d2_norm_trust_region | double precision | second derivative with respect to lambda of norm(x)^2 - Delta^2 | -! | ddot | double precision | blas dot product | - - -function d2_norm_trust_region(n,e_val,w,v_grad,lambda,delta) - - include 'pi.h' - - BEGIN_DOC - ! Compute the second derivative with respect to lambda of (||x(lambda)||^2 - Delta^2)^2 - END_DOC - - implicit none - - ! Variables - - ! in - integer, intent(in) :: n - double precision, intent(in) :: e_val(n) - double precision, intent(in) :: w(n,n) - double precision, intent(in) :: v_grad(n) - double precision, intent(in) :: lambda - double precision, intent(in) :: delta - - ! Functions - double precision :: d2_norm_trust_region - double precision :: ddot - - ! Internal - double precision :: wtg,accu1,accu2,accu3 - integer :: i, j - - ! Initialization - accu1 = 0d0 - accu2 = 0d0 - accu3 = 0d0 - - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - wtg = 0d0 - do j = 1, n - wtg = wtg + w(j,i) * v_grad(j) - enddo - !wtg = ddot(n,w(:,i),1,v_grad,1) - accu1 = accu1 + wtg**2 / (e_val(i) + lambda)**2 !4 - endif - enddo - - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - wtg = 0d0 - do j = 1, n - wtg = wtg + w(j,i) * v_grad(j) - enddo - !wtg = ddot(n,w(:,i),1,v_grad,1) - accu2 = accu2 - 2d0 * wtg**2 / (e_val(i) + lambda)**3 !2 - endif - enddo - - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - wtg = 0d0 - do j = 1, n - wtg = wtg + w(j,i) * v_grad(j) - enddo - !wtg = ddot(n,w(:,i),1,v_grad,1) - accu3 = accu3 + 6d0 * wtg**2 / (e_val(i) + lambda)**4 !3 - endif - enddo - - d2_norm_trust_region = 2d0 * (accu3 * (- delta**2 + accu1) + accu2**2) - -end function - -! Function value of ||x||^2 -! Version without OMP - -! *Compute the value of ||x||^2* - -! This function computes the value of ||x(lambda)||^2 - -! \begin{align*} -! ||\textbf{x}(\lambda)||^2 = \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^2} -! \end{align*} - -! Provided: -! | m_num | integer | number of MOs | - -! Input: -! | n | integer | mo_num*(mo_num-1)/2 | -! | e_val(n) | double precision | eigenvalues of the hessian | -! | W(n,n) | double precision | eigenvectors of the hessian | -! | v_grad(n) | double precision | gradient | -! | lambda | double precision | Lagrange multiplier | -! | delta | double precision | Delta of the trust region | - -! Internal: -! | wtg | double precision | temporary variable to store W^T.v_grad | -! | i,j | integer | indexes | - -! Function: -! | f_norm_trust_region | double precision | value of norm(x)^2 | -! | ddot | double precision | blas dot product | - - - -function f_norm_trust_region(n,e_val,tmp_wtg,lambda) - - include 'pi.h' - - BEGIN_DOC - ! Compute ||x(lambda)||^2 - END_DOC - - implicit none - - ! Variables - - ! in - integer, intent(in) :: n - double precision, intent(in) :: e_val(n) - double precision, intent(in) :: tmp_wtg(n) - double precision, intent(in) :: lambda - - ! function - double precision :: f_norm_trust_region - double precision :: ddot - - ! internal - integer :: i,j - - ! Initialization - f_norm_trust_region = 0d0 - - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - f_norm_trust_region = f_norm_trust_region + tmp_wtg(i)**2 / (e_val(i) + lambda)**2 - endif - enddo - -end function - -! OMP: First derivative of (1/||x||^2 - 1/Delta^2)^2 -! Version with OMP - -! *Compute the first derivative of (1/||x||^2 - 1/Delta^2)^2* - -! This function computes the value of (1/||x(lambda)||^2 - 1/Delta^2)^2 - -! \begin{align*} -! \frac{\partial}{\partial \lambda} (1/||\textbf{x}(\lambda)||^2 - 1/\Delta^2)^2 -! &= 4 \frac{\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3}} -! {(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^3} -! - \frac{4}{\Delta^2} \frac{\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3)}} -! {(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^2} \\ -! &= 4 \sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3} -! \left( \frac{1}{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^3} -! - \frac{1}{\Delta^2 (\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^2} \right) -! \end{align*} - -! \begin{align*} -! \text{accu1} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^2} \\ -! \text{accu2} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^3} -! \end{align*} - -! Provided: -! | m_num | integer | number of MOs | - -! Input: -! | n | integer | mo_num*(mo_num-1)/2 | -! | e_val(n) | double precision | eigenvalues of the hessian | -! | W(n,n) | double precision | eigenvectors of the hessian | -! | v_grad(n) | double precision | gradient | -! | lambda | double precision | Lagrange multiplier | -! | delta | double precision | Delta of the trust region | - -! Internal: -! | wtg | double precision | temporary variable to store W^T.v_grad | -! | tmp_accu1 | double precision | temporary array for the first sum | -! | tmp_accu2 | double precision | temporary array for the second sum | -! | tmp_wtg(n) | double precision | temporary array for W^t.v_grad | -! | i,j | integer | indexes | - -! Function: -! | d1_norm_inverse_trust_region | double precision | value of the first derivative | - - -function d1_norm_inverse_trust_region_omp(n,e_val,tmp_wtg,lambda,delta) - - use omp_lib - include 'pi.h' - - BEGIN_DOC - ! Compute the first derivative of (1/||x||^2 - 1/Delta^2)^2 - END_DOC - - implicit none - - ! Variables - - ! in - integer, intent(in) :: n - double precision, intent(in) :: e_val(n) - double precision, intent(in) :: tmp_wtg(n) - double precision, intent(in) :: lambda - double precision, intent(in) :: delta - - ! Internal - double precision :: accu1, accu2 - integer :: i,j - double precision, allocatable :: tmp_accu1(:), tmp_accu2(:) - - ! Functions - double precision :: d1_norm_inverse_trust_region_omp - - ! Allocation - allocate(tmp_accu1(n), tmp_accu2(n)) - - ! OMP - call omp_set_max_active_levels(1) - - ! OMP - !$OMP PARALLEL & - !$OMP PRIVATE(i,j) & - !$OMP SHARED(n,lambda, e_val, thresh_eig,& - !$OMP tmp_accu1, tmp_accu2, tmp_wtg, accu1, accu2) & - !$OMP DEFAULT(NONE) - - !$OMP MASTER - accu1 = 0d0 - accu2 = 0d0 - !$OMP END MASTER - - !$OMP DO - do i = 1, n - tmp_accu1(i) = 0d0 - enddo - !$OMP END DO - - !$OMP DO - do i = 1, n - tmp_accu2(i) = 0d0 - enddo - !$OMP END DO - -! !$OMP MASTER -! do i = 1, n -! if (ABS(e_val(i)+lambda) > 1d-12) then -! tmp_accu1(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**2 -! endif -! enddo -! !$OMP END MASTER - - !$OMP DO - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - tmp_accu1(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**2 - endif - enddo - !$OMP END DO - - !$OMP MASTER - do i = 1, n - accu1 = accu1 + tmp_accu1(i) - enddo - !$OMP END MASTER - -! !$OMP MASTER -! do i = 1, n -! if (ABS(e_val(i)+lambda) > 1d-12) then -! tmp_accu2(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**3 -! endif -! enddo -! !$OMP END MASTER - - !$OMP DO - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - tmp_accu2(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**3 - endif - enddo - !$OMP END DO - - !$OMP MASTER - do i = 1, n - accu2 = accu2 + tmp_accu2(i) - enddo - !$OMP END MASTER - - !$OMP END PARALLEL - - call omp_set_max_active_levels(4) - - d1_norm_inverse_trust_region_omp = 4d0 * accu2 * (1d0/accu1**3 - 1d0/(delta**2 * accu1**2)) - - deallocate(tmp_accu1, tmp_accu2) - -end - -! OMP: Second derivative of (1/||x||^2 - 1/Delta^2)^2 -! Version with OMP - -! *Compute the first derivative of (1/||x||^2 - 1/Delta^2)^2* - -! This function computes the value of (1/||x(lambda)||^2 - 1/Delta^2)^2 - -! \begin{align*} -! \frac{\partial^2}{\partial \lambda^2} (1/||\textbf{x}(\lambda)||^2 - 1/\Delta^2)^2 -! &= 4 \left[ \frac{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3)})^2}{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^4} -! - 3 \frac{\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^4}}{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^3} \right] \\ -! &- \frac{4}{\Delta^2} \left[ \frac{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3)})^2}{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^3} -! - 3 \frac{\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^4}}{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^2} \right] -! \end{align*} - - -! \begin{align*} -! \text{accu1} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^2} \\ -! \text{accu2} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^3} \\ -! \text{accu3} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^4} -! \end{align*} - -! Provided: -! | m_num | integer | number of MOs | - -! Input: -! | n | integer | mo_num*(mo_num-1)/2 | -! | e_val(n) | double precision | eigenvalues of the hessian | -! | W(n,n) | double precision | eigenvectors of the hessian | -! | v_grad(n) | double precision | gradient | -! | lambda | double precision | Lagrange multiplier | -! | delta | double precision | Delta of the trust region | - -! Internal: -! | wtg | double precision | temporary variable to store W^T.v_grad | -! | tmp_accu1 | double precision | temporary array for the first sum | -! | tmp_accu2 | double precision | temporary array for the second sum | -! | tmp_wtg(n) | double precision | temporary array for W^t.v_grad | -! | i,j | integer | indexes | - -! Function: -! | d1_norm_inverse_trust_region | double precision | value of the first derivative | - - -function d2_norm_inverse_trust_region_omp(n,e_val,tmp_wtg,lambda,delta) - - use omp_lib - include 'pi.h' - - BEGIN_DOC - ! Compute the second derivative of (1/||x||^2 - 1/Delta^2)^2 - END_DOC - - implicit none - - ! Variables - - ! in - integer, intent(in) :: n - double precision, intent(in) :: e_val(n) - double precision, intent(in) :: tmp_wtg(n) - double precision, intent(in) :: lambda - double precision, intent(in) :: delta - - ! Internal - double precision :: accu1, accu2, accu3 - integer :: i,j - double precision, allocatable :: tmp_accu1(:), tmp_accu2(:), tmp_accu3(:) - - ! Functions - double precision :: d2_norm_inverse_trust_region_omp - - ! Allocation - allocate(tmp_accu1(n), tmp_accu2(n), tmp_accu3(n)) - - ! OMP - call omp_set_max_active_levels(1) - - ! OMP - !$OMP PARALLEL & - !$OMP PRIVATE(i,j) & - !$OMP SHARED(n,lambda, e_val, thresh_eig,& - !$OMP tmp_accu1, tmp_accu2, tmp_accu3, tmp_wtg, & - !$OMP accu1, accu2, accu3) & - !$OMP DEFAULT(NONE) - - !$OMP MASTER - accu1 = 0d0 - accu2 = 0d0 - accu3 = 0d0 - !$OMP END MASTER - - !$OMP DO - do i = 1, n - tmp_accu1(i) = 0d0 - enddo - !$OMP END DO - - !$OMP DO - do i = 1, n - tmp_accu2(i) = 0d0 - enddo - !$OMP END DO - - !$OMP DO - do i = 1, n - tmp_accu3(i) = 0d0 - enddo - !$OMP END DO - - !$OMP DO - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - tmp_accu1(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**2 - endif - enddo - !$OMP END DO - - !$OMP MASTER - do i = 1, n - accu1 = accu1 + tmp_accu1(i) - enddo - !$OMP END MASTER - - !$OMP DO - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - tmp_accu2(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**3 - endif - enddo - !$OMP END DO - - !$OMP MASTER - do i = 1, n - accu2 = accu2 + tmp_accu2(i) - enddo - !$OMP END MASTER - - !$OMP DO - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - tmp_accu3(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**4 - endif - enddo - !$OMP END DO - - !$OMP MASTER - do i = 1, n - accu3 = accu3 + tmp_accu3(i) - enddo - !$OMP END MASTER - - !$OMP END PARALLEL - - call omp_set_max_active_levels(4) - - d2_norm_inverse_trust_region_omp = 4d0 * (6d0 * accu2**2/accu1**4 - 3d0 * accu3/accu1**3) & - - 4d0/delta**2 * (4d0 * accu2**2/accu1**3 - 3d0 * accu3/accu1**2) - - deallocate(tmp_accu1,tmp_accu2,tmp_accu3) - -end - -! First derivative of (1/||x||^2 - 1/Delta^2)^2 -! Version without OMP - -! *Compute the first derivative of (1/||x||^2 - 1/Delta^2)^2* - -! This function computes the value of (1/||x(lambda)||^2 - 1/Delta^2)^2 - -! \begin{align*} -! \frac{\partial}{\partial \lambda} (1/||\textbf{x}(\lambda)||^2 - 1/\Delta^2)^2 -! &= 4 \frac{\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3}} -! {(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^3} -! - \frac{4}{\Delta^2} \frac{\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3)}} -! {(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^2} \\ -! &= 4 \sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3} -! \left( \frac{1}{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^3} -! - \frac{1}{\Delta^2 (\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^2} \right) -! \end{align*} -! \begin{align*} -! \text{accu1} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^2} \\ -! \text{accu2} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^3} -! \end{align*} -! Provided: -! | m_num | integer | number of MOs | - -! Input: -! | n | integer | mo_num*(mo_num-1)/2 | -! | e_val(n) | double precision | eigenvalues of the hessian | -! | W(n,n) | double precision | eigenvectors of the hessian | -! | v_grad(n) | double precision | gradient | -! | lambda | double precision | Lagrange multiplier | -! | delta | double precision | Delta of the trust region | - -! Internal: -! | wtg | double precision | temporary variable to store W^T.v_grad | -! | i,j | integer | indexes | - -! Function: -! | d1_norm_inverse_trust_region | double precision | value of the first derivative | - - -function d1_norm_inverse_trust_region(n,e_val,w,v_grad,lambda,delta) - - include 'pi.h' - - BEGIN_DOC - ! Compute the first derivative of (1/||x||^2 - 1/Delta^2)^2 - END_DOC - - implicit none - - ! Variables - - ! in - integer, intent(in) :: n - double precision, intent(in) :: e_val(n) - double precision, intent(in) :: w(n,n) - double precision, intent(in) :: v_grad(n) - double precision, intent(in) :: lambda - double precision, intent(in) :: delta - - ! Internal - double precision :: wtg, accu1, accu2 - integer :: i,j - - ! Functions - double precision :: d1_norm_inverse_trust_region - - accu1 = 0d0 - accu2 = 0d0 - - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - wtg = 0d0 - do j = 1, n - wtg = wtg + w(j,i) * v_grad(j) - enddo - accu1 = accu1 + wtg**2 / (e_val(i) + lambda)**2 - endif - enddo - - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - wtg = 0d0 - do j = 1, n - wtg = wtg + w(j,i) * v_grad(j) - enddo - accu2 = accu2 + wtg**2 / (e_val(i) + lambda)**3 - endif - enddo - - d1_norm_inverse_trust_region = 4d0 * accu2 * (1d0/accu1**3 - 1d0/(delta**2 * accu1**2)) - -end - -! Second derivative of (1/||x||^2 - 1/Delta^2)^2 -! Version without OMP - -! *Compute the second derivative of (1/||x||^2 - 1/Delta^2)^2* - -! This function computes the value of (1/||x(lambda)||^2 - 1/Delta^2)^2 - -! \begin{align*} -! \frac{\partial^2}{\partial \lambda^2} (1/||\textbf{x}(\lambda)||^2 - 1/\Delta^2)^2 -! &= 4 \left[ \frac{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3)})^2}{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^4} -! - 3 \frac{\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^4}}{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^3} \right] \\ -! &- \frac{4}{\Delta^2} \left[ \frac{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3)})^2}{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^3} -! - 3 \frac{\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^4}}{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^2} \right] -! \end{align*} - -! \begin{align*} -! \text{accu1} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^2} \\ -! \text{accu2} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^3} \\ -! \text{accu3} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^4} -! \end{align*} - -! Provided: -! | m_num | integer | number of MOs | - -! Input: -! | n | integer | mo_num*(mo_num-1)/2 | -! | e_val(n) | double precision | eigenvalues of the hessian | -! | W(n,n) | double precision | eigenvectors of the hessian | -! | v_grad(n) | double precision | gradient | -! | lambda | double precision | Lagrange multiplier | -! | delta | double precision | Delta of the trust region | - -! Internal: -! | wtg | double precision | temporary variable to store W^T.v_grad | -! | i,j | integer | indexes | - -! Function: -! | d2_norm_inverse_trust_region | double precision | value of the first derivative | - - -function d2_norm_inverse_trust_region(n,e_val,w,v_grad,lambda,delta) - - include 'pi.h' - - BEGIN_DOC - ! Compute the second derivative of (1/||x||^2 - 1/Delta^2)^2 - END_DOC - - implicit none - - ! Variables - - ! in - integer, intent(in) :: n - double precision, intent(in) :: e_val(n) - double precision, intent(in) :: w(n,n) - double precision, intent(in) :: v_grad(n) - double precision, intent(in) :: lambda - double precision, intent(in) :: delta - - ! Internal - double precision :: wtg, accu1, accu2, accu3 - integer :: i,j - - ! Functions - double precision :: d2_norm_inverse_trust_region - - accu1 = 0d0 - accu2 = 0d0 - accu3 = 0d0 - - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - wtg = 0d0 - do j = 1, n - wtg = wtg + w(j,i) * v_grad(j) - enddo - accu1 = accu1 + wtg**2 / (e_val(i) + lambda)**2 - endif - enddo - - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - wtg = 0d0 - do j = 1, n - wtg = wtg + w(j,i) * v_grad(j) - enddo - accu2 = accu2 + wtg**2 / (e_val(i) + lambda)**3 - endif - enddo - - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - wtg = 0d0 - do j = 1, n - wtg = wtg + w(j,i) * v_grad(j) - enddo - accu3 = accu3 + wtg**2 / (e_val(i) + lambda)**4 - endif - enddo - - d2_norm_inverse_trust_region = 4d0 * (6d0 * accu2**2/accu1**4 - 3d0 * accu3/accu1**3) & - - 4d0/delta**2 * (4d0 * accu2**2/accu1**3 - 3d0 * accu3/accu1**2) - -end diff --git a/src/utils_trust_region/trust_region_optimal_lambda.org b/src/utils_trust_region/trust_region_optimal_lambda.org deleted file mode 100644 index b39c9a10..00000000 --- a/src/utils_trust_region/trust_region_optimal_lambda.org +++ /dev/null @@ -1,1665 +0,0 @@ -* Newton's method to find the optimal lambda - -*Compute the lambda value for the trust region* - -This subroutine uses the Newton method in order to find the optimal -lambda. This constant is added on the diagonal of the hessian to shift -the eiganvalues. It has a double role: -- ensure that the resulting hessian is positive definite for the - Newton method -- constrain the step in the trust region, i.e., - $||\textbf{x}(\lambda)|| \leq \Delta$, where $\Delta$ is the radius - of the trust region. -We search $\lambda$ which minimizes -\begin{align*} - f(\lambda) = (||\textbf{x}_{(k+1)}(\lambda)||^2 -\Delta^2)^2 -\end{align*} -or -\begin{align*} - \tilde{f}(\lambda) = (\frac{1}{||\textbf{x}_{(k+1)}(\lambda)||^2}-\frac{1}{\Delta^2})^2 -\end{align*} -and gives obviously 0 in both cases. \newline - -There are several cases: -- If $\textbf{H}$ is positive definite the interval containing the - solution is $\lambda \in (0, \infty)$ (and $-h_1 < 0$). -- If $\textbf{H}$ is indefinite ($h_1 < 0$) and $\textbf{w}_1^T \cdot - \textbf{g} \neq 0$ then the interval containing - the solution is $\lambda \in (-h_1, \infty)$. -- If $\textbf{H}$ is indefinite ($h_1 < 0$) and $\textbf{w}_1^T \cdot - \textbf{g} = 0$ then the interval containing the solution is - $\lambda \in (-h_1, \infty)$. The terms where $|h_i - \lambda| < - 10^{-12}$ are not computed, so the term where $i = 1$ is - automatically removed and this case becomes similar to the previous one. - -So to avoid numerical problems (cf. trust_region) we start the -algorithm at $\lambda=\max(0 + \epsilon,-h_1 + \epsilon)$, -with $\epsilon$ a little constant. -The research must be restricted to the interval containing the -solution. For that reason a little trust region in 1D is used. - -The Newton method to find the optimal $\lambda$ is : -\begin{align*} - \lambda_{(l+1)} &= \lambda_{(l)} - f^{''}(\lambda)_{(l)}^{-1} f^{'}(\lambda)_{(l)}^{} \\ -\end{align*} -$f^{'}(\lambda)_{(l)}$: the first derivative of $f$ with respect to -$\lambda$ at the l-th iteration, -$f^{''}(\lambda)_{(l)}$: the second derivative of $f$ with respect to -$\lambda$ at the l-th iteration.\newline - -Noting the Newton step $y = - f^{''}(\lambda)_{(l)}^{-1} -f^{'}(\lambda)_{(l)}^{}$ we constrain $y$ such as -\begin{align*} - y \leq \alpha -\end{align*} -with $\alpha$ a scalar representing the trust length (trust region in -1D) where the function $f$ or $\tilde{f}$ is correctly describe by the -Taylor series truncated at the second order. Thus, if $y > \alpha$, -the constraint is applied as -\begin{align*} - y^* = \alpha \frac{y}{|y|} -\end{align*} -with $y^*$ the solution in the trust region. - -The size of the trust region evolves in function of $\rho$ as for the -trust region seen previously cf. trust_region, rho_model. -The prediction of the value of $f$ or $\tilde{f}$ is done using the -Taylor series truncated at the second order cf. "trust_region", -"trust_e_model". - -The first and second derivatives of $f(\lambda) = (||\textbf{x}(\lambda)||^2 - -\Delta^2)^2$ with respect to $\lambda$ are: -\begin{align*} - \frac{\partial }{\partial \lambda} (||\textbf{x}(\lambda)||^2 - \Delta^2)^2 - = 2 \left(\sum_{i=1}^n \frac{-2(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^3} \right) - \left( - \Delta^2 + \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i+ \lambda)^2} \right) -\end{align*} -\begin{align*} -\frac{\partial^2}{\partial \lambda^2} (||\textbf{x}(\lambda)||^2 - \Delta^2)^2 -= 2 \left[ \left( \sum_{i=1}^n 6 \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^4} \right) \left( - \Delta^2 + \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^2} \right) + \left( \sum_{i=1}^n -2 \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^3} \right)^2 \right] -\end{align*} - -The first and second derivatives of $\tilde{f}(\lambda) = (1/||\textbf{x}(\lambda)||^2 - -1/\Delta^2)^2$ with respect to $\lambda$ are: -\begin{align*} - \frac{\partial}{\partial \lambda} (1/||\textbf{x}(\lambda)||^2 - 1/\Delta^2)^2 - &= 4 \frac{\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3}} - {(\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^3} - - \frac{4}{\Delta^2} \frac{\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3)}} - {(\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^2} \\ - &= 4 \sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3} - \left( \frac{1}{(\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^3} - - \frac{1}{\Delta^2 (\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^2} \right) -\end{align*} - -\begin{align*} - \frac{\partial^2}{\partial \lambda^2} (1/||\textbf{x}(\lambda)||^2 - 1/\Delta^2)^2 - &= 4 \left[ \frac{(\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3)})^2} - {(\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^4} - - 3 \frac{\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^4}} - {(\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^3} \right] \\ - &- \frac{4}{\Delta^2} \left[ \frac{(\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2} - {(h_i + \lambda)^3)})^2}{(\sum_ {i=1}^n\frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^3} - - 3 \frac{\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^4}} - {(\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^2} \right] -\end{align*} - -Provided in qp_edit: -| thresh_rho_2 | -| thresh_cc | -| nb_it_max_lambda | -| version_lambda_search | -| nb_it_max_pre_search | -see qp_edit for more details - -Input: -| n | integer | m*(m-1)/2 | -| e_val(n) | double precision | eigenvalues of the hessian | -| tmp_wtg(n) | double precision | w_i^T.v_grad(i) | -| delta | double precision | delta for the trust region | - -Output: -| lambda | double precision | Lagrange multiplier to constrain the norm of the size of the Newton step | -| | | lambda > 0 | - -Internal: -| d1_N | double precision | value of d1_norm_trust_region | -| d2_N | double precision | value of d2_norm_trust_region | -| f_N | double precision | value of f_norm_trust_region | -| prev_f_N | double precision | previous value of f_norm_trust_region | -| f_R | double precision | (norm(x)^2 - delta^2)^2 or (1/norm(x)^2 - 1/delta^2)^2 | -| prev_f_R | double precision | previous value of f_R | -| model | double precision | predicted value of f_R from prev_f_R and y | -| d_1 | double precision | value of the first derivative | -| d_2 | double precision | value of the second derivative | -| y | double precision | Newton's step, y = -f''^-1 . f' = lambda - prev_lambda | -| prev_lambda | double precision | previous value of lambda | -| t1,t2,t3 | double precision | wall time | -| i | integer | index | -| epsilon | double precision | little constant to avoid numerical problem | -| rho_2 | double precision | (prev_f_R - f_R)/(prev_f_R - model), agreement between model and f_R | -| version | integer | version of the root finding method | - -Function: -| d1_norm_trust_region | double precision | first derivative with respect to lambda of (norm(x)^2 - Delta^2)^2 | -| d2_norm_trust_region | double precision | first derivative with respect to lambda of (norm(x)^2 - Delta^2)^2 | -| d1_norm_inverse_trust_region | double precision | first derivative with respect to lambda of (1/norm(x)^2 - 1/Delta^2)^2 | -| d2_norm_inverse_trust_region | double precision | second derivative with respect to lambda of (1/norm(x)^2 - 1/Delta^2)^2 | -| f_norm_trust_region | double precision | value of norm(x)^2 | - - -#+BEGIN_SRC f90 :comments org :tangle trust_region_optimal_lambda.irp.f -subroutine trust_region_optimal_lambda(n,e_val,tmp_wtg,delta,lambda) - - include 'pi.h' - - BEGIN_DOC - ! Research the optimal lambda to constrain the step size in the trust region - END_DOC - - implicit none - - ! Variables - - ! in - integer, intent(in) :: n - double precision, intent(inout) :: e_val(n) - double precision, intent(in) :: delta - double precision, intent(in) :: tmp_wtg(n) - - ! out - double precision, intent(out) :: lambda - - ! Internal - double precision :: d1_N, d2_N, f_N, prev_f_N - double precision :: prev_f_R, f_R - double precision :: model - double precision :: d_1, d_2 - double precision :: t1,t2,t3 - integer :: i - double precision :: epsilon - double precision :: y - double precision :: prev_lambda - double precision :: rho_2 - double precision :: alpha - integer :: version - - ! Functions - double precision :: d1_norm_trust_region,d1_norm_trust_region_omp - double precision :: d2_norm_trust_region, d2_norm_trust_region_omp - double precision :: f_norm_trust_region, f_norm_trust_region_omp - double precision :: d1_norm_inverse_trust_region - double precision :: d2_norm_inverse_trust_region - double precision :: d1_norm_inverse_trust_region_omp - double precision :: d2_norm_inverse_trust_region_omp - - print*,'' - print*,'---Trust_newton---' - print*,'' - - call wall_time(t1) - - ! version_lambda_search - ! 1 -> ||x||^2 - delta^2 = 0, - ! 2 -> 1/||x||^2 - 1/delta^2 = 0 (better) - if (version_lambda_search == 1) then - print*, 'Research of the optimal lambda by solving ||x||^2 - delta^2 = 0' - else - print*, 'Research of the optimal lambda by solving 1/||x||^2 - 1/delta^2 = 0' - endif - ! Version 2 is normally better -#+END_SRC - -Resolution with the Newton method: - -#+BEGIN_SRC f90 :comments org :tangle trust_region_optimal_lambda.irp.f - ! Initialization - epsilon = 1d-4 - lambda =MAX(0d0, -e_val(1)) - - ! Pre research of lambda to start near the optimal lambda - ! by adding a constant epsilon and changing the constant to - ! have ||x(lambda + epsilon)|| ~ delta, before setting - ! lambda = lambda + epsilon - print*, 'Pre research of lambda:' - print*,'Initial lambda =', lambda - f_N = f_norm_trust_region_omp(n,e_val,tmp_wtg,lambda + epsilon) - print*,'||x(lambda)||=', dsqrt(f_N),'delta=',delta - i = 1 - - ! To increase lambda - if (f_N > delta**2) then - print*,'Increasing lambda...' - do while (f_N > delta**2 .and. i <= nb_it_max_pre_search) - - ! Update the previous norm - prev_f_N = f_N - ! New epsilon - epsilon = epsilon * 2d0 - ! New norm - f_N = f_norm_trust_region_omp(n,e_val,tmp_wtg,lambda + epsilon) - - print*, 'lambda', lambda + epsilon, '||x||', dsqrt(f_N), 'delta', delta - - ! Security - if (prev_f_N < f_N) then - print*,'WARNING, error: prev_f_N < f_N, exit' - epsilon = epsilon * 0.5d0 - i = nb_it_max_pre_search + 1 - endif - - i = i + 1 - enddo - - ! To reduce lambda - else - print*,'Reducing lambda...' - do while (f_N < delta**2 .and. i <= nb_it_max_pre_search) - - ! Update the previous norm - prev_f_N = f_N - ! New epsilon - epsilon = epsilon * 0.5d0 - ! New norm - f_N = f_norm_trust_region_omp(n,e_val,tmp_wtg,lambda + epsilon) - - print*, 'lambda', lambda + epsilon, '||x||', dsqrt(f_N), 'delta', delta - - ! Security - if (prev_f_N > f_N) then - print*,'WARNING, error: prev_f_N > f_N, exit' - epsilon = epsilon * 2d0 - i = nb_it_max_pre_search + 1 - endif - - i = i + 1 - enddo - endif - - print*,'End of the pre research of lambda' - - ! New value of lambda - lambda = lambda + epsilon - - print*, 'e_val(1):', e_val(1) - print*, 'Staring point, lambda =', lambda - - ! thresh_cc, threshold for the research of the optimal lambda - ! Leaves the loop when ABS(1d0-||x||^2/delta^2) > thresh_cc - ! thresh_rho_2, threshold to cancel the step in the research - ! of the optimal lambda, the step is cancelled if rho_2 < thresh_rho_2 - print*,'Threshold for the CC:', thresh_cc - print*,'Threshold for rho_2:', thresh_rho_2 - - print*, 'w_1^T . g =', tmp_wtg(1) - - ! Debug - !if (debug) then - ! print*, 'Iteration rho_2 lambda delta ||x|| |1-(||x||^2/delta^2)|' - !endif - - ! Initialization - i = 1 - f_N = f_norm_trust_region_omp(n,e_val,tmp_wtg,lambda) ! Value of the ||x(lambda)||^2 - model = 0d0 ! predicted value of (||x||^2 - delta^2)^2 - prev_f_N = 0d0 ! previous value of ||x||^2 - prev_f_R = 0d0 ! previous value of (||x||^2 - delta^2)^2 - f_R = 0d0 ! value of (||x||^2 - delta^2)^2 - rho_2 = 0d0 ! (prev_f_R - f_R)/(prev_f_R - m) - y = 0d0 ! step size - prev_lambda = 0d0 ! previous lambda - - ! Derivatives - if (version_lambda_search == 1) then - d_1 = d1_norm_trust_region_omp(n,e_val,tmp_wtg,lambda,delta) ! first derivative of (||x(lambda)||^2 - delta^2)^2 - d_2 = d2_norm_trust_region_omp(n,e_val,tmp_wtg,lambda,delta) ! second derivative of (||x(lambda)||^2 - delta^2)^2 - else - d_1 = d1_norm_inverse_trust_region_omp(n,e_val,tmp_wtg,lambda,delta) ! first derivative of (1/||x(lambda)||^2 - 1/delta^2)^2 - d_2 = d2_norm_inverse_trust_region_omp(n,e_val,tmp_wtg,lambda,delta) ! second derivative of (1/||x(lambda)||^2 - 1/delta^2)^2 - endif - - ! Trust length - alpha = DABS((1d0/d_2)*d_1) - - ! Newton's method - do while (i <= 100 .and. DABS(1d0-f_N/delta**2) > thresh_cc) - print*,'--------------------------------------' - print*,'Research of lambda, iteration:', i - print*,'--------------------------------------' - - ! Update of f_N, f_R and the derivatives - prev_f_N = f_N - if (version_lambda_search == 1) then - prev_f_R = (prev_f_N - delta**2)**2 - d_1 = d1_norm_trust_region_omp(n,e_val,tmp_wtg,lambda,delta) ! first derivative of (||x(lambda)||^2 - delta^2)^2 - d_2 = d2_norm_trust_region_omp(n,e_val,tmp_wtg,lambda,delta) ! second derivative of (||x(lambda)||^2 - delta^2)^2 - else - prev_f_R = (1d0/prev_f_N - 1d0/delta**2)**2 - d_1 = d1_norm_inverse_trust_region_omp(n,e_val,tmp_wtg,lambda,delta) ! first derivative of (1/||x(lambda)||^2 - 1/delta^2)^2 - d_2 = d2_norm_inverse_trust_region_omp(n,e_val,tmp_wtg,lambda,delta) ! second derivative of (1/||x(lambda)||^2 - 1/delta^2)^2 - endif - write(*,'(a,E12.5,a,E12.5)') ' 1st and 2nd derivative: ', d_1,', ', d_2 - - ! Newton's step - y = -(1d0/DABS(d_2))*d_1 - - ! Constraint on y (the newton step) - if (DABS(y) > alpha) then - y = alpha * (y/DABS(y)) ! preservation of the sign of y - endif - write(*,'(a,E12.5)') ' Step length: ', y - - ! Predicted value of (||x(lambda)||^2 - delta^2)^2, Taylor series - model = prev_f_R + d_1 * y + 0.5d0 * d_2 * y**2 - - ! Updates lambda - prev_lambda = lambda - lambda = prev_lambda + y - print*,'prev lambda:', prev_lambda - print*,'new lambda:', lambda - - ! Checks if lambda is in (-h_1, \infty) - if (lambda > MAX(0d0, -e_val(1))) then - ! New value of ||x(lambda)||^2 - f_N = f_norm_trust_region_omp(n,e_val,tmp_wtg,lambda) - - ! New f_R - if (version_lambda_search == 1) then - f_R = (f_N - delta**2)**2 ! new value of (||x(lambda)||^2 - delta^2)^2 - else - f_R = (1d0/f_N - 1d0/delta**2)**2 ! new value of (1/||x(lambda)||^2 -1/delta^2)^2 - endif - - if (version_lambda_search == 1) then - print*,'Previous value of (||x(lambda)||^2 - delta^2)^2:', prev_f_R - print*,'Actual value of (||x(lambda)||^2 - delta^2)^2:', f_R - print*,'Predicted value of (||x(lambda)||^2 - delta^2)^2:', model - else - print*,'Previous value of (1/||x(lambda)||^2 - 1/delta^2)^2:', prev_f_R - print*,'Actual value of (1/||x(lambda)||^2 - 1/delta^2)^2:', f_R - print*,'Predicted value of (1/||x(lambda)||^2 - 1/delta^2)^2:', model - endif - - print*,'previous - actual:', prev_f_R - f_R - print*,'previous - model:', prev_f_R - model - - ! Check the gain - if (DABS(prev_f_R - model) < thresh_model_2) then - print*,'' - print*,'WARNING: ABS(previous - model) <', thresh_model_2, 'rho_2 will tend toward infinity' - print*,'' - endif - - ! Will be deleted - !if (prev_f_R - f_R <= 1d-16 .or. prev_f_R - model <= 1d-16) then - ! print*,'' - ! print*,'WARNING: ABS(previous - model) <= 1d-16, exit' - ! print*,'' - ! exit - !endif - - ! Computes rho_2 - rho_2 = (prev_f_R - f_R)/(prev_f_R - model) - print*,'rho_2:', rho_2 - else - rho_2 = 0d0 ! in order to reduce the size of the trust region, alpha, until lambda is in (-h_1, \infty) - print*,'lambda < -e_val(1) ===> rho_2 = 0' - endif - - ! Evolution of the trust length, alpha - if (rho_2 >= 0.75d0) then - alpha = 2d0 * alpha - elseif (rho_2 >= 0.5d0) then - alpha = alpha - elseif (rho_2 >= 0.25d0) then - alpha = 0.5d0 * alpha - else - alpha = 0.25d0 * alpha - endif - write(*,'(a,E12.5)') ' New trust length alpha: ', alpha - - ! cancellaion of the step if rho < 0.1 - if (rho_2 < thresh_rho_2) then !0.1d0) then - lambda = prev_lambda - f_N = prev_f_N - print*,'Rho_2 <', thresh_rho_2,', cancellation of the step: lambda = prev_lambda' - endif - - print*,'' - print*,'lambda, ||x||, delta:' - print*, lambda, dsqrt(f_N), delta - print*,'CC:', DABS(1d0 - f_N/delta**2) - print*,'' - - i = i + 1 - enddo - - ! if trust newton failed - if (i > nb_it_max_lambda) then - print*,'' - print*,'######################################################' - print*,'WARNING: i >', nb_it_max_lambda,'for the trust Newton' - print*,'The research of the optimal lambda has failed' - print*,'######################################################' - print*,'' - endif - - print*,'Number of iterations :', i - print*,'Value of lambda :', lambda - print*,'Error on the trust region (1d0-f_N/delta**2) (Convergence criterion) :', 1d0-f_N/delta**2 - print*,'Error on the trust region (||x||^2 - delta^2)^2) :', (f_N - delta**2)**2 - print*,'Error on the trust region (1/||x||^2 - 1/delta^2)^2)', (1d0/f_N - 1d0/delta**2)**2 - - ! Time - call wall_time(t2) - t3 = t2 - t1 - print*,'Time in trust_newton:', t3 - - print*,'' - print*,'---End trust_newton---' - print*,'' - -end subroutine -#+END_SRC - -* OMP: First derivative of (||x||^2 - Delta^2)^2 - -*Function to compute the first derivative of (||x||^2 - Delta^2)^2* - -This function computes the first derivative of (||x||^2 - Delta^2)^2 -with respect to lambda. - -\begin{align*} -\frac{\partial }{\partial \lambda} (||\textbf{x}(\lambda)||^2 - \Delta^2)^2 -= -4 \left(\sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3} \right) -\left( - \Delta^2 + \sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i+ \lambda)^2} \right) -\end{align*} - -\begin{align*} - \text{accu1} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2} \\ - \text{accu2} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3} -\end{align*} - -Provided: -| mo_num | integer | number of MOs | - -Input: -| n | integer | mo_num*(mo_num-1)/2 | -| e_val(n) | double precision | eigenvalues of the hessian | -| W(n,n) | double precision | eigenvectors of the hessian | -| v_grad(n) | double precision | gradient | -| lambda | double precision | Lagrange multiplier | -| delta | double precision | Delta of the trust region | - -Internal: -| accu1 | double precision | first sum of the formula | -| accu2 | double precision | second sum of the formula | -| tmp_accu1 | double precision | temporary array for the first sum | -| tmp_accu2 | double precision | temporary array for the second sum | -| tmp_wtg(n) | double precision | temporary array for W^t.v_grad | -| i,j | integer | indexes | - -Function: -| d1_norm_trust_region | double precision | first derivative with respect to lambda of (norm(x)^2 - Delta^2)^2 | - -#+BEGIN_SRC f90 :comments org :tangle trust_region_optimal_lambda.irp.f -function d1_norm_trust_region_omp(n,e_val,tmp_wtg,lambda,delta) - - use omp_lib - include 'pi.h' - - BEGIN_DOC - ! Compute the first derivative with respect to lambda of (||x(lambda)||^2 - Delta^2)^2 - END_DOC - - implicit none - - ! in - integer, intent(in) :: n - double precision, intent(in) :: e_val(n) - double precision, intent(in) :: tmp_wtg(n) - double precision, intent(in) :: lambda - double precision, intent(in) :: delta - - ! Internal - double precision :: wtg,accu1,accu2 - integer :: i,j - double precision, allocatable :: tmp_accu1(:), tmp_accu2(:) - - ! Functions - double precision :: d1_norm_trust_region_omp - - ! Allocation - allocate(tmp_accu1(n), tmp_accu2(n)) - - ! OMP - call omp_set_max_active_levels(1) - - ! OMP - !$OMP PARALLEL & - !$OMP PRIVATE(i,j) & - !$OMP SHARED(n,lambda, e_val, thresh_eig,& - !$OMP tmp_accu1, tmp_accu2, tmp_wtg, accu1,accu2) & - !$OMP DEFAULT(NONE) - - !$OMP MASTER - accu1 = 0d0 - accu2 = 0d0 - !$OMP END MASTER - - !$OMP DO - do i = 1, n - tmp_accu1(i) = 0d0 - enddo - !$OMP END DO - - !$OMP DO - do i = 1, n - tmp_accu2(i) = 0d0 - enddo - !$OMP END DO - - !$OMP DO - do i = 1, n - if (ABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - tmp_accu1(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**2 - endif - enddo - !$OMP END DO - - !$OMP MASTER - do i = 1, n - accu1 = accu1 + tmp_accu1(i) - enddo - !$OMP END MASTER - - !$OMP DO - do i = 1, n - if (ABS(e_val(i)) > thresh_eig) then - tmp_accu2(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**3 - endif - enddo - !$OMP END DO - - !$OMP MASTER - do i = 1, n - accu2 = accu2 + tmp_accu2(i) - enddo - !$OMP END MASTER - - !$OMP END PARALLEL - - call omp_set_max_active_levels(4) - - d1_norm_trust_region_omp = -4d0 * accu2 * (accu1 - delta**2) - - deallocate(tmp_accu1, tmp_accu2) - -end function -#+END_SRC - -* OMP: Second derivative of (||x||^2 - Delta^2)^2 - -*Function to compute the second derivative of (||x||^2 - Delta^2)^2* - -This function computes the second derivative of (||x||^2 - Delta^2)^2 -with respect to lambda. -\begin{align*} -\frac{\partial^2 }{\partial \lambda^2} (||\textbf{x}(\lambda)||^2 - \Delta^2)^2 -= 2 \left[ \left( \sum_{i=1}^n 6 \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^4} \right) \left( - \Delta^2 + \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^2} \right) + \left( \sum_{i=1}^n -2 \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^3} \right)^2 \right] -\end{align*} - -\begin{align*} - \text{accu1} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^2} \\ - \text{accu2} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^3} \\ - \text{accu3} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^4} -\end{align*} - -Provided: -| m_num | integer | number of MOs | - -Input: -| n | integer | mo_num*(mo_num-1)/2 | -| e_val(n) | double precision | eigenvalues of the hessian | -| W(n,n) | double precision | eigenvectors of the hessian | -| v_grad(n) | double precision | gradient | -| lambda | double precision | Lagrange multiplier | -| delta | double precision | Delta of the trust region | - -Internal: -| accu1 | double precision | first sum of the formula | -| accu2 | double precision | second sum of the formula | -| accu3 | double precision | third sum of the formula | -| tmp_accu1 | double precision | temporary array for the first sum | -| tmp_accu2 | double precision | temporary array for the second sum | -| tmp_accu2 | double precision | temporary array for the third sum | -| tmp_wtg(n) | double precision | temporary array for W^t.v_grad | -| i,j | integer | indexes | - -Function: -| d2_norm_trust_region | double precision | second derivative with respect to lambda of (norm(x)^2 - Delta^2)^2 | - -#+BEGIN_SRC f90 :comments org :tangle trust_region_optimal_lambda.irp.f -function d2_norm_trust_region_omp(n,e_val,tmp_wtg,lambda,delta) - - use omp_lib - include 'pi.h' - - BEGIN_DOC - ! Compute the second derivative with respect to lambda of (||x(lambda)||^2 - Delta^2)^2 - END_DOC - - implicit none - - ! Variables - - ! in - integer, intent(in) :: n - double precision, intent(in) :: e_val(n) - double precision, intent(in) :: tmp_wtg(n) - double precision, intent(in) :: lambda - double precision, intent(in) :: delta - - ! Functions - double precision :: d2_norm_trust_region_omp - double precision :: ddot - - ! Internal - double precision :: accu1,accu2,accu3 - double precision, allocatable :: tmp_accu1(:), tmp_accu2(:), tmp_accu3(:) - integer :: i, j - - ! Allocation - allocate(tmp_accu1(n), tmp_accu2(n), tmp_accu3(n)) - - call omp_set_max_active_levels(1) - - ! OMP - !$OMP PARALLEL & - !$OMP PRIVATE(i,j) & - !$OMP SHARED(n,lambda, e_val, thresh_eig,& - !$OMP tmp_accu1, tmp_accu2, tmp_accu3, tmp_wtg, & - !$OMP accu1, accu2, accu3) & - !$OMP DEFAULT(NONE) - - ! Initialization - - !$OMP MASTER - accu1 = 0d0 - accu2 = 0d0 - accu3 = 0d0 - !$OMP END MASTER - - !$OMP DO - do i = 1, n - tmp_accu1(i) = 0d0 - enddo - !$OMP END DO - !$OMP DO - do i = 1, n - tmp_accu2(i) = 0d0 - enddo - !$OMP END DO - !$OMP DO - do i = 1, n - tmp_accu3(i) = 0d0 - enddo - !$OMP END DO - - ! Calculations - - ! accu1 - !$OMP DO - do i = 1, n - if (ABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - tmp_accu1(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**2 - endif - enddo - !$OMP END DO - - !$OMP MASTER - do i = 1, n - accu1 = accu1 + tmp_accu1(i) - enddo - !$OMP END MASTER - - ! accu2 - !$OMP DO - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - tmp_accu2(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**3 - endif - enddo - !$OMP END DO - - ! accu3 - !$OMP MASTER - do i = 1, n - accu2 = accu2 + tmp_accu2(i) - enddo - !$OMP END MASTER - - !$OMP DO - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - tmp_accu3(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**4 - endif - enddo - !$OMP END DO - - !$OMP MASTER - do i = 1, n - accu3 = accu3 + tmp_accu3(i) - enddo - !$OMP END MASTER - - !$OMP END PARALLEL - - d2_norm_trust_region_omp = 2d0 * (6d0 * accu3 * (- delta**2 + accu1) + (-2d0 * accu2)**2) - - deallocate(tmp_accu1, tmp_accu2, tmp_accu3) - -end function -#+END_SRC - -* OMP: Function value of ||x||^2 - -*Compute the value of ||x||^2* - -This function computes the value of ||x(lambda)||^2 - -\begin{align*} -||\textbf{x}(\lambda)||^2 = \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^2} -\end{align*} - -Provided: -| m_num | integer | number of MOs | - -Input: -| n | integer | mo_num*(mo_num-1)/2 | -| e_val(n) | double precision | eigenvalues of the hessian | -| W(n,n) | double precision | eigenvectors of the hessian | -| v_grad(n) | double precision | gradient | -| lambda | double precision | Lagrange multiplier | - -Internal: -| tmp_wtg(n) | double precision | temporary array for W^T.v_grad | -| tmp_fN | double precision | temporary array for the function | -| i,j | integer | indexes | - -#+BEGIN_SRC f90 :comments org :tangle trust_region_optimal_lambda.irp.f -function f_norm_trust_region_omp(n,e_val,tmp_wtg,lambda) - - use omp_lib - - include 'pi.h' - - BEGIN_DOC - ! Compute ||x(lambda)||^2 - END_DOC - - implicit none - - ! Variables - - ! in - integer, intent(in) :: n - double precision, intent(in) :: e_val(n) - double precision, intent(in) :: tmp_wtg(n) - double precision, intent(in) :: lambda - - ! functions - double precision :: f_norm_trust_region_omp - - ! internal - double precision, allocatable :: tmp_fN(:) - integer :: i,j - - ! Allocation - allocate(tmp_fN(n)) - - call omp_set_max_active_levels(1) - - ! OMP - !$OMP PARALLEL & - !$OMP PRIVATE(i,j) & - !$OMP SHARED(n,lambda, e_val, thresh_eig,& - !$OMP tmp_fN, tmp_wtg, f_norm_trust_region_omp) & - !$OMP DEFAULT(NONE) - - ! Initialization - - !$OMP MASTER - f_norm_trust_region_omp = 0d0 - !$OMP END MASTER - - !$OMP DO - do i = 1, n - tmp_fN(i) = 0d0 - enddo - !$OMP END DO - - ! Calculations - !$OMP DO - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - tmp_fN(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**2 - endif - enddo - !$OMP END DO - - !$OMP MASTER - do i = 1, n - f_norm_trust_region_omp = f_norm_trust_region_omp + tmp_fN(i) - enddo - !$OMP END MASTER - - !$OMP END PARALLEL - - deallocate(tmp_fN) - -end function -#+END_SRC - -* First derivative of (||x||^2 - Delta^2)^2 -Version without omp - -*Function to compute the first derivative of ||x||^2 - Delta* - -This function computes the first derivative of (||x||^2 - Delta^2)^2 -with respect to lambda. - -\begin{align*} -\frac{\partial }{\partial \lambda} (||\textbf{x}(\lambda)||^2 - \Delta^2)^2 -= 2 \left(-2\sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^3} \right) -\left( - \Delta^2 + \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i+ \lambda)^2} \right) -\end{align*} - -\begin{align*} -\text{accu1} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^2} \\ -\text{accu2} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^3} -\end{align*} - -Provided: -| m_num | integer | number of MOs | - -Input: -| n | integer | mo_num*(mo_num-1)/2 | -| e_val(n) | double precision | eigenvalues of the hessian | -| W(n,n) | double precision | eigenvectors of the hessian | -| v_grad(n) | double precision | gradient | -| lambda | double precision | Lagrange multiplier | -| delta | double precision | Delta of the trust region | - -Internal: -| accu1 | double precision | first sum of the formula | -| accu2 | double precision | second sum of the formula | -| wtg | double precision | temporary variable to store W^T.v_grad | -| i,j | integer | indexes | - -Function: -| d1_norm_trust_region | double precision | first derivative with respect to lambda of (norm(x)^2 - Delta^2)^2 | -| ddot | double precision | blas dot product | - -#+BEGIN_SRC f90 :comments org :tangle trust_region_optimal_lambda.irp.f -function d1_norm_trust_region(n,e_val,w,v_grad,lambda,delta) - - include 'pi.h' - - BEGIN_DOC - ! Compute the first derivative with respect to lambda of (||x(lambda)||^2 - Delta^2)^2 - END_DOC - - implicit none - - ! Variables - - ! in - integer, intent(in) :: n - double precision, intent(in) :: e_val(n) - double precision, intent(in) :: w(n,n) - double precision, intent(in) :: v_grad(n) - double precision, intent(in) :: lambda - double precision, intent(in) :: delta - - ! Internal - double precision :: wtg, accu1, accu2 - integer :: i, j - - ! Functions - double precision :: d1_norm_trust_region - double precision :: ddot - - ! Initialization - accu1 = 0d0 - accu2 = 0d0 - - do i = 1, n - wtg = 0d0 - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - do j = 1, n - wtg = wtg + w(j,i) * v_grad(j) - enddo - !wtg = ddot(n,w(:,i),1,v_grad,1) - accu1 = accu1 + wtg**2 / (e_val(i) + lambda)**2 - endif - enddo - - do i = 1, n - wtg = 0d0 - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - do j = 1, n - wtg = wtg + w(j,i) * v_grad(j) - enddo - !wtg = ddot(n,w(:,i),1,v_grad,1) - accu2 = accu2 - 2d0 * wtg**2 / (e_val(i) + lambda)**3 - endif - enddo - - d1_norm_trust_region = 2d0 * accu2 * (accu1 - delta**2) - -end function -#+END_SRC - -* Second derivative of (||x||^2 - Delta^2)^2 -Version without OMP - -*Function to compute the second derivative of ||x||^2 - Delta* - - -\begin{equation} -\frac{\partial^2 }{\partial \lambda^2} (||\textbf{x}(\lambda)||^2 - \Delta^2)^2 -= 2 \left[ \left( \sum_{i=1}^n 6 \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^4} \right) \left( - \Delta^2 + \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^2} \right) + \left( \sum_{i=1}^n -2 \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^3} \right)^2 \right] -\end{equation} - -\begin{align*} -\text{accu1} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^2} \\ -\text{accu2} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^3} \\ -\text{accu3} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^4} -\end{align*} -Provided: -| m_num | integer | number of MOs | - -Input: -| n | integer | mo_num*(mo_num-1)/2 | -| e_val(n) | double precision | eigenvalues of the hessian | -| W(n,n) | double precision | eigenvectors of the hessian | -| v_grad(n) | double precision | gradient | -| lambda | double precision | Lagrange multiplier | -| delta | double precision | Delta of the trust region | - -Internal: -| accu1 | double precision | first sum of the formula | -| accu2 | double precision | second sum of the formula | -| accu3 | double precision | third sum of the formula | -| wtg | double precision | temporary variable to store W^T.v_grad | -| i,j | integer | indexes | - -Function: -| d2_norm_trust_region | double precision | second derivative with respect to lambda of norm(x)^2 - Delta^2 | -| ddot | double precision | blas dot product | - -#+BEGIN_SRC f90 :comments org :tangle trust_region_optimal_lambda.irp.f -function d2_norm_trust_region(n,e_val,w,v_grad,lambda,delta) - - include 'pi.h' - - BEGIN_DOC - ! Compute the second derivative with respect to lambda of (||x(lambda)||^2 - Delta^2)^2 - END_DOC - - implicit none - - ! Variables - - ! in - integer, intent(in) :: n - double precision, intent(in) :: e_val(n) - double precision, intent(in) :: w(n,n) - double precision, intent(in) :: v_grad(n) - double precision, intent(in) :: lambda - double precision, intent(in) :: delta - - ! Functions - double precision :: d2_norm_trust_region - double precision :: ddot - - ! Internal - double precision :: wtg,accu1,accu2,accu3 - integer :: i, j - - ! Initialization - accu1 = 0d0 - accu2 = 0d0 - accu3 = 0d0 - - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - wtg = 0d0 - do j = 1, n - wtg = wtg + w(j,i) * v_grad(j) - enddo - !wtg = ddot(n,w(:,i),1,v_grad,1) - accu1 = accu1 + wtg**2 / (e_val(i) + lambda)**2 !4 - endif - enddo - - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - wtg = 0d0 - do j = 1, n - wtg = wtg + w(j,i) * v_grad(j) - enddo - !wtg = ddot(n,w(:,i),1,v_grad,1) - accu2 = accu2 - 2d0 * wtg**2 / (e_val(i) + lambda)**3 !2 - endif - enddo - - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - wtg = 0d0 - do j = 1, n - wtg = wtg + w(j,i) * v_grad(j) - enddo - !wtg = ddot(n,w(:,i),1,v_grad,1) - accu3 = accu3 + 6d0 * wtg**2 / (e_val(i) + lambda)**4 !3 - endif - enddo - - d2_norm_trust_region = 2d0 * (accu3 * (- delta**2 + accu1) + accu2**2) - -end function -#+END_SRC - -* Function value of ||x||^2 -Version without OMP - -*Compute the value of ||x||^2* - -This function computes the value of ||x(lambda)||^2 - -\begin{align*} -||\textbf{x}(\lambda)||^2 = \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^2} -\end{align*} - -Provided: -| m_num | integer | number of MOs | - -Input: -| n | integer | mo_num*(mo_num-1)/2 | -| e_val(n) | double precision | eigenvalues of the hessian | -| W(n,n) | double precision | eigenvectors of the hessian | -| v_grad(n) | double precision | gradient | -| lambda | double precision | Lagrange multiplier | -| delta | double precision | Delta of the trust region | - -Internal: -| wtg | double precision | temporary variable to store W^T.v_grad | -| i,j | integer | indexes | - -Function: -| f_norm_trust_region | double precision | value of norm(x)^2 | -| ddot | double precision | blas dot product | - - -#+BEGIN_SRC f90 :comments org :tangle trust_region_optimal_lambda.irp.f -function f_norm_trust_region(n,e_val,tmp_wtg,lambda) - - include 'pi.h' - - BEGIN_DOC - ! Compute ||x(lambda)||^2 - END_DOC - - implicit none - - ! Variables - - ! in - integer, intent(in) :: n - double precision, intent(in) :: e_val(n) - double precision, intent(in) :: tmp_wtg(n) - double precision, intent(in) :: lambda - - ! function - double precision :: f_norm_trust_region - double precision :: ddot - - ! internal - integer :: i,j - - ! Initialization - f_norm_trust_region = 0d0 - - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - f_norm_trust_region = f_norm_trust_region + tmp_wtg(i)**2 / (e_val(i) + lambda)**2 - endif - enddo - -end function -#+END_SRC - -* OMP: First derivative of (1/||x||^2 - 1/Delta^2)^2 -Version with OMP - -*Compute the first derivative of (1/||x||^2 - 1/Delta^2)^2* - -This function computes the value of (1/||x(lambda)||^2 - 1/Delta^2)^2 - -\begin{align*} - \frac{\partial}{\partial \lambda} (1/||\textbf{x}(\lambda)||^2 - 1/\Delta^2)^2 - &= 4 \frac{\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3}} - {(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^3} - - \frac{4}{\Delta^2} \frac{\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3)}} - {(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^2} \\ - &= 4 \sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3} - \left( \frac{1}{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^3} - - \frac{1}{\Delta^2 (\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^2} \right) -\end{align*} - -\begin{align*} -\text{accu1} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^2} \\ -\text{accu2} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^3} -\end{align*} - -Provided: -| m_num | integer | number of MOs | - -Input: -| n | integer | mo_num*(mo_num-1)/2 | -| e_val(n) | double precision | eigenvalues of the hessian | -| W(n,n) | double precision | eigenvectors of the hessian | -| v_grad(n) | double precision | gradient | -| lambda | double precision | Lagrange multiplier | -| delta | double precision | Delta of the trust region | - -Internal: -| wtg | double precision | temporary variable to store W^T.v_grad | -| tmp_accu1 | double precision | temporary array for the first sum | -| tmp_accu2 | double precision | temporary array for the second sum | -| tmp_wtg(n) | double precision | temporary array for W^t.v_grad | -| i,j | integer | indexes | - -Function: -| d1_norm_inverse_trust_region | double precision | value of the first derivative | - -#+BEGIN_SRC f90 :comments org :tangle trust_region_optimal_lambda.irp.f -function d1_norm_inverse_trust_region_omp(n,e_val,tmp_wtg,lambda,delta) - - use omp_lib - include 'pi.h' - - BEGIN_DOC - ! Compute the first derivative of (1/||x||^2 - 1/Delta^2)^2 - END_DOC - - implicit none - - ! Variables - - ! in - integer, intent(in) :: n - double precision, intent(in) :: e_val(n) - double precision, intent(in) :: tmp_wtg(n) - double precision, intent(in) :: lambda - double precision, intent(in) :: delta - - ! Internal - double precision :: accu1, accu2 - integer :: i,j - double precision, allocatable :: tmp_accu1(:), tmp_accu2(:) - - ! Functions - double precision :: d1_norm_inverse_trust_region_omp - - ! Allocation - allocate(tmp_accu1(n), tmp_accu2(n)) - - ! OMP - call omp_set_max_active_levels(1) - - ! OMP - !$OMP PARALLEL & - !$OMP PRIVATE(i,j) & - !$OMP SHARED(n,lambda, e_val, thresh_eig,& - !$OMP tmp_accu1, tmp_accu2, tmp_wtg, accu1, accu2) & - !$OMP DEFAULT(NONE) - - !$OMP MASTER - accu1 = 0d0 - accu2 = 0d0 - !$OMP END MASTER - - !$OMP DO - do i = 1, n - tmp_accu1(i) = 0d0 - enddo - !$OMP END DO - - !$OMP DO - do i = 1, n - tmp_accu2(i) = 0d0 - enddo - !$OMP END DO - -! !$OMP MASTER -! do i = 1, n -! if (ABS(e_val(i)+lambda) > 1d-12) then -! tmp_accu1(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**2 -! endif -! enddo -! !$OMP END MASTER - - !$OMP DO - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - tmp_accu1(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**2 - endif - enddo - !$OMP END DO - - !$OMP MASTER - do i = 1, n - accu1 = accu1 + tmp_accu1(i) - enddo - !$OMP END MASTER - -! !$OMP MASTER -! do i = 1, n -! if (ABS(e_val(i)+lambda) > 1d-12) then -! tmp_accu2(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**3 -! endif -! enddo -! !$OMP END MASTER - - !$OMP DO - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - tmp_accu2(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**3 - endif - enddo - !$OMP END DO - - !$OMP MASTER - do i = 1, n - accu2 = accu2 + tmp_accu2(i) - enddo - !$OMP END MASTER - - !$OMP END PARALLEL - - call omp_set_max_active_levels(4) - - d1_norm_inverse_trust_region_omp = 4d0 * accu2 * (1d0/accu1**3 - 1d0/(delta**2 * accu1**2)) - - deallocate(tmp_accu1, tmp_accu2) - -end -#+END_SRC - -* OMP: Second derivative of (1/||x||^2 - 1/Delta^2)^2 -Version with OMP - -*Compute the first derivative of (1/||x||^2 - 1/Delta^2)^2* - -This function computes the value of (1/||x(lambda)||^2 - 1/Delta^2)^2 - -\begin{align*} - \frac{\partial^2}{\partial \lambda^2} (1/||\textbf{x}(\lambda)||^2 - 1/\Delta^2)^2 - &= 4 \left[ \frac{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3)})^2}{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^4} - - 3 \frac{\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^4}}{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^3} \right] \\ - &- \frac{4}{\Delta^2} \left[ \frac{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3)})^2}{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^3} - - 3 \frac{\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^4}}{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^2} \right] -\end{align*} - - -\begin{align*} -\text{accu1} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^2} \\ -\text{accu2} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^3} \\ -\text{accu3} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^4} -\end{align*} - -Provided: -| m_num | integer | number of MOs | - -Input: -| n | integer | mo_num*(mo_num-1)/2 | -| e_val(n) | double precision | eigenvalues of the hessian | -| W(n,n) | double precision | eigenvectors of the hessian | -| v_grad(n) | double precision | gradient | -| lambda | double precision | Lagrange multiplier | -| delta | double precision | Delta of the trust region | - -Internal: -| wtg | double precision | temporary variable to store W^T.v_grad | -| tmp_accu1 | double precision | temporary array for the first sum | -| tmp_accu2 | double precision | temporary array for the second sum | -| tmp_wtg(n) | double precision | temporary array for W^t.v_grad | -| i,j | integer | indexes | - -Function: -| d1_norm_inverse_trust_region | double precision | value of the first derivative | - -#+BEGIN_SRC f90 :comments org :tangle trust_region_optimal_lambda.irp.f -function d2_norm_inverse_trust_region_omp(n,e_val,tmp_wtg,lambda,delta) - - use omp_lib - include 'pi.h' - - BEGIN_DOC - ! Compute the second derivative of (1/||x||^2 - 1/Delta^2)^2 - END_DOC - - implicit none - - ! Variables - - ! in - integer, intent(in) :: n - double precision, intent(in) :: e_val(n) - double precision, intent(in) :: tmp_wtg(n) - double precision, intent(in) :: lambda - double precision, intent(in) :: delta - - ! Internal - double precision :: accu1, accu2, accu3 - integer :: i,j - double precision, allocatable :: tmp_accu1(:), tmp_accu2(:), tmp_accu3(:) - - ! Functions - double precision :: d2_norm_inverse_trust_region_omp - - ! Allocation - allocate(tmp_accu1(n), tmp_accu2(n), tmp_accu3(n)) - - ! OMP - call omp_set_max_active_levels(1) - - ! OMP - !$OMP PARALLEL & - !$OMP PRIVATE(i,j) & - !$OMP SHARED(n,lambda, e_val, thresh_eig,& - !$OMP tmp_accu1, tmp_accu2, tmp_accu3, tmp_wtg, & - !$OMP accu1, accu2, accu3) & - !$OMP DEFAULT(NONE) - - !$OMP MASTER - accu1 = 0d0 - accu2 = 0d0 - accu3 = 0d0 - !$OMP END MASTER - - !$OMP DO - do i = 1, n - tmp_accu1(i) = 0d0 - enddo - !$OMP END DO - - !$OMP DO - do i = 1, n - tmp_accu2(i) = 0d0 - enddo - !$OMP END DO - - !$OMP DO - do i = 1, n - tmp_accu3(i) = 0d0 - enddo - !$OMP END DO - - !$OMP DO - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - tmp_accu1(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**2 - endif - enddo - !$OMP END DO - - !$OMP MASTER - do i = 1, n - accu1 = accu1 + tmp_accu1(i) - enddo - !$OMP END MASTER - - !$OMP DO - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - tmp_accu2(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**3 - endif - enddo - !$OMP END DO - - !$OMP MASTER - do i = 1, n - accu2 = accu2 + tmp_accu2(i) - enddo - !$OMP END MASTER - - !$OMP DO - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - tmp_accu3(i) = tmp_wtg(i)**2 / (e_val(i) + lambda)**4 - endif - enddo - !$OMP END DO - - !$OMP MASTER - do i = 1, n - accu3 = accu3 + tmp_accu3(i) - enddo - !$OMP END MASTER - - !$OMP END PARALLEL - - call omp_set_max_active_levels(4) - - d2_norm_inverse_trust_region_omp = 4d0 * (6d0 * accu2**2/accu1**4 - 3d0 * accu3/accu1**3) & - - 4d0/delta**2 * (4d0 * accu2**2/accu1**3 - 3d0 * accu3/accu1**2) - - deallocate(tmp_accu1,tmp_accu2,tmp_accu3) - -end -#+END_SRC - -* First derivative of (1/||x||^2 - 1/Delta^2)^2 -Version without OMP - -*Compute the first derivative of (1/||x||^2 - 1/Delta^2)^2* - -This function computes the value of (1/||x(lambda)||^2 - 1/Delta^2)^2 - -\begin{align*} - \frac{\partial}{\partial \lambda} (1/||\textbf{x}(\lambda)||^2 - 1/\Delta^2)^2 - &= 4 \frac{\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3}} - {(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^3} - - \frac{4}{\Delta^2} \frac{\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3)}} - {(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^2} \\ - &= 4 \sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3} - \left( \frac{1}{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^3} - - \frac{1}{\Delta^2 (\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^2} \right) -\end{align*} -\begin{align*} -\text{accu1} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^2} \\ -\text{accu2} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^3} -\end{align*} -Provided: -| m_num | integer | number of MOs | - -Input: -| n | integer | mo_num*(mo_num-1)/2 | -| e_val(n) | double precision | eigenvalues of the hessian | -| W(n,n) | double precision | eigenvectors of the hessian | -| v_grad(n) | double precision | gradient | -| lambda | double precision | Lagrange multiplier | -| delta | double precision | Delta of the trust region | - -Internal: -| wtg | double precision | temporary variable to store W^T.v_grad | -| i,j | integer | indexes | - -Function: -| d1_norm_inverse_trust_region | double precision | value of the first derivative | - -#+BEGIN_SRC f90 :comments org :tangle trust_region_optimal_lambda.irp.f -function d1_norm_inverse_trust_region(n,e_val,w,v_grad,lambda,delta) - - include 'pi.h' - - BEGIN_DOC - ! Compute the first derivative of (1/||x||^2 - 1/Delta^2)^2 - END_DOC - - implicit none - - ! Variables - - ! in - integer, intent(in) :: n - double precision, intent(in) :: e_val(n) - double precision, intent(in) :: w(n,n) - double precision, intent(in) :: v_grad(n) - double precision, intent(in) :: lambda - double precision, intent(in) :: delta - - ! Internal - double precision :: wtg, accu1, accu2 - integer :: i,j - - ! Functions - double precision :: d1_norm_inverse_trust_region - - accu1 = 0d0 - accu2 = 0d0 - - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - wtg = 0d0 - do j = 1, n - wtg = wtg + w(j,i) * v_grad(j) - enddo - accu1 = accu1 + wtg**2 / (e_val(i) + lambda)**2 - endif - enddo - - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - wtg = 0d0 - do j = 1, n - wtg = wtg + w(j,i) * v_grad(j) - enddo - accu2 = accu2 + wtg**2 / (e_val(i) + lambda)**3 - endif - enddo - - d1_norm_inverse_trust_region = 4d0 * accu2 * (1d0/accu1**3 - 1d0/(delta**2 * accu1**2)) - -end -#+END_SRC - -* Second derivative of (1/||x||^2 - 1/Delta^2)^2 -Version without OMP - -*Compute the second derivative of (1/||x||^2 - 1/Delta^2)^2* - -This function computes the value of (1/||x(lambda)||^2 - 1/Delta^2)^2 - -\begin{align*} - \frac{\partial^2}{\partial \lambda^2} (1/||\textbf{x}(\lambda)||^2 - 1/\Delta^2)^2 - &= 4 \left[ \frac{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3)})^2}{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^4} - - 3 \frac{\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^4}}{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^3} \right] \\ - &- \frac{4}{\Delta^2} \left[ \frac{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^3)})^2}{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^3} - - 3 \frac{\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^4}}{(\sum_i \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2})^2} \right] -\end{align*} - -\begin{align*} -\text{accu1} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^2} \\ -\text{accu2} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^3} \\ -\text{accu3} &= \sum_{i=1}^n \frac{(\textbf{w}_i^T \textbf{g})^2}{(h_i + \lambda)^4} -\end{align*} - -Provided: -| m_num | integer | number of MOs | - -Input: -| n | integer | mo_num*(mo_num-1)/2 | -| e_val(n) | double precision | eigenvalues of the hessian | -| W(n,n) | double precision | eigenvectors of the hessian | -| v_grad(n) | double precision | gradient | -| lambda | double precision | Lagrange multiplier | -| delta | double precision | Delta of the trust region | - -Internal: -| wtg | double precision | temporary variable to store W^T.v_grad | -| i,j | integer | indexes | - -Function: -| d2_norm_inverse_trust_region | double precision | value of the first derivative | - -#+BEGIN_SRC f90 :comments org :tangle trust_region_optimal_lambda.irp.f -function d2_norm_inverse_trust_region(n,e_val,w,v_grad,lambda,delta) - - include 'pi.h' - - BEGIN_DOC - ! Compute the second derivative of (1/||x||^2 - 1/Delta^2)^2 - END_DOC - - implicit none - - ! Variables - - ! in - integer, intent(in) :: n - double precision, intent(in) :: e_val(n) - double precision, intent(in) :: w(n,n) - double precision, intent(in) :: v_grad(n) - double precision, intent(in) :: lambda - double precision, intent(in) :: delta - - ! Internal - double precision :: wtg, accu1, accu2, accu3 - integer :: i,j - - ! Functions - double precision :: d2_norm_inverse_trust_region - - accu1 = 0d0 - accu2 = 0d0 - accu3 = 0d0 - - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - wtg = 0d0 - do j = 1, n - wtg = wtg + w(j,i) * v_grad(j) - enddo - accu1 = accu1 + wtg**2 / (e_val(i) + lambda)**2 - endif - enddo - - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - wtg = 0d0 - do j = 1, n - wtg = wtg + w(j,i) * v_grad(j) - enddo - accu2 = accu2 + wtg**2 / (e_val(i) + lambda)**3 - endif - enddo - - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - wtg = 0d0 - do j = 1, n - wtg = wtg + w(j,i) * v_grad(j) - enddo - accu3 = accu3 + wtg**2 / (e_val(i) + lambda)**4 - endif - enddo - - d2_norm_inverse_trust_region = 4d0 * (6d0 * accu2**2/accu1**4 - 3d0 * accu3/accu1**3) & - - 4d0/delta**2 * (4d0 * accu2**2/accu1**3 - 3d0 * accu3/accu1**2) - -end -#+END_SRC diff --git a/src/utils_trust_region/trust_region_rho.irp.f b/src/utils_trust_region/trust_region_rho.irp.f deleted file mode 100644 index 45738736..00000000 --- a/src/utils_trust_region/trust_region_rho.irp.f +++ /dev/null @@ -1,121 +0,0 @@ -! Agreement with the model: Rho - -! *Compute the ratio : rho = (prev_energy - energy) / (prev_energy - e_model)* - -! Rho represents the agreement between the model (the predicted energy -! by the Taylor expansion truncated at the 2nd order) and the real -! energy : - -! \begin{equation} -! \rho^{k+1} = \frac{E^{k} - E^{k+1}}{E^{k} - m^{k+1}} -! \end{equation} -! With : -! $E^{k}$ the energy at the previous iteration -! $E^{k+1}$ the energy at the actual iteration -! $m^{k+1}$ the predicted energy for the actual iteration -! (cf. trust_e_model) - -! If $\rho \approx 1$, the agreement is good, contrary to $\rho \approx 0$. -! If $\rho \leq 0$ the previous energy is lower than the actual -! energy. We have to cancel the last step and use a smaller trust -! region. -! Here we cancel the last step if $\rho < 0.1$, because even if -! the energy decreases, the agreement is bad, i.e., the Taylor expansion -! truncated at the second order doesn't represent correctly the energy -! landscape. So it's better to cancel the step and restart with a -! smaller trust region. - -! Provided in qp_edit: -! | thresh_rho | - -! Input: -! | prev_energy | double precision | previous energy (energy before the rotation) | -! | e_model | double precision | predicted energy after the rotation | - -! Output: -! | rho | double precision | the agreement between the model (predicted) and the real energy | -! | prev_energy | double precision | if rho >= 0.1 the actual energy becomes the previous energy | -! | | | else the previous energy doesn't change | - -! Internal: -! | energy | double precision | energy (real) after the rotation | -! | i | integer | index | -! | t* | double precision | time | - - -subroutine trust_region_rho(prev_energy, energy,e_model,rho) - - include 'pi.h' - - BEGIN_DOC - ! Compute rho, the agreement between the predicted criterion/energy and the real one - END_DOC - - implicit none - - ! Variables - - ! In - double precision, intent(inout) :: prev_energy - double precision, intent(in) :: e_model, energy - - ! Out - double precision, intent(out) :: rho - - ! Internal - double precision :: t1, t2, t3 - integer :: i - - print*,'' - print*,'---Rho_model---' - - call wall_time(t1) - -! Rho -! \begin{equation} -! \rho^{k+1} = \frac{E^{k} - E^{k+1}}{E^{k} - m^{k+1}} -! \end{equation} - -! In function of $\rho$ th step can be accepted or cancelled. - -! If we cancel the last step (k+1), the previous energy (k) doesn't -! change! -! If the step (k+1) is accepted, then the "previous energy" becomes E(k+1) - - -! Already done in an other subroutine - !if (ABS(prev_energy - e_model) < 1d-12) then - ! print*,'WARNING: prev_energy - e_model < 1d-12' - ! print*,'=> rho will tend toward infinity' - ! print*,'Check you convergence criterion !' - !endif - - rho = (prev_energy - energy) / (prev_energy - e_model) - - print*, 'previous energy, prev_energy :', prev_energy - print*, 'predicted energy, e_model :', e_model - print*, 'real energy, energy :', energy - print*, 'prev_energy - energy :', prev_energy - energy - print*, 'prev_energy - e_model :', prev_energy - e_model - print*, 'Rho :', rho - print*, 'Threshold for rho:', thresh_rho - - ! Modification of prev_energy in function of rho - if (rho < thresh_rho) then !0.1) then - ! the step is cancelled - print*, 'Rho <', thresh_rho,', the previous energy does not changed' - print*, 'prev_energy :', prev_energy - else - ! the step is accepted - prev_energy = energy - print*, 'Rho >=', thresh_rho,', energy -> prev_energy :', energy - endif - - call wall_time(t2) - t3 = t2 - t1 - print*,'Time in rho model:', t3 - - print*,'---End rho_model---' - print*,'' - -end subroutine diff --git a/src/utils_trust_region/trust_region_rho.org b/src/utils_trust_region/trust_region_rho.org deleted file mode 100644 index 9b25ee29..00000000 --- a/src/utils_trust_region/trust_region_rho.org +++ /dev/null @@ -1,123 +0,0 @@ -* Agreement with the model: Rho - -*Compute the ratio : rho = (prev_energy - energy) / (prev_energy - e_model)* - -Rho represents the agreement between the model (the predicted energy -by the Taylor expansion truncated at the 2nd order) and the real -energy : - -\begin{equation} -\rho^{k+1} = \frac{E^{k} - E^{k+1}}{E^{k} - m^{k+1}} -\end{equation} -With : -$E^{k}$ the energy at the previous iteration -$E^{k+1}$ the energy at the actual iteration -$m^{k+1}$ the predicted energy for the actual iteration -(cf. trust_e_model) - -If $\rho \approx 1$, the agreement is good, contrary to $\rho \approx 0$. -If $\rho \leq 0$ the previous energy is lower than the actual -energy. We have to cancel the last step and use a smaller trust -region. -Here we cancel the last step if $\rho < 0.1$, because even if -the energy decreases, the agreement is bad, i.e., the Taylor expansion -truncated at the second order doesn't represent correctly the energy -landscape. So it's better to cancel the step and restart with a -smaller trust region. - -Provided in qp_edit: -| thresh_rho | - -Input: -| prev_energy | double precision | previous energy (energy before the rotation) | -| e_model | double precision | predicted energy after the rotation | - -Output: -| rho | double precision | the agreement between the model (predicted) and the real energy | -| prev_energy | double precision | if rho >= 0.1 the actual energy becomes the previous energy | -| | | else the previous energy doesn't change | - -Internal: -| energy | double precision | energy (real) after the rotation | -| i | integer | index | -| t* | double precision | time | - -#+BEGIN_SRC f90 :comments org :tangle trust_region_rho.irp.f -subroutine trust_region_rho(prev_energy, energy,e_model,rho) - - include 'pi.h' - - BEGIN_DOC - ! Compute rho, the agreement between the predicted criterion/energy and the real one - END_DOC - - implicit none - - ! Variables - - ! In - double precision, intent(inout) :: prev_energy - double precision, intent(in) :: e_model, energy - - ! Out - double precision, intent(out) :: rho - - ! Internal - double precision :: t1, t2, t3 - integer :: i - - print*,'' - print*,'---Rho_model---' - - call wall_time(t1) -#+END_SRC - -** Rho -\begin{equation} -\rho^{k+1} = \frac{E^{k} - E^{k+1}}{E^{k} - m^{k+1}} -\end{equation} - -In function of $\rho$ th step can be accepted or cancelled. - -If we cancel the last step (k+1), the previous energy (k) doesn't -change! -If the step (k+1) is accepted, then the "previous energy" becomes E(k+1) - -#+BEGIN_SRC f90 :comments org :tangle trust_region_rho.irp.f - ! Already done in an other subroutine - !if (ABS(prev_energy - e_model) < 1d-12) then - ! print*,'WARNING: prev_energy - e_model < 1d-12' - ! print*,'=> rho will tend toward infinity' - ! print*,'Check you convergence criterion !' - !endif - - rho = (prev_energy - energy) / (prev_energy - e_model) - - print*, 'previous energy, prev_energy :', prev_energy - print*, 'predicted energy, e_model :', e_model - print*, 'real energy, energy :', energy - print*, 'prev_energy - energy :', prev_energy - energy - print*, 'prev_energy - e_model :', prev_energy - e_model - print*, 'Rho :', rho - print*, 'Threshold for rho:', thresh_rho - - ! Modification of prev_energy in function of rho - if (rho < thresh_rho) then !0.1) then - ! the step is cancelled - print*, 'Rho <', thresh_rho,', the previous energy does not changed' - print*, 'prev_energy :', prev_energy - else - ! the step is accepted - prev_energy = energy - print*, 'Rho >=', thresh_rho,', energy -> prev_energy :', energy - endif - - call wall_time(t2) - t3 = t2 - t1 - print*,'Time in rho model:', t3 - - print*,'---End rho_model---' - print*,'' - -end subroutine -#+END_SRC diff --git a/src/utils_trust_region/trust_region_step.irp.f b/src/utils_trust_region/trust_region_step.irp.f deleted file mode 100644 index 42aa6ed4..00000000 --- a/src/utils_trust_region/trust_region_step.irp.f +++ /dev/null @@ -1,716 +0,0 @@ -! Trust region - -! *Compute the next step with the trust region algorithm* - -! The Newton method is an iterative method to find a minimum of a given -! function. It uses a Taylor series truncated at the second order of the -! targeted function and gives its minimizer. The minimizer is taken as -! the new position and the same thing is done. And by doing so -! iteratively the method find a minimum, a local or global one depending -! of the starting point and the convexity/nonconvexity of the targeted -! function. - -! The goal of the trust region is to constrain the step size of the -! Newton method in a certain area around the actual position, where the -! Taylor series is a good approximation of the targeted function. This -! area is called the "trust region". - -! In addition, in function of the agreement between the Taylor -! development of the energy and the real energy, the size of the trust -! region will be updated at each iteration. By doing so, the step sizes -! are not too larges. In addition, since we add a criterion to cancel the -! step if the energy increases (more precisely if rho < 0.1), so it's -! impossible to diverge. \newline - -! References: \newline -! Nocedal & Wright, Numerical Optimization, chapter 4 (1999), \newline -! https://link.springer.com/book/10.1007/978-0-387-40065-5, \newline -! ISBN: 978-0-387-40065-5 \newline - -! By using the first and the second derivatives, the Newton method gives -! a step: -! \begin{align*} -! \textbf{x}_{(k+1)}^{\text{Newton}} = - \textbf{H}_{(k)}^{-1} \cdot -! \textbf{g}_{(k)} -! \end{align*} -! which leads to the minimizer of the Taylor series. -! !!! Warning: the Newton method gives the minimizer if and only if -! $\textbf{H}$ is positive definite, else it leads to a saddle point !!! -! But we want a step $\textbf{x}_{(k+1)}$ with a constraint on its (euclidian) norm: -! \begin{align*} -! ||\textbf{x}_{(k+1)}|| \leq \Delta_{(k+1)} -! \end{align*} -! which is equivalent to -! \begin{align*} -! \textbf{x}_{(k+1)}^T \cdot \textbf{x}_{(k+1)} \leq \Delta_{(k+1)}^2 -! \end{align*} - -! with: \newline -! $\textbf{x}_{(k+1)}$ is the step for the k+1-th iteration (vector of -! size n) \newline -! $\textbf{H}_{(k)}$ is the hessian at the k-th iteration (n by n -! matrix) \newline -! $\textbf{g}_{(k)}$ is the gradient at the k-th iteration (vector of -! size n) \newline -! $\Delta_{(k+1)}$ is the trust radius for the (k+1)-th iteration -! \newline - -! Thus we want to constrain the step size $\textbf{x}_{(k+1)}$ into a -! hypersphere of radius $\Delta_{(k+1)}$.\newline - -! So, if $||\textbf{x}_{(k+1)}^{\text{Newton}}|| \leq \Delta_{(k)}$ and -! $\textbf{H}$ is positive definite, the -! solution is the step given by the Newton method -! $\textbf{x}_{(k+1)} = \textbf{x}_{(k+1)}^{\text{Newton}}$. -! Else we have to constrain the step size. For simplicity we will remove -! the index $_{(k)}$ and $_{(k+1)}$. To restict the step size, we have -! to put a constraint on $\textbf{x}$ with a Lagrange multiplier. -! Starting from the Taylor series of a function E (here, the energy) -! truncated at the 2nd order, we have: -! \begin{align*} -! E(\textbf{x}) = E +\textbf{g}^T \cdot \textbf{x} + \frac{1}{2} -! \cdot \textbf{x}^T \cdot \textbf{H} \cdot \textbf{x} + -! \mathcal{O}(\textbf{x}^2) -! \end{align*} - -! With the constraint on the norm of $\textbf{x}$ we can write the -! Lagrangian -! \begin{align*} -! \mathcal{L}(\textbf{x},\lambda) = E + \textbf{g}^T \cdot \textbf{x} -! + \frac{1}{2} \cdot \textbf{x}^T \cdot \textbf{H} \cdot \textbf{x} -! + \frac{1}{2} \lambda (\textbf{x}^T \cdot \textbf{x} - \Delta^2) -! \end{align*} -! Where: \newline -! $\lambda$ is the Lagrange multiplier \newline -! $E$ is the energy at the k-th iteration $\Leftrightarrow -! E(\textbf{x} = \textbf{0})$ \newline - -! To solve this equation, we search a stationary point where the first -! derivative of $\mathcal{L}$ with respect to $\textbf{x}$ becomes 0, i.e. -! \begin{align*} -! \frac{\partial \mathcal{L}(\textbf{x},\lambda)}{\partial \textbf{x}}=0 -! \end{align*} - -! The derivative is: -! \begin{align*} -! \frac{\partial \mathcal{L}(\textbf{x},\lambda)}{\partial \textbf{x}} -! = \textbf{g} + \textbf{H} \cdot \textbf{x} + \lambda \cdot \textbf{x} -! \end{align*} - -! So, we search $\textbf{x}$ such as: -! \begin{align*} -! \frac{\partial \mathcal{L}(\textbf{x},\lambda)}{\partial \textbf{x}} -! = \textbf{g} + \textbf{H} \cdot \textbf{x} + \lambda \cdot \textbf{x} = 0 -! \end{align*} - -! We can rewrite that as: -! \begin{align*} -! \textbf{g} + \textbf{H} \cdot \textbf{x} + \lambda \cdot \textbf{x} -! = \textbf{g} + (\textbf{H} +\textbf{I} \lambda) \cdot \textbf{x} = 0 -! \end{align*} -! with $\textbf{I}$ is the identity matrix. - -! By doing so, the solution is: -! \begin{align*} -! (\textbf{H} +\textbf{I} \lambda) \cdot \textbf{x}= -\textbf{g} -! \end{align*} -! \begin{align*} -! \textbf{x}= - (\textbf{H} + \textbf{I} \lambda)^{-1} \cdot \textbf{g} -! \end{align*} -! with $\textbf{x}^T \textbf{x} = \Delta^2$. - -! We have to solve this previous equation to find this $\textbf{x}$ in the -! trust region, i.e. $||\textbf{x}|| = \Delta$. Now, this problem is -! just a one dimension problem because we can express $\textbf{x}$ as a -! function of $\lambda$: -! \begin{align*} -! \textbf{x}(\lambda) = - (\textbf{H} + \textbf{I} \lambda)^{-1} \cdot \textbf{g} -! \end{align*} - -! We start from the fact that the hessian is diagonalizable. So we have: -! \begin{align*} -! \textbf{H} = \textbf{W} \cdot \textbf{h} \cdot \textbf{W}^T -! \end{align*} -! with: \newline -! $\textbf{H}$, the hessian matrix \newline -! $\textbf{W}$, the matrix containing the eigenvectors \newline -! $\textbf{w}_i$, the i-th eigenvector, i.e. i-th column of $\textbf{W}$ \newline -! $\textbf{h}$, the matrix containing the eigenvalues in ascending order \newline -! $h_i$, the i-th eigenvalue in ascending order \newline - -! Now we use the fact that adding a constant on the diagonal just shifts -! the eigenvalues: -! \begin{align*} -! \textbf{H} + \textbf{I} \lambda = \textbf{W} \cdot (\textbf{h} -! +\textbf{I} \lambda) \cdot \textbf{W}^T -! \end{align*} - -! By doing so we can express $\textbf{x}$ as a function of $\lambda$ -! \begin{align*} -! \textbf{x}(\lambda) = - \sum_{i=1}^n \frac{\textbf{w}_i^T \cdot -! \textbf{g}}{h_i + \lambda} \cdot \textbf{w}_i -! \end{align*} -! with $\lambda \neq - h_i$. - -! An interesting thing in our case is the norm of $\textbf{x}$, -! because we want $||\textbf{x}|| = \Delta$. Due to the orthogonality of -! the eigenvectors $\left\{\textbf{w} \right\} _{i=1}^n$ we have: -! \begin{align*} -! ||\textbf{x}(\lambda)||^2 = \sum_{i=1}^n \frac{(\textbf{w}_i^T \cdot -! \textbf{g})^2}{(h_i + \lambda)^2} -! \end{align*} - -! So the $||\textbf{x}(\lambda)||^2$ is just a function of $\lambda$. -! And if we study the properties of this function we see that: -! \begin{align*} -! \lim_{\lambda\to\infty} ||\textbf{x}(\lambda)|| = 0 -! \end{align*} -! and if $\textbf{w}_i^T \cdot \textbf{g} \neq 0$: -! \begin{align*} -! \lim_{\lambda\to -h_i} ||\textbf{x}(\lambda)|| = + \infty -! \end{align*} - -! From these limits and knowing that $h_1$ is the lowest eigenvalue, we -! can conclude that $||\textbf{x}(\lambda)||$ is a continuous and -! strictly decreasing function on the interval $\lambda \in -! (-h_1;\infty)$. Thus, there is one $\lambda$ in this interval which -! gives $||\textbf{x}(\lambda)|| = \Delta$, consequently there is one -! solution. - -! Since $\textbf{x} = - (\textbf{H} + \lambda \textbf{I})^{-1} \cdot -! \textbf{g}$ and we want to reduce the norm of $\textbf{x}$, clearly, -! $\lambda > 0$ ($\lambda = 0$ is the unconstraint solution). But the -! Newton method is only defined for a positive definite hessian matrix, -! so $(\textbf{H} + \textbf{I} \lambda)$ must be positive -! definite. Consequently, in the case where $\textbf{H}$ is not positive -! definite, to ensure the positive definiteness, $\lambda$ must be -! greater than $- h_1$. -! \begin{align*} -! \lambda > 0 \quad \text{and} \quad \lambda \geq - h_1 -! \end{align*} - -! From that there are five cases: -! - if $\textbf{H}$ is positive definite, $-h_1 < 0$, $\lambda \in (0,\infty)$ -! - if $\textbf{H}$ is not positive definite and $\textbf{w}_1^T \cdot -! \textbf{g} \neq 0$, $(\textbf{H} + \textbf{I} -! \lambda)$ -! must be positve definite, $-h_1 > 0$, $\lambda \in (-h_1, \infty)$ -! - if $\textbf{H}$ is not positive definite , $\textbf{w}_1^T \cdot -! \textbf{g} = 0$ and $||\textbf{x}(-h_1)|| > \Delta$ by removing -! $j=1$ in the sum, $(\textbf{H} + \textbf{I} \lambda)$ must be -! positive definite, $-h_1 > 0$, $\lambda \in (-h_1, \infty$) -! - if $\textbf{H}$ is not positive definite , $\textbf{w}_1^T \cdot -! \textbf{g} = 0$ and $||\textbf{x}(-h_1)|| \leq \Delta$ by removing -! $j=1$ in the sum, $(\textbf{H} + \textbf{I} \lambda)$ must be -! positive definite, $-h_1 > 0$, $\lambda = -h_1$). This case is -! similar to the case where $\textbf{H}$ and $||\textbf{x}(\lambda = -! 0)|| \leq \Delta$ -! but we can also add to $\textbf{x}$, the first eigenvector $\textbf{W}_1$ -! time a constant to ensure the condition $||\textbf{x}(\lambda = -! -h_1)|| = \Delta$ and escape from the saddle point - -! Thus to find the solution, we can write: -! \begin{align*} -! ||\textbf{x}(\lambda)|| = \Delta -! \end{align*} -! \begin{align*} -! ||\textbf{x}(\lambda)|| - \Delta = 0 -! \end{align*} - -! Taking the square of this equation -! \begin{align*} -! (||\textbf{x}(\lambda)|| - \Delta)^2 = 0 -! \end{align*} -! we have a function with one minimum for the optimal $\lambda$. -! Since we have the formula of $||\textbf{x}(\lambda)||^2$, we solve -! \begin{align*} -! (||\textbf{x}(\lambda)||^2 - \Delta^2)^2 = 0 -! \end{align*} - -! But in practice, it is more effective to solve: -! \begin{align*} -! (\frac{1}{||\textbf{x}(\lambda)||^2} - \frac{1}{\Delta^2})^2 = 0 -! \end{align*} - -! To do that, we just use the Newton method with "trust_newton" using -! first and second derivative of $(||\textbf{x}(\lambda)||^2 - -! \Delta^2)^2$ with respect to $\textbf{x}$. -! This will give the optimal $\lambda$ to compute the -! solution $\textbf{x}$ with the formula seen previously: -! \begin{align*} -! \textbf{x}(\lambda) = - \sum_{i=1}^n \frac{\textbf{w}_i^T \cdot -! \textbf{g}}{h_i + \lambda} \cdot \textbf{w}_i -! \end{align*} - -! The solution $\textbf{x}(\lambda)$ with the optimal $\lambda$ is our -! step to go from the (k)-th to the (k+1)-th iteration, is noted $\textbf{x}^*$. - - - - -! Evolution of the trust region - -! We initialize the trust region at the first iteration using a radius -! \begin{align*} -! \Delta = ||\textbf{x}(\lambda=0)|| -! \end{align*} - -! And for the next iteration the trust region will evolves depending of -! the agreement of the energy prediction based on the Taylor series -! truncated at the 2nd order and the real energy. If the Taylor series -! truncated at the 2nd order represents correctly the energy landscape -! the trust region will be extent else it will be reduced. In order to -! mesure this agreement we use the ratio rho cf. "rho_model" and -! "trust_e_model". From that we use the following values: -! - if $\rho \geq 0.75$, then $\Delta = 2 \Delta$, -! - if $0.5 \geq \rho < 0.75$, then $\Delta = \Delta$, -! - if $0.25 \geq \rho < 0.5$, then $\Delta = 0.5 \Delta$, -! - if $\rho < 0.25$, then $\Delta = 0.25 \Delta$. - -! In addition, if $\rho < 0.1$ the iteration is cancelled, so it -! restarts with a smaller trust region until the energy decreases. - - - - -! Summary - -! To summarize, knowing the hessian (eigenvectors and eigenvalues), the -! gradient and the radius of the trust region we can compute the norm of -! the Newton step -! \begin{align*} -! ||\textbf{x}(\lambda = 0)||^2 = ||- \textbf{H}^{-1} \cdot \textbf{g}||^2 = \sum_{i=1}^n -! \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2}, \quad h_i \neq 0 -! \end{align*} - -! - if $h_1 \geq 0$, $||\textbf{x}(\lambda = 0)|| \leq \Delta$ and -! $\textbf{x}(\lambda=0)$ is in the trust region and it is not -! necessary to put a constraint on $\textbf{x}$, the solution is the -! unconstrained one, $\textbf{x}^* = \textbf{x}(\lambda = 0)$. -! - else if $h_1 < 0$, $\textbf{w}_1^T \cdot \textbf{g} = 0$ and -! $||\textbf{x}(\lambda = -h_1)|| \leq \Delta$ (by removing $j=1$ in -! the sum), the solution is $\textbf{x}^* = \textbf{x}(\lambda = -! -h_1)$, similarly to the previous case. -! But we can add to $\textbf{x}$, the first eigenvector $\textbf{W}_1$ -! time a constant to ensure the condition $||\textbf{x}(\lambda = -! -h_1)|| = \Delta$ and escape from the saddle point -! - else if $h_1 < 0$ and $\textbf{w}_1^T \cdot \textbf{g} \neq 0$ we -! have to search $\lambda \in (-h_1, \infty)$ such as -! $\textbf{x}(\lambda) = \Delta$ by solving with the Newton method -! \begin{align*} -! (||\textbf{x}(\lambda)||^2 - \Delta^2)^2 = 0 -! \end{align*} -! or -! \begin{align*} -! (\frac{1}{||\textbf{x}(\lambda)||^2} - \frac{1}{\Delta^2})^2 = 0 -! \end{align*} -! which is numerically more stable. And finally compute -! \begin{align*} -! \textbf{x}^* = \textbf{x}(\lambda) = - \sum_{i=1}^n \frac{\textbf{w}_i^T \cdot -! \textbf{g}}{h_i + \lambda} \cdot \textbf{w}_i -! \end{align*} -! - else if $h_1 \geq 0$ and $||\textbf{x}(\lambda = 0)|| > \Delta$ we -! do exactly the same thing that the previous case but we search -! $\lambda \in (0, \infty)$ -! - else if $h_1 < 0$ and $\textbf{w}_1^T \cdot \textbf{g} = 0$ and -! $||\textbf{x}(\lambda = -h_1)|| > \Delta$ (by removing $j=1$ in the -! sum), again we do exactly the same thing that the previous case -! searching $\lambda \in (-h_1, \infty)$. - - -! For the cases where $\textbf{w}_1^T \cdot \textbf{g} = 0$ it is not -! necessary in fact to remove the $j = 1$ in the sum since the term -! where $h_i - \lambda < 10^{-6}$ are not computed. - -! After that, we take this vector $\textbf{x}^*$, called "x", and we do -! the transformation to an antisymmetric matrix $\textbf{X}$, called -! m_x. This matrix $\textbf{X}$ will be used to compute a rotation -! matrix $\textbf{R}= \exp(\textbf{X})$ in "rotation_matrix". - -! NB: -! An improvement can be done using a elleptical trust region. - - - - -! Code - -! Provided: -! | mo_num | integer | number of MOs | - -! Cf. qp_edit in orbital optimization section, for some constants/thresholds - -! Input: -! | m | integer | number of MOs | -! | n | integer | m*(m-1)/2 | -! | H(n, n) | double precision | hessian | -! | v_grad(n) | double precision | gradient | -! | e_val(n) | double precision | eigenvalues of the hessian | -! | W(n, n) | double precision | eigenvectors of the hessian | -! | rho | double precision | agreement between the model and the reality, | -! | | | represents the quality of the energy prediction | -! | nb_iter | integer | number of iteration | - -! Input/Ouput: -! | delta | double precision | radius of the trust region | - -! Output: -! | x(n) | double precision | vector containing the step | - -! Internal: -! | accu | double precision | temporary variable to compute the step | -! | lambda | double precision | lagrange multiplier | -! | trust_radius2 | double precision | square of the radius of the trust region | -! | norm2_x | double precision | norm^2 of the vector x | -! | norm2_g | double precision | norm^2 of the vector containing the gradient | -! | tmp_wtg(n) | double precision | tmp_wtg(i) = w_i^T . g | -! | i, j, k | integer | indexes | - -! Function: -! | dnrm2 | double precision | Blas function computing the norm | -! | f_norm_trust_region_omp | double precision | compute the value of norm(x(lambda)^2) | - - -subroutine trust_region_step(n,nb_iter,v_grad,rho,e_val,w,x,delta) - - include 'pi.h' - - BEGIN_DOC - ! Compuet the step in the trust region - END_DOC - - implicit none - - ! Variables - - ! in - integer, intent(in) :: n - double precision, intent(in) :: v_grad(n), rho - integer, intent(inout) :: nb_iter - double precision, intent(in) :: e_val(n), w(n,n) - - ! inout - double precision, intent(inout) :: delta - - ! out - double precision, intent(out) :: x(n) - - ! Internal - double precision :: accu, lambda, trust_radius2 - double precision :: norm2_x, norm2_g - double precision, allocatable :: tmp_wtg(:) - integer :: i,j,k - double precision :: t1,t2,t3 - integer :: n_neg_eval - - - ! Functions - double precision :: ddot, dnrm2 - double precision :: f_norm_trust_region_omp - - print*,'' - print*,'==================' - print*,'---Trust_region---' - print*,'==================' - - call wall_time(t1) - - ! Allocation - allocate(tmp_wtg(n)) - -! Initialization and norm - -! The norm of the step size will be useful for the trust region -! algorithm. We start from a first guess and the radius of the trust -! region will evolve during the optimization. - -! avoid_saddle is actually a test to avoid saddle points - - -! Initialization of the Lagrange multiplier -lambda = 0d0 - -! List of w^T.g, to avoid the recomputation -tmp_wtg = 0d0 -do j = 1, n - do i = 1, n - tmp_wtg(j) = tmp_wtg(j) + w(i,j) * v_grad(i) - enddo -enddo - -! Replacement of the small tmp_wtg corresponding to a negative eigenvalue -! in the case of avoid_saddle -if (avoid_saddle .and. e_val(1) < - thresh_eig) then - i = 2 - ! Number of negative eigenvalues - do while (e_val(i) < - thresh_eig) - if (tmp_wtg(i) < thresh_wtg2) then - if (version_avoid_saddle == 1) then - tmp_wtg(i) = 1d0 - elseif (version_avoid_saddle == 2) then - tmp_wtg(i) = DABS(e_val(i)) - elseif (version_avoid_saddle == 3) then - tmp_wtg(i) = dsqrt(DABS(e_val(i))) - else - tmp_wtg(i) = thresh_wtg2 - endif - endif - i = i + 1 - enddo - - ! For the fist one it's a little bit different - if (tmp_wtg(1) < thresh_wtg2) then - tmp_wtg(1) = 0d0 - endif - -endif - -! Norm^2 of x, ||x||^2 -norm2_x = f_norm_trust_region_omp(n,e_val,tmp_wtg,0d0) -! We just use this norm for the nb_iter = 0 in order to initialize the trust radius delta -! We don't care about the sign of the eigenvalue we just want the size of the step in a normal Newton-Raphson algorithm -! Anyway if the step is too big it will be reduced -print*,'||x||^2 :', norm2_x - -! Norm^2 of the gradient, ||v_grad||^2 -norm2_g = (dnrm2(n,v_grad,1))**2 -print*,'||grad||^2 :', norm2_g - -! Trust radius initialization - -! At the first iteration (nb_iter = 0) we initialize the trust region -! with the norm of the step generate by the Newton's method ($\textbf{x}_1 = -! (\textbf{H}_0)^{-1} \cdot \textbf{g}_0$, -! we compute this norm using f_norm_trust_region_omp as explain just -! below) - - -! trust radius -if (nb_iter == 0) then - trust_radius2 = norm2_x - ! To avoid infinite loop of cancellation of this first step - ! without changing delta - nb_iter = 1 - - ! Compute delta, delta = sqrt(trust_radius) - delta = dsqrt(trust_radius2) -endif - -! Modification of the trust radius - -! In function of rho (which represents the agreement between the model -! and the reality, cf. rho_model) the trust region evolves. We update -! delta (the radius of the trust region). - -! To avoid too big trust region we put a maximum size. - - -! Modification of the trust radius in function of rho -if (rho >= 0.75d0) then - delta = 2d0 * delta -elseif (rho >= 0.5d0) then - delta = delta -elseif (rho >= 0.25d0) then - delta = 0.5d0 * delta -else - delta = 0.25d0 * delta -endif - -! Maximum size of the trust region -!if (delta > 0.5d0 * n * pi) then -! delta = 0.5d0 * n * pi -! print*,'Delta > delta_max, delta = 0.5d0 * n * pi' -!endif - -if (delta > 1d10) then - delta = 1d10 -endif - -print*, 'Delta :', delta - -! Calculation of the optimal lambda - -! We search the solution of $(||x||^2 - \Delta^2)^2 = 0$ -! - If $||\textbf{x}|| > \Delta$ or $h_1 < 0$ we have to add a constant -! $\lambda > 0 \quad \text{and} \quad \lambda > -h_1$ -! - If $||\textbf{x}|| \leq \Delta$ and $h_1 \geq 0$ the solution is the -! unconstrained one, $\lambda = 0$ - -! You will find more details at the beginning - - -! By giving delta, we search (||x||^2 - delta^2)^2 = 0 -! and not (||x||^2 - delta)^2 = 0 - -! Research of lambda to solve ||x(lambda)|| = Delta - -! Display -print*, 'e_val(1) = ', e_val(1) -print*, 'w_1^T.g =', tmp_wtg(1) - -! H positive definite -if (e_val(1) > - thresh_eig) then - norm2_x = f_norm_trust_region_omp(n,e_val,tmp_wtg,0d0) - print*, '||x(0)||=', dsqrt(norm2_x) - print*, 'Delta=', delta - - ! H positive definite, ||x(lambda = 0)|| <= Delta - if (dsqrt(norm2_x) <= delta) then - print*, 'H positive definite, ||x(lambda = 0)|| <= Delta' - print*, 'lambda = 0, no lambda optimization' - lambda = 0d0 - - ! H positive definite, ||x(lambda = 0)|| > Delta - else - ! Constraint solution - print*, 'H positive definite, ||x(lambda = 0)|| > Delta' - print*,'Computation of the optimal lambda...' - call trust_region_optimal_lambda(n,e_val,tmp_wtg,delta,lambda) - endif - -! H indefinite -else - if (DABS(tmp_wtg(1)) < thresh_wtg) then - norm2_x = f_norm_trust_region_omp(n,e_val,tmp_wtg, - e_val(1)) - print*, 'w_1^T.g <', thresh_wtg,', ||x(lambda = -e_val(1))|| =', dsqrt(norm2_x) - endif - - ! H indefinite, w_1^T.g = 0, ||x(lambda = -e_val(1))|| <= Delta - if (dsqrt(norm2_x) <= delta .and. DABS(tmp_wtg(1)) < thresh_wtg) then - ! Add e_val(1) in order to have (H - e_val(1) I) positive definite - print*, 'H indefinite, w_1^T.g = 0, ||x(lambda = -e_val(1))|| <= Delta' - print*, 'lambda = -e_val(1), no lambda optimization' - lambda = - e_val(1) - - ! H indefinite, w_1^T.g = 0, ||x(lambda = -e_val(1))|| > Delta - ! and - ! H indefinite, w_1^T.g =/= 0 - else - ! Constraint solution/ add lambda - if (DABS(tmp_wtg(1)) < thresh_wtg) then - print*, 'H indefinite, w_1^T.g = 0, ||x(lambda = -e_val(1))|| > Delta' - else - print*, 'H indefinite, w_1^T.g =/= 0' - endif - print*, 'Computation of the optimal lambda...' - call trust_region_optimal_lambda(n,e_val,tmp_wtg,delta,lambda) - endif - -endif - -! Recomputation of the norm^2 of the step x -norm2_x = f_norm_trust_region_omp(n,e_val,tmp_wtg,lambda) -print*,'' -print*,'Summary after the trust region:' -print*,'lambda:', lambda -print*,'||x||:', dsqrt(norm2_x) -print*,'delta:', delta - -! Calculation of the step x - -! x refers to $\textbf{x}^*$ -! We compute x in function of lambda using its formula : -! \begin{align*} -! \textbf{x}^* = \textbf{x}(\lambda) = - \sum_{i=1}^n \frac{\textbf{w}_i^T \cdot \textbf{g}}{h_i -! + \lambda} \cdot \textbf{w}_i -! \end{align*} - - -! Initialisation -x = 0d0 - -! Calculation of the step x - -! Normal version -if (.not. absolute_eig) then - - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - do j = 1, n - x(j) = x(j) - tmp_wtg(i) * W(j,i) / (e_val(i) + lambda) - enddo - endif - enddo - -! Version to use the absolute value of the eigenvalues -else - - do i = 1, n - if (DABS(e_val(i)) > thresh_eig) then - do j = 1, n - x(j) = x(j) - tmp_wtg(i) * W(j,i) / (DABS(e_val(i)) + lambda) - enddo - endif - enddo - -endif - -double precision :: beta, norm_x - -! Test -! If w_1^T.g = 0, the lim of ||x(lambda)|| when lambda tend to -e_val(1) -! is not + infinity. So ||x(lambda=-e_val(1))|| < delta, we add the first -! eigenvectors multiply by a constant to ensure the condition -! ||x(lambda=-e_val(1))|| = delta and escape the saddle point -if (avoid_saddle .and. e_val(1) < - thresh_eig) then - if (tmp_wtg(1) < 1d-15 .and. (1d0 - dsqrt(norm2_x)/delta) > 1d-3 ) then - - ! norm of x - norm_x = dnrm2(n,x,1) - - ! Computes the coefficient for the w_1 - beta = delta**2 - norm_x**2 - - ! Updates the step x - x = x + W(:,1) * dsqrt(beta) - - ! Recomputes the norm to check - norm_x = dnrm2(n,x,1) - - print*, 'Add w_1 * dsqrt(delta^2 - ||x||^2):' - print*, '||x||', norm_x - endif -endif - -! Transformation of x - -! x is a vector of size n, so it can be write as a m by m -! antisymmetric matrix m_x cf. "mat_to_vec_index" and "vec_to_mat_index". - - -! ! Step transformation vector -> matrix -! ! Vector with n element -> mo_num by mo_num matrix -! do j = 1, m -! do i = 1, m -! if (i>j) then -! call mat_to_vec_index(i,j,k) -! m_x(i,j) = x(k) -! else -! m_x(i,j) = 0d0 -! endif -! enddo -! enddo -! -! ! Antisymmetrization of the previous matrix -! do j = 1, m -! do i = 1, m -! if (i 0$ ($\lambda = 0$ is the unconstraint solution). But the -Newton method is only defined for a positive definite hessian matrix, -so $(\textbf{H} + \textbf{I} \lambda)$ must be positive -definite. Consequently, in the case where $\textbf{H}$ is not positive -definite, to ensure the positive definiteness, $\lambda$ must be -greater than $- h_1$. -\begin{align*} - \lambda > 0 \quad \text{and} \quad \lambda \geq - h_1 -\end{align*} - -From that there are five cases: -- if $\textbf{H}$ is positive definite, $-h_1 < 0$, $\lambda \in (0,\infty)$ -- if $\textbf{H}$ is not positive definite and $\textbf{w}_1^T \cdot - \textbf{g} \neq 0$, $(\textbf{H} + \textbf{I} - \lambda)$ - must be positve definite, $-h_1 > 0$, $\lambda \in (-h_1, \infty)$ -- if $\textbf{H}$ is not positive definite , $\textbf{w}_1^T \cdot - \textbf{g} = 0$ and $||\textbf{x}(-h_1)|| > \Delta$ by removing - $j=1$ in the sum, $(\textbf{H} + \textbf{I} \lambda)$ must be - positive definite, $-h_1 > 0$, $\lambda \in (-h_1, \infty$) -- if $\textbf{H}$ is not positive definite , $\textbf{w}_1^T \cdot - \textbf{g} = 0$ and $||\textbf{x}(-h_1)|| \leq \Delta$ by removing - $j=1$ in the sum, $(\textbf{H} + \textbf{I} \lambda)$ must be - positive definite, $-h_1 > 0$, $\lambda = -h_1$). This case is - similar to the case where $\textbf{H}$ and $||\textbf{x}(\lambda = - 0)|| \leq \Delta$ - but we can also add to $\textbf{x}$, the first eigenvector $\textbf{W}_1$ - time a constant to ensure the condition $||\textbf{x}(\lambda = - -h_1)|| = \Delta$ and escape from the saddle point - -Thus to find the solution, we can write: -\begin{align*} - ||\textbf{x}(\lambda)|| = \Delta -\end{align*} -\begin{align*} - ||\textbf{x}(\lambda)|| - \Delta = 0 -\end{align*} - -Taking the square of this equation -\begin{align*} - (||\textbf{x}(\lambda)|| - \Delta)^2 = 0 -\end{align*} -we have a function with one minimum for the optimal $\lambda$. -Since we have the formula of $||\textbf{x}(\lambda)||^2$, we solve -\begin{align*} - (||\textbf{x}(\lambda)||^2 - \Delta^2)^2 = 0 -\end{align*} - -But in practice, it is more effective to solve: -\begin{align*} - (\frac{1}{||\textbf{x}(\lambda)||^2} - \frac{1}{\Delta^2})^2 = 0 -\end{align*} - -To do that, we just use the Newton method with "trust_newton" using -first and second derivative of $(||\textbf{x}(\lambda)||^2 - -\Delta^2)^2$ with respect to $\textbf{x}$. -This will give the optimal $\lambda$ to compute the -solution $\textbf{x}$ with the formula seen previously: -\begin{align*} - \textbf{x}(\lambda) = - \sum_{i=1}^n \frac{\textbf{w}_i^T \cdot - \textbf{g}}{h_i + \lambda} \cdot \textbf{w}_i -\end{align*} - -The solution $\textbf{x}(\lambda)$ with the optimal $\lambda$ is our -step to go from the (k)-th to the (k+1)-th iteration, is noted $\textbf{x}^*$. - -#+BEGIN_SRC f90 :comments org :tangle trust_region_step.irp.f -#+END_SRC - -** Evolution of the trust region - -We initialize the trust region at the first iteration using a radius -\begin{align*} - \Delta = ||\textbf{x}(\lambda=0)|| -\end{align*} - -And for the next iteration the trust region will evolves depending of -the agreement of the energy prediction based on the Taylor series -truncated at the 2nd order and the real energy. If the Taylor series -truncated at the 2nd order represents correctly the energy landscape -the trust region will be extent else it will be reduced. In order to -mesure this agreement we use the ratio rho cf. "rho_model" and -"trust_e_model". From that we use the following values: -- if $\rho \geq 0.75$, then $\Delta = 2 \Delta$, -- if $0.5 \geq \rho < 0.75$, then $\Delta = \Delta$, -- if $0.25 \geq \rho < 0.5$, then $\Delta = 0.5 \Delta$, -- if $\rho < 0.25$, then $\Delta = 0.25 \Delta$. - -In addition, if $\rho < 0.1$ the iteration is cancelled, so it -restarts with a smaller trust region until the energy decreases. - -#+BEGIN_SRC f90 :comments org :tangle trust_region_step.irp.f -#+END_SRC - -** Summary - -To summarize, knowing the hessian (eigenvectors and eigenvalues), the -gradient and the radius of the trust region we can compute the norm of -the Newton step -\begin{align*} - ||\textbf{x}(\lambda = 0)||^2 = ||- \textbf{H}^{-1} \cdot \textbf{g}||^2 = \sum_{i=1}^n - \frac{(\textbf{w}_i^T \cdot \textbf{g})^2}{(h_i + \lambda)^2}, \quad h_i \neq 0 -\end{align*} - -- if $h_1 \geq 0$, $||\textbf{x}(\lambda = 0)|| \leq \Delta$ and - $\textbf{x}(\lambda=0)$ is in the trust region and it is not - necessary to put a constraint on $\textbf{x}$, the solution is the - unconstrained one, $\textbf{x}^* = \textbf{x}(\lambda = 0)$. -- else if $h_1 < 0$, $\textbf{w}_1^T \cdot \textbf{g} = 0$ and - $||\textbf{x}(\lambda = -h_1)|| \leq \Delta$ (by removing $j=1$ in - the sum), the solution is $\textbf{x}^* = \textbf{x}(\lambda = - -h_1)$, similarly to the previous case. - But we can add to $\textbf{x}$, the first eigenvector $\textbf{W}_1$ - time a constant to ensure the condition $||\textbf{x}(\lambda = - -h_1)|| = \Delta$ and escape from the saddle point -- else if $h_1 < 0$ and $\textbf{w}_1^T \cdot \textbf{g} \neq 0$ we - have to search $\lambda \in (-h_1, \infty)$ such as - $\textbf{x}(\lambda) = \Delta$ by solving with the Newton method - \begin{align*} - (||\textbf{x}(\lambda)||^2 - \Delta^2)^2 = 0 - \end{align*} - or - \begin{align*} - (\frac{1}{||\textbf{x}(\lambda)||^2} - \frac{1}{\Delta^2})^2 = 0 - \end{align*} - which is numerically more stable. And finally compute - \begin{align*} - \textbf{x}^* = \textbf{x}(\lambda) = - \sum_{i=1}^n \frac{\textbf{w}_i^T \cdot - \textbf{g}}{h_i + \lambda} \cdot \textbf{w}_i - \end{align*} -- else if $h_1 \geq 0$ and $||\textbf{x}(\lambda = 0)|| > \Delta$ we - do exactly the same thing that the previous case but we search - $\lambda \in (0, \infty)$ -- else if $h_1 < 0$ and $\textbf{w}_1^T \cdot \textbf{g} = 0$ and - $||\textbf{x}(\lambda = -h_1)|| > \Delta$ (by removing $j=1$ in the - sum), again we do exactly the same thing that the previous case - searching $\lambda \in (-h_1, \infty)$. - - -For the cases where $\textbf{w}_1^T \cdot \textbf{g} = 0$ it is not -necessary in fact to remove the $j = 1$ in the sum since the term -where $h_i - \lambda < 10^{-6}$ are not computed. - -After that, we take this vector $\textbf{x}^*$, called "x", and we do -the transformation to an antisymmetric matrix $\textbf{X}$, called -m_x. This matrix $\textbf{X}$ will be used to compute a rotation -matrix $\textbf{R}= \exp(\textbf{X})$ in "rotation_matrix". - -NB: -An improvement can be done using a elleptical trust region. - -#+BEGIN_SRC f90 :comments org :tangle trust_region_step.irp.f -#+END_SRC - -** Code - -Provided: -| mo_num | integer | number of MOs | - -Cf. qp_edit in orbital optimization section, for some constants/thresholds - -Input: -| m | integer | number of MOs | -| n | integer | m*(m-1)/2 | -| H(n, n) | double precision | hessian | -| v_grad(n) | double precision | gradient | -| e_val(n) | double precision | eigenvalues of the hessian | -| W(n, n) | double precision | eigenvectors of the hessian | -| rho | double precision | agreement between the model and the reality, | -| | | represents the quality of the energy prediction | -| nb_iter | integer | number of iteration | - -Input/Ouput: -| delta | double precision | radius of the trust region | - -Output: -| x(n) | double precision | vector containing the step | - -Internal: -| accu | double precision | temporary variable to compute the step | -| lambda | double precision | lagrange multiplier | -| trust_radius2 | double precision | square of the radius of the trust region | -| norm2_x | double precision | norm^2 of the vector x | -| norm2_g | double precision | norm^2 of the vector containing the gradient | -| tmp_wtg(n) | double precision | tmp_wtg(i) = w_i^T . g | -| i, j, k | integer | indexes | - -Function: -| dnrm2 | double precision | Blas function computing the norm | -| f_norm_trust_region_omp | double precision | compute the value of norm(x(lambda)^2) | - -#+BEGIN_SRC f90 :comments org :tangle trust_region_step.irp.f -subroutine trust_region_step(n,nb_iter,v_grad,rho,e_val,w,x,delta) - - include 'pi.h' - - BEGIN_DOC - ! Compuet the step in the trust region - END_DOC - - implicit none - - ! Variables - - ! in - integer, intent(in) :: n - double precision, intent(in) :: v_grad(n), rho - integer, intent(inout) :: nb_iter - double precision, intent(in) :: e_val(n), w(n,n) - - ! inout - double precision, intent(inout) :: delta - - ! out - double precision, intent(out) :: x(n) - - ! Internal - double precision :: accu, lambda, trust_radius2 - double precision :: norm2_x, norm2_g - double precision, allocatable :: tmp_wtg(:) - integer :: i,j,k - double precision :: t1,t2,t3 - integer :: n_neg_eval - - - ! Functions - double precision :: ddot, dnrm2 - double precision :: f_norm_trust_region_omp - - print*,'' - print*,'==================' - print*,'---Trust_region---' - print*,'==================' - - call wall_time(t1) - - ! Allocation - allocate(tmp_wtg(n)) -#+END_SRC - - -*** Initialization and norm - -The norm of the step size will be useful for the trust region -algorithm. We start from a first guess and the radius of the trust -region will evolve during the optimization. - -avoid_saddle is actually a test to avoid saddle points - -#+BEGIN_SRC f90 :comments org :tangle trust_region_step.irp.f - ! Initialization of the Lagrange multiplier - lambda = 0d0 - - ! List of w^T.g, to avoid the recomputation - tmp_wtg = 0d0 - do j = 1, n - do i = 1, n - tmp_wtg(j) = tmp_wtg(j) + w(i,j) * v_grad(i) - enddo - enddo - - ! Replacement of the small tmp_wtg corresponding to a negative eigenvalue - ! in the case of avoid_saddle - if (avoid_saddle .and. e_val(1) < - thresh_eig) then - i = 2 - ! Number of negative eigenvalues - do while (e_val(i) < - thresh_eig) - if (tmp_wtg(i) < thresh_wtg2) then - if (version_avoid_saddle == 1) then - tmp_wtg(i) = 1d0 - elseif (version_avoid_saddle == 2) then - tmp_wtg(i) = DABS(e_val(i)) - elseif (version_avoid_saddle == 3) then - tmp_wtg(i) = dsqrt(DABS(e_val(i))) - else - tmp_wtg(i) = thresh_wtg2 - endif - endif - i = i + 1 - enddo - - ! For the fist one it's a little bit different - if (tmp_wtg(1) < thresh_wtg2) then - tmp_wtg(1) = 0d0 - endif - - endif - - ! Norm^2 of x, ||x||^2 - norm2_x = f_norm_trust_region_omp(n,e_val,tmp_wtg,0d0) - ! We just use this norm for the nb_iter = 0 in order to initialize the trust radius delta - ! We don't care about the sign of the eigenvalue we just want the size of the step in a normal Newton-Raphson algorithm - ! Anyway if the step is too big it will be reduced - print*,'||x||^2 :', norm2_x - - ! Norm^2 of the gradient, ||v_grad||^2 - norm2_g = (dnrm2(n,v_grad,1))**2 - print*,'||grad||^2 :', norm2_g -#+END_SRC - -*** Trust radius initialization - - At the first iteration (nb_iter = 0) we initialize the trust region - with the norm of the step generate by the Newton's method ($\textbf{x}_1 = - (\textbf{H}_0)^{-1} \cdot \textbf{g}_0$, - we compute this norm using f_norm_trust_region_omp as explain just - below) - -#+BEGIN_SRC f90 :comments org :tangle trust_region_step.irp.f - ! trust radius - if (nb_iter == 0) then - trust_radius2 = norm2_x - ! To avoid infinite loop of cancellation of this first step - ! without changing delta - nb_iter = 1 - - ! Compute delta, delta = sqrt(trust_radius) - delta = dsqrt(trust_radius2) - endif -#+END_SRC - -*** Modification of the trust radius - -In function of rho (which represents the agreement between the model -and the reality, cf. rho_model) the trust region evolves. We update -delta (the radius of the trust region). - -To avoid too big trust region we put a maximum size. - -#+BEGIN_SRC f90 :comments org :tangle trust_region_step.irp.f - ! Modification of the trust radius in function of rho - if (rho >= 0.75d0) then - delta = 2d0 * delta - elseif (rho >= 0.5d0) then - delta = delta - elseif (rho >= 0.25d0) then - delta = 0.5d0 * delta - else - delta = 0.25d0 * delta - endif - - ! Maximum size of the trust region - !if (delta > 0.5d0 * n * pi) then - ! delta = 0.5d0 * n * pi - ! print*,'Delta > delta_max, delta = 0.5d0 * n * pi' - !endif - - if (delta > 1d10) then - delta = 1d10 - endif - - print*, 'Delta :', delta -#+END_SRC - -*** Calculation of the optimal lambda - -We search the solution of $(||x||^2 - \Delta^2)^2 = 0$ -- If $||\textbf{x}|| > \Delta$ or $h_1 < 0$ we have to add a constant - $\lambda > 0 \quad \text{and} \quad \lambda > -h_1$ -- If $||\textbf{x}|| \leq \Delta$ and $h_1 \geq 0$ the solution is the - unconstrained one, $\lambda = 0$ - -You will find more details at the beginning - -#+BEGIN_SRC f90 :comments org :tangle trust_region_step.irp.f - ! By giving delta, we search (||x||^2 - delta^2)^2 = 0 - ! and not (||x||^2 - delta)^2 = 0 - - ! Research of lambda to solve ||x(lambda)|| = Delta - - ! Display - print*, 'e_val(1) = ', e_val(1) - print*, 'w_1^T.g =', tmp_wtg(1) - - ! H positive definite - if (e_val(1) > - thresh_eig) then - norm2_x = f_norm_trust_region_omp(n,e_val,tmp_wtg,0d0) - print*, '||x(0)||=', dsqrt(norm2_x) - print*, 'Delta=', delta - - ! H positive definite, ||x(lambda = 0)|| <= Delta - if (dsqrt(norm2_x) <= delta) then - print*, 'H positive definite, ||x(lambda = 0)|| <= Delta' - print*, 'lambda = 0, no lambda optimization' - lambda = 0d0 - - ! H positive definite, ||x(lambda = 0)|| > Delta - else - ! Constraint solution - print*, 'H positive definite, ||x(lambda = 0)|| > Delta' - print*,'Computation of the optimal lambda...' - call trust_region_optimal_lambda(n,e_val,tmp_wtg,delta,lambda) - endif - - ! H indefinite - else - if (DABS(tmp_wtg(1)) < thresh_wtg) then - norm2_x = f_norm_trust_region_omp(n,e_val,tmp_wtg, - e_val(1)) - print*, 'w_1^T.g <', thresh_wtg,', ||x(lambda = -e_val(1))|| =', dsqrt(norm2_x) - endif - - ! H indefinite, w_1^T.g = 0, ||x(lambda = -e_val(1))|| <= Delta - if (dsqrt(norm2_x) <= delta .and. DABS(tmp_wtg(1)) < thresh_wtg) then - ! Add e_val(1) in order to have (H - e_val(1) I) positive definite - print*, 'H indefinite, w_1^T.g = 0, ||x(lambda = -e_val(1))|| <= Delta' - print*, 'lambda = -e_val(1), no lambda optimization' - lambda = - e_val(1) - - ! H indefinite, w_1^T.g = 0, ||x(lambda = -e_val(1))|| > Delta - ! and - ! H indefinite, w_1^T.g =/= 0 - else - ! Constraint solution/ add lambda - if (DABS(tmp_wtg(1)) < thresh_wtg) then - print*, 'H indefinite, w_1^T.g = 0, ||x(lambda = -e_val(1))|| > Delta' - else - print*, 'H indefinite, w_1^T.g =/= 0' - endif - print*, 'Computation of the optimal lambda...' - call trust_region_optimal_lambda(n,e_val,tmp_wtg,delta,lambda) - endif - - endif - - ! Recomputation of the norm^2 of the step x - norm2_x = f_norm_trust_region_omp(n,e_val,tmp_wtg,lambda) - print*,'' - print*,'Summary after the trust region:' - print*,'lambda:', lambda - print*,'||x||:', dsqrt(norm2_x) - print*,'delta:', delta -#+END_SRC - -*** Calculation of the step x - -x refers to $\textbf{x}^*$ -We compute x in function of lambda using its formula : -\begin{align*} -\textbf{x}^* = \textbf{x}(\lambda) = - \sum_{i=1}^n \frac{\textbf{w}_i^T \cdot \textbf{g}}{h_i -+ \lambda} \cdot \textbf{w}_i -\end{align*} - -#+BEGIN_SRC f90 :comments org :tangle trust_region_step.irp.f - ! Initialisation - x = 0d0 - - ! Calculation of the step x - - ! Normal version - if (.not. absolute_eig) then - - do i = 1, n - if (DABS(e_val(i)) > thresh_eig .and. DABS(e_val(i)+lambda) > thresh_eig) then - do j = 1, n - x(j) = x(j) - tmp_wtg(i) * W(j,i) / (e_val(i) + lambda) - enddo - endif - enddo - - ! Version to use the absolute value of the eigenvalues - else - - do i = 1, n - if (DABS(e_val(i)) > thresh_eig) then - do j = 1, n - x(j) = x(j) - tmp_wtg(i) * W(j,i) / (DABS(e_val(i)) + lambda) - enddo - endif - enddo - - endif - - double precision :: beta, norm_x - - ! Test - ! If w_1^T.g = 0, the lim of ||x(lambda)|| when lambda tend to -e_val(1) - ! is not + infinity. So ||x(lambda=-e_val(1))|| < delta, we add the first - ! eigenvectors multiply by a constant to ensure the condition - ! ||x(lambda=-e_val(1))|| = delta and escape the saddle point - if (avoid_saddle .and. e_val(1) < - thresh_eig) then - if (tmp_wtg(1) < 1d-15 .and. (1d0 - dsqrt(norm2_x)/delta) > 1d-3 ) then - - ! norm of x - norm_x = dnrm2(n,x,1) - - ! Computes the coefficient for the w_1 - beta = delta**2 - norm_x**2 - - ! Updates the step x - x = x + W(:,1) * dsqrt(beta) - - ! Recomputes the norm to check - norm_x = dnrm2(n,x,1) - - print*, 'Add w_1 * dsqrt(delta^2 - ||x||^2):' - print*, '||x||', norm_x - endif - endif -#+END_SRC - -*** Transformation of x - -x is a vector of size n, so it can be write as a m by m -antisymmetric matrix m_x cf. "mat_to_vec_index" and "vec_to_mat_index". - - #+BEGIN_SRC f90 :comments org :tangle trust_region_step.irp.f -! ! Step transformation vector -> matrix -! ! Vector with n element -> mo_num by mo_num matrix -! do j = 1, m -! do i = 1, m -! if (i>j) then -! call mat_to_vec_index(i,j,k) -! m_x(i,j) = x(k) -! else -! m_x(i,j) = 0d0 -! endif -! enddo -! enddo -! -! ! Antisymmetrization of the previous matrix -! do j = 1, m -! do i = 1, m -! if (i lower diagonal matrix (p,q), p > q - -! If a matrix is antisymmetric it can be reshaped as a vector. And the -! vector can be reshaped as an antisymmetric matrix - -! \begin{align*} -! \begin{pmatrix} -! 0 & -1 & -2 & -4 \\ -! 1 & 0 & -3 & -5 \\ -! 2 & 3 & 0 & -6 \\ -! 4 & 5 & 6 & 0 -! \end{pmatrix} -! \Leftrightarrow -! \begin{pmatrix} -! 1 & 2 & 3 & 4 & 5 & 6 -! \end{pmatrix} -! \end{align*} - -! !!! Here the algorithm only work for the lower diagonal !!! - -! Input: -! | i | integer | index in the vector | - -! Ouput: -! | p,q | integer | corresponding indexes in the lower diagonal of a matrix | -! | | | p > q, | -! | | | p -> row, | -! | | | q -> column | - - -subroutine vec_to_mat_index(i,p,q) - - include 'pi.h' - - BEGIN_DOC - ! Compute the indexes (p,q) of the element in the lower diagonal matrix knowing - ! its index i a vector - END_DOC - - implicit none - - ! Variables - - ! in - integer,intent(in) :: i - - ! out - integer, intent(out) :: p,q - - ! internal - integer :: a,b - double precision :: da - - da = 0.5d0*(1+ sqrt(1d0+8d0*DBLE(i))) - a = INT(da) - if ((a*(a-1))/2==i) then - p = a-1 - else - p = a - endif - b = p*(p-1)/2 - - ! Matrix element indexes - p = p + 1 - q = i - b - -end subroutine diff --git a/src/utils_trust_region/vec_to_mat_index.org b/src/utils_trust_region/vec_to_mat_index.org deleted file mode 100644 index 0a09fa86..00000000 --- a/src/utils_trust_region/vec_to_mat_index.org +++ /dev/null @@ -1,72 +0,0 @@ -* Vector to matrix indexes - -*Compute the indexes p,q of a matrix element with the vector index i* - -Vector (i) -> lower diagonal matrix (p,q), p > q - -If a matrix is antisymmetric it can be reshaped as a vector. And the -vector can be reshaped as an antisymmetric matrix - -\begin{align*} -\begin{pmatrix} -0 & -1 & -2 & -4 \\ -1 & 0 & -3 & -5 \\ -2 & 3 & 0 & -6 \\ -4 & 5 & 6 & 0 -\end{pmatrix} -\Leftrightarrow -\begin{pmatrix} -1 & 2 & 3 & 4 & 5 & 6 -\end{pmatrix} -\end{align*} - -!!! Here the algorithm only work for the lower diagonal !!! - -Input: -| i | integer | index in the vector | - -Ouput: -| p,q | integer | corresponding indexes in the lower diagonal of a matrix | -| | | p > q, | -| | | p -> row, | -| | | q -> column | - -#+BEGIN_SRC f90 :comments org :tangle vec_to_mat_index.irp.f -subroutine vec_to_mat_index(i,p,q) - - include 'pi.h' - - BEGIN_DOC - ! Compute the indexes (p,q) of the element in the lower diagonal matrix knowing - ! its index i a vector - END_DOC - - implicit none - - ! Variables - - ! in - integer,intent(in) :: i - - ! out - integer, intent(out) :: p,q - - ! internal - integer :: a,b - double precision :: da - - da = 0.5d0*(1+ sqrt(1d0+8d0*DBLE(i))) - a = INT(da) - if ((a*(a-1))/2==i) then - p = a-1 - else - p = a - endif - b = p*(p-1)/2 - - ! Matrix element indexes - p = p + 1 - q = i - b - -end subroutine -#+END_SRC diff --git a/src/utils_trust_region/vec_to_mat_v2.irp.f b/src/utils_trust_region/vec_to_mat_v2.irp.f deleted file mode 100644 index 9140b8d3..00000000 --- a/src/utils_trust_region/vec_to_mat_v2.irp.f +++ /dev/null @@ -1,39 +0,0 @@ -! Vect to antisymmetric matrix using mat_to_vec_index - -! Vector to antisymmetric matrix transformation using mat_to_vec_index -! subroutine. - -! Can be done in OMP (for the first part and with omp critical for the second) - - -subroutine vec_to_mat_v2(n,m,v_x,m_x) - - BEGIN_DOC - ! Vector to antisymmetric matrix - END_DOC - - implicit none - - integer, intent(in) :: n,m - double precision, intent(in) :: v_x(n) - double precision, intent(out) :: m_x(m,m) - - integer :: i,j,k - - ! 1D -> 2D lower diagonal - m_x = 0d0 - do j = 1, m - 1 - do i = j + 1, m - call mat_to_vec_index(i,j,k) - m_x(i,j) = v_x(k) - enddo - enddo - - ! Antisym - do i = 1, m - 1 - do j = i + 1, m - m_x(i,j) = - m_x(j,i) - enddo - enddo - -end diff --git a/src/utils_trust_region/vec_to_mat_v2.org b/src/utils_trust_region/vec_to_mat_v2.org deleted file mode 100644 index 4e358a88..00000000 --- a/src/utils_trust_region/vec_to_mat_v2.org +++ /dev/null @@ -1,40 +0,0 @@ -* Vect to antisymmetric matrix using mat_to_vec_index - -Vector to antisymmetric matrix transformation using mat_to_vec_index -subroutine. - -Can be done in OMP (for the first part and with omp critical for the second) - -#+BEGIN_SRC f90 :comments org :tangle vec_to_mat_v2.irp.f -subroutine vec_to_mat_v2(n,m,v_x,m_x) - - BEGIN_DOC - ! Vector to antisymmetric matrix - END_DOC - - implicit none - - integer, intent(in) :: n,m - double precision, intent(in) :: v_x(n) - double precision, intent(out) :: m_x(m,m) - - integer :: i,j,k - - ! 1D -> 2D lower diagonal - m_x = 0d0 - do j = 1, m - 1 - do i = j + 1, m - call mat_to_vec_index(i,j,k) - m_x(i,j) = v_x(k) - enddo - enddo - - ! Antisym - do i = 1, m - 1 - do j = i + 1, m - m_x(i,j) = - m_x(j,i) - enddo - enddo - -end -#+END_SRC