0001: /*
0002:  *  dverk.c
0003:  *  Dverk-Test
0004:  *
0005:  *  Created by ashley on 15/05/2009.
0006:  *  Copyright 2009 __MyCompanyName__. All rights reserved.
0007:  *
0008:  */
0009: 
0010: #include "control.h"
0011: #include "dverk.h"
0012: #include <stdio.h>
0013: #include <stdlib.h>
0014: #include <math.h>
0015: 
0016: void discrete(dverkcomm *c, void (*fcn)(), double x, double *y, double xend, double tol, int *ind, int nw, void *fdata);
0017: void dverk(dverkcomm *c, void (*fcn)(), double x, double *y, double xend, double tol, int *ind, int nw, void *fdata);
0018: 
0019: double dmax1(double x, double y)
0020: {
0021:         double z = (x>y)?x:y;
0022:         return z;
0023: }
0024: 
0025: double dmin1(double x, double y)
0026: {
0027:         double z = (x<y)?x:y;
0028:         return z;
0029: }
0030: 
0031: double dsign(double x, double y)
0032: {
0033:         double z;
0034:         if (y>0) {
0035:                 z = fabs(x);
0036:         } else {
0037:                 z = -fabs(x);
0038:         }
0039:         return z;
0040: }
0041: 
0042: dverkcomm *makeDverkCommunication(control *c)
0043: {
0044:         int n = c->nvariables;
0045:         dverkcomm *comm = (dverkcomm *)calloc(1, sizeof(dverkcomm));
0046:         comm->n = n;
0047:         if (c->integrator==discreteintegrator) {
0048:                 comm->solver = discrete;
0049:         } else if (c->integrator==dverkintegrator) {
0050:                 comm->floorvalues = calloc(n, sizeof(double));
0051:                 comm->solver = dverk;
0052:         } else {
0053:         }
0054:         return comm;
0055: }
0056: 
0057: void freeDverkCommunication(dverkcomm *comm)
0058: {
0059:         free(comm->floorvalues);
0060:         freecontrol(comm->ctrl);
0061:         free(comm);
0062: }
0063: 
0064: void discrete(dverkcomm *c, void (*fcn)(), double x, double *y, double xend, double tol, int *ind, int nw, void *fdata)
0065: {
0066: }
0067: 
0068: void euler(dverkcomm *c, void (*fcn)(), double x, double *y, double xend, double tol, int *ind, int nw, void *fdata)
0069: {
0070:         double t = x;
0071:         int n = c->n;
0072:         double tout = x + c->ctrl->dt;
0073:         double tstep = c->hstart;
0074:         double *w0 = (double *)calloc(n, sizeof(double));
0075: 
0076:         c->ctrl->output(c->ctrl, n, y, t);
0077: 
0078:         while (t<=xend) {
0079:                 fcn(n, x, y, w0, fdata);
0080:                 int j;
0081:                 for(j=0; j<n; j++) {
0082:                         y[j] += tstep*w0[j];
0083:                 }
0084:                 t += tstep;
0085: 
0086:                 while ((x>tout) && (c->ctrl->dt>0.0)) {
0087:                         c->ctrl->output(c->ctrl, n, y, x);
0088:                         tout += c->ctrl->dt;
0089:                 }
0090:         }
0091:         
0092:         free(w0);
0093: }
0094: 
0095: //      subroutine dverk (n, fcn, x, y, xend, tol, ind, c, nw, w)
0096: //      integer n, ind, nw, k
0097: //      double precision x, y(n), xend, tol, c(24), w(nw,9), temp
0098: 
0099: /*
0100: c***********************************************************************
0101: c                                                                      *
0102: c note added 11/14/85.                                                 *
0103: c                                                                      *
0104: c if you discover any errors in this subroutine, please contact        *
0105: c                                                                      *
0106: c        kenneth r. jackson                                            *
0107: c        department of computer science                                *
0108: c        university of toronto                                         *
0109: c        toronto, ontario,                                             *
0110: c        canada   m5s 1a4                                              *
0111: c                                                                      *
0112: c        phone: 416-978-7075                                           *
0113: c                                                                      *
0114: c        electronic mail:                                              *
0115: c        uucp:   {cornell,decvax,ihnp4,linus,uw-beaver}!utcsri!krj     *
0116: c        csnet:  krj@toronto                                           *
0117: c        arpa:   krj.toronto@csnet-relay                               *
0118: c        bitnet: krj%toronto@csnet-relay.arpa                          *
0119: c                                                                      *
0120: c dverk is written in fortran 66.                                      *
0121: c                                                                      *
0122: c the constants dwarf and rreb -- c(10) and c(11), respectively -- are *
0123: c set for a  vax  in  double  precision.  they  should  be  reset,  as *
0124: c described below, if this program is run on another machine.          *
0125: c                                                                      *
0126: c the c array is declared in this subroutine to have one element only, *
0127: c although  more  elements  are  referenced  in this subroutine.  this *
0128: c causes some compilers to issue warning messages.  there is,  though, *
0129: c no  error  provided  c is declared sufficiently large in the calling *
0130: c program, as described below.                                         *
0131: c                                                                      *
0132: c the following external statement  for  fcn  was  added  to  avoid  a *
0133: c warning  message  from  the  unix  f77 compiler.  the original dverk *
0134: c comments and code follow it.                                         *
0135: c                                                                      *
0136: c***********************************************************************
0137: c
0138:       external fcn
0139: c
0140: c***********************************************************************
0141: c                                                                      *
0142: c     purpose - this is a runge-kutta  subroutine  based  on  verner's *
0143: c fifth and sixth order pair of formulas for finding approximations to *
0144: c the solution of  a  system  of  first  order  ordinary  differential *
0145: c equations  with  initial  conditions. it attempts to keep the global *
0146: c error proportional to  a  tolerance  specified  by  the  user.  (the *
0147: c proportionality  depends  on the kind of error control that is used, *
0148: c as well as the differential equation and the range of integration.)  *
0149: c                                                                      *
0150: c     various options are available to the user,  including  different *
0151: c kinds  of  error control, restrictions on step sizes, and interrupts *
0152: c which permit the user to examine the state of the  calculation  (and *
0153: c perhaps make modifications) during intermediate stages.              *
0154: c                                                                      *
0155: c     the program is efficient for non-stiff systems.  however, a good *
0156: c variable-order-adams  method  will probably be more efficient if the *
0157: c function evaluations are very costly.  such a method would  also  be *
0158: c more suitable if one wanted to obtain a large number of intermediate *
0159: c solution values by interpolation, as might be the case  for  example *
0160: c with graphical output.                                               *
0161: c                                                                      *
0162: c                                    hull-enright-jackson   1/10/76    *
0163: c                                                                      *
0164: c***********************************************************************
0165: c                                                                      *
0166: c     use - the user must specify each of the following                *
0167: c                                                                      *
0168: c     n  number of equations                                           *
0169: c                                                                      *
0170: c   fcn  name of subroutine for evaluating functions - the  subroutine *
0171: c           itself must also be provided by the user - it should be of *
0172: c           the following form                                         *
0173: c              subroutine fcn(n, x, y, yprime)                         *
0174: c              integer n                                               *
0175: c              double precision x, y(n), yprime(n)                     *
0176: c                      *** etc ***                                     *
0177: c           and it should evaluate yprime, given n, x and y            *
0178: c                                                                      *
0179: c     x  independent variable - initial value supplied by user         *
0180: c                                                                      *
0181: c     y  dependent variable - initial values of components y(1), y(2), *
0182: c           ..., y(n) supplied by user                                 *
0183: c                                                                      *
0184: c  xend  value of x to which integration is to be carried out - it may *
0185: c           be less than the initial value of x                        *
0186: c                                                                      *
0187: c   tol  tolerance - the subroutine attempts to control a norm of  the *
0188: c           local  error  in  such  a  way  that  the  global error is *
0189: c           proportional to tol. in some problems there will be enough *
0190: c           damping  of  errors, as well as some cancellation, so that *
0191: c           the global error will be less than tol. alternatively, the *
0192: c           control   can   be  viewed  as  attempting  to  provide  a *
0193: c           calculated value of y at xend which is the exact  solution *
0194: c           to  the  problem y' = f(x,y) + e(x) where the norm of e(x) *
0195: c           is proportional to tol.  (the norm  is  a  max  norm  with *
0196: c           weights  that  depend on the error control strategy chosen *
0197: c           by the user.  the default weight for the k-th component is *
0198: c           1/max(1,abs(y(k))),  which therefore provides a mixture of *
0199: c           absolute and relative error control.)                      *
0200: c                                                                      *
0201: c   ind  indicator - on initial entry ind must be set equal to  either *
0202: c           1  or  2. if the user does not wish to use any options, he *
0203: c           should set ind to 1 - all that remains for the user to  do *
0204: c           then  is  to  declare c and w, and to specify nw. the user *
0205: c           may also  select  various  options  on  initial  entry  by *
0206: c           setting ind = 2 and initializing the first 9 components of *
0207: c           c as described in the next section.  he may also  re-enter *
0208: c           the  subroutine  with ind = 3 as mentioned again below. in *
0209: c           any event, the subroutine returns with ind equal to        *
0210: c              3 after a normal return                                 *
0211: c              4, 5, or 6 after an interrupt (see options c(8), c(9))  *
0212: c              -1, -2, or -3 after an error condition (see below)      *
0213: c                                                                      *
0214: c     c  communications vector - the dimension must be greater than or *
0215: c           equal to 24, unless option c(1) = 4 or 5 is used, in which *
0216: c           case the dimension must be greater than or equal to n+30   *
0217: c                                                                      *
0218: c    nw  first dimension of workspace w -  must  be  greater  than  or *
0219: c           equal to n                                                 *
0220: c                                                                      *
0221: c     w  workspace matrix - first dimension must be nw and second must *
0222: c           be greater than or equal to 9                              *
0223: c                                                                      *
0224: c     the subroutine  will  normally  return  with  ind  =  3,  having *
0225: c replaced the initial values of x and y with, respectively, the value *
0226: c of xend and an approximation to y at xend.  the  subroutine  can  be *
0227: c called  repeatedly  with new values of xend without having to change *
0228: c any other argument.  however, changes in tol, or any of the  options *
0229: c described below, may also be made on such a re-entry if desired.     *
0230: c                                                                      *
0231: c     three error returns are also possible, in which  case  x  and  y *
0232: c will be the most recently accepted values -                          *
0233: c     with ind = -3 the subroutine was unable  to  satisfy  the  error *
0234: c        requirement  with a particular step-size that is less than or *
0235: c        equal to hmin, which may mean that tol is too small           *
0236: c     with ind = -2 the value of hmin  is  greater  than  hmax,  which *
0237: c        probably  means  that the requested tol (which is used in the *
0238: c        calculation of hmin) is too small                             *
0239: c     with ind = -1 the allowed maximum number of fcn evaluations  has *
0240: c        been  exceeded,  but  this  can only occur if option c(7), as *
0241: c        described in the next section, has been used                  *
0242: c                                                                      *
0243: c     there are several circumstances that will cause the calculations *
0244: c to  be  terminated,  along with output of information that will help *
0245: c the user determine the cause of  the  trouble.  these  circumstances *
0246: c involve  entry with illegal or inconsistent values of the arguments, *
0247: c such as attempting a normal  re-entry  without  first  changing  the *
0248: c value of xend, or attempting to re-enter with ind less than zero.    *
0249: c                                                                      *
0250: c***********************************************************************
0251: c                                                                      *
0252: c     options - if the subroutine is entered with ind = 1, the first 9 *
0253: c components of the communications vector are initialized to zero, and *
0254: c the subroutine uses only default values  for  each  option.  if  the *
0255: c subroutine  is  entered  with ind = 2, the user must specify each of *
0256: c these 9 components - normally he would first set them all  to  zero, *
0257: c and  then  make  non-zero  those  that  correspond to the particular *
0258: c options he wishes to select. in any event, options may be changed on *
0259: c re-entry  to  the  subroutine  -  but if the user changes any of the *
0260: c options, or tol, in the course of a calculation he should be careful *
0261: c about  how  such changes affect the subroutine - it may be better to *
0262: c restart with ind = 1 or 2. (components 10 to 24 of c are used by the *
0263: c program  -  the information is available to the user, but should not *
0264: c normally be changed by him.)                                         *
0265: c                                                                      *
0266: c  c(1)  error control indicator - the norm of the local error is  the *
0267: c           max  norm  of  the  weighted  error  estimate  vector, the *
0268: c           weights being determined according to the value of c(1) -  *
0269: c              if c(1)=1 the weights are 1 (absolute error control)    *
0270: c              if c(1)=2 the weights are 1/abs(y(k))  (relative  error *
0271: c                 control)                                             *
0272: c              if c(1)=3 the  weights  are  1/max(abs(c(2)),abs(y(k))) *
0273: c                 (relative  error  control,  unless abs(y(k)) is less *
0274: c                 than the floor value, abs(c(2)) )                    *
0275: c              if c(1)=4 the weights are 1/max(abs(c(k+30)),abs(y(k))) *
0276: c                 (here individual floor values are used)              *
0277: c              if c(1)=5 the weights are 1/abs(c(k+30))                *
0278: c              for all other values of c(1), including  c(1) = 0,  the *
0279: c                 default  values  of  the  weights  are  taken  to be *
0280: c                 1/max(1,abs(y(k))), as mentioned earlier             *
0281: c           (in the two cases c(1) = 4 or 5 the user must declare  the *
0282: c           dimension of c to be at least n+30 and must initialize the *
0283: c           components c(31), c(32), ..., c(n+30).)                    *
0284: c                                                                      *
0285: c  c(2)  floor value - used when the indicator c(1) has the value 3    *
0286: c                                                                      *
0287: c  c(3)  hmin specification - if not zero, the subroutine chooses hmin *
0288: c           to be abs(c(3)) - otherwise it uses the default value      *
0289: c              10*max(dwarf,rreb*max(weighted norm y/tol,abs(x))),     *
0290: c           where dwarf is a very small positive  machine  number  and *
0291: c           rreb is the relative roundoff error bound                  *
0292: c                                                                      *
0293: c  c(4)  hstart specification - if not zero, the subroutine  will  use *
0294: c           an  initial  hmag equal to abs(c(4)), except of course for *
0295: c           the restrictions imposed by hmin and hmax  -  otherwise it *
0296: c           uses the default value of hmax*(tol)**(1/6)                *
0297: c                                                                      *
0298: c  c(5)  scale specification - this is intended to be a measure of the *
0299: c           scale of the problem - larger values of scale tend to make *
0300: c           the method more reliable, first  by  possibly  restricting *
0301: c           hmax  (as  described  below) and second, by tightening the *
0302: c           acceptance requirement - if c(5) is zero, a default  value *
0303: c           of  1  is  used.  for  linear  homogeneous  problems  with *
0304: c           constant coefficients, an appropriate value for scale is a *
0305: c           norm  of  the  associated  matrix.  for other problems, an *
0306: c           approximation to  an  average  value  of  a  norm  of  the *
0307: c           jacobian along the trajectory may be appropriate           *
0308: c                                                                      *
0309: c  c(6)  hmax specification - four cases are possible                  *
0310: c           if c(6).ne.0 and c(5).ne.0, hmax is taken to be            *
0311: c              min(abs(c(6)),2/abs(c(5)))                              *
0312: c           if c(6).ne.0 and c(5).eq.0, hmax is taken to be  abs(c(6)) *
0313: c           if c(6).eq.0 and c(5).ne.0, hmax is taken to be            *
0314: c              2/abs(c(5))                                             *
0315: c           if c(6).eq.0 and c(5).eq.0, hmax is given a default  value *
0316: c              of 2                                                    *
0317: c                                                                      *
0318: c  c(7)  maximum number of function evaluations  -  if  not  zero,  an *
0319: c           error  return with ind = -1 will be caused when the number *
0320: c           of function evaluations exceeds abs(c(7))                  *
0321: c                                                                      *
0322: c  c(8)  interrupt number  1  -  if  not  zero,  the  subroutine  will *
0323: c           interrupt   the  calculations  after  it  has  chosen  its *
0324: c           preliminary value of hmag, and just before choosing htrial *
0325: c           and  xtrial  in  preparation for taking a step (htrial may *
0326: c           differ from hmag in sign, and may  require  adjustment  if *
0327: c           xend  is  near) - the subroutine returns with ind = 4, and *
0328: c           will resume calculation at the point  of  interruption  if *
0329: c           re-entered with ind = 4                                    *
0330: c                                                                      *
0331: c  c(9)  interrupt number  2  -  if  not  zero,  the  subroutine  will *
0332: c           interrupt   the  calculations  immediately  after  it  has *
0333: c           decided whether or not to accept the result  of  the  most *
0334: c           recent  trial step, with ind = 5 if it plans to accept, or *
0335: c           ind = 6 if it plans to reject -  y(*)  is  the  previously *
0336: c           accepted  result, while w(*,9) is the newly computed trial *
0337: c           value, and w(*,2) is the unweighted error estimate vector. *
0338: c           the  subroutine  will  resume calculations at the point of *
0339: c           interruption on re-entry with ind = 5 or 6. (the user  may *
0340: c           change ind in this case if he wishes, for example to force *
0341: c           acceptance of a step that would otherwise be rejected,  or *
0342: c           vice versa. he can also restart with ind = 1 or 2.)        *
0343: c                                                                      *
0344: c***********************************************************************
0345: c                                                                      *
0346: c  summary of the components of the communications vector              *
0347: c                                                                      *
0348: c     prescribed at the option       determined by the program         *
0349: c           of the user                                                *
0350: c                                                                      *
0351: c                                    c(10) rreb(rel roundoff err bnd)  *
0352: c     c(1) error control indicator   c(11) dwarf (very small mach no)  *
0353: c     c(2) floor value               c(12) weighted norm y             *
0354: c     c(3) hmin specification        c(13) hmin                        *
0355: c     c(4) hstart specification      c(14) hmag                        *
0356: c     c(5) scale specification       c(15) scale                       *
0357: c     c(6) hmax specification        c(16) hmax                        *
0358: c     c(7) max no of fcn evals       c(17) xtrial                      *
0359: c     c(8) interrupt no 1            c(18) htrial                      *
0360: c     c(9) interrupt no 2            c(19) est                         *
0361: c                                    c(20) previous xend               *
0362: c                                    c(21) flag for xend               *
0363: c                                    c(22) no of successful steps      *
0364: c                                    c(23) no of successive failures   *
0365: c                                    c(24) no of fcn evals             *
0366: c                                                                      *
0367: c  if c(1) = 4 or 5, c(31), c(32), ... c(n+30) are floor values        *
0368: c                                                                      *
0369: c***********************************************************************
0370: c                                                                      *
0371: c  an overview of the program                                          *
0372: c                                                                      *
0373: c     begin initialization, parameter checking, interrupt re-entries   *
0374: c  ......abort if ind out of range 1 to 6                              *
0375: c  .     cases - initial entry, normal re-entry, interrupt re-entries  *
0376: c  .     case 1 - initial entry (ind .eq. 1 or 2)                      *
0377: c  v........abort if n.gt.nw or tol.le.0                               *
0378: c  .        if initial entry without options (ind .eq. 1)              *
0379: c  .           set c(1) to c(9) equal to zero                          *
0380: c  .        else initial entry with options (ind .eq. 2)               *
0381: c  .           make c(1) to c(9) non-negative                          *
0382: c  .           make floor values non-negative if they are to be used   *
0383: c  .        end if                                                     *
0384: c  .        initialize rreb, dwarf, prev xend, flag, counts            *
0385: c  .     case 2 - normal re-entry (ind .eq. 3)                         *
0386: c  .........abort if xend reached, and either x changed or xend not    *
0387: c  .        re-initialize flag                                         *
0388: c  .     case 3 - re-entry following an interrupt (ind .eq. 4 to 6)    *
0389: c  v        transfer control to the appropriate re-entry point.......  *
0390: c  .     end cases                                                  .  *
0391: c  .  end initialization, etc.                                      .  *
0392: c  .                                                                v  *
0393: c  .  loop through the following 4 stages, once for each trial step .  *
0394: c  .     stage 1 - prepare                                          .  *
0395: c***********error return (with ind=-1) if no of fcn evals too great .  *
0396: c  .        calc slope (adding 1 to no of fcn evals) if ind .ne. 6  .  *
0397: c  .        calc hmin, scale, hmax                                  .  *
0398: c***********error return (with ind=-2) if hmin .gt. hmax            .  *
0399: c  .        calc preliminary hmag                                   .  *
0400: c***********interrupt no 1 (with ind=4) if requested.......re-entry.v  *
0401: c  .        calc hmag, xtrial and htrial                            .  *
0402: c  .     end stage 1                                                .  *
0403: c  v     stage 2 - calc ytrial (adding 7 to no of fcn evals)        .  *
0404: c  .     stage 3 - calc the error estimate                          .  *
0405: c  .     stage 4 - make decisions                                   .  *
0406: c  .        set ind=5 if step acceptable, else set ind=6            .  *
0407: c***********interrupt no 2 if requested....................re-entry.v  *
0408: c  .        if step accepted (ind .eq. 5)                              *
0409: c  .           update x, y from xtrial, ytrial                         *
0410: c  .           add 1 to no of successful steps                         *
0411: c  .           set no of successive failures to zero                   *
0412: c**************return(with ind=3, xend saved, flag set) if x .eq. xend *
0413: c  .        else step not accepted (ind .eq. 6)                        *
0414: c  .           add 1 to no of successive failures                      *
0415: c**************error return (with ind=-3) if hmag .le. hmin            *
0416: c  .        end if                                                     *
0417: c  .     end stage 4                                                   *
0418: c  .  end loop                                                         *
0419: c  .                                                                   *
0420: c  begin abort action                                                  *
0421: c     output appropriate  message  about  stopping  the  calculations, *
0422: c        along with values of ind, n, nw, tol, hmin,  hmax,  x,  xend, *
0423: c        previous xend,  no of  successful  steps,  no  of  successive *
0424: c        failures, no of fcn evals, and the components of y            *
0425: c     stop                                                             *
0426: c  end abort action                                                    *
0427: c                                                                      *
0428: c***********************************************************************
0429: */
0430: 
0431: void dverk(dverkcomm *c, void (*fcn)(), double x, double *y, double xend, double tol, int *ind, int nw, void *fdata)
0432: {
0433: //     ******************************************************************
0434: //     * begin initialization, parameter checking, interrupt re-entries *
0435: //     ******************************************************************
0436: //
0437: //  ......abort if ind out of range 1 to 6
0438: 
0439:         int k;
0440:         double temp;
0441: //      int returnflag = 0;
0442:         int n = c->n;
0443:         double tout = x + c->ctrl->dt;
0444:         
0445:         double *w0 = (double *)calloc(n, sizeof(double));
0446:         double *w1 = (double *)calloc(n, sizeof(double));
0447:         double *w2 = (double *)calloc(n, sizeof(double));
0448:         double *w3 = (double *)calloc(n, sizeof(double));
0449:         double *w4 = (double *)calloc(n, sizeof(double));
0450:         double *w5 = (double *)calloc(n, sizeof(double));
0451:         double *w6 = (double *)calloc(n, sizeof(double));
0452:         double *w7 = (double *)calloc(n, sizeof(double));
0453:         double *w8 = (double *)calloc(n, sizeof(double));
0454: 
0455:         if ((*ind<1) || (*ind>6)) {
0456:                 goto L500;
0457:         }
0458: //        cases - initial entry, normal re-entry, interrupt re-entries
0459: //         go to (5, 5, 45, 1111, 2222, 2222), ind
0460:         switch (*ind) {
0461:                 case 1:
0462:                         goto L5;
0463:                         break;
0464:                 case 2:
0465:                         goto L5;
0466:                 case 3:
0467:                         goto L45;
0468:                 case 4:
0469:                         goto L1111;
0470:                 case 5:
0471:                         goto L2222;
0472:                 case 6:
0473:                         goto L2222;
0474:                 default:
0475:                         goto L500;
0476:         }
0477: 
0478: //        case 1 - initial entry (ind .eq. 1 or 2)
0479: //  .........abort if n.gt.nw or tol.le.0
0480: L5:     if ((n>nw) || (tol<0.0)) {
0481:                 goto L500;
0482:         }
0483:         if (*ind==2) {
0484:                 goto L15;
0485:         }
0486: //              initial entry without options (ind .eq. 1)
0487: //              set c(1) to c(9) equal to 0
0488: /*      for(k = 0; k<9; k++) {
0489:                 c[k] = 0.0e0;
0490:         } */
0491:         c->errorcontrol = 0;
0492:         c->floorvalue = 0.0;
0493:         c->hminspec = 0.0;
0494:         c->hstart = 0.0;
0495:         c->scalespec = 0.0;
0496:         c->hmaxspec = 0.0;
0497:         c->maxevals = 0;
0498:         c->interrupt1 = 0;
0499:         c->interrupt2 = 0;
0500:     goto L35;
0501: L15:
0502: //              initial entry with options (ind .eq. 2)
0503: //              make c(1) to c(9) non-negative
0504: //               do 20 k = 1, 9
0505: /*      for(k=0; k<9; k++) {
0506:         c[k] = fabs(c[k]);
0507:         } */
0508:         c->errorcontrol = abs(c->errorcontrol);
0509:         c->floorvalue =fabs(c->floorvalue);
0510:         c->hminspec = fabs(c->hminspec);
0511:         c->hstart = fabs(c->hstart);
0512:         c->scalespec = fabs(c->scalespec);
0513:         c->hmaxspec = fabs(c->hmaxspec);
0514:         c->maxevals = abs(c->maxevals);
0515:         c->interrupt1 = abs(c->interrupt1);
0516:         c->interrupt2 = abs(c->interrupt2);
0517: L20:
0518: //              make floor values non-negative if they are to be used
0519:         if ((c->errorcontrol!=4) && (c->errorcontrol!=5)) {
0520:                 goto L30;
0521:         }
0522: 
0523: //                  do 25 k = 1, n
0524:         for(k=0; k<n; k++) {
0525:                 c->floorvalues[k] = fabs(c->floorvalues[k]);
0526:         }
0527: L30:
0528: L35:
0529: //           initialize rreb, dwarf, prev xend, flag, counts
0530:         c->rreb = pow(2.0e0, -56);
0531:         c->dwarf = 1.0e-35;
0532: //           set previous xend initially to initial value of x
0533:         c->prevxend = x;
0534: //            do 40 k = 21, 24
0535: /*      for(k=20; k<24; k++) {
0536:                 c[k] = 0.0e0;
0537:         } */
0538:         c->xendflag = 0;
0539:         c->nsuccessful = 0;
0540:         c->nfailed = 0;
0541:         c->nevaluations = 0;
0542: 
0543: /*      if (c->ctrl->dt>0.0) {
0544:                 fcn(n, x, y, w0, fdata);
0545:                 c->nevaluations = c->nevaluations + 1;
0546:                 c->ctrl->output(n, y, x);
0547:                 tout += c->ctrl->dt;
0548:         } */
0549:                 
0550:         goto L50;
0551: //        case 2 - normal re-entry (ind .eq. 3)
0552: //  .........abort if xend reached, and either x changed or xend not
0553: L45:
0554:         if ((c->xendflag!=0) && (x!=c->prevxend || xend==c->prevxend)) {
0555:                 goto L500;
0556:         }
0557: //           re-initialize flag
0558:         c->xendflag = 0;
0559:         goto L50;
0560: //        case 3 - re-entry following an interrupt (ind .eq. 4 to 6)
0561: //           transfer control to the appropriate re-entry point..........
0562: //           this has already been handled by the computed go to        .
0563: //        end cases                                                     v
0564: L50:
0565: /*    end initialization, etc.
0566: c
0567: c     ******************************************************************
0568: c     * loop through the following 4 stages, once for each trial  step *
0569: c     * until the occurrence of one of the following                   *
0570: c     *    (a) the normal return (with ind .eq. 3) on reaching xend in *
0571: c     *        stage 4                                                 *
0572: c     *    (b) an error return (with ind .lt. 0) in stage 1 or stage 4 *
0573: c     *    (c) an interrupt return (with ind  .eq.  4,  5  or  6),  if *
0574: c     *        requested, in stage 1 or stage 4                        *
0575: c     ****************************************************************** */
0576: 
0577: L99999:
0578:         if ((x>tout) && (c->ctrl->dt>0.0)) {
0579:                 c->ctrl->output(c->ctrl, n, y, x);
0580:                 tout += c->ctrl->dt;
0581:         }
0582: 
0583: /*        ***************************************************************
0584: c        * stage 1 - prepare - do calculations of  hmin,  hmax,  etc., *
0585: c        * and some parameter  checking,  and  end  up  with  suitable *
0586: c        * values of hmag, xtrial and htrial in preparation for taking *
0587: c        * an integration step.                                        *
0588: c        *************************************************************** */
0589: 
0590: //      error return (with ind=-1) if no of fcn evals too great
0591:         if (c->maxevals==0 || c->nevaluations<c->maxevals) {
0592:                 goto L100;
0593:         }
0594:         *ind = -1;
0595:     goto L5000;
0596:         
0597: L100:
0598: 
0599: //      calculate slope (adding 1 to no of fcn evals) if ind .ne. 6
0600:     if (*ind==6) {
0601:                 goto L105;
0602:         }
0603:         fcn(n, x, y, w0, fdata);
0604:     c->nevaluations = c->nevaluations + 1;
0605: //      c->ctrl->output(n, y, x);
0606: //      tout += c->ctrl->dt;
0607: L105:
0608: 
0609: //      calculate hmin - use default unless value prescribed
0610:         c->hmin = c->hminspec;
0611:         if (c->hminspec != 0.0e0) goto L165;
0612: //              calculate default value of hmin
0613: //              first calculate weighted norm y - c(12) - as specified
0614: //              by the error control indicator c(1)
0615:         temp = 0.0e0;
0616:         if (c->errorcontrol != 1) goto L115;
0617: //                 absolute error control - weights are 1
0618: //                  do 110 k = 1, n
0619:         for(k=0; k<n; k++) {
0620:                 temp = (temp>fabs(y[k]))?temp:fabs(y[k]);
0621:         }
0622:         c->weightednormy = temp;
0623:         goto L160;
0624: L115:          
0625:         if (c->errorcontrol != 2) goto L120;
0626: //                 relative error control - weights are 1/dabs(y(k)) so
0627: //                 weighted norm y is 1
0628:     c->weightednormy = 1.0e0;
0629:     goto L160;
0630: L120:          
0631:         if (c->errorcontrol != 3) goto L130;
0632: //                 weights are 1/max(c(2),abs(y(k)))
0633: //                  do 125 k = 1, n
0634:         for(k=0; k<n; k++) {
0635:                 temp = dmax1(temp, fabs(y[k])/c->floorvalue);
0636:         }
0637:         c->weightednormy = dmin1(temp, 1.0e0);
0638:         goto L160;
0639: L130:          
0640:         if (c->errorcontrol != 4) goto L140;
0641: //                 weights are 1/max(c(k+30),abs(y(k)))
0642: //                  do 135 k = 1, n
0643:         for(k=0; k<n; k++) {
0644:                 temp = dmax1(temp, fabs(y[k])/c->floorvalues[k]);
0645:         }
0646:         c->weightednormy = dmin1(temp, 1.0e0);
0647:         goto L160;
0648: L140:          
0649:         if (c->errorcontrol != 5) goto L150;
0650: //                 weights are 1/c(k+30)
0651: //                  do 145 k = 1, n
0652:         for(k=0; k<n; k++) {
0653:                 temp = dmax1(temp, fabs(y[k])/c->floorvalues[k]);
0654:         }
0655:         c->weightednormy = temp;
0656:         goto L160;
0657: L150:
0658: //                 default case - weights are 1/max(1,abs(y(k)))
0659: //                  do 155 k = 1, n
0660:         for(k=0; k<n; k++) {
0661:                 temp = dmax1(temp, fabs(y[k]));
0662:         }
0663:         c->weightednormy = dmin1(temp, 1.0e0);
0664: L160:
0665:         c->hmin = 10.0e0*dmax1(c->dwarf, c->rreb*dmax1(c->weightednormy/tol,fabs(x)));
0666: L165:
0667: 
0668: //           calculate scale - use default unless value prescribed
0669:         c->scale = c->scalespec;
0670:         if (c->scalespec == 0.0e0) c->scale = 1.0e0;
0671: 
0672: //           calculate hmax - consider 4 cases
0673: //          case 1 both hmax and scale prescribed
0674:         if (c->hmaxspec!=0.0e0 && c->scalespec!=0.0e0) c->hmax = dmin1(c->hmaxspec, 2.0e0/c->scalespec);
0675: //          case 2 - hmax prescribed, but scale not
0676:         if (c->hmaxspec!=0.0e0 && c->scalespec==0.0e0) c->hmax = c->hmaxspec;
0677: //           case 3 - hmax not prescribed, but scale is
0678:         if (c->hmaxspec==0.0e0 && c->scalespec!=0.0e0) c->hmax = 2.0e0/c->scalespec;
0679: //           case 4 - neither hmax nor scale is provided
0680:         if (c->hmaxspec==0.0e0 && c->scalespec==0.0e0) c->hmax = 2.0e0;
0681: 
0682: //***********error return (with ind=-2) if hmin .gt. hmax
0683:         if (c->hmin < c->hmax) goto L170;
0684:         *ind = -2;
0685:         goto L5000;
0686: L170:
0687: 
0688: //           calculate preliminary hmag - consider 3 cases
0689:         if (*ind > 2) goto L175;
0690: //           case 1 - initial entry - use prescribed value of hstart, if
0691: //              any, else default
0692:         c->hmag = c->hstart;
0693:         if (c->hstart == 0.0e0) c->hmag = c->hmax*pow(tol, (1.0/6.0));
0694:         goto L185;
0695: L175:
0696:         if (c->nfailed > 1) goto L180;
0697: //           case 2 - after a successful step, or at most  one  failure,
0698: //              use min(2, .9*(tol/est)**(1/6))*hmag, but avoid possible
0699: //              overflow. then avoid reduction by more than half.
0700:         temp = 2.0e0*c->hmag;
0701:         if (tol < pow(2.0e0/0.9e0, 6)*c->est) temp = 0.9e0*tol/pow(c->est, (1.0/6.0))*c->hmag;
0702:         c->hmag = dmax1(temp, 0.5e0*c->hmag);
0703:         goto L185;
0704: L180:
0705: //           case 3 - after two or more successive failures
0706:         c->hmag = 0.5e0*c->hmag;
0707: L185:
0708: 
0709: //           check against hmax
0710:         c->hmag = dmin1(c->hmag, c->hmax);
0711: 
0712: //           check against hmin
0713:         c->hmag = dmax1(c->hmag, c->hmin);
0714: 
0715: //***********interrupt no 1 (with ind=4) if requested
0716:         if (c->interrupt1 == 0) goto L1111;
0717:         *ind = 4;
0718:         goto L5000;
0719: //           resume here on re-entry with ind .eq. 4   ........re-entry..
0720: L1111:
0721: 
0722: //           calculate hmag, xtrial - depending on preliminary hmag, xend
0723:         if (c->hmag >= fabs(xend - x)) goto L190;
0724: //              do not step more than half way to xend
0725:         c->hmag = dmin1(c->hmag, 0.5e0*fabs(xend - x));
0726:         c->xtrial = x + dsign(c->hmag, xend - x);
0727:     goto L195;
0728: L190:
0729: //              hit xend exactly
0730:         c->hmag = fabs(xend - x);
0731:     c->xtrial = xend;
0732: L195:
0733: 
0734: //           calculate htrial
0735:         c->htrial = c->xtrial - x;
0736: /*
0737: c        end stage 1
0738: c
0739: c        ***************************************************************
0740: c        * stage 2 - calculate ytrial (adding 7 to no of  fcn  evals). *
0741: c        * w(*,2), ... w(*,8)  hold  intermediate  results  needed  in *
0742: c        * stage 3. w(*,9) is temporary storage until finally it holds *
0743: c        * ytrial.                                                     *
0744: c        *************************************************************** */
0745: //      printf("%15.6e", c[16]);
0746: /*      printf("%8.5f", x);
0747:         for(k=0; k<n; k++) {
0748:                 printf(" %15.6e", y[k]);
0749:         }
0750:         printf("\n");  */
0751: 
0752: //      Here!
0753: /*      if ((x>tout) && (c->ctrl->dt>0.0)) {
0754:                 c->ctrl->output(c->ctrl, n, y, x);
0755:                 tout += c->ctrl->dt;
0756:         } */
0757: 
0758:         temp = c->htrial/1398169080000.0e0;
0759: 
0760: //            do 200 k = 1, n
0761:         for(k=0; k<n; k++) {
0762:                 w8[k] = y[k] + temp*w0[k]*233028180000.0e0;
0763:         }
0764:     fcn(n, x + c->htrial/6.0e0, w8, w0, fdata);
0765: 
0766: //            do 205 k = 1, n
0767:         for(k=0; k<n; k++) {
0768:                 w8[k] = y[k] + temp*(w0[k]*74569017600.0e0
0769:                                            + w1[k]*298276070400.0e0  );
0770:         }
0771:     fcn(n, x + c->htrial*(4.0e0/15.0e0), w8, w2, fdata);
0772: 
0773: //            do 210 k = 1, n
0774:         for(k=0; k<n; k++) {
0775:                 w8[k] = y[k] + temp*(w0[k]*1165140900000.0e0
0776:                                            - w1[k]*3728450880000.0e0
0777:                                            + w2[k]*3495422700000.0e0 );
0778:         }
0779:     fcn(n, x + c->htrial*(2.0e0/3.0e0), w8, w3, fdata);
0780: 
0781: //            do 215 k = 1, n
0782:         for(k=0; k<n; k++) {
0783:                 w8[k] = y[k] + temp*(-w0[k]*3604654659375.0e0
0784:                                                           + w1[k]*12816549900000.0e0 
0785:                                                           - w2[k]*9284716546875.0e0 
0786:                                                           + w3[k]*1237962206250.0e0 );
0787:         }
0788:     fcn(n, x + c->htrial*(5.0e0/6.0e0), w8, w4, fdata);
0789: 
0790: //            do 220 k = 1, n
0791:         for(k=0; k<n; k++) {
0792:                 w8[k] = y[k] + temp*(w0[k]*3355605792000.0e0 
0793:                                      - w1[k]*11185352640000.0e0 
0794:                                                          + w2[k]*9172628850000.0e0 
0795:                                                          - w3[k]*427218330000.0e0 
0796:                                                          + w4[k]*482505408000.0e0);
0797:         }
0798:     fcn(n, x + c->htrial, w8, w5, fdata);
0799: 
0800: //            do 225 k = 1, n
0801:         for(k=0; k<n; k++) {
0802:                 w8[k] = y[k] + temp*(-w0[k]*770204740536.0e0
0803:                                      + w1[k]*2311639545600.0e0
0804:                                      - w2[k]*1322092233000.0e0
0805:                                      - w3[k]*453006781920.0e0
0806:                                      + w4[k]*326875481856.0e0);
0807:         }
0808:         fcn(n, x + c->htrial/15.0e0, w8, w6, fdata);
0809: 
0810: //            do 230 k = 1, n
0811:         for(k=0; k<n; k++) {
0812:                 w8[k] = y[k] + temp*(w0[k]*2845924389000.0e0
0813:                                     - w1[k]*9754668000000.0e0
0814:                                     + w2[k]*7897110375000.0e0
0815:                                     - w3[k]*192082660000.0e0
0816:                                     + w4[k]*400298976000.0e0
0817:                                     + w6[k]*201586000000.0e0  );
0818:         }
0819:     fcn(n, x + c->htrial, w8, w7, fdata);
0820: 
0821: //           calculate ytrial, the extrapolated approximation and store
0822: //              in w(*,9)
0823: //            do 235 k = 1, n
0824:         for(k=0; k<n; k++) {
0825:                 w8[k] = y[k] + temp*(w0[k]*104862681000.0e0
0826:                                      + w2[k]*545186250000.0e0
0827:                                      + w3[k]*446637345000.0e0
0828:                                      + w4[k]*188806464000.0e0
0829:                                      + w6[k]*15076875000.0e0
0830:                                      + w7[k]*97599465000.0e0   );
0831:         }
0832: 
0833: //           add 7 to the no of fcn evals
0834:         c->nevaluations = c->nevaluations + 7;
0835: /*
0836: c        end stage 2
0837: c
0838: c        ***************************************************************
0839: c        * stage 3 - calculate the error estimate est. first calculate *
0840: c        * the  unweighted  absolute  error  estimate vector (per unit *
0841: c        * step) for the unextrapolated approximation and store it  in *
0842: c        * w(*,2).  then  calculate the weighted max norm of w(*,2) as *
0843: c        * specified by the error  control  indicator  c(1).  finally, *
0844: c        * modify  this result to produce est, the error estimate (per *
0845: c        * unit step) for the extrapolated approximation ytrial.       *
0846: c        ***************************************************************
0847: c
0848: c           calculate the unweighted absolute error estimate vector */
0849: //            do 300 k = 1, n
0850:         for(k=0; k<n; k++) {
0851:                 w1[k] = (w0[k]*8738556750.0e0
0852:                          + w2[k]*9735468750.0e0
0853:                          - w3[k]*9709507500.0e0
0854:                          + w4[k]*8582112000.0e0
0855:                          + w5[k]*95329710000.0e0
0856:                          - w6[k]*15076875000.0e0
0857:                          - w7[k]*97599465000.0e0)/1398169080000.0e0;
0858:         }
0859: 
0860: //           calculate the weighted max norm of w(*,2) as specified by
0861: //           the error control indicator c(1)
0862:         temp = 0.0e0;
0863:     if (c->errorcontrol != 1) goto L310;
0864: //              absolute error control
0865: //               do 305 k = 1, n
0866:         for(k=0; k<n; k++) {
0867:                 temp = dmax1(temp,fabs(w1[k]));
0868:         }
0869:     goto L360;
0870: L310:       
0871:         if (c->errorcontrol != 2) goto L320;
0872: //              relative error control
0873: //               do 315 k = 1, n
0874:         for(k=0; k<n; k++) {
0875:                 temp = dmax1(temp, fabs(w1[k]/y[k]));
0876:         }
0877:     goto L360;
0878: L320:
0879:         if (c->errorcontrol != 3) goto L330;
0880: //              weights are 1/max(c(2),abs(y(k)))
0881: //               do 325 k = 1, n
0882:         for(k=0; k<n; k++) {
0883:                 temp = dmax1(temp, fabs(w1[k]) / dmax1(c->floorvalue, fabs(y[k])) );
0884:         }
0885:     goto L360;
0886: L330:
0887:         if (c->errorcontrol != 4) goto L340;
0888: //              weights are 1/max(c(k+30),abs(y(k)))
0889: //               do 335 k = 1, n
0890:         for(k=0; k<n; k++) {
0891:                 temp = dmax1(temp, fabs(w1[k]) / dmax1(c->floorvalues[k], fabs(y[k])) );
0892:         }
0893:     goto L360;
0894: L340:
0895:     if (c->errorcontrol != 5) goto L350;
0896: //              weights are 1/c(k+30)
0897: //               do 345 k = 1, n
0898:         for(k=0; k<n; k++) {
0899:                 temp = dmax1(temp, fabs(w1[k]/c->floorvalues[k]));
0900:         }
0901:     goto L360;
0902: L350:
0903: //              default case - weights are 1/max(1,abs(y(k)))
0904: //               do 355 k = 1, n
0905:         for(k=0; k<n; k++) {
0906:                 temp = dmax1(temp, fabs(w1[k]) / dmax1(1.0e0, fabs(y[k])) );
0907:         }
0908: L360:
0909: 
0910: /*           calculate est - (the weighted max norm of w(*,2))*hmag*scale
0911: c              - est is intended to be a measure of the error  per  unit
0912: c              step in ytrial */
0913:         c->est = temp*c->hmag*c->scale;
0914: /*
0915: c        end stage 3
0916: c
0917: c        ***************************************************************
0918: c        * stage 4 - make decisions.                                   *
0919: c        ***************************************************************
0920: c
0921: c           set ind=5 if step acceptable, else set ind=6 */
0922:         *ind = 5;
0923:         if (c->est > tol) *ind = 6;
0924: 
0925: //**********interrupt no 2 if requested
0926:         if (c->interrupt2 == 0) goto L2222;
0927:         goto L5000;
0928: //          resume here on re-entry with ind .eq. 5 or 6   ...re-entry..
0929: L2222:
0930: 
0931:         if (*ind == 6) goto L410;
0932: //              step accepted (ind .eq. 5), so update x, y from xtrial,
0933: //                 ytrial, add 1 to the no of successful steps, and set
0934: //                 the no of successive failures to zero
0935:         x = c->xtrial;
0936: //               do 400 k = 1, n
0937:         for(k=0; k<n; k++) {
0938:                 y[k] = w8[k];
0939:         }
0940:         (c->nsuccessful)++;
0941:     c->nfailed = 0;
0942: //**************return(with ind=3, xend saved, flag set) if x .eq. xend
0943:         if (x != xend) goto L405;
0944:     *ind = 3;
0945:     c->prevxend = xend;
0946:     c->xendflag = 1;
0947:     goto L5000;
0948: L405:
0949:         goto L420;
0950: L410:
0951: //              step not accepted (ind .eq. 6), so add 1 to the no of
0952: //                 successive failures
0953:         (c->nfailed)++;// = c[22] + 1.0e0;
0954: //**************error return (with ind=-3) if hmag .le. hmin
0955:         if (c->hmag > c->hmin) goto L415;
0956:     *ind = -3;
0957:     goto L5000;
0958: L415:
0959: L420:
0960: 
0961: //        end stage 4
0962: 
0963:         goto L99999;
0964: //     end loop
0965: 
0966: //  begin abort action
0967: L500:
0968: 
0969: //      write(6,505) ind, tol, x, n, c(13), xend, nw, c(16), c(20),
0970: //     +      c(22), c(23), c(24), (y(k), k = 1, n)
0971: //
0972:         printf("Computation stopped on dverk: %4d %15.6e %15.6e %4d %15.6e %15.6e\n", *ind,tol,x,n,c->weightednormy,xend);
0973: 
0974: /*  505 format( /// 1h0, 58hcomputation stopped in dverk with the followin
0975: c     +g values -
0976: c     +   / 1h0, 5hind =, i4, 5x, 6htol  =, 1pd13.6, 5x, 11hx         =,
0977: c     +          1pd22.15
0978: c     +   / 1h , 5hn   =, i4, 5x, 6hhmin =, 1pd13.6, 5x, 11hxend      =,
0979: c     +          1pd22.15
0980: c     +   / 1h , 5hnw  =, i4, 5x, 6hhmax =, 1pd13.6, 5x, 11hprev xend =,
0981: c     +          1pd22.15
0982: c     +   / 1h0, 14x, 27hno of successful steps    =, 0pf8.0
0983: c     +   / 1h , 14x, 27hno of successive failures =, 0pf8.0
0984: c     +   / 1h , 14x, 27hno of function evals      =, 0pf8.0
0985: c     +   / 1h0, 23hthe components of y are
0986: c     +   // (1h , 1p5d24.15)                                           )*/
0987: 
0988: //505   format(///'Computation stopped on dverk:'/
0989: //     +                  i4,2d15.6,i4,2d14.5///)
0990: L5000:
0991:         if (c->ctrl->dt>0.0) {
0992:                 fcn(n, x, y, w0, fdata);
0993:                 c->nevaluations = c->nevaluations + 1;
0994:                 c->ctrl->output(c->ctrl, n, y, x);
0995: //              tout += c->ctrl->dt;
0996:         } 
0997:         
0998:         free(w0);
0999:         free(w1);
1000:         free(w2);
1001:         free(w3);
1002:         free(w4);
1003:         free(w5);
1004:         free(w6);
1005:         free(w7);
1006:         free(w8);
1007:         
1008:     return;
1009: 
1010: }