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
16 namespace plask { namespace optical { namespace effective {
17
18#define carg(x) realaxis? dcomplex(x, imag(start)) : dcomplex(real(start), x)
19#define fun(x) abs(valFunction(carg(x)))
20
21 double RootBrent::axisBrent(dcomplex start, double& fx, bool realaxis) const
22{
23 const double C = 0.5 * (3. - sqrt(5.));
24 const double G = 2. / (sqrt(5.) - 1.);
25
26 int i = 0;
27 double x, dist;
28
29 if (realaxis) {
30 x = real(start);
31 dist = real(params.initial_dist);
32 } else {
33 x = imag(start);
34 if (imag(params.initial_dist) == 0) dist = real(params.initial_dist);
35 else dist = imag(params.initial_dist);
36 }
37
38 // Exponential search
39 double a = x - dist,
40 b = x + dist;
41
42 writelog(LOG_DETAIL, "Searching for the root with Brent method between {0} and {1} in the {2} axis",
43 str(a), str(b), realaxis?"real":"imaginary");
44
45 if (isnan(fx)) fx = fun(x);
46
47 double fw = fun(a),
48 fv = fun(b);
49 log_value(carg(a), fw);
50 log_value(carg(x), fx);
51 log_value(carg(b), fv);
52
53 double d, u;
54
55 log_value.resetCounter();
56
57 bool bounded = false;
58 if (fw <= fx && fx <= fv) {
59 writelog(LOG_DETAIL, "Extending search range to lower values");
60 for (; i < params.maxiter; ++i) {
61 u = a - G * (x - a);
62 b = x; // fv = fx;
63 x = a; fx = fw;
64 a = u; fw = fun(a);
65 log_value.count(a, fw);
66 if (fw > fx) {
67 bounded = true;
68 break;
69 }
70 }
71 if (bounded) writelog(LOG_DETAIL, "Searching minimum in range between {0} and {1}", a, b);
72 } else if (fw >= fx && fx >= fv) {
73 writelog(LOG_DETAIL, "Extending search range to higher values");
74 for (; i < params.maxiter; ++i) {
75 u = b + G * (b - x);
76 a = x; // fw = fx;
77 x = b; fx = fv;
78 b = u; fv = fun(b);
79 log_value.count(carg(b), fv);
80 if (fv > fx) {
81 bounded = true;
82 break;
83 }
84 }
85 if (bounded) writelog(LOG_DETAIL, "Searching minimum in range between {0} and {1}", a, b);
86 } else bounded = true;
87
88 if (!bounded)
89 throw ComputationError(solver.getId(),
90 "Brent: {0}: minimum still unbounded after maximum number of iterations",
91 log_value.chartName());
92
93 double sa = a, sb = b, w = x, v = x, e = 0.0;
94 fw = fv = fx;
95
96 // Reference:
97 // Richard Brent,
98 // Algorithms for Minimization Without Derivatives,
99 // Dover, 2002,
100 // ISBN: 0-486-41998-3,
101 // LC: QA402.5.B74.
102
103 double tol = SMALL * abs(x) + params.tolx;
104 double tol2 = 2.0 * tol;
105
106 for (; i < params.maxiter; ++i) {
107 double m = 0.5 * (sa + sb) ;
108
109 // Check the stopping criterion.
110 if (abs(x - m) <= tol2 - 0.5 * (sb - sa)) return x;
111
112 // Fit a parabola.
113 double r = 0., q = 0., p = 0.;
114 if (tol < abs(e)) {
115 r = (x - w) * (fx - fv);
116 q = (x - v) * (fx - fw);
117 p = (x - v) * q - (x - w) * r;
118 q = 2.0 * (q - r);
119 if (0.0 < q) p = - p;
120 q = abs(q);
121 r = e;
122 e = d;
123 }
124 if (abs(p) < abs(0.5 * q * r) && q * (sa - x) < p && p < q * (sb - x)) {
125 // Take the parabolic interpolation step
126 d = p / q;
127 u = x + d;
128 // F must not be evaluated too close to a or b
129 if ((u - sa) < tol2 || (sb - u) < tol2) {
130 if (x < m) d = tol;
131 else d = - tol;
132 }
133 } else {
134 // A golden-section step.
135 if (x < m) e = sb - x;
136 else e = sa - x;
137 d = C * e;
138 }
139
140 // F must not be evaluated too close to x
141 if (tol <= abs(d)) u = x + d;
142 else if (0.0 < d) u = x + tol;
143 else u = x - tol;
144
145 // Update a, b, v, w, and x
146 double fu = fun(u);
147 if (fu <= fx) {
148 if (u < x) sb = x;
149 else sa = x;
150 v = w; fv = fw;
151 w = x; fw = fx;
152 x = u; fx = fu;
153 } else {
154 if (u < x) sa = u;
155 else sb = u;
156 if (fu <= fw || w == x) {
157 v = w; fv = fw;
158 w = u; fw = fu;
159 } else if (fu <= fv || v == x || v == w) {
160 v = u; fv = fu;
161 }
162 }
163 log_value.count(carg(x), fx);
164 }
165 throw ComputationError(solver.getId(), "Brent: {0}: maximum number of iterations reached", log_value.chartName());
166 return 0;
167}
168
169
170//**************************************************************************
172dcomplex RootBrent::find(dcomplex xstart) const
173{
174 double f0 = NAN;
175
176 xstart.real(axisBrent(xstart, f0, true));
177 for (int i = 0; i < params.stairs; ++i) {
178 xstart.imag(axisBrent(xstart, f0, false));
179 xstart.real(axisBrent(xstart, f0, true));
180 }
181
182 if (f0 > params.tolf_min)
183 ComputationError(solver.getId(),
184 "Brent: {0}: After real and imaginary minimum search, determinant still not small enough",
185 log_value.chartName());
186 return xstart;
187}
188
189}}} // namespace plask::optical::effective