roots.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: roots.f90 #
5!# #
6!# Copyright (C) 2006-2018 #
7!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8!# #
9!# This program is free software; you can redistribute it and/or modify #
10!# it under the terms of the GNU General Public License as published by #
11!# the Free Software Foundation; either version 2 of the License, or (at #
12!# your option) any later version. #
13!# #
14!# This program is distributed in the hope that it will be useful, but #
15!# WITHOUT ANY WARRANTY; without even the implied warranty of #
16!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
17!# NON INFRINGEMENT. See the GNU General Public License for more #
18!# details. #
19!# #
20!# You should have received a copy of the GNU General Public License #
21!# along with this program; if not, write to the Free Software #
22!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
23!# #
24!#############################################################################
25
26!----------------------------------------------------------------------------!
44!----------------------------------------------------------------------------!
45MODULE roots
46 IMPLICIT NONE
47 !--------------------------------------------------------------------------!
48 PRIVATE
50 TYPE roots_typ
51 REAL :: eps,tol,root,dx,dxold,df,d,e
52 REAL :: x(3),f(3)
53 INTEGER :: iter,max_iterations,error
54 LOGICAL :: iterate,do_modified_step
55 END TYPE roots_typ
56 ! exclude interface block from doxygen processing
58 INTERFACE getroot
59 MODULE PROCEDURE getroot_brentdekker
60 END INTERFACE
62 !--------------------------------------------------------------------------!
63 ! some constants
64! #define DEBUG_OUTPUT 1
65#undef DEBUG_OUTPUT
66 INTEGER, PARAMETER :: default_max_iterations = 1000
67 REAL, PARAMETER :: default_accuracy = 4*epsilon(default_accuracy)
68 CHARACTER(LEN=64), PARAMETER, DIMENSION(0:4) :: error_message = (/ &
69 "unknown error ", &
70 "root not bracketed ", &
71 "iteration exceeds maximum ", &
72 "requested accuracy smaller than machine precission ", &
73 "upper limit for iterations should be larger than 1 " /)
74 !--------------------------------------------------------------------------!
75 PUBLIC :: &
76 ! constants
77 default_accuracy, &
78 default_max_iterations, &
79 ! types
80 roots_typ, &
81 ! methods
82 initroots, &
83 geterrormessage, &
84 getroot, &
85 getroot_newton, &
86 getroot_regulafalsi, &
87 getroot_bisection, &
88 getroot_pegasus, &
89 getroot_king, &
90 getroot_andersonbjoerk, &
91 getroot_ridder, &
92 getroot_brentdekker
93 !--------------------------------------------------------------------------!
94
95CONTAINS
96
97 PURE SUBROUTINE initroots(this,x1,x2,f1,f2,dxacc,maxiter)
98 IMPLICIT NONE
99 !------------------------------------------------------------------------!
100 TYPE(roots_typ), INTENT(INOUT) :: this
101 REAL, INTENT(IN) :: x1,x2
102 REAL, INTENT(IN) :: f1,f2
103 REAL,OPTIONAL, INTENT(IN) :: dxacc
104 INTEGER,OPTIONAL, INTENT(IN) :: maxiter
105 !------------------------------------------------------------------------!
106 ! check if root is bracketed
107 IF (f1*f2.GT.0.0) THEN
108 ! f1 and f2 should have opposite signs
109 this%error = 1
110 ELSE
111 IF (x1.LT.x2) THEN
112 this%x(1) = x1
113 this%x(2) = x2
114 this%f(1) = f1
115 this%f(2) = f2
116 ELSE
117 this%x(1) = x2
118 this%x(2) = x1
119 this%f(1) = f2
120 this%f(2) = f1
121 END IF
122
123 this%dx = abs(this%x(2)-this%x(1))
124 this%dxold = this%dx
125
126 IF(PRESENT(dxacc)) THEN
127 IF (dxacc.GT.epsilon(dxacc)) THEN
128 this%eps = dxacc
129 ELSE
130 this%error = 3
131 END IF
132 ELSE
133 this%eps = default_accuracy
134 END IF
135
136 IF(PRESENT(maxiter)) THEN
137 IF (maxiter.GT.1) THEN
138 this%max_iterations = maxiter
139 ELSE
140 this%error = 4
141 END IF
142 ELSE
143 this%max_iterations = default_max_iterations
144 END IF
145
146 ! only used in Brent-Dekker method
147 this%d = this%x(2)-this%x(1)
148 this%e = this%d
149
150 this%error = 0
151 this%iter = 0
152 this%iterate = .true.
153 this%do_modified_step = .true.
154 END IF
155
156 END SUBROUTINE initroots
157
158
159 FUNCTION geterrormessage(error) RESULT(msg)
160 IMPLICIT NONE
161 !------------------------------------------------------------------------!
162 INTEGER, INTENT(IN) :: error
163 CHARACTER(LEN=64) :: msg
164 !------------------------------------------------------------------------!
165 SELECT CASE (error)
166 CASE(1,2)
167 msg = error_message(error)
168 CASE DEFAULT
169 msg = error_message(0)
170 END SELECT
171 END FUNCTION geterrormessage
172
173
174#ifndef DEBUG_OUTPUT
175 PURE &
176#endif
177 SUBROUTINE getroot_generic(this,Stepper,func,x1,x2,dxacc,maxiter,plist,xm)
178 IMPLICIT NONE
179 !------------------------------------------------------------------------!
180 TYPE(roots_typ) :: this
181 REAL :: x1,x2
182 REAL, OPTIONAL :: dxacc, plist(:), xm
183 INTEGER, OPTIONAL :: maxiter
184 !------------------------------------------------------------------------!
185 INTENT(IN) :: x1,x2,dxacc,maxiter,xm,plist
186 INTENT(INOUT) :: this
187 !------------------------------------------------------------------------!
188 INTERFACE
189 PURE SUBROUTINE stepper(root)
190 IMPORT roots_typ
191 IMPLICIT NONE
192 TYPE(roots_typ), INTENT(INOUT) :: root
193 END SUBROUTINE stepper
194 END INTERFACE
195 INTERFACE
196 PURE SUBROUTINE func(x,fx,plist)
197 IMPLICIT NONE
198 REAL, INTENT(IN) :: x
199 REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
200 REAL, INTENT(OUT) :: fx
201 END SUBROUTINE func
202 END INTERFACE
203 !------------------------------------------------------------------------!
204 REAL :: f1,f2
205 !------------------------------------------------------------------------!
206 ! compute function values at the boundaries
207 CALL func(x1,f1,plist)
208 CALL func(x2,f2,plist)
209
210 ! initialize root finding algorithm
211 CALL initroots(this,x1,x2,f1,f2)
212
213 IF (this%error.EQ.0) THEN
214 ! initial guess for the root
215 IF (PRESENT(xm)) THEN
216 this%x(3) = xm
217 ELSE
218 this%x(3) = arithmeticmean(this%x(1),this%x(2))
219 END IF
220#ifdef DEBUG_OUTPUT
221 print '(A3,A18,A10,A18,A18,2(A13,A10))', &
222 "#","x3","f3","dx","tol","x1","f1","x2","f2"
223#endif
224 ! main loop
225 DO
226 ! compute function value at x3
227 CALL func(this%x(3),this%f(3),plist)
228 ! print debug information if requested
229#ifdef DEBUG_OUTPUT
230 print '(I3,ES18.10,ES10.2,2(ES18.10),2(ES13.5,ES10.2))', &
231 this%iter,this%x(3),this%f(3),this%dx,this%tol,this%x(1),this%f(1),this%x(2),this%f(2)
232#endif
233 ! check convergence
234 CALL testconvergence(this)
235 IF (.NOT.this%iterate) EXIT
236
237 ! determine new estimate for x3
238 CALL stepper(this)
239 END DO
240 END IF
241 END SUBROUTINE getroot_generic
242
243#ifndef DEBUG_OUTPUT
244 PURE &
245#endif
246 SUBROUTINE getroot_newton(funcd,x1,x2,root,error,plist,xm,iterations)
247 IMPLICIT NONE
248 !------------------------------------------------------------------------!
249 REAL, INTENT(IN) :: x1,x2
250 REAL, INTENT(OUT) :: root
251 INTEGER, INTENT(OUT) :: error
252 REAL, DIMENSION(:), INTENT(IN), OPTIONAL :: plist
253 REAL, INTENT(IN), OPTIONAL :: xm
254 INTEGER, INTENT(OUT), OPTIONAL :: iterations
255 !------------------------------------------------------------------------!
256 INTERFACE
257 PURE SUBROUTINE funcd(x,fx,dfx,plist)
258 IMPLICIT NONE
259 REAL, INTENT(IN) :: x
260 REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
261 REAL, INTENT(OUT) :: fx,dfx
262 END SUBROUTINE funcd
263 END INTERFACE
264 !------------------------------------------------------------------------!
265 TYPE(roots_typ) :: this
266 REAL :: f1,f2,df1,df2
267 !------------------------------------------------------------------------!
268 ! compute function values and derivatives at the boundaries
269 CALL funcd(x1,f1,df1,plist)
270 CALL funcd(x2,f2,df2,plist)
271
272 ! initialize root finding algorithm
273 CALL initroots(this,x1,x2,f1,f2)
274
275 IF (this%error.EQ.0) THEN
276 ! initial guess for the root
277 IF (PRESENT(xm)) THEN
278 this%x(3) = xm
279 ELSE
280 this%x(3) = arithmeticmean(this%x(1),this%x(2))
281 END IF
282#ifdef DEBUG_OUTPUT
283 print '(A3,A18,A10,A18,A18,2(A13,A10))', &
284 "#","x3","f3","dx","tol","x1","f1","x2","f2"
285#endif
286 ! main loop
287 DO
288 ! compute function value and derivative at x3
289 CALL funcd(this%x(3),this%f(3),this%df,plist)
290 ! print debug information if requested
291#ifdef DEBUG_OUTPUT
292 print '(I3,ES18.10,ES10.2,2(ES18.10),2(ES13.5,ES10.2))', &
293 this%iter,this%x(3),this%f(3),this%dx,this%tol,this%x(1),this%f(1),this%x(2),this%f(2)
294#endif
295 ! check convergence
296 CALL testconvergence(this)
297 IF (.NOT.this%iterate) EXIT
298
299 ! determine new interval for enclosure and estimate x3
300 CALL step_newton(this)
301 END DO
302 END IF
303
304 root = this%root
305 error = this%error
306 IF (PRESENT(iterations)) iterations = this%iter
307
308 END SUBROUTINE getroot_newton
309
310
311#ifndef DEBUG_OUTPUT
312 PURE &
313#endif
314 SUBROUTINE getroot_regulafalsi(func,x1,x2,root,error,plist,xm,iterations)
315 IMPLICIT NONE
316 !------------------------------------------------------------------------!
317 REAL, INTENT(IN) :: x1,x2
318 REAL, INTENT(OUT) :: root
319 INTEGER, INTENT(OUT) :: error
320 REAL, DIMENSION(:), INTENT(IN), OPTIONAL :: plist
321 REAL, INTENT(IN), OPTIONAL :: xm
322 INTEGER, INTENT(OUT), OPTIONAL :: iterations
323 !------------------------------------------------------------------------!
324 INTERFACE
325 PURE SUBROUTINE func(x,fx,plist)
326 IMPLICIT NONE
327 REAL, INTENT(IN) :: x
328 REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
329 REAL, INTENT(OUT) :: fx
330 END SUBROUTINE func
331 END INTERFACE
332 !------------------------------------------------------------------------!
333 TYPE(roots_typ) :: this
334 !------------------------------------------------------------------------!
335 CALL getroot_generic(this,step_regulafalsi,func,x1,x2, &
336 default_accuracy,default_max_iterations,plist,xm)
337
338 root = this%root
339 error = this%error
340 IF (PRESENT(iterations)) iterations = this%iter
341
342 END SUBROUTINE getroot_regulafalsi
343
344
345#ifndef DEBUG_OUTPUT
346 PURE &
347#endif
348 SUBROUTINE getroot_pegasus(func,x1,x2,root,error,plist,xm,iterations)
349 IMPLICIT NONE
350 !------------------------------------------------------------------------!
351 REAL, INTENT(IN) :: x1,x2
352 REAL, INTENT(OUT) :: root
353 INTEGER, INTENT(OUT) :: error
354 REAL, DIMENSION(:), INTENT(IN), OPTIONAL :: plist
355 REAL, INTENT(IN), OPTIONAL :: xm
356 INTEGER, INTENT(OUT), OPTIONAL :: iterations
357 !------------------------------------------------------------------------!
358 INTERFACE
359 PURE SUBROUTINE func(x,fx,plist)
360 IMPLICIT NONE
361 REAL, INTENT(IN) :: x
362 REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
363 REAL, INTENT(OUT) :: fx
364 END SUBROUTINE func
365 END INTERFACE
366 !------------------------------------------------------------------------!
367 TYPE(roots_typ) :: this
368 !------------------------------------------------------------------------!
369 CALL getroot_generic(this,step_pegasus,func,x1,x2, &
370 default_accuracy,default_max_iterations,plist,xm)
371
372 root = this%root
373 error = this%error
374 IF (PRESENT(iterations)) iterations = this%iter
375
376 END SUBROUTINE getroot_pegasus
377
378
379#ifndef DEBUG_OUTPUT
380 PURE &
381#endif
382 SUBROUTINE getroot_king(func,x1,x2,root,error,plist,xm,iterations)
383 IMPLICIT NONE
384 !------------------------------------------------------------------------!
385 REAL, INTENT(IN) :: x1,x2
386 REAL, INTENT(OUT) :: root
387 INTEGER, INTENT(OUT) :: error
388 REAL, DIMENSION(:), INTENT(IN), OPTIONAL :: plist
389 REAL, INTENT(IN), OPTIONAL :: xm
390 INTEGER, INTENT(OUT), OPTIONAL :: iterations
391 !------------------------------------------------------------------------!
392 INTERFACE
393 PURE SUBROUTINE func(x,fx,plist)
394 IMPLICIT NONE
395 REAL, INTENT(IN) :: x
396 REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
397 REAL, INTENT(OUT) :: fx
398 END SUBROUTINE func
399 END INTERFACE
400 !------------------------------------------------------------------------!
401 TYPE(roots_typ) :: this
402 !------------------------------------------------------------------------!
403 CALL getroot_generic(this,step_king,func,x1,x2, &
404 default_accuracy,default_max_iterations,plist,xm)
405
406 root = this%root
407 error = this%error
408 IF (PRESENT(iterations)) iterations = this%iter
409
410 END SUBROUTINE getroot_king
411
412
413#ifndef DEBUG_OUTPUT
414 PURE &
415#endif
416 SUBROUTINE getroot_andersonbjoerk(func,x1,x2,root,error,plist,xm,iterations)
417 IMPLICIT NONE
418 !------------------------------------------------------------------------!
419 REAL, INTENT(IN) :: x1,x2
420 REAL, INTENT(OUT) :: root
421 INTEGER, INTENT(OUT) :: error
422 REAL, DIMENSION(:), INTENT(IN), OPTIONAL :: plist
423 REAL, INTENT(IN), OPTIONAL :: xm
424 INTEGER, INTENT(OUT), OPTIONAL :: iterations
425 !------------------------------------------------------------------------!
426 INTERFACE
427 PURE SUBROUTINE func(x,fx,plist)
428 IMPLICIT NONE
429 REAL, INTENT(IN) :: x
430 REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
431 REAL, INTENT(OUT) :: fx
432 END SUBROUTINE func
433 END INTERFACE
434 !------------------------------------------------------------------------!
435 TYPE(roots_typ) :: this
436 !------------------------------------------------------------------------!
437 CALL getroot_generic(this,step_andersonbjoerk,func,x1,x2, &
438 default_accuracy,default_max_iterations,plist,xm)
439
440 root = this%root
441 error = this%error
442 IF (PRESENT(iterations)) iterations = this%iter
443
444 END SUBROUTINE getroot_andersonbjoerk
445
446
447#ifndef DEBUG_OUTPUT
448 PURE &
449#endif
450 SUBROUTINE getroot_ridder(func,x1,x2,root,error,plist,xm,iterations)
451 IMPLICIT NONE
452 !------------------------------------------------------------------------!
453 REAL, INTENT(IN) :: x1,x2
454 REAL, INTENT(OUT) :: root
455 INTEGER, INTENT(OUT) :: error
456 REAL, DIMENSION(:), INTENT(IN), OPTIONAL :: plist
457 REAL, INTENT(IN), OPTIONAL :: xm
458 INTEGER, INTENT(OUT), OPTIONAL :: iterations
459 !------------------------------------------------------------------------!
460 INTERFACE
461 PURE SUBROUTINE func(x,fx,plist)
462 IMPLICIT NONE
463 REAL, INTENT(IN) :: x
464 REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
465 REAL, INTENT(OUT) :: fx
466 END SUBROUTINE func
467 END INTERFACE
468 !------------------------------------------------------------------------!
469 TYPE(roots_typ) :: this
470 !------------------------------------------------------------------------!
471 CALL getroot_generic(this,step_ridder,func,x1,x2, &
472 default_accuracy,default_max_iterations,plist,xm)
473
474 root = this%root
475 error = this%error
476 IF (PRESENT(iterations)) iterations = this%iter
477
478 END SUBROUTINE getroot_ridder
479
480
481#ifndef DEBUG_OUTPUT
482 PURE &
483#endif
484 SUBROUTINE getroot_brentdekker(func,x1,x2,root,error,plist,xm,iterations)
485 IMPLICIT NONE
486 !------------------------------------------------------------------------!
487 REAL, INTENT(IN) :: x1,x2
488 REAL, INTENT(OUT) :: root
489 INTEGER, INTENT(OUT) :: error
490 REAL, DIMENSION(:), INTENT(IN), OPTIONAL :: plist
491 REAL, INTENT(IN), OPTIONAL :: xm
492 INTEGER, INTENT(OUT), OPTIONAL :: iterations
493 !------------------------------------------------------------------------!
494 INTERFACE
495 PURE SUBROUTINE func(x,fx,plist)
496 IMPLICIT NONE
497 REAL, INTENT(IN) :: x
498 REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
499 REAL, INTENT(OUT) :: fx
500 END SUBROUTINE func
501 END INTERFACE
502 !------------------------------------------------------------------------!
503 TYPE(roots_typ) :: this
504 !------------------------------------------------------------------------!
505 CALL getroot_generic(this,step_brentdekker,func,x1,x2, &
506 default_accuracy,default_max_iterations,plist,xm)
507
508 root = this%root
509 error = this%error
510 IF (PRESENT(iterations)) iterations = this%iter
511
512 END SUBROUTINE getroot_brentdekker
513
514
515#ifndef DEBUG_OUTPUT
516 PURE &
517#endif
518 SUBROUTINE getroot_bisection(func,x1,x2,root,error,plist,xm,iterations)
519 IMPLICIT NONE
520 !------------------------------------------------------------------------!
521 REAL, INTENT(IN) :: x1,x2
522 REAL, INTENT(OUT) :: root
523 INTEGER, INTENT(OUT) :: error
524 REAL, DIMENSION(:), INTENT(IN), OPTIONAL :: plist
525 REAL, INTENT(IN), OPTIONAL :: xm
526 INTEGER, INTENT(OUT), OPTIONAL :: iterations
527 !------------------------------------------------------------------------!
528 INTERFACE
529 PURE SUBROUTINE func(x,fx,plist)
530 IMPLICIT NONE
531 REAL, INTENT(IN) :: x
532 REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
533 REAL, INTENT(OUT) :: fx
534 END SUBROUTINE func
535 END INTERFACE
536 !------------------------------------------------------------------------!
537 TYPE(roots_typ) :: this
538 !------------------------------------------------------------------------!
539 CALL getroot_generic(this,step_bisection,func,x1,x2, &
540 default_accuracy,default_max_iterations,plist,xm)
541
542 root = this%root
543 error = this%error
544 IF (PRESENT(iterations)) iterations = this%iter
545
546 END SUBROUTINE getroot_bisection
547
548
549 PURE SUBROUTINE updatebounds_bisection(this)
550 IMPLICIT NONE
551 !------------------------------------------------------------------------!
552 TYPE(roots_typ), INTENT(INOUT) :: this
553 !------------------------------------------------------------------------!
554 ! check sign
555 IF (this%f(1)*this%f(3).GT.0.0) THEN
556 ! sign(f1) == sign(f3) => x3 < root < x2
557 this%x(1)=this%x(3)
558 this%f(1)=this%f(3)
559 ELSE
560 ! sign(f1) != sign(f3) => x1 < root < x3
561 this%x(2)=this%x(3)
562 this%f(2)=this%f(3)
563 END IF
564 END SUBROUTINE updatebounds_bisection
565
566
567 PURE SUBROUTINE step_bisection(this)
568 IMPLICIT NONE
569 !------------------------------------------------------------------------!
570 TYPE(roots_typ), INTENT(INOUT) :: this
571 !------------------------------------------------------------------------!
572 CALL updatebounds_bisection(this)
573 CALL step_arithmeticmean(this)
574 END SUBROUTINE step_bisection
575
576
577 PURE SUBROUTINE step_newton(this)
578 IMPLICIT NONE
579 !------------------------------------------------------------------------!
580 TYPE(roots_typ), INTENT(INOUT) :: this
581 !------------------------------------------------------------------------!
582 ! determine new brackets for the root
583!CDIR IEXPAND
584 CALL updatebounds_bisection(this)
585 this%dxold = this%dx
586 ! check if we are out of bounds
587 ! or if the convergence is to slow
588 IF ( ( (this%df*(this%x(3)-this%x(2))-this%f(3).GE.0.0) &
589 .OR. (this%df*(this%x(3)-this%x(1))-this%f(3).LE.0.0) ) &
590 .OR. (abs(2*this%f(3)).GT.abs(this%dxold*this%df)) ) THEN
591 ! compute new estimate for the route using arithmetic mean value
592!CDIR IEXPAND
593 CALL step_arithmeticmean(this)
594 ELSE
595 ! compute new estimate for the root using Newton iteration
596 this%dx = this%f(3) / this%df
597 this%x(3) = this%x(3) - this%dx
598 this%dx = abs(this%dx)
599 END IF
600 END SUBROUTINE step_newton
601
602
603 PURE SUBROUTINE step_regulafalsi(this)
604 IMPLICIT NONE
605 !------------------------------------------------------------------------!
606 TYPE(roots_typ), INTENT(INOUT) :: this
607 !------------------------------------------------------------------------!
608 ! calculate new interval boundaries
609!CDIR IEXPAND
610 CALL updatebounds_bisection(this)
611 ! determine new approximation for the root
612!CDIR IEXPAND
613 CALL step_secantmethod(this)
614 END SUBROUTINE step_regulafalsi
615
616
617 PURE SUBROUTINE step_pegasus(this)
618 IMPLICIT NONE
619 !------------------------------------------------------------------------!
620 TYPE(roots_typ), INTENT(INOUT) :: this
621 !------------------------------------------------------------------------!
622 REAL :: g
623 !------------------------------------------------------------------------!
624 ! calculate new interval boundaries
625 IF (this%f(2)*this%f(3).LT.0.0) THEN
626 this%x(1) = this%x(2)
627 this%x(2) = this%x(3)
628 this%f(1) = this%f(2)
629 this%f(2) = this%f(3)
630 ELSE
631 g = this%f(2)/(this%f(3)+this%f(2))
632 this%x(2) = this%x(3)
633 this%f(1) = g*this%f(1)
634 this%f(2) = this%f(3)
635 END IF
636 ! determine new approximation for the root
637!CDIR IEXPAND
638 CALL step_secantmethod(this)
639 END SUBROUTINE step_pegasus
640
641
642 PURE SUBROUTINE step_king(this)
643 IMPLICIT NONE
644 !------------------------------------------------------------------------!
645 TYPE(roots_typ), INTENT(INOUT) :: this
646 REAL :: g
647 !------------------------------------------------------------------------!
648 IF (this%do_modified_step) THEN
649 IF (this%f(2)*this%f(3).LT.0.0) THEN
650 ! interchange (x1,f1) and (x2,f2)
651 g = this%x(1)
652 this%x(1) = this%x(2)
653 this%x(2) = g
654 g = this%f(1)
655 this%f(1) = this%f(2)
656 this%f(2) = g
657 END IF
658 ! Pegasus step
659 IF (this%f(2)*this%f(3).GT.0.0) THEN
660 g = this%f(2)/(this%f(3)+this%f(2))
661 this%x(2) = this%x(3)
662 this%f(1) = g*this%f(1)
663 this%f(2) = this%f(3)
664 END IF
665 this%do_modified_step = .false.
666 ELSE
667 IF (this%f(2)*this%f(3).LT.0.0) THEN
668 this%x(1) = this%x(2)
669 this%f(1) = this%f(2)
670 this%x(2) = this%x(3)
671 this%f(2) = this%f(3)
672 this%do_modified_step = .true.
673 ELSE
674 ! Pegasus step
675 IF (this%f(2)*this%f(3).GT.0.0) THEN
676 g = this%f(2)/(this%f(3)+this%f(2))
677 this%x(2) = this%x(3)
678 this%f(1) = g*this%f(1)
679 this%f(2) = this%f(3)
680 END IF
681 this%do_modified_step = .false.
682 END IF
683 END IF
684 ! determine new approximation for the root
685!CDIR IEXPAND
686 CALL step_secantmethod(this)
687 END SUBROUTINE step_king
688
689
690 PURE SUBROUTINE step_andersonbjoerk(this)
691 IMPLICIT NONE
692 !------------------------------------------------------------------------!
693 TYPE(roots_typ), INTENT(INOUT) :: this
694 !------------------------------------------------------------------------!
695 REAL :: g
696 !------------------------------------------------------------------------!
697 ! calculate new interval boundaries
698 IF (this%f(2)*this%f(3).LT.0.0) THEN
699 this%x(1) = this%x(2)
700 this%x(2) = this%x(3)
701 this%f(1) = this%f(2)
702 this%f(2) = this%f(3)
703 ELSE
704 g = 1.-this%f(3)/this%f(2)
705 IF(g.LE.0.) g = 0.5
706 this%x(2) = this%x(3)
707 this%f(1) = g*this%f(1)
708 this%f(2) = this%f(3)
709 END IF
710 ! determine new approximation for the root
711!CDIR IEXPAND
712 CALL step_secantmethod(this)
713 END SUBROUTINE step_andersonbjoerk
714
715
716 PURE SUBROUTINE step_ridder(this)
717 IMPLICIT NONE
718 !------------------------------------------------------------------------!
719 TYPE(roots_typ), INTENT(INOUT) :: this
720 !------------------------------------------------------------------------!
721 REAL :: denom,dx
722 !------------------------------------------------------------------------!
723 IF (this%do_modified_step) THEN
724 ! compute arithmetic mean of the new boundaries
725 CALL step_bisection(this)
726 ! next update uses Ridders 3 point formula
727 this%do_modified_step = .false.
728 ELSE
729 ! estimate the new step from the arithmetic mean value
730 denom = sqrt(this%f(3)*this%f(3) - this%f(1)*this%f(2)) + tiny(denom)
731 dx = (this%x(3)-this%x(1)) * sign(1.0,this%f(1)-this%f(2))*this%f(3) / denom
732 IF (dx.GT.0.0) THEN
733 this%x(1) = this%x(3)
734 this%f(1) = this%f(3)
735 ELSE
736 this%x(2) = this%x(3)
737 this%f(2) = this%f(3)
738 END IF
739 this%dx = 0.5*this%dx
740 this%x(3) = this%x(3) + dx
741 ! next step is a bisection step
742 this%do_modified_step = .true.
743 END IF
744 END SUBROUTINE step_ridder
745
746
747 PURE SUBROUTINE step_brentdekker(this)
748 IMPLICIT NONE
749 !------------------------------------------------------------------------!
750 TYPE(roots_typ), INTENT(INOUT) :: this
751 REAL :: m,p,q
752 !------------------------------------------------------------------------!
753 ! bracket the root
754 IF (this%f(3)*this%f(2).GT.0.0) THEN
755 ! x3, x2 on the same side of the root
756 this%x(2) = this%x(1)
757 this%f(2) = this%f(1)
758 this%d = this%x(1)-this%x(3)
759 this%e = this%d
760 END IF
761 ! root is now bracketed between x3 and x2
762
763 ! x3 should be closer to the root than x2, i.e. |f3|<|f2|
764 IF (abs(this%f(2)).LT.abs(this%f(3))) THEN
765 ! interchange (x3,f3) and (x2,f2) using (x1,f1) as temporary space
766 this%x(1) = this%x(3)
767 this%x(3) = this%x(2)
768 this%x(2) = this%x(1)
769 this%f(1) = this%f(3)
770 this%f(3) = this%f(2)
771 this%f(2) = this%f(1)
772 END IF
773
774 m = 0.5*(this%x(2)-this%x(3))
775
776 IF ((abs(this%e).LT.this%tol) .OR. abs(this%f(1)).LE.abs(this%f(3))) THEN
777 ! convergence too slow -> do bisection step
778 this%d = m
779 this%e = m
780 ELSE
781 ! try interpolation
782 CALL saveinversequadraticinterpolation(this%x(1),this%x(3),this%x(2), &
783 this%f(1),this%f(3),this%f(2),p,q)
784
785 IF (2*p.LT. min(3*m*q-abs(this%tol*q), abs(this%e*q))) THEN
786 ! interpolation worked store old correction d in e
787 ! and compute new correction d
788 this%e = this%d
789 this%d = p / q
790 ELSE
791 ! interpolation failed, do bisection
792 this%e = m
793 this%d = m
794 END IF
795 END IF
796
797 this%x(1) = this%x(3)
798 this%f(1) = this%f(3)
799
800 IF (abs(this%d).GT.this%tol) THEN
801 this%x(3) = this%x(3) + this%d
802 ELSE
803 this%x(3) = this%x(3) + sign(this%tol,m)
804 END IF
805
806 this%dx = abs(this%x(2)-this%x(3))
807 END SUBROUTINE step_brentdekker
808
809
810 PURE SUBROUTINE step_secantmethod(this)
811 IMPLICIT NONE
812 !------------------------------------------------------------------------!
813 TYPE(roots_typ), INTENT(INOUT) :: this
814 !------------------------------------------------------------------------!
815 INTEGER :: one,two
816 !------------------------------------------------------------------------!
817 IF (abs(this%f(1)).GT.abs(this%f(2))) THEN
818 ! x2 is probably the better approximation for the root
819 ! set indices as usual
820 one = 1
821 two = 2
822 ELSE
823 ! x1 is probably the better approximation for the root
824 ! interchange indices
825 one = 2
826 two = 1
827 END IF
828!CDIR IEXPAND
829 this%dx = linearinterpolation(this%x(one),this%x(two),this%f(one),this%f(two))
830 ! improves convergence
831 IF (abs(this%dx).LE.this%tol) THEN
832 this%dx = 0.8*this%tol*sign(1.0,this%x(one)-this%x(two))
833 END IF
834 ! new estimate for the root
835 this%x(3) = this%x(two) + this%dx
836 ! store this for error control
837 this%dx = abs(this%x(one)-this%x(two))
838 END SUBROUTINE step_secantmethod
839
840
841 PURE SUBROUTINE saveinversequadraticinterpolation(a,b,c,fa,fb,fc,db_numer,db_denom)
842 IMPLICIT NONE
843 !------------------------------------------------------------------------!
844 REAL, INTENT(IN) :: a,b,c,fa,fb,fc
845 REAL, INTENT(OUT) :: db_numer,db_denom
846 !------------------------------------------------------------------------!
847 REAL :: fbdivfa,fadivfc,fbdivfc,p,q
848 !------------------------------------------------------------------------!
849 fbdivfa = fb / fa
850 IF ((fa.EQ.fc).OR.(fb.EQ.fc)) THEN
851 ! fallback: linear interpolation
852 p = (c-b) * fbdivfa
853 q = 1.0 - fbdivfa
854 ELSE
855 ! try inverse quadratic interpolation
856 fadivfc = fa / fc
857 fbdivfc = fb / fc
858 p = fbdivfa * ((c-b)*fadivfc*(fadivfc-fbdivfc) - (b-a)*(fbdivfc-1.0) )
859 q = (fadivfc - 1.0) * (fbdivfc - 1.0) * (fbdivfa - 1.0)
860 END IF
861 ! move the sign of the ratio p/q to the denominator q;
862 ! thus the return value of the numerator is always > 0
863 IF (p.GT.0.0) THEN
864 db_numer = p
865 db_denom = -q
866 ELSE
867 db_numer = -p
868 db_denom = q
869 END IF
870 END SUBROUTINE saveinversequadraticinterpolation
871
872
873 PURE SUBROUTINE step_arithmeticmean(this)
874 IMPLICIT NONE
875 !------------------------------------------------------------------------!
876 TYPE(roots_typ), INTENT(INOUT) :: this
877 !------------------------------------------------------------------------!
878 ! new estimate for the root
879!CDIR IEXPAND
880 this%x(3) = arithmeticmean(this%x(1),this%x(2))
881 ! store this for error control
882 this%dx = abs(this%x(2)-this%x(1))
883 END SUBROUTINE step_arithmeticmean
884
885
886 PURE FUNCTION arithmeticmean(a,b) RESULT(s)
887 IMPLICIT NONE
888 !------------------------------------------------------------------------!
889 REAL, INTENT(IN) :: a,b
890 REAL :: s
891 !------------------------------------------------------------------------!
892 s = 0.5*(a+b)
893 END FUNCTION arithmeticmean
894
895
896 PURE FUNCTION linearinterpolation(x1,x2,f1,f2) RESULT(dx)
897 IMPLICIT NONE
898 !------------------------------------------------------------------------!
899 REAL, INTENT(IN) :: x1,x2,f1,f2
900 REAL :: dx
901 !------------------------------------------------------------------------!
902 dx = (x1-x2)*f2 / (f2 - f1 + tiny(f1))
903 END FUNCTION linearinterpolation
904
905
906 PURE SUBROUTINE testconvergence(this)
907 IMPLICIT NONE
908 !------------------------------------------------------------------------!
909 TYPE(roots_typ), INTENT(INOUT) :: this
910 !------------------------------------------------------------------------!
911 IF (this%error.NE.0) THEN
912 ! error occurred => abort iteration
913 this%iterate = .false.
914 ELSE
915 ! check convergence
916 this%tol = this%eps*abs(this%x(3))
917 IF (abs(this%f(3)).LE.tiny(this%f(3)) .OR. this%dx.LE.this%tol) THEN
918 ! converged => store current estimate xm in root and finish iteration
919 this%iterate = .false.
920 this%root = this%x(3)
921 ELSE
922 ! not converged => check iteration limit
923 IF (this%iter.GE.this%MAX_ITERATIONS) THEN
924 ! exceeded => store error code and abort iteration
925 this%error = 2
926 this%iterate = .false.
927 ELSE
928 ! continue iteration
929 this%iter = this%iter + 1
930 this%iterate = .true.
931 END IF
932 END IF
933 END IF
934 END SUBROUTINE testconvergence
935
936END MODULE roots
funcd
pure subroutine funcd(y, fy, dfy, plist)
Definition: bondi2d.f90:416
roots
root finding subroutines
Definition: roots.f90:45
roots::getroot_brentdekker
subroutine, public getroot_brentdekker(func, x1, x2, root, error, plist, xm, iterations)
Definition: roots.f90:485
roots::step_regulafalsi
pure subroutine step_regulafalsi(this)
Definition: roots.f90:604
roots::step_secantmethod
pure subroutine step_secantmethod(this)
Definition: roots.f90:811
roots::testconvergence
pure subroutine testconvergence(this)
Definition: roots.f90:907
roots::step_newton
pure subroutine step_newton(this)
Definition: roots.f90:578
roots::getroot_bisection
subroutine, public getroot_bisection(func, x1, x2, root, error, plist, xm, iterations)
Definition: roots.f90:519
roots::step_andersonbjoerk
pure subroutine step_andersonbjoerk(this)
Definition: roots.f90:691
roots::getroot_regulafalsi
subroutine, public getroot_regulafalsi(func, x1, x2, root, error, plist, xm, iterations)
Definition: roots.f90:315
roots::getroot_king
subroutine, public getroot_king(func, x1, x2, root, error, plist, xm, iterations)
Definition: roots.f90:383
roots::getroot_ridder
subroutine, public getroot_ridder(func, x1, x2, root, error, plist, xm, iterations)
Definition: roots.f90:451
roots::error_message
character(len=64), dimension(0:4), parameter error_message
Definition: roots.f90:68
roots::step_brentdekker
pure subroutine step_brentdekker(this)
Definition: roots.f90:748
roots::linearinterpolation
pure real function linearinterpolation(x1, x2, f1, f2)
Definition: roots.f90:897
roots::getroot_andersonbjoerk
subroutine, public getroot_andersonbjoerk(func, x1, x2, root, error, plist, xm, iterations)
Definition: roots.f90:417
roots::step_ridder
pure subroutine step_ridder(this)
Definition: roots.f90:717
roots::getroot_generic
subroutine getroot_generic(this, Stepper, func, x1, x2, dxacc, maxiter, plist, xm)
Definition: roots.f90:178
roots::getroot_newton
subroutine, public getroot_newton(funcd, x1, x2, root, error, plist, xm, iterations)
Definition: roots.f90:247
roots::arithmeticmean
pure real function arithmeticmean(a, b)
Definition: roots.f90:887
roots::default_max_iterations
integer, parameter, public default_max_iterations
Definition: roots.f90:66
roots::saveinversequadraticinterpolation
pure subroutine saveinversequadraticinterpolation(a, b, c, fa, fb, fc, db_numer, db_denom)
Definition: roots.f90:842
roots::initroots
pure subroutine, public initroots(this, x1, x2, f1, f2, dxacc, maxiter)
Definition: roots.f90:98
roots::step_bisection
pure subroutine step_bisection(this)
Definition: roots.f90:568
roots::getroot_pegasus
subroutine, public getroot_pegasus(func, x1, x2, root, error, plist, xm, iterations)
Definition: roots.f90:349
roots::geterrormessage
character(len=64) function, public geterrormessage(error)
Definition: roots.f90:160
roots::step_king
pure subroutine step_king(this)
Definition: roots.f90:643
roots::step_pegasus
pure subroutine step_pegasus(this)
Definition: roots.f90:618
roots::default_accuracy
real, parameter, public default_accuracy
Definition: roots.f90:67
roots::updatebounds_bisection
pure subroutine updatebounds_bisection(this)
Definition: roots.f90:550
roots::step_arithmeticmean
pure subroutine step_arithmeticmean(this)
Definition: roots.f90:874
func
pure subroutine func(x, fx, plist)
Definition: rootstest.f90:192
roots::roots_typ
basic type for root finding functions
Definition: roots.f90:50
Generated by   doxygen 1.9.4