GCC Code Coverage Report

Directory: src/solver/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 96.8% 60 / 0 / 62
Functions: 100.0% 2 / 0 / 2
Branches: 80.0% 40 / 0 / 50
findroot.c
Line Branch Exec Source
1 //-----------------------------------------------------------------------------
2 // findroot.c
3 //
4 // Finds solution of func(x) = 0 using either the Newton-Raphson
5 // method or Ridder's Method.
6 // Based on code from Numerical Recipes in C (Cambridge University
7 // Press, 1992).
8 //
9 // Date: 11/19/13
10 // Author: L. Rossman
11 //-----------------------------------------------------------------------------
12
13 #include <math.h>
14 #include "findroot.h"
15
16 #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
17 #define MAXIT 60
18
19
20 1833978 int findroot_Newton(double x1, double x2, double* rts, double xacc,
21 void (*func) (double x, double* f, double* df, void* p),
22 void* p)
23 //
24 // Using a combination of Newton-Raphson and bisection, find the root of a
25 // function func bracketed between x1 and x2. The root, returned in rts,
26 // will be refined until its accuracy is known within +/-xacc. func is a
27 // user-supplied routine, that returns both the function value and the first
28 // derivative of the function. p is a pointer to any auxilary data structure
29 // that func may require. It can be NULL if not needed. The function returns
30 // the number of function evaluations used or 0 if the maximum allowed
31 // iterations were exceeded.
32 //
33 // NOTES:
34 // 1. The calling program must insure that the signs of func(x1) and func(x2)
35 // are not the same, otherwise x1 and x2 do not bracket the root.
36 // 2. If func(x1) > func(x2) then the order of x1 and x2 should be
37 // switched in the call to Newton.
38 //
39 {
40 1833978 int j, n = 0;
41 double df, dx, dxold, f, x;
42 double temp, xhi, xlo;
43
44 // Initialize the "stepsize before last" and the last step.
45 1833978 x = *rts;
46 1833978 xlo = x1;
47 1833978 xhi = x2;
48 1833978 dxold = fabs(x2-x1);
49 1833978 dx = dxold;
50 1833978 func(x, &f, &df, p);
51 1833978 n++;
52
53 // Loop over allowed iterations.
54
1/2
✓ Branch 0 taken 8882060 times.
✗ Branch 1 not taken.
 8882060 for (j=1; j<=MAXIT; j++)
55 {
56 // Bisect if Newton out of range or not decreasing fast enough.
57
2/2
✓ Branch 0 taken 7849853 times.
✓ Branch 1 taken 1032207 times.
 8882060 if ( ( ( (x-xhi)*df-f)*((x-xlo)*df-f) >= 0.0
58
2/2
✓ Branch 0 taken 1069 times.
✓ Branch 1 taken 7848784 times.
 7849853 || (fabs(2.0*f) > fabs(dxold*df) ) ) )
59 {
60 1033276 dxold = dx;
61 1033276 dx = 0.5*(xhi-xlo);
62 1033276 x = xlo + dx;
63
2/2
✓ Branch 0 taken 252 times.
✓ Branch 1 taken 1033024 times.
 1033276 if ( xlo == x ) break;
64 }
65
66 // Newton step acceptable. Take it.
67 else
68 {
69 7848784 dxold = dx;
70 7848784 dx = f/df;
71 7848784 temp = x;
72 7848784 x -= dx;
73
2/2
✓ Branch 0 taken 13755 times.
✓ Branch 1 taken 7835029 times.
 7848784 if ( temp == x ) break;
74 }
75
76 // Convergence criterion.
77
2/2
✓ Branch 0 taken 1819971 times.
✓ Branch 1 taken 7048082 times.
 8868053 if ( fabs(dx) < xacc ) break;
78
79 // Evaluate function. Maintain bracket on the root.
80 7048082 func(x, &f, &df, p);
81 7048082 n++;
82
2/2
✓ Branch 0 taken 1003943 times.
✓ Branch 1 taken 6044139 times.
 7048082 if ( f < 0.0 ) xlo = x;
83 6044139 else xhi = x;
84 }
85 1833978 *rts = x;
86
1/2
✓ Branch 0 taken 1833978 times.
✗ Branch 1 not taken.
 1833978 if ( n <= MAXIT) return n;
87 ✗ else return 0;
88 };
89
90
91 322566 double findroot_Ridder(double x1, double x2, double xacc,
92 double (*func)(double, void* p), void* p)
93 {
94 int j;
95 double ans, fhi, flo, fm, fnew, s, xhi, xlo, xm, xnew;
96
97 322566 flo = func(x1, p);
98 322566 fhi = func(x2, p);
99
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 322566 times.
 322566 if ( flo == 0.0 ) return x1;
100
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 322566 times.
 322566 if ( fhi == 0.0 ) return x2;
101 322566 ans = 0.5*(x1+x2);
102
7/8
✓ Branch 0 taken 27617 times.
✓ Branch 1 taken 294949 times.
✓ Branch 2 taken 288 times.
✓ Branch 3 taken 27329 times.
✓ Branch 4 taken 294949 times.
✓ Branch 5 taken 288 times.
✓ Branch 6 taken 294949 times.
✗ Branch 7 not taken.
 322566 if ( (flo > 0.0 && fhi < 0.0) || (flo < 0.0 && fhi > 0.0) )
103 {
104 322278 xlo = x1;
105 322278 xhi = x2;
106
1/2
✓ Branch 0 taken 911937 times.
✗ Branch 1 not taken.
 911937 for (j=1; j<=MAXIT; j++) {
107 911937 xm = 0.5*(xlo + xhi);
108 911937 fm = func(xm, p);
109 911937 s = sqrt( fm*fm - flo*fhi );
110
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 911937 times.
 911937 if (s == 0.0) return ans;
111
2/2
✓ Branch 0 taken 559681 times.
✓ Branch 1 taken 352256 times.
 911937 xnew = xm + (xm-xlo)*( (flo >= fhi ? 1.0 : -1.0)*fm/s );
112
2/2
✓ Branch 0 taken 322278 times.
✓ Branch 1 taken 589659 times.
 911937 if ( fabs(xnew - ans) <= xacc ) break;
113 589659 ans = xnew;
114 589659 fnew = func(ans, p);
115
4/4
✓ Branch 0 taken 85617 times.
✓ Branch 1 taken 504042 times.
✓ Branch 2 taken 544207 times.
✓ Branch 3 taken 45452 times.
 589659 if ( SIGN(fm, fnew) != fm)
116 {
117 544207 xlo = xm;
118 544207 flo = fm;
119 544207 xhi = ans;
120 544207 fhi = fnew;
121 }
122
4/4
✓ Branch 0 taken 28310 times.
✓ Branch 1 taken 17142 times.
✓ Branch 2 taken 17142 times.
✓ Branch 3 taken 28310 times.
 45452 else if ( SIGN(flo, fnew) != flo )
123 {
124 17142 xhi = ans;
125 17142 fhi = fnew;
126 }
127
2/4
✓ Branch 0 taken 28310 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 28310 times.
✗ Branch 3 not taken.
 28310 else if ( SIGN(fhi, fnew) != fhi)
128 {
129 28310 xlo = ans;
130 28310 flo = fnew;
131 }
132 ✗ else return ans;
133
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 589659 times.
 589659 if ( fabs(xhi - xlo) <= xacc ) return ans;
134 }
135 322278 return ans;
136 }
137 288 return -1.e20;
138 }
139