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