share/colnew/fortran/colnew.f

   1 c From research!csnet!CSNET-RELAY!mit-multics.arpa!UBC.mailnet!USER=NBAF Tue, 3 Feb 87 15:25:36 PST
   2 C**********************************************************************
   3 C  this package solves boundary value problems for
   4 C  ordinary differential equations, as described below.
   5 C
   6 C  COLNEW is a modification of the package COLSYS by ascher,
   7 C  christiansen and russell [1]. It incorporates a new basis
   8 C  representation replacing b-splines, and improvements for
   9 C  the linear and nonlinear algebraic equation solvers.
  10 C  the package can be referenced as either COLNEW or COLSYS.
  11 C**********************************************************************
  12 C----------------------------------------------------------------------
  13 C                            p a r t  1
  14 C        main storage allocation and program control subroutines
  15 C----------------------------------------------------------------------
  16 C
  17       SUBROUTINE COLNEW (NCOMP, M, ALEFT, ARIGHT, ZETA, IPAR, LTOL,
  18      1                   TOL, FIXPNT, ISPACE, FSPACE, IFLAG,
  19      2                   FSUB, DFSUB, GSUB, DGSUB, GUESS)
  20 C
  21 C
  22 C**********************************************************************
  23 C
  24 C     written by
  25 C                  u. ascher,
  26 C                            department of computer science,
  27 C                            university of british columbia,
  28 C                            vancouver, b. c., canada   v6t 1w5
  29 C                  g. bader,
  30 C                            institut f. angewandte mathematik
  31 C                            university of heidelberg
  32 C                            im neuenheimer feld 294
  33 C                            d-6900 heidelberg 1
  34 C
  35 C**********************************************************************
  36 C
  37 C     purpose
  38 C
  39 C     this package solves a multi-point boundary value
  40 C     problem for a mixed order system of ode-s given by
  41 C
  42 C          (m(i))
  43 C         u       =  f  ( x; z(u(x)) )      i = 1, ... ,ncomp
  44 C          i          i
  45 C
  46 C                                          aleft .lt. x .lt. aright,
  47 C
  48 C
  49 C         g  ( zeta(j); z(u(zeta(j))) ) = 0   j = 1, ... ,mstar
  50 C          j
  51 C                                    mstar = m(1)+m(2)+...+m(ncomp),
  52 C
  53 C
  54 C         where                          t
  55 C               u = (u , u , ... ,u     )  is the exact solution vector
  56 C                     1   2        ncomp
  57 C
  58 C                (mi)
  59 C               u     is the mi=m(i) th  derivative of u
  60 C                i                                      i
  61 C
  62 C                                  (1)        (m1-1)       (mncomp-1)
  63 C               z(u(x)) = ( u (x),u  (x),...,u    (x),...,u      (x) )
  64 C                            1     1          1            ncomp
  65 C
  66 C                f (x,z(u))   is a (generally) nonlinear function of
  67 C                 i
  68 C                             z(u)=z(u(x)).
  69 C
  70 C                g (zeta(j);z(u))  is a (generally) nonlinear function
  71 C                 j
  72 C                               used to represent a boundary condition.
  73 C
  74 C         the boundary points satisfy
  75 C               aleft .le. zeta(1) .le. .. .le. zeta(mstar) .le. aright
  76 C
  77 C         the orders mi of the differential equations satisfy
  78 C                            1 .le. m(i) .le. 4.
  79 C
  80 C
  81 C**********************************************************************
  82 C
  83 C     method
  84 C
  85 C        the method used to approximate the solution u is
  86 C     collocation at gaussian points, requiring m(i)-1 continuous
  87 C     derivatives in the i-th component, i = 1, ..., ncomp.
  88 C     here, k is the number of collocation points (stages) per
  89 C     subinterval and is chosen such that k .ge. max m(i).
  90 C     a runge-kutta-monomial solution representation is utilized.
  91 C
  92 C     references
  93 C
  94 C     [1] u. ascher, j. christiansen and r.d. russell,
  95 C         collocation software for boundary-value odes,
  96 C         acm trans. math software 7 (1981), 209-222.
  97 C         this paper contains EXAMPLES where use of the code
  98 C         is demonstrated.
  99 C
 100 C     [2] g. bader and u. ascher,
 101 C         a new basis implementation for a mixed order
 102 C         boundary value ode solver,
 103 C         siam j. scient. stat. comput. (1987).
 104 C
 105 C     [3] u. ascher, j. christiansen and r.d. russell,
 106 C         a collocation solver for mixed order
 107 C         systems of boundary value problems,
 108 C         math. comp. 33 (1979), 659-679.
 109 C
 110 C     [4] u. ascher, j. christiansen and r.d. russell,
 111 C         colsys - a collocation code for boundary
 112 C         value problems,
 113 C         lecture notes comp.sc. 76, springer verlag,
 114 C         b. childs et. al. (eds.) (1979), 164-185.
 115 C
 116 C     [5] c. deboor and r. weiss,
 117 C         solveblok: a package for solving almost block diagonal
 118 C         linear systems,
 119 C         acm trans. math. software 6 (1980), 80-87.
 120 C
 121 C**********************************************************************
 122 C
 123 C     ***************     input to colnew     ***************
 124 C
 125 C     variables
 126 C
 127 C     ncomp - no. of differential equations   (ncomp .le. 20)
 128 C
 129 C     m(j) - order of the j-th differential equation
 130 C            ( mstar = m(1) + ... + m(ncomp) .le. 40 )
 131 C
 132 C     aleft - left end of interval
 133 C
 134 C     aright - right end of interval
 135 C
 136 C     zeta(j) - j-th side condition point (boundary point). must
 137 C               have  zeta(j) .le. zeta(j+1). all side condition
 138 C               points must be mesh points in all meshes used,
 139 C               see description of ipar(11) and fixpnt below.
 140 C
 141 C     ipar - an integer array dimensioned at least 11.
 142 C            a list of the parameters in ipar and their meaning follows
 143 C            some parameters are renamed in colnew; these new names are
 144 C            given in parentheses.
 145 C
 146 C     ipar(1)     ( = nonlin )
 147 C             = 0 if the problem is linear
 148 C             = 1 if the problem is nonlinear
 149 C
 150 C     ipar(2) = no. of collocation points per subinterval  (= k )
 151 C               where max m(i) .le.  k .le. 7 . if ipar(2)=0 then
 152 C               colnew sets  k = max ( max m(i)+1, 5-max m(i) )
 153 C
 154 C     ipar(3) = no. of subintervals in the initial mesh  ( = n ).
 155 C               if ipar(3) = 0 then colnew arbitrarily sets n = 5.
 156 C
 157 C     ipar(4) = no. of solution and derivative tolerances.  ( = ntol )
 158 C               we require  0 .lt. ntol .le. mstar.
 159 C
 160 C     ipar(5) = dimension of fspace.     ( = ndimf )
 161 C
 162 C     ipar(6) = dimension of ispace.     ( = ndimi )
 163 C
 164 C     ipar(7) -  output control ( = iprint )
 165 C              = -1 for full diagnostic printout
 166 C              = 0 for selected printout
 167 C              = 1 for no printout
 168 C
 169 C     ipar(8)     ( = iread )
 170 C             = 0 causes colnew to generate a uniform initial mesh.
 171 C             = 1 if the initial mesh is provided by the user.  it
 172 C                 is defined in fspace as follows:  the mesh
 173 C                 aleft=x(1).lt.x(2).lt. ... .lt.x(n).lt.x(n+1)=aright
 174 C                 will occupy  fspace(1), ..., fspace(n+1). the
 175 C                 user needs to supply only the interior mesh
 176 C                 points  fspace(j) = x(j), j = 2, ..., n.
 177 C             = 2 if the initial mesh is supplied by the user
 178 C                 as with ipar(8)=1, and in addition no adaptive
 179 C                 mesh selection is to be done.
 180 C
 181 C     ipar(9)     ( = iguess )
 182 C             = 0 if no initial guess for the solution is
 183 C                 provided.
 184 C             = 1 if an initial guess is provided by the user
 185 C                 in subroutine  guess.
 186 C             = 2 if an initial mesh and approximate solution
 187 C                 coefficients are provided by the user in  fspace.
 188 C                 (the former and new mesh are the same).
 189 C             = 3 if a former mesh and approximate solution
 190 C                 coefficients are provided by the user in fspace,
 191 C                 and the new mesh is to be taken twice as coarse;
 192 C                 i.e.,every second point from the former mesh.
 193 C             = 4 if in addition to a former initial mesh and
 194 C                 approximate solution coefficients, a new mesh
 195 C                 is provided in fspace as well.
 196 C                 (see description of output for further details
 197 C                 on iguess = 2, 3, and 4.)
 198 C
 199 C     ipar(10)= 0 if the problem is regular
 200 C             = 1 if the first relax factor is =rstart, and the
 201 C                 nonlinear iteration does not rely on past covergence
 202 C                 (use for an extra sensitive nonlinear problem only).
 203 C             = 2 if we are to return immediately upon  (a) two
 204 C                 successive nonconvergences, or  (b) after obtaining
 205 C                 error estimate for the first time.
 206 C
 207 C     ipar(11)= no. of fixed points in the mesh other than aleft
 208 C               and aright. ( = nfxpnt , the dimension of fixpnt)
 209 C               the code requires that all side condition points
 210 C               other than aleft and aright (see description of
 211 C               zeta ) be included as fixed points in fixpnt.
 212 C
 213 C     ltol  -  an array of dimension ipar(4). ltol(j) = l  specifies
 214 C              that the j-th tolerance in  tol  controls the error
 215 C              in the l-th component of z(u).   also require that
 216 C              1.le.ltol(1).lt.ltol(2).lt. ... .lt.ltol(ntol).le.mstar
 217 C
 218 C     tol    - an array of dimension ipar(4). tol(j) is the
 219 C              error tolerance on the ltol(j) -th component
 220 C              of z(u). thus, the code attempts to satisfy
 221 C              for j=1,...,ntol  on each subinterval
 222 C              abs(z(v)-z(u))       .le. tol(j)*abs(z(u))       +tol(j)
 223 C                            ltol(j)                     ltol(j)
 224 C
 225 C              if v(x) is the approximate solution vector.
 226 C
 227 C     fixpnt - an array of dimension ipar(11).   it contains
 228 C              the points, other than aleft and aright, which
 229 C              are to be included in every mesh.
 230 C
 231 C     ispace - an integer work array of dimension ipar(6).
 232 C              its size provides a constraint on nmax,
 233 C              the maximum number of subintervals. choose
 234 C              ipar(6) according to the formula
 235 C                      ipar(6)  .ge.  nmax*nsizei
 236 C                where
 237 C                      nsizei = 3 + kdm
 238 C                with
 239 C                      kdm = kd + mstar  ;  kd = k * ncomp ;
 240 C                      nrec = no. of right end boundary conditions.
 241 C
 242 C
 243 C     fspace - a real work array of dimension ipar(5).
 244 C              its size provides a constraint on nmax.
 245 C              choose ipar(5) according to the formula
 246 C                      ipar(5)  .ge.  nmax*nsizef
 247 C                where
 248 C                      nsizef = 4 + 3 * mstar + (5+kd) * kdm +
 249 C                              (2*mstar-nrec) * 2*mstar.
 250 C
 251 C
 252 C     iflag - the mode of return from colnew.
 253 C           = 1 for normal return
 254 C           = 0 if the collocation matrix is singular.
 255 C           =-1 if the expected no. of subintervals exceeds storage
 256 C               specifications.
 257 C           =-2 if the nonlinear iteration has not converged.
 258 C           =-3 if there is an input data error.
 259 C
 260 C
 261 C**********************************************************************
 262 C
 263 C     *************    user supplied subroutines   *************
 264 C
 265 C
 266 C     the following subroutines must be declared external in the
 267 C     main program which calls colnew.
 268 C
 269 C
 270 C     fsub  - name of subroutine for evaluating f(x,z(u(x))) =
 271 C                            t
 272 C             (f ,...,f     )  at a point x in (aleft,aright).  it
 273 C               1      ncomp
 274 C             should have the heading
 275 C
 276 C                       subroutine fsub (x , z , f)
 277 C
 278 C             where f is the vector containing the value of fi(x,z(u))
 279 C             in the i-th component and                            t
 280 C                                       z(u(x))=(z(1),...,z(mstar))
 281 C             is defined as above under  purpose .
 282 C
 283 C
 284 C     dfsub - name of subroutine for evaluating the jacobian of
 285 C             f(x,z(u)) at a point x.  it should have the heading
 286 C
 287 C                       subroutine dfsub (x , z , df)
 288 C
 289 C             where z(u(x)) is defined as for fsub and the (ncomp) by
 290 C             (mstar) array df should be filled by the partial deriv-
 291 C             atives of f, viz, for a particular call one calculates
 292 C                                df(i,j) = dfi / dzj, i=1,...,ncomp
 293 C                                                     j=1,...,mstar.
 294 C
 295 C
 296 C     gsub  - name of subroutine for evaluating the i-th component of
 297 C             g(x,z(u(x))) = g (zeta(i),z(u(zeta(i)))) at a point x =
 298 C                             i
 299 C             zeta(i) where 1.le.i.le.mstar. it should have the heading
 300 C
 301 C                       subroutine gsub (i , z , g)
 302 C
 303 C             where z(u) is as for fsub, and i and g=g  are as above.
 304 C                                                     i
 305 C             note that in contrast to f in  fsub , here
 306 C             only one value per call is returned in g.
 307 C
 308 C
 309 C     dgsub - name of subroutine for evaluating the i-th row of
 310 C             the jacobian of g(x,u(x)).  it should have the heading
 311 C
 312 C                       subroutine dgsub (i , z , dg)
 313 C
 314 C             where z(u) is as for fsub, i as for gsub and the mstar-
 315 C             vector dg should be filled with the partial derivatives
 316 C             of g, viz, for a particular call one calculates
 317 C                   dg(i,j) = dgi / dzj      j=1,...,mstar.
 318 C
 319 C
 320 C     guess - name of subroutine to evaluate the initial
 321 C             approximation for  z(u(x)) and for dmval(u(x))= vector
 322 C             of the mj-th derivatives of u(x). it should have the
 323 C             heading
 324 C
 325 C                       subroutine guess (x , z , dmval)
 326 C
 327 C             note that this subroutine is needed only if using
 328 C             ipar(9) = 1, and then all  mstar  components of z
 329 C             and  ncomp  components of  dmval  should be specified
 330 C             for any x,  aleft .le. x .le. aright .
 331 C
 332 C
 333 C**********************************************************************
 334 C
 335 C     ************   use of output from colnew   ************
 336 C
 337 C                 ***   solution evaluation   ***
 338 C
 339 C     on return from colnew, the arrays fspace and ispace
 340 C     contain information specifying the approximate solution.
 341 C     the user can produce the solution vector  z( u(x) )  at
 342 C     any point x, aleft .le. x .le. aright, by the statement,
 343 C
 344 C           call appsln (x, z, fspace, ispace)
 345 C
 346 C     when saving the coefficients for later reference, only
 347 C     ispace(1),...,ispace(7+ncomp)    and
 348 C     fspace(1),...,fspace(ispace(7))    need to be saved as
 349 C     these are the quantities used by appsln.
 350 C
 351 C
 352 C                 ***   simple continuation   ***
 353 C
 354 C
 355 C     a formerly obtained solution can easily be used as the
 356 C     first approximation for the nonlinear iteration for a
 357 C     new problem by setting   (iguess =) ipar(9) = 2, 3 or 4.
 358 C
 359 C     if the former solution has just been obtained then the
 360 C     values needed to define the first approximation are
 361 C     already in ispace and fspace.
 362 C     alternatively, if the former solution was obtained in a
 363 C     previous run and its coefficients were saved then those
 364 C     coefficients must be put back into
 365 C     ispace(1),..., ispace(7+ncomp)    and
 366 C     fspace(1),..., fspace(ispace(7)).
 367 C
 368 C     for ipar(9) = 2 or 3 set ipar(3) = ispace(1) ( = the
 369 C     size of the previous mesh ).
 370 C
 371 C     for ipar(9) = 4 the user specifies a new mesh of n subintervals
 372 C     as follows.
 373 C     the values in  fspace(1),...,fspace(ispace(7))  have to be
 374 C     shifted by n+1 locations to  fspace(n+2),..,fspace(ispace(7)+n+1)
 375 C     and the new mesh is then specified in fspace(1),..., fspace(n+1).
 376 C     also set ipar(3) = n.
 377 C
 378 C
 379 C**********************************************************************
 380 C
 381 C     ***************      package subroutines      ***************
 382 C
 383 C     the following description gives a brief overview of how the
 384 C     procedure is broken down into the subroutines which make up
 385 C     the package called  colnew . for further details the
 386 C     user should refer to documentation in the various subroutines
 387 C     and to the references cited above.
 388 C
 389 C     the subroutines fall into four groups:
 390 C
 391 C part 1 - the main storage allocation and program control subr
 392 C
 393 C     colnew - tests input values, does initialization and breaks up
 394 C              the work areas, fspace and ispace, into the arrays
 395 C              used by the program.
 396 C     colsys - another name for colnew
 397 C
 398 C     contrl - is the actual driver of the package. this routine
 399 C              contains the strategy for nonlinear equation solving.
 400 C
 401 C     skale  - provides scaling for the control
 402 C              of convergence in the nonlinear iteration.
 403 C
 404 C
 405 C part 2 - mesh selection and error estimation subroutines
 406 C
 407 C     consts - is called once by  colnew  to initialize constants
 408 C              which are used for error estimation and mesh selection.
 409 C
 410 C     newmsh - generates meshes. it contains the test to decide
 411 C              whether or not to redistribute a mesh.
 412 C
 413 C     errchk - produces error estimates and checks against the
 414 C              tolerances at each subinterval
 415 C
 416 C
 417 C part 3 - collocation system set-up subroutines
 418 C
 419 C     lsyslv - controls the set-up and solution of the linear
 420 C              algebraic systems of collocation equations which
 421 C              arise at each newton iteration.
 422 C
 423 C     gderiv - is used by lsyslv to set up the equation associated
 424 C              with a side condition point.
 425 C
 426 C     vwblok - is used by lsyslv to set up the equation(s) associated
 427 C              with a collocation point.
 428 C
 429 C     gblock - is used by lsyslv to construct a block of the global
 430 C              collocation matrix or the corresponding right hand
 431 C              side.
 432 C
 433 C
 434 C part 4 - service subroutines
 435 C
 436 C     appsln - sets up a standard call to  approx .
 437 C
 438 C     approx - evaluates a piecewise polynomial solution.
 439 C
 440 C     rkbas  - evaluates the mesh independent runge-kutta basis
 441 C
 442 C     vmonde - solves a vandermonde system for given right hand
 443 C              side
 444 C
 445 C     horder - evaluates the highest order derivatives of the
 446 C              current collocation solution used for mesh refinement.
 447 C
 448 C
 449 C part 5 - linear algebra  subroutines
 450 C
 451 C     to solve the global linear systems of collocation equations
 452 C     constructed in part 3,  colnew  uses a column oriented version
 453 C     of the package  solveblok originally due to de boor and weiss.
 454 C
 455 C     to solve the linear systems for static parameter condensation
 456 C     in each block of the collocation equations, the linpack
 457 C     routines  dgefa and  dgesl  are included. but these
 458 C     may be replaced when solving problems on vector processors
 459 C     or when solving large scale sparse jacobian problems.
 460 C
 461 C----------------------------------------------------------------------
 462       IMPLICIT REAL*8 (A-H,O-Z)
 463       DIMENSION M(1), ZETA(1), IPAR(1), LTOL(1), TOL(1), DUMMY(1),
 464      1          FIXPNT(1), ISPACE(1), FSPACE(1)
 465 C
 466       COMMON /COLOUT/ PRECIS, IOUT, IPRINT
 467       COMMON /COLLOC/ RHO(7), COEF(49)
 468       COMMON /COLORD/ K, NC, MSTAR, KD, MMAX, MT(20)
 469       COMMON /COLAPR/ N, NOLD, NMAX, NZ, NDMZ
 470       COMMON /COLMSH/ MSHFLG, MSHNUM, MSHLMT, MSHALT
 471       COMMON /COLSID/ TZETA(40), TLEFT, TRIGHT, IZETA, IDUM
 472       COMMON /COLNLN/ NONLIN, ITER, LIMIT, ICARE, IGUESS
 473       COMMON /COLEST/ TTL(40), WGTMSH(40), WGTERR(40), TOLIN(40),
 474      1                ROOT(40), JTOL(40), LTTOL(40), NTOL
 475 C
 476       EXTERNAL FSUB, DFSUB, GSUB, DGSUB, GUESS
 477 C
 478 C     this subroutine can be called either COLNEW or COLSYS
 479 C
 480       ENTRY      COLSYS (NCOMP, M, ALEFT, ARIGHT, ZETA, IPAR, LTOL,
 481      1                   TOL, FIXPNT, ISPACE, FSPACE, IFLAG,
 482      2                   FSUB, DFSUB, GSUB, DGSUB, GUESS)
 483 C
 484 C*********************************************************************
 485 C
 486 C     the actual subroutine colnew serves as an interface with
 487 C     the package of subroutines referred to collectively as
 488 C     colnew. the subroutine serves to test some of the input
 489 C     parameters, rename some of the parameters (to make under-
 490 C     standing of the coding easier), to do some initialization,
 491 C     and to break the work areas fspace and ispace up into the
 492 C     arrays needed by the program.
 493 C
 494 C**********************************************************************
 495 C
 496 C...  specify machine dependent output unit  iout  and compute machine
 497 C...  dependent constant  precis = 100 * machine unit roundoff
 498 C
 499       IF ( IPAR(7) .LE. 0 )  WRITE(6,99)
 500   99  FORMAT(//,33H VERSION *COLNEW* OF COLSYS .    ,//)
 501 C
 502       IOUT = 6
 503       PRECIS = 1.D0
 504    10 PRECIS = PRECIS / 2.D0
 505       PRECP1 = PRECIS + 1.D0
 506       IF ( PRECP1 .GT. 1.D0 )                       GO TO 10
 507       PRECIS = PRECIS * 100.D0
 508 C
 509 C...  in case incorrect input data is detected, the program returns
 510 C...  immediately with iflag=-3.
 511 C
 512       IFLAG = -3
 513       IF ( NCOMP .LT. 1 .OR. NCOMP .GT. 20 )        RETURN
 514       DO 20 I=1,NCOMP
 515          IF ( M(I) .LT. 1 .OR. M(I) .GT. 4 )        RETURN
 516    20 CONTINUE
 517 C
 518 C...  rename some of the parameters and set default values.
 519 C
 520       NONLIN = IPAR(1)
 521       K = IPAR(2)
 522       N = IPAR(3)
 523       IF ( N .EQ. 0 )  N = 5
 524       IREAD = IPAR(8)
 525       IGUESS = IPAR(9)
 526       IF ( NONLIN .EQ. 0 .AND. IGUESS .EQ. 1 )  IGUESS = 0
 527       IF ( IGUESS .GE. 2 .AND. IREAD .EQ. 0 )   IREAD = 1
 528       ICARE = IPAR(10)
 529       NTOL = IPAR(4)
 530       NDIMF = IPAR(5)
 531       NDIMI = IPAR(6)
 532       NFXPNT = IPAR(11)
 533       IPRINT = IPAR(7)
 534       MSTAR = 0
 535       MMAX = 0
 536       DO  30 I = 1, NCOMP
 537          MMAX = MAX0 ( MMAX, M(I) )
 538          MSTAR = MSTAR + M(I)
 539          MT(I) = M(I)
 540    30 CONTINUE
 541       IF ( K .EQ. 0 )   K = MAX0( MMAX + 1 , 5 - MMAX )
 542       DO 40 I = 1, MSTAR
 543    40 TZETA(I) = ZETA(I)
 544       DO 50 I = 1, NTOL
 545          LTTOL(I) = LTOL(I)
 546    50 TOLIN(I) = TOL(I)
 547       TLEFT = ALEFT
 548       TRIGHT = ARIGHT
 549       NC = NCOMP
 550       KD = K * NCOMP
 551 C
 552 C...  print the input data for checking.
 553 C
 554       IF ( IPRINT .GT. -1 )                         GO TO 80
 555       IF ( NONLIN .GT. 0 )                          GO TO 60
 556       WRITE (IOUT,260) NCOMP, (M(IP), IP=1,NCOMP)
 557       GO TO 70
 558    60 WRITE (IOUT,270) NCOMP, (M(IP), IP=1,NCOMP)
 559    70 WRITE (IOUT,280) (ZETA(IP), IP=1,MSTAR)
 560       IF ( NFXPNT .GT. 0 )
 561      1   WRITE (IOUT,340) NFXPNT, (FIXPNT(IP), IP=1,NFXPNT)
 562       WRITE (IOUT,290) K
 563       WRITE (IOUT,300) (LTOL(IP), IP=1,NTOL)
 564       WRITE (IOUT,310) (TOL(IP), IP=1,NTOL)
 565       IF (IGUESS .GE. 2) WRITE (IOUT,320)
 566       IF (IREAD .EQ. 2) WRITE (IOUT,330)
 567    80 CONTINUE
 568 C
 569 C...  check for correctness of data
 570 C
 571       IF ( K .LT. 0 .OR. K .GT. 7 )                 RETURN
 572       IF ( N .LT. 0 )                               RETURN
 573       IF ( IREAD .LT. 0 .OR. IREAD .GT. 2 )         RETURN
 574       IF ( IGUESS .LT. 0 .OR. IGUESS .GT. 4 )       RETURN
 575       IF ( ICARE .LT. 0 .OR. ICARE .GT. 2 )         RETURN
 576       IF ( NTOL .LT. 0 .OR. NTOL .GT. MSTAR )       RETURN
 577       IF ( NFXPNT .LT. 0 )                          RETURN
 578       IF ( IPRINT .LT. (-1) .OR. IPRINT .GT. 1 )    RETURN
 579       IF ( MSTAR .LT. 0 .OR. MSTAR .GT. 40 )        RETURN
 580       IP = 1
 581       DO 100 I = 1, MSTAR
 582       IF ( DABS(ZETA(I) - ALEFT) .LT. PRECIS .OR.
 583      1     DABS(ZETA(I) - ARIGHT) .LT. PRECIS )     GO TO 100
 584    90 IF ( IP .GT. NFXPNT )                         RETURN
 585         IF ( ZETA(I) - PRECIS .LT. FIXPNT(IP) )     GO TO 95
 586         IP = IP + 1
 587       GO TO 90
 588    95 IF ( ZETA(I) + PRECIS .LT. FIXPNT(IP) )       RETURN
 589   100 CONTINUE
 590 C
 591 C...  set limits on iterations and initialize counters.
 592 C...  limit = maximum number of newton iterations per mesh.
 593 C...  see subroutine  newmsh  for the roles of  mshlmt , mshflg ,
 594 C...  mshnum , and  mshalt .
 595 C
 596       MSHLMT = 3
 597       MSHFLG = 0
 598       MSHNUM = 1
 599       MSHALT = 1
 600       LIMIT = 40
 601 C
 602 C...  compute the maxium possible n for the given sizes of
 603 C...  ispace  and  fspace.
 604 C
 605       NREC = 0
 606       DO 110 I = 1, MSTAR
 607            IB = MSTAR + 1 - I
 608            IF ( ZETA(IB) .GE. ARIGHT )  NREC = I
 609   110 CONTINUE
 610       NFIXI = MSTAR
 611       NSIZEI = 3 + KD + MSTAR
 612       NFIXF = NREC * (2*MSTAR) + 5 * MSTAR + 3
 613       NSIZEF = 4 + 3 * MSTAR + (KD+5) * (KD+MSTAR) +
 614      1(2*MSTAR-NREC) * 2*MSTAR
 615       NMAXF = (NDIMF - NFIXF) / NSIZEF
 616       NMAXI = (NDIMI - NFIXI) / NSIZEI
 617       IF ( IPRINT .LT. 1 )  WRITE(IOUT,350) NMAXF, NMAXI
 618       NMAX = MIN0( NMAXF, NMAXI )
 619       IF ( NMAX .LT. N )                            RETURN
 620       IF ( NMAX .LT. NFXPNT+1 )                     RETURN
 621       IF (NMAX .LT. 2*NFXPNT+2 .AND. IPRINT .LT. 1)  WRITE(IOUT,360)
 622 C
 623 C...  generate pointers to break up  fspace  and  ispace .
 624 C
 625       LXI = 1
 626       LG = LXI + NMAX + 1
 627       LXIOLD = LG + 2*MSTAR * (NMAX * (2*MSTAR-NREC) + NREC)
 628       LW     = LXIOLD + NMAX + 1
 629       LV     = LW + KD**2 * NMAX
 630       LZ     = LV + MSTAR * KD * NMAX
 631       LDMZ   = LZ + MSTAR * (NMAX + 1)
 632       LDELZ  = LDMZ + KD * NMAX
 633       LDELDZ = LDELZ + MSTAR * (NMAX + 1)
 634       LDQZ   = LDELDZ + KD * NMAX
 635       LDQDMZ = LDQZ + MSTAR * (NMAX + 1)
 636       LRHS   = LDQDMZ + KD * NMAX
 637       LVALST = LRHS   + KD * NMAX + MSTAR
 638       LSLOPE = LVALST + 4 * MSTAR * NMAX
 639       LACCUM = LSLOPE + NMAX
 640       LSCL   = LACCUM + NMAX + 1
 641       LDSCL  = LSCL + MSTAR * (NMAX + 1)
 642       LPVTG = 1
 643       LPVTW = LPVTG + MSTAR * (NMAX + 1)
 644       LINTEG = LPVTW + KD * NMAX
 645 C
 646 C...  if  iguess .ge. 2, move  xiold, z, and  dmz  to their proper
 647 C...  locations in  fspace.
 648 C
 649       IF ( IGUESS .LT. 2 )                          GO TO 160
 650       NOLD = N
 651       IF (IGUESS .EQ. 4)  NOLD = ISPACE(1)
 652       NZ = MSTAR * (NOLD + 1)
 653       NDMZ = KD * NOLD
 654       NP1 = N + 1
 655       IF ( IGUESS .EQ. 4 )  NP1 = NP1 + NOLD + 1
 656       DO 120 I=1,NZ
 657   120 FSPACE( LZ+I-1 )  =  FSPACE( NP1+I )
 658       IDMZ = NP1 + NZ
 659       DO 125 I=1,NDMZ
 660   125 FSPACE( LDMZ+I-1 )  =  FSPACE( IDMZ+I )
 661       NP1 = NOLD + 1
 662       IF ( IGUESS .EQ. 4 )                          GO TO 140
 663       DO 130 I=1,NP1
 664   130 FSPACE( LXIOLD+I-1 )  =  FSPACE( LXI+I-1 )
 665       GO TO 160
 666   140 DO 150 I=1,NP1
 667   150 FSPACE( LXIOLD+I-1 )  =  FSPACE( N+1+I )
 668   160 CONTINUE
 669 C
 670 C...  initialize collocation points, constants, mesh.
 671 C
 672       CALL CONSTS ( K, RHO, COEF )
 673       CALL NEWMSH (3+IREAD, FSPACE(LXI), FSPACE(LXIOLD), DUMMY,
 674      1             DUMMY, DUMMY, DUMMY, DUMMY, NFXPNT, FIXPNT)
 675 C
 676 C...  determine first approximation, if the problem is nonlinear.
 677 C
 678       IF (IGUESS .GE. 2)                            GO TO 230
 679       NP1 = N + 1
 680       DO 210 I = 1, NP1
 681   210 FSPACE( I + LXIOLD - 1 ) = FSPACE( I + LXI - 1 )
 682       NOLD = N
 683       IF ( NONLIN .EQ. 0  .OR. IGUESS .EQ. 1 )      GO TO 230
 684 C
 685 C...  system provides first approximation of the solution.
 686 C...  choose z(j) = 0  for j=1,...,mstar.
 687 C
 688       DO 220 I=1, NZ
 689   220 FSPACE( LZ-1+I ) = 0.D0
 690       DO 225 I=1, NDMZ
 691   225 FSPACE( LDMZ-1+I ) = 0.D0
 692   230 CONTINUE
 693       IF (IGUESS .GE. 2)  IGUESS = 0
 694       CALL CONTRL (FSPACE(LXI),FSPACE(LXIOLD),FSPACE(LZ),FSPACE(LDMZ),
 695      1     FSPACE(LRHS),FSPACE(LDELZ),FSPACE(LDELDZ),FSPACE(LDQZ),
 696      2     FSPACE(LDQDMZ),FSPACE(LG),FSPACE(LW),FSPACE(LV),
 697      3     FSPACE(LVALST),FSPACE(LSLOPE),FSPACE(LSCL),FSPACE(LDSCL),
 698      4     FSPACE(LACCUM),ISPACE(LPVTG),ISPACE(LINTEG),ISPACE(LPVTW),
 699      5     NFXPNT,FIXPNT,IFLAG,FSUB,DFSUB,GSUB,DGSUB,GUESS )
 700 C
 701 C...  prepare output
 702 C
 703       ISPACE(1) = N
 704       ISPACE(2) = K
 705       ISPACE(3) = NCOMP
 706       ISPACE(4) = MSTAR
 707       ISPACE(5) = MMAX
 708       ISPACE(6) = NZ + NDMZ + N + 2
 709       K2 = K * K
 710       ISPACE(7) = ISPACE(6) + K2 - 1
 711       DO 240 I = 1, NCOMP
 712   240 ISPACE(7+I) = M(I)
 713       DO 250 I = 1, NZ
 714   250 FSPACE( N+1+I ) = FSPACE( LZ-1+I )
 715       IDMZ = N + 1 + NZ
 716       DO 255 I = 1, NDMZ
 717   255 FSPACE( IDMZ+I ) = FSPACE( LDMZ-1+I )
 718       IC = IDMZ + NDMZ
 719       DO 258 I = 1, K2
 720   258 FSPACE( IC+I ) = COEF(I)
 721       RETURN
 722 C----------------------------------------------------------------------
 723   260 FORMAT(/// 37H THE NUMBER OF (LINEAR) DIFF EQNS IS , I3/ 1X,
 724      1       16HTHEIR ORDERS ARE, 20I3)
 725   270 FORMAT(/// 40H THE NUMBER OF (NONLINEAR) DIFF EQNS IS , I3/ 1X,
 726      1       16HTHEIR ORDERS ARE, 20I3)
 727   280 FORMAT(27H SIDE CONDITION POINTS ZETA, 8F10.6, 4( / 27X, 8F10.6))
 728   290 FORMAT(37H NUMBER OF COLLOC PTS PER INTERVAL IS, I3)
 729   300 FORMAT(39H COMPONENTS OF Z REQUIRING TOLERANCES -,8(7X,I2,1X),
 730      1       4(/38X,8I10))
 731   310 FORMAT(33H CORRESPONDING ERROR TOLERANCES -,6X,8D10.2,
 732      1       4(/39X,8D10.2))
 733   320 FORMAT(44H INITIAL MESH(ES) AND Z,DMZ PROVIDED BY USER)
 734   330 FORMAT(27H NO ADAPTIVE MESH SELECTION)
 735   340 FORMAT(10H THERE ARE ,I5,27H FIXED POINTS IN THE MESH - ,
 736      1       10(6F10.6/))
 737   350 FORMAT(44H THE MAXIMUM NUMBER OF SUBINTERVALS IS MIN (, I4,
 738      1       23H (ALLOWED FROM FSPACE),,I4, 24H (ALLOWED FROM ISPACE) ))
 739   360 FORMAT(/53H INSUFFICIENT SPACE TO DOUBLE MESH FOR ERROR ESTIMATE)
 740       END