Commit 2d2a41b7 authored by Yingzhou Li's avatar Yingzhou Li

Tried cubic line search

parent 42ce8d6b
......@@ -8,6 +8,7 @@ LIST(APPEND elsi_rci_src
src/Davidson/elsi_rci_davidson.f90
src/OMM/elsi_rci_omm.f90
src/OMM/rci_omm_quartic.f90
src/OMM/rci_omm_cubic.f90
src/PPCG/elsi_rci_ppcg.f90)
ADD_LIBRARY(elsi_rci ${elsi_rci_src})
......
......@@ -572,12 +572,13 @@ contains
end do
n = n + next
next = n_state - int(sum(Vec_conv(1:n)), i4)
Vec_evold = Vec_ev
if (r_h%verbose > 0) write(*,'(a,i5,a,es15.7e3)') &
'Iter ', niter, ', Rel Err: ', sum(abs((Vec_ev(1:n_state) &
- Vec_evold(1:n_state))/Vec_ev(1)))
Vec_evold = Vec_ev
call rci_op_null(task)
if (next == 0 .or. niter >= max_iter) then
ijob = SID_FINISH
......
subroutine rci_omm_solve_cubic(coeff, x_min)
use ELSI_RCI_PRECISION, only:r8
implicit none
!**** INPUT ***********************************!
real(r8), intent(in) :: coeff(0:3) ! coeffs. of the cubic equation
!**** OUTPUT **********************************!
real(r8), intent(out) :: x_min ! position of minimum
!**** LOCAL ***********************************!
real(r8) :: a, b, c, d
real(r8) :: Pi = 3.141592653589793238462643383279502884197_r8
real(r8) :: p
real(r8) :: q
real(r8) :: p3
real(r8) :: q2
real(r8) :: delta, qrtd, gamma
! ax^3+bx^2+cx+d = 0
a = coeff(3)
b = coeff(2)
c = coeff(1)
d = coeff(0)
p = (3.0_r8*c/a-(b/a)**2.0_r8)/3.0_r8;
q = (2.0_r8*(b/a)**3.0_r8-9.0_r8*b/a*c/a+27.0_r8*d/a)/27.0_r8
p3 = p/3.0_r8
q2 = q/2.0_r8
delta = p3**3.0_r8 +q2**2.0_r8
if (delta >= 0.0_r8) then
qrtd = sqrt(delta)
x_min = sign(abs(-q2+qrtd)**(1.0_r8/3.0_r8), -q2+qrtd) &
+ sign(abs(-q2-qrtd)**(1.0_r8/3.0_r8), -q2-qrtd)
else
qrtd = sqrt(-delta)
if (q2 >= 0.0_r8) then
x_min = 2.0_r8*sqrt(-p3) &
*cos((atan2(-qrtd,-q2)-2.0_r8*pi)/3.0_r8)
else
x_min = 2.0_r8*sqrt(-p3)*cos(atan2(qrtd,-q2)/3.0_r8)
end if
end if
x_min = x_min - b/a/3.0_r8
end subroutine rci_omm_solve_cubic
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment