Commit c23265a7 authored by Yingzhou Li's avatar Yingzhou Li

Added csc, and omm cubic routine, also updated the convergence critia for PPCG

parent 52d540a1
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
......@@ -336,10 +336,12 @@ contains
! - check conv
if (ijob == SID_LOCK + 6) then
trG = resvec(1)
trdiff = abs(trG-trG_old)
trdiff = abs((trG-trG_old)/trG)
trG_old = trG
if (r_h%verbose > 0) write(*,'(a,i5,a,es15.7e3)') &
'Iter ', iter, ', Rel Err: ', trdiff/trG
if (r_h%verbose > 0 .and. iter > 0) then
write(*,'(a,i5,a,es15.7e3)') &
'Iter ', iter, ', Rel Err: ', trdiff
end if
if ((nact == 0) .or. (iter == max_iter)) then
ijob = SID_FINISH
......@@ -1329,13 +1331,12 @@ contains
! - Check conv
if (ijob == SID_NORR + 4) then
trG = resvec(1)
trdiff = abs(trG-trG_old)
trG_old = trG
trdiff = abs((trG-trG_old)/trG)
if (r_h%verbose > 0) write(*,'(a,i5,a,es15.7e3)') &
'Iter ', iter, ', Rel Err: ', trdiff/trG
'Iter ', iter, ', Rel Err: ', trdiff
if ((trdiff < tol_iter*sqrt(dble(nact))) &
if ((trdiff < tol_iter) &
.or. (iter == max_iter)) then
call rci_op_null(task)
ijob = SID_FINISH
......@@ -1372,7 +1373,7 @@ contains
write(*,*) 'Eigenvalues are:'
print *, resvec(1:n)
end if
call rci_op_converge(task, trdiff < tol_iter*sqrt(dble(nact)))
call rci_op_converge(task, trdiff < tol_iter)
return
end if
......
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