sturm_count Function

private pure function sturm_count(igu, igl, kc, lam) result(n_chg)

Number of Sturm-sequence sign changes (eigenvalues < lam) via the three-term determinant recursion with dynamic rescaling.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: igu(NZ_STACK_MAX+1)
real(kind=wp), intent(in) :: igl(NZ_STACK_MAX+1)
integer, intent(in) :: kc
real(kind=wp), intent(in) :: lam

Return Value integer


Called by

proc~~sturm_count~~CalledByGraph proc~sturm_count sturm_count proc~wavespeed_cg1_column wavespeed_cg1_column proc~wavespeed_cg1_column->proc~sturm_count proc~wavespeed_compute_impl wavespeed_compute_impl proc~wavespeed_compute_impl->proc~wavespeed_cg1_column proc~wavespeed_compute wavespeed_compute proc~wavespeed_compute->proc~wavespeed_compute_impl proc~ocean_dyn_step_split ocean_dyn_step_split proc~ocean_dyn_step_split->proc~wavespeed_compute proc~engine_step engine_step proc~engine_step->proc~ocean_dyn_step_split

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private :: absd
integer, private :: cur_s
real(kind=wp), private :: det
real(kind=wp), private :: detm1
real(kind=wp), private :: detm2
real(kind=wp), private :: dval
integer, private :: k
real(kind=wp), private :: offl
integer, private :: prev_s
real(kind=wp), private :: sfac

Source Code

   pure integer function sturm_count(igu, igl, kc, lam) result(n_chg)
      !! Number of Sturm-sequence sign changes (eigenvalues < lam) via
      !! the three-term determinant recursion with dynamic rescaling.
      !$acc routine seq
      real(wp), intent(in) :: igu(NZ_STACK_MAX + 1), igl(NZ_STACK_MAX + 1)
      integer, intent(in) :: kc
      real(wp), intent(in) :: lam
      integer :: k, prev_s, cur_s
      real(wp) :: det, detm1, detm2, dval, offl, absd, sfac

      n_chg = 0
      detm1 = 1.0_wp
      det = (igu(2) + igl(2)) - lam
      prev_s = 1
      cur_s = merge(1, -1, det >= 0.0_wp)
      if (cur_s /= prev_s) n_chg = n_chg + 1
      prev_s = cur_s
      do k = 3, kc
         absd = abs(det)
         sfac = 1.0_wp
         if (absd > RESCALE) then
            sfac = I_RESCALE
         else if (absd < I_RESCALE .and. absd > 0.0_wp) then
            sfac = RESCALE
         end if
         detm2 = C2_SCALE*detm1*sfac
         detm1 = C2_SCALE*det*sfac
         dval = (igu(k) + igl(k)) - lam
         offl = igu(k)*igl(k - 1)
         det = dval*detm1 - offl*detm2
         cur_s = merge(1, -1, det >= 0.0_wp)
         if (cur_s /= prev_s) n_chg = n_chg + 1
         prev_s = cur_s
      end do
   end function sturm_count