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', whic