subroutine mirkdc(method,tol,neqns,Nsub,MxNsub,
+ mesh_input,mesh,sol_input,Y,ldY,output_control,info,
+ iwork,work,init_de,init_Y,fsub,gsub,dfsub,dgsub)
c
c MIRKDC_4, Dec. 1999.
c
c Authors: Wayne Enright (enright@cs.toronto.edu),
c Paul Muir (muir@stmarys.ca).
c Additional programming assistance: Mark Adams, Nicole DeYoung.
c Thanks to Beta-Testers: Ian Gladwell, Larry Shampine, Richard Pancer.
c
c Note: This code makes use of the Level 1 BLAS routines. These must
c be linked with the MIRKDC package in order for it to work.
c
c***declaration of constants***
integer Mxs
parameter(Mxs=10)
c
c `Mxs' Maximum number of stages of the Runge-Kutta method.
c
c***declaration of parameters***
c imports:
integer method
double precision tol
integer neqns, Nsub, MxNsub
integer mesh_input
double precision mesh(0:MxNsub)
integer sol_input
double precision Y(ldY*(Nsub+1))
integer ldY, output_control
c
c `method' Defines the MIRK method to be used. Possible values
c for `method' and the corresponding MIRK schemes are:
c method | MIRK formula
c ---------------------------------------------------------
c Externally documented methods.
c 2 | MIRK221 scheme (trapezoidal rule, second order)
c 4 | MIRK343 scheme (Lobatto scheme, fourth order)
c 6 | MIRK563 scheme (6th order)
c Internally available methods.
c 121 | MIRK121 scheme (midpoint rule, second order)
c 221 | MIRK221 scheme (trapezoidal rule, second order)
c 232 | MIRK232 scheme (one-sided, abscissa 0, 2/3)
c 2321 | MIRK232 scheme (one-sided, abscissa 1, 1/3)
c 343 | MIRK343 scheme (Lobatto scheme, fourth order)
c 453 | MIRK453 scheme (one-sided, 5th order)
c 563 | MIRK563 scheme (6th order)
c
c `tol' Tolerance applied to an estimate of the defect of the approximate
c solution; the defect is the amount by which the continuous
c approximate solution, u(t), fails to satisfy the ODE system.
c We require |def(t)|/(|f(t,u(t))|+1) to be less than `tol', where
c y'(t)=f(t,y(t)) is the ODE and def(t)= u'(t)-f(t,u(t)).
c The same tolerance is applied to all defect components.
c
c `neqns' Number of differential equations. Also number of boundary
c conditions.
c `Nsub' On input, number of subintervals into which the problem
c interval is initially partitioned.
c `MxNsub' Maximum number of subintervals.
c `mesh_input' a) A value of 0 indicates that a uniform initial mesh
c should be set up by `mirkdc'.
c b) A value of 1 indicates that an initial mesh has
c already been provided by the user in the array `mesh'.
c This option should be chosen if the user has specific
c knowledge of the solution behavior.
c c) A value of 2 indicates that the user is employing
c continuation to solve the problem and that the mesh
c from a previous run will be stored in the array `mesh'.
c `mesh' On input, the set of points that partition the problem
c interval; initially empty if `mesh_input' = 0.
c `sol_input' a) Normal use of `mirkdc' requires that `sol_input'=1,
c indicating that an initial solution estimate is to be
c provided by the user through the `init_Y' routine. It
c is recommended that `Y' not be set to a constant.
c b) The option `sol_input' = 2 should only be chosen when
c using continuation to solve a problem, where `Y' takes
c its value from the solution of a previous run. A value
c of 2 indicates that an initial solution estimate at each
c mesh point is available in the array `Y'.
c `Y' On input, if `sol_input'=2, this is the initial discrete
c approximation to the solution, evaluated at each meshpoint.
c Externally, in the interfaces with the user main program and
c subroutines, `Y' is a two-dimensional array; the i-th column of `Y'
c stores the solution approximation for the i-th mesh point. Each
c solution approximation is of length `neqns'; the resulting
c dimensions of `Y' are `neqns x (Nsub+1)'. Inside MIRKDC it is
c convenient for `Y' to be stored as a column vector; in this
c one-dimensional version the solutions at each meshpoint come one
c after the other rather than side by side as they are in the
c two-dimensional array.
c `ldY' The leading dimension of Y. (`ldY' >= `neqns').
c `output_control' Controls the amount of information output. Provides
c profiling of code execution for activities such as the
c Newton iteration, mesh selection, or defect estimation.
c 0 - None. 1 - Intermediate. 2 - Full.
c exports:
c double precision mesh(0:Nsub)
c double precision Y(ldY*(Nsub+1))
integer info
c
c `mesh' On output, the set of points which define the final mesh.
c `Y' On output, if `info'=0, the converged discrete solution, evaluated
c on the final mesh.
c `info' Communication flag:
c 0 - implies a successful termination - solution obtained such
c that the defect is within user specified tolerance.
c -1 - implies an unsuccessful termination - the required size for
c the new mesh is be too large.
c Internal flag:
c The Newton iteration can fail for one of several reasons, each
c signaled by a non-zero `info' value. Even if the Newton
c iteration converges, the resultant solution may be
c unsuitable because the defect is too large. During the
c execution of `mirkdc', `info' may take on the values given
c below. None of these conditions lead directly to an
c unsuccessful termination. Rather, the code will react by
c attempting to double the size of the mesh and try again.
c 1 - too many iterations were taken without convergence being
c obtained.
c 2 - a singular coefficient matrix was encountered during the
c attempted solution of the Newton system.
c 3 - (a) it was impossible to obtain a suitable damping factor
c for the Newton update (indicative of an "effectively
c singular" Jacobian) or (b) evaluation of natural criterion
c function overflowed, indicative of a divergent iteration.
c 4 - the defect was too large and the solution from the Newton
c system could not be trusted.
c work space:
integer iwork(neqns*(MxNsub+1))
double precision work(2*MxNsub*Mxs*neqns +
+ 5*neqns*MxNsub + 2*neqns**2*MxNsub +
+ 6*neqns + 4*neqns**2 + 2*neqns**2*Mxs)
c
c `iwork' Integer work array provided to the `newiter' routine.
c `work' Double precision work array used within `mirkdc' and provided to
c `NewIter', `defect_estimate', `mesh_selector', and `interp_eval'.
c The largest amount of workspace used is during the call to
c `Newiter'. The space required then is at most
c 2 * MxNsub*Mxs*neqns + 5*neqns*MxNsub +
c 2*neqns^2*MxNsub + 6*neqns + 4*neqns^2 + 2*neqns^2*Mxs.
c
c user-supplied subroutines:
external init_de,init_Y,fsub,gsub,dfsub,dgsub
c
c `init_de' User-supplied subroutine to initialize the problem dependent
c variables, `leftbc', `a', and `b'.
c
c call init_de(leftbc, a, b)
c exports:
c integer leftbc
c double precision a,b
c `leftbc' Number of boundary conditions applied at the
c left endpoint of the problem interval.
c `a' Left endpoint of the problem interval.
c `b' Right endpoint of the problem interval.
c
c `init_Y' User-supplied subroutine to initialize the discrete solution
c approximation, `Y'. (Needed if `sol_input'=1).
c
c call init_Y(Nsub, neqns, mesh, Y, ldY)
c imports:
c integer Nsub, neqns
c double precision mesh(0:Nsub)
c integer ldY
c `Nsub' Number of subintervals in the initial mesh.
c `neqns' Number of differential equations in the
c first order system.
c `mesh' Points making up the partition of the problem
c interval.
c `ldY' Leading dimension of Y. (`ldY' >= `neqns').
c exports:
c double precision Y(ldY,Nsub+1)
c `Y' Initial approximation to the discrete
c solution at the points in the initial mesh.
c The i-th column of `Y' is the solution
c approximation vector at mesh point, `mesh(i-1)'.
c
c `fsub' User-supplied subroutine which defines f(t,y) for the first order
c system of differential equations, y' = f(t,y).
c
c call fsub(neqns, t, y, f)
c imports:
c integer neqns
c double precision t, y(neqns)
c `neqns' Number of differential equations in the
c first order system.
c `t' A point in the problem interval.
c `y' Current solution approximation at t.
c exports:
c double precision f(neqns)
c `f' Value of f(t,y).
c
c `gsub' User-supplied subroutine which defines the boundary condition
c equations.
c
c call gsub(neqns, ya, yb, bc)
c imports:
c integer neqns
c double precision ya(neqns), yb(neqns)
c `neqns' Number of differential equations in the
c first order system. Also the total number
c of boundary conditions.
c `ya' Current solution approximation at left endpoint.
c `yb' Current solution approximation at right endpoint.
c exports:
c double precision bc(neqns)
c `bc' Boundary condition equations evaluated at
c `ya' and `yb'. The first `leftbc' components
c of `bc' are `g_a(ya)'; the remaining
c `neqns-leftbc' components of `bc' are `g_b(yb)'.
c
c `dfsub' User-supplied subroutine which defines the Jacobian, df/dy,
c of the system of differential equations. Since the array
c `JJ' has already been set to zero, it is only necessary to
c set the non-zero elements.
c
c call dfsub(neqns, t, y, JJ)
c imports:
c integer neqns
c double precision t, y(neqns)
c `neqns' Number of differential equations in the
c first order system.
c `t' A point in the problem interval.
c `y' Current solution approximation at t.
c exports:
c double precision JJ(neqns,neqns)
c `JJ' Contains the Jacobian, df/dy.
c
c `dgsub' User-supplied subroutine which defines the Jacobian of
c the boundary conditions, that is d(bc)/dy. The array `dbc'
c has been set to zero prior to the call to `dgsub', so it
c is only necessary for the user to set the non-zero elements.
c
c call dgsub(neqns, ya, yb, dbc)
c imports:
c integer neqns
c double precision ya(neqns), yb(neqns)
c `neqns' Number of differential equations in the
c first order system. Also the total number
c of boundary conditions.
c `ya' Current solution approximation at left endpoint.
c `yb' Current solution approximation at right endpoint.
c exports:
c double precision dbc(neqns,neqns)
c `dbc' Contains the Jacobian of the boundary conditions.
c
c***declaration of local variables***
integer leftbc
double precision a,b,h
integer maxiter
double precision newton_tol
double precision defect_norm
integer Nsub_star
integer i,j
c
c `leftbc' The number of boundary conditions at the left end of the
c problem interval.
c `a' Left endpoint.
c `b' Right endpoint.
c `h' Initial subinterval size, when a uniform mesh is employed.
c `maxiter' Maximum number of Newton iterations.
c `newton_tol' Tolerance for each Newton iteration.
c `defect_norm' Max norm of the defect estimate.
c `Nsub_star' Number of subintervals of the new mesh.
c `i,j' Loop indices.
c
c***declaration of index related info for the work array***
c
c Arrays local to `mirkdc' to be represented within the `work' array:
c `k_discrete', `k_interp', `defect', `mesh_new', `Y_new'.
c Sizes:
integer s_k_discrete, s_k_interp, s_defect
integer s_mesh_new, s_Y_new
c
c The actual array declarations will be replaced by assignments giving
c the amount of space need within `work'.
c double precision k_discrete(Mxs*neqns*Nsub) =>
c `s_k_discrete' = `Mxs*neqns*Nsub'
c double precision k_interp(Mxs*neqns*Nsub) =>
c `s_k_interp' = `Mxs*neqns*Nsub'
c double precision defect(Nsub*neqns) => `s_defect' = `Nsub*neqns'
c double precision mesh_new(0:MxNsub) => `s_mesh_new' = `MxNsub+1'
c double precision Y_new((Nsub_star+1)*neqns) =>
c `s_Y_new' = `(Nsub_star+1)*neqns'
c
c `k_discrete' Stage information associated with the discrete
c Runge-Kutta method: the i-th set of `s*neqns' locations
c contain the `s' vectors of length `neqns' associated
c with the Runge-Kutta method on the i-th subinterval.
c (`Mxs' is the maximum number of stages; `s' is the
c actual number of stages.)
c `k_interp' Stage information associated with the continuous
c Runge-Kutta method; format is the same as `k_discrete'.
c `defect' Estimate of the (scaled) defect
c (u'(t) - f(t,u(t)))/(|f(t,u(t))|+1)
c on each subinterval (where u(t) is the continuous
c approximate solution). The defect estimate is a vector with
c `neqns' components and there is one such vector for each
c subinterval.
c `mesh_new' Points defining the new mesh.
c `Y_new' Discrete approximation to the solution on `mesh_new'
c
c Locations within `work':
integer i_k_discrete, i_k_interp
integer i_defect, i_mesh_new, i_Y_new
c
c `i_k_discrete' = 1;
c `i_k_interp' = `s_k_discrete' + 1
c `i_defect' = `s_k_discrete' + `s_k_interp' + 1
c `i_mesh_new' = `s_k_discrete' + `s_k_interp' + `s_defect' + 1
c `i_Y_new' = `s_k_discrete'+`s_k_interp'+`s_defect'+`s_mesh_new'+1
c
c Index related info for the representation within `work' of arrays local
c to `newiter'. See `newiter' for explanations of the use of the arrays.
c Sizes:
integer s_PHI, s_delta, s_top, s_bot, s_blocks
c
c `s_PHI' = `neqns * (Nsub+1)'
c `s_delta' = `neqns * (Nsub+1)'
c `s_top' = `neqns**2'
c `s_bot' = `neqns**2'
c `s_blocks' = `2 * neqns**2 * Nsub'
c
c Locations within `work' during call to `newiter'.
integer i_PHI, i_delta, i_top, i_bot, i_blocks
c
c `i_PHI' = `s_k_discrete + s_k_interp + 1'
c `i_delta' = `s_k_discrete + s_k_interp + s_PHI + 1'
c `i_top' = `s_k_discrete + s_k_interp + s_PHI + s_delta + 1'
c `i_bot' = `s_k_discrete + s_k_interp + s_PHI + s_delta + s_top + 1'
c `i_blocks' = `s_k_discrete+s_k_interp+s_PHI+s_delta+s_top+s_bot+1'
c
c Index related info for the work array for arrays local to `defect_estimate'.c See `defect_estimate' for explanations of the use of these arrays.
c Sizes:
c (All these arrays are of length `neqns'.)
c
c Locations within `work' during the call to `defect_estimate'.
integer i_z, i_z_prime, i_def_1, i_def_2
integer i_f_1, i_f_2, i_temp_1, i_temp_2
c
c `i_z = i_shift + 1; i_z_prime = i_shift + neqns + 1'
c `i_def_1 = i_shift + 2*neqns + 1; i_def_2 = i_shift + 3*neqns + 1'
c `i_f_1 = i_shift + 4*neqns + 1; i_f_2 = i_shift + 5*neqns + 1'
c `i_temp_1 = i_shift + 6*neqns + 1; i_temp_2 = i_shift + 7*neqns + 1'
c where `i_shift = s_k_discrete + s_k_interp + s_defect'
c
c Other work indexes.
integer i_work, i_shift
c
c `i_work' Index to free part of work array passed to various routines
c for use as a work array within those routines.
c `i_shift' Temporary variable used as a shift to simplify `work' index
c calculations.
c
c***declaration of variables for the /IOCNTRL/ common block***
c imports:
integer print_level_0, print_level_1, print_level_2
parameter (print_level_0 = 0, print_level_1 = 1,
* print_level_2 = 2)
integer profile
common /IOCNTRL/ profile
c
c `print_level_0' No output.
c `print_level_1' Intermediate output.
c `print_level_2' Full output.
c `profile' Controls output, to standard output, of profiling infor-
c mation such as Newton iteration counts, mesh selection,
c relative defect estimate. This variable is identical to
c the `output_control' parameter.
c
c
c The following two interfaces are included only to allow convenient loading
c of the integer work array upon successful exit. This work array is needed
c when the user calls the `sol_eval' routine, subsequent to the return from
c `mirkdc'.
c
c***declaration of variables from /RK_s/ common block***
c imports:
integer s
common /RK_s/ s
c `s' Number of discrete stages.
c
c***declaration of variables from /interp_s_star/ common block***
c imports:
integer s_star
common /interp_s_star/ s_star
c `s_star' Total number of stages required to form the
c interpolant on each subinterval. It includes all the
c stages of the discrete formula plus the additional
c stages required for the interpolant.
c----------------------------------------------------------------------------
c Called by : main
c Calls to : `rk_tableau',`interp_tableau',`init_de',`init_Y',`NewIter',
c `defect_estimate',`mesh_selector',`interp_eval',`half_mesh'
c----------------------------------------------------------------------------
c***initialization***
c
c Transfer output_control to common block variable `profile'.
profile = output_control
c
if (profile .GT. print_level_0) then
c Output required size of `work' and `iwork' based on given
c values for `neqns' and `MxNsub'.
write(6,*)
write(6,*) 'The input value for the order of the Runge-Kutta'
write(6,80) 'method is ',method,'.'
write(6,90) 'The input value for the number of ODEs is '
+ ,neqns,'.'
write(6,*) 'The input value for the maximum number of '
write(6,100) 'subintervals is ', MxNsub,'.'
write(6,*)
write(6,*) 'Based on these values the size of the double'
write(6,110) 'precision work array should be at least',
+ MxNsub*Mxs*neqns + 5*neqns*MxNsub + 2*neqns**2*MxNsub +
+ 6*neqns + 4*neqns**2 + 2*neqns**2*Mxs,'.'
write(6,*) 'and the size of the integer work array should'
write(6,120) 'be at least', neqns*(MxNsub+1),'.'
write(6,*)
80 format(1X,a10,i2,a1)
90 format(1x,a41,i3,a1)
100 format(1x,a15,i5,a1)
110 format(1x,a39,i7,a1)
120 format(1x,a11,i6,a1)
endif
c
c Set up formula dependent coefficients.
call rk_tableau(method)
call interp_tableau(method)
c
c Initialize remaining ODE dependent variables.
call init_de(leftbc,a,b)
c
c Setup initial mesh.
c If `mesh_input'=1 or 2 then there is already a mesh stored in `mesh'.
c If `mesh_input'=0 then setup uniform initial mesh of `Nsub' subintervals.
if (mesh_input.eq.0) then
mesh(0) = a
mesh(Nsub) = b
h = (b - a)/Nsub
do 10 i = 1, Nsub-1
mesh(i) = i*h + a
10 continue
endif
c
if (profile .GT. print_level_0) then
write(6,130) 'The size of the initial mesh is '
+ ,Nsub,' subintervals.'
write(6,*)
130 format(1x,a31,i5,a14)
end if
if (profile .GT. print_level_1) then
write(6,140)'The initial mesh of',Nsub,' subintervals:'
write(6,150)(mesh(i),i=0,Nsub)
write(6,*)
140 format(1x,a19,i5,a14)
150 format(7F10.6)
endif
c
c Initialize approximate solution.
c If `sol_input' = 2 then an approximate solution is already stored in `Y'.
c If `sol_input' = 1 then call `init_Y' to get an initial solution.
if (sol_input.eq.1) then
call init_Y(Nsub,neqns,mesh,Y,ldY)
endif
c
c Transfer the 2-D array, `Y', to its 1-D equivalent for use within MIRKDC.
do 15 i = 1,Nsub
do 20 j = 1,neqns
Y(i*neqns + j) = Y(i*ldY + j)
20 continue
15 continue
c
c Define the Newton tolerance and maximum iteration count
newton_tol = tol/100.0d0
maxiter = 40
c
c***main outer loop to control discrete problem sequence***
c REPEAT-UNTIL LOOP (Implemented here as a GOTO-IF pair)
1 continue
c
c Compute array sizes and indexes to provide access to the work
c array during the call to `newiter'. `s_defect' is also computed
c here since its value is required in more than one section of the code.
s_k_discrete = Mxs * neqns * Nsub
s_k_interp = Mxs * neqns * Nsub
i_k_discrete = 1
i_k_interp = s_k_discrete + 1
s_PHI = neqns * (Nsub+1)
s_delta = neqns * (Nsub+1)
s_top = neqns**2
s_bot = neqns**2
s_blocks = 2 * neqns**2 * Nsub
s_defect = Nsub*neqns
i_PHI = s_k_discrete + s_k_interp + 1
i_delta = s_k_discrete + s_k_interp + s_PHI + 1
i_top = s_k_discrete + s_k_interp + s_PHI +
+ s_delta + 1
i_bot = s_k_discrete + s_k_interp + s_PHI +
+ s_delta + s_top + 1
i_blocks = s_k_discrete + s_k_interp + s_PHI +
+ s_delta + s_top + s_bot + 1
i_work = s_k_discrete + s_k_interp + s_PHI +
+ s_delta + s_top + s_bot + s_blocks + 1
c
c Call `NewIter' to compute a solution for the discrete problem
c based on the current mesh.
c
if (profile .GT. print_level_0) then
write(6,*) 'Begin the Newton iteration:'
write(6,*)
end if
c
call NewIter(neqns,leftbc,Nsub,mesh,Y,
+ newton_tol,maxiter,info,work(i_PHI),work(i_delta),
+ work(i_top),work(i_bot),work(i_blocks),iwork,
+ work(i_k_discrete),work(i_work),fsub,gsub,dfsub,dgsub)
c
if ((info .NE. 0) .AND. (profile .EQ. print_level_1)) then
write(*,*)'the Newton iteration has failed.'
endif
if ((info .NE. 0 ) .AND. (profile .GT. print_level_1)) then
if (info .EQ. 1) then
write(6,*)'the Newton iteration failed because the',
+ ' maximum number of allowed iterations was exceeded.'
endif
if (info .EQ. 2) then
write(6,*)'the Newton iteration failed because a ',
+ 'singular matrix was encountered in the computation.'
endif
if (info .EQ. 3) then
write(6,*)'the Newton iteration failed because a ',
+ 'sufficiently large damping factor could not be found.'
endif
endif
c
c
c If the Newton iteration has converged, compute the defect estimate.
c
if (info .EQ. 0) then
c
c Set indexes to work array.
i_defect = s_k_discrete + s_k_interp + 1
i_shift = s_k_discrete + s_k_interp + s_defect
i_z = i_shift + 1
i_z_prime = i_shift + neqns + 1
i_def_1 = i_shift + 2*neqns + 1
i_def_2 = i_shift + 3*neqns + 1
i_f_1 = i_shift + 4*neqns + 1
i_f_2 = i_shift + 5*neqns + 1
i_temp_1 = i_shift + 6*neqns + 1
i_temp_2 = i_shift + 7*neqns + 1
i_work = i_shift + 8*neqns + 1
c
call defect_estimate(neqns,Nsub,mesh,Y,
+ work(i_defect), defect_norm,info,
+ work(i_z),work(i_z_prime),
+ work(i_def_1),work(i_def_2), work(i_f_1),
+ work(i_f_2),work(i_temp_1),work(i_temp_2),
+ work(i_k_discrete),work(i_k_interp),
+ work(i_work),fsub)
c
c The `defect_estimate' routine returns the defect estimate
c vector along with its max norm. It also checks the size of
c the defect norm to make sure it is less than a certain
c threshold value (see the `defect_estimate' routine for
c further details.) If the defect norm is too large the
c `defect_estimate' routine sets `info' to 4 to indicate that
c the defect is too large and that the solution should not be trusted.
c
if (profile .GT. print_level_0) then
write(6,*)'the Newton iteration was successful. ',
+ 'Build a continuous approximation to the solution and',
+ ' sample the defect.'
write(6,*)'The norm of defect is ',defect_norm
if (info .EQ. 4) then
write(6,*)'Since the defect is greater than ',
+ '10% the solution is not acceptable.'
endif
endif
c
endif
c End of `if-endif (info .EQ. 0)'
c
c
c If the Newton iteration failed or the defect was not acceptable,
c `info' will not be 0 and we will proceed to the `else' clause of the
c following `if-then-else' statement where mesh halving will be
c attempted. Otherwise, we will proceed to the `then' clause below
c and attempt mesh redistribution.
c
c Set work array indexes for `mesh_selector' and `half mesh'.
s_mesh_new = MxNsub + 1
i_mesh_new = s_k_discrete+s_k_interp
+ +s_defect+1
i_work = s_k_discrete+s_k_interp+s_defect
+ +s_mesh_new+1
i_Y_new = s_k_discrete+s_k_interp
+ +s_defect+s_mesh_new+1
c
c
if (info.EQ.0) then
c
c If we have not yet satisfied the tolerance, we must try again
c on a new mesh, if possible.
c
if (defect_norm .GT. tol) then
c
if (profile .GT. print_level_0) then
write(6,*) 'User defined tolerance',tol,' has not '
write(6,*) 'been satisfied.'
write(6,*) 'Construct a new mesh which ',
+ 'equidistributes the defect.'
endif
c
c Attempt to select a new mesh by calling `mesh_selector'.
call mesh_selector(neqns,Nsub,mesh,work(i_defect),
+ tol,Nsub_star,work(i_mesh_new),
+ info,MxNsub,work(i_work))
c
c If we were not successful in obtaining a new mesh, `info'
c will not equal 0 and we should terminate the computation.
c Otherwise we can compute a new solution estimate on the new
c mesh, and then update `Nsub', `mesh' and `Y'.
c
if (info .EQ. 0) then
if (profile .GT. print_level_0) then
write(6,*)
write(6,*)
write(6,160)'New mesh will be of size ',
+ Nsub_star,' subintervals.'
160 format(1x,a25,i5,a14)
end if
if (profile .GT. print_level_1) then
write(6,150)(work(i_mesh_new+i), i=0, Nsub_star)
write(6,*)
end if
c
c Use current computed solution to generate initial guess for next
c problem. Compute solution approximations at new mesh points.
c
c Set work array indexes.
s_Y_new = neqns*(Nsub_star+1)
i_Y_new = s_k_discrete+s_k_interp
+ +s_defect+s_mesh_new+1
i_work = s_k_discrete+s_k_interp
+ +s_defect+s_mesh_new+s_Y_new+1
c
do 30 i = 0, Nsub_star
call interp_eval(neqns,Nsub,mesh,Y,
+ work(i_mesh_new+i),
+ work(i_Y_new+i*neqns), work(i_work),
+ work(i_k_discrete),work(i_k_interp))
30 continue
c
c Copy `mesh_new' to `mesh', `Y_new' to `Y'; update `Nsub'.
call dcopy(Nsub_star+1,work(i_mesh_new),1,mesh,1)
call dcopy((Nsub_star+1)*neqns,work(i_Y_new),1,Y,1)
Nsub = Nsub_star
c
endif
c End of `if-endif (info .EQ. 0)'
c
endif
c End of `if-endif (defect_norm .GT. tol)'
c
if ( (info. EQ. 0) .AND. (defect_norm .LT. tol) .AND.
+ (profile .GT. print_level_0) ) then
write(6,*)'The user defined tolerance',tol,' has ',
+ 'been satisfied.'
write(6,*)'Successful computation. '
write(6,*)
endif
c
else
c Else clause associated with `if-then-else (info .EQ. 0)'
c
if (profile .GT. print_level_0) then
write(6,*)'Cannot obtain a solution for the current mesh.'
endif
c
c Simply half the current mesh, if possible, and try again.
if (2*Nsub .GT. MxNsub) then
c New mesh would be too large. Terminate the computation.
if (profile .GT. print_level_0) then
write(6,*)'New mesh would be too large.'
end if
info = -1
else
c Else clause associated with if_then_else `2*Nsub .GT. MxNsub'
c
call half_mesh(Nsub,mesh,work(i_mesh_new))
Nsub_star = 2*Nsub
c
c Copy `mesh_new' to `mesh'
call dcopy(Nsub_star+1,work(i_mesh_new),1,mesh,1)
Nsub = Nsub_star
c
if (profile .GT. print_level_0) then
write(6,*)'Double the mesh and try again.'
write(6,*)
write(6,*)
write(6,170)'New mesh will be of size', Nsub,'.'
170 format(1x,a24,i5,a1)
end if
if (profile .GT. print_level_1) then
write(6,150)(mesh(i), i = 0, Nsub)
write(6,*)
end if
c
c Initialize the approximate solution - store in `Y_new'.
if (sol_input .EQ. 1) then
call init_Y(Nsub,neqns,mesh,Y,ldY)
c Transferring the 2-D array, Y, to its 1-D equivalent for
c use with in MIRKDC.
do 35 i = 1,Nsub
do 40 j = 1,neqns
Y(i*neqns + j) = Y(i*ldY + j)
40 continue
35 continue
endif
c If sol_input = 2 then the initial guess on the initial mesh
c was provided through `Y' and `mesh', respectively. Since this
c initial guess is no longer usable, with the new mesh, we
c continue by simply setting `Y_new' to zero.
if (sol_input .EQ. 2) then
if (profile .GT. print_level_0) then
write(6,*)'Iteration failed and there is no ',
+ 'Init_Y routine; Y will be set to zero.'
write(6,*)'It would be better to retry the ',
+ 'continuation sequence with smaller steps'
endif
call dcopy((Nsub+1)*neqns,0.0d0,0,Y,1)
endif
c
c Reset `info' to 0 and `defect_norm' > `tol'
c to force a restart. (See UNTIL conditions below.)
info = 0
defect_norm = 2* tol
c
endif
c End of `if-then-else (2*Nsub .GT. MxNsub)'
c
endif
c End of `if-then-else (info .EQ. 0)', `Acceptable solution obtained'
c
c UNTIL( (info .NE. 0) .OR. (defect_norm .LE. tol))
if((info .EQ. 0) .AND.
+ (defect_norm .GT. tol))goto1
c
c***conclusion***
c
c Set values of integer work array to allow user convenient
c use of `sol_eval' routine, after a successful termination.
c See the `sol_eval' routine for further details.
if (info .EQ. 0) then
iwork(1) = neqns
iwork(2) = Nsub
iwork(3) = s
iwork(4) = s_star
iwork(5) = method
c
c Set values of double precision work array to allow user
c convenient use of `sol_eval' routine after a successful termination.
call dcopy(Nsub+1,mesh(0),1,work(2*Mxs*neqns*Nsub+1),1)
call dcopy((Nsub+1)*neqns,Y,1,
+ work(2*Mxs*neqns*Nsub+(Nsub+1)+1),1)
c
c Transfer the 1-D version of `Y' back to its 2-D equivalent.
do 45 i = 1,(Nsub-1),-1
do 50 j = 1,neqns
Y(i*ldY + j) = Y(i*neqns + j)
50 continue
45 continue
else
if (profile.GT.print_level_0) then
write(6,*)
write(6,*) 'The computation has failed.'
write(6,*)
write(6,*) '1) If the profiling information begins with and ',
+ 'simply proceeds through a series of failed Newton ',
+ 'iterations each followed by a doubling of the mesh, then',
+ 'there may be an error in the coding of the problem.',
+ 'It is recommended that the routines fsub, gsub, dfsub and',
+ ' dgsub be checked.'
write(6,*)
write(6,*) '2) If the profiling information shows a series',
+ 'of successful, converged Newton iterations with a halt',
+ 'simply because the value of MxNsub was exceeded, then it',
+ 'is recommended that the program be run again with a larger',
+ ' value of MxNsub.'
write(6,*)
write(6,*) '3) More information may be obtained by choosing',
+ 'a larger value for "output_control".'
endif
endif
c
return
end
SUBROUTINE COLROW(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,
* NCLBLK,NBLOKS,BOTBLK,NRWBOT,PIVOT,B,X,IFLAG)
C
C***************************************************************
C
C THIS PROGRAM SOLVES THE LINEAR SYSTEM A*X = B WHERE A IS
C AN ALMOST BLOCK DIAGONAL MATRIX OF THE FORM
C
C TOPBLK
C ARRAY(1)
C ARRAY(2)
C .
C .
C .
C .
C ARRAY(NBLOKS)
C BOTBLK
C
C WHERE
C TOPBLK IS NRWTOP BY NOVRLP
C ARRAY(K), K=1,NBLOKS, ARE NRWBLK BY NRWBLK+NOVRLP
C BOTBLK IS NRWBOT BY NOVRLP,
C AND
C NOVRLP = NRWTOP + NRWBOT
C WITH
C NOVRLP.LE.NRWBLK .
C
C THE LINEAR SYSTEM IS OF ORDER N = NBLOKS*NRWBLK + NOVRLP.
C
C THE METHOD IMPLEMENTED IS BASED ON GAUSS ELIMINATION WITH
C ALTERNATE ROW AND COLUMN ELIMINATION WITH PARTIAL PIVOTING,
C WHICH PRODUCES A STABLE DECOMPOSITION OF THE MATRIX A
C WITHOUT INTRODUCING FILL-IN.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C TO OBTAIN A SINGLE PRECISION VERSION OF THIS PACKAGE, REMOVE
C ALL DOUBLE PRECISION STATEMENTS. THERE IS ONE SUCH STATEMENT
C IN C O L R O W, THREE IN C R D C M P, AND TWO IN C R S O L V.
C IN ADDITION, REFERENCES TO BUILT-IN FUNCTIONS DABS AND DMAX1
C MUST BE REPLACED BY ABS AND AMAX1, RESPECTIVELY. DABS OCCURS
C NINE TIMES, IN C R D C M P. DMAX1 OCCURS FOUR TIMES, IN
C C R D C M P. FINALLY, ZERO IS INITIALISED TO 0.D0 IN A
C DATA STATEMENT IN C R D C M P. THIS MUST BE REPLACED BY:
C DATA ZERO/0.0/
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C ***** PARAMETERS *****
C
C *** ON ENTRY ...
C
C N - INTEGER
C THE ORDER OF THE LINEAR SYSTEM,
C GIVEN BY NBLOKS*NRWBLK + NOVRLP
C
C TOPBLK - DOUBLE PRECISION(NRWTOP,NOVRLP)
C THE FIRST BLOCK OF THE ALMOST BLOCK
C DIAGONAL MATRIX A
C
C NRWTOP - INTEGER
C NUMBER OF ROWS IN THE BLOCK TOPBLK
C
C NOVRLP - INTEGER
C THE NUMBER OF COLUMNS IN WHICH SUCC-
C ESSIVE BLOCKS OVERLAP, WHERE
C NOVRLP = NRWTOP + NRWBOT
C
C ARRAY - DOUBLE PRECISION(NRWBLK,NCLBLK,NBLOKS)
C ARRAY(,,K) CONTAINS THE K-TH NRWBLK
C BY NCLBLK BLOCK OF THE MATRIX A
C
C NRWBLK - INTEGER
C NUMBER OF ROWS IN K-TH BLOCK
C
C NCLBLK - INTEGER
C NUMBER OF COLUMNS IN K-TH BLOCK
C
C NBLOKS - INTEGER
C NUMBER OF NRWBLK BY NCLBLK BLOCKS IN
C THE MATRIX A
C
C BOTBLK - DOUBLE PRECISION(NRWBOT,NOVRLP)
C THE LAST BLOCK OF THE MATRIX A
C
C NRWBOT - INTEGER
C NUMBER OF ROWS IN THE BLOCK BOTBLK
C
C PIVOT - INTEGER(N)
C WORK SPACE
C
C B - DOUBLE PRECISION(N)
C THE RIGHT HAND SIDE VECTOR
C
C X - DOUBLE PRECISION(N)
C WORK SPACE
C
C *** ON RETURN ...
C
C TOPBLK,ARRAY,BOTBLK - ARRAYS CONTAINING THE
C DESIRED DECOMPOSITION OF THE MATRIX A
C (IF IFLAG = 0)
C
C PIVOT - INTEGER(N)
C RECORDS THE PIVOTING INDICES DETER-
C MINED IN THE DECOMPOSITION
C
C X - DOUBLE PRECISION(N)
C THE SOLUTION VECTOR (IF IFLAG = 0)
C
C IFLAG - INTEGER
C = 1, IF INPUT PARAMETERS ARE INVALID
C = -1, IF MATRIX IS SINGULAR
C = 0, OTHERWISE
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C ***** AUXILIARY PROGRAMS *****
C
C CRDCMP(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,NCLBLK,NBLOKS,
C * BOTBLK,NRWBOT,PIVOT,IFLAG)
C - DECOMPOSES THE MATRIX A USING MODIFIED
C ALTERNATE ROW AND COLUMN ELIMINATON WITH
C PARTIAL PIVOTING, AND IS USED FOR THIS
C PURPOSE IN C O L R O W.
C THE ARGUMENTS ARE AS IN C O L R O W.
C
C CRSLVE(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,NCLBLK,NBLOKS,
C * BOTBLK,NRWBOT,PIVOT,B,X)
C - SOLVES THE SYSTEM A*X = B ONCE A IS DECOMPOSED.
C THE ARGUMENTS ARE ALLAS IN C O L R O W.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C THE SUBROUTINE C O L R O W AUTOMATICALLY SOLVES THE
C INPUT SYSTEM WHEN IFLAG=0. C O L R O W IS CALLED ONLY ONCE
C FOR A GIVEN SYSTEM. THE SOLUTION FOR A SEQUENCE OF P RIGHT
C HAND SIDES CAN BE OBTAINED BY ONE CALL TO C O L R O W AND
C P-1 CALLS TO CRSLVE ONLY. SINCE THE ARRAYS TOPBLK,ARRAY,
C BOTBLK AND PIVOT CONTAIN THE DECOMPOSITION OF THE GIVEN
C COEFFICIENT MATRIX AND PIVOTING INFORMATION ON RETURN FROM
C C O L R O W , THEY MUST NOT BE ALTERED BETWEEN SUCCESSIVE
C CALLS TO CRSLVE WITH THE SAME LEFT HAND SIDES. FOR THE
C SAME REASON, IF THE USER WISHES TO SAVE THE COEFFICIENT
C MATRIX, THE ARRAYS TOPBLK,ARRAY,BOTBLK MUST BE COPIED
C BEFORE A CALL TO C O L R O W .
C
C*************************************************************************
C
C ***** SAMPLE CALLING PROGRAM *****
C
C THE FOLLOWING PROGRAM WILL EXERCISE COLROW, IN THE
C CASE WHEN THE COEFFICIENT MATRIX IS NON-SINGULAR.
C
C DOUBLE PRECISION TOP,AR,BOT,B,X
C DOUBLE PRECISION ERROR,ERR
C DIMENSION TOP(2,4),AR(4,8,2),BOT(2,4),B(12),X(12)
C INTEGER PIVOT(12)
C DATA N,NRWTOP,NOVRLP,NRWBLK,NCLBLK,NBLOKS,NRWBOT/12,2,4,4,8,2,2/
C DATA TOP(1,1),TOP(1,2),TOP(1,3),TOP(1,4),
C * TOP(2,1),TOP(2,2),TOP(2,3),TOP(2,4)/
C *0.0 D0,-0.98D0,-0.79D0,-0.15D0,
C *-1.00D0, 0.25D0,-0.87D0, 0.35D0/
C DATA AR(1,1,1),AR(1,2,1),AR(1,3,1),AR(1,4,1),
C * AR(1,5,1),AR(1,6,1),AR(1,7,1),AR(1,8,1)/
C *0.78D0, 0.31D0,-0.85D0, 0.89D0,-0.69D0,-0.98D0,-0.76D0,-0.82D0/
C DATA AR(2,1,1),AR(2,2,1),AR(2,3,1),AR(2,4,1),
C * AR(2,5,1),AR(2,6,1),AR(2,7,1),AR(2,8,1)/
C *0.12D0,-0.01D0, 0.75D0, 0.32D0,-1.00D0,-0.53D0,-0.83D0,-0.98D0/
C DATA AR(3,1,1),AR(3,2,1),AR(3,3,1),AR(3,4,1),
C * AR(3,5,1),AR(3,6,1),AR(3,7,1),AR(3,8,1)/
C *-0.58D0, 0.04D0, 0.87D0, 0.38D0,-1.00D0,-0.21D0,-0.93D0,-0.84D0/
C DATA AR(4,1,1),AR(4,2,1),AR(4,3,1),AR(4,4,1),
C * AR(4,5,1),AR(4,6,1),AR(4,7,1),AR(4,8,1)/
C *-0.21D0,-0.91D0,-0.09D0,-0.62D0,-1.99D0,-1.12D0,-1.21D0, 0.07D0/
C DATA AR(1,1,2),AR(1,2,2),AR(1,3,2),AR(1,4,2),
C * AR(1,5,2),AR(1,6,2),AR(1,7,2),AR(1,8,2)/
C *0.78D0,-0.93D0,-0.76D0, 0.48D0,-0.87D0,-0.14D0,-1.00D0,-0.59D0/
C DATA AR(2,1,2),AR(2,2,2),AR(2,3,2),AR(2,4,2),
C * AR(2,5,2),AR(2,6,2),AR(2,7,2),AR(2,8,2)/
C *-0.99D0, 0.21D0,-0.73D0,-0.48D0,-0.93D0,-0.91D0, 0.10D0,-0.89D0/
C DATA AR(3,1,2),AR(3,2,2),AR(3,3,2),AR(3,4,2),
C * AR(3,5,2),AR(3,6,2),AR(3,7,2),AR(3,8,2)/
C *-0.68D0,-0.09D0,-0.58D0,-0.21D0, 0.85D0,-0.39D0, 0.79D0,-0.71D0/
C DATA AR(4,1,2),AR(4,2,2),AR(4,3,2),AR(4,4,2),
C * AR(4,5,2),AR(4,6,2),AR(4,7,2),AR(4,8,2)/
C *0.39D0,-0.99D0,-0.12D0,-0.75D0,-0.68D0,-0.99D0, 0.50D0,-0.88D0/
C DATA BOT(1,1),BOT(1,2),BOT(1,3),BOT(1,4),
C * BOT(2,1),BOT(2,2),BOT(2,3),BOT(2,4)/
C *0.71D0,-0.64D0, 0.0 D0, 0.48D0,
C *0.08D0,100.0D0,50.00D0,15.00D0/
C DATA B(1),B(2),B(3),B(4),B(5),B(6),B(7),B(8),B(9),B(10),B(11),
C * B(12)/
C *-1.92D0,-1.27D0,-2.12D0,-2.16D0,-2.27D0,-6.08D0,-3.03D0,
C *-4.62D0,-1.02D0,-3.52D0,.55D0,165.08D0/
C
C*************************************************************************
C
C THE INPUT MATRIX IS GIVEN BY:
C
C 0.0 -0.98 -0.79 -0.15
C -1.00 0.25 -0.87 0.35
C 0.78 0.31 -0.85 0.89 -0.69 -0.98 -0.76 -0.82
C 0.12 -0.01 0.75 0.32 -1.00 -0.53 -0.83 -0.98
C -0.58 0.04 0.87 0.38 -1.00 -0.21 -0.93 -0.84
C -0.21 -0.91 -0.09 -0.62 -1.99 -1.12 -1.21 0.07
C 0.78 -0.93 -0.76 0.48 -0.87 -0.14 -1.00 -0.59
C -0.99 0.21 -0.73 -0.48 -0.93 -0.91 0.10 -0.89
C -0.68 -0.09 -0.58 -0.21 0.85 -0.39 0.79 -0.71
C 0.39 -0.99 -0.12 -0.75 -0.68 -0.99 0.50 -0.88
C 0.71 -0.64 0.0 0.48
C 0.08 100.0 50.00 15.00
C
C THE RIGHT HAND SIDE IS GIVEN BY:
C
C B = (-1.92,-1.27,-2.12,-2.16,-2.27,-6.08,-3.03,-4.62,
C -1.02,-3.52,0.55,165.08)
C
C THE SOLUTION OF THIS SYSTEM IS GIVEN BY;
C
C X = (1,1,1,1,1,1,1,1,1,1,1,1)
C
C*************************************************************************
C
C CALL COLROW(N,TOP,NRWTOP,NOVRLP,AR,NRWBLK,NCLBLK,NBLOKS,
C * BOT,NRWBOT,PIVOT,B,X,IFLAG)
C IF(IFLAG.NE.0)GO TO 1000
C ERROR = 0.D0
C DO 10 I=1,N
C ERR = 1.D0 - X(I)
C ERROR = DMAX1(ERROR,DABS(ERR))
C WRITE(6,100)X(I),ERR
C 10 CONTINUE
C WRITE(6,200)ERROR
C 200 FORMAT(12H MAX ERROR = ,D15.7)
C 100 FORMAT(1H ,F15.7,D15.7)
C RETURN
C1000 CONTINUE
C WRITE(6,300)IFLAG
C 300 FORMAT(9H IFLAG = ,I3)
C RETURN
C END
C
C***************************************************************
C
DOUBLE PRECISION TOPBLK,ARRAY,BOTBLK,B,X
INTEGER PIVOT(1)
DIMENSION TOPBLK(NRWTOP,1),ARRAY(NRWBLK,NCLBLK,1),
* BOTBLK(NRWBOT,1),B(1),X(1)
CALL CRDCMP(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,NCLBLK,NBLOKS,
* BOTBLK,NRWBOT,PIVOT,IFLAG)
IF(IFLAG.NE.0)RETURN
CALL CRSLVE(TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,NCLBLK,NBLOKS,
* BOTBLK,NRWBOT,PIVOT,B,X)
RETURN
END
SUBROUTINE CRDCMP(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,
* NCLBLK,NBLOKS,BOTBLK,NRWBOT,PIVOT,IFLAG)
C
C***************************************************************
C
C C R D C M P DECOMPOSES THE ALMOST BLOCK DIAGONAL MATRIX A
C USING MODIFIED ALTERNATE ROW AND COLUMN ELIMINATION WITH
C PARTIAL PIVOTING. THE MATRIX A IS STORED IN THE ARRAYS
C TOPBLK, ARRAY, AND BOTBLK.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C ***** PARAMETERS *****
C
C *** ON ENTRY ...
C
C N - INTEGER
C THE ORDER OF THE LINEAR SYSTEM,
C GIVEN BY NBLOKS*NRWBLK + NOVRLP
C
C TOPBLK - DOUBLE PRECISION(NRWTOP,NOVRLP)
C THE FIRST BLOCK OF THE ALMOST BLOCK
C DIAGONAL MATRIX A TO BE DECOMPOSED
C
C NRWTOP - INTEGER
C NUMBER OF ROWS IN THE BLOCK TOPBLK
C
C NOVRLP - INTEGER
C THE NUMBER OF COLUMNS IN WHICH SUCC-
C ESSIVE BLOCKS OVERLAP, WHERE
C NOVRLP = NRWTOP + NRWBOT
C
C ARRAY - DOUBLE PRECISION(NRWBLK,NCLBLK,NBLOKS)
C ARRAY(,,K) CONTAINS THE K-TH NRWBLK
C BY NCLBLK BLOCK OF THE MATRIX A
C
C NRWBLK - INTEGER
C NUMBER OF ROWS IN K-TH BLOCK
C
C NCLBLK - INTEGER
C NUMBER OF COLUMNS IN K-TH BLOCK
C
C NBLOKS - INTEGER
C NUMBER OF NRWBLK BY NCLBLK BLOCKS IN
C THE MATRIX A
C
C BOTBLK - DOUBLE PRECISION(NRWBOT,NOVRLP)
C THE LAST BLOCK OF THE MATRIX A
C
C NRWBOT - INTEGER
C NUMBER OF ROWS IN THE BLOCK BOTBLK
C
C PIVOT - INTEGER(N)
C WORK SPACE
C
C *** ON RETURN ...
C
C TOPBLK,ARRAY,BOTBLK - ARRAYS CONTAINING THE
C DESIRED DECOMPOSITION OF THE MATRIX A
C (IF IFLAG = 0)
C
C PIVOT - INTEGER(N)
C RECORDS THE PIVOTING INDICES DETER-
C MINED IN THE DECOMPOSITION
C
C IFLAG - INTEGER
C = 1, IF INPUT PARAMETERS ARE INVALID
C = -1, IF MATRIX IS SINGULAR
C = 0, OTHERWISE
C
C***************************************************************
C
DOUBLE PRECISION TOPBLK,ARRAY,BOTBLK
DOUBLE PRECISION ROWMAX,ROWPIV,ROWMLT,COLMAX,COLPIV
DOUBLE PRECISION SWAP,COLMLT,PIVMAX,ZERO,TEMPIV
INTEGER PIVOT(1)
DIMENSION TOPBLK(NRWTOP,1),ARRAY(NRWBLK,NCLBLK,1),
* BOTBLK(NRWBOT,1)
DATA ZERO/0.0D0/
C
C***************************************************************
C
C **** DEFINE THE CONSTANTS USED THROUGHOUT ****
C
C***************************************************************
C
IFLAG = 0
PIVMAX = ZERO
NRWTP1 = NRWTOP+1
NROWEL = NRWBLK-NRWTOP
NRWEL1 = NROWEL+1
NVRLP0 = NOVRLP-1
C
C***************************************************************
C
C **** CHECK VALIDITY OF THE INPUT PARAMETERS....
C
C IF PARAMETERS ARE INVALID THEN TERMINATE AT 10;
C ELSE CONTINUE AT 100.
C
C***************************************************************
C
IF(N.NE.NBLOKS*NRWBLK+NOVRLP)GO TO 10
IF(NOVRLP.NE.NRWTOP+NRWBOT)GO TO 10
IF(NCLBLK.NE.NOVRLP+NRWBLK)GO TO 10
IF(NOVRLP.GT.NRWBLK)GO TO 10
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C PARAMETERS ARE ACCEPTABLE - CONTINUE AT 100.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
GO TO 100
10 CONTINUE
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C PARAMETERS ARE INVALID. SET IFLAG = 1, AND TERMINATE
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
IFLAG = 1
RETURN
100 CONTINUE
C
C***************************************************************
C
C **** FIRST, IN TOPBLK....
C
C***************************************************************
C
C *** APPLY NRWTOP COLUMN ELIMINATIONS WITH COLUMN
C PIVOTING ....
C
C***************************************************************
C
DO 190 I = 1,NRWTOP
IPLUS1 = I+1
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C DETERMINE COLUMN PIVOT AND PIVOT INDEX
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
IPVT = I
COLMAX = DABS(TOPBLK(I,I))
DO 110 J = IPLUS1,NOVRLP
TEMPIV = DABS(TOPBLK(I,J))
IF(TEMPIV.LE.COLMAX)GO TO 110
IPVT = J
COLMAX = TEMPIV
110 CONTINUE
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C TEST FOR SINGULARITY:
C
C IF SINGULAR THEN TERMINATE AT 1000;
C ELSE CONTINUE.
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
IF(PIVMAX+COLMAX.EQ.PIVMAX)GO TO 1000
PIVMAX = DMAX1(COLMAX,PIVMAX)
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C IF NECESSARY INTERCHANGE COLUMNS
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
PIVOT(I) = IPVT
IF(IPVT.EQ.I)GO TO 140
DO 120 L = I,NRWTOP
SWAP = TOPBLK(L,IPVT)
TOPBLK(L,IPVT) = TOPBLK(L,I)
TOPBLK(L,I) = SWAP
120 CONTINUE
DO 130 L = 1,NRWBLK
SWAP = ARRAY(L,IPVT,1)
ARRAY(L,IPVT,1) = ARRAY(L,I,1)
ARRAY(L,I,1) = SWAP
130 CONTINUE
140 CONTINUE
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C COMPUTE MULTIPLIERS AND PERFORM COLUMN
C ELIMINATION
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
COLPIV = TOPBLK(I,I)
DO 180 J = IPLUS1,NOVRLP
COLMLT = TOPBLK(I,J)/COLPIV
TOPBLK(I,J) = COLMLT
IF(IPLUS1.GT.NRWTOP)GO TO 160
DO 150 L = IPLUS1,NRWTOP
TOPBLK(L,J) = TOPBLK(L,J)-COLMLT*TOPBLK(L,I)
150 CONTINUE
160 CONTINUE
DO 170 L = 1,NRWBLK
ARRAY(L,J,1) = ARRAY(L,J,1)-COLMLT*ARRAY(L,I,1)
170 CONTINUE
180 CONTINUE
190 CONTINUE
C
C***************************************************************
C
C **** IN EACH BLOCK ARRAY(,,K)....
C
C***************************************************************
C
INCR = 0
DO 395 K = 1,NBLOKS
KPLUS1 = K+1
C
C *****************************************************
C
C *** FIRST APPLY NRWBLK-NRWTOP ROW ELIMINATIONS WITH
C ROW PIVOTING....
C
C *****************************************************
C
DO 270 J = NRWTP1,NRWBLK
JPLUS1 = J+1
JMINN = J-NRWTOP
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C DETERMINE ROW PIVOT AND PIVOT INDEX
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
IPVT = JMINN
ROWMAX = DABS(ARRAY(JMINN,J,K))
LOOP = JMINN+1
DO 210 I = LOOP,NRWBLK
TEMPIV = DABS(ARRAY(I,J,K))
IF(TEMPIV.LE.ROWMAX)GO TO 210
IPVT = I
ROWMAX = TEMPIV
210 CONTINUE
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C TEST FOR SINGULARITY:
C
C IF SINGULAR THEN TERMINATE AT 1000;
C ELSE CONTINUE.
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
IF(PIVMAX+ROWMAX.EQ.PIVMAX)GO TO 1000
PIVMAX = DMAX1(ROWMAX,PIVMAX)
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C IF NECESSARY INTERCHANGE ROWS
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
INCRJ = INCR+J
PIVOT(INCRJ) = INCR+IPVT+NRWTOP
IF(IPVT.EQ.JMINN)GO TO 230
DO 220 L = J,NCLBLK
SWAP = ARRAY(IPVT,L,K)
ARRAY(IPVT,L,K) = ARRAY(JMINN,L,K)
ARRAY(JMINN,L,K) = SWAP
220 CONTINUE
230 CONTINUE
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C COMPUTE MULTIPLERS
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
ROWPIV = ARRAY(JMINN,J,K)
DO 240 I = LOOP,NRWBLK
ARRAY(I,J,K) = ARRAY(I,J,K)/ROWPIV
240 CONTINUE
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C PERFORM ROW ELIMINATION WITH COLUMN INDEXING
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
DO 260 L = JPLUS1,NCLBLK
ROWMLT = ARRAY(JMINN,L,K)
DO 250 I = LOOP,NRWBLK
ARRAY(I,L,K) = ARRAY(I,L,K)
* -ROWMLT*ARRAY(I,J,K)
250 CONTINUE
260 CONTINUE
270 CONTINUE
C
C *****************************************************
C
C *** NOW APPLY NRWTOP COLUMN ELIMINATIONS WITH
C COLUMN PIVOTING....
C
C *****************************************************
C
DO 390 I = NRWEL1,NRWBLK
IPLUSN = I+NRWTOP
IPLUS1 = I+1
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C DETERMINE COLUMN PIVOT AND PIVOT INDEX
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
IPVT = IPLUSN
COLMAX = DABS(ARRAY(I,IPVT,K))
LOOP = IPLUSN+1
DO 310 J = LOOP,NCLBLK
TEMPIV = DABS(ARRAY(I,J,K))
IF(TEMPIV.LE.COLMAX)GO TO 310
IPVT = J
COLMAX = TEMPIV
310 CONTINUE
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C TEST FOR SINGULARITY:
C
C IF SINGULAR THEN TERMINATE AT 1000;
C ELSE CONTINUE.
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
IF(PIVMAX+COLMAX.EQ.PIVMAX)GO TO 1000
PIVMAX = DMAX1(COLMAX,PIVMAX)
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C IF NECESSARY INTERCHANGE COLUMNS
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
INCRN = INCR+IPLUSN
PIVOT(INCRN) = INCR+IPVT
IRWBLK = IPLUSN-NRWBLK
IF(IPVT.EQ.IPLUSN)GO TO 340
DO 315 L = I,NRWBLK
SWAP = ARRAY(L,IPVT,K)
ARRAY(L,IPVT,K) = ARRAY(L,IPLUSN,K)
ARRAY(L,IPLUSN,K) = SWAP
315 CONTINUE
IPVBLK = IPVT-NRWBLK
IF(K.EQ.NBLOKS)GO TO 330
DO 320 L = 1,NRWBLK
SWAP = ARRAY(L,IPVBLK,KPLUS1)
ARRAY(L,IPVBLK,KPLUS1)
* = ARRAY(L,IRWBLK,KPLUS1)
ARRAY(L,IRWBLK,KPLUS1) = SWAP
320 CONTINUE
GO TO 340
330 CONTINUE
DO 335 L = 1,NRWBOT
SWAP = BOTBLK(L,IPVBLK)
BOTBLK(L,IPVBLK) = BOTBLK(L,IRWBLK)
BOTBLK(L,IRWBLK) = SWAP
335 CONTINUE
340 CONTINUE
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C COMPUTE MULTIPLIERS AND PERFORM COLUMN
C ELIMINATION
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
COLPIV = ARRAY(I,IPLUSN,K)
DO 380 J = LOOP,NCLBLK
COLMLT = ARRAY(I,J,K)/COLPIV
ARRAY(I,J,K) = COLMLT
IF(I.EQ.NRWBLK)GO TO 350
DO 345 L = IPLUS1,NRWBLK
ARRAY(L,J,K) = ARRAY(L,J,K)
* -COLMLT*ARRAY(L,IPLUSN,K)
345 CONTINUE
350 CONTINUE
JRWBLK = J-NRWBLK
IF(K.EQ.NBLOKS)GO TO 370
DO 360 L = 1,NRWBLK
ARRAY(L,JRWBLK,KPLUS1) =
* ARRAY(L,JRWBLK,KPLUS1)
* -COLMLT*ARRAY(L,IRWBLK,KPLUS1)
360 CONTINUE
GO TO 380
370 CONTINUE
DO 375 L = 1,NRWBOT
BOTBLK(L,JRWBLK) = BOTBLK(L,JRWBLK)
* -COLMLT*BOTBLK(L,IRWBLK)
375 CONTINUE
380 CONTINUE
390 CONTINUE
INCR = INCR + NRWBLK
395 CONTINUE
C
C***************************************************************
C
C **** FINALLY, IN BOTBLK....
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C *** APPLY NRWBOT ROW ELIMINATIONS WITH ROW
C PIVOTING....
C
C IF BOT HAS JUST ONE ROW GO TO 500
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
IF(NRWBOT.EQ.1)GO TO 500
DO 470 J = NRWTP1,NVRLP0
JPLUS1 = J+1
JMINN = J-NRWTOP
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C DETERMINE ROW PIVOT AND PIVOT INDEX
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
IPVT = JMINN
ROWMAX = DABS(BOTBLK(JMINN,J))
LOOP = JMINN+1
DO 410 I = LOOP,NRWBOT
TEMPIV = DABS(BOTBLK(I,J))
IF(TEMPIV.LE.ROWMAX) GO TO 410
IPVT = I
ROWMAX = TEMPIV
410 CONTINUE
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C TEST FOR SINGULARITY:
C
C IF SINGULAR THEN TERMINATE AT 1000;
C ELSE CONTINUE.
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
IF(PIVMAX+ROWMAX.EQ.PIVMAX)GO TO 1000
PIVMAX = DMAX1(ROWMAX,PIVMAX)
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C IF NECESSARY INTERCHANGE ROWS
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
INCRJ = INCR+J
PIVOT(INCRJ) = INCR+IPVT+NRWTOP
IF(IPVT.EQ.JMINN)GO TO 430
DO 420 L = J,NOVRLP
SWAP = BOTBLK(IPVT,L)
BOTBLK(IPVT,L) = BOTBLK(JMINN,L)
BOTBLK(JMINN,L) = SWAP
420 CONTINUE
430 CONTINUE
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C COMPUTE MULTIPLIERS
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
ROWPIV = BOTBLK(JMINN,J)
DO 440 I = LOOP,NRWBOT
BOTBLK(I,J) = BOTBLK(I,J)/ROWPIV
440 CONTINUE
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C PERFORM ROW ELIMINATION WITH COLUMN INDEXING
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
DO 460 L = JPLUS1,NOVRLP
ROWMLT = BOTBLK(JMINN,L)
DO 450 I = LOOP,NRWBOT
BOTBLK(I,L) = BOTBLK(I,L)-ROWMLT*BOTBLK(I,J)
450 CONTINUE
460 CONTINUE
470 CONTINUE
500 CONTINUE
C
C***************************************************************
C
C DONE PROVIDED THE LAST ELEMENT IS NOT ZERO
C
C***************************************************************
C
IF(PIVMAX+DABS(BOTBLK(NRWBOT,NOVRLP)).NE.PIVMAX) RETURN
C
C***************************************************************
C
C **** MATRIX IS SINGULAR - SET IFLAG = - 1.
C TERMINATE AT 1000.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
1000 CONTINUE
IFLAG = -1
RETURN
END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine criterion(neqns, leftbc, Nsub, mesh, Y, top,blocks,
+ bot, pivot, PHI, delta, g,k_discrete,work,fsub,gsub)
c
c***overview***
c
c This routine takes as input the current iterate, `Y' and computes
c a related function known as the natural criterion function, `g' = g(Y) =
c 0.5*|inv(J)*PHI(Y)|**2, where PHI(Y) is the residual function evaluated
c at `Y' and J is the Jacobian of the residual function evaluated at `Y',
c or possibly at a previous value of the iterate if a Fixed Jacobian
c iteration is employed. See [Ascher, Mattheij, and Russell, 1988,
c Pg. 335]. The intermediate quantities that arise during the calculation
c of `g', namely, `PHI = PHI(Y)' and `delta = inv(J)*PHI(Y)', are returned
c because they can be useful during the next step of the iteration.
c The mesh of `Nsub' subintervals is stored in `mesh'; `neqns' is the
c number of differential equations and the factored Jacobian is provided
c in the arrays `top', `blocks', `bot', and `pivot'. The number of boundary
c conditions imposed at the left endpoint is `leftbc'. The stages values
c generated by the `resid' routine are stored in the `k_discrete'
c array are passed out of this routine for use by other parts of the code.
c `work' is a work array passed to the `resid' routine and used here for
c temporary storage. The array, `work', is a work array.
c The subroutines, `fsub', and `gsub', are provided by the user to define
c the boundary value ODE and the boundary conditions.
c------------------------------------------------------------------------------
c***declaration of constants***
integer Mxs
parameter(Mxs=10)
c
c `Mxs' Maximum value for the total number of stages of the RK method.
c
c We wish to be able to detect when the value of the criterion
c function will overflow. Thus we set the overflow guard to
c sqrt(8.5d37)=9.22d18.
c
double precision overflow_guard
parameter(overflow_guard=9.22d18)
c
c***declaration of parameters***
c imports:
integer neqns,leftbc,Nsub
double precision mesh(0:Nsub)
double precision Y(neqns*(Nsub+1))
c
double precision top(neqns*leftbc)
double precision blocks(2*neqns**2 * Nsub)
double precision bot(neqns*(neqns-leftbc))
integer pivot(neqns*(Nsub+1))
c
c `neqns' Number of differential equations
c `leftbc' Number of boundary conditions at the left
c end of the problem interval
c `Nsub' Number of subintervals into which the problem
c interval is partitioned
c `mesh' Set of points which partition the problem interval
c `Y' Current discrete approximation to the solution,
c evaluated at each meshpoint.
c `top' Derivative information for left boundary conditions,
c in factored form.
c `blocks' Derivative information for residual function, in factored form.
c `bot' Derivative information for right boundary conditions,
c in factored form.
c `pivot' Contains pivoting information associated with the factorization
c of the Newton matrix. Together `top', `blocks', `bot', and
c `pivot' contain the almost block diagonal matrix representing
c the coefficient matrix of the Newton system in factored form
c subsequent to the call to `COLROW'.
c exports:
double precision PHI(neqns*(Nsub+1)), delta(neqns*(Nsub+1))
double precision g, k_discrete(Mxs*neqns*Nsub)
c
c `PHI' Value of the residual function evaluated at `Y' .
c `delta' Solution to the linear system with the usual Newton matrix as its
c coefficient matrix and the righthand side, `PHI'. That is,
c `delta = inv(J)*PHI', where J is the Jacobian of the residual
c function.
c `g' Value of natural criterion function at `Y'.
c i.e. g = 0.5*|inv(J)*PHI(Y)|^2
c `k_discrete' Vector containing the discrete stages for
c all subintervals. The ith set of `s*neqns'
c locations of `k_discrete' contains the `s'
c stages, each of length `neqns' corresponding
c to the ith subinterval. Passed from `resid'.
c
c***declaration of work space***
double precision work(neqns*(Nsub+1))
c
c `work' Work array passed to `resid' and also used here for
c a temporary copy of the contents of `PHI' provided
c to `crslve' because `crslve' will overwrite the
c contents of the `PHI' vector, which we need to export.
c
c***user-supplied subroutines***
external fsub,gsub
c
c `fsub' Defines f(t,y) for first order system of differential equations.
c `gsub' Defines the boundary condition equations.
c
c***declaration of variables for /IOCNTRL/ common block***
integer print_level_0, print_level_1, print_level_2
parameter (print_level_0 = 0, print_level_1 = 1,
* print_level_2 = 2)
integer profile
common /IOCNTRL/ profile
c
c `print_level_0' No output.
c `print_level_1' Intermediate output.
c `print_level_2' Full output.
c `profile' Controls output, to standard output, of profiling infor-
c mation such as Newton iteration counts, mesh selection,
c relative defect estimate. This variable is identical to
c the `output_control' parameter.
c----------------------------------------------------------------------------
c called by: `fixed_jacob', `damp_factor'
c calls to: `resid', `dcopy', `crslve', `idamax'
integer idamax
c-----------------------------------------------------------------------------
c***evaluate natural criterion function at Y ****
c First we compute PHI(Y) by calling `resid'; the value will
c be returned in the vector `PHI'.
c
call resid(neqns,leftbc,Nsub,mesh,Y,PHI,
+ k_discrete,work,fsub,gsub)
c
c We next compute `delta = inv(J)*PHI(Y)' by solving
c the corresponding linear system . J, stored in `top',
c `blocks', `bot', and `pivot' has already been
c factored so we can simply call `crslve' to get the solution
c corresponding to the right hand side, PHI(Y), stored in `PHI'.
c Since `crslve' would overwrite `PHI', we actually provide
c only a copy of `PHI' to `crslve'.
c
call dcopy((Nsub+1)*neqns,PHI,1,work,1)
c
call crslve(top, leftbc, neqns, blocks,
+ neqns,2*neqns,Nsub,bot,neqns-leftbc,pivot,work,delta)
c
c We complete this calculation by computing the value of the
c natural criterion function at `Y'. In terms of the quantities we
c have just computed, this is `g = 0.5*|delta|**2'.
c
g = dabs(delta(idamax((Nsub+1)*neqns,delta,1)))
c
c Here we insert a test here to guard against overflow. If the
c norm of the Newton correction is large enough that its square
c overflows, the iteration is in difficulty anyway.
c We signal difficulty by setting g to -1.0. The negative g
c condition will be trapped by the calling routine, and the
c iteration will be terminated.
c
if (g .GT. overflow_guard) then
g = -1.0d0
if (PROFILE .GT. PRINT_LEVEL_1) then
write(6,*)'Natural criterion function overflow'
end if
else
g = 0.5d0 * g**2
if (PROFILE .GT. PRINT_LEVEL_2) then
write(6,*)'Natural criterion function value',g
end if
endif
c
return
end
SUBROUTINE CRSLVE(TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,
* NCLBLK,NBLOKS,BOTBLK,NRWBOT,PIVOT,B,X)
C
C***************************************************************
C
C C R S L V E SOLVES THE LINEAR SYSTEM
C A*X = B
C USING THE DECOMPOSITION ALREADY GENERATED IN C R D C M P.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C ***** PARAMETERS *****
C
C *** ON ENTRY ...
C
C TOPBLK - DOUBLE PRECISION(NRWTOP,NOVRLP)
C OUTPUT FROM C R D C M P
C
C NOVRLP - INTEGER
C THE NUMBER OF COLUMNS IN WHICH SUCC-
C ESSIVE BLOCKS OVERLAP, WHERE
C NOVRLP = NRWTOP + NRWBOT
C
C NRWTOP - INTEGER
C NUMBER OF ROWS IN THE BLOCK TOPBLK
C
C ARRAY - DOUBLE PRECISION(NRWBLK,NCLBLK,NBLOKS)
C OUTPUT FROM C R D C M P
C
C NRWBLK - INTEGER
C NUMBER OF ROWS IN K-TH BLOCK
C
C NCLBLK - INTEGER
C NUMBER OF COLUMNS IN K-TH BLOCK
C
C NBLOKS - INTEGER
C NUMBER OF NRWBLK BY NCLBLK BLOCKS IN
C THE MATRIX A
C
C BOTBLK - DOUBLE PRECISION(NRWBOT,NOVRLP)
C OUTPUT FROM C R D C M P
C
C NRWBOT - INTEGER
C NUMBER OF ROWS IN THE BLOCK BOTBLK
C
C PIVOT - INTEGER(N)
C THE PIVOT VECTOR FROM C R D C M P
C
C B - DOUBLE PRECISION(N)
C THE RIGHT HAND SIDE VECTOR
C
C X - DOUBLE PRECISION(N)
C WORK SPACE
C
C *** ON RETURN ...
C
C
C X - DOUBLE PRECISION(N)
C THE SOLUTION VECTOR
C
C***************************************************************
C
DOUBLE PRECISION TOPBLK,ARRAY,BOTBLK,X,B
DOUBLE PRECISION DOTPRD,XJ,XINCRJ,BINCRJ,SWAP
INTEGER PIVOT(1)
DIMENSION TOPBLK(NRWTOP,1),ARRAY(NRWBLK,NCLBLK,1),
* BOTBLK(NRWBOT,1),B(1),X(1)
C
C***************************************************************
C
C **** DEFINE THE CONSTANTS USED THROUGHOUT ****
C
C***************************************************************
C
NRWTP1 = NRWTOP+1
NRWBK1 = NRWBLK+1
NVRLP1 = NOVRLP+1
NRWTP0 = NRWTOP-1
NRWBT1 = NRWBOT+1
NROWEL = NRWBLK-NRWTOP
NRWEL1 = NROWEL+1
NVRLP0 = NOVRLP-1
NBLKS1 = NBLOKS+1
NBKTOP = NRWBLK+NRWTOP
C
C***************************************************************
C
C **** FORWARD RECURSION ****
C
C***************************************************************
C
C *** FIRST, IN TOPBLK....
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C FORWARD SOLUTION
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
DO 130 J = 1,NRWTOP
X(J) = B(J)/TOPBLK(J,J)
IF(J.EQ.NRWTOP)GO TO 120
XJ = -X(J)
LOOP = J+1
DO 110 I = LOOP,NRWTOP
B(I) = B(I)+TOPBLK(I,J)*XJ
110 CONTINUE
120 CONTINUE
130 CONTINUE
C
C ********************************************************
C
C *** IN EACH BLOCK ARRAY(,,K)....
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
INCR = 0
DO 280 K = 1,NBLOKS
INCRTP = INCR+NRWTOP
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C FORWARD MODIFICATION
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
DO 220 J = 1,NRWTOP
INCRJ = INCR+J
XINCRJ = -X(INCRJ)
DO 210 I = 1,NRWBLK
INCRI = INCRTP+I
B(INCRI) = B(INCRI)+ARRAY(I,J,K)*XINCRJ
210 CONTINUE
220 CONTINUE
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C FORWARD ELIMINATION
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
DO 240 J = NRWTP1,NRWBLK
INCRJ = INCR+J
JPIVOT = PIVOT(INCRJ)
IF(JPIVOT.EQ.INCRJ)GO TO 225
SWAP = B(INCRJ)
B(INCRJ) = B(JPIVOT)
B(JPIVOT) = SWAP
225 CONTINUE
BINCRJ = -B(INCRJ)
LOOP = J-NRWTP0
DO 230 I = LOOP,NRWBLK
INCRI = INCRTP+I
B(INCRI) = B(INCRI)+ARRAY(I,J,K)*BINCRJ
230 CONTINUE
240 CONTINUE
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C FORWARD SOLUTION
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
DO 270 J = NRWBK1,NBKTOP
INCRJ = INCR+J
JRWTOP = J -NRWTOP
X(INCRJ) = B(INCRJ)/ARRAY(JRWTOP,J,K)
IF(J.EQ.NBKTOP)GO TO 260
XINCRJ = -X(INCRJ)
LOOP = J-NRWTP0
DO 250 I = LOOP,NRWBLK
INCRI = INCRTP+I
B(INCRI) = B(INCRI)+ARRAY(I,J,K)*XINCRJ
250 CONTINUE
260 CONTINUE
270 CONTINUE
INCR = INCR+NRWBLK
280 CONTINUE
C
C ********************************************************
C
C *** FINALLY, IN BOTBLK....
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C FORWARD MODIFICATION
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
INCRTP = INCR+NRWTOP
DO 320 J = 1,NRWTOP
INCRJ = INCR+J
XINCRJ = -X(INCRJ)
DO 310 I = 1,NRWBOT
INCRI = INCRTP+I
B(INCRI) = B(INCRI)+BOTBLK(I,J)*XINCRJ
310 CONTINUE
320 CONTINUE
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C FORWARD ELIMINATION
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
IF(NRWBOT.EQ.1)GO TO 350
DO 340 J = NRWTP1,NVRLP0
INCRJ = INCR+J
JPIVOT = PIVOT(INCRJ)
IF(JPIVOT.EQ.INCRJ)GO TO 325
SWAP = B(INCRJ)
B(INCRJ) = B(JPIVOT)
B(JPIVOT) = SWAP
325 CONTINUE
BINCRJ = -B(INCRJ)
LOOP = J-NRWTP0
DO 330 I = LOOP,NRWBOT
INCRI = INCRTP+I
B(INCRI) = B(INCRI)+BOTBLK(I,J)*BINCRJ
330 CONTINUE
340 CONTINUE
350 CONTINUE
C
C***************************************************************
C
C **** BACKWARD RECURSION ****
C
C***************************************************************
C
C *** FIRST IN BOTBLK....
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C BACKWARD SOLUTION
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
DO 430 LL = 1,NRWBOT
J = NVRLP1-LL
INCRJ = INCR+J
NRWBTL = NRWBT1-LL
X(INCRJ) = B(INCRJ)/BOTBLK(NRWBTL,J)
IF(LL.EQ.NRWBOT)GO TO 420
XINCRJ = -X(INCRJ)
LOOP = NRWBOT-LL
DO 410 I = 1,LOOP
INCRI = INCRTP+I
B(INCRI) = B(INCRI)+BOTBLK(I,J)*XINCRJ
410 CONTINUE
420 CONTINUE
430 CONTINUE
C
C ********************************************************
C
C *** THEN IN EACH BLOCK ARRAY(,,K)....
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
DO 490 L = 1,NBLOKS
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C BACKWARD ELIMINATION
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
K = NBLKS1-L
INCR = INCR-NRWBLK
DO 450 L1 = NRWEL1,NRWBLK
I = NRWBLK+NRWEL1-L1
IPLUSN = I+NRWTOP
LOOP = IPLUSN+1
INCRN = INCR+IPLUSN
DOTPRD = X(INCRN)
DO 440 J = LOOP,NCLBLK
INCRJ = INCR+J
DOTPRD = DOTPRD-ARRAY(I,J,K)*X(INCRJ)
440 CONTINUE
X(INCRN) = DOTPRD
IPVTN = PIVOT(INCRN)
IF(INCRN.EQ.IPVTN)GO TO 445
SWAP = X(INCRN)
X(INCRN) = X(IPVTN)
X(IPVTN) = SWAP
445 CONTINUE
450 CONTINUE
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C BACKWARD MODIFICATION
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
INCRTP = INCR+NRWTOP
DO 460 J = NRWBK1,NCLBLK
INCRJ = INCR+J
XINCRJ = -X(INCRJ)
DO 455 I = 1,NROWEL
INCRI = INCRTP+I
B(INCRI) = B(INCRI)+ARRAY(I,J,K)*XINCRJ
455 CONTINUE
460 CONTINUE
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C BACKWARD SOLUTION
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
DO 480 LL = 1,NROWEL
J = NRWBK1-LL
INCRJ = INCR+J
NRWELL = NRWEL1-LL
X(INCRJ) = B(INCRJ)/ARRAY(NRWELL,J,K)
IF(LL.EQ.NROWEL)GO TO 470
XINCRJ = -X(INCRJ)
LOOP = NROWEL-LL
DO 465 I = 1,LOOP
INCRI = INCRTP+I
B(INCRI) = B(INCRI)+ARRAY(I,J,K)*XINCRJ
465 CONTINUE
470 CONTINUE
480 CONTINUE
490 CONTINUE
C
C ********************************************************
C
C *** IN TOPBLK FINISH WITH....
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C BACKWARD ELIMINATION
C
C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
DO 520 L = 1,NRWTOP
I = NRWTP1-L
LOOP = I+1
DOTPRD = X(I)
DO 510 J = LOOP,NOVRLP
DOTPRD = DOTPRD-TOPBLK(I,J)*X(J)
510 CONTINUE
X(I) = DOTPRD
IPVTI = PIVOT(I)
IF(I.EQ.IPVTI)GO TO 515
SWAP = X(I)
X(I) = X(IPVTI)
X(IPVTI) = SWAP
515 CONTINUE
520 CONTINUE
RETURN
END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine damp_factor(neqns,leftbc,Nsub,mesh,Y,delta_0,
+ g_0,top,bot,blocks,pivot,lambda_min,lambda,
+ PHI,delta,g,Fixed_Jacobian,info,
+ k_discrete,Y_0,work,fsub,gsub)
c
c***overview***
c
c This routine takes a discrete solution approximation, `Y', provided
c on a mesh of `Nsub' subintervals (contained in `mesh') and uses
c an iteration to produce a suitable value for the damping
c factor `lambda', for use in a damped Newton update of `Y'. The number
c of differential equations is `neqns', while the number of boundary
c conditions at the lefthand endpoint of the problem interval is `leftbc'.
c
c With the iterate `Y' stored in `Y_0' and the corresponding Newton
c correction input in `delta_0', this routine predicts a new value
c for `lambda' and then computes a trial iterate value
c `Y = Y_0 - lambda*delta_0'. Its acceptability is measured by
c whether or not it produces a sufficient reduction in the natural
c criterion function `g = 0.5*|delta|^2', where `delta = inv(J)*PHI',
c where J is evaluated at `Y_0', the previous iterate, and PHI is the
c value of the residual function at the trial iterate `Y'. The arrays
c `top', `bot', `blocks', contain the derivative information associated
c with the Jacobian matrix J and `pivot' contains the pivoting information
c for the linear system solution. The new `g' value is compared with the
c previous value stored in `g_0'. If the reduction is not sufficient, a new
c damping factor is chosen, if possible, and then this evaluation
c step is repeated. If it is found that it is impossible to find a
c suitable damping factor, this is interpreted as a signal that the
c entire Newton iteration will not converge. This is indicated by
c setting `info' = 3. `lambda_min' is a minimum value for the
c damping factor. A successful return is indicated by `info'= 0. A large
c natural criterion function value, `g', signifies difficulty in the
c Newton iteration and is indicated by `info' = 3.
c
c At the end of a successful search for a damping factor, it may be
c appropriate to take one or more fixed_Jacobian steps. The criterion is
c that the current damped Newton step was a full Newton step i.e.
c `lambda' = 1. This is signaled by setting `Fixed_Jacobian' to TRUE.
c A by-product of the call to `resid' is the calculation of
c intermediate quantities associated with the underlying Runge-Kutta method
c called stages. These are returned in the array called `k_discrete'.
c This array is exported through the parameter list for later use
c during the construction and evaluation of the interpolant associated
c with the Runge-Kutta method. The array, `work', is a work array.
c The subroutines, `fsub', and `gsub', are provided by the user to define
c the boundary value ODE and the boundary conditions.
c------------------------------------------------------------------------------
c***declaration of constants***
integer Mxs
parameter (Mxs=10)
c
c `Mxs' Maximum number of stages of the RK method.
c
c***declaration of parameters***
c imports:
integer neqns, leftbc, Nsub
double precision mesh(0:Nsub), Y(neqns*(Nsub+1))
double precision delta_0(neqns*(Nsub+1)), g_0
double precision top(neqns*neqns), bot(neqns*neqns)
double precision blocks(2*neqns**2*Nsub)
integer pivot(neqns*(Nsub+1))
double precision lambda_min, lambda
c
c `neqns' Number of differential equations.
c `leftbc' Number of boundary conditions at the left
c end of the problem interval.
c `Nsub' Number of subintervals into which the problem
c interval is partitioned.
c `mesh' Set of points which partition the problem interval.
c `Y' On input, the current discrete approximation to the solution.
c `delta_0' Full Newton correction, equal to `inv(J)*PHI',
c where J and PHI are evaluated at the input iterate, `Y'.
c `g_0' Value of the natural criterion function at the input iterate.
c `top' Storage of derivative info for left boundary conditions.
c `bot' Storage of derivative information for right boundary conditions.
c `blocks' Storage of derivative information for residual function.
c `pivot' Pivoting information during linear system solution.
c Together `top', `blocks', and `bot' contain
c the almost block diagonal matrix representing the
c coefficient matrix of the Newton system, in factored form.
c `lambda_min' Minimum value for the damping factor. If the
c algorithm prescribes a value for lambda below this
c minimum, it is an indication that the Jacobian is
c "effectively singular" and that the Newton iteration
c will not converge.
c `lambda' On input, initial guess for the value of the damping
c factor used to perform the damped Newton update.
c
c exports:
c double precision Y(neqns*(Nsub+1)), lambda
double precision PHI(neqns*(Nsub+1)), delta(neqns*(Nsub+1)), g
logical Fixed_Jacobian
integer info
double precision k_discrete(Mxs*neqns*Nsub)
c
c `Y' On output, the updated discrete approximation to the
c solution, evaluated at each meshpoint.
c `lambda' On output, the updated damping factor used to perform
c the damped Newton update.
c `PHI' Value of the residual function associated with the updated iterate.
c `delta' Newton correction vector, equal to inv(J)*PHI,
c where PHI is evaluated at the new iterate and J is
c evaluated at the previous iterate.
c `g' Value of natural criterion function at the updated iterate.
c `Fixed_Jacobian' .TRUE. if the next step to be taken should be
c a fixed Jacobian step. .FALSE. if the next step
c to be taken should be a damped Newton step.
c `info' Communication flag:
c 0 - successful iteration converged to within user
c specified tolerance.
c 3 - unsuccessful termination - it was impossible to
c obtain a suitable damping factor for the Newton
c update (indicative of an "effectively singular"
c Jacobian) or evaluation of natural criterion
c function overflowed
c `k_discrete' Discrete stages for all subintervals. The i-th set of
c `s*neqns' locations of `k_discrete' contains the `s' stages,
c each of length `neqns' corresponding to the i-th subinterval.
c Computed by the `resid' routine within the call to
c `criterion'.
c work space:
double precision Y_0(neqns*(Nsub+1))
double precision work(neqns*(Nsub+1))
c `Y_0' Current iterate value is saved here during the search
c for a suitable damping factor and new iterate value.
c `work' Passed to the `criterion' routine.
c
c***user-supplied subroutines***
external fsub,gsub
c
c `fsub' Defines f(t,y) for first order system of differential equations.
c `gsub' Defines the boundary condition equations.
c
c***declaration of local variables***
logical accept
double precision sigma, tau
c
c `accept' Used to monitor the search for a new damping factor.
c `sigma, tau' Parameters used to the control selection
c of the damping factor - see descriptions below.
c
c***declaration of variables for /IOCNTRL/ common block***
c imports:
integer print_level_0,print_level_1,print_level_2
parameter (print_level_0 = 0,print_level_1 = 1,
+ print_level_2 = 2)
integer profile
common /IOCNTRL/ profile
c
c `print_level_0' No output.
c `print_level_1' Intermediate output.
c `print_level_2' Full output.
c `profile' Controls output, to standard output, of profiling infor-
c mation such as Newton iteration counts, mesh selection,
c relative defect estimate. This variable is identical to
c the `output_control' parameter.
c------------------------------------------------------------------------------
c Called by: `damped_newt'
c Calls to : `criterion', `dcopy', `daxpy'
c------------------------------------------------------------------------------
c***initialization***
c Two parameters which control the search for a suitable damping
c factor are defined here:
c
c (a) `sigma' ensures that the reduction in the size of the natural
c criterion function `g', evaluated at the new iterate, will be
c nonnegligible. This keeps the iteration from stalling.
c
c (b) `tau' controls how much the values for `lambda' are allowed to
c change from one step to the next in the iteration performed by
c this routine.
sigma = 0.01d0
tau = 0.1d0
c
c During the calculation of trial iterates for the determination
c of a suitable damping factor, `Y' will be overwritten, so
c save the current iterate, stored in `Y', in the vector, `Y_0'.
call dcopy((Nsub+1)*neqns,Y,1,Y_0,1)
c
c***iterative selection of new damping factor***
c
c Use a repeat-until loop to control the search for the new
c damping factor. The loop will be controlled by the flags,
c `accept' and `info': `accept' = .TRUE. or info .NE. 0
c will terminate the loop. (Implement the loop in FORTRAN
c using an if-goto pair.)
c
accept = .FALSE.
c REPEAT
2 continue
c
c ***compute trial iterate using current damping factor***
c
c Compute an updated iterate value, stored in `Y', using the
c current value for the damping factor, stored in `lambda'
c `Y = Y_0 - lambda*delta_0'.
c
call dcopy((Nsub+1)*neqns,Y_0,1,Y,1)
call daxpy((Nsub+1)*neqns,-lambda,delta_0,1,Y,1)
c
c ***evaluate natural criterion function at this new iterate***
c
c Measure the suitability of the current value of the
c damping factor by seeing how much the trial iterate
c reduces the value of the natural criterion function. Compute
c the value of the criterion function at `Y' by calling the
c `criterion' routine. The value is returned in `g'; the
c residual function, evaluated at `Y' is returned in `PHI'.
c The value `inv(J)*PHI', (with `J' evaluated at the previous
c iterate), is returned in `delta'.
c
call criterion(neqns,leftbc,Nsub,mesh,Y,top,
+ blocks, bot, pivot,PHI,delta,g,k_discrete,work,fsub,gsub)
c
c If the natural criterion function value, `g', has overflowed
c the `criterion' routine will trap that condition and set
c `g = -1.0' ( a value which cannot otherwise arise).
c Since the natural criterion function is related to the
c size of the Newton correction, a large `g' value signifies
c difficulty in the Newton iteration. Set `info'=3 and exit.
c
if (g .LT. 0.0d0) then
info = 3
else
c (g >= 0.0d0)
c Natural criterion function has not overflowed.
c
c ***accept current damping factor or compute a new one,
c depending on criterion function value***
c
c Check value of `g' to ensure that the reduction in size is
c sufficient. If so `lambda' and thus `Y' are acceptable.
c If not, select a new `lambda', if possible.
c
if (g .LE. (1.0d0 - 2.0d0*lambda*sigma)*g_0) then
accept = .TRUE.
c
c This damped Newton step is successful. The new iterate
c value has already been stored in `Y' so no update is
c needed at this point.
c
c Decide whether or not to take a fixed Jacobian step
c next time.
c
if (lambda .EQ. 1.0d0) Fixed_Jacobian = .TRUE.
else
c `(g > (1.0d0-2.0d0*lambda*sigma)*g_0')
c Compute a new damping factor.
lambda = max(tau*lambda,
+ (lambda**2 * g_0)/((2*lambda-1.0d0)*g_0+g))
c
if (profile .GT. print_level_2) then
write(6,*)'Damping factor lambda=', lambda
endif
c
c Check this new `lambda' value. If it is too small, then
c the Newton iteration will not converge, set `info' = 3
c to signal this.
c
if (lambda .LT. lambda_min) info = 3
endif
c End of `if (g.LE.(1.0d0-2.0d0*lambda*sigma)*g_0)'
end if
c End of `if (g.LT.0.0d0)'
c
c UNTIL `(accept .OR. (info .GT. 0))'
c
if ((.NOT. accept) .AND.
+ (info .LE. 0)) go to 2
c
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine damped_newt(neqns,leftbc,Nsub,mesh,Y,newtol,
+ lambda,PHI,top,bot,blocks,pivot,Fixed_Jacobian,
+ convrg,delta,delta_0_norm,g,info,k_discrete,delta_0,
+ work,fsub,gsub,dfsub,dgsub)
c
c***overview***
c
c This routine takes a mesh of `Nsub' subintervals, partitioning
c the problem interval and a discrete initial approximation provided
c at the meshpoints, (contained in the vectors `Y' and `mesh',
c respectively) and uses a full or damped Newton step to refine a
c solution approximation for the associated discrete nonlinear system.
c Convergence occurs if the Newton correction is less than the
c tolerance `newtol' and is indicated by `convrg' = TRUE. The
c number of differential equations is `neqns', while the number of
c boundary conditions at the lefthand endpoint of the problem interval
c is `leftbc'. A successful return is indicated by `info'= 0. If this
c routine encounters a singular matrix, it signals this be setting
c `info' = 2. A large natural criterion function value, `g', signifies
c difficulty and is indicated by `info' = 3. On the other hand, if the
c iteration appears to be converging well, this routine will signal that
c the iteration can switch over to a fixed Jacobian step (see below).
c
c At the beginning of each damped Newton step, it is assumed that
c the current iterate value, `Y', and the value of the residual function,
c `PHI', evaluated at `Y' are available. If the previous step was also
c a damped Newton step, this routine also expects `delta = inv(J)*PHI',
c where J, the Jacobian, is evaluated at the previous iterate and
c `delta_0_norm', the norm of the previous value of `delta_0 = inv(J)*PHI',
c where both J and `PHI' are associated with the previous iterate. The
c step begins with the evaluation of then Jacobian of the residual function
c at the current iterate, `Y'. The Newton system, consisting of the
c Jacobian as the coefficient matrix and `PHI', as the right hand side, is
c then solved to obtain the full Newton correction, which we store in the
c vector `delta_0'. The factored Jacobian is available in the vectors
c `top', `bot', `blocks'. Pivoting information for the solution of the
c linear system is stored in `pivot'. During the call to `damp_factor' the
c current iterate is copied into `Y_0', a damping factor, `lambda' is
c computed, and then a damped Newton step is taken to produce a new
c iterate, `Y = Y_0 - lambda*delta_0'.
c
c At the end of a successful damped Newton step, it may be
c appropriate to take one or more fixed Jacobian steps . The criterion
c is that the current damped Newton step was a full Newton step
c (`lambda'=1). This is signaled by setting `Fixed_Jacobian' to TRUE.
c The fixed Jacobian step, for reasons associated with efficient data
c management, expects to receive a number of quantities : (i) the new
c iterate, `Y', (ii) the residual function evaluated at this new iterate,
c `PHI', (iii) `delta = inv(J)*PHI', where J, the Jacobian, has been
c evaluated rather at the previous iterate and (iv) the value of the
c natural criterion function, `g = 0.5*|delta|**2'. All these values are
c to be naturally computed during the current damped Newton step and
c thus will be already set up for the fixed Jacobian step if appropriate.
c
c A by-product of the evaluation of the right hand side of the
c Newton system is the calculation of intermediate quantities associated
c with the underlying Runge-Kutta method called stages. These are
c returned to the `damped_newt' routine when it calls the `damp_factor'
c routine and are stored in the array called `k_discrete'.
c This array is exported through the parameter list for later use
c during the construction and evaluation of the interpolant associated
c with the Runge-Kutta method. The array, `work', is a work array provided
c to the `damp_factor' and `newmat' routines. The subroutines, `fsub',
c `gsub', `dfsub', and `dgsub', are provided by the user to define the
c boundary value ODE and the boundary conditions and their derivatives.
c------------------------------------------------------------------------------
c***declaration of constants***
integer Mxs
parameter (Mxs=10)
c
c `Mxs' is the maximum number of stages of the RK method.
c
c We wish to be able to detect when the criterion function value
c will overflow. Thus the overflow_guard to 9.22d18.
c
double precision overflow_guard
parameter (overflow_guard = 9.22d18)
c
c***declaration of parameters***
c imports:
integer neqns, leftbc, Nsub
double precision mesh(0:Nsub), Y(neqns*(Nsub+1))
double precision newtol, lambda, PHI(neqns*(Nsub+1))
double precision delta(neqns*(Nsub+1)), delta_0_norm
c
c `neqns' Number of differential equations.
c `leftbc' Number of boundary conditions at the left
c end of the problem interval.
c `Nsub' Number of subintervals into which the problem
c interval is partitioned.
c `mesh' Set of points which partition the problem interval.
c `Y' On input, the current discrete approximation to the solution.
c `newtol' Value to be applied to the Newton iteration,
c i.e. the iteration stops if the Newton correction is
c less than or equal to the user specified tolerance.
c A blended relative/absolute tolerance is applied.
c If `delta_j' is the Newton correction for, `Y_j', the jth
c component of `Y', then we require
c `(|delta_j| / |Y_j|) <= newtol'.
c `lambda' On input, the damping factor used in the previous damped
c Newton step. If its value is set to 0 on input, the
c previous step was not a damped Newton step.
c `PHI' On input, the value of the residual function on each
c subinterval evaluated at the current iterate `Y'.
c `delta' On input, it is `inv(J)*PHI' where J is evaluated at the
c previous iterate and `PHI' is evaluated at the input
c iterate, if the previous step was a damped Newton step.
c `delta_0_norm' On input, is |inv(J)*PHI| where J and `PHI' were
c evaluated at the previous iterate if the
c previous step was a damped Newton step.
c It is undefined if the previous step was not a
c damped Newton step.
c exports:
c double precision Y(neqns*(Nsub+1)), lambda
c double precision PHI(neqns*(Nsub+1))
double precision top(neqns*neqns)
double precision bot(neqns*neqns)
double precision blocks(2*neqns**2*Nsub)
integer pivot(neqns*(Nsub+1))
logical Fixed_Jacobian, convrg
c double precision delta(neqns*(Nsub+1))
c double precision delta_0_norm
double precision g
integer info
double precision k_discrete(Mxs*neqns*Nsub)
c
c `Y' On output, the updated discrete approximation to the solution.
c `lambda' On output, the value for the new damping factor.
c `PHI' On output, the value of the residual function on each
c subinterval evaluated at the updated iterate `Y'.
c `top' Storage of derivative information for left boundary
c conditions.
c `bot' Storage of derivative information for right boundary
c conditions.
c `blocks' Storage of derivative information for residual function.
c `pivot' Pivoting information during linear system solution.
c Together `top', `blocks', and `bot' contain
c the almost block diagonal matrix representing the
c coefficient matrix of the Newton system, in factored form.
c `Fixed_Jacobian' .TRUE. if the next step to be taken should be
c a fixed Jacobian step. .FALSE. if the next step
c to be taken should be a damped Newton step.
c `convrg' Used to check for convergence of the Newton iteration.
c `delta' Newton correction vector, equal to `inv(J)*PHI(Y)',
c where J is evaluated at the input iterate and `PHI'
c and `PHI' is evaluated at the new iterate `Y'.
c `delta_0_norm' on output, |inv(J)*PHI| where J and `PHI' are
c evaluated at the input iterate.
c `g' Value of natural criterion function at the updated iterate.
c `info' Internal communication flag:
c 0 - successful iteration converged to within user
c specified tolerance.
c 2 - unsuccessful termination - a singular coefficient
c matrix was encountered during the attempted solution
c of the Newton system.
c 3 - unsuccessful termination - it was impossible to obtain
c a suitable damping factor for the Newton update
c (indicative of an "effectively singular" Jacobian) or
c evaluation of natural criterion function overflowed.
c `k_discrete' Vector containing the discrete stages for all
c subintervals. The i-th set of `s*neqns' locations of
c `k_discrete' contains the `s' stages, each of length
c `neqns' corresponding to the i-th subinterval.
c work space:
double precision delta_0(neqns*(Nsub+1))
double precision
+ work(neqns*(2*Nsub+3+2*(Mxs+1)*neqns))
integer i_Y_0, i_work
c `delta_0' Full Newton correction, equal to `inv(J)*PHI',
c where J and `PHI' are evaluated at the input iterate.
c `work' Work array provided to the `damp_factor' and
c `newmat' routines.
c `i_Y_0' Index to the start of the part of `work' to be passed
c to `damp_factor' for use as the work array `Y_0'.
c `i_work' Index to the start of the remainder of `work' to be passed
c to `damp_factor' for use as the work array `work'.
c
c***user-supplied subroutines***
external fsub,gsub,dfsub,dgsub
c
c `fsub' Defines f(t,y) for the first order
c system of differential equations, y' = f(t,y).
c `gsub' Defines the boundary condition equations.
c `dfsub' Defines the Jacobian of the system
c of differential equations.
c `dgsub' Defines the Jacobian of the boundary conditions.
c
c***declaration of local variables***
integer factor
double precision g_0, lambda_min,mu
integer j
c
c `factor' Used by the `colrow' routine to indicate trouble; a
c return with `factor' equal to -1 indicates that a singular
c Jacobian was encountered.
c `g_0' Value of the natural criterion function at the input iterate.
c `lambda_min' Parameter used to the control selection of the
c damping factor - see descriptions below.
c `mu' Predicted value for the new damping factor.
c `j' Loop index, 1,...,`(Nsub+1)*neqns'.
c
c***declaration of variables for /IOCNTRL/ common block***
c imports:
integer print_level_0,print_level_1,print_level_2
parameter (print_level_0 = 0,print_level_1 = 1,
+ print_level_2 = 2)
integer profile
common /IOCNTRL/ profile
c
c `print_level_0' No output.
c `print_level_1' Intermediate output.
c `print_level_2' Full output.
c `profile' Controls output, to standard output, of profiling infor-
c mation such as Newton iteration counts, mesh selection,
c relative defect estimate. This variable is identical to
c the `output_control' parameter.
c------------------------------------------------------------------------------
c Called by : `NewIter'
c Calls to : `NewMat', `colrow', `damp_factor', `daxpy', `idamax'
integer idamax
c------------------------------------------------------------------------------
c***initialization of work array indexes***
i_Y_0 = 1
i_work = neqns*(Nsub+1)+1
c***perform a damped newton step
c The flag `info' is used to monitor for trouble; It is initialized to 0.
info = 0
c
c `lambda_min' is a parameter which controls the search for a
c suitable damping factor; it gives the minimum value
c for the damping factor. It is used here in the initial guess for the
c damping factor. It is also used with the `damp_factor' routine.
lambda_min = 0.01d0
c
c First construct the Newton system at the input iterate.
c It is assumed that the value of the residual function
c at the current iterate was obtained at the end of the previous
c iteration step and is available in the vector `PHI'. The
c coefficient matrix is computed by calling the `NewMat'
c routine and is returned in the vectors `top', `blocks', `bot'.
c
call NewMat(neqns,leftbc,Nsub,mesh,Y,top,blocks,bot,
+ k_discrete,work,dfsub,dgsub)
c
c The coefficient matrix constructed by `NewMat' has an
c "almost block diagonal" structure. `colrow', a package of
c subroutines designed to solve such systems is now employed.
c This routine will return with `factor' = -1 if the coefficient
c matrix of the Newton system is singular and `factor'=0 otherwise.
c The solution of the linear system, which gives the Newton correction,
c is returned in the vector `delta_0'.
c
call colrow((Nsub+1)*neqns,top,leftbc,neqns,blocks,
+ neqns,2*neqns,Nsub,bot,neqns-leftbc,
+ pivot,PHI,delta_0,factor)
c
if ( factor .EQ. -1 ) then
info = 2
c This is a signal that a singular Jacobian matrix was detected.
end if
c
c
if (info .EQ. 0) then
c
c `colrow' was successful in solving the linear system to obtain
c the Newton correction. Attempt to compute the corresponding
c new iterate using a damped Newton step. At this time, the value
c of the natural criterion function for the current iterate is
c also computed and stored in `g_0'. (It can be shown
c that, for this case, `g_0 = 0.5 * |delta_0|**2'.)
c
c Check here to see if the natural criterion function
c value will overflow. If so, trap the condition. Since the
c natural criterion function is related to the size of the
c Newton correction, a large `g' value signifies difficulty in
c the Newton iteration. Set `info' = 3 and exit.
c
g_0 = dabs(delta_0(idamax((Nsub+1)*neqns, delta_0,1)))
c
if (g_0 .GT. overflow_guard) then
if (profile .GT. print_level_1) then
write(6,*)'Natural criterion function overflow'
info = 3
end if
else
if (profile .GT. print_level_1) then
write(6,*)'Norm of Newton correction',g_0
endif
g_0 = 0.5d0 * g_0**2
endif
c
c Test for convergence now. If the norm of the scaled Newton
c correction, `delta_0', is sufficiently small, then the iteration
c has converged.
c
convrg = .TRUE.
do 20 j = 1, (Nsub+1)*neqns
if( dabs(delta_0(j)) .GT. newtol*(dabs(Y(j))+1.0d0))
+ convrg = .FALSE.
20 continue
c
if (convrg) then
c Update iterate with full Newton correction.
c The current iterate is still available in `Y' and
c the Newton correction is in `delta_0'.
call daxpy((Nsub+1)*neqns,-1.0d0,delta_0,1,Y,1)
c
c The stages of the RK method are updated in the next call to
c the `interp_setup' routine with this new iterate value.
c
else
c Attempt to update the Newton iterate with a damped Newton
c correction. Compute a prediction for the new damping factor.
c The value for `lambda' must be between `lambda_min' and 1.0d0.
c If the input `lambda' value is 0.0d0, this is a signal that the
c previous iteration was not a damped Newton iteration and
c appropriate information for the prediction of a new `lambda'
c value is not available, in which case set the predicted
c value to 1.0d0.
c
if (lambda .EQ. 0.0d0 ) then
lambda = 1.0d0
else
c For temporary storage overwrite `delta = delta - delta_0'.
call daxpy((Nsub+1)*neqns,-1.0d0,delta_0,1,delta,1)
c
c Use the norm of this difference in the prediction for
c `lambda'.
mu = (lambda*delta_0_norm)/
+ dabs(delta(idamax((Nsub+1)*neqns,delta,1)))
lambda = max(lambda_min, min(1.0d0, mu))
endif
c
c Print out current value for `lambda', if appropriate.
if (profile .GT. print_level_1) then
write(6,*)'Initial damping factor lambda =', lambda
endif
c
c The previous value of the norm of the Newton correction,
c stored in `delta_0_norm', was just used in the prediction
c of the new `lambda' value, and is no longer needed.
c Update the value of `delta_0_norm' using the current value
c of `delta_0' computed by `colrow' above (the norm of `delta_0'
c was already computed during the calculation of `g_0').
c
delta_0_norm = dsqrt(2.0d0*g_0)
c
c Compute a new iterate by updating the current iterate with
c the (possibly damped) Newton correction.
c
call damp_factor(neqns,leftbc,Nsub,mesh,Y,delta_0,
+ g_0,top,bot,blocks,pivot,lambda_min,lambda,
+ PHI,delta,g,Fixed_Jacobian,info,k_discrete,
+ work(i_Y_0),work(i_work),fsub,gsub)
if ((info .EQ. 0) .AND.
+ (profile .GT. print_level_2)) then
write(6,*)'Norm of damped Newton correction',
+ lambda*delta_0_norm
endif
c
c If successful, this routine will compute an appropriate
c damping factor, returned in `lambda', the updated iterate,
c returned in `Y', the corresponding residual function value,
c returned in `PHI', `delta = inv(J)*PHI', and `g', the natural
c criterion function, evaluated at the new iterate.
c
endif
c End of `if (convrg)'
c
endif
c End of `if (info.EQ.0)'
c
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine defect_estimate(neqns,Nsub,mesh,Y,defect,
+ defect_norm,info,z,z_prime,def_1,def_2,f_1,f_2,
+ temp_1,temp_2,k_discrete,k_interp,work,fsub)
c
c***overview***
c
c This routine uses the discrete solution approximation, `Y',
c obtained on the mesh stored in `mesh', plus the stages of the underlying
c RK method, stored in `k_discrete', plus some new stages, stored in
c `k_interp', to construct an interpolant for `Y'. `neqns' is the number of
c differential equations; `Nsub' is the the number of subintervals into
c which the mesh divides the problem interval. This interpolant is then used
c to compute an estimate of the defect on each subinterval, for each
c solution component. The i-th group of `neqns' locations of `defect' is
c used to store the relative defect information associated with the i-th
c subinterval. The max norm of the relative defect is returned in
c `defect_norm'. A successful return is indicated by `info = 0'. If the
c `defect_norm' is bigger than the constant parameter `defect_threshold'
c then `info' is set to 4.
c
c The primary information required to define the interpolant on the i-th
c subinterval is the current solution approximation at the (i-1)-st and
c i-th mesh points (which can be obtained from `Y'), and some
c stages associated with the i-th interval. The first `s' stages come
c from the discrete formula, and have already been computed. They are
c available in the `s' vector locations making up the i-th set of
c `s*neqns' locations of the vector,'k_discrete'. The remaining
c `s_star - s' stages are required only for the interpolant and must
c now be computed . For each subinterval, this is done through a call
c to the `interp_setup' routine, which returns the new stage evaluations
c for that subinterval. The resultant stages are then stored in the
c appropriate locations of the `k_interp' vector (which has the same
c format as the `k_discrete' vector).
c
c On the i-th subinterval the evaluation of the interpolant at
c `t = mesh(i-1)+tau*hi' is given by `z(t) = yim1 + hi sum br(tau)*kr',
c where `hi' is the subinterval size, 'yim1' is the solution approximation
c at `mesh(i-1)', `kr' is the r-th stage, `br(tau)' is the corresponding
c weight polynomial evaluated at `tau' and the sum runs from 1 to `s_star'.
c
c Once the interpolant is available, we can sample
c the defect at the point, `mesh(i-1)+tau*hi', by substituting the
c interpolant into the differential equation and subtracting the right
c side from the left side, i.e. the defect on the i-th subinterval is
c `delta = z'(mesh(i-1)+tau*hi)-f(mesh(i-1)+tau*hi,z(mesh(i-1)+tau*hi))'.
c A blended relative/absolute defect measure, `defect', equal to
c `|delta|/(|f|+1)' is computed, where `f' is the corresponding right hand
c side of the ODE.
c
c The sample point is available from the /sample_point/ common
c block and the subsequent evaluation of the weights is done by a call
c to the `interp_weights' routine. To improve the reliability of the
c defect estimate, we also sample at `1-tau' and choose the larger of
c the two defect values for the estimate, `delta'.
c
c `z','z_prime,'def_1','def_2','f_1','f_2','temp_1','temp_2'
c are work arrays used for temporary storage here while `work'
c is a work array passed to `interp_setup'. The subroutine, `fsub',
c is provided by the user to define the boundary value ODEs.
c------------------------------------------------------------------------------
c***declaration of constants***
integer Mxs
parameter(Mxs=10)
c
c `Mxs' Maximum value for the total number of stages of the interpolant.
c
double precision defect_threshold
parameter (defect_threshold = 1.0d-1)
c
c `defect_threshold' Used to monitor the size of the defect.
c
c***declaration of parameters***
c imports:
integer neqns, Nsub
double precision mesh(0:Nsub), Y(neqns*(Nsub+1))
double precision k_discrete(Mxs*neqns*Nsub)
c
c `neqns' Number of differential equations.
c `Nsub' Number of subintervals.
c `mesh' Array of points defining the subintervals.
c `Y' Contains the discrete solution at the mesh points
c the i-th set of `neqns' locations of `Y' contains
c the solution approximation at the i-th mesh point.
c `k_discrete' Contains the `s' discrete stage values for each subinterval.
c
c exports:
double precision defect(Nsub*neqns), defect_norm
integer info
double precision k_interp(Mxs*neqns*Nsub)
c
c `defect' Contains the corresponding blended relative/absolute
c defect measure.
c `defect_norm' Max norm of the defect measure.
c `info' Communication flag :
c 0 - Successful return.
c 4 - Size of the defect norm was too large and the
c solution from the current discrete system can not
c be trusted.
c `k_interp' Contains the `s_star-s' extra stage values for each
c subinterval need to form the interpolant.
c
c work space:
double precision z(neqns),z_prime(neqns)
double precision def_1(neqns),def_2(neqns)
double precision f_1(neqns),f_2(neqns)
double precision temp_1(neqns),temp_2(neqns)
double precision work(neqns)
c `z' Used to accumulate argument for fsub evaluation
c `z_prime' Used to save derivative of interpolant
c `def_1,def_2' Contain the values of the defect at the sample
c points `tau' and `1-tau' for each subinterval.
c `f_1,f_2' Contain the corresponding fsub value at the
c the two sample points.
c `temp_1', `temp_2' Temporary location for storage of scaled defect
c at `tau' and `1-tau'.
c `work' Work array passed to `interp_setup'.
c
c***user-supplied subroutines***
external fsub
c
c `fsub' Defines f(t,y) for the first order system of differential
c equations, y' = f(t,y).
c
c***declaration of local variables***
integer i,j
double precision hi, weights_1(Mxs), weights_2(Mxs)
double precision weights_1_prime(Mxs), weights_2_prime(Mxs)
double precision estimate_1, estimate_2
c
c `i' Loop index from 1 to Nsub.
c `j' Loop index from 1 to `neqns'.
c `hi' Size of the i-th subinterval.
c `f' Used to store `fsub' evaluation.
c `weights_1', `weights_2' Contains the weight polynomials of the
c interpolant evaluated at the sample point `tau'
c and `1-tau'.
c `weights_1_prime', `weights_2_prime' Contains the first derivatives of
c the weight polynomials of the
c interpolant evaluated at the sample
c points `tau' and `1-tau'.
c `estimate_1', `estimate_2' Maximum scaled defect at `tau' and `1-tau'
c
c***declaration of variables from /RK_s/ common block***
integer s
common /RK_s/ s
c `s' Number of discrete stages.
c
c***declaration of variables from /interp_s_star/ common block***
integer s_star
common /interp_s_star/ s_star
c `s_star' Total number of stages required to form the
c interpolant on each subinterval. It includes all the
c stages of the discrete formula plus the additional
c stages required for the interpolant.
c
c***declaration of variables from /sample_point/ common block***
double precision tau_star
common /sample_point/ tau_star
c
c `tau_star' Relative position of one of the sample points within each
c subinterval. The other sample point is `1-tau_star'.
c-----------------------------------------------------------------------------
c Called by: `mirkdc'
c Calls to : `interp_weights', `interp_setup', `sum_stages', `fsub'
c `daxpy', `dcopy'
integer idamax
c------------------------------------------------------------------------------
c Setup the weights, evaluated at the first sample point.
c
call interp_weights(s_star,tau_star,weights_1,weights_1_prime)
c
c Setup the weights, evaluated at the second sample point.
c
call interp_weights(s_star,1.0d0-tau_star,weights_2,
+ weights_2_prime)
c
do 30 i = 1, Nsub
hi = mesh(i) - mesh(i-1)
c
c We must first setup the new stages for this subinterval.
c
call interp_setup(neqns,mesh(i-1),hi,Y((i-1)*neqns+1),
+ Y(i*neqns+1),s,k_discrete((i-1)*s*neqns+1),
+ s_star,k_interp((i-1)*(s_star-s)*neqns+1),work,fsub)
c
c The new stages for the i-th subinterval have been computed
c and stored in the ith set of `(s_star-s)*neqns' locations of
c `k_interp'.
c
c ***Sample Point 1***
c
c With all the information needed to form the interpolant on the
c i-th subinterval already available, we can compute the value of
c the interpolant, and its first derivative at the sample point
c by calling the `sum_stages' routine.
c
call sum_stages(neqns,hi,Y((i-1)*neqns+1),
+ s,k_discrete((i-1)*s*neqns+1),s_star,
+ k_interp((i-1)*(s_star-s)*neqns+1),
+ weights_1,weights_1_prime,z,z_prime)
c
c This yields the value of the interpolant, in `z', and its first
c derivative, in `z_prime'. We must next evaluate the righthand
c side of the system of ode's at `z', i.e. compute f(z) and store
c in `f_1'.
c
call fsub(neqns,mesh(i-1)+tau_star*hi,z,f_1)
c
c Since the calculation of `z_prime' was already complete after
c the return from `sum_stages', we can now compute the defect which
c we store in `def_1'.
c
call daxpy(neqns,-1.0d0,f_1,1,z_prime,1)
call dcopy(neqns,z_prime,1,def_1,1)
c
c Compute the norm of the scaled defect for Sample Point 1.
c
do 10 j =1,neqns
temp_1(j) =
+ def_1(j)/(dabs(f_1(j))+1.0d0)
10 continue
c
estimate_1=dabs(temp_1(idamax(neqns,temp_1,1)))
c
c ***Sample Point 2***
c
c With all the information needed to form the interpolant on the
c i-th subinterval already available, we can compute the value of
c the interpolant, and its first derivative at the sample point
c by calling the `sum_stages' routine.
c
call sum_stages(neqns,hi,Y((i-1)*neqns+1),
+ s,k_discrete((i-1)*s*neqns+1),s_star,
+ k_interp((i-1)*(s_star-s)*neqns+1),
+ weights_2,weights_2_prime,z,z_prime)
c
c This yields the value of the interpolant, in `z', and its first
c derivative, in `z_prime'. We must evaluate the righthand side
c of the system of ode's at `z', i.e. compute f(z).
c
call fsub(neqns,mesh(i-1)+(1.0d0-tau_star)*hi,z,f_2)
c
c Since the calculation of `z_prime' was already complete after
c the return from `sum_stages', we can now compute the defect,
c which we store in `def_2'.
c
call daxpy(neqns,-1.0d0,f_2,1,z_prime,1)
call dcopy(neqns,z_prime,1,def_2,1)
c
do 20 j =1,neqns
temp_2(j) =
+ def_2(j)/(dabs(f_2(j))+1.0d0)
20 continue
c
estimate_2=dabs(temp_2(idamax(neqns,temp_2,1)))
c
c ***Compare defect estimates from the two sample points***
c
if (estimate_1 .GT. estimate_2) then
call dcopy(neqns,temp_1,1,defect((i-1)*neqns+1),1)
else
call dcopy(neqns,temp_2,1,defect((i-1)*neqns+1),1)
endif
c
c The scaled defect for the i-th subinterval has now been copied
c into the ith group of `neqns' locations of the vector `defect'.
c
30 continue
c
c Compute the max norm of the defect estimate.
defect_norm = dabs(defect(idamax(Nsub*neqns,defect,1)))
c
c If this defect_norm is bigger than `defect_threshold' then
c this is an indication that, the defect is no longer a
c small relative perturbation of the original system of ODEs
c and the solution from the current discrete system should
c not be trusted. It is probably better to start over as if the
c Newton iteration had failed. Here, we check the size of `defect_norm',
c and set `info' = 4 (to cause a mesh doubling), if it is too large.
if (defect_norm .GT. defect_threshold) info = 4
c
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine fixed_jacob(neqns,leftbc,Nsub,mesh,Y,newtol,
+ delta,g,PHI,top,bot,blocks,pivot,Fixed_Jacobian,lambda,
+ convrg,info,k_discrete,Y_hat,PHI_hat,work,fsub,gsub)
c
c***overview***
c
c This routine takes a mesh of `Nsub' subintervals, partitioning
c the problem interval and a discrete initial approximation provided
c at the meshpoints, (contained in the vectors `Y' and `mesh',
c respectively) and uses a Fixed Jacobian step, as part of a Newton
c iteration, to produce an updated solution approximation. Convergence
c occurs if the Newton correction is less than the tolerance `newtol'
c and is indicated by `convrg' = TRUE. The number of differential equations
c is `neqns', while the number of boundary conditions at the lefthand
c endpoint of the problem interval is `leftbc'. A successful return is
c indicated by `info'= 0. A large natural criterion function value,'g',
c signifies difficulty in the Newton iteration and is indicated by
c `info' = 3.
c
c When a fixed Jacobian step is to be taken, the following
c quantities are assumed to be available from the previous iteration
c step (i) the iterate, `Y', (ii) the residual function evaluated at this
c iterate, `PHI', (iii) `delta = inv(J)*PHI', where J, the Jacobian ,
c available in factored form in the arrays `top', `bot', `blocks', and
c `pivot', has not been evaluated at `Y', but rather at some previous
c iterate and (iv) the value of the natural criterion function, `g', at
c the iterate, `Y'.
c
c The fixed Jacobian step consists of first computing a (potential)
c new iterate, `Y_hat = Y - delta'. The acceptability of `Y_hat' is
c determined by evaluating the natural criterion function at `Y_hat',
c giving `g_hat'. During the calculation of `g_hat, it is necessary to
c evaluate the residual function at `Y_hat' and store it in the vector,
c `PHI_hat', and a new value for `delta = inv(J)*PHI_hat', both of which
c will then be available for the next step of the iteration. The
c acceptability of the new iterate value is determined by comparing `g_hat'
c with the value of the natural criterion function for the previous iterate,
c stored in `g'. If `g_hat < rho*g' then the new iterate value is accepted
c and another fixed Jacobian step is attempted. If `g_hat > rho*g', the
c code returns to taking damped Newton steps (the value of `rho' is set in
c this routine). This is signaled by setting `Fixed_Jacobian' to FALSE.
c The value of `lambda', the damping factor, is set to 0 as well. This is
c to indicate appropriate information for the prediction of a new damping
c factor is not available for the upcoming damped Newton step.
c
c A by-product of the evaluation of the right hand side of the
c Newton system is the calculation of intermediate quantities associated
c with the underlying Runge-Kutta method called stages. These are
c returned in the array called `k_discrete', which is exported through
c the parameter list for later use during the construction and evaluation
c of the interpolant associated with the Runge-Kutta method. The array,
c `work', is a work array. The subroutines, `fsub', and `gsub', are provided
c by the user to define the boundary value ODE and the boundary conditions.
c------------------------------------------------------------------------------
c***declaration of constants***
integer Mxs
parameter (Mxs=10)
c
c `Mxs' is the maximum number of stages of the RK method.
c
c***declaration of parameters***
c imports:
integer neqns, leftbc, Nsub
double precision mesh(0:Nsub), Y(neqns*(Nsub+1))
double precision newtol, delta(neqns*(Nsub+1)), g
double precision PHI(neqns*(Nsub+1)), top(neqns*neqns)
double precision bot(neqns*neqns), blocks(2*neqns**2*Nsub)
integer pivot(neqns*(Nsub+1))
c
c `neqns' the number of differential equations
c `leftbc' the number of boundary conditions at the left
c end of the problem interval
c `Nsub' the number of subintervals into which the problem
c interval is partitioned
c `mesh' the set of points which partition the problem interval.
c `Y' on input, the current discrete approximation to the solution.
c `newtol' the tolerance value to be applied to the Newton iteration,
c i.e. the iteration stops if the Newton correction
c is less than or equal to the user specified
c tolerance. A blended relative/absolute tolerance
c is applied. If `delta_j' is the Newton correction
c for the jth component of Y, Y_j, then we require
c |delta_j| / (|Y_j|+1) <= newtol.
c `delta' on input, the Newton correction vector, equal to
c inv(J)*PHI(Y), where J is evaluated at a previous iterate.
c `g' on input, the value of natural criterion function at the
c current iterate, `Y'.
c `PHI' on input, is the value of the residual function on each
c subinterval, evaluated at `Y'.
c `top' storage of derivative information for left boundary conditions.
c `bot' storage of derivative information for right boundary conditions.
c `blocks' storage of derivative information for residual function.
c `pivot' pivoting information associated with the factorization
c of J. Together `top', `blocks', and `bot' contain
c the almost block diagonal matrix representing the
c coefficient matrix of the Newton system, in factored form.
c exports:
c double precision Y(neqns*(Nsub+1)),delta(neqns*(Nsub+1))
c double precision g, PHI(neqns*(Nsub+1))
logical Fixed_Jacobian
double precision lambda
logical convrg
integer info
double precision k_discrete(Mxs*neqns*Nsub)
c
c `Y' On output, the updated discrete approximation to the solution.
c `delta' On output, equal to inv(J)*PHI where the
c Jacobian J is evaluated at the previous iterate
c and PHI is evaluated at the updated iterate `Y_hat'
c `g' On output, value of the criterion function evaluated
c at the updated iterate.
c `PHI' On output, is value of the residual function
c on each subinterval evaluated at the updated iterate.
c `Fixed_Jacobian' .TRUE. if the next step to be taken should be a
c fixed Jacobian step. .FALSE. if the next step
c to be taken should be a damped Newton step.
c `lambda' Set to 0 to indicate appropriate information for the
c prediction of a new lambda is not available during the
c upcoming damped Newton step
c `convrg' Used to check for convergence of the Newton iteration.
c `info' Communication flag:
c 0 - successful iteration converged to
c within user specified tolerance.
c 3 - unsuccessful termination - A large natural
c criterion function value, g, signifies
c difficulty in the Newton iteration.
c `k_discrete' Discrete stages for all subintervals. The ith set of
c `s*neqns' locations of `k_discrete' contains the `s'
c stages, each of length `neqns', corresponding
c to the ith subinterval. Computed by `resid'
c from within `criterion`.
c work space:
double precision Y_hat(neqns*(Nsub+1))
double precision PHI_hat(neqns*(Nsub+1))
double precision work(neqns*(Nsub+1))
c `Y_hat' Trial iterate value.
c `PHI_hat' Residual function value associated with the trial iterate.
c `work' Work space passed to `criterion'.
c
c***user - supplied subroutines***
external fsub,gsub
c `fsub' Defines f(t,y) for the first order
c system of differential equations, y' = f(t,y).
c `gsub' Defines the boundary condition equations.
c
c***declaration of local variables***
double precision g_hat, rho
integer j
c
c `g_hat' Value of the natural criterion function at the trial iterate.
c `rho' Used in measuring performance of fixed Jacobian iteration.
c `j' Loop index, 1,...,(Nsub+1)*neqns.
c
c***declaration of variables for /IOCNTRL/ common block***
c imports:
integer print_level_0,print_level_1,print_level_2
parameter (print_level_0 = 0,print_level_1 = 1,
+ print_level_2 = 2)
integer profile
common /IOCNTRL/ profile
c
c `print_level_0' No output.
c `print_level_1' Intermediate output.
c `print_level_2' Full output.
c `profile' Controls output, to standard output, of profiling infor-
c mation such as Newton iteration counts, mesh selection,
c relative defect estimate. This variable is identical to
c the `output_control' parameter.
c------------------------------------------------------------------------------
c Called by : `Newiter'
c Calls to : `criterion', `dcopy', `daxpy', `idamax'
integer idamax
c------------------------------------------------------------------------------
c The flag `info' is used to monitor for trouble; it is initialized to 0.
info = 0
c
c When a fixed Jacobian step has been taken, the decision about
c whether or not to follow with another fixed Jacobian step or not
c will be determined by the size of the decrease in the value of the
c natural criterion function. This decrease must be at least a factor
c of `rho'. Here `rho' is initialized.
rho = 0.5d0
c
c In addition to the new iterate `Y' being available, we also
c have the corresponding evaluation of the residual function,
c `PHI', the Newton correction, `delta', for the current fixed
c Jacobian step, and the evaluation of the criterion function
c for the this iterate, `g = 0.5*|inv(J)*PHI|**2 = 0.5*|delta|**2'.
c
c If the Newton correction, `delta', is less than the tolerance, the
c iteration has converged. A blended relative/absolute tolerance is
c employed.
c
if (profile .GT. print_level_1) then
write(6,*)'Norm of Newton correction',
+ dabs(delta(idamax((Nsub+1)*neqns,delta,1)))
endif
c
convrg = .TRUE.
do 10 j = 1, (Nsub+1)*neqns
if(dabs(delta(j)).GT.newtol*(dabs(Y(j))+1.0d0)) convrg=.FALSE.
10 continue
c
if (convrg) then
c Update current iterate before exiting.
call daxpy((Nsub+1)*neqns,-1.0d0,delta,1,Y,1)
c
c The stages of the RK method are updated in the next
c call to the `interp_setup' routine with this new
c iterate value.
c
else
c If the iteration has not converged, some care must be taken
c in using the fixed Newton step `delta' to update the
c current iterate. The correction will not be useful if it
c does not lead to an iterate, `Y', which provides a sufficient
c reduction in the size of the natural criterion function.
c The acceptability of the new iterate is determined by
c examining the resultant reduction in the size of the natural
c criterion function. The following calls compute the trial iterate
c `Y_hat = Y - delta'.
c
call dcopy((Nsub+1)*neqns,Y,1,Y_hat,1)
call daxpy((Nsub+1)*neqns,-1.0d0,delta,1,Y_hat,1)
c
c The natural criterion function is now evaluated to provide a
c measure of how good this trial iterate is. By-products
c of this calculation, are the quantities `PHI_hat', the evaluation
c of the residual function at `Y_hat', and `delta = 0.5*
c |inv(J)*PHI_hat|' for the NEXT iterate. These MAY be useful for the
c next iteration step. If `Y_hat' turns out to be not acceptable `PHI',
c the residual function at `Y', will be needed, so here it is not
c overwritten. Also, `g', the value of the natural criterion function
c at `Y', will be needed for comparison with `g_hat', the value of
c the natural criterion function at `Y_hat', so `g' is not overwritten
c either. The previous value of `delta = inv(J)*PHI' is no longer
c needed in either case, so it is overwritten with the new value
c `inv(J)*PHI_hat'.
c
call criterion(neqns,leftbc,Nsub,mesh,Y_hat,top,blocks,bot,
+ pivot,PHI_hat,delta,g_hat,k_discrete,work,fsub,gsub)
c
c If the natural criterion function value is about to overflow
c the `criterion' routine will trap that condition and set
c g_hat = -1.0 ( a value which cannot otherwise arise).
c Since the natural criterion function is related to the
c size of the Newton correction, a large value signifies
c difficulty in the Newton iteration. Set `info'=3 and exit
c the loop.
c
if (g_hat .LT. 0.0d0) then
info = 3
else
c The Natural criterion function value has not overflowed.
c If the new iterate results in a sufficient reduction in
c the size of the natural criterion function, accept it.
c
if ( g_hat .LE. rho * g ) then
c The new iterate is acceptable.
c Update the appropriate values `Y','PHI','g' in
c preparation for another fixed Jacobian step.
c (The flag `Fixed_Jacobian' remains TRUE.)
c
call dcopy((Nsub+1)*neqns,Y_hat,1,Y,1)
call dcopy((Nsub+1)*neqns,PHI_hat,1,PHI,1)
g = g_hat
else
c The new iterate is not acceptable. Set the flag so that
c next iteration step is a damped Newton step. This also
c requires setting `lambda=0'.
Fixed_Jacobian = .FALSE.
lambda = 0.0d0
c
c If `g_hat' is greater than `rho' * `g' but less than `g'
c itself, it may be reasonable to use the new iterate
c and corresponding residual function value to begin the
c subsequent full Newton step. These new values, `Y_hat'
c and `PHI_hat', are copied to `Y' and `PHI'.
c
if (g_hat .LT. g) then
call dcopy((Nsub+1)*neqns,Y_hat,1,Y,1)
call dcopy((Nsub+1)*neqns,PHI_hat,1,PHI,1)
else
c Begin the full Newton step using the values in
c `Y' and `PHI' from the beginning of this step.
c A call to `resid' is made here to obtain the `k_discrete'
c values, evaluated at `Y', as these must be up-to-date
c with the `Y' and `PHI' values before the damped Newton
c step begins.
c
call resid(neqns,leftbc,Nsub,mesh,Y,PHI,
+ k_discrete,work,fsub,gsub)
endif
c `if (g_hat .LT. g)'
endif
c `if (g_hat .LE. rho*g)'
endif
c `if (g_hat .LT. 0.0d0)'
endif
c `if (convrg)'
c
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine half_mesh(Nsub,mesh_current,mesh_new)
c
c***overview***
c
c This routine takes as input a mesh of `Nsub' subintervals which
c partitions the problem interval, stored in the vector `mesh_current',
c The routine returns a new mesh, stored in `mesh_new', which is exactly
c twice as fine as the given one, i.e. has `2*Nsub' subintervals, obtained
c by splitting each subinterval of the given mesh in half.
c------------------------------------------------------------------------------
c***declaration of parameters***
c imports:
integer Nsub
double precision mesh_current(0:Nsub)
c
c `Nsub' Number of subintervals in the input mesh.
c `mesh_current' Input meshpoints.
c
c exports:
double precision mesh_new(0:2*Nsub)
c
c `mesh_new' Meshpoints making up the new mesh.
c
c***declaration of local variables***
integer i
c
c `i' loop index from 1 to `Nsub'.
c
c------------------------------------------------------------------------------
c Called by: `mirkdc', `mesh_selector'.
c------------------------------------------------------------------------------
c
mesh_new(0) = mesh_current(0)
c
do 10 i = 1, Nsub
c
mesh_new(2*i) = mesh_current(i)
c
c Generate a new mesh point at location 2*i-1 using the average of
c the adjacent meshpoint values.
c
mesh_new(2*i-1) = (mesh_current(i) + mesh_current(i-1))/2.0d0
c
10 continue
c
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine interp_eval(neqns,Nsub,mesh,Y,t,z,z_prime,
+ k_discrete,k_interp)
c***overview***
c
c This routine uses the current discrete solution approximation, `Y',
c obtained on a mesh of `Nsub' subintervals, stored in `mesh', plus
c additional stage information associated with the underlying discrete RK
c method, to compute an evaluation of an interpolant to the discrete
c solution. This interpolant is evaluated at `t', with the value returned
c in `z', and the value of the first derivative returned in `z_prime'.
c
c The interpolant is computed as follows: Given `t', assumed to
c be somewhere in the problem interval, the subinterval of `mesh' that
c contains `t' is determined first, i.e. we find i such that
c `mesh(i-1) =< t < mesh(i) for t < mesh(Nsub)' or
c `i = Nsub for t = mesh(Nsub)'.
c The fraction of the way `t' is through the i-th subinterval, a quantity
c which we call `tau', is also determined during this step. Then the
c interpolant on this i-th subinterval is evaluated - the interpolant has
c the form `z(mesh(i-1) + tau*hi) = yim1 + hi sum br(tau)*kr'.
c The weights, `br(tau)', can be computed with a call to `interp_weights'.
c
c The remaining information required to define the interpolant
c on the i-th subinterval is the current solution approximation at the
c (i-1)st mesh point, `yim1', which can be obtained from `Y', and the
c stages associated with the i-th subinterval. The first `s' stages,
c contained in `k_discrete', come from the discrete formula and
c remaining `s_star - s' stages, contained in `k_interp', are necessary
c for the interpolant. Each stage is a vector of length `neqns' where
c `neqns' is the number of differential equations.
c
c `k_discrete' is the vector containing the discrete stages for
c all subintervals. The i-th set of `s*neqns' locations of `k_discrete'
c contains the `s' stages, each of length `neqns' corresponding to the
c i-th subinterval. `k_interp' plays a similar role for the `s_star-s'
c stages associated with the interpolant.
c------------------------------------------------------------------------------
c***declaration of constants***
integer Mxs
parameter(Mxs=10)
c
c `Mxs' Maximum value for the total number of stages of the interpolant.
c
c***declaration of parameters***
c imports:
integer neqns, Nsub
double precision mesh(0:Nsub), Y(neqns*(Nsub+1)), t
double precision k_discrete(Mxs*neqns*Nsub)
double precision k_interp(Mxs*neqns*Nsub)
c
c `neqns' Number of differential equations.
c `Nsub' Number of subintervals.
c `mesh' Array of points defining the subintervals.
c `Y' Contains the discrete solution at the mesh points;
c the i-th set of `neqns' locations of `Y' contains
c the solution approximation at the i-th mesh point.
c `t' Point at which the interpolant is to be evaluated.
c `k_discrete' Contains the `s' discrete stage values for each subinterval.
c `k_interp' Contains the `s_star-s' extra stage values for each
c subinterval need to form the interpolant.
c exports:
double precision z(neqns), z_prime(neqns)
c
c `z' Value of the interpolant at `t'.
c `z_prime' Value of the first derivative of the interpolant at `t'.
c
c***declaration of local variables***
integer i
double precision hi
double precision tau, weights(Mxs), weights_prime(Mxs)
c
c `i' Index for subinterval containing `t', i ranges from 1 to `Nsub'.
c `hi' Size of the i-th subinterval.
c `tau' Equal to (t-mesh(i-1))/hi.
c `weights' Contains the weight polynomials of the interpolant
c evaluated at `tau'.
c `weights_prime' Contains the first derivatives of the weight
c polynomials of the interpolant evaluated at `tau'.
c
c***declaration of variables from /RK_s/ common block***
integer s
common /RK_s/ s
c `s' Number of discrete stages.
c
c***declaration of variables from /interp_s_star/ common block***
integer s_star
common /interp_s_star/ s_star
c `s_star' Total number of stages required to form the
c interpolant on each subinterval. It includes all the
c stages of the discrete formula plus the additional
c stages required for the interpolant.
c-----------------------------------------------------------------------------
c Called by: `mirkdc'.
c Calls to: `interval', `interp_weights', `sum_stages'.
c------------------------------------------------------------------------------
c Call `interval' to determine the subinterval containing t.
c The index of the subinterval is returned in `i'.
c
call interval(Nsub,mesh,t,i)
c
c Compute `tau' for use in the evaluation of the weight polynomials.
c
hi = mesh(i) - mesh(i-1)
tau = (t - mesh(i-1))/hi
c
c Setup the weights, evaluated at `tau'.
c
call interp_weights(s_star, tau, weights,weights_prime)
c
c We need to evaluate the interpolant and its derivative at `tau'.
c This is done by taking weighted sums of the stages, using the weights
c we have just computed, and the stages corresponding to the i-th
c subinterval, and calling the `sum_stages' routine.
c
call sum_stages(neqns,hi,Y((i-1)*neqns+1),s,
+ k_discrete((i-1)*s*neqns+1),s_star,
+ k_interp((i-1)*(s_star-s)*neqns+1),
+ weights, weights_prime, z, z_prime)
c
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine interp_setup(neqns,tim1,hi,yim1,yi,s,ki_discrete,
+ s_star,ki_interp,y,fsub)
c***overview***
c
c The purpose of this routine is to perform some preliminary
c calculations to set up stages for later use in the evaluation of the
c interpolant. The interpolant is defined on the i-th subinterval by,
c `Z_i(tim1 + hi*tau) = yim1 + hi sum b_r(tau)*k_r',
c where the stages, `k_r', depend on `yim1' and `yi', the solution
c approximations at the left and right endpoints, `hi', the size of the
c subinterval, 'tim1' the left endpoint, and the coefficients of the
c underlying Runge-Kutta method. The weight polynomials, `b_r(tau)', can be
c evaluated through a call to the `interp_weights' routine.
c The first `s' stages of each subinterval were already obtained during the
c computation of the discrete solution, so only the last `s_star - s' stages
c need to be computed now.
c
c This routine will construct the `s_star-s' extra stages needed by
c the interpolant for a single subinterval. This routine takes as input
c the vector, `ki_discrete', containing the `s' stages associated with
c the discrete formula on the i-th subinterval, and then uses the
c Runge-Kutta coefficients for the extra stages, available from the
c `/interp_coeffs/' common block, to construct the `s_star' new stages.
c These are returned in the vector `ki_interp'. Each stage is a vector of
c length `neqns', where `neqns' is the number of differential equations.
c `y' is a work array. The subroutine, `fsub', is provided by the user to
c define the boundary value ODE.
c-------------------------------------------------------------------------------
c***declaration of constants***
integer Mxs
parameter (Mxs=10)
c
c `Mxs' Maximum value for the total number of stages of the interpolant.
c
c***declaration of parameters***
c imports:
integer neqns
double precision tim1, hi, yim1(neqns), yi(neqns)
integer s, s_star
double precision ki_discrete(s*neqns)
c
c `neqns' Number of ordinary differential equations.
c `tim1' Left endpoint of the ith subinterval.
c `hi' Size of the current subinterval.
c `yim1' Solution approximation at the left endpoint of the current
c subinterval.
c `yi' Solution approximation at the right endpoint of the current
c subinterval.
c `s' Number of discrete stages.
c `s_star' Total number of stages required to form the
c interpolant on each subinterval. It includes all the
c stages of the discrete formula plus the additional
c stages required for the interpolant.
c `ki_discrete' Vector of `s*neqns' components. The j-th
c component is of length `neqns' and corresponds
c to the j-th of the `s' stage values associated
c with the discrete formula on the i-th subinterval.
c exports:
double precision ki_interp((s_star-s)*neqns)
c
c `ki_interp' Vector with same structure as `ki_discrete'.
c The j-th vector component corresponds to the j-th
c of the new stages needed to form the interpolant.
c
c work space:
double precision y(neqns)
c
c `y' Used to hold the argument for the function evaluation
c for each new stage as it is being computed.
c
c***user-supplied subroutines***
external fsub
c
c `fsub' Defines f(t,y) for the first order system
c of differential equations, y' = f(t,y).
c
c***declaration of local variables***
integer r, j
c `r,j' Loop indexes from 1 to `s_star-s'.
c
c***declaration of variables from common block /interp_coeff/***
c imports:
double precision c_star(Mxs), v_star(Mxs), x_star(Mxs*Mxs)
common /interp_coeff/ c_star,v_star,x_star
c
c `c_star' First component contains the abscissa for the first new stage.
c The second component contains the abscissa for the second new
c stage, and so on.
c
c `v_star' First component contains the blending coefficient for the first
c new stage. The second component contains the blending coefficient
c for the second new stage, and so on.
c
c `x_star' Contains the coupling coefficients for each of the new stages.
c The Runge-Kutta coefficient matrix, X_star, is related to the
c data structure, `x_star' as follows (the matrix is stored
c columnwise in the vector):
c
c X_star(i,j) = x_star( (j-1)*s + i )
c
c------------------------------------------------------------------------------
c Called by: `defect_estimate'.
c Calls to : `fsub', `dcopy', `daxpy'.
c------------------------------------------------------------------------------
c***construction of the new stages***
c
do 20 r = 1 , s_star-s
c
c initialize `y' to zero
call dcopy(neqns,0.0d0,0,y,1)
c
c Accumulate the argument for the rth new stage in `y'.
c
c The contributions from the discrete stages.
do 5 j = 1 , s
c
c The contribution from the j-th discrete stage.
call daxpy(neqns,x_star((j-1)*(s_star-s)+r),
+ ki_discrete((j-1)*neqns+1),1,y,1)
5 continue
c
c The contributions from the previous new stages.
c
do 10 j = 1 , r-1
c
c The contribution from the j-th new stage.
call daxpy(neqns,x_star((j+s-1)*(s_star-s)+r),
+ ki_interp((j-1)*neqns+1),1,y,1)
10 continue
c
c Multiply the sum of the stages by `hi'.
call dscal(neqns,hi,y,1)
c
c Add on the contribution from `yim1'.
call daxpy(neqns,(1-v_star(r)),yim1,1,y,1)
c
c Add on the contribution from `yi'.
call daxpy(neqns,v_star(r),yi,1,y,1)
c
c All contributions to the argument for r-th new stage have
c been made. Simply evaluate `fsub' at this argument to get the
c new stage. The result is copied into the appropriate locations
c of `ki_interp'.
c
call fsub(neqns,tim1+c_star(r)*hi,y,
+ ki_interp((r-1)*neqns+1))
20 continue
c
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine interp_tableau (method)
c
c***overview***
c
c This routine defines all the extra coefficients needed to augment
c the discrete RK formula with an interpolant. `method' defines the MIRK
c method to be used. The possible values for `method' and the corresponding
c MIRK scheme associated with each of these values are listed below.
c
c Value of `method' | MIRK formula
c -------------------------------------
c 121 | MIRK121 scheme
c 221 | MIRK221 scheme
c 232 | MIRK232 scheme (abscissa = 0,2/3)
c 2321 | MIRK232 scheme (abscissa = 1,1/3)
c 343 | MIRK343 scheme
c 453 | MIRK453 scheme
c 563 | MIRK563 scheme
c 5631 | An improved MIRK563 scheme
c
c The interpolant is defined on a subinterval by subinterval basis, as
c follows: When t can be expressed as t_i + hi*tau for some i and tau,
c with hi = t_i+1 - t_i then the interpolant, Z(t) = Zi(t_i+ hi*tau) =
c yi + hi sum br(tau)*kr, where the kr depend on information from the ith
c subinterval as well as certain blending and coupling coefficients that
c are specific to the MIRK formula and interpolant. br(tau) is the rth
c weight polynomial of the interpolant evaluated at tau.
c
c The main purpose of this routine is to define the parameter which
c gives,`s_star', the total number of stages for the interpolant and the
c coefficients for the `s_star-s' extra stages needed by the interpolant.
c Once defined by this routine, the coefficients are then
c available from the /interp_coeff/ common block, while `s_star' is
c available from the /interp_s_star/ common block. The `order' of the
c interpolant is available through the /interp_order/ common block.
c
c This routine also defines the sample point, relative to the size
c of a subinterval, at which the defect should be sampled. This is the
c point at which the interpolant must be evaluated and is thus specific
c to the particular interpolant being implemented. The defect sample
c point, which we call `tau_star', is made available through the
c /sample_point/ common block.
c
c The final portion of formula dependent information involves
c the expressions for the weight polynomials, arising as part of the
c definition of the interpolant for the discrete formula. This information
c is provided through the `interp_weights' routine which evaluates all
c the weight polynomials and their derivatives at a given point. This
c routine can be called by other parts of the code whenever the weight
c polynomials and derivatives are required - i.e. whenever an evaluation
c of the interpolant is required.
c-----------------------------------------------------------------------------
c***declaration of constants***
integer Mxs
parameter (Mxs=10)
c
c `Mxs' is the maximum number of stages of the RK method.
c
c***declaration of parameters***
c imports:
integer method
c
c `method' defines the method to be used
c
c***declaration of variables from common block /interp_s_star/***
integer s_star
common /interp_s_star/ s_star
c
c `s_star' is the total number of stages required to form the
c interpolant on each subinterval. It includes all the
c stages of the discrete formula plus the additional
c stages required for the interpolant.
c
c***declaration of variables from common block /interp_coeff/***
c exports:
double precision c_star(Mxs), v_star(Mxs), x_star(Mxs*Mxs)
common /interp_coeff/ c_star,v_star,x_star
c
c The first component of `c_star' contains the abscissa
c for the first new stage. The second component of `c_star' contains
c the abscissa for the second new stage, and so on.
c
c The first component of `v_star' contains the blending coefficient
c for the first new stage. The second component of `v_star' contains
c the blending coefficient for the second new stage, and so on.
c
c The matrix, X_star, is related to the data structure, `x_star'
c as follows (the matrix is stored columnwise in the vector):
c
c X_star(i,j) = x_star( (j-1)*s + i )
c
c `x_star' contains the coupling coefficients for each of the new
c stages.
c
c***declaration of variables for the /sample_point/ common block***
c exports:
double precision tau_star
common /sample_point/ tau_star
c
c `tau_star' is the point within each subinterval at which the
c defect is to be sampled.
c
c***declaration of variables for the /rk_method/ common block***
c exports:
integer amethod
common /rk_method/ amethod
c
c `amethod' the same as `method', defines the method to be
c used. It is defined twice here so that it may be
c passed into a common block.
c
c***declaration of variables for the /interp_order/ common block***
c exports:
integer p
common /interp_order/ p
c
c `p' is the order of the interpolant => local error is O(h**(p+1)).
c------------------------------------------------------------------------------
c called by: `mirkdc'
c calls to: `dcopy'
c------------------------------------------------------------------------------
c Translate short form of method values.
if (method .EQ. 2) method = 221
if (method .EQ. 4) method = 343
if (method .EQ. 6) method = 563
c Assign the value of `method' to `amethod' so that it may be
c passed into a common block.
amethod = method
c Beginning method definition based on value of `method'.
c
if (method. EQ. 121) then
c
c SETUP OF THE COEFFICIENTS FOR THE NEW STAGES
c
c Here we give coefficients defining the interpolant for a 1-stage,
c 2nd order MIRK, stage order 1, with a 2nd order interpolant.
c
s_star = 3
c
c_star(1) = 0.0d0
c_star(2) = 1.0d0
c
v_star(1) = 0.0d0
v_star(2) = 1.0d0
c
c Column 1 of x_star
c
x_star(1) = 0.0d0
x_star(2) = 0.0d0
c
c Column 2 of x_star
c
x_star(3) = 0.0d0
x_star(4) = 0.0d0
c
c Column 3 of x_star
c
x_star(5) = 0.0d0
x_star(6) = 0.0d0
c
c Define the sample point for the interpolant of the 1-stage, 2nd
c order MIRK discrete formula.
c
tau_star = 0.25d0
c
c Define the value of the order of the interpolant, `p'.
c
p = 2
c
elseif (method. EQ. 221) then
c
c SETUP OF THE COEFFICIENTS FOR THE NEW STAGES
c
c Here we give coefficients defining the interpolant for a 2-stage,
c 2nd order MIRK, stage order 1, with a 2nd order interpolant.
c
s_star = 2
c
c Define the sample point for the interpolant of the 2-stage, 2nd
c order MIRK discrete formula.
c
tau_star = 0.25d0
c
c Define the value of the order of the interpolant, `p'.
c
p = 2
c
elseif (method. EQ. 232) then
c
c SETUP OF THE COEFFICIENTS FOR THE NEW STAGES
c
c Here we give coefficients defining the interpolant for a 2-stage,
c 3rd order MIRK, stage order 2, with a 3rd order interpolant.
c
s_star = 3
c
c_star(1) = 1.0d0
c
v_star(1) = 1.0d0
c
c Column 1 of x_star
c
x_star(1) = 0.0d0
c
c Column 2 of x_star
c
x_star(2) = 0.0d0
c
c Column 3 of x_star
c
x_star(3) = 0.0d0
c
c Define the sample point for the interpolant of the 2-stage, 3rd
c order MIRK discrete formula.
c
tau_star = 0.25d0
c
c Define the value of the order of the interpolant, `p'.
c
p = 3
c
elseif (method. EQ. 2321) then
c
c SETUP OF THE COEFFICIENTS FOR THE NEW STAGES
c
c Here we give coefficients defining the interpolant for the
c reflection of the above 2-stage, 3rd order MIRK, stage order 2,
c with a 3rd order interpolant.
c
s_star = 3
c
c_star(1) = 0.0d0
c
v_star(1) = 0.0d0
c
c Column 1 of x_star
c
x_star(1) = 0.0d0
c
c Column 2 of x_star
c
x_star(2) = 0.0d0
c
c Column 3 of x_star
c
x_star(3) = 0.0d0
c
c Define the sample point for the interpolant of the 2-stage, 3rd
c order MIRK discrete formula.
c
tau_star = 0.25d0
c
c Define the value of the order of the interpolant, `p'.
c
p = 3
c
elseif (method. EQ. 343) then
c
c SETUP OF THE COEFFICIENTS FOR THE NEW STAGES
c
c Here we give coefficients defining the interpolant for a 3-stage,
c 4th order MIRK, stage order 3, with a 4th order interpolant.
c
c There is one extra stage. The corresponding row of the RK tableau
c is 3/4 | 27/32 | 3/64 -9/64 0 0
c
s_star = 4
c
c_star(1) = 3.0d0/4.0d0
c
v_star(1) = 27.0d0/32.0d0
c
c Column 1 of x_star
c
x_star(1) = 3.0d0/64.0d0
c
c Column 2 of x_star
c
x_star(2) = -9.0d0/64.0d0
c
c Column 3 of x_star
c
x_star(3) = 0.0d0
c
c Column 4 of x_star
c
x_star(4) = 0.0d0
c
c Define the sample point for the 3-stage, 4th order MIRK
c discrete formula.
c
tau_star = 0.226d0
c
c Define the value of the order of the interpolant, `p'.
c
p = 4
c
elseif (method .EQ. 453) then
c
c SETUP OF THE COEFFICIENTS FOR THE NEW STAGES
c
c Here we give coefficients defining the interpolant for a 4-stage,
c 5th order MIRK, stage order 3, with a 5th order interpolant.
c
s_star = 6
c
c_star(1) = 4.0d0/5.0d0
c_star(2) = 13.0d0/23.0d0
c
v_star(1) = 4.0d0/5.0d0
v_star(2) = 13.0d0/23.0d0
c
c Column 1 of x_star
c
x_star(1) = 14.0d0/1125.0d0
x_star(2) = 1.0d0/2.0d0
c
c Column 2 of x_star
c
x_star(3) = -74.0d0/875.0d0
x_star(4) = 4508233.0d0/1958887.0d0
c
c Column 3 of x_star
c
x_star(5) = -128.0d0/3375.0d0
x_star(6) = 48720832.0d0/2518569.0d0
c
c Column 4 of x_star
c
x_star(7) = 104.0d0/945.0d0
x_star(8) = -27646420.0d0/17629983.0d0
c
c Column 5 of x_star
c
x_star(9) = 0.0d0
x_star(10) = -11517095.0d0/559682.0d0
c
c Column 6 of x_star
c
x_star(11) = 0.0d0
x_star(12) = 0.0d0
c
c Define the sample point for the interpolant of the 4-stage, 5th
c order MIRK discrete formula.
c
tau_star = 0.30d0
c
c Define the value of the order of the interpolant, `p'.
c
p = 5
c
elseif (method .EQ. 4531) then
c
c SETUP OF THE COEFFICIENTS FOR THE NEW STAGES
c
c Here we give coefficients defining the interpolant for a 4-stage,
c 5th order MIRK, stage order 3, with a 5th order interpolant.
c
s_star = 6
c_star(1) = 7.0d0/13.0d0
c_star(2) = 13.0d0/23.0d0
c
v_star(1) = 4.0d0/5.0d0
v_star(2) = 23.0d0/33.0d0
c
c Column 1 of x_star
c
x_star(1) = 36445.0d0/1542294.0d0
x_star(2) = 1.0d0/2.0d0
c
c Column 2 of x_star
c
x_star(3) = -125759.0d0/1999270.0d0
x_star(4) = 1499810527.0d0/2327157756.0d0
c
c Column 3 of x_star
c
x_star(5) = -3224608.0d0/11567205.0d0
x_star(6) = -9339357584.0d0/2742721641.0d0
c
c Column 4 of x_star
c
x_star(7) = 915050.0d0/16194087.0d0
x_star(8) = -11845932250.0d0/4918765257.0d0
c
c Column 5 of x_star
c
x_star(9) = 0.0d0
x_star(10) = 514365992257.0d0/113365827828.0d0
c
c Column 6 of x_star
c
x_star(11) = 0.0d0
x_star(12) = 0.0d0
c
c Define the sample point for the interpolant of the 4-stage, 5th
c order MIRK discrete formula.
c
tau_star = 0.30d0
c
c Define the value of the order of the interpolant, `p'.
c
p = 5
c
elseif (method .EQ. 563) then
c
c SETUP OF THE COEFFICIENTS FOR THE NEW STAGES
c
c Here we give coefficients defining the interpolant for a 5-stage,
c 6th order MIRK, stage order 3, with a 6th order interpolant.
c
c There are four extra stages. The corresponding rows of the RK tableau:
c 7/16 | 7/16 | 1547/32768 -1225/32768 749/4096 -287/2048 -861/16384
c 3/8 | 3/8 | 83/1536 -13/384 283/1536 -167/1536 -49/512
c 9/16 | 9/16 | 1225/32768 -1547/32768 287/2048 -749/4096 861/16384
c 1/8 | 1/8 | 233/3456 -19/1152 0 0 0
c ... -5/72 7/72 -17/216
c
s_star = 9
c
c_star(1) = 7.0d0/16.0d0
c_star(2) = 3.0d0/8.0d0
c_star(3) = 9.0d0/16.0d0
c_star(4) = 1.0d0/8.0d0
c
v_star(1) = 7.0d0/16.0d0
v_star(2) = 3.0d0/8.0d0
v_star(3) = 9.0d0/16.0d0
v_star(4) = 1.0d0/8.0d0
c
c Column 1 of x_star
c
x_star(1) = 1547.0d0/32768.0d0
x_star(2) = 83.0d0/1536.0d0
x_star(3) = 1225.0d0/32768.0d0
x_star(4) = 233.0d0/3456.0d0
c
c Column 2 of x_star
c
x_star(5) = -1225.0d0/32768.0d0
x_star(6) = -13.0d0/384.0d0
x_star(7) = -1547.0d0/32768.0d0
x_star(8) = -19.0d0/1152.0d0
c
c Column 3 of x_star
c
x_star(9) = 749.0d0/4096.0d0
x_star(10) = 283.0d0/1536.0d0
x_star(11) = 287.0d0/2048.0d0
x_star(12) = 0.0d0
c
c Column 4 of x_star
c
x_star(13) = -287.0d0/2048.0d0
x_star(14) = -167.0d0/1536.0d0
x_star(15) = -749.0d0/4096.0d0
x_star(16) = 0.0d0
c
c Column 5 of x_star
c
x_star(17) = -861.0d0/16384.0d0
x_star(18) = -49.0d0/512.0d0
x_star(19) = 861.0d0/16384.0d0
x_star(20) = 0.0d0
c
c Column 6 of x_star
c
x_star(21) = 0.0d0
x_star(22) = 0.0d0
x_star(23) = 0.0d0
x_star(24) = -5.0d0/72.0d0
c
c Column 7 of x_star
c
x_star(25) = 0.0d0
x_star(26) = 0.0d0
x_star(27) = 0.0d0
x_star(28) = 7.0d0/72.0d0
c
c Column 8 of x_star
c
x_star(29) = 0.0d0
x_star(30) = 0.0d0
x_star(31) = 0.0d0
x_star(32) = -17.0d0/216.0d0
c
c Column 9 of x_star
c
x_star(33) = 0.0d0
x_star(34) = 0.0d0
x_star(35) = 0.0d0
x_star(36) = 0.0d0
c
c Define the sample point for the interpolant of the 5-stage, 6th
c order MIRK discrete formula.
c
tau_star = 0.7156d0
c
c Define the value of the order of the interpolant, `p'.
c
p = 6
c
elseif (method .EQ. 5631) then
c
c SETUP OF THE COEFFICIENTS FOR THE NEW STAGES
c
c Here we give coefficients defining the interpolant for a 5-stage,
c 6th order MIRK, stage order 3, with a 6th order interpolant.
c
c There are four extra stages. The corresponding rows of the RK tableau:
c 7/16 | 7/16 | 1547/32768 -1225/32768 749/4096 -287/2048 -861/16384
c 3/8 | 3/8 | 83/1536 -13/384 283/1536 -167/1536 -49/512
c 9/16 | 9/16 | 1225/32768 -1547/32768 287/2048 -749/4096 861/16384
c 1/8 | 1/8 | 233/3456 -19/1152 0 0 0
c ... -5/72 7/72 -17/216
c
s_star = 8
c
c_star(1) = 1.0d0/2.0d0
c_star(2) = 1.0d0/2.0d0-sqrt(7.0d0)/14.0d0
c_star(3) = 87.0d0/100.0d0
c
v_star(1) = 1.0d0/2.0d0
v_star(2) = 1.0d0/2.0d0-sqrt(7.0d0)/14.0d0
v_star(3) = 87.0d0/100.0d0
c
c (The elements of the matrix containing the stage coefficients are
c stored columnwise in the vector `x_star'.
c Column 1 of x_star
c
x_star(1) = 1.D0/64.D0
x_star(2) = 3.D0/112.D0+9.D0/1960.D0*dsqrt(7.D0)
x_star(3) = 2707592511.D0/1000000000000.D0-1006699707.D0/
#1000000000000.D0*dsqrt(7.D0)
c
c Column 2 of x_star
c
x_star(4) = -1.D0/64.D0
x_star(5) = -3.D0/112.D0+9.D0/1960.D0*dsqrt(7.D0)
x_star(6) = -51527976591.D0/1000000000000.D0-1006699707.D0/
#1000000000000.D0*dsqrt(7.D0)
c
c Column 3 of x_star
c
x_star(7) = 7.D0/192.D0*dsqrt(21.D0)
x_star(8) = 11.D0/840.D0*dsqrt(7.D0)+3.D0/112.D0*dsqrt(7.D0)
#*dsqrt(3.D0)
x_star(9) = -610366393.D0/75000000000.D0+7046897949.D0/
#1000000000000.D0*dsqrt(7.D0)+14508670449.D0/1000000000000.D0
#*dsqrt(7.D0)*dsqrt(3.D0)
c
c Column 4 of x_star
c
x_star(10) = -7.D0/192.D0*dsqrt(21.D0)
x_star(11) = 11.D0/840.D0*dsqrt(7.D0)-3.D0/112.D0*dsqrt(7.D0)
#*dsqrt(3.D0)
x_star(12) = -610366393.D0/75000000000.D0+7046897949.D0/
#1000000000000.D0*dsqrt(7.D0)-14508670449.D0/1000000000000.D0
#*dsqrt(7.D0)*dsqrt(3.D0)
c
c Column 5 of x_star
c
x_star(13) = 0
x_star(14) = 88.D0/5145.D0*dsqrt(7.D0)
x_star(15) = -12456457.D0/1171875000.D0+1006699707.D0/
#109375000000.D0*dsqrt(7.D0)
c
c Column 6 of x_star
c
x_star(16) = 0
x_star(17) = -18.D0/343.D0*dsqrt(7.D0)
x_star(18) = 3020099121.D0/437500000000.D0*dsqrt(7.D0)+
#47328957.D0/625000000.D0
c
c Column 7 of x_star
c
x_star(19) = 0
x_star(20) = 0
x_star(21) = -7046897949.D0/250000000000.D0*dsqrt(7.D0)
c
c Column 8 of x_star
c
x_star(22) = 0
x_star(23) = 0
x_star(24) = 0
c
c Define the sample point for the interpolant of the 5-stage, 6th
c order MIRK discrete formula.
c
tau_star = 0.4d0
c
c Define the value of the order of the interpolant, `p'.
c
p = 6
c
endif
c
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine interp_weights(s_star,tau,w,wp)
c
c***overview***
c
c This routine evaluates the weight polynomials (and their first
c derivatives) corresponding to the `s_star' stages used to form the
c interpolant for the discrete formula defined by the `RK_tableau'
c routine. The evaluation takes place at `tau'. The results are returned
c in the `s_star' locations of the vectors `w' and `wp'. The
c interp_weights routine will be dependent on the parameter, `method'.
c `method' defines the method to be used and will be accessed via the
c /rk_method/ common block. The possible values for `method' and the
c corresponding MIRK method associated with each of these values are
c listed below.
c
c Value of `method' | MIRK formula
c -------------------------------------
c 121 | MIRK121 scheme
c 221 | MIRK221 scheme
c 232 | MIRK232 scheme (abscissa = 0,2/3)
c 2321 | MIRK232 scheme (abscissa = 1,1/3)
c 343 | MIRK343 scheme
c 453 | MIRK453 scheme
c 563 | MIRK563 scheme
c 5631 | An improved MIRK563 scheme
c
c The presence of the interp_weights routine demonstrates the use
c of an "information hiding" principle. The more common approach would be
c to simply provide the coefficients of the polynomial weights in
c an array made available to the rest of the software package through
c the /interp_coeffs/ common block. In the approach employed here we
c "hide" the polynomial expression for each weight inside the
c `interp_weights' routine whose function is to provide an evaluation
c of the weights at a given point. An advantage of this approach is
c that changes can be made in the method by which the weights are
c determined with no subsequent change needed in the rest of the
c software.
c----------------------------------------------------------------------------
c***declaration of parameters***
c imports:
integer s_star
double precision tau
c
c `s_star' is the number of weight polynomials.
c `tau' is the point at which the weight polynomials
c are to be evaluated.
c exports:
double precision w(s_star), wp(s_star)
c
c `w' contains the values of the weight polynomials
c of the interpolant, evaluated at `tau'.
c
c `wp' contains the values of the first derivative
c of the weight polynomials of the interpolant,
c evaluated at `tau'.
c
c***declaration of local variables***
double precision t2, tm1, t4m3, t2m1
c
c Factors appearing in the expressions for the weight polynomials
c and derivatives, defined below.
c
c***declaration of variables for the /rk_method/ common block***
c exports:
integer method
common /rk_method/ method
c
c `method' defines the method to be used
c------------------------------------------------------------------------------
c called by: `defect_estimate', `interp_eval'
c------------------------------------------------------------------------------
c The polynomials in `tau' that define the interpolant are very
c dependent on the underlying discrete formula and the interpolant
c derived. The following part of this routine will be different for
c each formula implemented.
c
if (method. EQ. 121) then
c
c WEIGHT POLYNOMIAL IMPLEMENTATION for a 1-stage, 2nd order MIRK formula
c
w(1) = 0.0d0
c
w(2) = tau*(1.0d0-tau/2.0d0)
c
w(3) = tau**2/2.0d0
c
c Derivative polynomials.
c
wp(1) = 0.0d0
c
wp(2) = 1.0d0 - tau
c
wp(3) = tau
c
elseif (method. EQ. 221) then
c
c WEIGHT POLYNOMIAL IMPLEMENTATION for a 2-stage, 2nd order MIRK formula
c
w(1) = tau*(1.0d0-tau/2.0d0)
c
w(2) = tau**2/2.0d0
c
c Derivative polynomials.
c
wp(1) = 1.0d0 - tau
c
wp(2) = tau
c
elseif (method .EQ. 232) then
c
c WEIGHT POLYNOMIAL IMPLEMENTATION for a 2-stage, 3rd order MIRK formula
c
w(1) = tau*(2.0d0*tau**2 - 5.0d0*tau + 4.0d0)/4.0d0
c
w(2) = -3.0d0*tau**2*(2.0d0*tau - 3.0d0)/4.0d0
c
w(3) = tau**2*(tau - 1.0d0)
c
c Derivative polynomials.
c
wp(1) = (tau - 2.0d0/3.0d0)*(tau - 1.0d0)/(2.0d0/3.0d0)
c
wp(2) = tau*(tau - 1.0d0)/(-2.0d0/9.0d0)
c
wp(3) = tau*(tau - 2.0d0/3.0d0)/(1.0d0/3.0d0)
c
elseif (method .EQ. 2321) then
c
c WEIGHT POLYNOMIAL IMPLEMENTATION for the reflection of the above
c 2-stage, 3rd order MIRK formula
c
w(1) = 1.0d0/2.0d0*tau**3 - 1.0d0/4.0d0*tau**2
c
w(2) = -3.0d0/2.0d0*tau**3 + 9.0d0/4.0d0*tau**2
c
w(3) = tau**3 - 2.0d0*tau**2 + tau
c
c Derivative polynomials.
c
wp(1) = 3.0d0*tau*(tau-1.0d0/3.0d0)/2.0d0
c
wp(2) = -9.0d0*tau*(tau-1.0d0)/2.0d0
c
wp(3) = 3.0d0*(tau-1.0d0)*(tau-1.0d0/3.0d0)
c
elseif (method .EQ. 343) then
c
c WEIGHT POLYNOMIAL IMPLEMENTATION for a 3-stage, 4th order MIRK formula
c
t2 = tau*tau
tm1 = tau - 1.0d0
t4m3 = tau*4.0d0 - 3.0d0
t2m1 = tau*2.0d0 - 1.0d0
c
w(1) = -tau*(2.0d0*tau-3.0d0)*(2.0d0*t2-3.0d0*tau+2.0d0)/6.0d0
c
w(2) = t2*(12.0d0*t2-20.0d0*tau + 9.0d0)/6.0d0
c
w(3) = 2.0d0*t2*(6.0d0*t2-14.0d0*tau+9.0d0)/3.0d0
c
w(4) = -16.0d0*t2*tm1*tm1/3.0d0
c
c Derivative polynomials.
c
wp(1) = -tm1*t4m3*t2m1/3.0d0
c
wp(2) = tau*t2m1*t4m3
c
wp(3) = 4.0d0*tau*t4m3*tm1
c
wp(4) = -32.0d0*tau*t2m1*tm1/3.0d0
c
elseif (method .EQ. 453) then
c
w(1) = tau*(22464.0d0-83910.0d0*tau+143041.0d0*tau**2
+ -113808.0d0*tau**3+33256.0d0*tau**4)/22464.0d0
c
w(2) = tau**2*(-2418.0d0+12303.0d0*tau-19512.0d0*tau**2
++10904.0d0*tau**3)/3360.0d0
c
w(3) = -8.0d0/81.0d0*tau**2*(-78.0d0+209.0d0*tau-204.0d0
+*tau**2+8.0d0*tau**3)
c
w(4) = -25.0d0/1134.0d0*tau**2*(-390.0d0+1045.0d0*tau
+-1020.0d0*tau**2+328.0d0*tau**3)
c
w(5) = -25.0d0/5184.0d0*tau**2*(390.0d0+255.0d0*tau-
+1680.0d0*tau**2+2072.0d0*tau**3)
c
w(6) = 279841.0d0/168480.0d0*tau**2*(-6.0d0+21.0d0*tau
+-24.0d0*tau**2+8.0d0*tau**3)
c
c Derivative polynomials.
c
wp(1) = 1.0d0-13985.0d0/1872.0d0*tau+143041.0d0/7488.0d0
+*tau**2-2371.0d0/117.0d0*tau**3+20785.0d0/2808.0d0*tau**4
c
wp(2) = -403.0d0/280.0d0*tau+12303.0d0/1120.0d0*tau**2
+-813.0d0/35.0d0*tau**3+1363.0d0/84.0d0*tau**4
c
wp(3) = 416.0d0/27.0d0*tau-1672.0d0/27.0d0*tau**2+2176.0d0
+/27.0d0*tau**3-320.0d0/81.0d0*tau**4
c
wp(4) = 3250.0d0/189.0d0*tau-26125.0d0/378.0d0*tau**2+
+17000.0d0/189.0d0*tau**3-20500.0d0/567.0d0*tau**4
c
wp(5) = -1625.0d0/432.0d0*tau-2125.0d0/576.0d0*tau**2
++875.0d0/27.0d0*tau**3-32375.0d0/648.0d0*tau**4
c
wp(6) = -279841.0d0/14040.0d0*tau+1958887.0d0/18720.0d0*
+tau**2-279841.0d0/1755.0d0*tau**3+279841.0d0/4212.0d0*tau**4
c
elseif (method .EQ. 4531) then
c
w(1) = tau*(3857490.0d0-14075970.0d0*tau+24780520.0d0
+*tau**2-20977215.0d0*tau**3+7479584.0d0*tau**4)/3857490.0d0
c
w(2) = tau**2*(-959595.0d0+3459675.0d0*tau-4676415.0d0
+*tau**2+1893368.0d0*tau**3)/989100.0d0
c
w(3) = -32.0d0/63585.0d0*tau**2*(-13650.0d0+42100.0d0*tau
+-47175.0d0*tau**2+10376.0d0*tau**3)
c
w(4) = -250.0d0/89019.0d0*tau**2*(-2730.0d0+8420.0d0*tau
+-9435.0d0*tau**2+4336.0d0*tau**3)
c
w(5) = 28561.0d0/791280.0d0*tau**2*(390.0d0+255.0d0*tau
+-1680.0d0*tau**2+2072.0d0*tau**3)
c
w(6) = -279841.0d0/2449200.0d0*tau**2*(210.0d0-225.0d0*
+tau-180.0d0*tau**2+536.0d0*tau**3)
c
wp(1) = 1.0d0-938398.0d0/128583.0d0*tau+2478052.0d0/128583.0d0
+*tau**2-399566.0d0/18369.0d0*tau**3+534256.0d0/55107.0d0*tau**4
c
wp(2) = -9139.0d0/4710.0d0*tau+46129.0d0/4396.0d0*tau**2-
+311761.0d0/16485.0d0*tau**3+473342.0d0/49455.0d0*tau**4
c
wp(3) = 58240.0d0/4239.0d0*tau-269440.0d0/4239.0d0*tau**2+
+402560.0d0/4239.0d0*tau**3-332032.0d0/12717.0d0*tau**4
c
wp(4) = 65000.0d0/4239.0d0*tau-2105000.0d0/29673.0d0*tau**2
++3145000.0d0/29673.0d0*tau**3-5420000.0d0/89019.0d0*tau**4
c
wp(5) = 371293.0d0/13188.0d0*tau+485537.0d0/17584.0d0*tau**2
+-114244.0d0/471.0d0*tau**3+1056757.0d0/2826.0d0*tau**4
c
wp(6) = -1958887.0d0/40820.0d0*tau+2518569.0d0/32656.0d0*
+tau**2+839523.0d0/10205.0d0*tau**3-18749347.0d0/61230.0d0*tau**4
c
elseif (method .EQ. 563) then
c
c WEIGHT POLYNOMIAL IMPLEMENTATION for a 5-stage, 6th order MIRK formula
c
w(1) =
+tau-28607.D0/7434.D0*tau**2-166210.D0/33453.D0*tau**3+334780.D0/1
+1151.D0*tau**4-1911296.D0/55755.D0*tau**5+406528.D0/33453.D0*tau*
+*6
c
w(2) =
+777.D0/590.D0*tau**2-2534158.D0/234171.D0*tau**3+0.208858D7/78057
+.D0*tau**4-10479104.D0/390285.D0*tau**5+0.11328512D8/0.1170855D7*
+tau**6
c
w(3) =
+-1008.D0/59.D0*tau**2+222176.D0/1593.D0*tau**3-180032.D0/531.D0
+*tau**4+876544.D0/2655.D0*tau**5-180224.D0/1593.D0*tau**6
c
w(4) =
+-1008.D0/59.D0*tau**2+222176.D0/1593.D0*tau**3-180032.D0
+/531.D0*tau**4+876544.D0/2655.D0*tau**5-180224.D0/1593.D0*tau**6
c
w(5) =
+-378.D0/59.D0*tau**2+27772.D0/531.D0*tau**3-22504.D0/177.D0*tau*
+*4+109568.D0/885.D0*tau**5-22528.D0/531.D0*tau**6
c
w(6) =
+-95232.D0/413.D0*tau**2+0.62384128D8/33453.D0*tau**3-49429504.D0
+/11151.D0*tau**4+0.46759936D8/11151.D0*tau**5-46661632.D0/33453.
+D0*tau**6
c
w(7) =
+896.D0/5.D0*tau**2-4352.D0/3.D0*tau**3+3456*tau**4-16384.D0/5.D0
+*tau**5+16384.D0/15.D0*tau**6
c
w(8) =
+50176.D0/531.D0*tau**2-179554304.D0/234171.D0*tau**3+0.143363072D9
+/78057.D0*tau**4-136675328.D0/78057.D0*tau**5+0.137363456D9
+/234171.D0*tau**6
c
w(9) =
+16384.D0/441.D0*tau**3-16384.D0/147.D0*tau**4+16384.D0/147.D0
+*tau**5-16384.D0/441.D0*tau**6
c
c Derivative polynomials.
c
wp(1) =
+1-28607.D0/3717.D0*tau-166210.D0/11151.D0*tau**2+0.133912D7/
+11151.D0*tau**3-1911296.D0/11151.D0*tau**4+813056.D0/11151.
+D0*tau**5
c
wp(2) =
+777.D0/295.D0*tau-2534158.D0/78057.D0*tau**2+0.835432D7
+/78057.D0*tau**3-10479104.D0/78057.D0*tau**4+0.22657024D8
+/390285.D0*tau**5
c
wp(3) =
+-2016.D0/59.D0*tau+222176.D0/531.D0*tau**2-720128.D0/531.D0
+*tau**3+876544.D0/531.D0*tau**4-360448.D0/531.D0*tau**5
c
wp(4) =
+-2016.D0/59.D0*tau+222176.D0/531.D0*tau**2-720128.D0/531.D0
+*tau**3+876544.D0/531.D0*tau**4-360448.D0/531.D0*tau**5
c
wp(5) =
+-756.D0/59.D0*tau+27772.D0/177.D0*tau**2-90016.D0/177.D0*
+tau**3+109568.D0/177.D0*tau**4-45056.D0/177.D0*tau**5
c
wp(6) =
+-190464.D0/413.D0*tau+0.62384128D8/11151.D0*tau**2-1977
+18016.D0/11151.D0*tau**3+0.23379968D9/11151.D0*tau**4-
+93323264.D0/11151.D0*tau**5
c
wp(7) =
+1792.D0/5.D0*tau-4352*tau**2+13824*tau**3-16384*tau**4
++32768.D0/5.D0*tau**5
c
wp(8) =
+100352.D0/531.D0*tau-179554304.D0/78057.D0*tau**2+0.57345
+2288D9/78057.D0*tau**3-683376640.D0/78057.D0*tau**4+0.274
+726912D9/78057.D0*tau**5
c
wp(9) =
+16384.D0/147.D0*tau**2-65536.D0/147.D0*tau**3+81920.D0/147
+.D0*tau**4-32768.D0/147.D0*tau**5
c
elseif (method .EQ. 5631) then
c
c WEIGHT POLYNOMIAL IMPLEMENTATION for an improved 5-stage, 6th order
c MIRK formula
c
w(1) =
#-(12233+1450*dsqrt(7.D0))*(800086000*tau**5+63579600*dsqrt
#(7.D0)*tau**4-2936650584.D0*tau**4+4235152620.D0*tau**3
#-201404565*dsqrt(7.D0)*tau**3+232506630*dsqrt(7.D0)*tau**2
#-3033109390.D0*tau**2+1116511695*tau-116253315*dsqrt(7.D0)
#*tau+22707000*dsqrt(7.D0)-191568780)*tau/2112984835740.D0
c
w(2) =
#-(-10799+650*dsqrt(7.D0))*(24962000*tau**4
#+473200*dsqrt(7.
#D0)*tau**3-67024328*tau**3-751855*dsqrt(7.D0)*tau**2
#+66629600*tau**2-29507
#250*tau+236210*dsqrt(7.D0)*tau+5080365+50895*dsqrt(7.D0))
#*tau**2/29551834260.D0
c
w(3) =
#7.D0/1274940.D0*(259+50*dsqrt(7.D0))
#*(14000*tau**4-48216*tau
#**3+1200*dsqrt(7.D0)*tau**3-3555*dsqrt(7.D0)
#*tau**2+62790*tau**2+3610*ds
#qrt(7.D0)*tau-37450*tau+9135-1305*dsqrt(7.D0))*tau**2
c
w(4) =
#7.D0/1274940.D0*(259+50*dsqrt(7.D0))*(14000*tau**4
#-48216*tau
#**3+1200*dsqrt(7.D0)*tau**3-3555*dsqrt(7.D0)
#*tau**2+62790*tau**2+3610*ds
#qrt(7.D0)*tau-37450*tau+9135-1305*dsqrt(7.D0))*tau**2
c
w(5) =
#16.D0/2231145.D0*(259+50*dsqrt(7.D0))
#*(14000*tau**4-48216*
#tau**3+1200*dsqrt(7.D0)*tau**3-3555*dsqrt(7.D0)*tau**2
#+62790*tau**2+3610*d
#sqrt(7.D0)*tau-37450*tau+9135
#-1305*dsqrt(7.D0))*tau**2
c
w(6) =
#4.D0/1227278493.D0*(740*dsqrt(7.D0)-6083)*(1561000*tau**2-
#2461284*tau-109520*dsqrt(7.D0)*tau+979272+86913*
#dsqrt(7.D0))*(tau-1)**2*tau**2
c
w(7) =
#-49.D0/63747.D0*dsqrt(7.D0)*(20000*tau**2
#-20000*tau+3393)*(tau-1)**2*tau**2
c
w(8) =
#-1250000000.D0/889206903.D0*(28*tau**2-28*tau+9)
#*(tau-1)**2*tau**2
c
c Derivative polynomials.
c
wp(1) =
#(1450*dsqrt(7.D0)+12233)*(14*tau-7
#+dsqrt(7.D0))*(tau-1)*(-4
#00043*tau+75481+2083*dsqrt(7.D0))*(100*tau-87)
#*(2*tau-1)/493029795006.D0
c
wp(2) =
# -(650*dsqrt(7.D0)-10799)*(14*tau-7+dsqrt(7.D0))*(37443*tau-
#13762-2083*dsqrt(7.D0))*(100*tau-87)
#*(2*tau-1)*tau/20686283982.D0
c
wp(3) =
# 7.D0/42498.D0*(259+50*dsqrt(7.D0))*(14*tau-7+dsqrt(7.D0))
#*(tau-1)*(100*tau-87)*(2*tau-1)*tau
c
wp(4) =
# 7.D0/42498.D0*(259+50*dsqrt(7.D0))*(14*tau-7+dsqrt(7.D0))
#*(tau-1)*(100*tau-87)*(2*tau-1)*tau
c
wp(5) =
# 32.D0/148743.D0*(259+50*dsqrt(7.D0))*(14*tau-7+dsqrt(7.D0
#))*(tau-1)*(100*tau-87)*(2*tau-1)*tau
c
wp(6) =
# 4.D0/1227278493.D0*(740*dsqrt(7.D0)-6083)*(14*tau-7+dsqrt
#(7.D0))*(tau-1)*(100*tau-87)*(6690*tau-4085
#-869*dsqrt(7.D0))*tau
c
wp(7) =
# -98.D0/21249.D0*dsqrt(7.D0)*(tau-1)*(100*tau-13)
#*(100*tau-87)*(2*tau-1)*tau
c
wp(8) =
# -1250000000.D0/2074816107.D0*(14*tau-7+dsqrt(7.D0))
#*(tau-1)*(14*tau-7-dsqrt(7.D0))*(2*tau-1)*tau
c
endif
c
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine interval(Nsub, mesh, t, i)
c
c***overview***
c
c This routine finds `i' such that `mesh(i-1) =< t < mesh(i)' for
c `t < mesh(Nsub)'. The value of `Nsub' is returned when `t = mesh(Nsub)'.
c `Nsub' is the number of subintervals and `mesh' is the list of points
c defining the subintervals.
c----------------------------------------------------------------------------
c***declaration of parameters***
c imports:
integer Nsub
double precision mesh(0:Nsub), t
c
c `Nsub' Number of subintervals.
c `mesh' Array of points defining the subintervals.
c `t' Point to be located within `mesh'.
c
c exports:
integer i
c
c `i' Mesh index such that `t' satisfies `mesh(i-1) =< t < mesh(i)'
c for `t < mesh(Nsub)' or `i = Nsub for t = mesh(Nsub)'.
c
c***declaration of local variables***
double precision meshhalfnew
integer inew, iold, iright, ileft
c
c `meshhalfnew' New mesh value resulting from halving the current interval.
c `inew' Current subinterval containing `meshhalfnew'.
c `iold' Previous subinterval containing `meshhalfnew'.
c `ileft' Left hand mesh point bounding `meshhalfnew'.
c `iright' Right hand mesh point bounding `meshhalfnew'.
c---------------------------------------------------------------------------
c Called by: `interp_eval' or `sol_eval'.
c---------------------------------------------------------------------------
c Handle special cases for the location of `t' first
c
if ((t .LT. mesh(0)) .OR. (t .GT. mesh(Nsub))) then
write(6,*)'Error-t is outside of problem interval'
endif
c
if (t .EQ. mesh(Nsub)) then
i = Nsub
return
end if
c
if (t .EQ. mesh(0)) then
i = 1
return
end if
c
c At this point we know that `mesh(0) < t < mesh(Nsub)'.
c Binary search to determine subinterval containing `t'.
c
inew = (Nsub/2)
meshhalfnew = mesh(inew)
c
c Determine if `t' is equal to `meshhalfnew'.
c
if (t .EQ. meshhalfnew) then
i = inew + 1
return
endif
c
c `t' is not equal to the meshpoint.
c
if (t .LT. meshhalfnew) then
c
c `t' is on the left hand side of the interval.
c
dowhile (t .LT. meshhalfnew)
iold = inew
inew = iold/2
meshhalfnew = mesh(inew)
enddo
c
iright = iold
ileft = inew
c
else
c
c `t' is on the right hand side of the interval.
c
dowhile (t .GT. meshhalfnew)
iold = inew
inew = nint(real((iold + (Nsub - iold)/2)) + 0.5)
meshhalfnew = mesh(inew)
enddo
c
ileft = iold
iright = inew
c
endif
c
i = iright
c
c Determine if `t' is equal to `meshhalfnew'.
c
if (t .EQ. meshhalfnew) then
i = inew + 1
return
endif
c
c `t' is not equal to the meshpoint.
c
c `t' is bounded by `ileft' and `iright'.
c
dowhile (abs(iright - ileft) .NE. 1)
inew = (ileft + iright)/2
if(mesh(inew) .LT. t) then
ileft = inew
elseif (mesh(inew) .GT. t) then
iright = inew
else
i = inew + 1
return
endif
enddo
c
i = iright
c
c `t' is in the `i'-th subinterval.
c
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine jacblk(neqns,hi,ti,yi,yip1,k,LLi,RRi,
+ yr,dfdy,ur,dkdyi,dkdyip1,dfsub)
c
c***overview***
c
c ****Evaluation of the Jacobian for one subinterval****
c
c This routine computes the left and right blocks (of the Jacobian of
c the residual function) associated with the i-th subinterval. The left
c and right block values are returned in the parameters `LLi' and `RRi'.
c Each block is of size `neqns**2', where `neqns' is the number of
c differential equations. We have,
c
c s s
c -- --
c LLi = -I - hi * > br * dkr/dyi, RRi = I - hi * > br * dkr/dyi+1
c -- --
c r=1 r=1
c
c Each partial derivative, dkr/dyi, is computed as follows,
c
c r-1
c --
c dkr/dyi = df/dy(tr,yr)*((1-vr)*I+hi*> xrj * dkj/dyi)
c --
c j=1,
c
c with a similar expression for each dkr/dyi+1.
c
c The solution value, yr, at which df/dy is evaluated depends on
c the input parameters `yi' and `yip1', the solution approximations
c at each end of the ith subinterval, and the previously computed
c stages provided by the input parameter `k'. The size of the ith
c subinterval is `hi' and the endpoint is `ti'. The computation of
c `yr' also depends on the coefficients of the RK method provided
c through the /RK_s/,/RK_coeff/ common blocks.
c
c (All matrices employed in this routine are actually stored
c columnwise within vectors.)
c--------------------------------------------------------------------------------
c***declaration of constants***
integer Mxs
parameter(Mxs=10)
c
c `Mxs' is the maximum value for the total number of stages
c of the interpolant.
c
c***declaration of parameters***
c imports:
integer neqns
double precision hi, ti, yi(neqns),yip1(neqns)
double precision k(Mxs*neqns)
c
c `neqns' is the number of differential equations.
c
c `hi' is the size of the current subinterval.
c
c `ti' is the value of the lefthand endpoint of the
c ith subinterval.
c
c `yi' is the current solution approximation at `ti'.
c
c `yip1' is the current solution approximation at `ti+h', i.e.
c at the right hand endpoint of the subinterval.
c
c `k' contains the `s' stage evaluations, each of size
c `neqns'associated with the `ith' subinterval.
c exports:
double precision LLi(neqns**2),RRi(neqns**2)
c
c `LLi' is the partial derivative of the ith component of
c the residual function with respect to yi
c
c `RRi' is the partial derivative of the ith component of the
c residual function with respect to yi+1
c
c***declaration of work arrays***
double precision yr(neqns)
double precision dfdy(neqns**2)
double precision ur(neqns**2)
double precision dkdyi(neqns**2*Mxs)
double precision dkdyip1(neqns**2*Mxs)
c `yr' the argument for the evaluation of the `rth' stage.
c
c `dfdy' during the `rth' time through the main loop below, it
c contains the derivative of the right hand side of the
c system of ODE's, evaluated at the arguments of the
c 'rth' stage of the RK scheme.
c
c `ur' used to store the intermediate values arising in the
c computation of each partial derivative.
c
c `dkdyi' The `rth' set of `neqns**2' locations of this vector
c will contain the value of the rth partial derivative
c with respect to `yi'.
c
c `dkdyip1' The `rth' set of `neqns**2' locations of this vector
c will contain the value of the rth partial derivative
c with respect to `yip1'.
c
c***user - supplied subroutines***
external dfsub
c
c dfsub - This routine defines the Jacobian of the system of
c differential equations with respect to y, i.e.
c df/dy.
c
c***declaration of local variables***
integer l,i,m,r, j,i1,i2
double precision sum
integer neqnssq
integer alpha,beta,sigma
double precision tr
logical dkdyi_zero(Mxs)
logical dkdyip1_zero(Mxs)
logical matrix_mult_dkdyi
logical matrix_mult_dkdyip1
c
c `l,i,m,r,j' are loop index variables
c
c `sum' is an accumulator
c
c `i1',`i2' are indexes used in calculating ur and yr
c
c `neqnssq' is the value of neqns squared
c
c `alpha',`beta',`sigma' are indexes used in the
c multiplication of the two matrices `ur' and `dfdy'
c
c `tr' the abscissa for the evaluation of the `rth' stage.
c
c `dkdyi_zero' is a logical array. If the computation of the rth
c stage has been completed then a value of `.TRUE.' in
c `dkdyi_zero(r)' indicates that every component of
c the rth partial derivative with respect to `yi'
c is zero. If the computation of the rth stage has
c not been completed then a value of `.TRUE.'
c indicates that the `ur' components of the rth
c partial derivative are currently zero. A value of
c `.FALSE.' in `dkdyi_zero(r)' implies the rth
c partial derivative is non-zero in at least one of
c its components.
c
c `dkdyip1_zero' is a logical array. If the computation of the rth
c stage has been completed then a value of `.TRUE.'
c in `dkdyip1_zero(r)' indicates that every component
c of the rth partial derivative with respect to
c `yip1' is zero. If the computation of the rth
c stage has not been completed then a value of
c `.TRUE.' indicates that the `ur' components of
c the rth partial derivative are currently zero.
c A value of `.FALSE.' in `dkdyip1_zero(r)'
c implies the rth partial derivative is non-zero
c in at least one of its components. `dkdyi_zero'
c and `dkdyip1_zero' are both used to avoid
c multiplications by zero.
c
c `matrix_mult_dkdyi' is a logical variable. For the compuation
c of a given stage, a value of `.FALSE.'
c indicates that the partial derivative
c with respect to `yi'
c has the form `(1-v(r))*dfdy' and that
c multiplying the matrix `ur' by the matrix
c `dfdy' is unnecessary. A value of `.TRUE.'
c indicates that the partial derivative
c does not have the form `(1-v(r))*dfdy' and
c that multiplying the matrix `ur' by the
c matrix `dfdy' may be necessary.
c
c `matrix_mult_dkdyip1' is a logical variable. For the compuation
c of a given stage, a value of `.FALSE.'
c indicates that the partial derivative
c with respect to `yip1'
c has the form `v(r)*dfdy' and that
c multiplying the matrix `ur' by the matrix
c `dfdy' is unnecessary. A value of `.TRUE.'
c indicates that the partial derivative
c does not have the form `v(r)*dfdy' and
c that multiplying the matrix `ur' by the
c matrix `dfdy' may be necessary.
c
c***declaration of variables in common block /RK_s/***
integer s
common /RK_s/ s
c
c `s' is the number of discrete stages of the Runge-Kutta formula.
c
c***declaration of variables in common block /RK_coeff/***
double precision v(Mxs), b(Mxs)
double precision c(Mxs), x(Mxs**2)
common /RK_coeff/ v, b, c, x
c
c `v', `b', `c',and `x' are coefficients that define
c the discrete Runge-Kutta formula.
c
c The first component of `c' contains the abscissa for the first
c new stage. The second component of `c' contains the abscissa for
c the second new stage, and so on.
c
c The first component of `v' contains the blending coefficient
c for the first new stage. The second component of `v' contains
c the blending coefficient for the second new stage, and so on.
c
c The matrix, X, is related to the data structure, `x'
c as follows (the matrix is stored columnwise in the vector):
c
c X(i,j) = x( (j-1)*s + i )
c
c `x' contains the coupling coefficients for each of the new
c stages
c
c------------------------------------------------------------------------------
c called by: `newmat'
c calls to: `dfsub'
c------------------------------------------------------------------------------
c***COMPUTATION OF THE PARTIAL DERIVATIVES***
c
c We compute the partial derivatives in pairs: dkr/dyi and dkr/dyi+1
c during the `rth' pass of the following loop. The results are stored
c in the appropriate locations of the vectors `dkdyi' and `dkdyip1',
c respectively.
c
c Initialize the components `dkdyi_zero' and `dkdyi_zero' to `.TRUE.'
c indicating that all of the partial derivatives are currently zero.
c
do 1 j = 1,s
dkdyi_zero(j) = .TRUE.
dkdyip1_zero(j) = .TRUE.
1 continue
c
neqnssq = neqns**2
c
do 40 r = 1, s
c
c Initialize `matrix_mult_dkdyi' to `.TRUE.' indicating that the
c partial derivative is not of the form `(1-v(r))*dfdy' and that
c multiplying the matrix `ur' by `dfdy' may still be necessary.
c
matrix_mult_dkdyi = .TRUE.
c
c Construct the arguments for the `rth' stage, needed here for the
c evaluation of `df/dy'
c
c Compute the abscissa
c
tr = ti + c(r)*hi
c
c Begin the construction of `yr' by setting it to zero.
c
do 2 j=1, neqns
yr(j) = 0.0d0
2 continue
c
c Add on weighted contributions from first `r-1' stages.
c The `rth' stage of the current subinterval is stored in
c `k((r-1)*neqns+1)'.
c
do 5 j = 1, r-1
c
i1 = (j-1)*s + r
i2 = (j-1)*neqns
c
if ( x(i1) .NE. 0.0d0 ) then
c
c If `x(i1) = 0' then there is no need to perform the following
c multiplication.
do 3 l = 1, neqns
yr(l) = yr(l) + x(i1)*k(i2+l)
3 continue
c
endif
c
5 continue
c
c Multiply yr by hi.
c
do 7 l =1, neqns
yr(l) = hi*yr(l)
7 continue
c
c Add on weighted contributions from `yi' and `yip1'.
c
do 10 l=1,neqns
yr(l) = (1.0d0 - v(r))*yi(l) + v(r)*yip1(l) + yr(l)
10 continue
c
c The arguments for the `rth' stage, `tr' and `yr' have been
c constructed. Evaluate df/dy, the derivative of the system of ODE's
c by first setting it equal to zero then calling the user defined
c subroutine `dfsub' at these arguments to fill in the non-zero
c entries.
c
do 12 i = 1,neqnssq
dfdy(i) = 0.0d0
12 continue
c
call dfsub(neqns,tr,yr,dfdy)
c
c The values of the `rth' partial derivatives with respect to `yi'
c and `yip1' are obtained by multiplying `dfdy' by a weighted linear
c combination of the partial derivatives of the previous `r-1'
c partial derivatives. We will store these linear combinations, as we
c construct them, in the matrix, `ur'.
c
c Construct the `ur' matrix for dkr/dyi
c
do 14 j = 1, neqnssq
ur(j) = 0.0d0
14 continue
c
c Add on weighted multiples of the previous partial derivatives.
c The `jth' previous partial derivative is stored in
c `dkdyi((j-1)*neqns**2+1)'. This can be done through a call to
c `daxpy' since we are storing all matricies in vector form.
c
do 15 j = 1, r-1
c
i1 = (j-1)*s + r
i2 = (j-1)*neqnssq
c
if ( x(i1) .NE. 0 .AND. ( .NOT. dkdyi_zero(j) ) ) then
c
c If `x(i1) = 0' or the jth partial derivative equals zero then
c there is no need to perform the following multiplication.
c
dkdyi_zero(r) = .FALSE.
c Indicates the matrix `ur' is non zero in at least one of the
c components and we will require a matrix multiply later.
c It also indicates that the rth partial derivative of
c `dkdyi' will be non zero in one of the components.
c
do 16 l= 1, neqnssq
ur(l) = ur(l) + x(i1)*dkdyi(i2+l)
16 continue
c
endif
c
15 continue
c
c Multiply ur by hi
c
if ( .NOT. dkdyi_zero(r) ) then
c
c If the `ur' matrix for the rth partial derivative equals zero
c then there is no need to perform the following multiplication.
do 18 j = 1, neqnssq
ur(j) = hi*ur(j)
18 continue
c
endif
c
c Add on weighted contribution (1.0-v(r))*I.
c
alpha = (r-1)*neqns**2
if ( v(r) .NE. 1.0d0 ) then
c
c If `1-v(r) = 0' then there is no need to perform the following
c multiplication.
c
if ( dkdyi_zero(r) ) then
matrix_mult_dkdyi = .FALSE.
c Signifies that multiplying the matrix `ur' by the
c matrix `dfdy' is not necessary later on
c because the matrix dkdyi can be obtained simply
c by multiplying the scalar (1-v(r)) times the matrix dfdy.
c We proceed to perform the scalar multiplication here.
c
do j = 1,neqnssq
dkdyi(alpha + j) = ( 1.0d0-v(r) )*dfdy(j)
enddo
c
else
do 20 j =1, neqnssq, neqns+1
ur(j) = 1.0d0-v(r) + ur(j)
20 continue
c
endif
c
dkdyi_zero(r) = .FALSE.
c
c Indicates the matrix `ur' is non zero in at least one of the
c components and we will require a matrix multiply later.
c It also indicates that the rth partial derivative of
c `dkdyi' will be non zero in one of the components.
c
endif
c
c `ur' now contains the proper linear combination of previously
c computed partial derivatives with respect to `yi'.
c The value of dkr/dyi is obtained simply by multiplying `ur' and
c `dfdy' together. The resultant matrix value is stored in
c the `neqns*neqns' locations of `dkdyi((r-1)*neqns**2+1)'.
c
c Multiply `ur' and `dfdy'
c
if ( matrix_mult_dkdyi ) then
c
if ( .NOT. dkdyi_zero(r) ) then
c
c If the `ur' matrix for the rth partial derivative equals zero
c or rth partial derivative could be obtained by multiplying
c the scalar `(1-v(r)' by `dfdy' then there is no need to perform
c the following matrix multiplication.
c
do 26 i =1, neqns
beta = alpha + i
c
do 27 j = 1, neqns
sigma = (j-1)*neqns
sum = 0.0d0
c
do 28 m = 1,neqns
sum = sum + dfdy((m-1)*neqns+i)*ur(sigma+m)
28 continue
c
dkdyi(sigma+beta) = sum
27 continue
c
26 continue
c
else
c If the `ur' matrix is equal to zero then the corresponding rth
c partial derivative is set to zero. As well, the value of
c `dkdyi_zero(r)' will left as `.TRUE.' to avoid
c multiplying by zero in subsequent stage computations.
c
do 25 j = 1,neqnssq
dkdyi(alpha + j) = 0.0d0
25 continue
c
endif
c `(.NOT. dkdyi_zero(r) )'
c
endif
c `( matrix_mult_dkdyi) '
c
c Initialize `matrix_mult_dkdyip1' to `.TRUE.' indicating that the
c partial derivative is currently not of the form `v(r)*dfdy' and
c that multiplying the matrix `ur' by `dfdy' may still be necessary.
c
matrix_mult_dkdyip1 = .TRUE.
c
c Construct the `ur' matrix for dkr/dyip1
c
do 24 j = 1, neqnssq
ur(j) = 0.0d0
24 continue
c
c Add on weighted multiples of the previous partial derivatives.
c The `jth' previous partial derivative is stored in
c `dkdyip1((j-1)*neqns**2+1)'.
c
do 30 j = 1, r-1
c
i1 = (j-1)*s+ r
i2 = (j-1)*neqnssq
c
if ( x(i1) .NE. 0 .AND.( .NOT. dkdyip1_zero(j)) ) then
dkdyip1_zero(r) = .FALSE.
c Indicates the matrix `ur' is non zero in at least one of the
c components and we will require a matrix multiply later.
c It also indicates that the rth partial derivative of
c `dkdyip1' will be non zero in one of the components.
c
do 31 l = 1, neqnssq
ur(l) = ur(l) + x(i1)*dkdyip1(i2+l)
31 continue
c
endif
c
30 continue
c
c Multiply ur by hi
c
if ( .NOT. dkdyip1_zero(r) ) then
c
c If the `ur' matrix for the rth partial derivative equals zero
c then there is no need to perform the following multiplication.
do 35 j = 1, neqnssq
ur(j) = hi*ur(j)
35 continue
c
endif
c
c Add on weighted contribution v(r)*I
c
if ( v(r) .NE. 0.0d0 ) then
c
c If `v(r) = 0' then there is no need to perform the following
c multiplication.
c
if ( dkdyip1_zero(r) ) then
c
matrix_mult_dkdyip1 = .FALSE.
c Signifies that multiplying the matrix `ur' by the
c matrix `dfdy' is not necessary later on
c because the matrix dkdyi can be obtained simply
c by multiplying the scalar v(r) times the matrix dfdy.
c We proceed to perform the scalar multiplication here.
c
do j = 1,neqnssq
dkdyip1(alpha + j) = v(r) * dfdy(j)
enddo
c
else
do 36 j = 1, neqnssq, neqns+1
ur(j) = v(r) + ur(j)
36 continue
c
endif
c
dkdyip1_zero(r) = .FALSE.
c
c Indicates the matrix `ur' is non zero in at least one of the
c components and we will require a matrix multiply later.
c It also indicates that the rth partial derivative of
c `dkdyip1' will be non zero in one of the components.
c
endif
c
c `ur' now contains the proper linear combination of previously
c computed partial derivatives with respect to `yip1'.
c The value of dkr/dyi+1 is obtained simply by multiplying `ur'
c and `dfdy' together. The resultant matrix value is stored in
c the `neqns*neqns' locations of `dkdyip1((r-1)*neqns**2+1)'.
c
c Multiply `ur' by `dfdy'
c
if ( matrix_mult_dkdyip1 ) then
if ( .NOT. dkdyip1_zero(r) ) then
c If the `ur' matrix for the rth partial derivative equals zero
c or rth partial derivative could be obtained by multiplying
c the scalar `v(r)' by `dfdy' then there is no need to perform
c the following matrix multiplication.
c
do 39 i = 1, neqns
beta = alpha + i
c
do 38 j = 1 ,neqns
sum=0.0d0
sigma = (j-1)*neqns
c
do 37 m = 1, neqns
sum = sum + dfdy((m-1)*neqns+i)*ur(sigma+m)
37 continue
c
dkdyip1(sigma+beta) = sum
38 continue
c
39 continue
c
else
c
c If the `ur' matrix is equal to zero then the corresponding rth
c partial derivative is set to zero. As well, the value of
c `dkdyip1_zero(r)' will left as `.TRUE.' to avoid
c multiplying by zero in subsequent stage computations.
c
do 41 j = 1,neqnssq
dkdyip1(alpha + j) = 0.0d0
41 continue
c
endif
c `( .NOT. dkdyip1_zero(r) )'
c
endif
c `( matrix_mult_dkdyip1 )'
c
40 continue
c
c All partial derivatives have now been computed.
c
c COMPUTATION OF THE BLOCKS, LLi AND RRi.
c
c The left and right blocks are now constructed from linear combinations
c of the partial derivatives.
c
c Initialize LLi and RRi to 0
c
do 44 j = 1, neqnssq
LLi(j) = 0.0d0
RRi(j) = 0.0d0
44 continue
c
c Add appropriately weighted partial derivative values to LLi and RRi.
c
do 50 r = 1, s
c
i1 = (r-1)*neqnssq
c
do 51 l = 1, neqnssq
LLi(l) = LLi(l) - b(r)*dkdyi(i1+l)
RRi(l) = RRi(l) - b(r)*dkdyip1(i1+l)
51 continue
c
50 continue
c
c Multiply LLi and RRi by hi
c
do 53 j = 1, neqnssq
LLi(j) = hi*LLi(j)
RRi(j) = hi*RRi(j)
53 continue
c
c Add on -I and I, respectively
c
do 55 j= 1, neqnssq, neqns+1
LLi(j) = -1.0d0 + LLi(j)
RRi(j) = 1.0d0 + RRi(j)
55 continue
c
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine mesh_selector(neqns,Nsub,mesh_current,defect,tol,
* Nsub_star,mesh_new,info,MxNsub,s_hat)
c
c***overview***
c
c The purpose of this routine is to examine the relative defect
c associated with a continuous solution approximation in order to decide
c upon a good mesh. The `Nsub+1' points describing the current mesh are
c passed in through the vector `mesh_current'. An estimate of the blended
c relative/absolute defect, is passed in through the vector `defect',
c which has `Nsub' components, each of size `neqns', where `neqns' is
c the number of differential equations. The user defined tolerance, which
c will be used to assess the relative defect, is passed in through the
c parameter, `tol'.
c
c After examining the blended relative/absolute defect, this routine
c will call either the `half_mesh' routine, to generate a new mesh by
c halving each subinterval of the current mesh, or the routine
c `redistribute'. `redistribute' generates a new mesh with a different
c number of points than the current mesh and with a different distribution
c of points, depending how the relative defect is distributed over the
c current mesh.
c
c The redistribution process involves monitoring the mesh function
c values which are based on the defect. A mesh redistribution is performed
c if the maximum mesh function value, is sufficiently greater than
c the average value. If so, a new mesh is determined that approximately
c equidistributes the mesh function over the subintervals of the new mesh.
c The number of points in the new mesh is determined from accuracy
c considerations.
c
c The algorithm implemented here is based on [Ascher, Mattheij, and
c Russell, Chapter 9]. The number of subintervals in the new mesh is
c returned in `Nsub_star', while the points of the new mesh are returned
c in the vector `mesh_new'. This routine will fail if the prescribed number
c of new points is greater than the maximum number allowed, `MxNsub', in
c which case we set `info = -1' to signal this.
c
c `s_hat' is used here for temporary storage.
c-------------------------------------------------------------------------------
c***declaration of constants***
integer Mxs
parameter (Mxs=10)
c `Mxs' is the maximum value for the total number of stages
c of the interpolant.
c
c***declaration of factors controlling mesh redistribution***
double precision safety_factor
double precision rho
double precision upper_new_mesh
double precision lower_new_mesh
c
parameter(safety_factor = 1.3d0)
parameter(rho = 1.0d0)
parameter(upper_new_mesh = 4.0d0)
parameter(lower_new_mesh = 0.5d0)
c
c `safety_factor' a factor to increase the prediction of the appropriate
c number of subintervals needed to achieve the given
c tolerance in the new mesh.
c
c `rho' a factor which determines how much bigger the maximum mesh
c function value must be, compared to the average value, in
c order for a mesh redistribution to occur.
c
c Setting `rho=1' implies that a mesh redistribution will take
c place everytime.
c
c `upper_new_mesh', `lower_new_mesh'
c
c are upper and lower bounds on the ratio between the number
c of subintervals in current mesh and the number of subintervals
c in the new mesh.
c
c***declaration of parameters***
c imports:
integer neqns
integer Nsub
double precision mesh_current(0:Nsub)
double precision defect(neqns*Nsub)
double precision tol
integer MxNsub
c
c `neqns' the number of differential equations.
c
c `Nsub' the number of subintervals in the current mesh.
c
c `mesh_current' a vector of points defining the current mesh.
c
c `defect' a vector giving an estimate of the blended
c relative/absolute defect, with `Nsub' components,
c each of size `neqns'.
c
c `tol' the user defined tolerance.
c
c `MxNsub' is the maximum number of subintervals.
c exports:
integer Nsub_star
double precision mesh_new(0:MxNsub)
integer info
c
c `Nsub_star' the number of subintervals of the new mesh.
c
c `mesh_new' the points defining the new mesh.
c
c `info' communication flag
c
c info = -1 size of the new mesh would be too
c large
c
c info = 0 successful generation of a new mesh
c
c***declaration of work array***
double precision s_hat(Nsub)
c
c `s_hat' a vector of mesh function values, with one component per
c subinterval. The mesh function value is based on the
c value of the norm of the defect, the tolerance `tol',
c the size of the subinterval, and the order of the Runge-Kutta
c scheme, `p'.
c***declaration of local variables***
integer i
double precision hi
double precision r1, r2, r3
integer offset
integer Nsub_pred
double precision norm
c
c `i' loop index from 1 to `Nsub'.
c
c `hi' size of ith subinterval of current mesh.
c
c 1/(p+1)
c s_hat(i) = (|| defect(i) ||/tol )
c ----------------------
c hi
c
c `r1' maximum size of weighted mesh function values for any
c subinterval of current mesh.
c
c `r2' weighted sum of mesh function values, `s_hat(i)', with
c weights, `hi'. `r2' is the predicted number of points
c that should be in the new mesh.
c
c `r3' average of weighted mesh function values.
c
c `offset' offset used in the vector `defect'
c for the `neqns' solution components of the ith
c subinterval.
c
c `Nsub_pred' predicted value for number of subintervals in
c new mesh if a redistribution is needed.
c
c `norm' max norm of the defect estimate for the ith subinterval.
c
c***declaration of variables in common block /interp_order/***
c imports:
integer p
common /interp_order/ p
c
c `p' the order of the interpolant associated with the underlying
c discrete RK formula.
c
c***declaration of variables in common block /IOCNTRL/***
c imports:
integer print_level_0, print_level_1, print_level_2
parameter (print_level_0 = 0, print_level_1 = 1,
* print_level_2 = 2)
integer profile
common /IOCNTRL/ profile
c
c `print_level_0' No output.
c `print_level_1' Intermediate output.
c `print_level_2' Full output.
c `profile' Controls output, to standard output, of profiling infor-
c mation such as Newton iteration counts, mesh selection,
c relative defect estimate. This variable is identical to
c the `output_control' parameter.
c-----------------------------------------------------------------------------
c called by: `mirkdc'
c calls to : `half_mesh', `redistribute', `idamax'
integer idamax
c------------------------------------------------------------------------------
c Compute weighted mesh function values as well as r1 and r2.
c `r1' will be the maximum value of `s_hat(i)*hi'; `r2' will be the sum
c of these same values. Initialize each to zero. Also set `info' to 0.
c
info = 0
r1 = 0.0d0
r2 = 0.0d0
c
do 10 i = 1 , Nsub
c
hi = mesh_current(i)-mesh_current(i-1)
c
c Compute offset to beginning of defect estimate vector for
c ith subinterval. This vector is of length `neqns'.
c
offset = (i-1)*neqns
c
c Compute norm of the relative defect estimate vector for the
c ith subinterval.
c
norm = dabs(defect(offset+idamax(neqns,defect(offset+1),1)))
c
c Compute corresponding mesh function value, `s_hat(i)'.
c
s_hat(i) = (norm/tol)**(1.0d0/(p+1))/hi
c
c Update maximum weighted mesh function value, if necessary.
c
if ((s_hat(i)*hi) .GT. r1) r1 = s_hat(i)*hi
c
c Accumulate ith weighted mesh function value in `r2'
c
r2 = r2 + s_hat(i)*hi
c
10 continue
c
c `r3' is the average of the weighted mesh function values
c
r3 = r2/Nsub
c
c `r2' gives an estimate of the number of subintervals that the
c new mesh should have. We employ a safety factor to increase
c the number of points slightly to increase the probability
c of achieving the tolerance.
c
Nsub_pred = (safety_factor*r2)+1
c
c To avoid cycling, make sure new mesh has a few more points than
c the previous one.
c
if ( dabs( (Nsub_pred - Nsub)/(1.0d0*Nsub)) .lt. 0.1d0)
+ Nsub_pred = 1.1d0 * Nsub
c
c Decide whether to half the current mesh or produce a new mesh by
c redistributing the points of the current mesh, while adding or deleting
c a few. Redistribute if the maximum weighted mesh function value
c is more than `rho' times the size of the average value.
c
if (r1 .LE. rho*r3) then
c
c The maximum value is not sufficiently large. That is, the
c distribution of mesh function values is sufficiently uniform
c and mesh redistribution is unwarrented. In this case, the current
c mesh is halved.
c
Nsub_star = 2*Nsub
c
if (Nsub_star .GT. MxNsub) then
c
c mesh halving fails; signal this by setting info = -1
info = -1
c
if (profile .GT. print_level_0) then
write(6,*)' New mesh would be too large'
endif
c
else
c
if (profile .GT. print_level_0) then
write(6,*)' Half the current mesh'
end if
c
call half_mesh(Nsub,mesh_current,mesh_new)
c A new mesh twice as fine as the current one is returned
c in `mesh_new'.
c
endif
c
else
c
c The distribution of mesh function values is sufficiently
c non-uniform to warrent a mesh redistribution.
c
c The value of `Nsub_pred' gives the required number of
c subintervals for the new mesh. However, if this is
c too big or too small we will not use it.
c
Nsub_star = Nsub_pred
c
if ( Nsub_star .GT. upper_new_mesh * Nsub )
+ Nsub_star = upper_new_mesh * Nsub
if ( Nsub_star .LT. lower_new_mesh * Nsub )
+ Nsub_star = lower_new_mesh * Nsub
c
c Also, check that `Nsub_star' is less than the maximum value.
c
if (Nsub_star .GT. MxNsub) then
c
c mesh redistribution fails; signal this by setting info = -1
info = -1
c
if (profile .GT. print_level_0) then
write(6,*)' New mesh would be too large'
endif
c
else
c
c if (profile .GT. print_level_0) then
c write(6,*)'Call to redistibute with Nsub_star=',
c + Nsub_star
c end if
c
call redistribute(Nsub,mesh_current,s_hat,
+ Nsub_star,mesh_new)
c The new mesh, having `Nsub_star' subintervals, is returned in
c `mesh_new'
c
endif
c
endif
c `r1. LE. rho*r3'
c
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine NewIter(neqns,leftbc,Nsub,mesh,Y,newtol,maxiter,
+ info,PHI,delta,top,bot,blocks,pivot,k_discrete,
+ work,fsub,gsub,dfsub,dgsub)
c
c***overview***
c
c This routine takes a mesh of `Nsub' subintervals, partitioning
c the problem interval and a discrete initial approximation provided
c at the meshpoints, (contained in the vectors `Y' and `mesh',
c respectively) and uses a modified Newton method to produce a
c solution approximation to within a tolerance of `newtol', for
c the associated discrete nonlinear system `PHI(Y) = 0', where `PHI(Y)' is
c the residual function. The number of differential equations
c is `neqns', while the number of boundary conditions at the lefthand
c endpoint of the problem interval is `leftbc'. The iteration completes
c successfully if it produces a converged solution within `maxiter'
c iterations. A successful return is indicated by `info'= 0. The
c updated solution is returned in `Y'. The arrays `PHI', `delta', `top',
c `bot', `blocks', `pivot', and `work' provide working space for the
c `resid', `damped_newt', and `fixed_jacob' routines. This routine
c provides, as a by-product, the stage values associated with the
c application of a Runge-Kutta method on each subinterval of the mesh.
c The array `k_discrete' returns these stages. The subroutines,
c `fsub', `gsub', `dfsub', and `dgsub', are provided by the user to
c define the boundary value ODE and the boundary conditions and their
c derivatives.
c
c The modified Newton iteration implemented here proceeds
c through a sequence of steps, each of which must be either a damped
c Newton step or a fixed Jacobian step. The code decides, depending
c on how well the iterates appear to be converging, which type of step
c to take each time. This routine is organized as follows. After some
c initialization, we begin a loop which executes the modified Newton
c iteration. The kind of step to be taken is controlled by the value of
c the flag, `Fixed_Jacobian', which is initially set to FALSE since the
c first step must be a full Newton step.
c
c During a damped Newton step an iterative process is used to
c choose a damping factor, `lambda', which is used to provide a damped
c update to the current iterate `Y = Y - lambda*delta' where `delta' is
c the Newton correction; `delta = inv(J)*PHI' where `PHI' is the residual
c function evaluated at `Y' and J is the corresponding Jacobian also
c evaluated at `Y'. The use of a damping factor can improve the
c convergence properties of the Newton iteration as well as allowing for
c the early detection of a divergent iteration. The update is accepted
c if it produces a sufficient reduction in the value of the natural
c criterion function (which is related to the size of the residual
c function). At the end of a successful damped Newton step, the routine
c takes a fixed Jacobian step if the step was a full Newton step
c (i.e. `lambda' = 1).
c
c The fixed Jacobian step consists of first computing a potential
c new iterate, `Y_hat = Y - delta' where `delta' is the Newton correction
c for this fixed Jacobian step. Here `delta = inv(J)*PHI', where `PHI' is
c the residual function evaluated at `Y' and where J, the Jacobian, has
c not been evaluated at `Y', but rather at some previous iterate.
c The new iterate, `Y_hat' will be acceptable if it leads to a sufficient
c reduction in the size of the natural criterion function. If the new
c iterate value is accepted then another fixed Jacobian step is attempted.
c Otherwise, the code returns to taking full Newton steps.
c
c This routine can fail in three ways. If the Newton iteration takes
c too many steps, the iteration halts with `info'=1. If a singular matrix is
c encountered during the factorization of the Newton matrix (during the call
c to `damped_newt'), the iteration halts with `info'=2. If it is impossible
c for the `damped_newt' routine to obtain a suitable damping factor for the
c Newton update (indicative of an ``effectively singular" Jacobian) the
c iteration halts with `info'=3.
c-----------------------------------------------------------------------------
c***declaration of constants***
integer Mxs
parameter(Mxs=10)
c
c `Mxs' Maximum number of stages of the Runge-Kutta method.
c
c***declaration of parameters***
c imports:
integer neqns, leftbc, Nsub
double precision mesh(0:Nsub)
double precision Y(neqns*(Nsub+1))
double precision newtol
integer maxiter
c
c `neqns' Number of differential equations.
c `leftbc' Number of boundary conditions at the left
c end of the problem interval.
c `Nsub' Number of subintervals into which the problem
c interval is partitioned.
c `mesh' Set of points which partition the problem interval.
c `Y' On input, the initial approximation for the Newton iteration.
c `newtol' Tolerance value to be applied to the Newton iteration,
c i.e. the iteration stops if the Newton correction
c is less than or equal to the user specified
c tolerance. A blended relative/absolute tolerance
c is applied. If `delta_j' is the Newton correction
c for, `Y_j', the jth component of `Y', then we require
c `(|delta_j| / (|Y_j|+1)) <= newtol'.
c `maxiter' Maximum number of damped Newton iterations
c to be allowed during each call to this routine.
c exports:
c double precision Y(neqns*(Nsub+1))
integer info
double precision k_discrete(Mxs*neqns*Nsub)
c
c `Y' On output, when `info' = 0, the converged solution to
c the Newton iteration.
c `info' Internal communication flag:
c 0 - successful iteration converged to
c within user specified tolerance.
c 1 - unsuccessful termination - more than `maxiter' iterations
c were taken without convergence being achieved.
c 2 - unsuccessful termination - a singular coefficient matrix
c was encountered during the attempted solution of the
c Newton system.
c 3 - unsuccessful termination - it was impossible to obtain a
c suitable damping factor for the Newton update (indicative
c of an ``effectively singular" Jacobian) or evaluation
c of natural criterion function overflowed.
c `k_discrete' The discrete stages for all subintervals. The i-th set of
c `s*neqns' locations of `k_discrete' contains the `s
c stages, each of length `neqns' corresponding to the i-th
c subinterval.
c work space:
double precision PHI(neqns*(Nsub+1))
double precision delta(neqns*(Nsub+1))
double precision top(neqns*neqns)
double precision bot(neqns*neqns)
double precision blocks(2*neqns**2 * Nsub)
integer pivot(neqns*(Nsub+1))
double precision
+ work(3*neqns*(Nsub+1)+neqns+2*neqns**2*(Mxs+1))
c
c `PHI' Residual function; it has one set of `neqns' locations per
c subinterval. The first `leftbc' locations correspond to the
c left boundary conditions; the final `neqns-leftbc' locations
c correspond to the right boundary conditions. Used as the
c right hand side for the Newton systems.
c `delta' Newton correction vector, equal to `inv(J)*PHI(Y)',
c where J is the Jacobian evaluated at the current iterate,
c `Y', when a damped Newton step is taken.
c `top' Storage of derivative information for left boundary
c conditions.
c `bot' Storage of derivative information for right boundary
c conditions.
c `blocks' Storage of derivative information for residual function.
c `pivot' Pivoting information associated with the factorization
c of J. Together `top', `blocks', and `bot' contain
c the almost block diagonal matrix representing the
c coefficient matrix of the Newton system.
c `work' Work array passed to `resid' and `fixed_jacob'.
c
c***user-supplied subroutines***
external fsub,gsub,dfsub,dgsub
c
c `fsub' Defines f(t,y) for the first order system
c of differential equations, y' = f(t,y).
c `gsub' Defines the boundary condition equations.
c `dfsub' Defines the Jacobian of the system of
c differential equations.
c `dgsub' Defines the Jacobian of the boundary conditions.
c
c***declaration of local variables***
c
logical convrg
integer count
double precision delta_0_norm, lambda, g
logical Fixed_Jacobian
integer i_delta_0, i_work_d
integer i_Y_hat, i_PHI_hat, i_work_f
c
c convrg' Used to check for convergence of the Newton iteration.
c count' Counts number of damped Newton iterations.
c delta_0_norm' Norm of the Newton correction vector `delta_0'.
c lambda' Damping factor used to perform a damped Newton update.
c g' Value of natural criterion function.
c Fixed_Jacobian' .TRUE. if the next step to be taken should be
c a fixed Jacobian step. .FALSE. if the next step
c to be taken should be a damped Newton step.
c i_delta_0' Index into `work' for the start of the part of
c `work' that will correspond to the temporary
c storage array called `delta_0' within `damped_newt'.
c i_work_d' Index to `work' for the start of the part of
c `work' that will correspond to the work space used
c within `damped_newt'.
c i_Y_hat' Index to `work' for the start of the part of
c `work' that will correspond to the temporary storage
c array `Y_hat' used within `fixed_jacob'.
c i_PHI_hat' Index to `work' for the start of the part of
c `work' that will correspond to the temporary storage
c array `PHI_hat' used within `fixed_jacob'.
c i_work_f' Index to `work' for the start of the part of
c `work' that will correspond to the work space used
c within `fixed_jacob'.
c
c***declaration of variables for /IOCNTRL/ common block***
c imports:
integer print_level_0, print_level_1, print_level_2
parameter (print_level_0 = 0, print_level_1 = 1,
* print_level_2 = 2)
integer profile
common /IOCNTRL/ profile
c
c `print_level_0' No output.
c `print_level_1' Intermediate output.
c `print_level_2' Full output.
c `profile' Controls output, to standard output, of profiling infor-
c mation such as Newton iteration counts, mesh selection,
c and defect estimates.
c-----------------------------------------------------------------------------
c Called by: `mirkdc'
c Calls to: `fixed_jacob',`damped_newt', `resid'
c-----------------------------------------------------------------------------
c***initialization of work array indexes
i_delta_0 = 1
i_work_d = neqns*(Nsub+1) + 1
i_Y_hat = 1
i_PHI_hat = neqns*(Nsub+1) + 1
i_work_f = 2*neqns*(Nsub+1) + 1
c
c***initialization for Newton iteration***
c The number of damped Newton iterations (equal to the number of
c Jacobian evaluations and factorizations) will be counted using the
c variable `count' which is initialized to 0. The variable `convrg',
c monitors whether the iteration has converged or not and is
c initialized to .FALSE..
c
count = 0
convrg = .FALSE.
c
c The first iteration step must be a damped Newton step.
c This is handled by setting `Fixed_Jacobian' to .FALSE.
c
Fixed_Jacobian = .FALSE.
c
c When a damped Newton step is used, it is assumed that the current
c iterate, `Y', and also the value of, `PHI', the residual function
c at `Y' will be available from the previous iteration step. For the
c first step, which will be a damped Newton step, `Y' will be
c available, but the residual function must be computed here in
c preparation for the first time through the damped Newton iteration
c step. This is done by calling `resid' which will return with the
c residual function stored in the vector `PHI'.
c
call resid(neqns, leftbc, Nsub, mesh, Y, PHI,
+ k_discrete, work, fsub, gsub)
c
c In order to predict the new damping factor, each damped Newton step
c also assumes several other quantities are available from the previous
c iteration step. For the first step, these will not be available. In
c order for the prediction algorithm to detect this situation we set
c `lambda=0' (The value of `lambda' can never reach 0 during the
c usual execution of the algorithm).
c
lambda = 0.0d0
c
c***Newton iteration loop***
c REPEAT-UNTIL LOOP (Implemented here using an if-goto pair)
1 continue
if (Fixed_Jacobian) then
c
if (profile .GT. print_level_1) then
write(6,*)'Fixed Jacobian step.'
end if
c
call fixed_jacob(neqns,leftbc,Nsub,mesh,Y,newtol,
+ delta,g,PHI,top,bot,blocks,pivot,Fixed_Jacobian,
+ lambda, convrg,info, k_discrete,work(i_Y_hat),
+ work(i_PHI_hat),work(i_work_f), fsub,gsub)
c
else
c
c Update `count' to monitor the number of damped Newton iterations
count = count + 1
c
if (profile .GT. print_level_1) then
write(6,100)'Damped Newton step, count =',count,'.'
100 format(1x,a27,i3,a1)
end if
c
call damped_newt(neqns,leftbc,Nsub,mesh,Y,newtol,
+ lambda,PHI,top,bot,blocks,pivot,Fixed_Jacobian,
+ convrg,delta,delta_0_norm,g,info,k_discrete,
+ work(i_delta_0),work(i_work_d),fsub,gsub,
+ dfsub,dgsub)
endif
c End of `If (Fixed Jacobian)'
c
if ((.NOT. convrg) .AND. (info .LE. 0)) then
c Check to see that the maximum number of iterations has
c not been exceeded.
if ( count .GT. maxiter) info = 1
endif
c
c Repeat the loop unless convergence has been obtained, the allowed
c number of iterations has been exceeded, a singular matrix was
c encountered during the Newton iteration, or a suitable damping
c factor has not been found (the latter three signaled by `info' > 0).
c
c UNTIL ( convrg .OR. (info .GT. 0))
if ((.NOT. convrg) .AND.
+ (info.LE. 0)) go to 1
c***end of Newton iteration loop***
c
c***conclusion***
if (profile .GT. print_level_0) then
if (count .GT. 1) then
write(6,200)'After',count,' Newton steps'
else
write(6,300)'After',count,' Newton step'
200 format(1x,a5,i3,a13)
300 format(1x,a5,i3,a12)
endif
end if
c
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine newmat(neqns,leftbc,Nsub,mesh,Y,top,blocks,bot,
+ k_discrete,work,dfsub,dgsub)
c
c***overview***
c
c This routine computes the Newton matrix for the discrete system
c modeling the BVP. The Newton matrix is the Jacobian of the residual
c function. It is evaluated at the current iterate `Y'. The Jacobian is
c returned in the three vectors `top',`bot', and `blocks'; `top' contains
c the elements of the Jacobian associated with the boundary conditions
c evaluated at the left endpoint - it consists of `leftbc' rows, each
c of size `neqns', where `leftbc' is the number of boundary conditions
c imposed at the left endpoint, and `neqns' is the number of differential
c equations; `bot' contains the corresponding information for the right
c endpoint. It consists of `neqns-leftbc' rows, each of length `neqns'.
c The vector `blocks' contains, for the i-th subinterval of the mesh which
c partitions the problem interval, as prescribed by the vector, `mesh', the
c blocks L_i and R_i, i=1,...,Nsub where `Nsub' is the number of subinter-
c vals defined by the mesh.
c
c L_i is the derivative of the i-th component of the residual function
c with respect to y_i; R_i is the derivative of the i-th component of the
c residual function with respect to y_i+1. These blocks, stored column
c by column, are placed in `blocks' one after the other in the order
c L_i, R_i, i=1,...,Nsub. Each such block is `neqns' rows, each of length
c `neqns'.
c
c This routine accesses the discrete stages for each subinterval,
c which are computed by the `subcom' routine, and saved in the array
c `k_discrete' provided to this routine through the parameter list. Thus
c a precondition for the use of this routine is that the `resid' routine
c has already been called, with the same `Y' input. `work' is a work
c array provided to `jacblk' and also used here.
c-------------------------------------------------------------------------------
c***declaration of constants***
integer Mxs
parameter(Mxs=10)
c
c `Mxs' Maximum value for the number of stages of the RK method.
c
c***declaration of parameters***
c imports:
integer neqns,leftbc,Nsub
double precision mesh(0:Nsub), Y(neqns*(Nsub+1))
double precision k_discrete(Mxs*neqns*Nsub)
c
c `neqns' Number of differential equations.
c `leftbc' Number of boundary conditions at the left end of the problem
c interval.
c `Nsub' Number of subintervals into which the problem
c interval is partitioned.
c `mesh' Set of points which partition the problem interval.
c `Y' Current discrete approximation to the solution,
c evaluated at each meshpoint.
c `k_discrete' Vector containing the discrete stages for
c all subintervals. The i-th set of `s*neqns'
c locations of `k_discrete' contains the `s'
c stages, each of length `neqns' corresponding
c to the i-th subinterval.
c exports:
double precision top(leftbc*neqns)
double precision blocks(2*neqns**2*(Nsub+1))
double precision bot((neqns-leftbc)*neqns)
c
c `top' Derivative information for left boundary conditions.
c `blocks' Derivative information for residual function.
c `bot' Derivative information for right boundary conditions.
c
c***declaration of work arrays***
double precision work(neqns+2*(Mxs+1)*neqns**2)
c
c `work' Passed for working storage to the `jacblk' routine and
c used here as the vector where the values of the
c derivatives of the boundary conditions are returned from
c the user defined routine `dgsub'.
c
c***user-supplied subroutines***
external dfsub, dgsub
c
c `dfsub' Defines the Jacobian of the system of differential equations.
c `dgsub' Defines the Jacobian of the boundary condition equations.
c
c***declaration of local variables***
integer i, j, k, neqnsq
integer i1,i2,i3,i4,n1
integer toploc, botloc, workloc
double precision h, t
c
c `i,j,k' Loop variables; `i' ranges from 1 to `Nsub'; `j,k' range from 1
c to `neqns'.
c `i1',`i2',`i3',`i4',`n1' Indexes for the variables `Y', `k_discrete', and
c `blocks' .
c `toploc', `botloc', `workloc' Offset locations in the vectors `top',
c `bot', and `work'.
c `neqnsq' Equals `neqns*neqns'.
c
c `h' Size of the i-th subinterval, `mesh(i) - mesh(i-1)'.
c `t' Value of the (i-1)st mesh point, `mesh(i-1)'.
c
c***declarations of indices for work array passing to `jacblk'***
integer i_yr, i_dfdy, i_ur, i_dkdyi, i_dkdyip1
c `i_yr' Starting location in `work': work(1), for `y_r', the argument for
c the evaluation of the `rth' stage.
c `i_dfdy' Starting location in `work' : work(neqns+1), for `dfdy', which
c during the `rth' time through the main loop below, contains the
c derivative of the right hand side of the system of ODEs,
c evaluated at the arguments of the `rth' stage of the RK scheme.
c `i_ur' Starting location in `work' : work(neqns+neqns^2+1), for `ur',
c which is used to store the intermediate values arising in the
c computation of each partial derivative.
c `i_dkdyi' Starting location in `work' : work(neqns+2*neqns^2+1), for
c `dkdyi'; the `rth' set of `neqns**2' locations of this vector
c will contain the value of the `rth' partial derivative
c with respect to `y_i'.
c `i_dkdyip1' Starting location in `work' : work(neqns+2*neqns^2+
c Mxs*neqns^2+1),
c for `dkdyip1'; the `rth' set of `neqns**2' locations of this
c vector will contain the value of the `rth' partial derivative
c with respect to `y_ip1'.
c
c***declaration of variables in common block /RK_s/***
integer s
common /RK_s/ s
c
c `s' Number of discrete stages of the Runge-Kutta formula.
c----------------------------------------------------------------------------
c Called by: `damped_newt'.
c Calls to: `jacblk', `dgsub'.
c-----------------------------------------------------------------------------
c***computation of `blocks'***
c
neqnsq = neqns**2
i_yr = 1
i_dfdy = neqns+1
i_ur = neqns+neqnsq+1
i_dkdyi = neqns+2*neqnsq+1
i_dkdyip1 = neqns+2*neqnsq+Mxs*neqnsq+1
c
c We need to compute the Jacobian blocks L and R for each subinterval.
c
c For each subinterval ...
c
i1 = s*neqns
n1= Nsub*neqns
c
do 5 i=1, Nsub
i2 = (i-1)*neqns
i3 = (i-1)*i1
i4 = 2*(i-1)*neqnsq
h = mesh(i) - mesh(i-1)
t = mesh(i-1)
c
c y_i and y_i+1, the current solution approximations at the left and
c right endpoints of the i-th subinterval, respectively, are required
c for the computation of the i-th Jacobian blocks. Y((i-1)*neqns+1)
c is y_i; Y(i*neqns+1) is y_i+1. We will also need the `s' stage
c evaluations associated with the `ith' subinterval. These are already
c available in the `ith' set of `s*neqns' locations of `k_discrete',
c `k_discrete((i-1)*s*neqns+1)'.
c
c The Jacobian blocks are passed to the appropriate places
c in `blocks.' L is stored in the neqnsq locations starting at
c location 2*(i-1)*neqnsq+1; R is stored in the neqnsq locations
c (2*i-1)*neqnsq+1.
c
call jacblk(neqns,h,t,Y(i2+1),Y(i2+neqns+1),
+ k_discrete(i3+1), blocks(i4+1), blocks(i4+neqnsq+1),
+ work(i_yr),work(i_dfdy),work(i_ur),
+ work(i_dkdyi), work(i_dkdyip1),dfsub)
c
5 continue
c
c***computation of `top' and `bot'***
c
c We now compute the parts of the Jacobian of the residual
c function associated with the boundary conditions.
c
c `Y(1)' is y_0, `Y(Nsub*neqns+1)' is y_Nsub. We call `dgsub' to get
c the boundary condition derivative information. The information is
c returned in the vector `work'. Although, `work' is viewed inside
c `dgsub' as a square matrix of size `neqns' by `neqns', it is more
c convenient to view it here in the usual columnwise fashion, i.e.
c with the columns of the square matrix laid out one after the
c other in the column vector, `work'.
c
c Set the Jacobian to zero then call `dgsub' to obtain the
c non-zero entries.
c
do 7 i = 1,neqns*neqns
work(i) = 0.0d0
7 continue
c
call dgsub(neqns,Y(1),Y(n1+1),work)
c
c The vector `top' must be assigned the values corresponding to
c the first `leftbc' rows of the square version of the boundary condition
c matrix. Since `work' contains the elements of the boundary condition
c matrix stored in a columnwise fashion, this corresponds to selecting
c the first `leftbc' elements from each column stored within `work', and
c copying them into the appropriate locations of the vector `top'.
c
c The vector `top' will also store the first `leftbc' rows of the
c boundary condition Jacobian in a columnwise fashion. Thus it will
c consist of `neqns' vectors, each of length `leftbc', representing the
c columns of the `leftbc' by `neqns' matrix corresponding to the
c derivatives of the boundary conditions at the left end of the problem
c interval.
c
c A similar situation holds for the vector `bot', which will
c contain, in a columnwise fashion, the `neqns-leftbc' by `neqns'
c matrix, corresponding to the derivatives of the boundary conditions
c at the right end of the problem interval. The variables `toploc',
c `workloc', and `botloc' are used to identify, within their respective
c vectors, the appropriate locations for the jth column, for j=1 to
c `neqns'.
c
do 25 j = 1, neqns
c
toploc = (j-1)*leftbc
workloc = (j-1)*neqns
botloc = (j-1)*(neqns-leftbc)
c
c Copy the first leftbc locations of the jth "column" of `work'
c into the leftbc locations corresponding to the jth "column" of `top'.
c
do 10 k = 1, leftbc
top(toploc+k) = work(workloc+k)
10 continue
c
c Copy the remaining `neqns-leftbc' locations of the jth "column"
c of `work' into the `neqns-leftbc' locations corresponding to the
c jth "column" of `bot'.
c
do 20 k = leftbc+1, neqns
bot(botloc+k-leftbc) = work(workloc+k)
20 continue
c
25 continue
c
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine redistribute(Nsub,mesh_current,s_hat,Nsub_star,
+ mesh_new)
c
c***overview***
c
c This routine takes as input a mesh of `Nsub' subintervals which
c partitions the problem interval, stored in the vector `mesh_current',
c and, for each subinterval of this mesh, the value of the mesh function,
c stored in the vector, `s_hat'. The mesh function, which is piecewise
c constant, i.e. has a constant value on each subinterval of the current
c mesh, associated with some measure of the relative defect on each
c subinterval (see `mesh_selector' routine for further details). The
c routine returns a new mesh, of `Nsub_star' points, which equidistributes
c the integral of the mesh function over the subintervals of the new mesh
c in the vector `mesh_new'.
c
c The equidistribution is performed as follows. The idea is to select
c new mesh points so that the area under the mesh function is equally
c distributed over the subintervals of the new mesh. The total area under
c the mesh function is first computed and then divided by the total number
c of subintervals in the new mesh, giving a quantity we call `zeta'.
c The value of the integral of the mesh function over each subinterval of
c the new mesh must equal `zeta'. The area under the mesh function for
c each subinterval of the old mesh is computed and added to a sum,
c proceeding from left to right. Whenever the sum reaches a value beyond
c `zeta', only the area of that fraction of the subinterval needed to
c reach the value of `zeta' is added. This fraction is used to decide
c where to place the new mesh point. The procedure beginning with the new
c mesh point is repeated and continued until every subinterval of the
c current mesh is processed. See [Ascher, Mattheij, Russell, Chap.9] for
c details
c
c This routine assumes that `Nsub_star' is less than or equal to the
c maximum allowable number of subintervals.
c--------------------------------------------------------------------------------
c***declaration of parameters***
c imports:
integer Nsub
double precision mesh_current(0:Nsub)
double precision s_hat(Nsub)
integer Nsub_star
c
c `Nsub' the number of subintervals in the input mesh
c
c `mesh_current' the vector containing the current mesh
c
c `s_hat' the vector containing the mesh function
c value for each subinterval.
c
c `Nsub_star' the number of points in the new mesh.
c exports:
double precision mesh_new(0:Nsub_star)
c
c `mesh_new' the meshpoints making up the new mesh
c
c***declaration of local variables***
integer i,k
double precision sum
double precision zeta
double precision t
double precision Integral
double precision Next_Piece
c
c `i' loop index from 1 to Nsub_star.
c
c `k' loop index from 1 to Nsub.
c
c `sum' temporary variable used to compute `zeta'.
c
c `zeta' weighted sum of mesh function values, averaged over
c number of subintervals in the new mesh, `Nsub_star'.
c
c `t' used as marker to indicate the point within the
c problem interval to which the integration for the
c determination of new mesh points has proceeded so far.
c
c `Integral' used to store accumulation of contributions
c to the area under the mesh function for each
c subinterval of the new mesh.
c
c `Next_Piece' used to store the value of each contribution
c to `Integral' as it arises during the determination
c of the position of each new mesh point.
c------------------------------------------------------------------------------
c called by: `mesh_selector'
c------------------------------------------------------------------------------
c Compute equidistribution constant, `zeta'. It equals the
c weighted sum of the mesh function values divided by the desired
c number of mesh points in the new mesh, `Nsub_star'.
c
sum = 0.0d0
c
do 10 k = 1, Nsub
sum = sum + s_hat(k)*(mesh_current(k)-mesh_current(k-1))
10 continue
c
zeta = sum/Nsub_star
c
c Loop to determine points of new mesh.
c
c k gives index of the subinterval of the current mesh whose mesh
c function contribution is being added to the integral for the ith
c new mesh point.
c
k = 1
i = 0
mesh_new(0) = mesh_current(0)
t = mesh_current(0)
Integral = 0.0d0
c
c REPEAT
c --------
1 continue
c --------
c
Next_Piece = s_hat(k)*(mesh_current(k)-t)
c
if ((Integral+Next_Piece) .GT. zeta) then
c
c The next piece is too big. Set the next new meshpoint
c so that it allows us to include only the appropriate
c fraction of the current piece necessary to make the
c integral for this new subinterval equal to `zeta'.
c
mesh_new(i+1) = (zeta - Integral)/s_hat(k) + t
c
c Set the "point of integration" marker, `t', to this new mesh
c point value, set Integral to zero, and increment `i' in
c preparation for the determination of the next new mesh point.
c
t = mesh_new(i+1)
i = i + 1
Integral = 0
c
else
c
c The next piece to be added to the integral is still not
c large enough to make the value of the integral `zeta'
c or bigger, so just add it on to the current integral
c value and adjust the "point of integration" marker, `t'.
c
Integral = Integral + Next_Piece
t = mesh_current(k)
c
c The contribution from the kth subinterval of the current
c mesh is complete. Increment k.
c
k = k + 1
c
endif
c
c UNTIL ( k .GT. Nsub)
c ----------------------
if (k. LE. Nsub) goto1
c ----------------------
c
c Finish up the process by setting the last new mesh point value
c to the value of the old last mesh point.
c
mesh_new(Nsub_star) = mesh_current(Nsub)
c
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine resid(neqns,leftbc,Nsub,mesh,Y,PHI,
+ k_discrete,work,fsub,gsub)
c
c***overview***
c
c This routine evaluates the residual function at the current
c iterate, `Y', on the current mesh, `mesh'. The residual function value
c is returned in the vector `PHI'. The mesh has `Nsub' subintervals.
c
c The evaluation of the residual function is done in two phases:
c one involving the calculation of the components of `PHI' associated
c with the boundary conditions (`leftbc' conditions at the left boundary
c and `neqns-leftbc' at the right boundary, where `neqns' is the number of
c differential equations), and the other associated with the
c other `Nsub' components, each of size `neqns', corresponding to each
c subinterval of the current mesh. The components associated with
c the boundary conditions are obtained directly from the subroutine
c `gsub'. Each of the other components is obtained through a call
c to the `subcom' routine.
c
c In the first part of this routine we perform the necessary
c initializations. We next treat the components of `PHI' related to the
c subintervals, and lastly deal with the components corresponding to the
c boundary conditions.
c
c A by-product of the call to `subcom' routine is the return
c of the `s' stages that arise as intermediate quantities during the
c computation of the part of `PHI' associated with the i-th subinterval.
c These stage values will be used later during the call to `Newmat' and the
c defect estimation phases of this code, and are stored here in the
c `k_discrete' vector which is returned from this routine through the
c parameter list. `work' is a work array provided to the `subcom' routine.
c The subroutines, `fsub' and `gsub' are provided by the user to define the
c boundary value ODE and the boundary conditions.
c-------------------------------------------------------------------------------
c***declaration of constants***
integer Mxs
parameter (Mxs=10)
c
c `Mxs' is the maximum number of stages.
c
c***declaration of parameters***
c imports:
integer neqns,leftbc,Nsub
double precision mesh(0:Nsub), Y(neqns*(Nsub+1))
c
c `neqns' Number of differential equations.
c `leftbc' Number of boundary conditions at the left end of the problem
c interval.
c `Nsub' Number of subintervals into which the problem interval is
c partitioned.
c `mesh' Set of points which partition the problem interval.
c `Y' Current discrete approximation to the solution,
c evaluated at each meshpoint.
c
c exports:
double precision PHI(neqns*(Nsub+1))
double precision k_discrete(Mxs*neqns*Nsub)
c
c `PHI' Value of the residual function on each subinterval.
c It has one set of `neqns' locations per subinterval.
c In addition, the first `leftbc' locations correspond to the
c left boundary conditions; the final `neqns-leftbc' locations
c correspond to the right boundary conditions.
c `k_discrete' Vector containing the discrete stages for
c all subintervals. The i-th set of `s*neqns'
c locations of `k_discrete' contains the `s'
c stages, each of length `neqns' corresponding
c to the i-th subinterval.
c work arrays:
double precision work(neqns)
c
c `work' Temporary work array required by the `subcom' routine
c and also used here to receive the values of the residual
c function corresponding to the boundary conditions as
c returned from the user defined routine `gsub'.
c
c***declaration of variables in common block /RK_s/***
integer s
common /RK_s/ s
c
c `s' Number of discrete stages of the Runge-Kutta formula.
c
c***user-supplied subroutines***
external fsub, gsub
c
c `fsub' Defines f(t,y) for the first order
c system of differential equations, y' = f(t,y).
c `gsub' Defines the boundary condition equations.
c
c***declaration of local variables***
integer i, j, i1, n1
double precision h, t
c
c `i,j' Loop variables; `i' ranges from 1 to `Nsub';
c `j' ranges from 1 to `neqns'.
c `i1','n1' Indexes used for the variables `Y', `PHI', and `k_discrete'.
c `h' Size of the i-th subinterval, `mesh(i) - mesh(i-1)'.
c `t' Value of the i-th mesh point, `mesh(i-1)'.
c-------------------------------------------------------------------------------
c Called by: `newiter','criterion', `fixed_jacob'.
c Calls: `subcom', `gsub'.
c-------------------------------------------------------------------------------
c***computation of the main part of `PHI' associated with the subintervals***
c
c For each subinterval ...
c
do 5 i=1, Nsub
i1 = (i-1)*neqns
h = mesh(i) - mesh(i-1)
t = mesh(i-1)
c
c y_i and y_ip1, the current solution approximations at the endpoints
c of the i-th subinterval, are required for the computation of the i-th
c residual component. Y((i-1)*neqns+1) is y_i, Y(i*neqns+1) is y_ip1.
c
c We now compute the component of `PHI' associated with the i-th sub-
c interval. Note that `PHI(leftbc+(i-1)*neqns+1)' is simply
c the place in `PHI' where this information is stored. The stages
c of the discrete formula, for the i-th subinterval, are returned
c by the `subcom' routine to the appropriate part of the `k_discrete'
c vector.
c
call subcom(neqns,h,t,Y(i1+1),Y(i1+neqns+1),
+ PHI(leftbc+i1+1), k_discrete(i1*s+1),work,fsub)
c
5 continue
c
c***boundary condition case***
c
n1=Nsub*neqns
c
c Y(1) is y_0, Y(Nsub*neqns+1) is y_Nsub. We call `gsub' to get the
c boundary condition information.
c
call gsub(neqns,Y(1),Y(n1+1),work)
c
c We next assign the boundary condition info, stored in the array
c `work' to the appropriate part of `PHI'.
c
do 10 j = 1, leftbc
PHI(j) = work(j)
10 continue
c
do 15 j = leftbc+1, neqns
PHI(n1 + j) = work(j)
15 continue
c
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine RK_tableau (method)
c
c***overview****
c
c This routine sets up all coefficients related to the discrete
c RK formula. `method' defines the MIRK method to be used. The possible
c values for `method' and the corresponding MIRK scheme associated with
c each of these values are listed below.
c
c Value of `method' | MIRK formula
c -------------------------------------
c 121 | MIRK121 scheme
c 221 | MIRK221 scheme
c 232 | MIRK232 scheme (abscissa = 0,2/3)
c 2321 | MIRK232 scheme (abscissa = 1,1/3)
c 343 | MIRK343 scheme
c 453 | MIRK453 scheme
c 563 | MIRK563 scheme
c 5631 | An improved MIRK563 scheme
c
c This routine assigns values to the parameters `v',`b',`c',`x' and `s';
c thus defining the MIRK formula to be used. The tableau
c of the method is
c
c c(1) | v(1) | 0 0 ... 0
c c(2) | v(2) | x(2,1) 0 ... 0
c c(3) | v(3) | x(3,1) x(3,2) ... 0
c . . . . .
c . . . . .
c . . . . .
c c(s) | v(s) | x(s,1) x(s,2) ... 0
c --------------------------------------------------------
c b(1) b(2) ... b(s)
c
c These coefficients are then available through the
c /RK_s/,/RK_coeff/ common blocks.
c------------------------------------------------------------------------------
c***declaration of constants***
integer Mxs
parameter (Mxs=10)
c
c `Mxs' maximum number of stages in the discrete formula.
c***declaration of parameters***
c imports:
integer method
c
c `method' defines the method to be used
c
c***declaration of variables in common block /RK_s/***
integer s
common /RK_s/ s
c
c `s' is the number of discrete stages of the Runge-Kutta formula.
c***declaration of variables in common block /RK_coeff/***
double precision v(Mxs), b(Mxs)
double precision c(Mxs), x(Mxs**2)
common /RK_coeff/ v, b, c, x
c
c `v', `b', `c',and `x' are coefficients that define
c the discrete Runge-Kutta formula. See tableau above.
c
c The first component of `c' contains the abscissa for the first
c new stage. The second component of `c' contains the abscissa for
c the second new stage, and so on.
c
c The first component of `v' contains the blending coefficient
c for the first new stage. The second component of `v' contains
c the blending coefficient for the second new stage, and so on.
c
c The matrix, X, is related to the data structure, `x'
c as follows (the matrix is stored columnwise in the vector):
c
c X(i,j) = x( (j-1)*s + i )
c
c `x' contains the coupling coefficients for each of the new
c stages
c-----------------------------------------------------------------------------
c called by: `mirkdc'
c------------------------------------------------------------------------------
c Translate short form of method values.
if (method .EQ. 2) method = 221
if (method .EQ. 4) method = 343
if (method .EQ. 6) method = 563
c Beginning method definition based on value of `method'.
if (method .EQ. 121) then
c
c Define the 1-stage, 2nd order, stage order 1 MIRK formula
c
s = 1
c
c(1) = 1.0d0/2.0d0
c
v(1) = 1.0d0/2.0d0
c
b(1) = 1.0d0
c
c Column 1 of x
c
x(1) = 0.0d0
c
elseif (method .EQ. 221) then
c
c Define the 2-stage, 2nd order, stage order 1 MIRK formula
c
s = 2
c
c(1) = 0.0d0
c(2) = 1.0d0
c
v(1) = 0.0d0
v(2) = 1.0d0
c
b(1) = 1.0d0/2.0d0
b(2) = 1.0d0/2.0d0
c
c Column 1 of x
c
x(1) = 0.0d0
x(2) = 0.0d0
c
c Column 2 of x
c
x(3) = 0.0d0
x(4) = 0.0d0
c
elseif (method .EQ. 232) then
c
c Define the 2-stage, 3rd order, stage order 2 MIRK formula
c
s = 2
c
c(1) = 0.0d0
c(2) = 2.0d0/3.0d0
c
v(1) = 0.0d0
v(2) = 4.0d0/9.0d0
c
b(1) = 1.0d0/4.0d0
b(2) = 3.0d0/4.0d0
c
c Column 1 of x
c
x(1) = 0.0d0
x(2) = 2.0d0/9.0d0
c
c Column 2 of x
c
x(3) = 0.0d0
x(4) = 0.0d0
c
elseif (method. EQ. 2321) then
c
c Define the reflection of the above 2-stage, 3rd order, stage order 2
c MIRK formula
c
s = 2
c
c(1) = 1.0d0
c(2) = 1.0d0/3.0d0
v(1) = 1.0d0
v(2) = 5.0d0/9.0d0
b(1) = 1.0d0/4.0d0
b(2) = 3.0d0/4.0d0
c
c Column 1 of x
c
x(1) = 0.0d0
x(2) = -2.0d0/9.0d0
c
c Column 2 of x
c
x(3) = 0.0d0
x(4) = 0.0d0
c
elseif (method .EQ. 343) then
c
c Define the 3-stage, 4th order, stage order 3 MIRK formula
c
s = 3
c
c(1) = 0.0d0
c(2) = 1.0d0
c(3) = 1.0d0/2.0d0
c
v(1) = 0.0d0
v(2) = 1.0d0
v(3) = 1.0d0 / 2.0d0
c
b(1) = 1.0d0 / 6.0d0
b(2) = 1.0d0 / 6.0d0
b(3) = 2.0d0 / 3.0d0
c
c Column 1 of x
c
x(1) = 0.0d0
x(2) = 0.0d0
x(3) = 1.0d0 / 8.0d0
c
c Column 2 of x
c
x(4) = 0.0d0
x(5) = 0.0d0
x(6) = -1.0d0 / 8.0d0
c
c Column 3 of x
c
x(7) = 0.0d0
x(8) = 0.0d0
x(9) = 0.0d0
c
elseif (method. EQ. 453) then
c
c Define the 4-stage, 5th order, stage order 3 MIRK formula
c
s = 4
c
c(1) = 0.0d0
c(2) = 1.0d0
c(3) = 3.0d0/4.0d0
c(4) = 3.0d0/10.0d0
c
v(1) = 0.0d0
v(2) = 1.0d0
v(3) = 27.0d0/32.0d0
v(4) = 837.0d0/1250.0d0
c
b(1) = 5.0d0/54.0d0
b(2) = 1.0d0/14.0d0
b(3) = 32.0d0/81.0d0
b(4) = 250.0d0/567.0d0
c
c Column 1 of x
c
x(1) = 0.0d0
x(2) = 0.0d0
x(3) = 3.0d0/64.0d0
x(4) = 21.0d0/1000.0d0
c
c Column 2 of x
c
x(5) = 0.0d0
x(6) = 0.0d0
x(7) = -9.0d0/64.0d0
x(8) = 63.0d0/5000.0d0
c
c Column 3 of x
c
x(9) = 0.0d0
x(10) = 0.0d0
x(11) = 0.0d0
x(12) = -252.0d0/625.0d0
c
c Column 4 of x
c
x(13) = 0.0d0
x(14) = 0.0d0
x(15) = 0.0d0
x(16) = 0.0d0
c
elseif (method. EQ. 4531) then
c
s = 4
c
c(1) = 0.0d0
c(2) = 1.0d0
c(3) = 3.0d0/4.0d0
c(4) = 3.0d0/10.0d0
c
v(1) = 0.0d0
v(2) = 1.0d0
v(3) = 27.0d0/32.0d0
v(4) = 837.0d0/1250.0d0
c
b(1) = 5.0d0/54.0d0
b(2) = 1.0d0/14.0d0
b(3) = 32.0d0/81.0d0
b(4) = 250.0d0/567.0d0
c
c Column 1 of x
c
x(1) = 0.0d0
x(2) = 0.0d0
x(3) = 3.0d0/64.0d0
x(4) = 21.0d0/1000.0d0
c
c Column 2 of x
c
x(5) = 0.0d0
x(6) = 0.0d0
x(7) = -9.0d0/64.0d0
x(8) = 63.0d0/5000.0d0
c
c Column 3 of x
c
x(9) = 0.0d0
x(10) = 0.0d0
x(11) = 0.0d0
x(12) = -252.0d0/625.0d0
c
c Column 4 of x
c
x(13) = 0.0d0
x(14) = 0.0d0
x(15) = 0.0d0
x(16) = 0.0d0
c
elseif (method. EQ. 563) then
c
c Define the 5-stage, 6th order, stage order 3 MIRK formula
c
s = 5
c
c(1) = 0.0d0
c(2) = 1.0d0
c(3) = 1.0d0/4.0d0
c(4) = 3.0d0/4.0d0
c(5) = 1.0d0/2.0d0
c
v(1) = 0.0d0
v(2) = 1.0d0
v(3) = 5.0d0 / 32.0d0
v(4) = 27.0d0 / 32.0d0
v(5) = 1.0d0 / 2.0d0
c
b(1) = 7.0d0 / 90.0d0
b(2) = 7.0d0 / 90.0d0
b(3) = 32.0d0 / 90.0d0
b(4) = 32.0d0 / 90.0d0
b(5) = 12.0d0 / 90.0d0
c
c Column 1 of x
c
x(1) = 0.0d0
x(2) = 0.0d0
x(3) = 9.0d0 / 64.0d0
x(4) = 3.0d0 / 64.0d0
x(5) = -5.0d0 / 24.0d0
c
c Column 2 of x
c
x(6) = 0.0d0
x(7) = 0.0d0
x(8) = -3.0d0 / 64.0d0
x(9) = -9.0d0 / 64.0d0
x(10) = 5.0d0 / 24.0d0
c
c Column 3 of x
c
x(11) = 0.0d0
x(12) = 0.0d0
x(13) = 0.0d0
x(14) = 0.0d0
x(15) = 2.0d0 / 3.0d0
c
c Column 4 of x
c
x(16) = 0.0d0
x(17) = 0.0d0
x(18) = 0.0d0
x(19) = 0.0d0
x(20) = -2.0d0 / 3.0d0
c
c Column 5 of x
c
x(21) = 0.0d0
x(22) = 0.0d0
x(23) = 0.0d0
x(24) = 0.0d0
x(25) = 0.0d0
c
elseif (method. EQ. 5631) then
c
c Define the improved 5-stage, 6th order, stage order 3 MIRK formula
c
s = 5
c
c(1) = 0.0d0
c(2) = 1.0d0
c(3) = 1.0d0/2.0d0-sqrt(21.0d0)/14.0d0
c(4) = 1.0d0/2.0d0+sqrt(21.0d0)/14.0d0
c(5) = 1.0d0/2.0d0
c
v(1) = 0.0d0
v(2) = 1.0d0
v(3) = 1.0d0/2.0d0-9.0d0*sqrt(21.0d0)/98.0d0
v(4) = 1.0d0/2.0d0+9.0d0*sqrt(21.0d0)/98.0d0
v(5) = 1.0d0 / 2.0d0
c
b(1) = 1.0d0 / 20.0d0
b(2) = 1.0d0 / 20.0d0
b(3) = 49.0d0 / 180.0d0
b(4) = 49.0d0 / 180.0d0
b(5) = 16.0d0 / 45.0d0
c
c Column 1 of x
c
x(1) = 0.0d0
x(2) = 0.0d0
x(3) = 1.0d0/14.0d0+sqrt(21.0d0)/98.0d0
x(4) = 1.0d0/14.0d0-sqrt(21.0d0)/98.0d0
x(5) = -5.0d0 / 128.0d0
c
c Column 2 of x
c
x(6) = 0.0d0
x(7) = 0.0d0
x(8) = -1.0d0/14.0d0+sqrt(21.0d0)/98.0d0
x(9) = -1.0d0/14.0d0-sqrt(21.0d0)/98.0d0
x(10) = 5.0d0 / 128.0d0
c
c Column 3 of x
c
x(11) = 0.0d0
x(12) = 0.0d0
x(13) = 0.0d0
x(14) = 0.0d0
x(15) = 7.0d0*sqrt(21.0d0)/128.0d0
c
c Column 4 of x
c
x(16) = 0.0d0
x(17) = 0.0d0
x(18) = 0.0d0
x(19) = 0.0d0
x(20) = -7.0d0*sqrt(21.0d0)/128.0d0
c
c Column 5 of x
c
x(21) = 0.0d0
x(22) = 0.0d0
x(23) = 0.0d0
x(24) = 0.0d0
x(25) = 0.0d0
c
endif
c
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine sol_eval(t,z,z_prime,fail,iwork,work)
c
c***overview***
c
c This routine uses the current solution approximation, `Y',
c obtained on the mesh stored in `mesh', plus the stages of the
c underlying discrete RK method plus some new stages associated with
c the interpolant in order to evaluate the interpolant. This interpolant
c is evaluated at `t', with the value returned in `z', and the value of
c the first derivative returned in `z_prime'.
c
c The interpolant is computed as follows: Given `t', assumed to
c be somewhere in the problem interval, the subinterval of `mesh' that
c contains `t' is determined first, i.e. find i such that mesh(i-1)
c =< t < mesh(i) for t < mesh(Nsub) or i = Nsub for t = mesh(Nsub).
c The fraction of the way `t' is through the ith subinterval, a quantity
c which we call `tau', is also determined during this step. Then the
c interpolant on this ith subinterval is evaluated - the interpolant has
c the form z(mesh(i-1) + tau*hi) = yim1 + hi sum br(tau)*kr. The weights,
c br(tau), can be computed with a call to `interp_weights'. `Nsub' is
c the number of subintervals.
c
c The remaining information required to define the interpolant
c on the ith subinterval is the current solution approximation at the
c (i-1)st mesh point, `yim1', which can be obtained from `Y', and the
c stages associated with the ith subinterval. The first `s' stages,
c contained in `k_discrete', come from the discrete formula and
c remaining `s_star - s' stages, contained in `k_interp', are necessary
c for the interpolant. Each stage is a vector of length `neqns' where
c `neqns' is the number of differential equations.
c
c `k_discrete' is the vector containing the discrete stages for
c all subintervals. The ith set of `s*neqns' locations of `k_discrete'
c contains the `s' stages, each of length `neqns' corresponding to the
c ith subinterval. `k_interp' plays a similar role for the `s_star-s'
c stages associated with the interpolant.
c------------------------------------------------------------------------------
c***declaration of constants***
integer Mxs
parameter(Mxs=10)
c
c `Mxs' is the maximum value for the total number of stages
c of the interpolant.
c
c***declaration of parameters***
c imports:
double precision t
c `t' the point at which the interpolant is to be evaluated.
c exports:
c double precision z(neqns)
double precision z(*)
c
c double precision z_prime(neqns)
double precision z_prime(*)
c
integer fail
c
c `z' is the value of the interpolant at x.
c
c `z_prime' is the value of the first derivative of the
c interpolant at x.
c
c `fail' if fail = 1 then the point `t', at which the interpolant
c is to be evaluated, is outside the problem interval and
c therefore `sol_eval' will not produce a solution
c approximation. If fail = 0 then `t' is a valid meshpoint.
c
c***LAYOUT OF WORK ARRAYS***
c
c All of the input info to this routine except for the evaluation
c point `t' is contained within the `iwork' and `work' arrays.
c
integer iwork(5)
c
c iwork(1) = neqns
integer neqns
c `neqns' is the number of differential equations.
c
c iwork(2) = Nsub
integer Nsub
c `Nsub' is the number of subintervals.
c
c iwork(3) = s
integer s
c `s' is the number of discrete stages.
c
c iwork(4) = s_star
integer s_star
c `s_star' is the total number of stages required to form the
c interpolant on each subinterval. It includes all the
c stages of the discrete formula plus the additional
c stages required for the interpolant.
c
c iwork(5) = method
integer method
c `method' defines the Runge-Kutta method that was used to compute the
c solution. It is passed via the /rk_method/ common block to
c the `interp_weights' routine.
c
c
c let MnN = Mxs*neqns*Nsub.
c
c double precision work(2*Mxs*neqns*Nsub+(Nsub+1)+(Nsub+1)*neqns)
double precision work(*)
c
c work(1..MnN) = k_discrete(Mxs*neqns*Nsub)
c `k_discrete' contains the `s' discrete stage values for each
c subinterval.
c work(MnN+1..2*MnN) = k_interp(Mxs*neqns*Nsub)
c `k_interp' contains the `s_star-s' extra stage values for each
c subinterval need to form the interpolant.
c
c work(2*MnN+1..2*MnN+(Nsub+1)) = mesh(0:Nsub)
c `mesh' is the list of points defining the subintervals
c
c work(2*MnN+(Nsub+1)+1..2*MnN+(Nsub+1)+(Nsub+1)*neqns)) =
c Y(neqns*(Nsub+1))
c `Y' contains the discrete solution at the mesh points;
c the ith set of `neqns' locations of `Y' contains
c the solution approximation at the ith mesh point.
c
c***declaration of local variables***
integer i
double precision hi
double precision tau
double precision weights(Mxs)
double precision weights_prime(Mxs)
double precision a, b
c
c `i' index for subinterval containing `t', i ranges
c from 1 to `Nsub'
c
c `hi' the size of the ith subinterval.
c
c `tau' equal to (t-mesh(i-1))/hi
c
c `weights' contains the weight polynomials of the interpolant
c evaluated at `tau'.
c
c `weights_prime' contains the first derivatives of the weight
c polynomials of the interpolant evaluated at `tau'.
c
c `a' the left hand endpoint of the interval.
c `b' the right hand endpoint of the interval.
c
c***declaration of work array indexes***
integer i_k_discrete
integer i_k_interp
integer i_mesh, i_Y
c `i_k_discrete' = 1; `i_k_interp' = Mxs*neqns*Nsub+1;
c `i_mesh' = 2*Mxs*neqns*Nsub+1; `i_Y' = 2*Mxs*neqns*Nsub+(Nsub+1)+1
c***declaration of variables in /rk_method/ common block***
c exports:
c integer method
common /rk_method/ method
c `method' defines the method to be used
c
c This value is passed through the rk_method common block
c to the `interp_weights' routine.
c-----------------------------------------------------------------------------
c called by:`main', user's main program subsequent to the
c completion of the call to mirkdc.
c
c calls to: `interval', `interp_weights', `sum_stages'
c------------------------------------------------------------------------------
c Setup access to info in `iwork' and `work'
neqns = iwork(1)
Nsub = iwork(2)
s = iwork(3)
s_star = iwork(4)
method = iwork(5)
i_k_discrete = 1
i_k_interp = Mxs*neqns*Nsub+1
i_mesh = 2*Mxs*neqns*Nsub+1
i_Y = 2*Mxs*neqns*Nsub+(Nsub+1)+1
a = work(i_mesh)
b = work(i_mesh+(Nsub))
c
fail = 0
if ((t.LT.a).or.(t.GT.b)) then
fail = 1
else
c Call `interval' to determine the subinterval containing t.
c The index of the subinterval is returned in `i'.
c
call interval(Nsub,work(i_mesh),t,i)
c
c Compute `tau' for use in the evaluation of the weight polynomials.
c
hi = work(i_mesh+i) - work(i_mesh+(i-1))
tau = (t - work(i_mesh+(i-1)))/hi
c
c Setup the weights, evaluated at `tau'.
c
call interp_weights(s_star, tau, weights,weights_prime)
c
c We need to evaluate the interpolant and its derivative at `tau'.
c
c This is done by taking weighted sums of the stages, using the weights
c we have just computed, and the stages corresponding to the ith
c subinterval, and calling the `sum_stages' routine.
c
call sum_stages(neqns,hi,work(i_Y+(i-1)*neqns),s,
+ work(i_k_discrete+(i-1)*s*neqns),s_star,
+ work(i_k_interp+(i-1)*(s_star-s)*neqns),
+ weights, weights_prime, z, z_prime)
c
end if
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine subcom(neqns,hi,ti,yi,yip1,phii,k,y,fsub)
c
c***overview***
c
c ***Subcom Routine(Calculation of subcomponents of residual function)***
c
c The purpose of this routine is to evaluate the component
c of PHI, the residual function, associated with the i-th sub-
c interval. This component of PHI is returned through the parameter
c `phii'.
c s
c --
c phii = yip1-yi- hi * < br*kr
c --
c r=1
c
c The evaluation takes place for the current iterate `Y'.
c However, the i-th component of PHI depends only on the solution
c information directly associated with the i-th subinterval, namely
c yi and yi+1, which have been assigned to `yi' and `yip1'. Also
c the size of the i-th subinterval is assigned to `hi' and the left
c hand endpoint of the ith subinterval is assigned to `ti'.
c
c The calculation uses the s-stage MIRK formula defined by the
c values of `s', `c',`v', `b' and `x' (as passed in from the
c /RK_s/,/RK_coeff/ commom blocks), The solution approximation vectors,
c `yi' and `yip1', as well as `phii', are of size `neqns', the number
c of differential equations.
c
c During the computation of `phii', we will compute `s'
c intermediate vectors, kr, r=1..s, each of size `neqns', and called
c the stages of the formula. The stages are stored in the vector `k'
c in order, and are returned by this routine for later use during the
c call to `newmat' and during the interpolation phase.
c `y' provides a work array of size `neqns'.
c-----------------------------------------------------------------------------
c***declaration of constants***
integer Mxs
parameter (Mxs=10)
c
c `Mxs' maximum number of stages.
c
c***declaration of parameters***
c imports :
integer neqns
double precision hi, ti, yi(neqns), yip1(neqns)
c
c `neqns' is the number of differential equations
c
c `hi' is the size of the ith subinterval
c
c `ti' is the lefthand endpoint of the ith subinterval
c
c `yi' is the solution approximation value at the left endpoint
c
c `yip1' is the solution approxiamtion at the right endpoint
c exports :
double precision phii(neqns)
double precision k(Mxs*neqns)
c
c `phii' is the value of the component of the residual function
c associated with the ith subinterval
c
c `k' stores the `s' stage values k1, k2, ..., ks, each of size
c `neqns'
c
c***declaration of work array***
double precision y(neqns)
c `y' a work array for temporary storage.
c
c***user-supplied subroutines***
external fsub
c
c `fsub' - This routine defines f(t,y) for the first order system
c of differential equations, y' = f(t,y).
c
c***declaration of local variables***
integer j, r, l
integer i1,i2
double precision tr
c
c `j,r,l' are loop index variables; `j,r' range from 1 to `s'
c `l' ranges from 1 to `neqns'.
c
c `i1',`i2' are indexes used in calculating ur.
c
c `tr' point at which rth stage of Runge-Kutta method is
c evaluated.
c
c***declaration of variables in common block /RK_s/***
integer s
common /RK_s/ s
c
c `s' is the number of discrete stages of the Runge-Kutta formula.
c***declaration of variables in common block /RK_coeff/***
double precision v(Mxs), b(Mxs)
double precision c(Mxs), x(Mxs**2)
common /RK_coeff/ v, b, c, x
c
c `v', `b', `c',and `x' are coefficients that define
c the discrete Runge-Kutta formula.
c
c The first component of `c' contains the abscissa for the first
c new stage. The second component of `c' contains the abscissa for
c the second new stage, and so on.
c
c The first component of `v' contains the blending coefficient
c for the first new stage. The second component of `v' contains
c the blending coefficient for the second new stage, and so on.
c
c The matrix, X, is related to the data structure, `x'
c as follows (the matrix is stored columnwise in the vector):
c
c X(i,j) = x( (j-1)*s + i )
c
c `x' contains the coupling coefficients for each of the new
c stages
c------------------------------------------------------------------------------
c called by: `resid'
c calls to: `fsub'
c------------------------------------------------------------------------------
c***COMPUTATION OF `phii'***
c
c For each stage of the Runge-Kutta method, we compute the
c argument of the function f, the right hand side of the ODE system,
c and then evaluate f there by calling `fsub', with the results being
c stored in the appropiate part of the vector `k'. The vector `y'
c provides storage for accumulation of the argument for f.
c
do 7 r = 1, s
c
do 2 l = 1, neqns
y(l) = 0.0d0
2 continue
c
c Add the weighted sum of the previous stages
c
do 4 j = 1, r-1
i1 = (j-1)*s + r
i2 = (j-1)*neqns
c
do 3 l = 1, neqns
y(l) = y(l) + x(i1)*k(i2+l)
3 continue
c
4 continue
c
c Multiply by hi
c
do 5 l = 1, neqns
y(l) = hi*y(l)
5 continue
c Add (1-v(r))*yi + v(r)*yip1
c
do 6 l = 1, neqns
y(l) = (1.0d0 - v(r))*yi(l) + v(r)*yip1(l) + y(l)
6 continue
c
c The second argument to f has been accumulated in the vector `y'.
c Compute the abscissa for the argument to fsub.
c
tr = ti + c(r)*hi
c
c Evaluate the function f to get kr, the rth stage.
c The user routine which evaluates the right hand side of the
c system of ODE's must be called `fsub'
c
call fsub(neqns, tr, y, k((r-1)*neqns+1))
c
c The rth stage has now been computed and stored in the rth
c vector location of k
c
7 continue
c
c All `s' stages of the formula have now been computed. We can now
c take their weighted sum in order to complete the calculation of `phii'.
c
do 8 l = 1, neqns
phii(l) = 0.0d0
8 continue
c
do 9 r = 1, s
i1 = (r-1)*neqns
do 10 l = 1, neqns
phii(l) = phii(l) - b(r)*k(i1+l)
10 continue
9 continue
c
c Multiply phii by hi
c
do 11 l = 1, neqns
phii(l) = hi*phii(l)
11 continue
c
c Add (yip1 - yi) to phii
c
do 12 l = 1, neqns
phii(l) = yip1(l) - yi(l) + phii(l)
12 continue
c
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine sum_stages(neqns,hi,yim1,s,ki_discrete,
+ s_star,ki_interp,weights,weights_prime,z,z_prime)
c
c***overview***
c
c Given the stages and corresponding weights, this routine
c computes the value of the interpolant, `z', and its first derivative,
c `z_prime', as follows:
c z = yim1 + hi sum weights(r)*ki(r),
c z_prime = sum weights_prime(r)*ki(r),
c where the ki(r)'s are the stages. The first `s' stages, contained
c in `ki_discrete', come from the discrete formula; the remaining
c `s_star - s' stages, contained in `ki_interp', are necessary for
c the interpolant only. Each stage is a vector of length `neqns' where
c `neqns' is the number of differential equations. `hi' is the size of
c the ith subinterval and `yim1' is the discrete solution at the left
c endpoint of the ith subinterval. `weights' contains the values of the
c weight polynomials of the interpolant and `weights_prime' contains the
c values of the first derivative of the weight polynomials.
c------------------------------------------------------------------------------
c***declaration of parameters***
c imports:
integer neqns
double precision hi
double precision yim1(neqns)
integer s
double precision ki_discrete(s*neqns)
integer s_star
double precision ki_interp((s_star-s)*neqns)
double precision weights(s_star)
double precision weights_prime(s_star)
c
c `neqns' is the number of differential equations.
c
c `hi' is the size of the current (ith) subinterval.
c
c `yim1' is the discrete solution at the left endpoint
c of the ith subinterval.
c
c `s' is the number of discrete stages.
c
c `ki_discrete' the set of `s' discrete stages associated with
c the ith subinterval.
c
c `s_star' is the total number of stages.
c
c `ki_interp' the set of `s_star-s' interpolation stages
c associated with the ith subinterval.
c
c `weights' the weights associated with the stages needed to
c evaluate the interpolant at the appropriate
c point. The weights have aleady been evaluated at
c that point.
c
c `weights_prime' the weights associated with the stages needed to
c form the first derivative of the interpolant
c at the appropriate point. The weights have already
c been evaluated at this point.
c exports:
double precision z(neqns)
double precision z_prime(neqns)
c
c `z' is the value of the interpolant.
c
c `z_prime' is the value of the first derivative of the
c interpolant.
c
c***declaration of local variables***
integer j
c
c `j' loop index from 1 to 's', or 1 to 's_star-s'
c---------------------------------------------------------------------------
c Called by: `defect_estimate', `interp_eval'
c Calls to : `dcopy', `daxpy'
c---------------------------------------------------------------------------
c Initialize all vectors to be used in this computation to zero.
c
call dcopy(neqns,0.0d0,0,z,1)
call dcopy(neqns,0.0d0,0,z_prime,1)
c
c We can now sum the appropriately weighted stages to compute
c the interpolant `z' and the first derivative of the
c interpolant, `z_prime'.
c
c Accumulate weighted sum of discrete stages for `z' and `z_prime'.
c
do 10 j = 1, s
c
call daxpy(neqns,weights(j),
+ ki_discrete((j-1)*neqns+1),1,z,1)
call daxpy(neqns,weights_prime(j),
+ ki_discrete((j-1)*neqns+1),1,z_prime,1)
10 continue
c
c Accumulate weighted sum of new interpolant stages for `z'
c and `z_prime'.
c
do 20 j = 1,(s_star-s)
c
call daxpy(neqns,weights(s+j),
+ ki_interp((j-1)*neqns+1),1,z,1)
call daxpy(neqns,weights_prime(s+j),
+ ki_interp((j-1)*neqns+1),1,z_prime,1)
20 continue
c
c The calculation of `z_prime' is now complete.
c Scale the weighted sum of stages by hi, and then add
c on yim1. This will complete the calculation of `z'.
c
call dscal(neqns,hi,z,1)
call daxpy(neqns,1.0d0,yim1,1,z,1)
c
return
end