639 lines
30 KiB
C++
639 lines
30 KiB
C++
module;
|
|
|
|
#include <algorithm>
|
|
#include <chrono>
|
|
#include <cmath>
|
|
#include <complex>
|
|
#include <cstdint>
|
|
#include <limits>
|
|
#include <memory>
|
|
#include <ranges>
|
|
#include <stdexcept>
|
|
#include <string>
|
|
#include <utility>
|
|
#include <vector>
|
|
|
|
#include <Eigen/Dense>
|
|
#include <Eigen/Eigenvalues>
|
|
#include <Eigen/SVD>
|
|
#include <mfem.hpp>
|
|
#include <mpi.h>
|
|
|
|
module mean_field;
|
|
|
|
import :solver.preconditioning_diagnostics;
|
|
|
|
namespace {
|
|
using Clock = std::chrono::steady_clock;
|
|
|
|
[[nodiscard]] double seconds_between(
|
|
const Clock::time_point start,
|
|
const Clock::time_point finish
|
|
) {
|
|
return std::chrono::duration<double>(finish - start).count();
|
|
}
|
|
|
|
void verify_finite_vector(
|
|
const mfem::Vector &vector,
|
|
const char *message
|
|
) {
|
|
for (int index = 0; index < vector.Size(); ++index) {
|
|
if (!std::isfinite(vector(index))) {
|
|
throw std::invalid_argument(message);
|
|
}
|
|
}
|
|
}
|
|
|
|
[[nodiscard]] double global_dot(
|
|
const mfem::Vector &left,
|
|
const mfem::Vector &right,
|
|
const MPI_Comm communicator
|
|
) {
|
|
if (communicator == MPI_COMM_NULL) {
|
|
throw std::invalid_argument("Preconditioning diagnostics require a valid MPI communicator.");
|
|
}
|
|
if (left.Size() != right.Size()) {
|
|
throw std::invalid_argument("A distributed inner product received vectors with different sizes.");
|
|
}
|
|
|
|
const double localValue = left * right;
|
|
double globalValue = 0.0;
|
|
MPI_Allreduce(&localValue, &globalValue, 1, MPI_DOUBLE, MPI_SUM, communicator);
|
|
return globalValue;
|
|
}
|
|
|
|
[[nodiscard]] double global_norm(
|
|
const mfem::Vector &vector,
|
|
const MPI_Comm communicator
|
|
) {
|
|
return std::sqrt(std::max(global_dot(vector, vector, communicator), 0.0));
|
|
}
|
|
|
|
[[nodiscard]] mean_field::solver::OperatorApplicationStatistics maximum_rank_statistics(
|
|
const mean_field::solver::OperatorApplicationStatistics &local,
|
|
const MPI_Comm communicator
|
|
) {
|
|
unsigned long long localApplications = static_cast<unsigned long long>(local.applications);
|
|
unsigned long long maximumApplications{0};
|
|
MPI_Allreduce(&localApplications, &maximumApplications, 1, MPI_UNSIGNED_LONG_LONG, MPI_MAX, communicator);
|
|
|
|
mean_field::solver::OperatorApplicationStatistics result;
|
|
result.applications = static_cast<std::uint64_t>(maximumApplications);
|
|
MPI_Allreduce(&local.totalSeconds, &result.totalSeconds, 1, MPI_DOUBLE, MPI_MAX, communicator);
|
|
MPI_Allreduce(&local.maximumSeconds, &result.maximumSeconds, 1, MPI_DOUBLE, MPI_MAX, communicator);
|
|
return result;
|
|
}
|
|
|
|
[[nodiscard]] double maximum_rank_value(
|
|
const double localValue,
|
|
const MPI_Comm communicator
|
|
) {
|
|
double result = 0.0;
|
|
MPI_Allreduce(&localValue, &result, 1, MPI_DOUBLE, MPI_MAX, communicator);
|
|
return result;
|
|
}
|
|
|
|
[[nodiscard]] mean_field::solver::PreconditionerLifecycleStatistics maximum_rank_lifecycle_statistics(
|
|
const mean_field::solver::PreconditionerLifecycleStatistics &local,
|
|
const MPI_Comm communicator
|
|
) {
|
|
unsigned long long localSetups = static_cast<unsigned long long>(local.setups);
|
|
unsigned long long localRefreshes = static_cast<unsigned long long>(local.refreshes);
|
|
unsigned long long maximumSetups{0};
|
|
unsigned long long maximumRefreshes{0};
|
|
MPI_Allreduce(&localSetups, &maximumSetups, 1, MPI_UNSIGNED_LONG_LONG, MPI_MAX, communicator);
|
|
MPI_Allreduce(&localRefreshes, &maximumRefreshes, 1, MPI_UNSIGNED_LONG_LONG, MPI_MAX, communicator);
|
|
|
|
mean_field::solver::PreconditionerLifecycleStatistics result;
|
|
result.setups = static_cast<std::uint64_t>(maximumSetups);
|
|
result.refreshes = static_cast<std::uint64_t>(maximumRefreshes);
|
|
MPI_Allreduce(&local.setupSeconds, &result.setupSeconds, 1, MPI_DOUBLE, MPI_MAX, communicator);
|
|
MPI_Allreduce(&local.refreshSeconds, &result.refreshSeconds, 1, MPI_DOUBLE, MPI_MAX, communicator);
|
|
return result;
|
|
}
|
|
|
|
[[nodiscard]] Eigen::MatrixXd copy_hessenberg(
|
|
const Eigen::MatrixXd &source,
|
|
const int rowCount,
|
|
const int columnCount
|
|
) {
|
|
return source.topLeftCorner(rowCount, columnCount);
|
|
}
|
|
} // namespace
|
|
|
|
namespace mean_field::solver {
|
|
InstrumentedOperator::InstrumentedOperator(const mfem::Operator &operation)
|
|
: mfem::Operator(
|
|
operation.Height(),
|
|
operation.Width()
|
|
),
|
|
m_operation(std::addressof(operation)) {
|
|
}
|
|
|
|
void InstrumentedOperator::Mult(
|
|
const mfem::Vector &input,
|
|
mfem::Vector &output
|
|
) const {
|
|
const Clock::time_point start = Clock::now();
|
|
m_operation->Mult(input, output);
|
|
const double elapsed = seconds_between(start, Clock::now());
|
|
|
|
++m_statistics.applications;
|
|
m_statistics.totalSeconds += elapsed;
|
|
m_statistics.maximumSeconds = std::max(m_statistics.maximumSeconds, elapsed);
|
|
}
|
|
|
|
void InstrumentedOperator::ResetStatistics() const noexcept {
|
|
m_statistics = {};
|
|
}
|
|
|
|
const OperatorApplicationStatistics &InstrumentedOperator::GetStatistics() const noexcept {
|
|
return m_statistics;
|
|
}
|
|
|
|
const mfem::Operator &InstrumentedOperator::GetOperation() const noexcept {
|
|
return *m_operation;
|
|
}
|
|
|
|
InstrumentedPreconditioner::InstrumentedPreconditioner(mfem::Solver &preconditioner)
|
|
: mfem::Solver(
|
|
preconditioner.Height(),
|
|
preconditioner.Width(),
|
|
preconditioner.iterative_mode
|
|
),
|
|
m_preconditioner(std::addressof(preconditioner)) {
|
|
}
|
|
|
|
void InstrumentedPreconditioner::SetOperator(const mfem::Operator &operation) {
|
|
const Clock::time_point start = Clock::now();
|
|
m_preconditioner->SetOperator(operation);
|
|
m_lifecycleStatistics.setupSeconds += seconds_between(start, Clock::now());
|
|
++m_lifecycleStatistics.setups;
|
|
if (m_preconditioner->Height() != Height() || m_preconditioner->Width() != Width()) {
|
|
throw std::invalid_argument("An instrumented preconditioner changed dimensions during SetOperator.");
|
|
}
|
|
}
|
|
|
|
void InstrumentedPreconditioner::Mult(
|
|
const mfem::Vector &input,
|
|
mfem::Vector &output
|
|
) const {
|
|
const Clock::time_point start = Clock::now();
|
|
m_preconditioner->Mult(input, output);
|
|
const double elapsed = seconds_between(start, Clock::now());
|
|
|
|
++m_statistics.applications;
|
|
m_statistics.totalSeconds += elapsed;
|
|
m_statistics.maximumSeconds = std::max(m_statistics.maximumSeconds, elapsed);
|
|
}
|
|
|
|
void InstrumentedPreconditioner::ResetStatistics() const noexcept {
|
|
m_statistics = {};
|
|
}
|
|
|
|
const OperatorApplicationStatistics &InstrumentedPreconditioner::GetStatistics() const noexcept {
|
|
return m_statistics;
|
|
}
|
|
|
|
const PreconditionerLifecycleStatistics &InstrumentedPreconditioner::GetLifecycleStatistics() const noexcept {
|
|
return m_lifecycleStatistics;
|
|
}
|
|
|
|
const mfem::Solver &InstrumentedPreconditioner::GetPreconditioner() const noexcept {
|
|
return *m_preconditioner;
|
|
}
|
|
|
|
IdentityPreconditioner::IdentityPreconditioner(const int size) : mfem::Solver(size) {
|
|
if (size <= 0) {
|
|
throw std::invalid_argument("An identity preconditioner requires a positive dimension.");
|
|
}
|
|
}
|
|
|
|
void IdentityPreconditioner::SetOperator(const mfem::Operator &operation) {
|
|
if (operation.Height() != Height() || operation.Width() != Width()) {
|
|
throw std::invalid_argument("The identity preconditioner received an incompatible operator.");
|
|
}
|
|
}
|
|
|
|
void IdentityPreconditioner::Mult(
|
|
const mfem::Vector &input,
|
|
mfem::Vector &output
|
|
) const {
|
|
if (input.Size() != Width()) {
|
|
throw std::invalid_argument("The identity preconditioner received an input with the wrong size.");
|
|
}
|
|
output = input;
|
|
}
|
|
|
|
FixedRightPreconditionedOperator::FixedRightPreconditionedOperator(
|
|
const mfem::Operator &jacobian,
|
|
const mfem::Solver &inversePreconditioner
|
|
)
|
|
: mfem::Operator(
|
|
jacobian.Height(),
|
|
inversePreconditioner.Width()
|
|
),
|
|
m_jacobian(std::addressof(jacobian)),
|
|
m_inversePreconditioner(std::addressof(inversePreconditioner)),
|
|
m_preconditionedDirection(inversePreconditioner.Height()) {
|
|
if (jacobian.Height() != jacobian.Width()) {
|
|
throw std::invalid_argument("A preconditioned stellar Jacobian must be square.");
|
|
}
|
|
if (inversePreconditioner.Height() != jacobian.Width() || inversePreconditioner.Width() != jacobian.Height()) {
|
|
throw std::invalid_argument("The inverse preconditioner does not map residuals into Jacobian states.");
|
|
}
|
|
if (Height() != Width()) {
|
|
throw std::invalid_argument("The fixed right-preconditioned product must be square.");
|
|
}
|
|
}
|
|
|
|
void FixedRightPreconditionedOperator::Mult(
|
|
const mfem::Vector &input,
|
|
mfem::Vector &output
|
|
) const {
|
|
if (input.Size() != Width()) {
|
|
throw std::invalid_argument("The right-preconditioned operator received an input with the wrong size.");
|
|
}
|
|
m_inversePreconditioner->Mult(input, m_preconditionedDirection);
|
|
m_jacobian->Mult(m_preconditionedDirection, output);
|
|
}
|
|
|
|
const mfem::Operator &FixedRightPreconditionedOperator::GetJacobian() const noexcept {
|
|
return *m_jacobian;
|
|
}
|
|
|
|
const mfem::Solver &FixedRightPreconditionedOperator::GetInversePreconditioner() const noexcept {
|
|
return *m_inversePreconditioner;
|
|
}
|
|
|
|
void ResidualHistoryMonitor::Reset() {
|
|
mfem::IterativeSolverMonitor::Reset();
|
|
m_history.clear();
|
|
}
|
|
|
|
void ResidualHistoryMonitor::MonitorResidual(
|
|
const int iteration,
|
|
const double norm,
|
|
const mfem::Vector &,
|
|
const bool final
|
|
) {
|
|
m_history.push_back({.iteration = iteration, .reportedNorm = norm, .final = final});
|
|
}
|
|
|
|
const std::vector<IterationResidualMeasurement> &ResidualHistoryMonitor::GetHistory() const noexcept {
|
|
return m_history;
|
|
}
|
|
|
|
DirectResidualMeasurement measureDirectResidual(
|
|
const mfem::Operator &jacobian,
|
|
const mfem::Vector &rightHandSide,
|
|
const mfem::Vector &solution,
|
|
const std::span<const operators::RootBlockDescriptor> residualBlocks,
|
|
const MPI_Comm communicator,
|
|
const double denominatorFloor
|
|
) {
|
|
if (jacobian.Height() != jacobian.Width() || rightHandSide.Size() != jacobian.Height() ||
|
|
solution.Size() != jacobian.Width()) {
|
|
throw std::invalid_argument("Direct residual measurement received incompatible linear-system dimensions.");
|
|
}
|
|
if (!std::isfinite(denominatorFloor) || denominatorFloor <= 0.0) {
|
|
throw std::invalid_argument("The direct-residual denominator floor must be finite and positive.");
|
|
}
|
|
verify_finite_vector(rightHandSide, "Direct residual measurement received a non-finite right-hand side.");
|
|
verify_finite_vector(solution, "Direct residual measurement received a non-finite solution.");
|
|
|
|
int expectedOffset = 0;
|
|
for (const operators::RootBlockDescriptor &block : residualBlocks) {
|
|
if (block.kind != operators::RootBlockKind::residual || block.offset != expectedOffset || block.size < 0 ||
|
|
block.offset + block.size > jacobian.Height() || !std::isfinite(block.scale) || block.scale <= 0.0) {
|
|
throw std::invalid_argument("Residual block descriptors do not form the canonical equation layout.");
|
|
}
|
|
expectedOffset += block.size;
|
|
}
|
|
if (expectedOffset != jacobian.Height()) {
|
|
throw std::invalid_argument("Residual block descriptors do not cover the complete equation vector.");
|
|
}
|
|
|
|
mfem::Vector action(jacobian.Height());
|
|
jacobian.Mult(solution, action);
|
|
if (action.Size() != rightHandSide.Size()) {
|
|
throw std::runtime_error("The Jacobian returned an action with the wrong size.");
|
|
}
|
|
mfem::Vector trueResidual(rightHandSide);
|
|
trueResidual -= action;
|
|
verify_finite_vector(trueResidual, "Direct residual measurement produced a non-finite residual.");
|
|
|
|
DirectResidualMeasurement measurement;
|
|
measurement.rightHandSideNorm = global_norm(rightHandSide, communicator);
|
|
measurement.trueResidualNorm = global_norm(trueResidual, communicator);
|
|
const double denominator = std::max(measurement.rightHandSideNorm, denominatorFloor);
|
|
measurement.relativeResidual = measurement.trueResidualNorm / denominator;
|
|
measurement.blocks.reserve(residualBlocks.size());
|
|
|
|
for (const operators::RootBlockDescriptor &block : residualBlocks) {
|
|
const mfem::Vector blockRightHandSide(
|
|
const_cast<mfem::real_t *>(rightHandSide.GetData()) + block.offset, block.size
|
|
);
|
|
const mfem::Vector blockResidual(trueResidual.GetData() + block.offset, block.size);
|
|
const double blockRightHandSideNorm = global_norm(blockRightHandSide, communicator);
|
|
const double blockResidualNorm = global_norm(blockResidual, communicator);
|
|
const double blockDenominator = std::max(blockRightHandSideNorm, denominatorFloor);
|
|
const double globalResidualFraction =
|
|
measurement.trueResidualNorm > denominatorFloor
|
|
? blockResidualNorm * blockResidualNorm /
|
|
(measurement.trueResidualNorm * measurement.trueResidualNorm)
|
|
: 0.0;
|
|
measurement.blocks.push_back(
|
|
{.stableId = std::string(block.stableId),
|
|
.size = block.size,
|
|
.descriptorScale = block.scale,
|
|
.rightHandSideNorm = blockRightHandSideNorm,
|
|
.trueResidualNorm = blockResidualNorm,
|
|
.blockRelativeResidual = blockResidualNorm / blockDenominator,
|
|
.scaledRightHandSideNorm = blockRightHandSideNorm / block.scale,
|
|
.scaledTrueResidualNorm = blockResidualNorm / block.scale,
|
|
.contributionToGlobalRelativeResidual = blockResidualNorm / denominator,
|
|
.fractionOfGlobalSquaredResidualNorm = globalResidualFraction}
|
|
);
|
|
}
|
|
return measurement;
|
|
}
|
|
|
|
LinearSolveMeasurement measureLinearSolve(
|
|
const mfem::IterativeSolver &iterativeSolver,
|
|
const mfem::Operator &jacobian,
|
|
const mfem::Vector &rightHandSide,
|
|
const mfem::Vector &solution,
|
|
const std::span<const operators::RootBlockDescriptor> residualBlocks,
|
|
const OperatorApplicationStatistics &jacobianStatistics,
|
|
const OperatorApplicationStatistics &inversePreconditionerStatistics,
|
|
const PreconditionerLifecycleStatistics &inversePreconditionerLifecycle,
|
|
const ResidualHistoryMonitor &monitor,
|
|
const double localSolveSeconds,
|
|
const MPI_Comm communicator,
|
|
const double denominatorFloor
|
|
) {
|
|
if (!std::isfinite(localSolveSeconds) || localSolveSeconds < 0.0) {
|
|
throw std::invalid_argument("A linear-solve duration must be finite and nonnegative.");
|
|
}
|
|
|
|
const DirectResidualMeasurement directResidual =
|
|
measureDirectResidual(jacobian, rightHandSide, solution, residualBlocks, communicator, denominatorFloor);
|
|
const double reportedInitial = iterativeSolver.GetInitialNorm();
|
|
const double reportedFinal = iterativeSolver.GetFinalNorm();
|
|
const double reportedReduction =
|
|
std::abs(reportedInitial) > denominatorFloor ? std::abs(reportedFinal) / std::abs(reportedInitial) : 0.0;
|
|
double digitsPerJacobianApplication = 0.0;
|
|
if (jacobianStatistics.applications > 0 && directResidual.relativeResidual >= 0.0 &&
|
|
std::isfinite(directResidual.relativeResidual)) {
|
|
digitsPerJacobianApplication = -std::log10(std::max(directResidual.relativeResidual, denominatorFloor)) /
|
|
static_cast<double>(jacobianStatistics.applications);
|
|
}
|
|
|
|
return {
|
|
.solverConverged = iterativeSolver.GetConverged(),
|
|
.outerIterations = iterativeSolver.GetNumIterations(),
|
|
.solverReportedInitialNorm = reportedInitial,
|
|
.solverReportedFinalNorm = reportedFinal,
|
|
.solverReportedResidualReduction = reportedReduction,
|
|
.trueResidualDigitsReducedPerJacobianApplication = digitsPerJacobianApplication,
|
|
.solveSecondsMaximumRank = maximum_rank_value(localSolveSeconds, communicator),
|
|
.jacobian = maximum_rank_statistics(jacobianStatistics, communicator),
|
|
.inversePreconditioner = maximum_rank_statistics(inversePreconditionerStatistics, communicator),
|
|
.inversePreconditionerLifecycle =
|
|
maximum_rank_lifecycle_statistics(inversePreconditionerLifecycle, communicator),
|
|
.directResidual = directResidual,
|
|
.reportedResidualHistory = monitor.GetHistory()
|
|
};
|
|
}
|
|
|
|
ArnoldiSpectralMeasurement measureArnoldiSpectrum(
|
|
const mfem::Operator &operation,
|
|
const mfem::Vector &initialDirection,
|
|
const MPI_Comm communicator,
|
|
const ArnoldiOptions &options
|
|
) {
|
|
if (operation.Height() != operation.Width() || operation.Width() <= 0) {
|
|
throw std::invalid_argument("Arnoldi diagnostics require a nonempty square operator.");
|
|
}
|
|
if (initialDirection.Size() != operation.Width()) {
|
|
throw std::invalid_argument("The Arnoldi initial direction has the wrong size.");
|
|
}
|
|
if (options.krylovDimension <= 0 || !std::isfinite(options.breakdownRelativeTolerance) ||
|
|
options.breakdownRelativeTolerance < 0.0 || !std::isfinite(options.ritzConvergenceRelativeTolerance) ||
|
|
options.ritzConvergenceRelativeTolerance < 0.0) {
|
|
throw std::invalid_argument("Arnoldi diagnostic options are invalid.");
|
|
}
|
|
verify_finite_vector(initialDirection, "Arnoldi diagnostics received a non-finite initial direction.");
|
|
|
|
const Clock::time_point measurementStart = Clock::now();
|
|
OperatorApplicationStatistics localApplicationStatistics;
|
|
|
|
const double initialNorm = global_norm(initialDirection, communicator);
|
|
if (!std::isfinite(initialNorm) || initialNorm <= 0.0) {
|
|
throw std::invalid_argument("Arnoldi diagnostics require a nonzero initial direction.");
|
|
}
|
|
|
|
const int requestedDimension = std::min(options.krylovDimension, operation.Width());
|
|
Eigen::MatrixXd hessenberg = Eigen::MatrixXd::Zero(requestedDimension + 1, requestedDimension);
|
|
std::vector<mfem::Vector> basis;
|
|
basis.reserve(static_cast<std::size_t>(requestedDimension + 1));
|
|
basis.emplace_back(initialDirection);
|
|
basis.back() /= initialNorm;
|
|
|
|
int achievedDimension{0};
|
|
bool invariantSubspaceFound{false};
|
|
|
|
for (int column = 0; column < requestedDimension; ++column) {
|
|
mfem::Vector candidate(operation.Height());
|
|
const Clock::time_point applicationStart = Clock::now();
|
|
operation.Mult(basis[static_cast<std::size_t>(column)], candidate);
|
|
const double applicationSeconds = seconds_between(applicationStart, Clock::now());
|
|
++localApplicationStatistics.applications;
|
|
localApplicationStatistics.totalSeconds += applicationSeconds;
|
|
localApplicationStatistics.maximumSeconds =
|
|
std::max(localApplicationStatistics.maximumSeconds, applicationSeconds);
|
|
if (candidate.Size() != operation.Height()) {
|
|
throw std::runtime_error("The Arnoldi operator returned a vector with the wrong size.");
|
|
}
|
|
verify_finite_vector(candidate, "The Arnoldi operator produced a non-finite vector.");
|
|
const double unorthogonalizedNorm = global_norm(candidate, communicator);
|
|
|
|
const int passCount = options.reorthogonalize ? 2 : 1;
|
|
for (int pass = 0; pass < passCount; ++pass) {
|
|
for (int row = 0; row <= column; ++row) {
|
|
const double projection = global_dot(basis[static_cast<std::size_t>(row)], candidate, communicator);
|
|
hessenberg(row, column) += projection;
|
|
candidate.Add(-projection, basis[static_cast<std::size_t>(row)]);
|
|
}
|
|
}
|
|
|
|
const double nextNorm = global_norm(candidate, communicator);
|
|
hessenberg(column + 1, column) = nextNorm;
|
|
achievedDimension = column + 1;
|
|
const double breakdownScale = std::max(unorthogonalizedNorm, 1.0);
|
|
if (nextNorm <= options.breakdownRelativeTolerance * breakdownScale) {
|
|
invariantSubspaceFound = true;
|
|
break;
|
|
}
|
|
if (column + 1 < requestedDimension) {
|
|
candidate /= nextNorm;
|
|
basis.push_back(std::move(candidate));
|
|
}
|
|
}
|
|
|
|
if (achievedDimension <= 0) {
|
|
throw std::runtime_error("Arnoldi diagnostics did not construct a Krylov projection.");
|
|
}
|
|
|
|
const Eigen::MatrixXd projected = copy_hessenberg(hessenberg, achievedDimension, achievedDimension);
|
|
const Eigen::MatrixXd projectedRectangular =
|
|
copy_hessenberg(hessenberg, achievedDimension + 1, achievedDimension);
|
|
|
|
Eigen::EigenSolver<Eigen::MatrixXd> eigenSolver(projected, true);
|
|
if (eigenSolver.info() != Eigen::Success) {
|
|
throw std::runtime_error("The projected Arnoldi eigenproblem did not converge.");
|
|
}
|
|
Eigen::JacobiSVD<Eigen::MatrixXd> singularValueDecomposition(projectedRectangular);
|
|
if (singularValueDecomposition.info() != Eigen::Success) {
|
|
throw std::runtime_error("The projected Arnoldi singular-value problem did not converge.");
|
|
}
|
|
|
|
ArnoldiSpectralMeasurement measurement;
|
|
measurement.requestedDimension = requestedDimension;
|
|
measurement.achievedDimension = achievedDimension;
|
|
measurement.invariantSubspaceFound = invariantSubspaceFound;
|
|
const OperatorApplicationStatistics globalApplicationStatistics =
|
|
maximum_rank_statistics(localApplicationStatistics, communicator);
|
|
measurement.operatorApplications = globalApplicationStatistics.applications;
|
|
measurement.operatorApplicationSecondsMaximumRank = globalApplicationStatistics.totalSeconds;
|
|
measurement.operatorMaximumApplicationSecondsMaximumRank = globalApplicationStatistics.maximumSeconds;
|
|
measurement.ritzValues.reserve(static_cast<std::size_t>(achievedDimension));
|
|
|
|
const Eigen::VectorXd singularValues = singularValueDecomposition.singularValues();
|
|
measurement.projectedLargestSingularValue = singularValues(0);
|
|
measurement.projectedSmallestSingularValue = singularValues(singularValues.size() - 1);
|
|
measurement.projectedConditionProxy =
|
|
measurement.projectedSmallestSingularValue > 0.0
|
|
? measurement.projectedLargestSingularValue / measurement.projectedSmallestSingularValue
|
|
: std::numeric_limits<double>::infinity();
|
|
|
|
const double finalSubdiagonal = hessenberg(achievedDimension, achievedDimension - 1);
|
|
std::complex<double> centroid{0.0, 0.0};
|
|
const auto eigenvalues = eigenSolver.eigenvalues();
|
|
const auto eigenvectors = eigenSolver.eigenvectors();
|
|
for (int index = 0; index < achievedDimension; ++index) {
|
|
const std::complex<double> eigenvalue = eigenvalues(index);
|
|
const double eigenvectorNorm = eigenvectors.col(index).norm();
|
|
const double residualEstimate =
|
|
eigenvectorNorm > 0.0
|
|
? std::abs(finalSubdiagonal * eigenvectors(achievedDimension - 1, index)) / eigenvectorNorm
|
|
: std::numeric_limits<double>::infinity();
|
|
const double convergenceScale = std::max(std::abs(eigenvalue), 1.0);
|
|
const double relativeResidualEstimate = residualEstimate / convergenceScale;
|
|
const bool converged = relativeResidualEstimate <= options.ritzConvergenceRelativeTolerance;
|
|
|
|
measurement.ritzValues.push_back(
|
|
{.realPart = eigenvalue.real(),
|
|
.imaginaryPart = eigenvalue.imag(),
|
|
.magnitude = std::abs(eigenvalue),
|
|
.distanceFromOne = std::abs(eigenvalue - std::complex<double>{1.0, 0.0}),
|
|
.residualEstimate = residualEstimate,
|
|
.relativeResidualEstimate = relativeResidualEstimate,
|
|
.converged = converged}
|
|
);
|
|
centroid += eigenvalue;
|
|
measurement.convergedRitzValueCount += converged ? 1 : 0;
|
|
measurement.negativeRealPartCount += eigenvalue.real() < 0.0 ? 1 : 0;
|
|
}
|
|
centroid /= static_cast<double>(achievedDimension);
|
|
measurement.centroidRealPart = centroid.real();
|
|
measurement.centroidImaginaryPart = centroid.imag();
|
|
|
|
measurement.minimumMagnitude = std::numeric_limits<double>::infinity();
|
|
measurement.minimumRealPart = std::numeric_limits<double>::infinity();
|
|
measurement.maximumRealPart = -std::numeric_limits<double>::infinity();
|
|
double squaredDistanceFromOne{0.0};
|
|
double squaredClusterRadius{0.0};
|
|
for (const RitzValueMeasurement &ritz : measurement.ritzValues) {
|
|
const std::complex<double> value{ritz.realPart, ritz.imaginaryPart};
|
|
measurement.minimumMagnitude = std::min(measurement.minimumMagnitude, ritz.magnitude);
|
|
measurement.maximumMagnitude = std::max(measurement.maximumMagnitude, ritz.magnitude);
|
|
measurement.minimumRealPart = std::min(measurement.minimumRealPart, ritz.realPart);
|
|
measurement.maximumRealPart = std::max(measurement.maximumRealPart, ritz.realPart);
|
|
measurement.maximumAbsoluteImaginaryPart =
|
|
std::max(measurement.maximumAbsoluteImaginaryPart, std::abs(ritz.imaginaryPart));
|
|
squaredDistanceFromOne += ritz.distanceFromOne * ritz.distanceFromOne;
|
|
squaredClusterRadius += std::norm(value - centroid);
|
|
|
|
double pairDefect = std::numeric_limits<double>::infinity();
|
|
for (const RitzValueMeasurement &candidate : measurement.ritzValues) {
|
|
pairDefect = std::min(
|
|
pairDefect,
|
|
std::abs(std::complex<double>{candidate.realPart, candidate.imaginaryPart} - std::conj(value))
|
|
);
|
|
}
|
|
measurement.conjugatePairDefect = std::max(measurement.conjugatePairDefect, pairDefect);
|
|
}
|
|
measurement.rmsDistanceFromOne = std::sqrt(squaredDistanceFromOne / achievedDimension);
|
|
measurement.rmsClusterRadius = std::sqrt(squaredClusterRadius / achievedDimension);
|
|
|
|
const double projectedFrobeniusSquared = projected.squaredNorm();
|
|
if (projectedFrobeniusSquared > 0.0) {
|
|
const Eigen::MatrixXd normalityCommutator =
|
|
projected.transpose() * projected - projected * projected.transpose();
|
|
measurement.projectedDepartureFromNormality = normalityCommutator.norm() / projectedFrobeniusSquared;
|
|
}
|
|
|
|
const Eigen::MatrixXd hermitianPart = 0.5 * (projected + projected.transpose());
|
|
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> fieldOfValuesSolver(hermitianPart);
|
|
if (fieldOfValuesSolver.info() != Eigen::Success) {
|
|
throw std::runtime_error("The projected field-of-values problem did not converge.");
|
|
}
|
|
measurement.projectedFieldOfValuesMinimumRealPart = fieldOfValuesSolver.eigenvalues().minCoeff();
|
|
measurement.projectedFieldOfValuesMaximumRealPart = fieldOfValuesSolver.eigenvalues().maxCoeff();
|
|
const double localMeasurementSeconds = seconds_between(measurementStart, Clock::now());
|
|
const double localNonApplicationSeconds =
|
|
std::max(localMeasurementSeconds - localApplicationStatistics.totalSeconds, 0.0);
|
|
measurement.measurementSecondsMaximumRank = maximum_rank_value(localMeasurementSeconds, communicator);
|
|
measurement.nonApplicationSecondsMaximumRank = maximum_rank_value(localNonApplicationSeconds, communicator);
|
|
return measurement;
|
|
}
|
|
|
|
std::vector<RitzValueMeasurement> selectRitzValues(
|
|
const ArnoldiSpectralMeasurement &measurement,
|
|
const RitzValueOrdering ordering,
|
|
const int count
|
|
) {
|
|
if (count < 0) {
|
|
throw std::invalid_argument("The requested Ritz-value count must be nonnegative.");
|
|
}
|
|
|
|
std::vector<RitzValueMeasurement> selected;
|
|
selected.reserve(measurement.ritzValues.size());
|
|
for (const RitzValueMeasurement &value : measurement.ritzValues) {
|
|
if (value.converged) {
|
|
selected.push_back(value);
|
|
}
|
|
}
|
|
|
|
std::ranges::sort(selected, [ordering](const RitzValueMeasurement &left, const RitzValueMeasurement &right) {
|
|
switch (ordering) {
|
|
case RitzValueOrdering::closest_to_zero:
|
|
return left.magnitude < right.magnitude;
|
|
case RitzValueOrdering::farthest_from_one:
|
|
return left.distanceFromOne > right.distanceFromOne;
|
|
case RitzValueOrdering::smallest_real_part:
|
|
return left.realPart < right.realPart;
|
|
case RitzValueOrdering::largest_magnitude:
|
|
return left.magnitude > right.magnitude;
|
|
}
|
|
return false;
|
|
});
|
|
if (static_cast<int>(selected.size()) > count) {
|
|
selected.resize(static_cast<std::size_t>(count));
|
|
}
|
|
return selected;
|
|
}
|
|
} // namespace mean_field::solver
|