perf(assembly): use MKL entry products with ordered load diagnostics
This commit is contained in:
@@ -11,6 +11,7 @@
|
||||
#include "fesa/constraints/essential_constraint_policy.h"
|
||||
#include "fesa/loads/load.h"
|
||||
#include "fesa/model/source_target_resolver.h"
|
||||
#include "math/dense_blas_internal.h"
|
||||
|
||||
namespace fesa {
|
||||
namespace {
|
||||
@@ -279,13 +280,35 @@ Result<Vector> LoadAssembler::EffectiveFreeRhs(const Vector& full_load,
|
||||
return Result<Vector>::Failure(prescribed_status);
|
||||
}
|
||||
|
||||
Vector entry_coefficients{kfc.Values().size()};
|
||||
Vector entry_prescribed{kfc.Values().size()};
|
||||
for (std::size_t position = 0U; position < kfc.Values().size(); ++position) {
|
||||
entry_coefficients[position] = kfc.Values()[position];
|
||||
entry_prescribed[position] =
|
||||
prescribed_values[kfc.ColumnIndices()[position]];
|
||||
}
|
||||
Vector entry_products{0U};
|
||||
try {
|
||||
entry_products = dense_blas_internal::MultiplyDiagonal(entry_coefficients,
|
||||
entry_prescribed);
|
||||
} catch (const std::length_error&) {
|
||||
return Result<Vector>::Failure(LoadFailure(
|
||||
"load-backend-dimension-overflow", location, "LOAD_ASSEMBLER", "Kfc",
|
||||
"Prescribed displacement products exceed the backend integer range."));
|
||||
} catch (const std::invalid_argument&) {
|
||||
return Result<Vector>::Failure(LoadFailure(
|
||||
"invalid-load-product", location, "LOAD_ASSEMBLER", "Kfc",
|
||||
"Prescribed displacement product inputs must be finite and compatible."));
|
||||
}
|
||||
|
||||
// MKL computes independent entries; ordered CSR reduction below retains
|
||||
// cancellation and the first product/row-sum overflow diagnostic.
|
||||
Vector correction{kfc.Rows()};
|
||||
for (std::size_t row = 0U; row < kfc.Rows(); ++row) {
|
||||
double sum = 0.0;
|
||||
for (std::size_t position = kfc.RowOffsets()[row];
|
||||
position < kfc.RowOffsets()[row + 1U]; ++position) {
|
||||
const double product = kfc.Values()[position] *
|
||||
prescribed_values[kfc.ColumnIndices()[position]];
|
||||
const double product = entry_products[position];
|
||||
if (!std::isfinite(product)) {
|
||||
return Result<Vector>::Failure(LoadFailure(
|
||||
"nonfinite-load-accumulation", location, "LOAD_ASSEMBLER",
|
||||
|
||||
Reference in New Issue
Block a user