Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80601078
dfsane_solver.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sun, Sep 1, 14:52
Size
5 KB
Mime Type
text/x-c
Expires
Tue, Sep 3, 14:52 (2 d)
Engine
blob
Format
Raw Data
Handle
20404637
Attached To
rTAMAAS tamaas
dfsane_solver.cpp
View Options
/**
* @file
* LICENSE
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "dfsane_solver.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
DFSANESolver::DFSANESolver(Residual& residual)
: EPSolver(residual),
search_direction(residual.getVector().dataSize(),
residual.getVector().getNbComponents()),
previous_residual(residual.getVector().dataSize(),
residual.getVector().getNbComponents()),
current_x(residual.getVector().dataSize(),
residual.getVector().getNbComponents()),
delta_x(residual.getVector().dataSize(),
residual.getVector().getNbComponents()),
delta_residual(residual.getVector().dataSize(),
residual.getVector().getNbComponents()) {}
void DFSANESolver::solve() {
EPSolver::beforeSolve();
auto& x_prev = *_x;
UInt k = 0, nmax = 500;
Real F_norm = 0;
previous_merits.clear();
_residual.computeResidual(x_prev);
previous_residual = _residual.getVector();
const auto F_norm0 = previous_residual.l2norm();
eta = [F_norm0](UInt k) { return F_norm0 / ((1 + k) * (1 + k)); };
const std::pair<Real, Real> sigma_extrema{1e-10, 1e10};
Real sigma = 1;
do {
computeSearchDirection(sigma);
lineSearch(eta(k));
// No need for a solution update
// it was done automatically in line search
sigma = computeSpectralCoeff(sigma_extrema);
x_prev = current_x;
previous_residual = _residual.getVector();
F_norm = previous_residual.l2norm();
Logger().get(LogLevel::debug) << TAMAAS_DEBUG_MSG(k << ' ' << F_norm);
} while (F_norm > getTolerance() and k++ < nmax);
if (k >= nmax)
TAMAAS_EXCEPTION("DF-SANE did not converge");
Logger().get(LogLevel::info) << "DF-SANE: sucessful convergence (" << k
<< " iterations, " << getTolerance() << ")\n";
_residual.computeResidualDisplacement(*_x);
}
Real DFSANESolver::computeSpectralCoeff(const std::pair<Real, Real>& bounds) {
delta_x = current_x;
delta_x -= *_x;
delta_residual = _residual.getVector();
delta_residual -= previous_residual;
auto trial = delta_x.dot(delta_x) / delta_x.dot(delta_residual);
if (std::abs(trial) >= bounds.first and std::abs(trial) <= bounds.second)
return trial;
else {
auto F = previous_residual.l2norm();
// Rule from Cruz et al. (2006)
if (F > 1)
return 1;
else if (F < 1e-5)
return 1e5;
else
return 1 / F;
}
}
void DFSANESolver::computeSearchDirection(Real sigma) {
_residual.computeResidual(*_x);
search_direction = _residual.getVector();
search_direction *= -sigma;
}
void DFSANESolver::lineSearch(Real eta_k) {
const UInt M = 10;
const Real gamma = 1e-4;
Real alphap = 1, alpham = 1;
const std::pair<Real, Real> tau_extrema{0.1, 0.5};
// Merit function (nexp = 2)
const auto f = [&](Real alpha, Real sign) {
current_x = search_direction;
current_x *= alpha * sign;
current_x += *_x;
_residual.computeResidual(current_x);
return _residual.getVector().dot(_residual.getVector());
};
// Manage previous merits
const auto fk = _residual.getVector().dot(_residual.getVector());
previous_merits.push_front(fk);
if (previous_merits.size() > M)
previous_merits.pop_back();
const auto fbar =
*std::max_element(previous_merits.begin(), previous_merits.end());
do {
Real fp = f(alphap, +1);
if (fp < fbar + eta_k - gamma * alphap * alphap * fk)
return;
Real fm = f(alpham, -1);
if (fm < fbar + eta_k - gamma * alpham * alpham * fk)
return;
alphap = computeAlpha(alphap, fp, fk, tau_extrema);
alpham = computeAlpha(alpham, fm, fk, tau_extrema);
if (std::isnan(alphap) or std::isnan(alpham))
TAMAAS_EXCEPTION("NaN error in line search");
} while (true);
}
Real DFSANESolver::computeAlpha(const Real alpha, Real f, Real fk,
const std::pair<Real, Real>& bounds) {
const auto alphat = alpha * alpha * fk / (f + (2 * alpha - 1) * fk);
if (alphat < bounds.first * alpha)
return bounds.first * alpha;
else if (alphat > bounds.second * alpha)
return bounds.second * alpha;
else
return alphat;
}
} // namespace tamaas
Event Timeline
Log In to Comment