LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-21 11:03:29 +01:00
Code Releases Activity

Added utils_trust_region directory
Some checks failed
continuous-integration/drone/push Build is failing
Details

Browse Source
This commit is contained in:
Anthony Scemama 2023-02-08 16:44:40 +01:00
parent f947412a16
commit d84092a29d
27 changed files with 7734 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

89
src/utils_trust_region/EZFIO.cfg Normal file
View File

@ -0,0 +1,89 @@
[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

1
src/utils_trust_region/NEED Normal file
View File

@ -0,0 +1 @@
hartree_fock

5
src/utils_trust_region/README.rst Normal file
View File

@ -0,0 +1,5 @@
============
trust_region
============
The documentation can be found in the org files.

7
src/utils_trust_region/TANGLE_org_mode.sh Executable file
View File

@ -0,0 +1,7 @@
#!/bin/sh
list='ls *.org'
for element in $list
do
emacs --batch $element -f org-babel-tangle
done

248
src/utils_trust_region/algo_trust.irp.f Normal file
View File

@ -0,0 +1,248 @@
! 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

593
src/utils_trust_region/algo_trust.org Normal file
View File

@ -0,0 +1,593 @@
* 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

85
src/utils_trust_region/apply_mo_rotation.irp.f Normal file
View File

@ -0,0 +1,85 @@
! 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

86
src/utils_trust_region/apply_mo_rotation.org Normal file
View File

@ -0,0 +1,86 @@
* 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

61
src/utils_trust_region/mat_to_vec_index.irp.f Normal file
View File

@ -0,0 +1,61 @@
! 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

63
src/utils_trust_region/mat_to_vec_index.org Normal file
View File

@ -0,0 +1,63 @@
* 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

2
src/utils_trust_region/pi.h Normal file
View File

@ -0,0 +1,2 @@
!logical, parameter :: debug=.False.
double precision, parameter :: pi = 3.1415926535897932d0

443
src/utils_trust_region/rotation_matrix.irp.f Normal file
View File

@ -0,0 +1,443 @@
! 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

454
src/utils_trust_region/rotation_matrix.org Normal file
View File

@ -0,0 +1,454 @@
* 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

64
src/utils_trust_region/sub_to_full_rotation_matrix.irp.f Normal file
View File

@ -0,0 +1,64 @@
! 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

65
src/utils_trust_region/sub_to_full_rotation_matrix.org Normal file
View File

@ -0,0 +1,65 @@
* 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

119
src/utils_trust_region/trust_region_expected_e.irp.f Normal file
View File

@ -0,0 +1,119 @@
! 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

121
src/utils_trust_region/trust_region_expected_e.org Normal file
View File

@ -0,0 +1,121 @@
* 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

1655
src/utils_trust_region/trust_region_optimal_lambda.irp.f Normal file
View File

File diff suppressed because it is too large Load Diff

1665
src/utils_trust_region/trust_region_optimal_lambda.org Normal file
View File

File diff suppressed because it is too large Load Diff

121
src/utils_trust_region/trust_region_rho.irp.f Normal file
View File

@ -0,0 +1,121 @@
! 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

123
src/utils_trust_region/trust_region_rho.org Normal file
View File

@ -0,0 +1,123 @@
* 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

716
src/utils_trust_region/trust_region_step.irp.f Normal file
View File

@ -0,0 +1,716 @@
! 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<j) then
! m_x(i,j) = - m_x(j,i)
! endif
! enddo
! enddo
! Deallocation, end
deallocate(tmp_wtg)
call wall_time(t2)
t3 = t2 - t1
print*,'Time in trust_region:', t3
print*,'======================'
print*,'---End trust_region---'
print*,'======================'
print*,''
end

726
src/utils_trust_region/trust_region_step.org Normal file
View File

@ -0,0 +1,726 @@
* 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}^*$.
#+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<j) then
! m_x(i,j) = - m_x(j,i)
! endif
! enddo
! enddo
#+END_SRC
*** Deallocation, end
#+BEGIN_SRC f90 :comments org :tangle trust_region_step.irp.f
deallocate(tmp_wtg)
call wall_time(t2)
t3 = t2 - t1
print*,'Time in trust_region:', t3
print*,'======================'
print*,'---End trust_region---'
print*,'======================'
print*,''
end
#+END_SRC

71
src/utils_trust_region/vec_to_mat_index.irp.f Normal file
View File

@ -0,0 +1,71 @@
! 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 |
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

72
src/utils_trust_region/vec_to_mat_index.org Normal file
View File

@ -0,0 +1,72 @@
* 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

39
src/utils_trust_region/vec_to_mat_v2.irp.f Normal file
View File

@ -0,0 +1,39 @@
! 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

40
src/utils_trust_region/vec_to_mat_v2.org Normal file
View File

@ -0,0 +1,40 @@
* 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