PLaSK library

Loading...
Searching...
No Matches
brent.cpp
Go to the documentation of this file.
1/*
2 * This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
3 * Copyright (c) 2022 Lodz University of Technology
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, version 3.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 */
14#include "brent.hpp"
15#include "rootdigger-impl.hpp"
16
17 namespace plask { namespace optical { namespace modal {
18
19#define carg(x) realaxis? dcomplex(x, imag(start)) : dcomplex(real(start), x)
20#define fun(x) abs(valFunction(carg(x)))
21
22 double RootBrent::axisBrent(dcomplex start, double& fx, bool realaxis, int& counter)
23{
24 const double C = 0.5 * (3. - sqrt(5.));
25 const double G = 2. / (sqrt(5.) - 1.);
26
27 int i = 0;
28 double x, dist;
29
30 if (realaxis) {
31 x = real(start);
32 dist = real(params.initial_dist);
33 } else {
34 x = imag(start);
35 if (imag(params.initial_dist) == 0) dist = real(params.initial_dist);
36 else dist = imag(params.initial_dist);
37 }
38
39 // Exponential search
40 double a = x - dist,
41 b = x + dist;
42
43 writelog(LOG_DETAIL, "Searching for the root with Brent method between {0} and {1} in the {2} axis",
44 str(a), str(b), realaxis?"real":"imaginary");
45
46 if (isnan(fx)) fx = fun(x);
47
48 double fw = fun(a),
49 fv = fun(b);
50 log_value(carg(a), fw);
51 log_value(carg(x), fx);
52 log_value(carg(b), fv);
53
54 double d, u;
55
56 log_value.resetCounter();
57
58 bool bounded = false;
59 if (fw <= fx && fx <= fv) {
60 writelog(LOG_DETAIL, "Extending search range to lower values");
61 for (; i < params.maxiter; ++i) {
62 u = a - G * (x - a);
63 b = x; // fv = fx;
64 x = a; fx = fw;
65 a = u; fw = fun(a);
66 counter++;
67 log_value.count(carg(b), fw);
68 if (fw > fx) {
69 bounded = true;
70 break;
71 }
72 }
73 if (bounded) writelog(LOG_DETAIL, "Searching minimum in range between {0} and {1}", a, b);
74 } else if (fw >= fx && fx >= fv) {
75 writelog(LOG_DETAIL, "Extending search range to higher values");
76 for (; i < params.maxiter; ++i) {
77 u = b + G * (b - x);
78 a = x; // fw = fx;
79 x = b; fx = fv;
80 b = u; fv = fun(b);
81 counter++;
82 log_value.count(carg(b), fv);
83 if (fv > fx) {
84 bounded = true;
85 break;
86 }
87 }
88 if (bounded) writelog(LOG_DETAIL, "Searching minimum in range between {0} and {1}", a, b);
89 } else bounded = true;
90
91 if (!bounded)
92 throw ComputationError(solver.getId(),
93 "Brent: {0}: minimum still unbounded after maximum number of iterations",
94 log_value.chartName());
95
96 double sa = a, sb = b, w = x, v = x, e = 0.0;
97 fw = fv = fx;
98
99 // Reference:
100 // Richard Brent,
101 // Algorithms for Minimization Without Derivatives,
102 // Dover, 2002,
103 // ISBN: 0-486-41998-3,
104 // LC: QA402.5.B74.
105
106 double tol = SMALL * abs(x) + params.tolx;
107 double tol2 = 2.0 * tol;
108
109 for (; i < params.maxiter; ++i) {
110 double m = 0.5 * (sa + sb) ;
111
112 // Check the stopping criterion.
113 if (abs(x - m) <= tol2 - 0.5 * (sb - sa)) return x;
114
115 // Fit a parabola.
116 double r = 0., q = 0., p = 0.;
117 if (tol < abs(e)) {
118 r = (x - w) * (fx - fv);
119 q = (x - v) * (fx - fw);
120 p = (x - v) * q - (x - w) * r;
121 q = 2.0 * (q - r);
122 if (0.0 < q) p = - p;
123 q = abs(q);
124 r = e;
125 e = d;
126 }
127 if (abs(p) < abs(0.5 * q * r) && q * (sa - x) < p && p < q * (sb - x)) {
128 // Take the parabolic interpolation step
129 d = p / q;
130 u = x + d;
131 // F must not be evaluated too close to a or b
132 if ((u - sa) < tol2 || (sb - u) < tol2) {
133 if (x < m) d = tol;
134 else d = - tol;
135 }
136 } else {
137 // A golden-section step.
138 if (x < m) e = sb - x;
139 else e = sa - x;
140 d = C * e;
141 }
142
143 // F must not be evaluated too close to x
144 if (tol <= abs(d)) u = x + d;
145 else if (0.0 < d) u = x + tol;
146 else u = x - tol;
147
148 // Update a, b, v, w, and x
149 double fu = fun(u);
150 if (fu <= fx) {
151 if (u < x) sb = x;
152 else sa = x;
153 v = w; fv = fw;
154 w = x; fw = fx;
155 x = u; fx = fu;
156 } else {
157 if (u < x) sa = u;
158 else sb = u;
159 if (fu <= fw || w == x) {
160 v = w; fv = fw;
161 w = u; fw = fu;
162 } else if (fu <= fv || v == x || v == w) {
163 v = u; fv = fu;
164 }
165 }
166 counter++;
167 log_value.count(carg(x), fx);
168 }
169 throw ComputationError(solver.getId(), "Brent: {0}: maximum number of iterations reached", log_value.chartName());
170 return 0;
171}
172
173
174//**************************************************************************
176dcomplex RootBrent::find(dcomplex xstart)
177{
178 double f0 = NAN;
179 dcomplex xprev = NAN;
180 double tolx2 = params.tolx * params.tolx;
181
182 int counter = 0;
183
184 while (counter < params.maxiter && !(f0 <= params.tolf_min || abs2(xstart - xprev) <= tolx2)) {
185 xprev = xstart;
186 xstart.real(axisBrent(xstart, f0, true, counter));
187 xstart.imag(axisBrent(xstart, f0, false, counter));
188 }
189
190 if (f0 > params.tolf_max)
191 ComputationError(solver.getId(),
192 "Brent: {0}: After real and imaginary minimum search, determinant still not small enough",
193 log_value.chartName());
194 return xstart;
195}
196
197}}} // namespace plask::optical::modal