51 REAL :: eps,tol,root,dx,dxold,df,d,e
53 INTEGER :: iter,max_iterations,error
54 LOGICAL :: iterate,do_modified_step
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 " /)
97 PURE SUBROUTINE initroots(this,x1,x2,f1,f2,dxacc,maxiter)
101 REAL,
INTENT(IN) :: x1,x2
102 REAL,
INTENT(IN) :: f1,f2
103 REAL,
OPTIONAL,
INTENT(IN) :: dxacc
104 INTEGER,
OPTIONAL,
INTENT(IN) :: maxiter
107 IF (f1*f2.GT.0.0)
THEN
123 this%dx = abs(this%x(2)-this%x(1))
126 IF(
PRESENT(dxacc))
THEN
127 IF (dxacc.GT.epsilon(dxacc))
THEN
136 IF(
PRESENT(maxiter))
THEN
137 IF (maxiter.GT.1)
THEN
138 this%max_iterations = maxiter
147 this%d = this%x(2)-this%x(1)
152 this%iterate = .true.
153 this%do_modified_step = .true.
162 INTEGER,
INTENT(IN) :: error
163 CHARACTER(LEN=64) :: msg
182 REAL,
OPTIONAL :: dxacc, plist(:), xm
183 INTEGER,
OPTIONAL :: maxiter
185 INTENT(IN) :: x1,x2,dxacc,maxiter,xm,plist
186 INTENT(INOUT) :: this
189 PURE SUBROUTINE stepper(root)
193 END SUBROUTINE stepper
196 PURE SUBROUTINE func(x,fx,plist)
198 REAL,
INTENT(IN) :: x
199 REAL,
INTENT(IN),
DIMENSION(:),
OPTIONAL :: plist
200 REAL,
INTENT(OUT) :: fx
207 CALL func(x1,f1,plist)
208 CALL func(x2,f2,plist)
213 IF (this%error.EQ.0)
THEN
215 IF (
PRESENT(xm))
THEN
221 print
'(A3,A18,A10,A18,A18,2(A13,A10))', &
222 "#",
"x3",
"f3",
"dx",
"tol",
"x1",
"f1",
"x2",
"f2"
227 CALL func(this%x(3),this%f(3),plist)
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)
235 IF (.NOT.this%iterate)
EXIT
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
257 PURE SUBROUTINE funcd(x,fx,dfx,plist)
259 REAL,
INTENT(IN) :: x
260 REAL,
INTENT(IN),
DIMENSION(:),
OPTIONAL :: plist
261 REAL,
INTENT(OUT) :: fx,dfx
266 REAL :: f1,f2,df1,df2
269 CALL funcd(x1,f1,df1,plist)
270 CALL funcd(x2,f2,df2,plist)
275 IF (this%error.EQ.0)
THEN
277 IF (
PRESENT(xm))
THEN
283 print
'(A3,A18,A10,A18,A18,2(A13,A10))', &
284 "#",
"x3",
"f3",
"dx",
"tol",
"x1",
"f1",
"x2",
"f2"
289 CALL funcd(this%x(3),this%f(3),this%df,plist)
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)
297 IF (.NOT.this%iterate)
EXIT
306 IF (
PRESENT(iterations)) iterations = this%iter
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
325 PURE SUBROUTINE func(x,fx,plist)
327 REAL,
INTENT(IN) :: x
328 REAL,
INTENT(IN),
DIMENSION(:),
OPTIONAL :: plist
329 REAL,
INTENT(OUT) :: fx
340 IF (
PRESENT(iterations)) iterations = this%iter
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
359 PURE SUBROUTINE func(x,fx,plist)
361 REAL,
INTENT(IN) :: x
362 REAL,
INTENT(IN),
DIMENSION(:),
OPTIONAL :: plist
363 REAL,
INTENT(OUT) :: fx
374 IF (
PRESENT(iterations)) iterations = this%iter
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
393 PURE SUBROUTINE func(x,fx,plist)
395 REAL,
INTENT(IN) :: x
396 REAL,
INTENT(IN),
DIMENSION(:),
OPTIONAL :: plist
397 REAL,
INTENT(OUT) :: fx
408 IF (
PRESENT(iterations)) iterations = this%iter
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
427 PURE SUBROUTINE func(x,fx,plist)
429 REAL,
INTENT(IN) :: x
430 REAL,
INTENT(IN),
DIMENSION(:),
OPTIONAL :: plist
431 REAL,
INTENT(OUT) :: fx
442 IF (
PRESENT(iterations)) iterations = this%iter
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
461 PURE SUBROUTINE func(x,fx,plist)
463 REAL,
INTENT(IN) :: x
464 REAL,
INTENT(IN),
DIMENSION(:),
OPTIONAL :: plist
465 REAL,
INTENT(OUT) :: fx
476 IF (
PRESENT(iterations)) iterations = this%iter
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
495 PURE SUBROUTINE func(x,fx,plist)
497 REAL,
INTENT(IN) :: x
498 REAL,
INTENT(IN),
DIMENSION(:),
OPTIONAL :: plist
499 REAL,
INTENT(OUT) :: fx
510 IF (
PRESENT(iterations)) iterations = this%iter
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
529 PURE SUBROUTINE func(x,fx,plist)
531 REAL,
INTENT(IN) :: x
532 REAL,
INTENT(IN),
DIMENSION(:),
OPTIONAL :: plist
533 REAL,
INTENT(OUT) :: fx
544 IF (
PRESENT(iterations)) iterations = this%iter
555 IF (this%f(1)*this%f(3).GT.0.0)
THEN
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
596 this%dx = this%f(3) / this%df
597 this%x(3) = this%x(3) - this%dx
598 this%dx = abs(this%dx)
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)
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)
648 IF (this%do_modified_step)
THEN
649 IF (this%f(2)*this%f(3).LT.0.0)
THEN
652 this%x(1) = this%x(2)
655 this%f(1) = this%f(2)
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)
665 this%do_modified_step = .false.
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.
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)
681 this%do_modified_step = .false.
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)
704 g = 1.-this%f(3)/this%f(2)
706 this%x(2) = this%x(3)
707 this%f(1) = g*this%f(1)
708 this%f(2) = this%f(3)
723 IF (this%do_modified_step)
THEN
727 this%do_modified_step = .false.
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
733 this%x(1) = this%x(3)
734 this%f(1) = this%f(3)
736 this%x(2) = this%x(3)
737 this%f(2) = this%f(3)
739 this%dx = 0.5*this%dx
740 this%x(3) = this%x(3) + dx
742 this%do_modified_step = .true.
754 IF (this%f(3)*this%f(2).GT.0.0)
THEN
756 this%x(2) = this%x(1)
757 this%f(2) = this%f(1)
758 this%d = this%x(1)-this%x(3)
764 IF (abs(this%f(2)).LT.abs(this%f(3)))
THEN
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)
774 m = 0.5*(this%x(2)-this%x(3))
776 IF ((abs(this%e).LT.this%tol) .OR. abs(this%f(1)).LE.abs(this%f(3)))
THEN
783 this%f(1),this%f(3),this%f(2),p,q)
785 IF (2*p.LT. min(3*m*q-abs(this%tol*q), abs(this%e*q)))
THEN
797 this%x(1) = this%x(3)
798 this%f(1) = this%f(3)
800 IF (abs(this%d).GT.this%tol)
THEN
801 this%x(3) = this%x(3) + this%d
803 this%x(3) = this%x(3) + sign(this%tol,m)
806 this%dx = abs(this%x(2)-this%x(3))
817 IF (abs(this%f(1)).GT.abs(this%f(2)))
THEN
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))
835 this%x(3) = this%x(two) + this%dx
837 this%dx = abs(this%x(one)-this%x(two))
844 REAL,
INTENT(IN) :: a,b,c,fa,fb,fc
845 REAL,
INTENT(OUT) :: db_numer,db_denom
847 REAL :: fbdivfa,fadivfc,fbdivfc,p,q
850 IF ((fa.EQ.fc).OR.(fb.EQ.fc))
THEN
858 p = fbdivfa * ((c-b)*fadivfc*(fadivfc-fbdivfc) - (b-a)*(fbdivfc-1.0) )
859 q = (fadivfc - 1.0) * (fbdivfc - 1.0) * (fbdivfa - 1.0)
882 this%dx = abs(this%x(2)-this%x(1))
889 REAL,
INTENT(IN) :: a,b
899 REAL,
INTENT(IN) :: x1,x2,f1,f2
902 dx = (x1-x2)*f2 / (f2 - f1 + tiny(f1))
911 IF (this%error.NE.0)
THEN
913 this%iterate = .false.
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
919 this%iterate = .false.
920 this%root = this%x(3)
923 IF (this%iter.GE.this%MAX_ITERATIONS)
THEN
926 this%iterate = .false.
929 this%iter = this%iter + 1
930 this%iterate = .true.
pure subroutine funcd(y, fy, dfy, plist)
subroutine, public getroot_brentdekker(func, x1, x2, root, error, plist, xm, iterations)
pure subroutine step_regulafalsi(this)
pure subroutine step_secantmethod(this)
pure subroutine testconvergence(this)
pure subroutine step_newton(this)
subroutine, public getroot_bisection(func, x1, x2, root, error, plist, xm, iterations)
pure subroutine step_andersonbjoerk(this)
subroutine, public getroot_regulafalsi(func, x1, x2, root, error, plist, xm, iterations)
subroutine, public getroot_king(func, x1, x2, root, error, plist, xm, iterations)
subroutine, public getroot_ridder(func, x1, x2, root, error, plist, xm, iterations)
character(len=64), dimension(0:4), parameter error_message
pure subroutine step_brentdekker(this)
pure real function linearinterpolation(x1, x2, f1, f2)
subroutine, public getroot_andersonbjoerk(func, x1, x2, root, error, plist, xm, iterations)
pure subroutine step_ridder(this)
subroutine getroot_generic(this, Stepper, func, x1, x2, dxacc, maxiter, plist, xm)
subroutine, public getroot_newton(funcd, x1, x2, root, error, plist, xm, iterations)
pure real function arithmeticmean(a, b)
integer, parameter, public default_max_iterations
pure subroutine saveinversequadraticinterpolation(a, b, c, fa, fb, fc, db_numer, db_denom)
pure subroutine, public initroots(this, x1, x2, f1, f2, dxacc, maxiter)
pure subroutine step_bisection(this)
subroutine, public getroot_pegasus(func, x1, x2, root, error, plist, xm, iterations)
character(len=64) function, public geterrormessage(error)
pure subroutine step_king(this)
pure subroutine step_pegasus(this)
real, parameter, public default_accuracy
pure subroutine updatebounds_bisection(this)
pure subroutine step_arithmeticmean(this)
pure subroutine func(x, fx, plist)
basic type for root finding functions