subroutine setup_r1(comm,t_start,y_start,t_end,tolerance,thresholds, &
method,task,error_assess,h_start,message)
!
! 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
!
real(kind=wp), intent(in) :: t_end, t_start !indep!
real(kind=wp), intent(in) :: tolerance
real(kind=wp), dimension(:), intent(in) :: y_start !dep!
real(kind=wp), dimension(:), intent(in) :: thresholds !shp-dep!
type(rk_comm_real_1d) :: comm
real(kind=wp), intent(in), optional :: h_start !indep!
logical, intent(in), optional :: error_assess, message
character(len=*), intent(in), optional :: task, method
!
character(len=*), parameter :: srname="SETUP"
!
real(kind=wp) :: hmin !indep!
real(kind=wp) :: cdiff
integer :: ier, nrec, tr_dim_of_stages
logical :: legalt
character(len=1) :: task1, method1
!
integer, parameter :: not_ready=-1, fatal=911, just_fine=1
real(kind=wp), parameter :: zero=0.0_wp, pt01=0.01_wp, fivepc=0.05_wp, &
third=1.0_wp/3.0_wp, one=1.0_wp, two=2.0_wp, ten=10.0_wp, fifty=50.0_wp
!
ier = just_fine; nrec = 0
!
! Clear previous state of the suite.
!
call setup_global_stuff
nullify(comm%thresh, comm%err_estimates, comm%weights, comm%y_old, &
comm%scratch, &
comm%y, comm%yp, comm%y_new, comm%yp_old, comm%stages, comm%ge_y, &
comm%ge_yp, comm%ge_err_estimates, comm%ge_assess, comm%ge_y_new, &
comm%ge_stages, comm%v0, comm%v1, comm%v2, comm%v3, comm%vtemp, &
comm%xstage, comm%ytemp, comm%p)
!
! Fetch output channel and machine constants;
!
call machine_const(comm%round_off,comm%sqrrmc,comm%cubrmc,comm%sqtiny, &
comm%outch)
!
body: do
!
! Check for valid shape
if (size(shape(y_start))>0) then
if (any(shape(y_start)==0)) then
ier = fatal; nrec = 2; write (comm%rec,"(a)") &
" ** An extent of Y_START has zero length. This is not permitted."
exit body
end if
end if
!
! Check and process non-trivial optional arguments
if (present(task)) then
task1 = task(1:1); comm%use_range = task1 == "R" .or. task1 == "r"
legalt = comm%use_range .or. task1 == "S" .or. task1 == "s"
if (.not.legalt) then
ier = fatal; nrec = 2; write (comm%rec,"(a,a,a/a)") &
" ** You have set the first character of TASK to be '",TASK1,"'.", &
" ** It must be one of 'R','r','S' or 's'."
exit body
end if
end if
if (present(method)) then
method1 = method(1:1)
select case (method1)
case("L","l"); comm%rk_method = 1
case("M","m"); comm%rk_method = 2
case("H","h"); comm%rk_method = 3
case default
ier = fatal; nrec = 2; write (comm%rec,"(a,a,a/a)") &
" ** You have set the first character of METHOD to be '",METHOD1,"'.", &
" ** It must be one of 'L','l','M','m','H' or 'h'."
exit body
end select
end if
if (present(message)) comm%print_message = message
!
! Check consistency of array arguments
!
if (any(shape(y_start)/=shape(thresholds))) then
ier = fatal; nrec = 1; write (comm%rec,"(a)") &
" ** The shapes of Y_START and THRESHOLDS are not consistent."
exit body
end if
!
! Check and process compulsory arguments
if (t_start == t_end) then
ier = fatal; nrec = 1; write (comm%rec,"(a,e13.5,a)") &
" ** You have set T_START = T_END = ",T_START,"."
exit body
else
comm%t_end = t_end; comm%t_start = t_start
comm%t_old = t_start; comm%t = t_start
comm%dir = sign(one,t_end-t_start)
end if
if ((tolerance > pt01) .or. (tolerance < comm%round_off)) then
ier = fatal; nrec = 2; write (comm%rec,"(a,e13.5,a/a,e13.5,a)") &
" ** You have set TOLERANCE = ",tolerance," which is not permitted. The", &
" ** range of permitted values is (",comm%round_off,",0.01)."
exit body
else
comm%tol = tolerance
end if
if (minval(thresholds) < comm%sqtiny) then !spec-ar!
ier = fatal; nrec = 2; write (comm%rec,"(a/a,e13.5,a)") &
" ** You have set a component of THRESHOLDS to be less than the permitted", &
" ** minimum,",comm%sqtiny,"."
exit body
end if
!
! Set formula definitions and characteristics
call method_const(comm%rk_method, comm%a, comm%b, comm%c, comm%bhat, &
comm%r, comm%e, comm%ptr, comm%no_of_stages, comm%intrp_degree, &
comm%intrp_able, comm%intrp_needs_stages, comm%cost, &
comm%safety, comm%expon, comm%stability_radius, comm%tan_angle, &
comm%rs, comm%rs1, comm%rs2, comm%rs3, comm%rs4, comm%order, &
comm%last_stage, comm%max_stiff_iters, comm%no_of_ge_steps, comm%fsal,&
cdiff)
!
tr_dim_of_stages = maxval(comm%ptr(2:comm%no_of_stages))
comm%toosml = comm%round_off/cdiff
!
! In STEP_INTEGRATE the first step taken will have magnitude H. If
! H_START = ABS(H_START) is not equal to zero, H = H_START. If H_START is
! equal to zero, the code is to find an on-scale initial step size H. To
! start this process, H is set here to an upper bound on the first step
! size that reflects the scale of the independent variable.
! RANGE_INTEGRATE has some additional information, namely the first output
! point, that is used to refine this bound in STEP_INTEGRATE when called
! from RANGE_INTEGRATE. If H_START is not zero, but it is either too big
! or too small, the input H_START is ignored and H_START is set to zero to
! activate the automatic determination of an on-scale initial step size.
!
hmin = max(comm%sqtiny,comm%toosml*max(abs(t_start),abs(t_end)))
if (abs(t_end-t_start) < hmin) then
ier = fatal; nrec = 4; write (comm%rec,"(a/a/a,e13.5,a/a,e13.5,a)") &
" ** You have set values for T_END and T_START that are not clearly", &
" ** distinguishable for the method and the precision of the computer", &
" ** being used. ABS(T_END-T_START) is ",ABS(T_END-T_START)," but should be", &
" ** at least ",hmin,"."
exit body
end if
if (present(h_start)) comm%h_start = abs(h_start)
if (comm%h_start > abs(t_end-t_start) .or. comm%h_start < hmin) &
comm%h_start = zero
if (comm%h_start == zero) then
comm%h = max(abs(t_end-t_start)/comm%rs3,hmin)
else
comm%h = comm%h_start
end if
!
! Allocate a number of arrays using pointers.
!
allocate(comm%thresh(size(y_start,1)), & !alloc!
comm%err_estimates(size(y_start,1)), & !alloc!
comm%weights(size(y_start,1)), & !alloc!
comm%y_old(size(y_start,1)), & !alloc!
comm%scratch(size(y_start,1)), & !alloc!
comm%y(size(y_start,1)), & !alloc!
comm%yp(size(y_start,1)), & !alloc!
comm%stages(size(y_start,1),tr_dim_of_stages), & !alloc!
comm%ymax(size(y_start,1)),stat=ier) !alloc!
if (ier /= 0) then
ier = fatal; nrec = 1 ; write (comm%rec,"(a)") &
" ** Not enough storage available to create workspace required internally."
exit body
else
comm%y = y_start; comm%ymax = abs(y_start)
comm%thresh = thresholds;
comm%y_new => comm%scratch; comm%yp_old => comm%scratch
comm%v0 => comm%err_estimates; comm%vtemp => comm%scratch
comm%v1 => comm%stages(:,1); comm%v2 => comm%stages(:,2)
comm%v3 => comm%stages(:,3)
end if
!
! Pre-allocate storage for interpolation if the TASK = `R' was specified.
!
if (comm%use_range) then
if (comm%intrp_able) then
if (comm%rk_method==1) then
comm%p => comm%stages(:,1:2)
else if (comm%rk_method==2) then
allocate(comm%p(size(y_start,1),5), & !alloc!
comm%ytemp(size(y_start,1)), & !alloc!
comm%xstage(size(y_start,1)),stat=ier) !alloc!
if (ier /= 0) then
ier = fatal; nrec = 1 ; write (comm%rec,"(a)") &
" ** Not enough storage available to create workspace required internally."
exit body
end if
end if
end if
end if
!
! Initialise state and allocate storage for global error assessment
!
comm%t_ge_max_contrib = t_start
if (present(error_assess)) comm%erason = error_assess
if (comm%erason) then
!
! Storage is required for the stages of a secondary integration. The
! stages of the primary intergration can only be overwritten in the
! cases where there is no interpolant or the interpolant does not
! require information about the stages (e.g. METHOD 'H' and METHOD 'L',
! respectively).
if (.not.comm%intrp_needs_stages) then
comm%ge_stages => comm%stages
else
allocate(comm%ge_stages(size(y_start,1),tr_dim_of_stages),stat=ier) !alloc!
if (ier /= 0) then
ier = fatal; nrec = 1 ; write (comm%rec,"(a)") &
" ** Not enough storage available to create workspace required internally."
exit body
end if
end if
allocate(comm%ge_y(size(y_start,1)), & !alloc!
comm%ge_yp(size(y_start,1)), & !alloc!
comm%ge_err_estimates(size(y_start,1)), & !alloc!
comm%ge_assess(size(y_start,1)), & !alloc!
comm%ge_y_new(size(y_start,1)),stat=ier) !alloc!
if (ier /= 0) then
ier = fatal; nrec = 1 ; write (comm%rec,"(a)") &
" ** Not enough storage available to create workspace required internally."
exit body
else
comm%ge_assess = 0.0_wp; comm%ge_y = y_start
end if
end if
exit body
end do body
!
call rkmsg_r1(ier,srname,nrec,comm)
!
contains
subroutine setup_global_stuff
!
comm%h_start=0.0_wp; comm%h_old=0.0_wp
comm%f_count=0; comm%full_f_count=0; comm%step_count=0; comm%bad_step_count=0
comm%at_t_start=.true.; comm%at_t_end=.false.
comm%rk_method=2;
comm%ge_max_contrib=0.0_wp; comm%ge_f_count=0
comm%erason=.false.; comm%erasfl=.false.
comm%print_message=.true.; comm%use_range=.true.
comm%stiff_bad_step_count=0; comm%hit_t_end_count=0
comm%errold=0.0_wp; comm%h_average=0.0_wp
comm%chkeff=.false.; comm%phase2=.true.
comm%save_states(1:7)=not_ready
comm%stop_on_fatal=.true.; comm%saved_fatal_err=.false.
!
end subroutine setup_global_stuff
end subroutine setup_r1