subroutine stiff_r1(comm,f,toomch,sure_stiff)
!
! Part of rksuite_90 v1.0 (Aug 1994)
! software for initial value problems in ODEs
!
! Authors: R.W. Brankin (NAG Ltd., Oxford, England)
! I. Gladwell (Math Dept., SMU, Dallas, TX, USA)
! see main doc for contact details
!
type(rk_comm_real_1d), intent(inout), target :: comm
logical, intent(in) :: toomch
logical, intent(out) :: sure_stiff
!
interface
function f(t,y)
use rksuite_90_prec, only:wp
real(kind=wp), intent(in) :: t !indep!
real(kind=wp), dimension(:), intent(in) :: y !dep!
real(kind=wp), dimension(size(y,1)) :: f !dep!
end function f
end interface
!
logical :: maybe_stiff, lots_of_fails
!
real(kind=wp) :: alpha1, alpha2, beta1, beta2 !dep!
real(kind=wp) :: rold, v1v0, v2v0, v2v1, v3v1, v3v2 !dep!
real(kind=wp) :: dist, res2, scale, v0nrm, v3nrm, ynrm, rho, v0v0, v1v1, &
v2v2, v3v3, yy, det1, det2
integer :: ntry
complex(kind=wp), dimension(2) :: root_pair, prv_root_pair
integer, parameter :: bigr=1, smlr=2
real(kind=wp), dimension(:), pointer :: v0, v1, v2, v3, y, y_old !dep!
real(kind=wp), dimension(:), pointer :: weights, thresh !shp-dep!
!
integer, parameter :: max_f_count=5000
real(kind=wp),parameter :: zero=0.0_wp, pt001=0.001_wp, pt9=0.9_wp, &
fifth=0.2_wp, half=0.5_wp, one=1.0_wp, two=2.0_wp, five=5.0_wp, &
large=1.0e+10_wp
!
v0 => comm%v0
v1 => comm%v1
v2 => comm%v2
v3 => comm%v3
weights => comm%weights
thresh => comm%thresh
y => comm%y
y_old => comm%y_old
!
sure_stiff = .false.
lots_of_fails = .false.
!
if (mod(comm%step_count-10,40)==0) then
lots_of_fails = comm%stiff_bad_step_count >= 10
comm%stiff_bad_step_count = 0
end if
!
! If either too much work has been done or there are lots of failed steps,
! test for stiffness.
!
maybe_stiff = toomch .or. lots_of_fails
if (maybe_stiff) then
!
! Regenerate weight vector
!
weights = max(half*(abs(y)+abs(y_old)),thresh)
maybe_stiff = fifth < abs(comm%h/comm%h_average) .and. &
abs(comm%h/comm%h_average) < five
if (maybe_stiff) then
!
! The average step size is used to predict the cost in function evaluations
! of finishing the integration to T_END. If this cost is no more than
! MAX_F_COUNT, the problem is declared not stiff. If the step size is
! being restricted on grounds of stability, it will stay close to H_AVERAGE.
! The prediction will then be good, but the cost is too low to consider
! the problem stiff. If the step size is not close to H_AVERAGE, the
! problem is not stiff. Either way there is no point to testing for a step
! size restriction due to stability.
!
maybe_stiff = comm%cost*abs((comm%t_end-comm%t)/comm%h_average) > &
real(max_f_count,kind=wp)
if (maybe_stiff) then
!
! There have been many step failures or a lot of work has been done. Now
! we must determine if this is due to the stability characteristics of the
! formula. This is done by calculating the dominant eigenvalues of the
! local Jacobian and then testing whether H_AVERAGE corresponds to being
! on the boundary of the stability region.
! The size of Y provides scale information needed to approximate
! the Jacobian by differences.
!
v0v0 = wt_inner_prod(v0,v0)
yy = wt_inner_prod(y,y)
ynrm = sqrt(yy)
scale = ynrm*comm%sqrrmc
if (scale==zero) then
!
! Degenerate case. Y is (almost) the zero vector so the scale is not
! defined. The input vector V0 is the difference between Y and a
! lower order approximation to the solution that is within the error
! tolerance. When Y vanishes, V0 is itself an acceptable approximate
! solution, so we take SCALE from it, if this is possible.
!
scale = v0v0*comm%sqrrmc
maybe_stiff = scale > zero
end if
end if
end if
end if
!
if (.not. maybe_stiff) return
!
if (v0v0==zero) then
!
! Degenerate case. V0 is (almost) the zero vector so cannot
! be used to define a direction for an increment to Y. Try a
! "random" direction.
!
v0 = one; v0v0 = wt_inner_prod(v0,v0)
end if
!
v0nrm = sqrt(v0v0)
v0 = v0/v0nrm; v0v0 = one
!
! Use a nonlinear power method to estimate the two dominant eigenvalues.
! V0 is often very rich in the two associated eigenvectors. For this
! reason the computation is organized with the expectation that a minimal
! number of iterations will suffice. Indeed, it is necessary to recognize
! a kind of degeneracy when there is a dominant eigenvalue. The function
! DOMINANT_EIGENVALUE does this. In the first try, NTRY = 1, a Rayleigh
! quotient for such an eigenvalue is initialized as ROLD. After each
! iteration, DOMINANT_EIGENVALUE computes a new Rayleigh quotient and
! tests whether the two approximations agree to one tenth of one per cent
! and the eigenvalue, eigenvector pair satisfy a stringent test on the
! residual.
!
ntry = 1
do
!
v1 = approx_jacobian(f,v0,v0v0)
v1v1 = wt_inner_prod(v1,v1)
!
! The quantity SQRT(V1V1/V0V0) is a lower bound for the product of H_AVERAGE
! and a Lipschitz constant. If it should be LARGE, stiffness is not
! restricting the step size to the stability region. The principle is
! clear enough, but the real reason for this test is to recognize an
! extremely inaccurate computation of V1V1 due to finite precision
! arithmetic in certain degenerate circumstances.
!
if (sqrt(v1v1) > large*sqrt(v0v0)) return
!
v1v0 = wt_inner_prod(v1,v0)
if (ntry==1) then
rold = v1v0/v0v0
!
! This is the first Rayleigh quotient approximating the product of H_AVERAGE
! and a dominant eigenvalue. If it should be very small, the
! problem is not stiff. It is important to test for this possibility so
! as to prevent underflow and degeneracies in the subsequent iteration.
!
if (abs(rold) < comm%cubrmc) return
else
!
if (dominant_eigenvalue(v1v1,v1v0,v0v0)) exit
end if
!
v2 = approx_jacobian(f,v1,v1v1)
v2v2 = wt_inner_prod(v2,v2)
v2v0 = wt_inner_prod(v2,v0)
v2v1 = wt_inner_prod(v2,v1)
if (dominant_eigenvalue(v2v2,v2v1,v1v1)) exit
!
! Fit a quadratic in the eigenvalue to the three successive iterates
! V0,V1,V2 of the power method to get a first approximation to
! a pair of eigenvalues. A test made earlier in DOMINANT_EIGENVALUE
! implies that the quantity DET1 here will not be too small.
!
det1 = v0v0*v1v1 - v1v0*rev_wt_inner_prod(v1v0)
alpha1 = (-v0v0*v2v1 + rev_wt_inner_prod(v1v0)*v2v0)/det1
beta1 = (v1v0*v2v1 - v1v1*v2v0)/det1
!
! Iterate again to get V3, test again for degeneracy, and then fit a
! quadratic to V1,V2,V3 to get a second approximation to a pair
! of eigenvalues.
!
v3 = approx_jacobian(f,v2,v2v2)
v3v3 = wt_inner_prod(v3,v3)
v3v1 = wt_inner_prod(v3,v1)
v3v2 = wt_inner_prod(v3,v2)
if (dominant_eigenvalue(v3v3,v3v2,v2v2)) exit
!
det2 = v1v1*v2v2 - v2v1*rev_wt_inner_prod(v2v1)
alpha2 = (-v1v1*v3v2 + rev_wt_inner_prod(v2v1)*v3v1)/det2
beta2 = (v2v1*v3v2 - v2v2*v3v1)/det2
!
! First test the residual of the quadratic fit to see if we might
! have determined a pair of eigenvalues.
!
res2 = abs( v3v3 + rev_wt_inner_prod(alpha2)*v3v2 + &
rev_wt_inner_prod(beta2)*v3v1 + &
alpha2*rev_wt_inner_prod(v3v2) + &
alpha2*rev_wt_inner_prod(alpha2)*v2v2 + &
alpha2*rev_wt_inner_prod(beta2)*v2v1 + &
beta2*rev_wt_inner_prod(v3v1) + &
beta2*rev_wt_inner_prod(alpha2)*rev_wt_inner_prod(v2v1) + &
beta2*rev_wt_inner_prod(beta2)*v1v1 )
if (res2 <= abs(v3v3)*pt001**2) then
!
! Calculate the two approximate pairs of eigenvalues.
!
prv_root_pair(1:2) = quadratic_roots(alpha1,beta1)
root_pair(1:2) = quadratic_roots(alpha2,beta2)
!
! The test for convergence is done on the larger root of the second
! approximation. It is complicated by the fact that one pair of roots
! might be real and the other complex. First calculate the spectral
! radius RHO of HAVG*J as the magnitude of ROOT1. Then see if one of
! the roots R1,R2 is within one per cent of ROOT1. A subdominant root
! may be very poorly approximated if its magnitude is much smaller than
! RHO -- this does not matter in our use of these eigenvalues.
!
rho = abs(prv_root_pair(bigr))
dist = min( abs(root_pair(bigr) - prv_root_pair(bigr)), &
abs(root_pair(bigr) - prv_root_pair(smlr)) )
if (dist <= pt001*rho) exit
end if
!
! Do not have convergence yet. Because the iterations are cheap, and
! because the convergence criterion is stringent, we are willing to try
! a few iterations.
!
ntry = ntry + 1
if (ntry > comm%max_stiff_iters) return
v3nrm = sqrt(v3v3)
v0 = v3/v3nrm
v0v0 = one
!
end do
!
! We now have the dominant eigenvalues. Decide if the average step
! size is being restricted on grounds of stability. Check the real
! parts of the eigenvalues. First see if the dominant eigenvalue is
! in the left half plane -- there won't be a stability restriction
! unless it is. If there is another eigenvalue of comparable magnitude
! with a positive real part, the problem is not stiff. If the dominant
! eigenvalue is too close to the imaginary axis, we cannot diagnose
! stiffness.
!
if ( real(root_pair(bigr)) < zero) then
if ( .not. ( abs(root_pair(smlr)) >= pt9*rho .and. &
real(root_pair(smlr)) > zero ) ) then
if ( abs(aimag(root_pair(bigr))) <= &
abs(real(root_pair(bigr)))*comm%tan_angle) then
!
! If the average step size corresponds to being well within the
! stability region, the step size is not being restricted because
! of stability.
!
sure_stiff = rho >= pt9*comm%stability_radius
end if
end if
end if
!
contains
function approx_jacobian(f,v,vdotv)
!
real(kind=wp), intent(in) :: vdotv
real(kind=wp), dimension(:), intent(in) :: v !dep!
real(kind=wp), dimension(size(v,1)) :: approx_jacobian !dep!
!
interface
function f(t,y)
use rksuite_90_prec, only:wp
real(kind=wp), intent(in) :: t !indep!
real(kind=wp), dimension(:), intent(in) :: y !dep!
real(kind=wp), dimension(size(y,1)) :: f !dep!
end function f
end interface
!
real(kind=wp) :: temp1
!
! Scale V so that it can be used as an increment to Y
! for an accurate difference approximation to the Jacobian.
!
temp1 = scale/sqrt(vdotv)
comm%vtemp = y + temp1*v
!
approx_jacobian = f(comm%t,comm%vtemp)
comm%f_count = comm%f_count + 1
!
! Form the difference approximation. At the same time undo
! the scaling of V and introduce the factor of H_AVERAGE.
!
approx_jacobian = &
(comm%h_average/temp1)*(approx_jacobian-comm%yp)
!
end function approx_jacobian
function quadratic_roots(alpha,beta)
!
real(kind=wp), intent(in) :: alpha, beta !dep!
complex(kind=wp), dimension(2) :: quadratic_roots
!
complex(kind=wp) :: temp, sqdisc, r1, r2
!
! For types other than real/complex, this procedure must be constructed
! such that temp and sqdisc are evaluated as compelx quantities
!
temp = alpha/two; sqdisc = sqrt(temp**2 - beta)
!
! Do we have double root?
!
if (sqdisc==zero) then
quadratic_roots = (/ -temp, -temp /)
!
! Distinct roots
!
else
r1 = -temp + sqdisc; r2 = -temp + sqdisc
if (abs(r1) > abs(r2)) then
quadratic_roots = (/ r1, r2 /)
else
quadratic_roots = (/ r2, r1 /)
end if
end if
!
end function quadratic_roots
function dominant_eigenvalue(v1v1,v1v0,v0v0)
!
real(kind=wp), intent(in) :: v0v0, v1v1
real(kind=wp), intent(in) :: v1v0 !dep!
logical :: dominant_eigenvalue
!
real(kind=wp) :: ratio !dep!
real(kind=wp) :: res, det
logical :: big
!
ratio = v1v0/v0v0; rho = abs(ratio)
det = v0v0*v1v1 - v1v0*rev_wt_inner_prod(v1v0); res = abs(det/v0v0)
!
big = det == zero .or. &
(res<=abs(v1v1)*pt001**2 .and. abs(ratio-rold)<=pt001*rho)
!
if (big) then
root_pair(bigr) = cmplx(ratio)
root_pair(smlr) = cmplx(zero)
end if
!
rold = ratio
dominant_eigenvalue = big
!
end function dominant_eigenvalue
function wt_inner_prod(vec_1,vec_2)
!
real(kind=wp), dimension(:), intent(in) :: vec_1, vec_2 !dep!
real(kind=wp) :: wt_inner_prod !dep!
!
!
wt_inner_prod = sum ( (vec_1/weights) * (vec_2/weights) ) !spec-ar!
!
end function wt_inner_prod
function rev_wt_inner_prod(value)
!
real(kind=wp), intent(in) :: value !dep!
real(kind=wp) :: rev_wt_inner_prod !dep!
!
! given result of inner product value = v1.v0
! must return the reverse, ie v0.v1
!
! for real variables the value is the same
! for complex need to conjugate
!
rev_wt_inner_prod = value !spec-line!
!
end function rev_wt_inner_prod
end subroutine stiff_r1