@@ -28,10 +28,10 @@ program test_hybrd
28
28
1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ]
29
29
30
30
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
33
33
real (wp) :: factor, fnorm1, fnorm2
34
- real (wp),allocatable :: fvec(:) , wa(:) , x(:)
34
+ real (wp),dimension (:), allocatable :: fvec, wa, x
35
35
36
36
integer , parameter :: nwrite = output_unit ! ! logical output unit
37
37
real (wp), parameter :: one = 1.0_wp
@@ -84,24 +84,50 @@ program test_hybrd
84
84
' EXIT PARAMETER' , info, &
85
85
' FINAL APPROXIMATE SOLUTION' , x(1 :n)
86
86
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)
98
88
end do
99
89
end if
100
90
end do
101
91
102
92
contains
103
93
! *****************************************************************************************
104
94
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
+
105
131
! *****************************************************************************************
106
132
! >
107
133
! The calling sequence of fcn should be identical to the
@@ -388,14 +414,20 @@ subroutine vecfcn(n, x, Fvec, Nprob)
388
414
Fvec(4 ) = c6* temp2 + c4* (x(4 ) - one) + c5* (x(2 ) - one)
389
415
case (5 )
390
416
! HELICAL VALLEY FUNCTION.
417
+ write (* ,* ) ' 1'
391
418
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
395
426
temp2 = sqrt (x(1 )** 2 + x(2 )** 2 )
396
427
Fvec(1 ) = ten* (x(3 ) - ten* temp1)
397
428
Fvec(2 ) = ten* (temp2 - one)
398
429
Fvec(3 ) = x(3 )
430
+ write (* ,* ) ' 2'
399
431
case (6 )
400
432
! WATSON FUNCTION.
401
433
do k = 1 , n
0 commit comments