PLaSK library
Loading...
Searching...
No Matches
broyden.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 "broyden.hpp"
15#include "rootdigger-impl.hpp"
16using namespace std;
17
18namespace plask { namespace optical { namespace modal {
19
20//**************************************************************************
22dcomplex RootBroyden::find(dcomplex start)
23{
24 writelog(LOG_DETAIL, "Searching for the root with Broyden method starting from " + str(start));
25 log_value.resetCounter();
26 dcomplex x = Broyden(start);
27 writelog(LOG_RESULT, "Found root at " + str(x));
28 return x;
29}
30
31//**************************************************************************
32//**************************************************************************
33//******** The Broyden globally convergent method for root finding *********
34//**************************************************************************
35
36//**************************************************************************
37// Return Jacobian of F(x)
38void RootBroyden::fdjac(dcomplex x, dcomplex F, dcomplex& Jr, dcomplex& Ji)
39{
40 double xr0 = real(x), xi0 = imag(x);
41 double hr = EPS*abs(xr0), hi = EPS*abs(xi0);
42 if (hr == 0.0) hr = EPS;
43 if (hi == 0.0) hi = EPS;
44
45 double xr1 = xr0 + hr, xi1 = xi0 + hi;
46 hr = xr1 - xr0; hi = xi1 - xi0; // trick to reduce finite precision error
47
48 dcomplex xr = dcomplex(xr1, xi0), xi = dcomplex(xr0, xi1);
49 dcomplex Fr = valFunction(xr); log_value(xr, Fr);
50 dcomplex Fi = valFunction(xi); log_value(xi, Fi);
51
52 Jr = (Fr - F) / hr;
53 Ji = (Fi - F) / hi;
54}
55
56//**************************************************************************
57// Search for the new point x along direction p for which
58// functional f decreased sufficiently
59// g - (approximate) gradient of 1/2(F*F), stpmax - maximum allowed step
60// return true if performed step or false if could not find sufficient function decrease
61bool RootBroyden::lnsearch(dcomplex& x, dcomplex& F, dcomplex g, dcomplex p, double stpmax)
62{
63 if (double absp=abs(p) > stpmax) p *= stpmax/absp; // Ensure step <= stpmax
64
65 double slope = real(g)*real(p) + imag(g)*imag(p); // slope = grad(f)*p
66
67 // Compute the functional
68 double f = 0.5 * (real(F)*real(F) + imag(F)*imag(F));
69
70 // Remember original values
71 dcomplex x0 = x;
72 double f0 = f;
73
74 double lambda = 1.0; // lambda parameter x = x0 + lambda*p
75 double lambda1, lambda2 = 0., f2 = 0.;
76
77 bool first = true;
78
79 while(true) {
80 if (lambda < params.lambda_min) { // we have (possible) convergence of x
81 x = x0; // f = f0;
82 return false;
83 }
84
85 x = x0 + lambda*p;
86 F = valFunction(x);
87 log_value.count(x, F);
88
89 f = 0.5 * (real(F)*real(F) + imag(F)*imag(F));
90 if (std::isnan(f)) throw ComputationError(solver.getId(), "computed value is NaN");
91
92 if (f < f0 + params.alpha*lambda*slope) // sufficient function decrease
93 return true;
94
95 lambda1 = lambda;
96
97 if (first) { // first backtrack
98 lambda = -slope / (2. * (f-f0-slope));
99 first = false;
100 } else { // subsequent backtracks
101 double rsh1 = f - f0 - lambda1*slope;
102 double rsh2 = f2 - f0 - lambda2*slope;
103 double a = (rsh1 / (lambda1*lambda1) - rsh2 / (lambda2*lambda2)) / (lambda1-lambda2);
104 double b = (-lambda2 * rsh1 / (lambda1*lambda1) + lambda1 * rsh2 / (lambda2*lambda2))
105 / (lambda1-lambda2);
106
107 if (a == 0.0)
108 lambda = -slope/(2.*b);
109 else {
110 double delta = b*b - 3.*a*slope;
111 if (delta < 0.0) throw ComputationError(solver.getId(), "broyden lnsearch: roundoff problem");
112 lambda = (-b + sqrt(delta)) / (3.0*a);
113 }
114 }
115
116 lambda2 = lambda1; f2 = f; // store the second last parameters
117
118 lambda = max(lambda, 0.1*lambda1); // guard against too fast decrease of lambda
119
120 writelog(LOG_DETAIL, "Broyden step decreased to the fraction " + str(lambda) + " of the original step");
121 }
122
123}
124
125//**************************************************************************
126// Search for the root of char_val using globally convergent Broyden method
127// starting from point x
128dcomplex RootBroyden::Broyden(dcomplex x)
129{
130 // Compute the initial guess of the function (and check for the root)
131 dcomplex F = valFunction(x);
132 double absF = abs(F);
133 log_value.count(x, F);
134 if (absF < params.tolf_min) return x;
135
136 bool restart = true; // do we have to recompute Jacobian?
137 bool trueJacobian; // did we recently update Jacobian?
138
139 dcomplex Br, Bi; // Broyden matrix columns
140 dcomplex dF, dx; // Computed dist
141
142 dcomplex oldx, oldF;
143
144 // Main loop
145 for (int i = 0; i < params.maxiter; i++) {
146 oldx = x; oldF = F;
147
148 if (restart) { // compute Broyden matrix as a Jacobian
149 fdjac(x, F, Br, Bi);
150 restart = false;
151 trueJacobian = true;
152 } else { // update Broyden matrix
153 dcomplex dB = dF - dcomplex(real(Br)*real(dx)+real(Bi)*imag(dx), imag(Br)*real(dx)+imag(Bi)*imag(dx));
154 double m = (real(dx)*real(dx) + imag(dx)*imag(dx));
155 Br += (dB * real(dx)) / m;
156 Bi += (dB * imag(dx)) / m;
157 trueJacobian = false;
158 }
159
160 // compute g ~= B^T * F
161 dcomplex g = dcomplex(real(Br)*real(F)+imag(Br)*imag(F), real(Bi)*real(F)+imag(Bi)*imag(F));
162
163 // compute p = - B**(-1) * F
164 double M = real(Br)*imag(Bi) - imag(Br)*real(Bi);
165 if (M == 0) throw ComputationError(solver.getId(), "singular Jacobian in Broyden method");
166 dcomplex p = - dcomplex(real(F)*imag(Bi)-imag(F)*real(Bi), real(Br)*imag(F)-imag(Br)*real(F)) / M;
167
168 // find the right step
169 if (lnsearch(x, F, g, p, params.maxstep)) { // found sufficient functional decrease
170 dx = x - oldx;
171 dF = F - oldF;
172 if ((abs(dx) < params.tolx && abs(F) < params.tolf_max) || abs(F) < params.tolf_min)
173 return x; // convergence!
174 } else {
175 if (abs(F) < params.tolf_max) // convergence!
176 return x;
177 else if (!trueJacobian) { // first try reinitializing the Jacobian
178 writelog(LOG_DETAIL, "Reinitializing Jacobian");
179 restart = true;
180 } else { // either spurious convergence (local minimum) or failure
181 throw ComputationError(solver.getId(), "Broyden method failed to converge");
182 }
183 }
184 }
185
186 throw ComputationError(solver.getId(), "Broyden: maximum number of iterations reached");
187}
188
189}}} // namespace plask::optical::modal