KOKO\Mimi
1269 lines
54 KiB
C++
1269 lines
54 KiB
C++
#include "fesa/solvers/eigen/mkl_feast_eigen_solver.h"
|
|
|
|
#include <mkl.h>
|
|
|
|
#include <algorithm>
|
|
#include <array>
|
|
#include <cmath>
|
|
#include <cstddef>
|
|
#include <cstdint>
|
|
#include <iomanip>
|
|
#include <limits>
|
|
#include <new>
|
|
#include <numeric>
|
|
#include <sstream>
|
|
#include <stdexcept>
|
|
#include <string>
|
|
#include <utility>
|
|
#include <vector>
|
|
|
|
#include "fesa/math/sparse_matrix.h"
|
|
#include "fesa/fem/dof_manager.h"
|
|
#include "fesa/math/vector.h"
|
|
#include "fesa/solvers/linear/mkl_pardiso_solver.h"
|
|
#include "solvers/eigen/mkl_feast_eigen_solver_test_seam.h"
|
|
#include "math/dense_blas_internal.h"
|
|
|
|
namespace fesa {
|
|
|
|
namespace internal {
|
|
namespace {
|
|
constexpr std::size_t kMaximumCountCalls = 512U;
|
|
thread_local FeastBackendTestHook* g_feast_backend_test_hook = nullptr;
|
|
thread_local FeastWorkspaceTestHook* g_feast_workspace_test_hook = nullptr;
|
|
}
|
|
|
|
bool CanCallFeast(const std::size_t used) {
|
|
return used < kMaximumCountCalls;
|
|
}
|
|
|
|
bool CanStartUpperCandidate(const std::size_t used) {
|
|
return used <= 384U - 21U;
|
|
}
|
|
|
|
bool ProductFits(const std::size_t left, const std::size_t right) {
|
|
return left == 0U ||
|
|
right <= (std::numeric_limits<std::size_t>::max)() / left;
|
|
}
|
|
|
|
bool FitsDoubleStorage(const std::size_t count) {
|
|
return count <= std::vector<double>{}.max_size() &&
|
|
count <= static_cast<std::size_t>(
|
|
(std::numeric_limits<std::ptrdiff_t>::max)()) /
|
|
sizeof(double);
|
|
}
|
|
|
|
ScopedFeastBackendTestHook::ScopedFeastBackendTestHook(
|
|
FeastBackendTestHook hook)
|
|
: hook_(std::move(hook)), previous_(g_feast_backend_test_hook) {
|
|
g_feast_backend_test_hook = &hook_;
|
|
}
|
|
|
|
ScopedFeastBackendTestHook::~ScopedFeastBackendTestHook() {
|
|
g_feast_backend_test_hook = previous_;
|
|
}
|
|
|
|
ScopedFeastWorkspaceTestHook::ScopedFeastWorkspaceTestHook(
|
|
FeastWorkspaceTestHook hook)
|
|
: hook_(std::move(hook)), previous_(g_feast_workspace_test_hook) {
|
|
g_feast_workspace_test_hook = &hook_;
|
|
}
|
|
|
|
ScopedFeastWorkspaceTestHook::~ScopedFeastWorkspaceTestHook() {
|
|
g_feast_workspace_test_hook = previous_;
|
|
}
|
|
|
|
void ApplyFeastWorkspaceTestHook(const bool count_workspace) {
|
|
if (g_feast_workspace_test_hook != nullptr)
|
|
(*g_feast_workspace_test_hook)(count_workspace);
|
|
}
|
|
|
|
void ApplyFeastBackendTestHook(const FeastBackendTestCall& call, MKL_INT* m0,
|
|
MKL_INT* returned, MKL_INT* loop,
|
|
double* epsout, MKL_INT* info,
|
|
double* eigenvalues,
|
|
const std::size_t eigenvalue_count,
|
|
double* vectors,
|
|
const std::size_t vector_count,
|
|
double* residuals,
|
|
const std::size_t residual_count) {
|
|
if (g_feast_backend_test_hook == nullptr) return;
|
|
FeastBackendTestOutputs outputs{
|
|
static_cast<int>(*info), static_cast<std::ptrdiff_t>(*m0),
|
|
static_cast<std::ptrdiff_t>(*returned), static_cast<int>(*loop), *epsout,
|
|
std::vector<double>(eigenvalues, eigenvalues + eigenvalue_count),
|
|
std::vector<double>(vectors, vectors + vector_count),
|
|
std::vector<double>(residuals, residuals + residual_count)};
|
|
(*g_feast_backend_test_hook)(call, outputs);
|
|
*info = static_cast<MKL_INT>(outputs.info);
|
|
*m0 = static_cast<MKL_INT>(outputs.m0);
|
|
*returned = static_cast<MKL_INT>(outputs.returned);
|
|
*loop = static_cast<MKL_INT>(outputs.loop);
|
|
*epsout = outputs.epsout;
|
|
if (outputs.eigenvalues.size() == eigenvalue_count)
|
|
std::copy(outputs.eigenvalues.begin(), outputs.eigenvalues.end(), eigenvalues);
|
|
if (outputs.vectors.size() == vector_count)
|
|
std::copy(outputs.vectors.begin(), outputs.vectors.end(), vectors);
|
|
if (outputs.residuals.size() == residual_count)
|
|
std::copy(outputs.residuals.begin(), outputs.residuals.end(), residuals);
|
|
}
|
|
} // namespace internal
|
|
|
|
namespace {
|
|
|
|
constexpr std::uint64_t kSeed = 0x4D495443344D4645ULL;
|
|
|
|
Status Failure(const std::string& code, const std::string& identity,
|
|
const std::string& message) {
|
|
return Status::Failure(
|
|
FailureCategory::kSolver,
|
|
{{Severity::kError, code, {{}, 0U}, "FEAST", identity, message}});
|
|
}
|
|
|
|
bool FitsMkl(const std::size_t value) {
|
|
return value <=
|
|
static_cast<std::size_t>((std::numeric_limits<MKL_INT>::max)());
|
|
}
|
|
|
|
double NextUp(const double value) {
|
|
return std::nextafter(value, (std::numeric_limits<double>::infinity)());
|
|
}
|
|
|
|
double Delta(const double value, const double beta) {
|
|
const double up =
|
|
std::nextafter(value, (std::numeric_limits<double>::infinity)());
|
|
const double down =
|
|
std::nextafter(value, -(std::numeric_limits<double>::infinity)());
|
|
const double ulp = (std::max)(std::abs(up - value), std::abs(value - down));
|
|
return 1.0e-8 * beta + 64.0 * ulp;
|
|
}
|
|
|
|
/// @brief Restores the calling thread's oneMKL local-thread setting.
|
|
class LocalThreadGuard {
|
|
public:
|
|
LocalThreadGuard() : previous_(mkl_set_num_threads_local(1)) {}
|
|
~LocalThreadGuard() {
|
|
static_cast<void>(mkl_set_num_threads_local(previous_));
|
|
}
|
|
|
|
private:
|
|
int previous_;
|
|
};
|
|
|
|
struct Csr {
|
|
std::vector<MKL_INT> offsets;
|
|
std::vector<MKL_INT> columns;
|
|
std::vector<double> values;
|
|
};
|
|
|
|
Result<Csr> CopyUpperOneBased(const SparseMatrix& matrix, const double scale) {
|
|
Csr copy;
|
|
copy.offsets.reserve(matrix.Rows() + 1U);
|
|
copy.offsets.push_back(1);
|
|
const auto& offsets = matrix.RowOffsets();
|
|
const auto& columns = matrix.ColumnIndices();
|
|
const auto& values = matrix.Values();
|
|
for (std::size_t row = 0U; row < matrix.Rows(); ++row) {
|
|
for (std::size_t position = offsets[row]; position < offsets[row + 1U];
|
|
++position) {
|
|
const std::size_t column = columns[position];
|
|
if (column < row) {
|
|
continue;
|
|
}
|
|
if (!FitsMkl(column + 1U) || !FitsMkl(copy.values.size() + 1U)) {
|
|
return Result<Csr>::Failure(Failure(
|
|
"feast-integer-overflow", "csr-index",
|
|
"The public CSR index exceeds the oneMKL ILP64 integer range."));
|
|
}
|
|
copy.columns.push_back(static_cast<MKL_INT>(column + 1U));
|
|
copy.values.push_back(values[position] / scale);
|
|
}
|
|
if (!FitsMkl(copy.values.size() + 1U)) {
|
|
return Result<Csr>::Failure(Failure("feast-integer-overflow",
|
|
"csr-offset",
|
|
"The converted one-based CSR offset "
|
|
"exceeds the oneMKL integer range."));
|
|
}
|
|
copy.offsets.push_back(static_cast<MKL_INT>(copy.values.size() + 1U));
|
|
}
|
|
return Result<Csr>::Success(std::move(copy));
|
|
}
|
|
|
|
/// @brief Keeps the independently scaled full or upper-mirrored operator.
|
|
Result<SparseMatrix> ScaledActionMatrix(const SparseMatrix& matrix,
|
|
const double scale,
|
|
const bool mirror_upper = false) {
|
|
std::vector<CooContribution> contributions;
|
|
std::vector<std::vector<std::size_t>> columns(matrix.Rows());
|
|
for (std::size_t row = 0; row < matrix.Rows(); ++row) {
|
|
for (auto position = matrix.RowOffsets()[row];
|
|
position < matrix.RowOffsets()[row + 1U]; ++position) {
|
|
const auto column = matrix.ColumnIndices()[position];
|
|
if (mirror_upper && column < row) continue;
|
|
const auto value = matrix.Values()[position] / scale;
|
|
columns[row].push_back(column);
|
|
contributions.push_back({row, column, value, row, position});
|
|
if (mirror_upper && column != row) {
|
|
columns[column].push_back(row);
|
|
contributions.push_back({column, row, value, row, position});
|
|
}
|
|
}
|
|
}
|
|
SparsePattern pattern{{0U}, {}};
|
|
for (auto& row : columns) {
|
|
std::sort(row.begin(), row.end());
|
|
pattern.column_indices.insert(pattern.column_indices.end(), row.begin(),
|
|
row.end());
|
|
pattern.row_offsets.push_back(pattern.column_indices.size());
|
|
}
|
|
return SparseMatrix::FromCoo(matrix.Rows(), matrix.Columns(),
|
|
std::move(contributions), pattern);
|
|
}
|
|
|
|
Result<Csr> IdentityOneBased(const std::size_t size) {
|
|
if (!FitsMkl(size) || !FitsMkl(size + 1U)) {
|
|
return Result<Csr>::Failure(
|
|
Failure("feast-integer-overflow", "identity",
|
|
"The identity CSR exceeds oneMKL range."));
|
|
}
|
|
Csr identity;
|
|
identity.offsets.reserve(size + 1U);
|
|
identity.columns.reserve(size);
|
|
identity.values.assign(size, 1.0);
|
|
for (std::size_t index = 0U; index <= size; ++index) {
|
|
identity.offsets.push_back(static_cast<MKL_INT>(index + 1U));
|
|
if (index < size) {
|
|
identity.columns.push_back(static_cast<MKL_INT>(index + 1U));
|
|
}
|
|
}
|
|
return Result<Csr>::Success(std::move(identity));
|
|
}
|
|
|
|
Status ValidateSymmetry(const SparseMatrix& matrix) {
|
|
const auto& offsets = matrix.RowOffsets();
|
|
const auto& columns = matrix.ColumnIndices();
|
|
const auto& values = matrix.Values();
|
|
double scale = 0.0;
|
|
for (const double value : values) scale = (std::max)(scale, std::abs(value));
|
|
double norm = 0.0;
|
|
double asymmetry = 0.0;
|
|
for (std::size_t row = 0U; row < matrix.Rows(); ++row) {
|
|
for (std::size_t p = offsets[row]; p < offsets[row + 1U]; ++p) {
|
|
const std::size_t column = columns[p];
|
|
const auto begin =
|
|
columns.begin() + static_cast<std::ptrdiff_t>(offsets[column]);
|
|
const auto end =
|
|
columns.begin() + static_cast<std::ptrdiff_t>(offsets[column + 1U]);
|
|
const auto found = std::lower_bound(begin, end, row);
|
|
if (found == end || *found != row) {
|
|
return Failure("feast-nonsymmetric-csr",
|
|
std::to_string(row) + ":" + std::to_string(column),
|
|
"FEAST requires a full symmetric public CSR before "
|
|
"upper extraction.");
|
|
}
|
|
const std::size_t q = static_cast<std::size_t>(found - columns.begin());
|
|
if (scale != 0.0) {
|
|
norm = std::hypot(norm, values[p] / scale);
|
|
asymmetry =
|
|
std::hypot(asymmetry, values[p] / scale - values[q] / scale);
|
|
}
|
|
}
|
|
}
|
|
if (norm > 0.0 && asymmetry / norm > 1e-12)
|
|
return Failure("feast-nonsymmetric-csr", "whitened-operator",
|
|
"The public CSR violates the Frobenius symmetry contract.");
|
|
return Status::Ok();
|
|
}
|
|
|
|
void SetControls(std::array<MKL_INT, 128>& fpm, const bool count) {
|
|
feastinit(fpm.data());
|
|
fpm[0] = 0;
|
|
fpm[1] = 16;
|
|
fpm[2] = 12;
|
|
fpm[3] = 20;
|
|
fpm[4] = 1;
|
|
fpm[5] = 0; // Trace convergence alone can leave individual pairs inaccurate.
|
|
fpm[13] = 0;
|
|
fpm[26] = 1;
|
|
fpm[27] = 1;
|
|
fpm[49] = count ? 1 : 0;
|
|
fpm[63] = 0;
|
|
}
|
|
|
|
std::string FormatDouble(const double value) {
|
|
std::ostringstream stream;
|
|
stream << std::scientific << std::setprecision(17) << value;
|
|
return stream.str();
|
|
}
|
|
|
|
std::string CallHistory(const std::vector<FeastBackendCall>& calls) {
|
|
std::string history;
|
|
for (const FeastBackendCall& call : calls) {
|
|
if (!history.empty()) history += ";";
|
|
history += std::to_string(call.ordinal) + ":" + call.kind +
|
|
":info=" + std::to_string(call.info) +
|
|
":m0=" + std::to_string(call.allocated_m0) + ":(" +
|
|
FormatDouble(call.lower) + "," + FormatDouble(call.upper) + ")";
|
|
}
|
|
return history;
|
|
}
|
|
|
|
bool FillInitialBasis(std::vector<double>& vectors, const std::size_t n,
|
|
const std::size_t columns) {
|
|
std::fill(vectors.begin(), vectors.end(), 0.0);
|
|
if (columns == n) {
|
|
for (std::size_t i = 0U; i < n; ++i) vectors[i + i * n] = 1.0;
|
|
return true;
|
|
}
|
|
std::uint64_t state = kSeed;
|
|
for (double& value : vectors) {
|
|
state = 6364136223846793005ULL * state + 1442695040888963407ULL;
|
|
value = 2.0 * static_cast<double>(state >> 11U) * 0x1.0p-53 - 1.0;
|
|
}
|
|
std::vector<std::vector<double>> reflectors;
|
|
std::vector<double> signs;
|
|
reflectors.reserve(columns);
|
|
signs.reserve(columns);
|
|
for (std::size_t column = 0U; column < columns; ++column) {
|
|
double* current = vectors.data() + column * n;
|
|
double norm = 0.0;
|
|
double column_scale = 0.0;
|
|
for (std::size_t row = 0U; row < n; ++row) {
|
|
column_scale = std::hypot(column_scale, current[row]);
|
|
}
|
|
for (std::size_t row = column; row < n; ++row) {
|
|
norm = std::hypot(norm, current[row]);
|
|
}
|
|
if (!(column_scale > 0.0) || !(norm > 1.0e-12 * column_scale) ||
|
|
!std::isfinite(norm)) {
|
|
return false;
|
|
}
|
|
// Avoid cancellation in the reflector. Adjust Q's column sign below so
|
|
// R has a positive diagonal without subtracting nearly equal numbers.
|
|
const double alpha = -std::copysign(norm, current[column]);
|
|
signs.push_back(std::copysign(1.0, alpha));
|
|
std::vector<double> reflector(n - column, 0.0);
|
|
reflector[0U] = current[column] - alpha;
|
|
for (std::size_t row = column + 1U; row < n; ++row) {
|
|
reflector[row - column] = current[row];
|
|
}
|
|
double reflector_norm = 0.0;
|
|
for (const double value : reflector)
|
|
reflector_norm = std::hypot(reflector_norm, value);
|
|
if (!(reflector_norm > 0.0) || !std::isfinite(reflector_norm)) {
|
|
return false;
|
|
}
|
|
for (double& value : reflector) value /= reflector_norm;
|
|
for (std::size_t target = column; target < columns; ++target) {
|
|
const std::size_t length = n - column;
|
|
auto* target_column = vectors.data() + column + target * n;
|
|
const double dot = dense_blas_internal::DotValues(
|
|
reflector.data(), length, target_column);
|
|
dense_blas_internal::AxpyValues(-2.0 * dot, reflector.data(), length,
|
|
target_column);
|
|
}
|
|
reflectors.push_back(std::move(reflector));
|
|
}
|
|
std::fill(vectors.begin(), vectors.end(), 0.0);
|
|
for (std::size_t column = 0U; column < columns; ++column)
|
|
vectors[column + column * n] = 1.0;
|
|
for (std::size_t index = reflectors.size(); index > 0U; --index) {
|
|
const std::size_t column = index - 1U;
|
|
const auto& reflector = reflectors[column];
|
|
for (std::size_t target = 0U; target < columns; ++target) {
|
|
const std::size_t length = n - column;
|
|
auto* target_column = vectors.data() + column + target * n;
|
|
const double dot = dense_blas_internal::DotValues(
|
|
reflector.data(), length, target_column);
|
|
dense_blas_internal::AxpyValues(-2.0 * dot, reflector.data(), length,
|
|
target_column);
|
|
}
|
|
}
|
|
for (std::size_t column = 0U; column < columns; ++column)
|
|
for (std::size_t row = 0U; row < n; ++row)
|
|
vectors[row + column * n] *= signs[column];
|
|
double gram_error = 0.0;
|
|
for (std::size_t left = 0U; left < columns; ++left) {
|
|
for (std::size_t right = 0U; right < columns; ++right) {
|
|
const double dot = dense_blas_internal::DotValues(
|
|
vectors.data() + left * n, n, vectors.data() + right * n);
|
|
gram_error = std::hypot(gram_error, dot - (left == right ? 1.0 : 0.0));
|
|
}
|
|
}
|
|
return std::isfinite(gram_error) &&
|
|
gram_error / std::sqrt(static_cast<double>(columns)) <= 1e-12;
|
|
}
|
|
|
|
} // namespace
|
|
|
|
Result<FeastEigenResult> MklFeastEigenSolver::SolveGeneralized(
|
|
const SparseMatrix& stiffness, const SparseMatrix& mass,
|
|
const std::size_t requested_modes, const double mass_lower_bound,
|
|
const double classification_scale) const {
|
|
try {
|
|
return SolveImpl(stiffness, &mass, requested_modes, classification_scale,
|
|
mass_lower_bound);
|
|
} catch (const std::bad_alloc&) {
|
|
return Result<FeastEigenResult>::Failure(
|
|
Failure("eigen-allocation-failed", "generalized-workspace",
|
|
"Allocation of generalized sparse-pencil storage failed."));
|
|
}
|
|
}
|
|
|
|
Result<FeastEigenResult> MklFeastEigenSolver::Solve(
|
|
const SparseMatrix& whitened_operator, const std::size_t requested_modes,
|
|
const double zero_operator_scale) const {
|
|
return SolveImpl(whitened_operator, nullptr, requested_modes,
|
|
zero_operator_scale, 1.0);
|
|
}
|
|
|
|
Result<FeastEigenResult> MklFeastEigenSolver::SolveImpl(
|
|
const SparseMatrix& whitened_operator, const SparseMatrix* mass,
|
|
const std::size_t requested_modes, const double zero_operator_scale,
|
|
const double mass_lower_bound) const try {
|
|
static_assert(sizeof(MKL_INT) == 8U,
|
|
"FEAST requires the installed ILP64 ABI.");
|
|
const Status csr_status = whitened_operator.Validate();
|
|
if (!csr_status.IsOk())
|
|
return Result<FeastEigenResult>::Failure(
|
|
Failure("feast-invalid-csr", "public-csr",
|
|
"The public sparse matrix failed CSR validation."));
|
|
const std::size_t n = whitened_operator.Rows();
|
|
if (mass != nullptr) {
|
|
if (!mass->Validate().IsOk() || mass->Rows() != n || mass->Columns() != n ||
|
|
!std::isfinite(mass_lower_bound) || !(mass_lower_bound > 0.0) ||
|
|
!std::isfinite(zero_operator_scale) || !(zero_operator_scale > 0.0)) {
|
|
return Result<FeastEigenResult>::Failure(Failure(
|
|
"feast-invalid-mass-bound", "generalized-pencil",
|
|
"A matching finite mass and certified positive bound are required."));
|
|
}
|
|
const Status mass_symmetry = ValidateSymmetry(*mass);
|
|
if (!mass_symmetry.IsOk())
|
|
return Result<FeastEigenResult>::Failure(mass_symmetry);
|
|
MklPardisoSolver spd_check;
|
|
const Status spd = spd_check.Factorize(*mass);
|
|
if (!spd.IsOk())
|
|
return Result<FeastEigenResult>::Failure(
|
|
Failure("feast-mass-not-positive-definite", "mass-spd",
|
|
"Sparse SPD mass factorization failed."));
|
|
}
|
|
if (n == 0U || whitened_operator.Columns() != n || requested_modes == 0U ||
|
|
requested_modes > n) {
|
|
return Result<FeastEigenResult>::Failure(Failure(
|
|
"feast-invalid-dimension", "matrix-or-request",
|
|
"FEAST requires 1 <= requested modes <= square matrix dimension."));
|
|
}
|
|
if (!FitsMkl(n) || !internal::ProductFits(n, n)) {
|
|
return Result<FeastEigenResult>::Failure(Failure(
|
|
"feast-integer-overflow", "matrix-shape",
|
|
"FEAST workspace dimensions exceed the private backend range."));
|
|
}
|
|
if (!internal::FitsDoubleStorage(n) || !internal::FitsDoubleStorage(n * n)) {
|
|
return Result<FeastEigenResult>::Failure(Failure(
|
|
"eigen-allocation-failed", "count-workspace",
|
|
"The required count workspace exceeds addressable double storage."));
|
|
}
|
|
const MKL_INT mkl_n = static_cast<MKL_INT>(n);
|
|
const Status symmetry = ValidateSymmetry(whitened_operator);
|
|
if (!symmetry.IsOk()) return Result<FeastEigenResult>::Failure(symmetry);
|
|
|
|
double infinity_norm = 0.0;
|
|
const auto& offsets = whitened_operator.RowOffsets();
|
|
const auto& values = whitened_operator.Values();
|
|
for (std::size_t row = 0U; row < n; ++row) {
|
|
double sum = 0.0;
|
|
for (std::size_t p = offsets[row]; p < offsets[row + 1U]; ++p) {
|
|
sum = mass != nullptr ? NextUp(sum + std::abs(values[p]))
|
|
: sum + std::abs(values[p]);
|
|
}
|
|
if (!std::isfinite(sum)) {
|
|
return Result<FeastEigenResult>::Failure(Failure(
|
|
"feast-nonfinite-operator-scale", "operator-scale",
|
|
"The operator infinity norm overflowed while computing its scale."));
|
|
}
|
|
infinity_norm = (std::max)(infinity_norm, sum);
|
|
}
|
|
const bool zero_operator = infinity_norm == 0.0;
|
|
if (zero_operator &&
|
|
(!(zero_operator_scale > 0.0) || !std::isfinite(zero_operator_scale))) {
|
|
return Result<FeastEigenResult>::Failure(Failure(
|
|
"feast-invalid-zero-scale", "zero-operator-scale",
|
|
"An exactly zero operator requires a finite positive physical scale."));
|
|
}
|
|
const double operator_scale =
|
|
mass != nullptr || zero_operator ? zero_operator_scale : infinity_norm;
|
|
if (!std::isfinite(operator_scale) || !(operator_scale > 0.0)) {
|
|
return Result<FeastEigenResult>::Failure(
|
|
Failure("feast-nonfinite-operator-scale", "operator-scale",
|
|
"The operator infinity norm is not a finite positive scale."));
|
|
}
|
|
auto a = CopyUpperOneBased(whitened_operator, operator_scale);
|
|
if (!a.HasValue()) return Result<FeastEigenResult>::Failure(a.GetStatus());
|
|
auto b =
|
|
mass != nullptr ? CopyUpperOneBased(*mass, 1.0) : IdentityOneBased(n);
|
|
if (!b.HasValue()) return Result<FeastEigenResult>::Failure(b.GetStatus());
|
|
auto action_matrix = ScaledActionMatrix(whitened_operator, operator_scale);
|
|
if (!action_matrix.HasValue())
|
|
return Result<FeastEigenResult>::Failure(action_matrix.GetStatus());
|
|
|
|
double lower = -1.0;
|
|
double upper = 1.0;
|
|
double g_lower = 0.0;
|
|
double g_upper = 0.0;
|
|
double beta = 1.0;
|
|
if (mass != nullptr) {
|
|
// ||K||2 <= ||K||inf and M >= m_lower*I bound the genuine pencil.
|
|
const double bound =
|
|
NextUp(NextUp(infinity_norm / mass_lower_bound) / operator_scale);
|
|
g_lower = -bound;
|
|
g_upper = bound;
|
|
beta = (std::max)(1.0, bound);
|
|
lower = std::nextafter(g_lower - 16.0 * Delta(g_lower, beta),
|
|
-(std::numeric_limits<double>::infinity)());
|
|
upper = NextUp(g_upper + 16.0 * Delta(g_upper, beta));
|
|
} else if (!zero_operator) {
|
|
g_lower = (std::numeric_limits<double>::infinity)();
|
|
g_upper = -(std::numeric_limits<double>::infinity)();
|
|
const auto& columns = whitened_operator.ColumnIndices();
|
|
for (std::size_t row = 0U; row < n; ++row) {
|
|
double off_diagonal = 0.0;
|
|
double diagonal = 0.0;
|
|
for (std::size_t p = offsets[row]; p < offsets[row + 1U]; ++p) {
|
|
const double value = values[p] / operator_scale;
|
|
if (columns[p] == row)
|
|
diagonal = value;
|
|
else
|
|
off_diagonal = NextUp(off_diagonal + std::abs(value));
|
|
}
|
|
g_lower =
|
|
(std::min)(g_lower, std::nextafter(
|
|
diagonal - off_diagonal,
|
|
-(std::numeric_limits<double>::infinity)()));
|
|
g_upper = (std::max)(g_upper, NextUp(diagonal + off_diagonal));
|
|
beta = (std::max)(beta, NextUp(std::abs(diagonal) + off_diagonal));
|
|
}
|
|
const double lower_delta = Delta(g_lower, beta);
|
|
const double upper_delta = Delta(g_upper, beta);
|
|
lower = std::nextafter(g_lower - 16.0 * lower_delta,
|
|
-(std::numeric_limits<double>::infinity)());
|
|
upper = NextUp(g_upper + 16.0 * upper_delta);
|
|
}
|
|
if (!std::isfinite(lower) || !std::isfinite(upper) || !(lower < upper)) {
|
|
return Result<FeastEigenResult>::Failure(
|
|
Failure("feast-invalid-enclosure", "gershgorin",
|
|
"A strict finite FEAST enclosure could not be constructed."));
|
|
}
|
|
|
|
LocalThreadGuard threads;
|
|
std::vector<FeastBackendCall> calls;
|
|
if (mass != nullptr)
|
|
calls.push_back({0U, "mass-spd", "original-mass-factorization", 0.0,
|
|
mass_lower_bound, 0U, 0, false, 0U, false, 0, false, 0.0,
|
|
"accepted",
|
|
"sparse-SPD; caller-certified-roundoff-bound"});
|
|
auto TraceFailure = [&](const std::string& code, const std::string& identity,
|
|
const std::string& message) {
|
|
return Failure(code, identity,
|
|
mass != nullptr
|
|
? message + " Backend calls=" + CallHistory(calls)
|
|
: message);
|
|
};
|
|
std::array<MKL_INT, 128> snapshot{};
|
|
SetControls(snapshot, false);
|
|
std::vector<std::int64_t> base_fpm;
|
|
base_fpm.reserve(snapshot.size());
|
|
for (const MKL_INT value : snapshot)
|
|
base_fpm.push_back(static_cast<std::int64_t>(value));
|
|
std::vector<double> count_e;
|
|
std::vector<double> count_res;
|
|
std::vector<double> count_x;
|
|
try {
|
|
internal::ApplyFeastWorkspaceTestHook(true);
|
|
count_e.assign(n, 0.0);
|
|
count_res.assign(n, 0.0);
|
|
count_x.assign(n * n, 0.0);
|
|
} catch (const std::bad_alloc&) {
|
|
return Result<FeastEigenResult>::Failure(TraceFailure(
|
|
"eigen-allocation-failed", "count-workspace",
|
|
"Allocation of the required FEAST count workspace failed."));
|
|
}
|
|
|
|
auto Count = [&](const double call_lower, const double call_upper,
|
|
const std::string& purpose, std::size_t* count) -> Status {
|
|
if (!internal::CanCallFeast(calls.size()))
|
|
return TraceFailure("eigen-search-budget-exhausted", "count-calls",
|
|
"The bounded FEAST count-call budget was exhausted.");
|
|
std::array<MKL_INT, 128> fpm{};
|
|
SetControls(fpm, true);
|
|
MKL_INT m0 = static_cast<MKL_INT>(n);
|
|
MKL_INT returned = 0;
|
|
MKL_INT loop = 0;
|
|
MKL_INT info = 0;
|
|
double epsout = 0.0;
|
|
if (!FillInitialBasis(count_x, n, n)) {
|
|
return TraceFailure(
|
|
"eigenpair-validation-failed", "count-initial-basis",
|
|
"The prescribed count initial basis failed QR validation.");
|
|
}
|
|
const char uplo = 'U';
|
|
dfeast_scsrgv(&uplo, &mkl_n, a.Value().values.data(),
|
|
a.Value().offsets.data(), a.Value().columns.data(),
|
|
b.Value().values.data(), b.Value().offsets.data(),
|
|
b.Value().columns.data(), fpm.data(), &epsout, &loop,
|
|
&call_lower, &call_upper, &m0, count_e.data(), count_x.data(),
|
|
&returned, count_res.data(), &info);
|
|
internal::ApplyFeastBackendTestHook(
|
|
{true, calls.size(), call_lower, call_upper, n}, &m0, &returned,
|
|
&loop, &epsout, &info, count_e.data(), count_e.size(), count_x.data(),
|
|
count_x.size(), count_res.data(), count_res.size());
|
|
const bool count_available =
|
|
info == 0 && m0 >= 0 && static_cast<std::size_t>(m0) <= n;
|
|
FeastBackendCall event{calls.size(),
|
|
"count",
|
|
purpose,
|
|
call_lower,
|
|
call_upper,
|
|
n,
|
|
static_cast<int>(info),
|
|
count_available,
|
|
0U,
|
|
false,
|
|
0,
|
|
false,
|
|
0.0,
|
|
"",
|
|
""};
|
|
if (count_available) {
|
|
event.count = static_cast<std::size_t>(m0);
|
|
}
|
|
event.outcome = count_available ? "accepted" : "failed";
|
|
event.decision_reason =
|
|
count_available ? "accepted" : "vendor-info=" + std::to_string(info);
|
|
calls.push_back(std::move(event));
|
|
if (info != 0)
|
|
return TraceFailure(
|
|
"eigen-feast-failed", "info=" + std::to_string(info),
|
|
"FEAST did not certify the requested count interval.");
|
|
if (m0 < 0 || static_cast<std::size_t>(m0) > n) {
|
|
return TraceFailure("eigen-count-inconsistent", "count-output",
|
|
"FEAST count output is outside [0,n].");
|
|
}
|
|
*count = static_cast<std::size_t>(m0);
|
|
for (std::size_t i = 0U; i + 1U < calls.size(); ++i) {
|
|
const auto& prior = calls[i];
|
|
if (prior.kind != "count" || !prior.count_available ||
|
|
prior.lower != call_lower)
|
|
continue;
|
|
if ((prior.upper == call_upper && prior.count != *count) ||
|
|
(prior.upper < call_upper && prior.count > *count) ||
|
|
(prior.upper > call_upper && prior.count < *count))
|
|
return TraceFailure(
|
|
"eigen-count-inconsistent", "count-history",
|
|
"Successful prefix counts violate identity or monotonicity.");
|
|
}
|
|
return Status::Ok();
|
|
};
|
|
|
|
std::size_t count = 0U;
|
|
Status count_status = Status::Ok();
|
|
bool enclosure_certified = false;
|
|
for (std::size_t attempt = 0U; attempt < 3U; ++attempt) {
|
|
if (attempt > 0U) {
|
|
if (zero_operator) {
|
|
const double extent = std::ldexp(1.0, static_cast<int>(attempt));
|
|
lower = -extent;
|
|
upper = extent;
|
|
} else {
|
|
const double multiplier = std::ldexp(1.0, static_cast<int>(attempt));
|
|
lower =
|
|
std::nextafter(g_lower - 16.0 * multiplier * Delta(g_lower, beta),
|
|
-(std::numeric_limits<double>::infinity)());
|
|
upper = NextUp(g_upper + 16.0 * multiplier * Delta(g_upper, beta));
|
|
}
|
|
}
|
|
count_status = Count(lower, upper, "total-enclosure", &count);
|
|
if (count_status.IsOk() && count == n) {
|
|
enclosure_certified = true;
|
|
break;
|
|
}
|
|
const int info = calls.back().info;
|
|
const bool may_retry =
|
|
info == 1 || info == -5 || (count_status.IsOk() && count != n);
|
|
if (may_retry && attempt + 1U < 3U) {
|
|
calls.back().outcome = "retry";
|
|
calls.back().decision_reason =
|
|
count_status.IsOk() ? "count-mismatch"
|
|
: "retryable-vendor-info=" + std::to_string(info);
|
|
}
|
|
if (!may_retry) break;
|
|
}
|
|
if (!enclosure_certified) {
|
|
return Result<FeastEigenResult>::Failure(
|
|
!count_status.IsOk()
|
|
? count_status
|
|
: TraceFailure(
|
|
"eigen-count-inconsistent", "total-enclosure",
|
|
"The enclosing FEAST count did not equal matrix dimension."));
|
|
}
|
|
if (mass != nullptr) {
|
|
// Count the actual pencil before removing any part of the conservative
|
|
// negative enclosure. A far-negative endpoint degrades contour filtering
|
|
// of low positive modes on a stiff mesh; an empty guarded prefix proves
|
|
// that moving it cannot remove a mode, including a repeated rigid cluster.
|
|
const double margin = Delta(0.0, beta);
|
|
const double candidate = -16.0 * margin;
|
|
const std::array<double, 3U> guards{candidate - margin, candidate,
|
|
candidate + margin};
|
|
if (lower < guards.front()) {
|
|
for (const double guard : guards) {
|
|
std::size_t negative_count = 0U;
|
|
const Status prefix_status =
|
|
Count(lower, guard, "lower-empty-prefix", &negative_count);
|
|
if (!prefix_status.IsOk())
|
|
return Result<FeastEigenResult>::Failure(prefix_status);
|
|
if (negative_count != 0U)
|
|
return Result<FeastEigenResult>::Failure(TraceFailure(
|
|
"negative-eigenvalue-spectrum", "certified-negative-prefix",
|
|
"The generalized pencil contains significant negative "
|
|
"eigenvalues."));
|
|
}
|
|
lower = candidate;
|
|
calls.back().decision_reason = "guarded-empty-prefix; lowest-spectrum-preserved";
|
|
}
|
|
}
|
|
const double fixed_lower = lower;
|
|
const double enclosure_upper = upper;
|
|
double solve_upper = enclosure_upper;
|
|
std::string upper_search_termination =
|
|
n <= (std::min)(n, 2U * requested_modes + 8U) ? "count-target"
|
|
: "optimization-budget";
|
|
double search_lower = fixed_lower;
|
|
std::size_t certified_count = n;
|
|
std::vector<double> certified_uppers{enclosure_upper};
|
|
for (std::size_t midpoint_attempt = 0U;
|
|
upper_search_termination != "count-target" && midpoint_attempt < 32U;
|
|
++midpoint_attempt) {
|
|
if (!internal::CanStartUpperCandidate(calls.size())) {
|
|
upper_search_termination = "optimization-budget";
|
|
break;
|
|
}
|
|
const double midpoint = 0.5 * (search_lower + solve_upper);
|
|
if (!std::isfinite(midpoint) || midpoint == search_lower ||
|
|
midpoint == solve_upper) {
|
|
upper_search_termination = "no-midpoint";
|
|
break;
|
|
}
|
|
const double midpoint_delta = Delta(midpoint, beta);
|
|
const std::array<double, 7U> collision_offsets{0.0, 1.0, -1.0, 2.0,
|
|
-2.0, 4.0, -4.0};
|
|
bool certified = false;
|
|
bool optional_abandoned = false;
|
|
std::size_t candidate_count = 0U;
|
|
for (const double offset : collision_offsets) {
|
|
const double candidate = midpoint + offset * midpoint_delta;
|
|
const double margin = Delta(candidate, beta);
|
|
const std::array<double, 3U> guards{candidate - margin, candidate,
|
|
candidate + margin};
|
|
if (!(fixed_lower < guards[0U] && guards[2U] < enclosure_upper)) continue;
|
|
std::array<std::size_t, 3U> guard_counts{};
|
|
bool successful = true;
|
|
for (std::size_t guard = 0U; guard < guards.size(); ++guard) {
|
|
const Status guard_status = Count(fixed_lower, guards[guard],
|
|
"upper-guard", &guard_counts[guard]);
|
|
if (!guard_status.IsOk()) {
|
|
calls.back().outcome = "optional-abandoned";
|
|
const int optional_info = calls.back().info;
|
|
calls.back().decision_reason =
|
|
"optional-status=" + std::to_string(optional_info);
|
|
if (optional_info != 1 && optional_info != 2 && optional_info != 3 &&
|
|
optional_info != -5) {
|
|
return Result<FeastEigenResult>::Failure(guard_status);
|
|
}
|
|
optional_abandoned = true;
|
|
successful = false;
|
|
break;
|
|
}
|
|
}
|
|
if (optional_abandoned) break;
|
|
if (!successful || guard_counts[0U] != guard_counts[1U] ||
|
|
guard_counts[1U] != guard_counts[2U])
|
|
continue;
|
|
candidate_count = guard_counts[1U];
|
|
certified = true;
|
|
if (candidate_count >= requested_modes) {
|
|
solve_upper = candidate;
|
|
certified_uppers.push_back(candidate);
|
|
}
|
|
break;
|
|
}
|
|
if (optional_abandoned) {
|
|
upper_search_termination = "optional-status";
|
|
break;
|
|
}
|
|
if (!certified) {
|
|
upper_search_termination = "collision-exhausted";
|
|
break;
|
|
}
|
|
if (candidate_count >= requested_modes) {
|
|
certified_count = candidate_count;
|
|
if (candidate_count <= (std::min)(n, 2U * requested_modes + 8U)) {
|
|
upper_search_termination = "count-target";
|
|
break;
|
|
}
|
|
} else {
|
|
search_lower = midpoint;
|
|
}
|
|
if (solve_upper - search_lower <= 8.0 * Delta(midpoint, beta)) {
|
|
upper_search_termination = "width-limit";
|
|
break;
|
|
}
|
|
}
|
|
count = 0U;
|
|
count_status = Count(fixed_lower, solve_upper, "final-count", &count);
|
|
if (!count_status.IsOk() || count < requested_modes ||
|
|
count != certified_count) {
|
|
return Result<FeastEigenResult>::Failure(TraceFailure(
|
|
"eigen-count-inconsistent", "final-count",
|
|
"The final FEAST interval count is not the certified guarded count."));
|
|
}
|
|
if (solve_upper != enclosure_upper) {
|
|
const double margin = Delta(solve_upper, beta);
|
|
std::size_t below = 0U;
|
|
std::size_t above = 0U;
|
|
const Status below_status =
|
|
Count(fixed_lower, solve_upper - margin, "final-count", &below);
|
|
const Status above_status =
|
|
Count(fixed_lower, solve_upper + margin, "final-count", &above);
|
|
if (!below_status.IsOk() || !above_status.IsOk() || below != count ||
|
|
above != count) {
|
|
return Result<FeastEigenResult>::Failure(
|
|
TraceFailure("eigen-count-inconsistent", "final-guard",
|
|
"The final interior endpoint is not guarded by equal "
|
|
"successful counts."));
|
|
}
|
|
}
|
|
if (count > (std::numeric_limits<std::size_t>::max)() - 8U ||
|
|
count > ((std::numeric_limits<std::size_t>::max)() - 1U) / 3U) {
|
|
return Result<FeastEigenResult>::Failure(
|
|
TraceFailure("eigen-allocation-failed", "solve-capacity",
|
|
"FEAST solve-capacity padding overflows size_t."));
|
|
}
|
|
const std::size_t padded_count =
|
|
(std::max)(count + 8U, (3U * count + 1U) / 2U);
|
|
std::size_t capacity = count == n ? n : (std::min)(n, padded_count);
|
|
if (!internal::ProductFits(n, capacity))
|
|
return Result<FeastEigenResult>::Failure(
|
|
TraceFailure("feast-workspace-overflow", "solve-workspace",
|
|
"The FEAST solve workspace size overflows."));
|
|
if (!internal::FitsDoubleStorage(n * capacity) ||
|
|
!internal::FitsDoubleStorage(capacity)) {
|
|
return Result<FeastEigenResult>::Failure(TraceFailure(
|
|
"eigen-allocation-failed", "solve-workspace",
|
|
"The required solve workspace exceeds addressable double storage."));
|
|
}
|
|
std::vector<double> eigenvalues;
|
|
std::vector<double> residuals;
|
|
std::vector<double> vectors;
|
|
MKL_INT returned = 0;
|
|
MKL_INT info = 0;
|
|
bool solve_accepted = false;
|
|
std::size_t solve_capacity = capacity;
|
|
for (std::size_t attempt = 0U; attempt < 3U; ++attempt) {
|
|
if (!internal::CanCallFeast(calls.size()))
|
|
return Result<FeastEigenResult>::Failure(TraceFailure(
|
|
"eigen-search-budget-exhausted", "solve-call",
|
|
"The total backend call budget is exhausted before solve."));
|
|
try {
|
|
internal::ApplyFeastWorkspaceTestHook(false);
|
|
eigenvalues.assign(solve_capacity, 0.0);
|
|
residuals.assign(solve_capacity, 0.0);
|
|
vectors.assign(n * solve_capacity, 0.0);
|
|
} catch (const std::bad_alloc&) {
|
|
return Result<FeastEigenResult>::Failure(TraceFailure(
|
|
"eigen-allocation-failed", "solve-workspace",
|
|
"Allocation of the required FEAST solve workspace failed."));
|
|
}
|
|
if (!FillInitialBasis(vectors, n, solve_capacity)) {
|
|
return Result<FeastEigenResult>::Failure(TraceFailure(
|
|
"eigenpair-validation-failed", "solve-initial-basis",
|
|
"The prescribed solve initial basis failed QR validation."));
|
|
}
|
|
std::array<MKL_INT, 128> fpm{};
|
|
SetControls(fpm, false);
|
|
MKL_INT m0 = static_cast<MKL_INT>(solve_capacity);
|
|
returned = 0;
|
|
MKL_INT loop = 0;
|
|
info = 0;
|
|
double epsout = 0.0;
|
|
const char uplo = 'U';
|
|
dfeast_scsrgv(&uplo, &mkl_n, a.Value().values.data(),
|
|
a.Value().offsets.data(), a.Value().columns.data(),
|
|
b.Value().values.data(), b.Value().offsets.data(),
|
|
b.Value().columns.data(), fpm.data(), &epsout, &loop,
|
|
&fixed_lower, &solve_upper, &m0, eigenvalues.data(),
|
|
vectors.data(), &returned, residuals.data(), &info);
|
|
internal::ApplyFeastBackendTestHook(
|
|
{false, calls.size(), fixed_lower, solve_upper, solve_capacity}, &m0,
|
|
&returned, &loop, &epsout, &info, eigenvalues.data(),
|
|
eigenvalues.size(), vectors.data(), vectors.size(), residuals.data(),
|
|
residuals.size());
|
|
const bool solve_outputs_available =
|
|
info == 0 && loop >= 0 && std::isfinite(epsout) && returned >= 0 &&
|
|
static_cast<std::size_t>(returned) <= solve_capacity;
|
|
calls.push_back(
|
|
{calls.size(), "solve", "solve", fixed_lower, solve_upper,
|
|
solve_capacity, static_cast<int>(info), solve_outputs_available,
|
|
solve_outputs_available ? static_cast<std::size_t>(returned) : 0U,
|
|
solve_outputs_available,
|
|
solve_outputs_available ? static_cast<int>(loop) : 0,
|
|
solve_outputs_available, solve_outputs_available ? epsout : 0.0,
|
|
solve_outputs_available ? "accepted" : "failed",
|
|
solve_outputs_available ? "accepted"
|
|
: "vendor-info=" + std::to_string(info)});
|
|
if (info == 0 && (!std::isfinite(epsout) || loop < 0)) {
|
|
return Result<FeastEigenResult>::Failure(
|
|
TraceFailure("invalid-eigenvalue-result", "solve-output",
|
|
"A successful FEAST solve returned nonfinite epsout or "
|
|
"invalid loop."));
|
|
}
|
|
if (info == 0 && returned >= 0 &&
|
|
static_cast<std::size_t>(returned) <= solve_capacity) {
|
|
for (std::size_t column = 0U; column < static_cast<std::size_t>(returned);
|
|
++column) {
|
|
if (!std::isfinite(eigenvalues[column]))
|
|
return Result<FeastEigenResult>::Failure(TraceFailure(
|
|
"invalid-eigenvalue-result", std::to_string(column),
|
|
"A nonfinite eigenvalue cannot trigger an interval retry."));
|
|
for (std::size_t row = 0U; row < n; ++row)
|
|
if (!std::isfinite(vectors[row + column * n]))
|
|
return Result<FeastEigenResult>::Failure(TraceFailure(
|
|
"invalid-eigenvalue-result", std::to_string(column),
|
|
"A nonfinite eigenvector cannot trigger an interval retry."));
|
|
}
|
|
}
|
|
bool endpoint_separated = info == 0 && returned >= 0 &&
|
|
static_cast<std::size_t>(returned) == count;
|
|
if (endpoint_separated) {
|
|
for (std::size_t index = 0U; index < count; ++index) {
|
|
if (!std::isfinite(eigenvalues[index])) {
|
|
return Result<FeastEigenResult>::Failure(TraceFailure(
|
|
"invalid-eigenvalue-result", std::to_string(index),
|
|
"FEAST returned a nonfinite eigenvalue before endpoint "
|
|
"validation."));
|
|
}
|
|
endpoint_separated =
|
|
endpoint_separated &&
|
|
eigenvalues[index] - fixed_lower >= Delta(fixed_lower, beta) &&
|
|
solve_upper - eigenvalues[index] >= Delta(solve_upper, beta);
|
|
}
|
|
}
|
|
if (endpoint_separated) {
|
|
solve_accepted = true;
|
|
break;
|
|
}
|
|
const bool retryable =
|
|
!endpoint_separated &&
|
|
(info == 1 || info == 2 || info == 3 ||
|
|
(info == 0 &&
|
|
(returned < 0 || static_cast<std::size_t>(returned) != count ||
|
|
!endpoint_separated)));
|
|
if (!retryable || attempt == 2U) break;
|
|
calls.back().outcome = "retry";
|
|
calls.back().decision_reason =
|
|
info == 0 ? "endpoint-or-count-mismatch"
|
|
: "retryable-vendor-info=" + std::to_string(info);
|
|
double larger_upper = enclosure_upper;
|
|
for (const double certified : certified_uppers) {
|
|
if (certified > solve_upper && certified < larger_upper)
|
|
larger_upper = certified;
|
|
}
|
|
const bool interval_changed = larger_upper > solve_upper;
|
|
if (interval_changed) {
|
|
solve_upper = larger_upper;
|
|
std::size_t recounted = 0U;
|
|
const Status recount_status =
|
|
Count(fixed_lower, solve_upper, "final-count", &recounted);
|
|
if (!recount_status.IsOk() || recounted < requested_modes) {
|
|
return Result<FeastEigenResult>::Failure(TraceFailure(
|
|
"eigen-count-inconsistent", "solve-retry",
|
|
"A larger certified solve interval could not be recounted."));
|
|
}
|
|
count = recounted;
|
|
if (solve_upper != enclosure_upper) {
|
|
const double retry_margin = Delta(solve_upper, beta);
|
|
std::size_t below = 0U;
|
|
std::size_t above = 0U;
|
|
const Status below_status = Count(
|
|
fixed_lower, solve_upper - retry_margin, "final-count", &below);
|
|
const Status above_status = Count(
|
|
fixed_lower, solve_upper + retry_margin, "final-count", &above);
|
|
if (!below_status.IsOk() || !above_status.IsOk() || below != count ||
|
|
above != count) {
|
|
return Result<FeastEigenResult>::Failure(TraceFailure(
|
|
"eigen-count-inconsistent", "solve-retry-guard",
|
|
"The larger retry interval is not guarded by matching counts."));
|
|
}
|
|
}
|
|
if (count > (std::numeric_limits<std::size_t>::max)() - 8U ||
|
|
count > ((std::numeric_limits<std::size_t>::max)() - 1U) / 3U) {
|
|
return Result<FeastEigenResult>::Failure(
|
|
TraceFailure("eigen-allocation-failed", "solve-capacity",
|
|
"FEAST retry capacity padding overflows size_t."));
|
|
}
|
|
const std::size_t retry_padding =
|
|
(std::max)(count + 8U, (3U * count + 1U) / 2U);
|
|
capacity = count == n ? n : (std::min)(n, retry_padding);
|
|
}
|
|
const std::size_t previous_capacity = solve_capacity;
|
|
const std::size_t doubled =
|
|
solve_capacity > n / 2U ? n : 2U * solve_capacity;
|
|
solve_capacity = (std::max)(capacity, doubled);
|
|
if (solve_capacity == previous_capacity && !interval_changed) break;
|
|
}
|
|
if (!solve_accepted) {
|
|
return Result<FeastEigenResult>::Failure(
|
|
TraceFailure("eigen-feast-failed", "info=" + std::to_string(info),
|
|
"FEAST solve did not return the certified interval count "
|
|
"within its bounded retries."));
|
|
}
|
|
|
|
for (std::size_t index = 0U; index < count; ++index) {
|
|
if (!std::isfinite(eigenvalues[index])) {
|
|
return Result<FeastEigenResult>::Failure(TraceFailure(
|
|
"invalid-eigenvalue-result", std::to_string(index),
|
|
"FEAST returned a nonfinite eigenvalue before sorting."));
|
|
}
|
|
for (std::size_t row = 0U; row < n; ++row) {
|
|
if (!std::isfinite(vectors[row + index * n])) {
|
|
return Result<FeastEigenResult>::Failure(TraceFailure(
|
|
"invalid-eigenvalue-result", std::to_string(index),
|
|
"FEAST returned a nonfinite eigenvector before sorting."));
|
|
}
|
|
}
|
|
}
|
|
|
|
std::vector<std::size_t> order(count);
|
|
std::iota(order.begin(), order.end(), 0U);
|
|
std::stable_sort(order.begin(), order.end(),
|
|
[&](const std::size_t left, const std::size_t right) {
|
|
return eigenvalues[left] < eigenvalues[right];
|
|
});
|
|
Matrix ordered_vectors(n, count);
|
|
std::vector<double> ordered_values;
|
|
std::vector<bool> ordered_backend_available;
|
|
std::vector<double> ordered_backend_res;
|
|
std::vector<double> ordered_mu;
|
|
std::vector<double> ordered_rho_b;
|
|
ordered_values.reserve(count);
|
|
ordered_backend_available.reserve(count);
|
|
ordered_backend_res.reserve(count);
|
|
ordered_mu.reserve(count);
|
|
ordered_rho_b.reserve(count);
|
|
for (std::size_t output = 0U; output < count; ++output) {
|
|
const std::size_t input = order[output];
|
|
const double mu = eigenvalues[input];
|
|
if (!std::isfinite(mu)) {
|
|
return Result<FeastEigenResult>::Failure(
|
|
TraceFailure("invalid-eigenvalue-result", std::to_string(input),
|
|
"FEAST returned a nonfinite eigenvalue."));
|
|
}
|
|
if (mu < -1.0e-9) {
|
|
return Result<FeastEigenResult>::Failure(TraceFailure(
|
|
"negative-eigenvalue-spectrum", std::to_string(input),
|
|
"The raw normalized FEAST eigenvalue is significantly negative."));
|
|
}
|
|
if (mu - fixed_lower < Delta(fixed_lower, beta) ||
|
|
solve_upper - mu < Delta(solve_upper, beta)) {
|
|
return Result<FeastEigenResult>::Failure(TraceFailure(
|
|
"eigen-endpoint-separation", std::to_string(input),
|
|
"A returned FEAST eigenvalue is too close to a solve endpoint."));
|
|
}
|
|
Vector raw(n);
|
|
for (std::size_t row = 0U; row < n; ++row)
|
|
raw[row] = vectors[row + input * n];
|
|
double norm = raw.Dot(raw);
|
|
if (mass != nullptr) {
|
|
norm = raw.Dot(mass->Multiply(raw));
|
|
}
|
|
norm = std::sqrt(norm);
|
|
if (!(norm > 0.0) || !std::isfinite(norm))
|
|
return Result<FeastEigenResult>::Failure(
|
|
TraceFailure("eigenpair-validation-failed", std::to_string(input),
|
|
"FEAST returned a non-normalizable eigenvector."));
|
|
double residual_norm = 0.0;
|
|
const auto action = action_matrix.Value().Multiply(raw);
|
|
const auto mass_action = mass != nullptr ? mass->Multiply(raw) : raw;
|
|
for (std::size_t row = 0U; row < n; ++row) {
|
|
const double difference = action[row] - mu * mass_action[row];
|
|
residual_norm += difference * difference;
|
|
ordered_vectors(row, output) = vectors[row + input * n] / norm;
|
|
if (!std::isfinite(ordered_vectors(row, output)))
|
|
return Result<FeastEigenResult>::Failure(
|
|
TraceFailure("eigenpair-validation-failed", std::to_string(input),
|
|
"FEAST returned a non-normalizable eigenvector."));
|
|
}
|
|
double residual_denominator = norm;
|
|
if (mass != nullptr) {
|
|
double norm_k = 0.0;
|
|
double norm_m = 0.0;
|
|
double euclidean = 0.0;
|
|
for (const double value : values)
|
|
norm_k = std::hypot(norm_k, value / operator_scale);
|
|
for (const double value : mass->Values())
|
|
norm_m = std::hypot(norm_m, value);
|
|
for (std::size_t row = 0U; row < n; ++row)
|
|
euclidean = std::hypot(euclidean, vectors[row + input * n]);
|
|
residual_denominator = (norm_k + std::abs(mu) * norm_m) * euclidean;
|
|
}
|
|
if (!std::isfinite(residual_norm) || !std::isfinite(residual_denominator))
|
|
return Result<FeastEigenResult>::Failure(
|
|
TraceFailure("eigenpair-validation-failed", std::to_string(input),
|
|
"A nonfinite independent residual numerator or scale "
|
|
"cannot certify a pair."));
|
|
const double rho_b =
|
|
residual_denominator == 0.0
|
|
? (residual_norm == 0.0 ? 0.0
|
|
: (std::numeric_limits<double>::infinity)())
|
|
: std::sqrt(residual_norm) / residual_denominator;
|
|
if (!std::isfinite(rho_b) || rho_b > 1.0e-10) {
|
|
if (mass != nullptr)
|
|
return Result<FeastEigenResult>::Failure(TraceFailure(
|
|
"eigenpair-validation-failed", std::to_string(input),
|
|
"Independent original-pencil residual exceeds 1e-10: eta=" +
|
|
FormatDouble(rho_b) + ", mu=" + FormatDouble(mu) +
|
|
"; calls=" + CallHistory(calls)));
|
|
const FeastBackendCall& last_call = calls.back();
|
|
auto mirrored_matrix =
|
|
ScaledActionMatrix(whitened_operator, operator_scale, true);
|
|
if (!mirrored_matrix.HasValue())
|
|
return Result<FeastEigenResult>::Failure(mirrored_matrix.GetStatus());
|
|
const auto mirrored_action = mirrored_matrix.Value().Multiply(raw);
|
|
double mirrored_residual_squared = 0.0;
|
|
double action_difference_squared = 0.0;
|
|
for (std::size_t row = 0U; row < n; ++row) {
|
|
const double mirrored_residual =
|
|
mirrored_action[row] - mu * vectors[row + input * n];
|
|
mirrored_residual_squared += mirrored_residual * mirrored_residual;
|
|
const double action_difference =
|
|
action[row] - mirrored_action[row];
|
|
action_difference_squared += action_difference * action_difference;
|
|
}
|
|
const double rho_b_upper = std::sqrt(mirrored_residual_squared) / norm;
|
|
const double action_difference_norm =
|
|
std::sqrt(action_difference_squared);
|
|
return Result<FeastEigenResult>::Failure(TraceFailure(
|
|
"eigenpair-validation-failed", std::to_string(input),
|
|
"Independent whitened residual exceeds 1e-10: rho_b=" +
|
|
FormatDouble(rho_b) + ", mu=" + FormatDouble(mu) +
|
|
", norm=" + FormatDouble(norm) +
|
|
", backend_res=" + FormatDouble(residuals[input]) +
|
|
", info=" + std::to_string(last_call.info) +
|
|
", epsout=" + FormatDouble(last_call.epsout) +
|
|
", loop=" + std::to_string(last_call.loop) +
|
|
", count=" + std::to_string(count) + ", m0=" +
|
|
std::to_string(last_call.allocated_m0) + ", interval=(" +
|
|
FormatDouble(fixed_lower) + "," + FormatDouble(solve_upper) +
|
|
"), calls=" + std::to_string(calls.size()) +
|
|
", history=" + CallHistory(calls) +
|
|
", rho_b_upper=" + FormatDouble(rho_b_upper) +
|
|
", full_upper_action_difference_norm=" +
|
|
FormatDouble(action_difference_norm) + "."));
|
|
}
|
|
const double lambda = mu * operator_scale;
|
|
if (!std::isfinite(lambda) || (mu != 0.0 && lambda == 0.0)) {
|
|
return Result<FeastEigenResult>::Failure(TraceFailure(
|
|
"invalid-eigenvalue-result", std::to_string(input),
|
|
"Physical eigenvalue conversion overflowed or underflowed."));
|
|
}
|
|
ordered_values.push_back(lambda);
|
|
const bool backend_res_finite = std::isfinite(residuals[input]);
|
|
if (!backend_res_finite && mu != 0.0) {
|
|
return Result<FeastEigenResult>::Failure(TraceFailure(
|
|
"invalid-eigenvalue-result", std::to_string(input),
|
|
"A nonzero FEAST mode has a nonfinite backend residual."));
|
|
}
|
|
ordered_backend_available.push_back(backend_res_finite);
|
|
ordered_backend_res.push_back(backend_res_finite ? residuals[input] : 0.0);
|
|
ordered_mu.push_back(mu);
|
|
ordered_rho_b.push_back(rho_b);
|
|
}
|
|
double gram_error = 0.0;
|
|
double gram_squared = 0.0;
|
|
Matrix mass_modes{n, count};
|
|
for (std::size_t right = 0; right < count; ++right) {
|
|
Vector mode(n);
|
|
for (std::size_t row = 0; row < n; ++row)
|
|
mode[row] = ordered_vectors(row, right);
|
|
const auto action = mass != nullptr ? mass->Multiply(mode) : std::move(mode);
|
|
for (std::size_t row = 0; row < n; ++row) {
|
|
if (!std::isfinite(action[row]))
|
|
return Result<FeastEigenResult>::Failure(
|
|
TraceFailure("eigenpair-validation-failed", "gram",
|
|
"The independent mass action is nonfinite."));
|
|
mass_modes(row, right) = action[row];
|
|
}
|
|
}
|
|
const auto gram = dense_blas_internal::MultiplyTransposeLeft(ordered_vectors,
|
|
mass_modes);
|
|
for (std::size_t left = 0U; left < count; ++left) {
|
|
for (std::size_t right = 0U; right < count; ++right) {
|
|
const double dot = gram(left, right);
|
|
const double expected = left == right ? 1.0 : 0.0;
|
|
const double difference = std::abs(dot - expected);
|
|
gram_squared += difference * difference;
|
|
gram_error = (std::max)(gram_error, difference);
|
|
if (std::abs(dot - expected) > 1.0e-10)
|
|
return Result<FeastEigenResult>::Failure(TraceFailure(
|
|
"eigenpair-validation-failed",
|
|
std::to_string(left) + ":" + std::to_string(right),
|
|
"The all-c FEAST vectors fail the independent Gram check."));
|
|
}
|
|
}
|
|
gram_error = std::sqrt(gram_squared /
|
|
(mass != nullptr ? 1.0 : static_cast<double>(count)));
|
|
if (!std::isfinite(gram_error) || gram_error > 1.0e-10) {
|
|
return Result<FeastEigenResult>::Failure(TraceFailure(
|
|
"eigenpair-validation-failed", "gram",
|
|
"The normalized all-c Frobenius Gram error exceeds 1e-10: " +
|
|
FormatDouble(gram_error) + "."));
|
|
}
|
|
char version[256]{};
|
|
mkl_get_version_string(version, static_cast<int>(sizeof(version)));
|
|
FeastEigenResult result{requested_modes,
|
|
std::move(ordered_values),
|
|
std::move(ordered_vectors),
|
|
operator_scale,
|
|
order,
|
|
std::move(ordered_backend_available),
|
|
std::move(ordered_backend_res),
|
|
std::move(ordered_mu),
|
|
std::move(ordered_rho_b),
|
|
gram_error,
|
|
std::move(calls),
|
|
std::move(base_fpm),
|
|
version,
|
|
"ILP64",
|
|
1U,
|
|
kSeed,
|
|
fixed_lower,
|
|
enclosure_upper,
|
|
fixed_lower,
|
|
solve_upper,
|
|
upper_search_termination};
|
|
return Result<FeastEigenResult>::Success(std::move(result));
|
|
} catch (const std::invalid_argument&) {
|
|
return Result<FeastEigenResult>::Failure(
|
|
Failure("eigenpair-validation-failed", "independent-action",
|
|
"Modal product inputs must be finite and shape compatible."));
|
|
} catch (const std::length_error&) {
|
|
return Result<FeastEigenResult>::Failure(
|
|
Failure("feast-integer-overflow", "independent-action",
|
|
"Independent sparse action exceeds the backend range."));
|
|
} catch (const std::runtime_error&) {
|
|
return Result<FeastEigenResult>::Failure(
|
|
Failure("eigen-multiplication-failed", "independent-action",
|
|
"Independent sparse action failed."));
|
|
}
|
|
|
|
} // namespace fesa
|