@@ -28,10 +28,10 @@ program test_hybrd
2828 1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ]
2929
3030 integer :: i, ic, info, k, lwa, n, NFEv, NPRob, ntries, icase
31- integer :: na( 55 ) , nf( 55 ) , np( 55 ) , nx( 55 )
32- real (wp) :: fnm( 55 )
31+ integer , dimension ( 55 ) :: na, nf, np, nx
32+ real (wp), dimension ( 55 ) :: fnm
3333 real (wp) :: factor, fnorm1, fnorm2
34- real (wp),allocatable :: fvec(:) , wa(:) , x(:)
34+ real (wp),dimension (:), allocatable :: fvec, wa, x
3535
3636 integer , parameter :: nwrite = output_unit ! ! logical output unit
3737 real (wp), parameter :: one = 1.0_wp
@@ -84,24 +84,50 @@ program test_hybrd
8484 ' EXIT PARAMETER' , info, &
8585 ' FINAL APPROXIMATE SOLUTION' , x(1 :n)
8686 factor = ten* factor
87-
88- ! compare with previously generated solutions:
89- if (info_original(ic)<5 .and. & ! ignore any where the original minpack failed
90- any (abs ( solution(ic) - x)>tol .and. &
91- abs ((solution(ic) - x)/ (solution(ic))) > solution_reltol)) then
92- write (nwrite,' (A)' ) ' Failed case'
93- write (nwrite, ' (//5x, a//(5x, 5d15.7))' ) ' Expected x: ' , solution(ic)
94- write (nwrite, ' (/5x, a//(5x, 5d15.7))' ) ' Computed x: ' , x
95- error stop
96- end if
97-
87+ call compare_solutions(ic, x, solution_reltol, tol)
9888 end do
9989 end if
10090 end do
10191
10292 contains
10393! *****************************************************************************************
10494
95+ ! *****************************************************************************************
96+ ! >
97+ ! Compare with previously generated solutions.
98+
99+ subroutine compare_solutions (ic , x , reltol , abstol )
100+
101+ implicit none
102+
103+ integer ,intent (in ) :: ic ! ! problem number (index is `solution` vector)
104+ real (wp),dimension (:),intent (in ) :: x ! ! computed `x` vector from the method
105+ real (wp),intent (in ) :: reltol ! ! relative tolerance for `x` to pass
106+ real (wp),intent (in ) :: abstol ! ! absolute tolerance for `x` to pass
107+
108+ real (wp),dimension (size (x)) :: diff, absdiff, reldiff
109+
110+ if (info_original(ic)<5 ) then ! ignore any where the original minpack failed
111+ diff = solution(ic) - x
112+ absdiff = abs (diff)
113+ if (any (absdiff> abstol)) then ! first do an absolute diff
114+ ! also do a rel diff if the abs diff fails (also protect for divide by zero)
115+ reldiff = absdiff
116+ where (solution(ic) /= 0.0_wp ) reldiff = absdiff / abs (solution(ic))
117+ if (any (reldiff > reltol)) then
118+ write (nwrite,' (A)' ) ' Failed case'
119+ write (nwrite, ' (//5x, a//(5x, 5d15.7))' ) ' Expected x: ' , solution(ic)
120+ write (nwrite, ' (/5x, a//(5x, 5d15.7))' ) ' Computed x: ' , x
121+ write (nwrite, ' (/5x, a//(5x, 5d15.7))' ) ' absdiff: ' , absdiff
122+ write (nwrite, ' (/5x, a//(5x, 5d15.7))' ) ' reldiff: ' , reldiff
123+ error stop ! test failed
124+ end if
125+ end if
126+ end if
127+
128+ end subroutine compare_solutions
129+ ! *****************************************************************************************
130+
105131! *****************************************************************************************
106132! >
107133! The calling sequence of fcn should be identical to the
@@ -388,14 +414,20 @@ subroutine vecfcn(n, x, Fvec, Nprob)
388414 Fvec(4 ) = c6* temp2 + c4* (x(4 ) - one) + c5* (x(2 ) - one)
389415 case (5 )
390416 ! HELICAL VALLEY FUNCTION.
417+ write (* ,* ) ' 1'
391418 tpi = eight* atan (one)
392- temp1 = sign (c7, x(2 ))
393- if (x(1 ) > zero) temp1 = atan (x(2 )/ x(1 ))/ tpi
394- if (x(1 ) < zero) temp1 = atan (x(2 )/ x(1 ))/ tpi + c8
419+ if (x(1 ) > zero) then
420+ temp1 = atan (x(2 )/ x(1 ))/ tpi
421+ else if (x(1 ) < zero) then
422+ temp1 = atan (x(2 )/ x(1 ))/ tpi + c8
423+ else
424+ temp1 = sign (c7, x(2 )) ! does this ever happen?
425+ end if
395426 temp2 = sqrt (x(1 )** 2 + x(2 )** 2 )
396427 Fvec(1 ) = ten* (x(3 ) - ten* temp1)
397428 Fvec(2 ) = ten* (temp2 - one)
398429 Fvec(3 ) = x(3 )
430+ write (* ,* ) ' 2'
399431 case (6 )
400432 ! WATSON FUNCTION.
401433 do k = 1 , n
0 commit comments