Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : ! Ernst Hairer's copyright for rodas can be found at the end of this file.
21 :
22 : module mod_rosenbrock
23 : use mod_dc_decsol
24 : use utils_lib
25 : use const_def, only: dp
26 : use math_lib
27 :
28 : logical, parameter :: dbg = .false.
29 :
30 : integer, parameter :: ns_max = 8 ! current max allowed value for number of stages
31 : ! okay to increase this if necessary.
32 :
33 : contains
34 :
35 0 : subroutine null_mas(n,am,lmas,lrpar,rpar,lipar,ipar)
36 : integer, intent(in) :: n, lmas, lrpar, lipar
37 : real(dp), intent(inout) :: am(lmas,n)
38 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
39 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
40 0 : am = 0
41 0 : end subroutine null_mas
42 :
43 :
44 3 : subroutine do_ros2(
45 : & n,fcn,ifcn,x,y,xend,
46 : & h,max_step_size,max_steps,
47 : & rtol,atol,itol,y_min,y_max,
48 : & jac,ijac,sjac,nzmax,isparse,
49 : & mljac_in,mujac_in,dfx,idfx,
50 : & mas,imas,mlmas,mumas,
51 : & solout,iout,
52 : & decsol, decsols, decsolblk,
53 : & lrd, rpar_decsol, lid, ipar_decsol,
54 : & caller_id, nvar, nz, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk,
55 : & fcn_blk_dble, jac_blk_dble,
56 : & work,lwork,iwork,liwork,
57 : & lrpar,rpar,lipar,ipar,
58 : & lout,idid)
59 : implicit real(dp) (a-h,o-z)
60 : #include "rodas_args.dek"
61 : integer, parameter :: ns = 2 ! number of stages
62 : call do_rodas(
63 : & ns,contro3,ros2_coeffs,n,fcn,ifcn,x,y,xend,
64 : & h,max_step_size,max_steps,
65 : & rtol,atol,itol,y_min,y_max,
66 : & jac,ijac,sjac,nzmax,isparse,
67 : & mljac_in,mujac_in,dfx,idfx,
68 : & mas,imas,mlmas,mumas,
69 : & solout,iout,
70 : & decsol, decsols, decsolblk,
71 : & lrd, rpar_decsol, lid, ipar_decsol,
72 : & caller_id, nvar, nz, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk,
73 : & fcn_blk_dble, jac_blk_dble,
74 : & work,lwork,iwork,liwork,
75 : & lrpar,rpar,lipar,ipar,
76 3 : & lout,idid)
77 3 : end subroutine do_ros2
78 :
79 :
80 3 : subroutine do_rose2(
81 : & n,fcn,ifcn,x,y,xend,
82 : & h,max_step_size,max_steps,
83 : & rtol,atol,itol,y_min,y_max,
84 : & jac,ijac,sjac,nzmax,isparse,
85 : & mljac_in,mujac_in,dfx,idfx,
86 : & mas,imas,mlmas,mumas,
87 : & solout,iout,
88 : & decsol, decsols, decsolblk,
89 : & lrd, rpar_decsol, lid, ipar_decsol,
90 : & caller_id, nvar, nz, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk,
91 : & fcn_blk_dble, jac_blk_dble,
92 : & work,lwork,iwork,liwork,
93 : & lrpar,rpar,lipar,ipar,
94 : & lout,idid)
95 : implicit real(dp) (a-h,o-z)
96 : #include "rodas_args.dek"
97 : integer, parameter :: ns = 3 ! number of stages
98 : call do_rodas(
99 : & ns,contro3,rose2_coeffs,n,fcn,ifcn,x,y,xend,
100 : & h,max_step_size,max_steps,
101 : & rtol,atol,itol,y_min,y_max,
102 : & jac,ijac,sjac,nzmax,isparse,
103 : & mljac_in,mujac_in,dfx,idfx,
104 : & mas,imas,mlmas,mumas,
105 : & solout,iout,
106 : & decsol, decsols, decsolblk,
107 : & lrd, rpar_decsol, lid, ipar_decsol,
108 : & caller_id, nvar, nz, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk,
109 : & fcn_blk_dble, jac_blk_dble,
110 : & work,lwork,iwork,liwork,
111 : & lrpar,rpar,lipar,ipar,
112 3 : & lout,idid)
113 3 : end subroutine do_rose2
114 :
115 :
116 3 : subroutine do_ros3p(
117 : & n,fcn,ifcn,x,y,xend,
118 : & h,max_step_size,max_steps,
119 : & rtol,atol,itol,y_min,y_max,
120 : & jac,ijac,sjac,nzmax,isparse,
121 : & mljac_in,mujac_in,dfx,idfx,
122 : & mas,imas,mlmas,mumas,
123 : & solout,iout,
124 : & decsol, decsols, decsolblk,
125 : & lrd, rpar_decsol, lid, ipar_decsol,
126 : & caller_id, nvar, nz, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk,
127 : & fcn_blk_dble, jac_blk_dble,
128 : & work,lwork,iwork,liwork,
129 : & lrpar,rpar,lipar,ipar,
130 : & lout,idid)
131 : implicit real(dp) (a-h,o-z)
132 : #include "rodas_args.dek"
133 : integer, parameter :: ns = 3 ! number of stages
134 : call do_rodas(
135 : & ns,contro3,ros3p_coeffs,n,fcn,ifcn,x,y,xend,
136 : & h,max_step_size,max_steps,
137 : & rtol,atol,itol,y_min,y_max,
138 : & jac,ijac,sjac,nzmax,isparse,
139 : & mljac_in,mujac_in,dfx,idfx,
140 : & mas,imas,mlmas,mumas,
141 : & solout,iout,
142 : & decsol, decsols, decsolblk,
143 : & lrd, rpar_decsol, lid, ipar_decsol,
144 : & caller_id, nvar, nz, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk,
145 : & fcn_blk_dble, jac_blk_dble,
146 : & work,lwork,iwork,liwork,
147 : & lrpar,rpar,lipar,ipar,
148 3 : & lout,idid)
149 3 : end subroutine do_ros3p
150 :
151 :
152 7 : subroutine do_ros3pl(
153 : & n,fcn,ifcn,x,y,xend,
154 : & h,max_step_size,max_steps,
155 : & rtol,atol,itol,y_min,y_max,
156 : & jac,ijac,sjac,nzmax,isparse,
157 : & mljac_in,mujac_in,dfx,idfx,
158 : & mas,imas,mlmas,mumas,
159 : & solout,iout,
160 : & decsol, decsols, decsolblk,
161 : & lrd, rpar_decsol, lid, ipar_decsol,
162 : & caller_id, nvar, nz, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk,
163 : & fcn_blk_dble, jac_blk_dble,
164 : & work,lwork,iwork,liwork,
165 : & lrpar,rpar,lipar,ipar,
166 : & lout,idid)
167 : implicit real(dp) (a-h,o-z)
168 : #include "rodas_args.dek"
169 : integer, parameter :: ns = 4 ! number of stages
170 : call do_rodas(
171 : & ns,contro3,ros3pl_coeffs,n,fcn,ifcn,x,y,xend,
172 : & h,max_step_size,max_steps,
173 : & rtol,atol,itol,y_min,y_max,
174 : & jac,ijac,sjac,nzmax,isparse,
175 : & mljac_in,mujac_in,dfx,idfx,
176 : & mas,imas,mlmas,mumas,
177 : & solout,iout,
178 : & decsol, decsols, decsolblk,
179 : & lrd, rpar_decsol, lid, ipar_decsol,
180 : & caller_id, nvar, nz, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk,
181 : & fcn_blk_dble, jac_blk_dble,
182 : & work,lwork,iwork,liwork,
183 : & lrpar,rpar,lipar,ipar,
184 7 : & lout,idid)
185 7 : end subroutine do_ros3pl
186 :
187 :
188 7 : subroutine do_rodas3(
189 : & n,fcn,ifcn,x,y,xend,
190 : & h,max_step_size,max_steps,
191 : & rtol,atol,itol,y_min,y_max,
192 : & jac,ijac,sjac,nzmax,isparse,
193 : & mljac_in,mujac_in,dfx,idfx,
194 : & mas,imas,mlmas,mumas,
195 : & solout,iout,
196 : & decsol, decsols, decsolblk,
197 : & lrd, rpar_decsol, lid, ipar_decsol,
198 : & caller_id, nvar, nz, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk,
199 : & fcn_blk_dble, jac_blk_dble,
200 : & work,lwork,iwork,liwork,
201 : & lrpar,rpar,lipar,ipar,
202 : & lout,idid)
203 : implicit real(dp) (a-h,o-z)
204 : #include "rodas_args.dek"
205 : integer, parameter :: ns = 4 ! number of stages
206 : call do_rodas(
207 : & ns,contro3,rodas3_coeffs,n,fcn,ifcn,x,y,xend,
208 : & h,max_step_size,max_steps,
209 : & rtol,atol,itol,y_min,y_max,
210 : & jac,ijac,sjac,nzmax,isparse,
211 : & mljac_in,mujac_in,dfx,idfx,
212 : & mas,imas,mlmas,mumas,
213 : & solout,iout,
214 : & decsol, decsols, decsolblk,
215 : & lrd, rpar_decsol, lid, ipar_decsol,
216 : & caller_id, nvar, nz, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk,
217 : & fcn_blk_dble, jac_blk_dble,
218 : & work,lwork,iwork,liwork,
219 : & lrpar,rpar,lipar,ipar,
220 7 : & lout,idid)
221 7 : end subroutine do_rodas3
222 :
223 :
224 9 : subroutine do_rodas4(
225 : & n,fcn,ifcn,x,y,xend,
226 : & h,max_step_size,max_steps,
227 : & rtol,atol,itol,y_min,y_max,
228 : & jac,ijac,sjac,nzmax,isparse,
229 : & mljac_in,mujac_in,dfx,idfx,
230 : & mas,imas,mlmas,mumas,
231 : & solout,iout,
232 : & decsol, decsols, decsolblk,
233 : & lrd, rpar_decsol, lid, ipar_decsol,
234 : & caller_id, nvar, nz, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk,
235 : & fcn_blk_dble, jac_blk_dble,
236 : & work,lwork,iwork,liwork,
237 : & lrpar,rpar,lipar,ipar,
238 : & lout,idid)
239 : implicit real(dp) (a-h,o-z)
240 : #include "rodas_args.dek"
241 : integer, parameter :: ns = 6 ! number of stages
242 : call do_rodas(
243 : & ns,contro4,rodas4_coeffs,n,fcn,ifcn,x,y,xend,
244 : & h,max_step_size,max_steps,
245 : & rtol,atol,itol,y_min,y_max,
246 : & jac,ijac,sjac,nzmax,isparse,
247 : & mljac_in,mujac_in,dfx,idfx,
248 : & mas,imas,mlmas,mumas,
249 : & solout,iout,
250 : & decsol, decsols, decsolblk,
251 : & lrd, rpar_decsol, lid, ipar_decsol,
252 : & caller_id, nvar, nz, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk,
253 : & fcn_blk_dble, jac_blk_dble,
254 : & work,lwork,iwork,liwork,
255 : & lrpar,rpar,lipar,ipar,
256 9 : & lout,idid)
257 9 : end subroutine do_rodas4
258 :
259 :
260 7 : subroutine do_rodasp(
261 : & n,fcn,ifcn,x,y,xend,
262 : & h,max_step_size,max_steps,
263 : & rtol,atol,itol,y_min,y_max,
264 : & jac,ijac,sjac,nzmax,isparse,
265 : & mljac_in,mujac_in,dfx,idfx,
266 : & mas,imas,mlmas,mumas,
267 : & solout,iout,
268 : & decsol, decsols, decsolblk,
269 : & lrd, rpar_decsol, lid, ipar_decsol,
270 : & caller_id, nvar, nz, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk,
271 : & fcn_blk_dble, jac_blk_dble,
272 : & work,lwork,iwork,liwork,
273 : & lrpar,rpar,lipar,ipar,
274 : & lout,idid)
275 : implicit real(dp) (a-h,o-z)
276 : #include "rodas_args.dek"
277 : integer, parameter :: ns = 6 ! number of stages
278 : call do_rodas(
279 : & ns,contro4,rodasp_coeffs,n,fcn,ifcn,x,y,xend,
280 : & h,max_step_size,max_steps,
281 : & rtol,atol,itol,y_min,y_max,
282 : & jac,ijac,sjac,nzmax,isparse,
283 : & mljac_in,mujac_in,dfx,idfx,
284 : & mas,imas,mlmas,mumas,
285 : & solout,iout,
286 : & decsol, decsols, decsolblk,
287 : & lrd, rpar_decsol, lid, ipar_decsol,
288 : & caller_id, nvar, nz, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk,
289 : & fcn_blk_dble, jac_blk_dble,
290 : & work,lwork,iwork,liwork,
291 : & lrpar,rpar,lipar,ipar,
292 7 : & lout,idid)
293 7 : end subroutine do_rodasp
294 :
295 :
296 39 : subroutine do_rodas(
297 : & ns,contro,coeffs,n,fcn,ifcn,x,y,xend,
298 : & h,max_step_size,max_steps,
299 : & rtol,atol,itol,y_min,y_max,
300 : & jac,ijac,sjac,nzmax,isparse,
301 : & mljac_in,mujac_in,dfx,idfx,
302 : & mas,imas,mlmas,mumas,
303 : & solout,iout,
304 : & decsol, decsols, decsolblk,
305 : & lrd, rpar_decsol, lid, ipar_decsol,
306 : & caller_id, nvar, nz, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk,
307 : & fcn_blk_dble, jac_blk_dble,
308 : & work,lwork,iwork,liwork,
309 : & lrpar,rpar,lipar,ipar,
310 : & lout,idid)
311 : implicit real(dp) (a-h,o-z)
312 : integer, intent(in) :: ns ! number of stages
313 : interface
314 : real(dp) function contro(i,x,rwork,iwork,ierr)
315 : use const_def, only: dp
316 : integer, intent(in) :: i
317 : real(dp), intent(in) :: x
318 : real(dp), intent(inout), target :: rwork(*)
319 : integer, intent(inout), target :: iwork(*)
320 : integer, intent(out) :: ierr
321 : end function contro
322 : subroutine coeffs (ns,ra,rc,rd,ros_d,ros_m,ros_e,ros_alpha,
323 : & ros_gamma,ros_newf,ros_elo,no_aux_in_error,ros_name)
324 : use const_def, only: dp
325 : integer, intent(in) :: ns
326 : real(dp), intent(inout) ::
327 : & ros_m(ns), ros_e(ns), ros_d(ns,ns), ra(ns,ns), rc(ns,ns), rd(ns,ns)
328 : real(dp), intent(inout) :: ros_alpha(ns), ros_gamma(ns)
329 : integer, intent(out) :: ros_elo
330 : logical, intent(out) :: no_aux_in_error, ros_newf(ns)
331 : character (len=12), intent(out) :: ros_name
332 : end subroutine coeffs
333 : end interface
334 : #include "rodas_args.dek"
335 :
336 : logical :: autnms,implct,jband,arret,pred
337 : integer :: mljac, mujac, needed_lwork, needed_liwork
338 : real(dp), intent(in) :: y_min, y_max
339 39 : real(dp), pointer, dimension(:) :: p1, p2, p3, p4, p5
340 39 : integer, pointer, dimension(:) :: ip1
341 : integer :: i, iedy, iedy1, iecon, ieak
342 : integer :: ieia, iefx, iee, iemas, iejac, ieja, ieip, ieynew, iesj, iesa, ierr
343 : integer :: ijob, lde, ldjac, ldmas, ldmas2, m1, m2, meth
344 : integer :: nm1, njac, nfcn, nerror, ndec, naccpt, nsol, nstep, nmax, nrejct
345 39 : mljac = mljac_in; mujac = mujac_in
346 : ! *** *** *** *** *** *** ***
347 : ! setting the parameters
348 : ! *** *** *** *** *** *** ***
349 39 : nfcn=0
350 39 : naccpt=0
351 39 : nrejct=0
352 39 : nstep=0
353 39 : njac=0
354 39 : ndec=0
355 39 : nsol=0
356 39 : arret=.false.
357 : ! -------- nmax , the maximal number of steps -----
358 39 : if(max_steps == 0)then
359 0 : nmax=100000
360 : else
361 39 : nmax=max_steps
362 39 : if(nmax <= 0)then
363 0 : if (lout > 0) write(lout,*)' wrong input max_steps=',max_steps
364 : arret=.true.
365 : end if
366 : end if
367 : ! -------- meth coefficients of the method
368 39 : if(iwork(2) == 0)then
369 39 : meth=1
370 : else
371 0 : meth=iwork(2)
372 0 : if(meth <= 0.or.meth >= 4)then
373 0 : if (lout > 0) write(lout,*)' curious input iwork(2)=',iwork(2)
374 : arret=.true.
375 : end if
376 : end if
377 : ! -------- pred step size control
378 39 : if(iwork(3) <= 1)then
379 39 : pred=.true.
380 : else
381 0 : pred=.false.
382 : end if
383 : ! -------- parameter for second order equations
384 39 : m1=iwork(9)
385 39 : m2=iwork(10)
386 39 : nm1=n-m1
387 39 : if (m1 == 0) m2=n
388 39 : if (m2 == 0) m2=m1
389 39 : if (m1 < 0.or.m2 < 0.or.m1+m2 > n) then
390 0 : if (lout > 0) write(lout,*)' curious input for iwork(9,10)=',m1,m2
391 : arret=.true.
392 : end if
393 39 : nerror=iwork(11) ! number of variables to use for tolerances
394 39 : if (nerror == 0) nerror=n
395 : ! -------- uround smallest number satisfying 1.d0+uround>1.d0
396 39 : if(work(1) == 0.d0)then
397 39 : uround=1.d-16
398 : else
399 0 : uround=work(1)
400 0 : if(uround < 1.d-16.or.uround >= 1.d0)then
401 0 : if (lout > 0) write(lout,*)' coefficients have 16 digits, uround=',work(1)
402 : arret=.true.
403 : end if
404 : end if
405 : ! -------- maximal step size
406 39 : if(max_step_size == 0.d0)then
407 37 : hmax=xend-x
408 : else
409 2 : hmax=max_step_size
410 : end if
411 : ! ------- fac1,fac2 parameters for step size selection
412 39 : if(work(3) == 0.d0)then
413 39 : fac1=5.d0
414 : else
415 0 : fac1=1.d0/work(3)
416 : end if
417 39 : if(work(4) == 0.d0)then
418 39 : fac2=1.d0/3d0 ! 6.0d0
419 : ! originally default was 1/6, but that gave poor results on the diffusion test
420 : else
421 0 : fac2=1.d0/work(4)
422 : end if
423 39 : if (fac1 < 1.0d0.or.fac2 > 1.0d0) then
424 0 : if (lout > 0) write(lout,*)' curious input work(3,4)=',work(3),work(4)
425 : arret=.true.
426 : end if
427 : ! --------- safe safety factor in step size prediction
428 39 : if (work(5) == 0.0d0) then
429 39 : safe=0.9d0
430 : else
431 0 : safe=work(5)
432 0 : if (safe <= 0.001d0.or.safe >= 1.0d0) then
433 0 : if (lout > 0) write(lout,*)' curious input for work(5)=',work(5)
434 : arret=.true.
435 : end if
436 : end if
437 : ! --------- check if tolerances are o.k.
438 39 : if (itol == 0) then
439 39 : if (atol(1) <= 0.d0.or.rtol(1) <= 10.d0*uround) then
440 0 : if (lout > 0) write(lout,*) ' tolerances are too small'
441 : arret=.true.
442 : end if
443 : else
444 0 : do i=1,n
445 0 : if (atol(i) <= 0.d0.or.rtol(i) <= 10.d0*uround) then
446 0 : if (lout > 0) write(lout,*) ' tolerances(',i,') are too small'
447 : arret=.true.
448 : end if
449 : end do
450 : end if
451 : ! *** *** *** *** *** *** *** *** *** *** *** *** ***
452 : ! computation of array entries
453 : ! *** *** *** *** *** *** *** *** *** *** *** *** ***
454 : ! ---- autonomous, implicit, banded or not ?
455 39 : autnms=ifcn == 0
456 39 : implct=imas /= 0
457 39 : jband=mljac < nm1 .and. nzmax == 0
458 39 : if ((nzmax > 0) .and. (jband .or. ijac==0 .or. m1 /= 0)) then
459 0 : if (lout > 0) write(lout,*) 'sparse matrix -- nzmax > 0 -- requires ijac=1, m1=0, mljac=n'
460 : arret=.true.
461 : end if
462 : ! -------- computation of the row-dimensions of the 2-arrays ---
463 : ! -- jacobian and matrix e
464 39 : if(jband)then
465 23 : ldjac=mljac+mujac+1
466 23 : lde=mljac+ldjac
467 : else
468 16 : mljac=nm1
469 16 : mujac=nm1
470 16 : ldjac=nm1
471 16 : lde=nm1
472 : end if
473 : ! -- mass matrix
474 39 : if (implct) then
475 0 : if (mlmas /= nm1) then
476 0 : ldmas=mlmas+mumas+1
477 0 : if (nzmax > 0) then ! sparse jacobian
478 0 : ijob=9
479 0 : else if (jband) then
480 0 : ijob=4
481 : else
482 0 : ijob=3
483 : end if
484 : else
485 0 : ldmas=nm1
486 0 : ijob=5
487 : end if
488 : ! ------ bandwidth of "mas" not larger than bandwidth of "jac"
489 0 : if (mlmas > mljac.or.mumas > mujac) then
490 0 : if (lout > 0) then
491 0 : write(lout,*) 'bandwidth of "mas" must not be larger than bandwidth of "jac"'
492 0 : write(lout,*) 'mlmas', mlmas
493 0 : write(lout,*) 'mljac', mljac
494 0 : write(lout,*) 'mumas', mumas
495 0 : write(lout,*) 'mujac', mujac
496 : end if
497 : arret=.true.
498 : end if
499 : else
500 39 : ldmas=0
501 39 : if (nzmax > 0) then ! sparse jacobian
502 0 : ijob=8
503 39 : else if (jband) then
504 23 : ijob=2
505 : else
506 16 : ijob=1
507 : end if
508 : end if
509 39 : ldmas2=max(1,ldmas)
510 :
511 : call calculate_work_sizes(
512 : & n, ns_max, ldjac, nm1, ldmas, lde, nzmax,
513 : & needed_lwork, needed_liwork, ieynew, iedy1, iedy, ieak, iefx, iecon,
514 39 : & iejac, iemas, iee, iesj, iesa, ieip, ieia, ieja)
515 :
516 39 : if(needed_lwork > lwork)then
517 : ierr = 0
518 39 : call realloc_double(work,needed_lwork,ierr)
519 39 : if (ierr /= 0) then
520 0 : write(lout,*) ' insufficient storage for work, min. lwork=',needed_lwork
521 0 : arret=.true.
522 : end if
523 : end if
524 :
525 39 : if(needed_liwork > liwork)then
526 : ierr = 0
527 0 : call realloc_integer(iwork,needed_liwork,ierr)
528 0 : if (ierr /= 0) then
529 0 : write(lout,*) ' insufficient storage for iwork, min. liwork=',needed_liwork
530 : arret=.true.
531 : end if
532 : end if
533 :
534 : ! ------ when a fail has occurred, we return with idid=-1
535 39 : if (arret) then
536 0 : idid=-1
537 0 : return
538 : end if
539 :
540 : ! -------- call to core integrator ------------
541 39 : p1(1:n*ns) => work(ieak:ieak+n*ns-1) ! arg ak1
542 39 : p2(1:lde*nm1) => work(iee:iee+lde*nm1-1) ! arg e1
543 39 : p3(1:n) => work(ieynew:ieynew+n-1) ! ynew
544 39 : p4(1:n) => work(iedy1:iedy1+n-1) ! dy1
545 39 : p5(1:n) => work(iedy:iedy+n-1) ! dy
546 :
547 39 : ip1(1:nm1) => iwork(ieip:ieip+nm1-1)
548 : call roscor(
549 : & ns,contro,coeffs,n,fcn,x,y,xend,hmax,h,rtol,atol,itol,y_min,y_max,
550 : & jac,ijac,sjac,nzmax,isparse,
551 : & mljac,mujac,dfx,idfx,mas,mlmas,mumas,solout,iout,idid,nmax,
552 : & uround,meth,ijob,fac1,fac2,safe,autnms,implct,jband,pred,ldjac,
553 : & lde,ldmas2,p3,p4,p5,p1,
554 : & work(iefx:lwork),work(iejac:lwork),p2,work(iemas:lwork),
555 : & ip1,work(iecon:lwork),
556 : & decsol,decsols,decsolblk,
557 : & caller_id, nvar, nz, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk,
558 : & fcn_blk_dble, jac_blk_dble,
559 : & iwork(ieia:liwork),iwork(ieja:liwork),
560 : & work(iesj:lwork),work(iesa:lwork),lrd,rpar_decsol,lid,ipar_decsol,m1,m2,nm1,nerror,
561 39 : & nfcn,njac,nstep,naccpt,nrejct,ndec,nsol,lout,lrpar,rpar,lipar,ipar)
562 39 : iwork(14)=nfcn
563 39 : iwork(15)=njac
564 39 : iwork(16)=nstep
565 39 : iwork(17)=naccpt
566 39 : iwork(18)=nrejct
567 39 : iwork(19)=ndec
568 39 : iwork(20)=nsol
569 : ! ----------- return -----------
570 39 : return
571 39 : end subroutine do_rodas
572 : !
573 :
574 :
575 39 : subroutine calculate_work_sizes(
576 : & n, ns_max, ldjac, nm1, ldmas, lde, nzmax,
577 : & lwork, liwork, ieynew, iedy1, iedy, ieak, iefx, iecon,
578 : & iejac, iemas, iee, iesj, iesa, ieip, ieia, ieja)
579 : integer, intent(in) :: n, ns_max, ldjac, nm1, ldmas, lde, nzmax
580 : integer, intent(out) ::
581 : & lwork, liwork, ieynew, iedy1, iedy, ieak, iefx, iecon,
582 : & iejac, iemas, iee, iesj, iesa, ieip, ieia, ieja
583 39 : ieynew=21
584 39 : iedy1=ieynew+n
585 39 : iedy=iedy1+n
586 39 : ieak=iedy+n
587 39 : iefx=ieak+n*ns_max
588 39 : iecon=iefx+n
589 39 : iejac=iecon+2+5*n
590 39 : iemas=iejac+n*ldjac
591 39 : iee=iemas+nm1*ldmas
592 39 : iesj=iee+nm1*lde
593 39 : iesa=iesj+nzmax
594 39 : lwork=iesa+nzmax-1
595 39 : ieip=21
596 39 : ieia=ieip+nm1
597 39 : ieja=ieia+n+1
598 39 : liwork=ieja+nzmax-1
599 39 : end subroutine calculate_work_sizes
600 :
601 :
602 : !
603 : !
604 : ! ----- ... and here is the core integrator ----------
605 : !
606 39 : subroutine roscor(
607 : & ns,contro,coeffs,n,fcn,x,y,xend,hmax,h,rtol,atol,itol,y_min,y_max,
608 : & jac,ijac,sjac,nzmax,isparse,mljac,mujac,dfx,idfx,mas,mlmas,mumas,
609 : & solout,iout,idid,nmax,uround,meth,ijob,fac1,fac2,safe,autnms,implct,banded,
610 39 : & pred,ldjac,lde,ldmas,ynew,dy1,dy,ak1,fx,fjac,e1,fmas,ip,rwork,
611 : & decsol,decsols,decsolblk,
612 : & caller_id, nvar, nz, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk,
613 : & fcn_blk_dble, jac_blk_dble,
614 39 : & ia,ja,sparse_jac,sa,
615 : & lrd,rpar_decsol,lid,ipar_decsol,m1,m2,nm1,nerror,
616 : & nfcn,njac,nstep,naccpt,nrejct,ndec,nsol,lout,lrpar,rpar,lipar,ipar)
617 : ! ----------------------------------------------------------
618 : ! core integrator for rodas4
619 : ! parameters same as in rodas4 with workspace added
620 : ! ----------------------------------------------------------
621 : ! declarations
622 : ! ----------------------------------------------------------
623 : implicit real(dp) (a-h,o-z)
624 : integer :: n, itol, ijac, isparse, mljac, mujac, idfx, mlmas, mumas, iout, idid, nmax, meth, ijob,
625 : & ldjac, lde, ldmas, m1, m2, nm1, nerror, nfcn, njac, nstep, naccpt, nrejct, ndec, nsol, lout
626 : interface
627 : real(dp) function contro(i,x,rwork,iwork,ierr)
628 : use const_def, only: dp
629 : integer, intent(in) :: i
630 : real(dp), intent(in) :: x
631 : real(dp), intent(inout), target :: rwork(*)
632 : integer, intent(inout), target :: iwork(*)
633 : integer, intent(out) :: ierr
634 : end function contro
635 : subroutine coeffs(ns,ra,rc,rd,ros_d,ros_m,ros_e,ros_alpha,ros_gamma,ros_newf,ros_elo,no_aux_in_error,ros_name)
636 : use const_def, only: dp
637 : implicit none
638 : integer, intent(in) :: ns
639 : real(dp), intent(inout) :: ros_m(ns), ros_e(ns), ros_d(ns,ns), ra(ns,ns), rc(ns,ns), rd(ns,ns)
640 : real(dp), intent(inout) :: ros_alpha(ns), ros_gamma(ns)
641 : integer, intent(out) :: ros_elo
642 : logical, intent(out) :: no_aux_in_error, ros_newf(ns)
643 : character (len=12), intent(out) :: ros_name
644 : end subroutine coeffs
645 : #include "num_solout.dek"
646 : #include "num_mas.dek"
647 : #include "num_dfx.dek"
648 : #include "num_fcn.dek"
649 : #include "num_fcn_blk_dble.dek"
650 : #include "num_jac.dek"
651 : ! for double block tridiagonal matrix
652 : #include "num_jac_blk_dble.dek"
653 : #include "num_sjac.dek"
654 : #include "mtx_decsol.dek"
655 : #include "mtx_decsols.dek"
656 : #include "mtx_decsolblk.dek"
657 : end interface
658 : integer, intent(in) :: nzmax, lrpar, lipar, lrd, lid
659 : integer, intent(out) :: ia(n+1), ja(nzmax)
660 : real(dp), intent(inout) :: sparse_jac(nzmax), sa(nzmax)
661 : integer, intent(in) :: caller_id, nvar, nz
662 : real(dp), dimension(:), pointer, intent(inout) :: lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk ! =(nvar,nvar,nz)
663 : ! the rosenbrock method parameters
664 : integer, intent(in) :: ns ! number of stages
665 430 : real(dp) :: ros_m(ns), ros_e(ns)
666 430 : real(dp) :: ros_alpha(ns), ros_gamma(ns)
667 : integer :: ros_elo
668 78 : logical :: ros_newf(ns)
669 : character (len=12) :: ros_name
670 :
671 : ! args
672 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
673 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
674 : real(dp), intent(inout), pointer :: y(:) ! (n)
675 : real(dp), intent(inout), pointer :: dy(:) ! (n)
676 : real(dp), intent(inout), pointer :: ynew(:) ! (n)
677 : real(dp), intent(inout), pointer :: dy1(:) ! (n)
678 :
679 : dimension fx(n),fjac(ldjac,n),fmas(ldmas,nm1),atol(*),rtol(*)
680 :
681 : integer :: iwork(2)
682 : real(dp), target :: rwork(2+5*n)
683 : real(dp), intent(inout), pointer :: rpar_decsol(:) ! (lrd)
684 : integer, intent(inout), pointer :: ipar_decsol(:) ! (lid)
685 : real(dp), intent(in) :: y_min, y_max
686 : real(dp), pointer :: ak1(:) ! =(n,ns)
687 : real(dp), pointer :: e1(:) ! =(lde,nm1)
688 : integer, pointer :: ip(:) ! (nm1)
689 :
690 : ! locals
691 5542 : dimension hd(ns), ra(ns,ns), rhc(ns,ns), rc(ns,ns), rd(ns,ns), ros_d(ns,ns)
692 : logical :: reject,autnms,implct,banded,last,pred,not_stage1,no_aux_in_error,need_free
693 39 : real(dp), pointer :: cont(:), dy2(:)
694 215 : real(dp) :: hprev, rd32
695 : integer :: i, j, k, l, j1, is
696 : integer :: lrc, lbeg, lend, irtrn
697 : integer :: mbb, mbdiag, mbjac, md, mdiag, mdiff, mle, mm, mue, mujacj, mujacp
698 : integer :: nsing, nn, nn2, nn3, nn4
699 : integer :: ierr, ier
700 :
701 39 : real(dp), pointer :: ak(:,:), e(:,:), p1(:)
702 39 : ak(1:n,1:ns) => ak1(1:n*ns)
703 39 : e(1:lde,1:nm1) => e1(1:lde*nm1)
704 :
705 : ! *** *** *** *** *** *** ***
706 : ! initialisations
707 : ! *** *** *** *** *** *** ***
708 1081 : rc = 0
709 39 : md = 0 !n/2-1 ! for dbg
710 39 : need_free = .false.
711 39 : cont => rwork(1:5*n)
712 39 : dy2 => cont(n+1:2*n)
713 39 : nn=n
714 39 : nn2=2*n
715 39 : nn3=3*n
716 39 : nn4=4*n
717 39 : lrc=5*n
718 : ! ------- compute mass matrix for implicit case ----------
719 39 : if (implct) call mas (nm1,fmas,ldmas,lrpar,rpar,lipar,ipar)
720 : ! ------ set the parameters of the method -----
721 : call coeffs(ns,ra,rc,rd,ros_d,ros_m,ros_e,ros_alpha,
722 39 : & ros_gamma,ros_newf,ros_elo,no_aux_in_error,ros_name)
723 78 : gamma = ros_gamma(1)
724 39 : if (ns >=3) then
725 36 : rd32 = rd(3,2)
726 : else
727 : rd32 = 0
728 : end if
729 39 : if (no_aux_in_error .and. nerror > m2) nerror = m2
730 : ! --- initial preparations
731 39 : if (m1 > 0) ijob=ijob+10
732 :
733 : if (dbg) write(*,*) trim(ros_name) // ' ijob', ijob
734 :
735 39 : posneg=sign(1.d0,xend-x)
736 78 : hmaxn=min(abs(hmax),abs(xend-x))
737 39 : if (abs(h) <= 10.d0*uround) h=1.0d-6
738 39 : h=min(abs(h),hmaxn)
739 39 : h=sign(h,posneg)
740 39 : reject=.false.
741 39 : last=.false.
742 39 : nsing=0
743 39 : ierr=0
744 39 : irtrn=1
745 215 : hd=0
746 : ! -------- prepare band-widths --------
747 39 : mbdiag=mumas+1
748 39 : if (banded) then
749 23 : mle=mljac
750 23 : mue=mujac
751 23 : mbjac=mljac+mujac+1
752 23 : mbb=mlmas+mumas+1
753 23 : mdiag=mle+mue+1
754 23 : mdiff=mle+mue-mumas
755 : end if
756 39 : if (iout /= 0) then
757 64 : xold=x
758 : irtrn=1
759 64 : hout=h
760 25 : iwork(1) = lrc
761 25 : iwork(2) = n
762 25 : j=1+lrc
763 25 : rwork(j) = xold; j=j+1
764 25 : rwork(j) = h
765 25 : call solout(naccpt,xold,x,n,y,rwork,iwork,contro,lrpar,rpar,lipar,ipar,irtrn)
766 25 : if (irtrn < 0) GOTO 179
767 : end if
768 :
769 : ! --- basic integration step
770 25211 : 1 if (nstep > nmax) GOTO 178
771 25211 : hprev = h
772 25211 : if (0.1d0*abs(h) <= abs(x)*uround) GOTO 177
773 25211 : if (last) then
774 78 : h=hopt
775 39 : idid=1
776 39 : GOTO 191
777 : end if
778 25172 : hopt=h
779 25172 : if ((x+h*1.0001d0-xend)*posneg >= 0.d0) then
780 39 : h=xend-x
781 39 : last=.true.
782 : end if
783 : ! *** *** *** *** *** *** ***
784 : ! computation of the jacobian
785 : ! *** *** *** *** *** *** ***
786 25172 : njac=njac+1
787 25172 : if (ijac == 0) then
788 10652 : if (nvar > 0) then
789 0 : call fcn_blk_dble(n,caller_id,nvar,nz,x,h,y,dy1,lrpar,rpar,lipar,ipar,ierr)
790 : else
791 10652 : call fcn(n,x,h,y,dy1,lrpar,rpar,lipar,ipar,ierr)
792 : end if
793 10652 : if (ierr /= 0) GOTO 180
794 10652 : nfcn=nfcn+1
795 : ! --- compute jacobian matrix numerically
796 10652 : if (banded) then
797 : ! --- jacobian is banded
798 2430 : mujacp=mujac+1
799 2430 : md=min(mbjac,n)
800 4860 : do mm=1,m1/m2+1
801 17010 : do k=1,md
802 12150 : j=k+(mm-1)*m2
803 972000 : 12 ak(j,2)=y(j)
804 972000 : ak(j,3)=dsqrt(uround*max(1.d-5,abs(y(j))))
805 972000 : y(j)=y(j)+ak(j,3)
806 972000 : j=j+md
807 972000 : if (j <= mm*m2) GOTO 12
808 :
809 12150 : p1(1:n) => ak(1:n,1)
810 12150 : if (nvar > 0) then
811 0 : call fcn_blk_dble(n,caller_id,nvar,nz,x,h,y,p1,lrpar,rpar,lipar,ipar,ierr)
812 : else
813 12150 : call fcn(n,x,h,y,p1,lrpar,rpar,lipar,ipar,ierr)
814 : end if
815 :
816 12150 : if (ierr /= 0) GOTO 180
817 12150 : j=k+(mm-1)*m2
818 12150 : j1=k
819 972000 : lbeg=max(1,j1-mujac)+m1
820 972000 : 14 lend=min(m2,j1+mljac)+m1
821 972000 : y(j)=ak(j,2)
822 972000 : mujacj=mujacp-j1-m1
823 5817420 : do l=lbeg,lend
824 5817420 : fjac(l+mujacj,j)=(ak(l,1)-dy1(l))/ak(j,3)
825 : end do
826 972000 : j=j+md
827 972000 : j1=j1+md
828 972000 : lbeg=lend+1
829 974430 : if (j <= mm*m2) GOTO 14
830 : end do
831 : end do
832 : else
833 : ! --- jacobian is full
834 24666 : do i=1,n
835 16483 : ysafe=y(i)
836 16483 : delt=dsqrt(uround*max(1.d-5,abs(ysafe)))
837 16444 : y(i)=ysafe+delt
838 16444 : p1(1:n) => ak(1:n,1)
839 16444 : if (nvar > 0) then
840 0 : call fcn_blk_dble(n,caller_id,nvar,nz,x,h,y,p1,lrpar,rpar,lipar,ipar,ierr)
841 : else
842 16444 : call fcn(n,x,h,y,p1,lrpar,rpar,lipar,ipar,ierr)
843 : end if
844 :
845 16444 : if (ierr /= 0) GOTO 180
846 49332 : do j=m1+1,n
847 49332 : fjac(j-m1,i)=(ak(j,1)-dy1(j))/delt
848 : end do
849 24666 : y(i)=ysafe
850 : end do
851 : end if
852 : else
853 : ! --- compute jacobian matrix analytically
854 14520 : if (nzmax == 0) then
855 : if (dbg) write(*,11) 'jac y(:)', y(1:min(4,n))
856 :
857 14520 : if (nvar > 0) then
858 0 : call jac_blk_dble(n,caller_id,nvar,nz,x,h,y,dy1,uf_lblk,uf_dblk,uf_ublk,lrpar,rpar,lipar,ipar,ierr)
859 : else
860 14520 : call jac(n,x,h,y,dy1,fjac,ldjac,lrpar,rpar,lipar,ipar,ierr)
861 : end if
862 :
863 : if (dbg) write(*,11) 'jac dy1(:)', dy1(1:min(4,n))
864 : else
865 0 : call sjac(n,x,h,y,dy1,nzmax,ia,ja,sparse_jac,lrpar,rpar,lipar,ipar,ierr)
866 : end if
867 14520 : if (ierr /= 0) GOTO 180
868 : end if
869 25172 : if (.not.autnms) then
870 0 : if (idfx == 0) then
871 : ! --- compute numerically the derivative with respect to x
872 0 : delt=sqrt(uround*max(1.d-5,abs(x)))
873 39 : xdelt=x+delt
874 0 : p1(1:n) => ak(1:n,1)
875 0 : if (nvar > 0) then
876 0 : call fcn_blk_dble(n,caller_id,nvar,nz,xdelt,h,y,p1,lrpar,rpar,lipar,ipar,ierr)
877 : else
878 0 : call fcn(n,xdelt,h,y,p1,lrpar,rpar,lipar,ipar,ierr)
879 : end if
880 :
881 0 : if (ierr /= 0) GOTO 180
882 0 : do j=1,n
883 0 : fx(j)=(ak(j,1)-dy1(j))/delt
884 : end do
885 : else
886 : ! --- compute analytically the derivative with respect to x
887 0 : call dfx(n,x,y,fx,lrpar,rpar,lipar,ipar,ierr)
888 0 : if (ierr /= 0) GOTO 180
889 : end if
890 : end if
891 : 2 continue
892 25299 : if (h <= hprev*1d-30) GOTO 190
893 25299 : ierr=0
894 : ! *** *** *** *** *** *** ***
895 : ! compute the stages
896 : ! *** *** *** *** *** *** ***
897 25338 : fac=1.d0/(h*gamma)
898 25299 : if (need_free) then
899 : call decsol_done(n,fjac,ldjac,fmas,ldmas,mlmas,mumas,
900 : & m1,m2,nm1,fac,e1,lde,ip,ak1,ier,ijob,implct,ip,
901 : & mle,mue,mbjac,mbb,mdiag,mdiff,mbdiag,
902 : & decsol,decsols,decsolblk,
903 : & caller_id, nvar, nz, lblk, dblk, ublk,
904 25260 : & sparse_jac,nzmax,isparse,ia,ja,sa,lrd,rpar_decsol,lid,ipar_decsol)
905 : end if
906 : call decomr(n,fjac,ldjac,fmas,ldmas,mlmas,mumas,
907 : & m1,m2,nm1,fac,e1,lde,ip,ak1,ier,ijob,implct,ip,
908 : & mle,mue,mbjac,mbb,mdiag,mdiff,mbdiag,
909 : & decsol,decsols,decsolblk,
910 : & caller_id, nvar, nz, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk,
911 25299 : & sparse_jac,nzmax,isparse,ia,ja,sa,lrd,rpar_decsol,lid,ipar_decsol)
912 25299 : need_free = .true.
913 25299 : if (ier /= 0) GOTO 80
914 : if (dbg .and. .false.) then
915 : write(*,11) 'fac', fac
916 : write(*,*) 'n', n
917 : write(*,*) 'ldjac', ldjac
918 : write(*,*) 'ldmas', ldmas
919 : write(*,*) 'mlmas', mlmas
920 : write(*,*) 'mumas', mumas
921 : write(*,*) 'm1', m1
922 : write(*,*) 'm2', m2
923 : write(*,*) 'nm1', nm1
924 : write(*,*) 'lde', lde
925 : do i=1,lde
926 : write(*,11) 'e', e(i,1), e(i,2), e(i,min(3,n)), e(i,min(4,n))
927 : end do
928 : write(*,*)
929 : end if
930 25299 : ndec=ndec+1
931 : ! --- prepare for the computation of the stages
932 25299 : if (.not.autnms) hd = h*ros_gamma
933 25299 : if (h <= 0d0 .or. is_bad(h)) GOTO 177
934 407921 : rhc = rc/h
935 :
936 : ! ------------ stages ---------------
937 : if (dbg) then
938 : write(*,*)
939 : write(*,*) 'nstep', nstep
940 : write(*,11) 'x', x
941 : write(*,11) 'h', h
942 : if (nstep > 30) stop 'rosenbrock nstep > 1'
943 : end if
944 :
945 25299 : if (is_bad(h)) GOTO 177
946 :
947 106471 : do is=1,ns
948 : if (dbg) write(*,*) 'stage', is
949 81174 : if (is == 1) then
950 4307539 : dy = dy1
951 25299 : not_stage1 = .false.
952 : else
953 55875 : if (ros_newf(is)) then
954 5912391 : do j=1,n
955 20307900 : ynew(j) = y(j) + sum(ra(is,1:is-1)*ak(j,1:is-1))
956 5912391 : if (ynew(j) < y_min .or. ynew(j) > y_max) then
957 : if (dbg) write(*,*) 'stage ynew(j) < y_min .or. ynew(j) > y_max',
958 : & is, j, ynew(j), y(j), sum(ra(is,1:is-1)*ak(j,1:is-1))
959 : GOTO 82
960 : end if
961 : end do
962 : if (dbg) write(*,11) 'fcn y(:)', ynew(1:min(4,n))
963 :
964 49715 : if (nvar > 0) then
965 0 : call fcn_blk_dble(n,caller_id,nvar,nz,x+ros_alpha(is)*h,h,ynew,dy,lrpar,rpar,lipar,ipar,ierr)
966 : else
967 49715 : call fcn(n,x+ros_alpha(is)*h,h,ynew,dy,lrpar,rpar,lipar,ipar,ierr)
968 : end if
969 :
970 : if (dbg) write(*,11) 'fcn dy(:)', dy(1:min(4,n))
971 49715 : if (ierr /= 0) GOTO 81
972 311619 : if (rd32 /= 0) dy2(1:n) = dy(1:n)
973 : end if
974 7346681 : do j=1,n
975 24558559 : cont(j) = sum(rhc(is,1:is-1)*ak(j,1:is-1))
976 : end do
977 55873 : if (is == 3 .and. rd32 /= 0) then
978 72362 : do j=1,n
979 72362 : cont(j) = cont(j) + rd32*dy2(j) + rd(3,1)*dy1(j)
980 : end do
981 : end if
982 55873 : not_stage1 = .true.
983 : end if
984 81172 : p1(1:n) => ak1(1+(is-1)*n:n*is)
985 : call slvrod(n,fjac,ldjac,mljac,mujac,fmas,ldmas,mlmas,mumas,
986 : & m1,m2,nm1,fac,e1,lde,ip,dy,p1,fx,cont,hd(is),ijob,not_stage1,
987 : & mle,mue,mbjac,mbb,mdiag,mdiff,mbdiag,
988 : & decsol,decsols,decsolblk,
989 : & caller_id, nvar, nz, lblk, dblk, ublk, uf_lblk, uf_dblk, uf_ublk,
990 81172 : & nzmax,isparse,ia,ja,sa,lrd,rpar_decsol,lid,ipar_decsol,ierr)
991 81172 : if (ierr /= 0) GOTO 81
992 : if (dbg) write(*,11) 'ak(:,is)', ak(1,is), ak(2,is), ak(min(3,n),is), ak(min(4,n),is)
993 25297 : if (dbg) write(*,*)
994 : end do
995 :
996 : ! ------------ new solution ---------------
997 2166373 : do j=1,n
998 11572828 : ynew(j) = y(j) + sum(ros_m(1:ns)*ak(j,1:ns))
999 2166373 : if (ynew(j) < y_min .or. ynew(j) > y_max) then
1000 : if (dbg) write(*,*) 'new solution ynew(j) < y_min .or. ynew(j) > y_max', j, ynew(j)
1001 : GOTO 82
1002 : end if
1003 : end do
1004 25297 : nsol=nsol+ns
1005 25297 : nfcn=nfcn+ns-1
1006 :
1007 : ! ------------ error estimation ----------
1008 25297 : nstep=nstep+1
1009 : ! put the error vector in dy
1010 2166373 : do i=1,n
1011 11598125 : dy(i) = sum(ros_e(1:ns)*ak(i,1:ns))
1012 : end do
1013 :
1014 25336 : err=0.d0
1015 2166373 : do i=1,nerror
1016 2141076 : if (itol == 0) then
1017 2141076 : sk=atol(1)+rtol(1)*max(abs(y(i)),abs(ynew(i)))
1018 : else
1019 0 : sk=atol(i)+rtol(i)*max(abs(y(i)),abs(ynew(i)))
1020 : !if (dbg) write(*,*) 'sk', i, sk
1021 : end if
1022 2166373 : err=err+(dy(i)/sk)**2
1023 : !if (dbg) write(*,*) 'err', i, err
1024 : end do
1025 :
1026 25297 : err=sqrt(err/nerror)
1027 25297 : if (is_bad(err)) GOTO 81
1028 :
1029 :
1030 : 11 format(a20,2x,99(e26.16,1x))
1031 :
1032 : if (dbg) then
1033 :
1034 : write(*,11) 'y(:)', y(1:min(4,n))
1035 : write(*,11) 'ynew(:)', ynew(1:min(4,n))
1036 : write(*,11) 'err(:)', dy(1:min(4,n))
1037 : write(*,11) 'err=', err
1038 : write(*,*)
1039 : write(*,*)
1040 :
1041 : if (nstep > 20) stop
1042 :
1043 : if (is_bad(y(1)) .or. is_bad(y(2))) stop
1044 :
1045 : end if
1046 :
1047 :
1048 : ! --- computation of hnew
1049 : ! --- we require .2<=hnew/h<=6.
1050 25336 : eloi = 1d0/ros_elo ! inverse of estimated local order
1051 25297 : fac=max(fac2,min(fac1,pow(err,eloi)/safe))
1052 :
1053 2191670 : if (minval(ynew(1:n)) < y_min) then
1054 : if (dbg) write(*,*) 'reject < y_min', minval(ynew(1:n)), y_min
1055 0 : fac = 100
1056 0 : err = 100
1057 2191670 : else if (maxval(ynew(1:n)) > y_max) then
1058 : if (dbg) write(*,*) 'reject > y_max', maxval(ynew(1:n)), y_max
1059 0 : fac = 100
1060 0 : err = 100
1061 : end if
1062 :
1063 :
1064 39 : hnew=h/fac
1065 25297 : if (is_bad(hnew)) GOTO 81
1066 :
1067 25297 : hopt=hnew
1068 :
1069 : ! *** *** *** *** *** *** ***
1070 : ! is the error small enough ?
1071 : ! *** *** *** *** *** *** ***
1072 25297 : if (err <= 1.d0) then
1073 : ! --- step is accepted
1074 25172 : naccpt=naccpt+1
1075 25172 : if (pred) then
1076 : ! --- predictive controller of gustafsson
1077 25172 : if (naccpt > 1) then
1078 39 : facgus=(hacc/h)*pow(err**2/erracc,eloi)/safe
1079 25133 : facgus=max(fac2,min(fac1,facgus))
1080 25133 : fac=max(fac,facgus)
1081 25133 : hnew=h/fac
1082 : end if
1083 25172 : hacc=h
1084 25172 : erracc=max(1.0d-2,err)
1085 : end if
1086 25172 : if (iout /= 0) then
1087 2116549 : do i=1,n
1088 2116549 : cont(i)=y(i)
1089 : end do
1090 : end if
1091 2165878 : do i=1,n
1092 2165878 : y(i)=ynew(i)
1093 : end do
1094 25172 : xold=x
1095 25172 : x=x+h
1096 25172 : if (iout /= 0) then
1097 2116549 : do i=1,n
1098 2116549 : cont(i+nn)=ynew(i)
1099 : end do
1100 8729 : if (ros_elo >= 4) then
1101 558282 : do i=1,n
1102 3335976 : cont(i+nn2)=sum(ros_d(2,1:ns-1)*ak(i,1:ns-1))
1103 3338262 : cont(i+nn3)=sum(ros_d(3,1:ns-1)*ak(i,1:ns-1))
1104 : end do
1105 2286 : if (ros_elo == 5) then
1106 0 : do i=1,n
1107 0 : cont(i+nn4)=sum(ros_d(4,1:ns-1)*ak(i,1:ns-1))
1108 : end do
1109 : end if
1110 : else
1111 1558267 : do i=1,n
1112 1558267 : cont(i+nn2) = dy1(i)
1113 : end do
1114 : end if
1115 : irtrn=1
1116 8729 : hout=h
1117 8729 : iwork(1) = lrc
1118 8729 : iwork(2) = n
1119 8729 : j=1+lrc
1120 8729 : rwork(j) = xold; j=j+1
1121 8729 : rwork(j) = h
1122 8729 : call solout(naccpt,xold,x,n,y,rwork,iwork,contro,lrpar,rpar,lipar,ipar,irtrn)
1123 8729 : if (irtrn < 0) GOTO 179
1124 : end if
1125 25172 : if (abs(hnew) > hmaxn) hnew=posneg*hmaxn
1126 25172 : if (reject) hnew=posneg*min(abs(hnew),abs(h))
1127 25172 : reject=.false.
1128 25172 : h=hnew
1129 25172 : GOTO 1
1130 : else
1131 : ! --- step is rejected
1132 125 : reject=.true.
1133 125 : last=.false.
1134 125 : h=hnew
1135 125 : if (naccpt >= 1) nrejct=nrejct+1
1136 0 : GOTO 2
1137 : end if
1138 : ! --- singular matrix
1139 0 : 80 nsing=nsing+1
1140 0 : if (nsing >= 5) GOTO 176
1141 0 : h=h*0.5d0
1142 0 : reject=.true.
1143 0 : last=.false.
1144 0 : GOTO 2
1145 : ! --- step rejected
1146 0 : 81 h=h*0.5d0
1147 0 : reject=.true.
1148 0 : last=.false.
1149 0 : GOTO 2
1150 : ! --- step rejected for y_min or y_max
1151 2 : 82 h=h*0.1d0
1152 2 : reject=.true.
1153 2 : last=.false.
1154 2 : GOTO 2
1155 : ! --- fail exit
1156 : 176 continue
1157 0 : if (lout > 0) write(lout,979)x
1158 0 : if (lout > 0) write(lout,*) ' matrix is repeatedly singular, ier=',ier
1159 0 : idid=-4
1160 0 : GOTO 191
1161 : 177 continue
1162 0 : if (lout > 0) write(lout,979)x
1163 0 : if (lout > 0) write(lout,*) ' step size too small, h=',h
1164 0 : idid=-3
1165 0 : GOTO 191
1166 : 178 continue
1167 0 : if (lout > 0) write(lout,979)x
1168 0 : if (lout > 0) write(lout,*) ' more than nmax =',nmax,'steps are needed'
1169 0 : idid=-2
1170 0 : GOTO 191
1171 : ! --- solout exit
1172 : 179 continue
1173 0 : if (lout > 0) write(lout,979)x
1174 : 979 format(' exit at x=',e18.4)
1175 0 : idid=2
1176 0 : GOTO 191
1177 : ! --- forced exit because of ierr /= 0
1178 : 180 continue
1179 0 : idid=-5
1180 0 : GOTO 191
1181 : ! --- too many reductions in stepsize and still not okay error
1182 : 190 continue
1183 0 : idid=-8
1184 :
1185 :
1186 : 191 continue
1187 :
1188 39 : if (need_free) then
1189 39 : p1(1:n) => ak(1:n,1)
1190 : call decsol_done(n,fjac,ldjac,fmas,ldmas,mlmas,mumas,
1191 : & m1,m2,nm1,fac,e1,lde,ip,p1,ier,ijob,implct,ip,
1192 : & mle,mue,mbjac,mbb,mdiag,mdiff,mbdiag,
1193 : & decsol,decsols,decsolblk,
1194 : & caller_id, nvar, nz, lblk, dblk, ublk,
1195 39 : & sparse_jac,nzmax,isparse,ia,ja,sa,lrd,rpar_decsol,lid,ipar_decsol)
1196 : end if
1197 :
1198 39 : return
1199 39 : end subroutine roscor
1200 :
1201 :
1202 : !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1203 3 : subroutine ros2_coeffs (ns,ra,rc,rd,ros_d,ros_m,ros_e,ros_alpha,
1204 3 : & ros_gamma,ros_newf,ros_elo,no_aux_in_error,ros_name)
1205 :
1206 : ! Rosenbrock Method "ROS2"
1207 : ! CWI, MAS-R9717
1208 : ! J.G.Verwer, E.J.Spee, J.G.Blom, W.H.Hundsdorfer
1209 :
1210 : ! --- an L-stable method, 2 stages, order 2
1211 : !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1212 : implicit none
1213 : integer, intent(in) :: ns
1214 : real(dp), intent(inout) :: ros_m(ns), ros_e(ns)
1215 : real(dp), intent(inout) :: ros_d(ns,ns), ra(ns,ns), rc(ns,ns), rd(ns,ns)
1216 : real(dp), intent(inout) :: ros_alpha(ns), ros_gamma(ns)
1217 : integer, intent(out) :: ros_elo
1218 : logical, intent(out) :: no_aux_in_error, ros_newf(ns)
1219 : character (len=12), intent(out) :: ros_name
1220 :
1221 3 : real(dp) :: g
1222 3 : no_aux_in_error = .true.
1223 3 : g = 1.0d0 + 1.0d0/sqrt(2.0d0)
1224 :
1225 : !~~~> name of the method
1226 3 : ros_name = 'ros2'
1227 : !~~~> number of stages
1228 3 : if (ns /= 2) stop 'bad ns arg for ros2_coeffs'
1229 :
1230 21 : rd = 0
1231 21 : ros_d = 0
1232 :
1233 3 : ra(2,1) = (1.d0)/g
1234 3 : rc(2,1) = (-2.d0)/g
1235 :
1236 : !~~~> does the stage i require a new function evaluation (ros_newf(i)=true)
1237 : ! or does it reuse the function evaluation from stage i-1 (ros_newf(i)=false)
1238 3 : ros_newf(1) = .true.
1239 3 : ros_newf(2) = .true.
1240 : !~~~> m_i = coefficients for new step solution
1241 3 : ros_m(1)= (3.d0)/(2.d0*g)
1242 3 : ros_m(2)= (1.d0)/(2.d0*g)
1243 : ! e_i = coefficients for error estimator
1244 3 : ros_e(1) = 1.d0/(2.d0*g)
1245 3 : ros_e(2) = 1.d0/(2.d0*g)
1246 : !~~~> ros_elo = estimator of local order
1247 3 : ros_elo = 2
1248 : !~~~> y_stage_i ~ y( t + h*alpha_i )
1249 3 : ros_alpha(1) = 0.0d0
1250 3 : ros_alpha(2) = 1.0d0
1251 : !~~~> gamma_i = \sum_j gamma_{i,j}
1252 3 : ros_gamma(1) = g
1253 3 : ros_gamma(2) =-g
1254 :
1255 3 : return
1256 : end subroutine ros2_coeffs
1257 :
1258 :
1259 : !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1260 3 : subroutine rose2_coeffs (ns,ra,rc,rd,ros_d,ros_m,ros_e,ros_alpha,
1261 3 : & ros_gamma,ros_newf,ros_elo,no_aux_in_error,ros_name)
1262 :
1263 : ! Rosenbrock Method "ROSE2"
1264 : ! Shampine & Reichelt, SIAM J Sci. Comput., 18, (1997) 1-22.
1265 : ! ode23s from matlab ode suite.
1266 : ! --- an l-stable method, 3 stages, order 2
1267 : ! --- final function evaluation uses the new state.
1268 : !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1269 : implicit none
1270 : integer, intent(in) :: ns
1271 : real(dp), intent(inout) :: ros_m(ns), ros_e(ns), ros_d(ns,ns), ra(ns,ns), rc(ns,ns), rd(ns,ns)
1272 : real(dp), intent(inout) :: ros_alpha(ns), ros_gamma(ns)
1273 : integer, intent(out) :: ros_elo
1274 : logical, intent(out) :: no_aux_in_error, ros_newf(ns)
1275 : character (len=12), intent(out) :: ros_name
1276 3 : real(dp) :: g, e32
1277 : real(dp), parameter :: sqrt2 = 1.4142135623731d0 ! sqrt(2d0)
1278 3 : no_aux_in_error = .true.
1279 3 : g = 1.0d0/(2d0 + sqrt2)
1280 3 : e32 = 6d0 + sqrt2
1281 : !~~~> name of the method
1282 3 : ros_name = 'rose2'
1283 : !~~~> number of stages
1284 3 : if (ns /= 3) stop 'bad ns arg for rose2_coeffs'
1285 39 : ros_d = 0
1286 :
1287 : !~~~> the coefficient matrices a and c are strictly lower triangular.
1288 :
1289 3 : ra(2,1)= 0.5d0/g
1290 3 : ra(3,1)= 1.0d0/g
1291 3 : ra(3,2)= 1.0d0/g
1292 :
1293 3 : rc(2,1) = -1.0d0/g
1294 3 : rc(3,1) = -(e32+2d0)/g
1295 3 : rc(3,2) = -e32/g
1296 :
1297 3 : rd(2,1) = 0d0
1298 3 : rd(3,1) = 2d0
1299 3 : rd(3,2) = e32
1300 :
1301 : !~~~> does the stage i require a new function evaluation (ros_newf(i)=true)
1302 : ! or does it reuse the function evaluation from stage i-1 (ros_newf(i)=false)
1303 3 : ros_newf(1) = .true.
1304 3 : ros_newf(2) = .true.
1305 3 : ros_newf(3) = .true.
1306 :
1307 : !~~~> m_i = coefficients for new step solution
1308 3 : ros_m(1) = 1/g
1309 3 : ros_m(2) = 1/g
1310 3 : ros_m(3) = 0
1311 :
1312 : ! e_i = coefficients for error estimator
1313 3 : ros_e(1) = -1/(6*g)
1314 3 : ros_e(2) = -1/(3*g)
1315 3 : ros_e(3) = 1/(6*g)
1316 :
1317 : !~~~> ros_elo = estimator of local order
1318 3 : ros_elo = 2
1319 :
1320 : !~~~> y_stage_i ~ y( t + h*alpha_i )
1321 3 : ros_alpha(1)= 0.0d0
1322 3 : ros_alpha(2)= 0.5d0
1323 3 : ros_alpha(3)= 1.0d0
1324 :
1325 : !~~~> gamma_i = \sum_j gamma_{i,j}
1326 3 : ros_gamma(1)= g
1327 3 : ros_gamma(2)= 0
1328 3 : ros_gamma(3)= g
1329 :
1330 3 : return
1331 : end subroutine rose2_coeffs
1332 :
1333 :
1334 : !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1335 3 : subroutine ros3p_coeffs (ns,ra,rc,rd,ros_d,ros_m,ros_e,ros_alpha,
1336 3 : & ros_gamma,ros_newf,ros_elo,no_aux_in_error,ros_name)
1337 : !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1338 :
1339 : ! Rosenbrock Method "ROS3P"
1340 : ! Visit at CWI, January 2000
1341 : ! J.Lang, J.G.Verwer
1342 :
1343 : ! --- A-stable, 3 stages, order 3, 2 function evaluations
1344 : !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1345 : implicit none
1346 : integer, intent(in) :: ns
1347 : real(dp), intent(inout) :: ros_m(ns), ros_e(ns), ros_d(ns,ns), ra(ns,ns), rc(ns,ns), rd(ns,ns)
1348 : real(dp), intent(inout) :: ros_alpha(ns), ros_gamma(ns)
1349 : integer, intent(out) :: ros_elo
1350 : logical, intent(out) :: no_aux_in_error, ros_newf(ns)
1351 : character (len=12), intent(out) :: ros_name
1352 3 : no_aux_in_error = .true.
1353 : !~~~> name of the method
1354 3 : ros_name = 'ros3p'
1355 : !~~~> number of stages
1356 3 : if (ns /= 3) stop 'bad ns arg for ros3p_coeffs'
1357 39 : ra = 0
1358 39 : rc = 0
1359 39 : rd = 0
1360 39 : ros_d = 0
1361 :
1362 : !~~~> the coefficient matrices a and c are strictly lower triangular.
1363 :
1364 3 : ra(2,1)= 1.26794919243112273d0
1365 3 : ra(3,1)= 1.26794919243112273d0
1366 3 : ra(3,2)= 0.d0
1367 :
1368 3 : rc(2,1) = -1.60769515458673630d0
1369 3 : rc(3,1) = -3.46410161513775469d0
1370 3 : rc(3,2) = -1.73205080756887734d0
1371 :
1372 : !~~~> does the stage i require a new function evaluation (ros_newf(i)=true)
1373 : ! or does it reuse the function evaluation from stage i-1 (ros_newf(i)=false)
1374 3 : ros_newf(1) = .true.
1375 3 : ros_newf(2) = .true.
1376 3 : ros_newf(3) = .false.
1377 : !~~~> m_i = coefficients for new step solution
1378 3 : ros_m(1) = 2.0d0
1379 3 : ros_m(2) = 0.57735026918962578d0
1380 3 : ros_m(3) = 0.42264973081037424d0
1381 : ! e_i = coefficients for error estimator
1382 3 : ros_e(1) = ros_m(1) - 2.11324865405187124d0
1383 3 : ros_e(2) = ros_m(2) - 1.0d0
1384 3 : ros_e(3) = ros_m(3) - 0.42264973081037424d0
1385 : !~~~> ros_elo = estimator of local order
1386 3 : ros_elo = 3
1387 : !~~~> y_stage_i ~ y( t + h*alpha_i )
1388 3 : ros_alpha(1)= 0.0d+00
1389 3 : ros_alpha(2)= 1d+00
1390 3 : ros_alpha(3)= 1d+00
1391 : !~~~> gamma_i = \sum_j gamma_{i,j}
1392 3 : ros_gamma(1)= 0.78867513459481286d0
1393 3 : ros_gamma(2)= -0.21132486540518713d0
1394 3 : ros_gamma(3)= -1.07735026918962573d0
1395 3 : return
1396 : end subroutine ros3p_coeffs
1397 :
1398 : !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1399 :
1400 : !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1401 7 : subroutine ros3pl_coeffs (ns,ra,rc,rd,ros_d,ros_m,ros_e,ros_alpha,
1402 7 : & ros_gamma,ros_newf,ros_elo,no_aux_in_error,ros_name)
1403 : !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1404 : ! --- a stiffly-stable method for parabolic equations; 4 stages, order 3, 3 function evaluations.
1405 : !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1406 : implicit none
1407 : integer, intent(in) :: ns
1408 : real(dp), intent(inout) :: ros_m(ns), ros_e(ns), ros_d(ns,ns), ra(ns,ns), rc(ns,ns), rd(ns,ns)
1409 : real(dp), intent(inout) :: ros_alpha(ns), ros_gamma(ns)
1410 : integer, intent(out) :: ros_elo
1411 : logical, intent(out) :: no_aux_in_error, ros_newf(ns)
1412 : character (len=12), intent(out) :: ros_name
1413 7 : no_aux_in_error = .true.
1414 : !~~~> name of the method
1415 7 : ros_name = 'ros3pl'
1416 : !~~~> number of stages
1417 7 : if (ns /= 4) stop 'bad ns arg for ros3pl_coeffs'
1418 147 : ra = 0
1419 147 : rc = 0
1420 147 : rd = 0
1421 147 : ros_d = 0
1422 :
1423 : !~~~> the coefficient matrices a and c are strictly lower triangular.
1424 :
1425 : ! matA (alpha21, alpha31, ...)
1426 7 : ra(2,1) = 1.147140180139521d0
1427 7 : ra(3,1) = 2.463070773030053d0
1428 7 : ra(3,2) = ra(2,1)
1429 7 : ra(4,1) = ra(3,1)
1430 7 : ra(4,2) = ra(2,1)
1431 7 : ra(4,3) = 0.d0
1432 :
1433 : ! -matC (c21, c31, ...)
1434 7 : rc(2,1) = -2.631861185781065d0
1435 7 : rc(3,1) = -1.302364158113095d0
1436 7 : rc(3,2) = 2.769432022251304d0
1437 7 : rc(4,1) = -1.552568958732400d0
1438 7 : rc(4,2) = 2.587743501215153d0
1439 7 : rc(4,3) = -1.416993298352020d0
1440 :
1441 : !~~~> does the stage i require a new function evaluation (ros_newf(i)=true)
1442 : ! or does it reuse the function evaluation from stage i-1 (ros_newf(i)=false)
1443 7 : ros_newf(1) = .true.
1444 7 : ros_newf(2) = .true.
1445 7 : ros_newf(3) = .true.
1446 7 : ros_newf(4) = .false.
1447 :
1448 : !~~~> m_i = coefficients for new step solution
1449 : ! vecB1 (m1, m2, ...)
1450 7 : ros_m(1) = 2.463070773030053d0
1451 7 : ros_m(2) = 1.147140180139521d0
1452 7 : ros_m(3) = 0
1453 7 : ros_m(4) = 1
1454 :
1455 : ! e_i = coefficients for error estimator
1456 : ! vecB1-vecB2 (m1-mhat1, ...)
1457 7 : ros_e(1) = ros_m(1) - 2.346947683513665d0
1458 7 : ros_e(2) = ros_m(2) - 0.456530569451895d0
1459 7 : ros_e(3) = ros_m(3) - 0.056949243945495d0
1460 7 : ros_e(4) = ros_m(4) - 0.738684936166224d0
1461 :
1462 : !~~~> ros_elo = estimator of local order
1463 7 : ros_elo = 3
1464 :
1465 : !~~~> y_stage_i ~ y( t + h*alpha_i )
1466 : ! vecA (alpha1, alpha2, ...)
1467 7 : ros_alpha(1)= 0
1468 7 : ros_alpha(2)= 0.5d0
1469 7 : ros_alpha(3)= 1
1470 7 : ros_alpha(4)= 1
1471 :
1472 : !~~~> gamma_i = \sum_j gamma_{i,j}
1473 : ! vecG (gamma1, gamma2,...)
1474 7 : ros_gamma(1)= 0.435866521508459d0
1475 7 : ros_gamma(2)= -0.064133478491541d0
1476 7 : ros_gamma(3)= 0.111028172512505d0
1477 7 : ros_gamma(4)= 0
1478 :
1479 7 : return
1480 : end subroutine ros3pl_coeffs
1481 :
1482 : !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1483 :
1484 : !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1485 7 : subroutine rodas3_coeffs (ns,ra,rc,rd,ros_d,ros_m,ros_e,ros_alpha,
1486 7 : & ros_gamma,ros_newf,ros_elo,no_aux_in_error,ros_name)
1487 : !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1488 : ! --- a stiffly-stable method; 4 stages, order 3, 3 function evaluations.
1489 : !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1490 : implicit none
1491 : integer, intent(in) :: ns
1492 : real(dp), intent(inout) :: ros_m(ns), ros_e(ns), ros_d(ns,ns), ra(ns,ns), rc(ns,ns), rd(ns,ns)
1493 : real(dp), intent(inout) :: ros_alpha(ns), ros_gamma(ns)
1494 : integer, intent(out) :: ros_elo
1495 : logical, intent(out) :: no_aux_in_error, ros_newf(ns)
1496 : character (len=12), intent(out) :: ros_name
1497 147 : ra = 0
1498 147 : rc = 0
1499 147 : rd = 0
1500 147 : ros_d = 0
1501 7 : no_aux_in_error = .false.
1502 : !~~~> name of the method
1503 7 : ros_name = 'rodas3'
1504 : !~~~> number of stages
1505 7 : if (ns /= 4) stop 'bad ns arg for rodas3_coeffs'
1506 :
1507 : !~~~> the coefficient matrices a and c are strictly lower triangular.
1508 7 : ra(2,1) = 0.0d+00
1509 7 : ra(3,1) = 2.0d+00
1510 7 : ra(3,2) = 0.0d+00
1511 7 : ra(4,1) = 2.0d+00
1512 7 : ra(4,2) = 0.0d+00
1513 7 : ra(4,3) = 1.0d+00
1514 :
1515 7 : rc(2,1) = 4.0d+00
1516 7 : rc(3,1) = 1.0d+00
1517 7 : rc(3,2) =-1.0d+00
1518 7 : rc(4,1) = 1.0d+00
1519 7 : rc(4,2) =-1.0d+00
1520 7 : rc(4,3) =-(8.0d+00/3.0d+00)
1521 :
1522 : !~~~> does the stage i require a new function evaluation (ros_newf(i)=true)
1523 : ! or does it reuse the function evaluation from stage i-1 (ros_newf(i)=false)
1524 7 : ros_newf(1) = .true.
1525 7 : ros_newf(2) = .false.
1526 7 : ros_newf(3) = .true.
1527 7 : ros_newf(4) = .true.
1528 : !~~~> m_i = coefficients for new step solution
1529 7 : ros_m(1) = 2.0d+00
1530 7 : ros_m(2) = 0.0d+00
1531 7 : ros_m(3) = 1.0d+00
1532 7 : ros_m(4) = 1.0d+00
1533 : !~~~> e_i = coefficients for error estimator
1534 7 : ros_e(1) = 0.0d+00
1535 7 : ros_e(2) = 0.0d+00
1536 7 : ros_e(3) = 0.0d+00
1537 7 : ros_e(4) = 1.0d+00
1538 : !~~~> ros_elo = estimator of local order
1539 7 : ros_elo = 3
1540 : !~~~> y_stage_i ~ y( t + h*alpha_i )
1541 7 : ros_alpha(1) = 0.0d+00
1542 7 : ros_alpha(2) = 0.0d+00
1543 7 : ros_alpha(3) = 1.0d+00
1544 7 : ros_alpha(4) = 1.0d+00
1545 : !~~~> gamma_i = \sum_j gamma_{i,j}
1546 7 : ros_gamma(1) = 0.5d+00
1547 7 : ros_gamma(2) = 1.5d+00
1548 7 : ros_gamma(3) = 0.0d+00
1549 7 : ros_gamma(4) = 0.0d+00
1550 7 : return
1551 : end subroutine rodas3_coeffs
1552 :
1553 :
1554 : !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1555 9 : subroutine rodas4_coeffs (ns,ra,rc,rd,ros_d,ros_m,ros_e,ros_alpha,
1556 9 : & ros_gamma,ros_newf,ros_elo,no_aux_in_error,ros_name)
1557 : !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1558 : ! stiffly-stable rosenbrock method of order 4, with 6 stages
1559 : !
1560 : ! e. hairer and g. wanner, solving ordinary differential
1561 : ! equations ii. stiff and differential-algebraic problems.
1562 : ! springer series in computational mathematics,
1563 : ! springer-verlag (1996)
1564 : !
1565 : !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1566 :
1567 : implicit none
1568 : integer, intent(in) :: ns
1569 : real(dp), intent(inout) ::
1570 : & ros_m(ns), ros_e(ns), ros_d(ns,ns), ra(ns,ns), rc(ns,ns), rd(ns,ns)
1571 : real(dp), intent(inout) :: ros_alpha(ns), ros_gamma(ns)
1572 : integer, intent(out) :: ros_elo
1573 : logical, intent(out) :: no_aux_in_error, ros_newf(ns)
1574 : character (len=12), intent(out) :: ros_name
1575 :
1576 9 : no_aux_in_error = .false.
1577 387 : rd = 0
1578 : !~~~> name of the method
1579 9 : ros_name = 'rodas4'
1580 : !~~~> number of stages
1581 9 : if (ns /= 6) stop 'bad ns arg for rodas4_coeffs'
1582 :
1583 : !~~~> y_stage_i ~ y( t + h*alpha_i )
1584 9 : ros_alpha(1) = 0.000d0
1585 9 : ros_alpha(2) = 0.386d0
1586 9 : ros_alpha(3) = 0.210d0
1587 9 : ros_alpha(4) = 0.630d0
1588 9 : ros_alpha(5) = 1.000d0
1589 9 : ros_alpha(6) = 1.000d0
1590 :
1591 : !~~~> gamma_i = \sum_j gamma_{i,j}
1592 9 : ros_gamma(1) = 0.2500000000000000d+00
1593 9 : ros_gamma(2) =-0.1043000000000000d+00
1594 9 : ros_gamma(3) = 0.1035000000000000d+00
1595 9 : ros_gamma(4) =-0.3620000000000023d-01
1596 9 : ros_gamma(5) = 0.0d0
1597 9 : ros_gamma(6) = 0.0d0
1598 :
1599 : !~~~> the coefficient matrices a and c are strictly lower triangular.
1600 :
1601 9 : ra(2,1) = 0.1544000000000000d+01
1602 9 : ra(3,1) = 0.9466785280815826d+00
1603 9 : ra(3,2) = 0.2557011698983284d+00
1604 9 : ra(4,1) = 0.3314825187068521d+01
1605 9 : ra(4,2) = 0.2896124015972201d+01
1606 9 : ra(4,3) = 0.9986419139977817d+00
1607 9 : ra(5,1) = 0.1221224509226641d+01
1608 9 : ra(5,2) = 0.6019134481288629d+01
1609 9 : ra(5,3) = 0.1253708332932087d+02
1610 9 : ra(5,4) =-0.6878860361058950d+00
1611 :
1612 45 : ra(6,1:4) = ra(5,1:4)
1613 9 : ra(6,5) = 1
1614 :
1615 9 : rc(2,1) =-0.5668800000000000d+01
1616 9 : rc(3,1) =-0.2430093356833875d+01
1617 9 : rc(3,2) =-0.2063599157091915d+00
1618 9 : rc(4,1) =-0.1073529058151375d+00
1619 9 : rc(4,2) =-0.9594562251023355d+01
1620 9 : rc(4,3) =-0.2047028614809616d+02
1621 9 : rc(5,1) = 0.7496443313967647d+01
1622 9 : rc(5,2) =-0.1024680431464352d+02
1623 9 : rc(5,3) =-0.3399990352819905d+02
1624 9 : rc(5,4) = 0.1170890893206160d+02
1625 9 : rc(6,1) = 0.8083246795921522d+01
1626 9 : rc(6,2) =-0.7981132988064893d+01
1627 9 : rc(6,3) =-0.3152159432874371d+02
1628 9 : rc(6,4) = 0.1631930543123136d+02
1629 9 : rc(6,5) =-0.6058818238834054d+01
1630 :
1631 : !~~~> m_i = coefficients for new step solution
1632 45 : ros_m(1:4) = ra(5,1:4)
1633 9 : ros_m(5) = 1
1634 9 : ros_m(6) = 1
1635 :
1636 : !~~~> e_i = coefficients for error estimator
1637 9 : ros_e(1) = 0.0d+00
1638 9 : ros_e(2) = 0.0d+00
1639 9 : ros_e(3) = 0.0d+00
1640 9 : ros_e(4) = 0.0d+00
1641 9 : ros_e(5) = 0.0d+00
1642 9 : ros_e(6) = 1.0d+00
1643 :
1644 : !~~~> does the stage i require a new function evaluation (ros_newf(i)=true)
1645 : ! or does it reuse the function evaluation from stage i-1 (ros_newf(i)=false)
1646 9 : ros_newf(1) = .true.
1647 9 : ros_newf(2) = .true.
1648 9 : ros_newf(3) = .true.
1649 9 : ros_newf(4) = .true.
1650 9 : ros_newf(5) = .true.
1651 9 : ros_newf(6) = .true.
1652 :
1653 : ! dense output
1654 9 : ros_d(2,1)= 0.1012623508344586d+02
1655 9 : ros_d(2,2)=-0.7487995877610167d+01
1656 9 : ros_d(2,3)=-0.3480091861555747d+02
1657 9 : ros_d(2,4)=-0.7992771707568823d+01
1658 9 : ros_d(2,5)= 0.1025137723295662d+01
1659 9 : ros_d(3,1)=-0.6762803392801253d+00
1660 9 : ros_d(3,2)= 0.6087714651680015d+01
1661 9 : ros_d(3,3)= 0.1643084320892478d+02
1662 9 : ros_d(3,4)= 0.2476722511418386d+02
1663 9 : ros_d(3,5)=-0.6594389125716872d+01
1664 :
1665 : !~~~> ros_elo = estimator of local order
1666 9 : ros_elo = 4
1667 :
1668 9 : return
1669 : end subroutine rodas4_coeffs
1670 :
1671 :
1672 : !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1673 7 : subroutine rodasp_coeffs (ns,ra,rc,rd,ros_d,ros_m,ros_e,ros_alpha,
1674 7 : & ros_gamma,ros_newf,ros_elo,no_aux_in_error,ros_name)
1675 : !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1676 : ! stiffly-stable rosenbrock method of order 4, with 6 stages
1677 : ! for parabolic equations (G.Steinbach,1993).
1678 : !
1679 : ! e. hairer and g. wanner, solving ordinary differential
1680 : ! equations ii. stiff and differential-algebraic problems.
1681 : ! springer series in computational mathematics,
1682 : ! springer-verlag (1996)
1683 : !
1684 : !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1685 :
1686 : implicit none
1687 : integer, intent(in) :: ns
1688 : real(dp), intent(inout) :: ros_m(ns), ros_e(ns), ros_d(ns,ns), ra(ns,ns), rc(ns,ns), rd(ns,ns)
1689 : real(dp), intent(inout) :: ros_alpha(ns), ros_gamma(ns)
1690 : integer, intent(out) :: ros_elo
1691 : logical, intent(out) :: no_aux_in_error, ros_newf(ns)
1692 : character (len=12), intent(out) :: ros_name
1693 7 : no_aux_in_error = .false.
1694 301 : rd = 0
1695 : !~~~> name of the method
1696 7 : ros_name = 'rodasp'
1697 : !~~~> number of stages
1698 7 : if (ns /= 6) stop 'bad ns arg for rodasp_coeffs'
1699 :
1700 : !~~~> y_stage_i ~ y( t + h*alpha_i )
1701 7 : ros_alpha(1) = 0.000d0
1702 7 : ros_alpha(2) = 0.750d0
1703 7 : ros_alpha(3) = 0.210d0
1704 7 : ros_alpha(4) = 0.630d0
1705 7 : ros_alpha(5) = 1.000d0
1706 7 : ros_alpha(6) = 1.000d0
1707 :
1708 : !~~~> gamma_i = \sum_j gamma_{i,j}
1709 7 : ros_gamma(1) = 0.2500000000000000d+00
1710 7 : ros_gamma(2) =-0.5000000000000000d+00
1711 7 : ros_gamma(3) =-0.2350400000000000d-01
1712 7 : ros_gamma(4) =-0.3620000000000000d-01
1713 7 : ros_gamma(5) = 0.0d0
1714 7 : ros_gamma(6) = 0.0d0
1715 :
1716 : !~~~> the coefficient matrices a and c are strictly lower triangular.
1717 7 : ra(2,1) = 0.3000000000000000d+01
1718 7 : ra(3,1) = 0.1831036793486759d+01
1719 7 : ra(3,2) = 0.4955183967433795d+00
1720 7 : ra(4,1) = 0.2304376582692669d+01
1721 7 : ra(4,2) =-0.5249275245743001d-01
1722 7 : ra(4,3) =-0.1176798761832782d+01
1723 7 : ra(5,1) =-0.7170454962423024d+01
1724 7 : ra(5,2) =-0.4741636671481785d+01
1725 7 : ra(5,3) =-0.1631002631330971d+02
1726 7 : ra(5,4) =-0.1062004044111401d+01
1727 :
1728 35 : ra(6,1:4) = ra(5,1:4)
1729 7 : ra(6,5) = 1
1730 :
1731 7 : rc(2,1)=-0.1200000000000000d+02
1732 7 : rc(3,1)=-0.8791795173947035d+01
1733 7 : rc(3,2)=-0.2207865586973518d+01
1734 7 : rc(4,1)= 0.1081793056857153d+02
1735 7 : rc(4,2)= 0.6780270611428266d+01
1736 7 : rc(4,3)= 0.1953485944642410d+02
1737 7 : rc(5,1)= 0.3419095006749676d+02
1738 7 : rc(5,2)= 0.1549671153725963d+02
1739 7 : rc(5,3)= 0.5474760875964130d+02
1740 7 : rc(5,4)= 0.1416005392148534d+02
1741 7 : rc(6,1)= 0.3462605830930532d+02
1742 7 : rc(6,2)= 0.1530084976114473d+02
1743 7 : rc(6,3)= 0.5699955578662667d+02
1744 7 : rc(6,4)= 0.1840807009793095d+02
1745 7 : rc(6,5)=-0.5714285714285717d+01
1746 :
1747 : !~~~> m_i = coefficients for new step solution
1748 35 : ros_m(1:4) = ra(5,1:4)
1749 7 : ros_m(5) = 1
1750 7 : ros_m(6) = 1
1751 :
1752 : !~~~> e_i = coefficients for error estimator
1753 7 : ros_e(1) = 0.0d+00
1754 7 : ros_e(2) = 0.0d+00
1755 7 : ros_e(3) = 0.0d+00
1756 7 : ros_e(4) = 0.0d+00
1757 7 : ros_e(5) = 0.0d+00
1758 7 : ros_e(6) = 1.0d+00
1759 :
1760 : !~~~> does the stage i require a new function evaluation (ros_newf(i)=true)
1761 : ! or does it reuse the function evaluation from stage i-1 (ros_newf(i)=false)
1762 7 : ros_newf(1) = .true.
1763 7 : ros_newf(2) = .true.
1764 7 : ros_newf(3) = .true.
1765 7 : ros_newf(4) = .true.
1766 7 : ros_newf(5) = .true.
1767 7 : ros_newf(6) = .true.
1768 :
1769 : ! dense output
1770 7 : ros_d(2,1)= 0.2509876703708589d+02
1771 7 : ros_d(2,2)= 0.1162013104361867d+02
1772 7 : ros_d(2,3)= 0.2849148307714626d+02
1773 7 : ros_d(2,4)=-0.5664021568594133d+01
1774 7 : ros_d(2,5)= 0.0000000000000000d+00
1775 7 : ros_d(3,1)= 0.1638054557396973d+01
1776 7 : ros_d(3,2)=-0.7373619806678748d+00
1777 7 : ros_d(3,3)= 0.8477918219238990d+01
1778 7 : ros_d(3,4)= 0.1599253148779520d+02
1779 7 : ros_d(3,5)=-0.1882352941176471d+01
1780 :
1781 : !~~~> ros_elo = estimator of local order
1782 7 : ros_elo = 4
1783 :
1784 7 : return
1785 : end subroutine rodasp_coeffs
1786 :
1787 : ! continuous output routines
1788 :
1789 :
1790 0 : real(dp) function contro3(i,x,rwork,iwork,ierr)
1791 : integer, intent(in) :: i
1792 : real(dp), intent(in) :: x
1793 : real(dp), intent(inout), target :: rwork(*)
1794 : integer, intent(inout), target :: iwork(*)
1795 : integer, intent(out) :: ierr
1796 :
1797 0 : real(dp) :: xold, h, dx, s, s2
1798 0 : real(dp), pointer :: cont(:)
1799 : integer :: lrc, n, j
1800 0 : ierr=0
1801 0 : lrc = iwork(1)
1802 0 : n = iwork(2)
1803 0 : cont => rwork(1:lrc); j=lrc+1
1804 0 : xold = rwork(j); j=j+1
1805 0 : h = rwork(j)
1806 0 : dx=x-xold
1807 0 : s = dx/h
1808 0 : s2 = s**2
1809 0 : contro3 = (1-s2)*cont(i) + dx*(1-s)*cont(i+2*n) + s2*cont(i+n)
1810 : return
1811 0 : end function contro3
1812 :
1813 :
1814 0 : real(dp) function contro4(i,x,rwork,iwork,ierr)
1815 : integer, intent(in) :: i
1816 : real(dp), intent(in) :: x
1817 : real(dp), intent(inout), target :: rwork(*)
1818 : integer, intent(inout), target :: iwork(*)
1819 : integer, intent(out) :: ierr
1820 :
1821 0 : real(dp) :: xold, h, s
1822 0 : real(dp), pointer :: cont(:)
1823 : integer :: lrc, n, j
1824 0 : ierr=0
1825 0 : lrc = iwork(1)
1826 0 : n = iwork(2)
1827 :
1828 0 : cont => rwork(1:lrc); j=1+lrc
1829 0 : xold = rwork(j); j=j+1
1830 0 : h = rwork(j)
1831 0 : s=(x-xold)/h
1832 0 : contro4=cont(i)*(1-s)+s*(cont(i+n)+(1-s)*(cont(i+n*2)+s*cont(i+n*3)))
1833 : return
1834 0 : end function contro4
1835 :
1836 : end module mod_rosenbrock
1837 :
1838 :
1839 : ! copyright (c) 2004, ernst hairer
1840 :
1841 : ! redistribution and use in source and binary forms, with or without
1842 : ! modification, are permitted provided that the following conditions are
1843 : ! met:
1844 :
1845 : ! - redistributions of source code must retain the above copyright
1846 : ! notice, this list of conditions and the following disclaimer.
1847 :
1848 : ! - redistributions in binary form must reproduce the above copyright
1849 : ! notice, this list of conditions and the following disclaimer in the
1850 : ! documentation and/or other materials provided with the distribution.
1851 :
1852 : ! this software is provided by the copyright holders and contributors “as
1853 : ! is” and any express or implied warranties, including, but not limited
1854 : ! to, the implied warranties of merchantability and fitness for a
1855 : ! particular purpose are disclaimed. in no event shall the regents or
1856 : ! contributors be liable for any direct, indirect, incidental, special,
1857 : ! exemplary, or consequential damages (including, but not limited to,
1858 : ! procurement of substitute goods or services; loss of use, data, or
1859 : ! profits; or business interruption) however caused and on any theory of
1860 : ! liability, whether in contract, strict liability, or tort (including
1861 : ! negligence or otherwise) arising in any way out of the use of this
1862 : ! software, even if advised of the possibility of such damage.
|
