diff --git a/docs/calculation-performance/t4-shell-optimization.md b/docs/calculation-performance/t4-shell-optimization.md new file mode 100644 index 0000000..4e0d940 --- /dev/null +++ b/docs/calculation-performance/t4-shell-optimization.md @@ -0,0 +1,162 @@ +# Shell 계산 경로의 성능 변경 + +이 기록은 reduced-coordinate shell의 local 계산 변경과 근거를 구분한다. +전체 시간 목표와 최종 원본 비교는 [PRD](PRD.md)의 계약이며 아래의 단계 검사는 +그 계약의 최종 합격을 대신하지 않는다. + +## 계측과 변경의 연결 + +원본 meshfine shell의 기준 전체 실행은 6372.3256209초다. VTune 소프트웨어 +샘플링의 별도 150초 부분 수집은 accepted increment 2 / time 0.2까지 포함한다. +`Jet*`는 sampled CPU 133.466초(27.6%), `malloc`은 30.289초(6.3%)였다. +`EvaluateAssembly`의 inclusive 257.862958 CPU초와 중첩 selector의 1.600014초는 +겹치는 호출 경로이므로 더하거나 전체 경과 시간으로 바꾸지 않는다. + +`nonlinear_mitc4_shell.cpp`의 모든 Jet는 기존에 값이 상수인 경우도 576-double +Hessian을 영 초기화한 `shared_ptr`로 할당했다. 상수의 영 Hessian은 이제 +읽기 전용 정적 영 배열을 참조하고 처음 쓰는 때만 기존 공유 storage를 할당한다. +읽기 전용 영 배열에는 trial 값이 없으며 동시 호출은 그 배열을 수정하지 않는다. +이미 할당된 Hessian의 복사/공유 수명은 그대로다. 현재 trial 이외의 형상·접선· +도함수 응답을 저장하거나 주소에 따라 재사용하는 cache는 추가하지 않았다. + +후속 후보는 576-entry 반복문 앞에서 입력 read pointer와 fresh output pointer를 +한 번씩 취득한다. fresh output은 단일 `make_shared` 할당을 유지하고 영 초기화를 +생략한 뒤 576개 전체를 기록한다. `operator+/-/*`, scalar scale, Compose, Atan2의 +모든 complete write와 MKL 거부 시의 NaN fill 뒤에만 결과를 읽거나 복사한다. +일반 부분 mutable write는 영 초기화를 유지한다. 상수의 nullptr 복사에는 별도 +mutable 공유 객체가 없지만 실제 consumer는 복사된 상수 Hessian을 수정하지 않는다. +이미 할당된 결과의 shared copy와 입력 alias는 수정하지 않는다. + +## 일반 행렬 연산의 MKL 경계 + +| 파일 / 함수 | 연산과 크기 | 변경 | +| --- | --- | --- | +| `mitc4_shell.cpp::Stiffness`, `Congruence` | 5×20 strain의 BᵀCB, 20/4 local→24 congruence | `Matrix::Multiply`와 private transpose-left MKL GEMM, 상삼각 결과를 그대로 대칭 복사 | +| `nonlinear_mitc4_inertia.cpp::CheckPhysicalMass` | 24×20 complement Gram·mass congruence, 20×20 Cholesky·inverse·reconstruction | GEMM 및 기존 private LAPACK helper; 같은 positivity, reconstruction, 2^-32 conditioning gate 유지 | +| `nonlinear_mitc4_inertia.cpp::Evaluate` | 각 lumped mass-node의 6×3 coefficient Gram, velocity/position bilinear, 3방향 mass variation의 6×6 block | 무할당 row-major transpose-right GEMM; 원래 point/node/row/column scatter와 가중 누적 유지 | +| `nonlinear_mitc4_shell.cpp::operator*`, `Compose`, `Atan2`, `DrillingEnergy` | ordinary 24×24 derivative outer products | private `OuterProduct`의 row-major GEMM k=1, 기존 Jet storage에 직접 기록 | +| `hht_dynamic_algebra.cpp::BuildPhysicalSpace` | scaled free mass×analytic null 및 일반 complement congruence의 sparse contributions | entrywise compact-diagonal MKL products 뒤 원래 CSR/COO 순서 누적; full global dense matrix 없음 | +| `hht_dynamic_algebra.cpp::CheckLocalMass` | nonempty nodal physical block, 크기1–6 | MKL Cholesky/inverse; 동일 inverse infinity norm과 2^-32 conditioning 검사 | +| `hht_dynamic_algebra.cpp::CheckRegularity`, `hht_reduced_operators.cpp::CheckGraphRegularity` | scalar null stiffness metric와 gradient×null projection | MKL GEMV/dot 또는 entry products, 기존 fresh factorization/모든 unit RHS/residual/condition gate 유지 | +| `hht_reduced_operators.cpp::ContractConstraintHessian` | failure-only bilinear replay의 sparse Hessian entries | 두 MKL diagonal product 뒤 기존 left association/entry order sum | + +Outer product는 배열 copy나 Matrix 객체를 만들지 않는다. 크기·integer range· +nonempty null storage를 검사하며 비유한 입력은 MKL 호출/출력 쓰기 전에 거부한다. +유한 입력의 곱 overflow는 후속 owner finite gate가 판단한다. Jet consumer는 +거부된 operand의 결과를 NaN으로 표시하고 원래 equilibrium/constraint finite +실패를 반환한다. 직접 scalar product로 계산하는 fallback은 없다. + +`CholeskyLowerAndInverse`의 optional private failure reason은 validation 전에 +`kNone`으로 초기화된다. factor pivot/finite 실패와 inverse/finite 실패를 구분하여 +기존 positive-mass와 ill-conditioned-mass 진단을 유지한다. 실패에서는 lower와 +inverse output candidate를 교체하지 않는다. 기존 default-nullptr consumer의 +bool 의미는 그대로이며 live inverse 실패를 발생시킨 검사는 아직 없다. + +Element-local dense work의 `ScopedLocalDenseThreads`는 caller TLS 요청이 default +0일 때만 1을 적용한다. caller의 양수 override는 유지하고 정상 반환과 unwind 뒤 +원래 default를 복원한다. 전역 MKL/OpenMP 정책이나 이후 global solve에는 영향을 +주지 않는다. 각 scalar metric scope는 MKL 호출 뒤 끝나며 global PARDISO를 포함하지 +않는다. 현재 후보의 실제 비용과 복원 검사는 아래 단계 검증 기록으로 구분한다. +TLS default 0, caller override 우선순위와 저장한 값을 복원하는 동작은 +[Intel 2026 thread-local reference](https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2026-0/mkl-set-num-threads-local.html)의 계약을 따른다. + +Hessian의 합산 순서와 ordinary product는 별도다. 곱셈 Jet는 원래의 +`left.second * right.value + outer1 + outer2 + left.value * right.second`를 +유지한다. Compose와 Atan2는 계수×gradient를 먼저 계산하여 원래 left association을 +유지한 뒤 outer product를 계산한다. GEMM의 연산과 signed-zero 결과는 원래 scalar +구현의 bit 동일성을 주장하지 않으며 접선/방향/혼합 도함수의 수치 검토 대상이다. +일반 24×24 연산을 AD라는 이유로 특수 정밀도 예외로 처리하지 않는다. + +## 유지하는 특수 연산 + +| 파일 / 함수 | 크기·정밀도·이유 | 검증 범위 | +| --- | --- | --- | +| `nonlinear_mitc4_inertia.cpp::Jet`, `ExpSkew`, `LeftJacobian` | 3개의 nodal rotation coordinate에 대한 최대 3차 analytic derivative; 3×3×3 tensor | finite rotation, inertia force/position/velocity tangent와 null·momentum oracle | +| `hht_reduced_operators.cpp::BasisValue`, fixed basis 생성 | normalized high+low 및 outward radius; 3-vector basis certificate | fixed chart/radius replay·partial-boundary·selector transversality | +| `hht_reduced_operators.cpp::RateSum`, `GradientDot`, reduced block 누적 | sparse support의 원래 순서 Neumaier correction | cancellation, full defect, mixed rates·channel acceptance | +| `nonlinear_mitc4_inertia.cpp::InertiaSum` | kinetic force/energy/momentum의 원래 순서 compensated contraction | 독립 inertia oracle와 omission 검출, 좌표/단위 변화 | +| shell의 fixed3 coordinate tensor/rotation 변환 | 좌표·director·chart의 3-vector 또는 3×3 연산 | rigid motion/objectivity·principal chart 실패 | +| `hht_dynamic_algebra.cpp::BuildPhysicalSpace`의 analytic nodal basis Gram, `InitializeAcceleration`의 force projection/lift | translation 1계수 또는 단일 절점의 최대3 rotation components; 현재 assembly의 nodal frame에 대한 Float64 좌표 변환 | 서로 다른 node의 basis support는 disjoint, 같은 node의 최대3-vector inner product만 사용; excluded fixed rows·원래 association·초기 acceleration 진단 순서 유지 | + +보상 누적과 interval/basis 인증은 일반 Float64 GEMM으로 접어 수치 의미를 잃지 않는다. +`hht_reduced_operators.cpp::GradientDot`, `Lift`, `ProjectPhysical`, +`ReconstructPhysicalRates`, `AssembleAugmentedReducedSystem`의 support contractions는 +원래 entry 순서의 Neumaier correction과 high+low/outward fixed-basis certificate를 +유지하는 예외다. Plain sum과 달리 partition cancellation에서 소실되는 1을 유지하는 +`ReconstructsRatesAcrossSelectedPartitionCancellation` 및 augmented full-model defect +검사가 해당 의미를 검증한다. 일반 24×24 outer와 mass matrix 제품은 이 예외로 묶지 않는다. +`BuildPhysicalSpace`의 좌표 basis는 현재 assembly의 analytic null로 각 호출에서 +생성한다. 초기화는 initial assembly를 사용하고 후속 질량 인증은 해당 현재 assembly를 +사용한다. 실행 전체의 fixed initial chart와 같은 자료라고 판단하지 않으며 +다른 trial에 basis/응답을 재사용하지 않는다. + +## 단계 검증 상태 + +변경 전 clean source `8c7c41729eb9eaa904f6739b54d34c37b1bd555f`의 fresh MSVC +Debug configure/build는 native 0, local shell/inertia/reduced/math 70개는 모두 PASS다. +별도 MITC4 linear kernel·mass·recovery 20개도 native 0으로 모두 PASS다. +Mixed-time actual shell fixture는 144.84초에 정상 통과했으며 중단되지 않았다. +[원래 실행 기록](../../.harness/performance/t4/baseline-native.json)과 +[linear 실행 기록](../../.harness/performance/t4/baseline-linear-native.json)을 보존한다. + +첫 후보 build의 `Status` 기본 생성자 오류(native 2)는 별도로 보존하고 +`Status::Ok` 초기화로 수정했다. 수정 후보의 configure/build는 native 0이고 +91개 중 90 PASS, `HhtReducedEmbedding.ReconstructsActualShellRatesWithMixedTimeLoad`는 +180초 timeout으로 미통과다. 전체 CTest native 8과 실제 종료 상태를 보존한다. +이 Debug 비용 악화를 Release 개선으로 표시하지 않는다. 이는 +`candidate-local-fixed` snapshot의 결과다. Read pointer hoist, fresh output의 +단일 allocation/full overwrite, 별도 이름의 local TLS scope와 추가 coverage를 +적용한 `packed-output-mkl-coverage` 결과는 별도로 기록한다. + +기존 Korean compiler 환경의 Ninja header prefix가 mojibake이고 actual dependency +graph가 0개인 사실을 확인했다. 그 build의 직접 변경 source/검사 결과는 보존하되 +header consumer 전체 재빌드를 입증한다고 해석하지 않는다. 설치된 compiler의 +locale fallback 때문에 `VSLANG=1033`만으로 해결되지 않는다는 root의 후속 근거가 +있다. 정확한 include prefix를 지정한 ignored metadata helper의 별도 probe 뒤 +full fresh build에서 actual Ninja header deps를 확인했다. 공용 preset/flag나 +source object를 바꾸어 혼합하지 않는다. + +후속 `packed-output-mkl-coverage` snapshot의 source SHA-256은 +[source manifest](../../.harness/performance/t4/candidate-packed-source.json)에 있다. +Corrected-prefix fresh configure/build는 각각 native 0이고 +`eigen_reduction.cpp.obj`의 실제 dependency inventory에 +`dense_blas_internal.h`가 포함되어 있다. +기존 Cholesky/eigen consumer, local tangent/objectivity/patch·inertia, +fixed chart·channel/failure·boundary·MKL scope/shape/finite 검사 **105/105 PASS**, +native 0, CTest 36.94초다. 별도 mixed-time fixture도 **1/1 PASS**, native 0, +fixture 60.69초다. [실행 기록](../../.harness/performance/t4/packed-native.json), +[검사 로그](../../.harness/performance/t4/packed-4.txt), +[mixed-time 기록](../../.harness/performance/t4/packed-mixed-native.json)을 보존한다. +이 106개 결과는 전체 Debug/원본 비교/시간수렴을 통과했다는 뜻이 아니다. + +첫 Release 보조 configure는 runtime DLL PATH를 compiler 탐색보다 앞에 두어 +icx가 선택되었고 `libircmt.lib` 없음으로 native 1 실패했다. 수치 실행은 없었으며 +실패 자료/cache를 보존한다. 동일 설치 MSVC cl을 명시한 별도 fresh directory의 +Release configure/build 97 edges는 native 0이다. [보조 실행 기록](../../.harness/performance/t4/release-micro-fixed-native.json)을 보존한다. + +Jet 보조 측정은 frozen `8c7c417`의 실제 Jet source span과 현재 production +`nonlinear_mitc4_shell.cpp`를 같은 translation unit에 두고 현재 source의 +Release library와 연결했다. 동일 입력·flags로 각 variant의 3000-product loop를 +워밍업 없이 한 번씩 측정했다. 모든 value/gradient/Hessian/directional/mixed +coefficient 및 checksum이 같고 실제 native 0, scope 안 effective threads 1, +종료 뒤 caller TLS 복원 true다. 정확성 대조는 모든 타이밍 후 수행했다. + +| 보조 variant | 3000-product loop 경과 시간 | +| --- | ---: | +| frozen original | 0.0025384초 | +| MKL candidate / default 8 | 0.0082295초 | +| MKL candidate / element scope 1 | 0.0029779초 | + +Scope는 새 MKL kernel의 default-thread 비용을 줄이지만 이 dense nonzero +kernel 자체가 원래 구현보다 빠르다는 근거는 아니다. Scope 후보도 original보다 +약17% 느리다. 전체 모델의 준비·조립·allocator·선형 해법·출력 비용과 최종 +100초 목표는 이 보조 값으로 판단하지 않는다. 실제 production의 Constant/Jet +allocation 빈도 감소는 최종 VTune에서 확인해야 한다. + +Standalone include의 vendor C4819 경고는 보존했다. 실제 production과 같은 +`/W4 /WX /external:W0`로 compile/link만 추가 실행하여 native 0/경고 없음이고 +수치 타이밍은 반복하지 않았다. 측정된 exe는 별도 원본 copy로 보존했으며 +추가 exe의 `.text` SHA-256이 동일한 것을 확인했다. +[기계어/실행 파일 보존 기록](../../.harness/performance/t4/micro-warning-policy-text-comparison.json)을 따른다. +전체 Debug·원본 비교·C08 current grid·VTune 최종 수집·공식 단일 시간 측정은 +이 단계에서 실행하지 않았다. 변경되지 않은 reference와 tolerance를 사용해야 한다. diff --git a/src/fesa/analysis/hht_dynamic_algebra.cpp b/src/fesa/analysis/hht_dynamic_algebra.cpp index 22b02f6..0332161 100644 --- a/src/fesa/analysis/hht_dynamic_algebra.cpp +++ b/src/fesa/analysis/hht_dynamic_algebra.cpp @@ -2,9 +2,11 @@ #include #include #include +#include #include #include #include +#include #include #include #include @@ -13,6 +15,7 @@ #include "fesa/math/matrix.h" #include "fesa/solvers/linear/linear_solver.h" #include "fesa/solvers/linear/tangent_solver.h" +#include "math/dense_blas_internal.h" namespace fesa { namespace hht_internal { @@ -132,13 +135,13 @@ Status CheckLocalMass(const Domain& domain, const SparseMatrix& matrix, const std::size_t begin, const std::size_t end) { if (begin == end) return Status::Ok(); const std::size_t size = end - begin; - std::array, 6> block{}; - std::array, 6> lower{}; + const dense_blas_internal::ScopedLocalDenseThreads local_threads{}; + Matrix block{size, size}; double scale = 0.0; for (std::size_t i = 0; i < size; ++i) { for (std::size_t j = 0; j < size; ++j) { - block[i][j] = Entry(matrix, begin + i, begin + j); - scale = std::max(scale, std::abs(block[i][j])); + block(i, j) = Entry(matrix, begin + i, begin + j); + scale = std::max(scale, std::abs(block(i, j))); } } if (!std::isfinite(scale) || !(scale > 0.0)) { @@ -149,36 +152,34 @@ Status CheckLocalMass(const Domain& domain, const SparseMatrix& matrix, for (std::size_t i = 0; i < size; ++i) { double sum = 0.0; for (std::size_t j = 0; j < size; ++j) { - block[i][j] /= scale; - sum += std::abs(block[i][j]); + block(i, j) /= scale; + sum += std::abs(block(i, j)); } norm = std::max(norm, sum); - for (std::size_t j = 0; j <= i; ++j) { - double value = block[i][j]; - for (std::size_t k = 0; k < j; ++k) value -= lower[i][k] * lower[j][k]; - if (!std::isfinite(value) || (i == j && !(value > 0.0))) { - return Failure( - domain, "dynamic-positive-mass", - "Physical mass requires strictly positive Cholesky pivots."); - } - lower[i][j] = i == j ? std::sqrt(value) : value / lower[j][j]; - } } + Matrix lower{size, size}; + Matrix inverse_lower{size, size}; + auto factor_failure = dense_blas_internal::CholeskyFailure::kNone; + try { + if (!dense_blas_internal::CholeskyLowerAndInverse( + block, lower, inverse_lower, &factor_failure)) { + if (factor_failure == dense_blas_internal::CholeskyFailure::kInverse) + return Failure( + domain, "dynamic-positive-mass-conditioning", + "A nodal free physical block reciprocal condition is below 2^-32."); + return Failure( + domain, "dynamic-positive-mass", + "Physical mass requires strictly positive Cholesky pivots."); + } + } catch (const std::exception& error) { + return Failure(domain, "dynamic-positive-mass", error.what()); + } + const Matrix inverse = + dense_blas_internal::MultiplyTransposeLeft(inverse_lower, inverse_lower); std::array sums{}; for (std::size_t column = 0; column < size; ++column) { - std::array solution{}; for (std::size_t i = 0; i < size; ++i) { - double value = i == column ? 1.0 : 0.0; - for (std::size_t k = 0; k < i; ++k) value -= lower[i][k] * solution[k]; - solution[i] = value / lower[i][i]; - } - for (std::size_t reverse = size; reverse > 0U; --reverse) { - const std::size_t i = reverse - 1U; - double value = solution[i]; - for (std::size_t k = i + 1U; k < size; ++k) - value -= lower[k][i] * solution[k]; - solution[i] = value / lower[i][i]; - sums[i] += std::abs(solution[i]); + sums[i] += std::abs(inverse(i, column)); } } const double inverse_norm = *std::max_element(sums.begin(), sums.end()); @@ -212,19 +213,33 @@ Result BuildPhysicalSpace( double symmetry_error = 0.0; double null_error = 0.0; if (mass_scale > 0.0) { + Vector normalized_mass{mass.Values().size()}; + Vector null_coefficients{mass.Values().size()}; + for (std::size_t index = 0U; index < mass.Values().size(); ++index) { + normalized_mass[index] = mass.Values()[index] / mass_scale; + const auto full = dofs.FreeDofs()[mass.ColumnIndices()[index]]; + if (full % 6U >= 3U && node_null[full / 6U] != no_null) + null_coefficients[index] = + assembly.rotations[full / 6U].drill_null[full % 6U - 3U]; + } + Vector null_products{0U}; + { + const dense_blas_internal::ScopedLocalDenseThreads local_threads{}; + null_products = dense_blas_internal::MultiplyDiagonal(normalized_mass, + null_coefficients); + } for (std::size_t row = 0; row < mass.Rows(); ++row) { std::vector action(null_nodes.size(), 0.0); for (std::size_t index = mass.RowOffsets()[row]; index < mass.RowOffsets()[row + 1U]; ++index) { const std::size_t column = mass.ColumnIndices()[index]; - const double entry = mass.Values()[index] / mass_scale; + const double entry = normalized_mass[index]; mass_norm = std::hypot(mass_norm, entry); symmetry_error = std::hypot( symmetry_error, entry - Entry(mass, column, row) / mass_scale); const std::size_t full = dofs.FreeDofs()[column]; if (full % 6U >= 3U && node_null[full / 6U] != no_null) { - action[node_null[full / 6U]] += - entry * assembly.rotations[full / 6U].drill_null[full % 6U - 3U]; + action[node_null[full / 6U]] += null_products[index]; } } for (const double value : action) @@ -291,6 +306,7 @@ Result BuildPhysicalSpace( } } std::vector entries; + std::vector> mass_factors; node_starts.push_back(dimension); double complement_gram_error = 0.0; double null_gram_error = 0.0; @@ -350,8 +366,10 @@ Result BuildPhysicalSpace( domain, "dynamic-invalid-analytic-null", "The scaled physical and analytic null Gram identities failed.")); } - for (std::size_t i = 0; i < dimension; ++i) + for (std::size_t i = 0; i < dimension; ++i) { entries.push_back({i, i, 0.0, 0U, entries.size()}); + mass_factors.push_back({0.0, 0.0}); + } for (std::size_t row = 0; row < mass.Rows(); ++row) { const auto& left = rows[dofs.FreeDofs()[row]]; for (std::size_t index = mass.RowOffsets()[row]; @@ -359,13 +377,54 @@ Result BuildPhysicalSpace( const auto& right = rows[dofs.FreeDofs()[mass.ColumnIndices()[index]]]; for (const auto& a : left) { for (const auto& b : right) { - entries.push_back({a.first, b.first, - a.second * mass.Values()[index] * b.second, 0U, - entries.size()}); + entries.push_back({a.first, b.first, a.second, 0U, entries.size()}); + mass_factors.push_back({mass.Values()[index], b.second}); } } } } + Vector coefficients{entries.size()}; + Vector mass_values{entries.size()}; + Vector right_values{entries.size()}; + for (std::size_t i = 0U; i < entries.size(); ++i) { + coefficients[i] = entries[i].value; + mass_values[i] = mass_factors[i].first; + right_values[i] = mass_factors[i].second; + } + Vector projected{0U}; + { + const dense_blas_internal::ScopedLocalDenseThreads local_threads{}; + projected = + dense_blas_internal::MultiplyDiagonal(coefficients, mass_values); + bool finite = true; + for (std::size_t i = 0U; i < projected.Size(); ++i) { + entries[i].value = projected[i]; + finite = finite && std::isfinite(projected[i]); + } + // A rejected intermediate cannot produce a finite original contribution. + // Still compute the remaining finite entries so an earlier final-product + // overflow keeps its diagnostic priority. No invalid operand reaches MKL. + if (!finite) { + std::vector valid_positions; + for (std::size_t i = 0U; i < projected.Size(); ++i) + if (std::isfinite(projected[i])) valid_positions.push_back(i); + Vector valid_left{valid_positions.size()}; + Vector valid_right{valid_positions.size()}; + for (std::size_t i = 0U; i < valid_positions.size(); ++i) { + valid_left[i] = projected[valid_positions[i]]; + valid_right[i] = right_values[valid_positions[i]]; + } + const Vector valid_products = + dense_blas_internal::MultiplyDiagonal(valid_left, valid_right); + for (std::size_t i = 0U; i < valid_positions.size(); ++i) + entries[valid_positions[i]].value = valid_products[i]; + auto rejected = MakeSparse(dimension, dimension, std::move(entries)); + return Result::Failure(rejected.GetStatus()); + } + projected = dense_blas_internal::MultiplyDiagonal(projected, right_values); + } + for (std::size_t i = 0U; i < entries.size(); ++i) + entries[i].value = projected[i]; auto block = MakeSparse(dimension, dimension, std::move(entries)); if (!block.HasValue()) return Result::Failure(block.GetStatus()); @@ -436,15 +495,20 @@ Result CheckRegularity(const Domain& domain, for (std::size_t i = 0; i < nodes.size(); ++i) { node_null[nodes[i]] = i; const auto& n = assembly.rotations[nodes[i]].drill_null; - double stiffness = 0.0; + Matrix metric{3U, 3U}; + Vector direction{3U}; for (std::size_t a = 0; a < 3; ++a) { + direction[a] = n[a]; for (std::size_t b = 0; b < 3; ++b) { - stiffness += n[a] * - Entry(assembly.drilling_metric, 6U * nodes[i] + 3U + a, - 6U * nodes[i] + 3U + b) * - n[b]; + metric(a, b) = Entry(assembly.drilling_metric, 6U * nodes[i] + 3U + a, + 6U * nodes[i] + 3U + b); } } + double stiffness = 0.0; + { + const dense_blas_internal::ScopedLocalDenseThreads local_threads{}; + stiffness = direction.Dot(metric.Multiply(direction)); + } if (!std::isfinite(stiffness) || !(stiffness > 0.0)) { return Result::Failure( Failure(domain, "dynamic-algebraic-regularity", @@ -454,22 +518,37 @@ Result CheckRegularity(const Domain& domain, weights[i] = std::sqrt(stiffness); } std::vector entries; + std::vector null_coefficients; for (std::size_t row = 0; row < nodes.size(); ++row) { entries.push_back({row, row, 0.0, 0U, entries.size()}); + null_coefficients.push_back(0.0); for (std::size_t index = gradient.RowOffsets()[row]; index < gradient.RowOffsets()[row + 1U]; ++index) { const std::size_t full = gradient.ColumnIndices()[index]; if (full % 6U >= 3U && node_null[full / 6U] != no_null) { const std::size_t column = node_null[full / 6U]; entries.push_back( - {row, column, - gradient.Values()[index] * - assembly.rotations[full / 6U].drill_null[full % 6U - 3U] / - weights[row] / weights[column], - 0U, entries.size()}); + {row, column, gradient.Values()[index], 0U, entries.size()}); + null_coefficients.push_back( + assembly.rotations[full / 6U].drill_null[full % 6U - 3U]); } } } + Vector gradient_values{entries.size()}; + Vector null_values{entries.size()}; + for (std::size_t i = 0U; i < entries.size(); ++i) { + gradient_values[i] = entries[i].value; + null_values[i] = null_coefficients[i]; + } + Vector products{0U}; + { + const dense_blas_internal::ScopedLocalDenseThreads local_threads{}; + products = + dense_blas_internal::MultiplyDiagonal(gradient_values, null_values); + } + for (std::size_t i = 0U; i < entries.size(); ++i) + entries[i].value = + products[i] / weights[entries[i].row] / weights[entries[i].column]; auto matrix = MakeSparse(nodes.size(), nodes.size(), std::move(entries)); if (!matrix.HasValue()) return Result::Failure(matrix.GetStatus()); diff --git a/src/fesa/analysis/hht_reduced_operators.cpp b/src/fesa/analysis/hht_reduced_operators.cpp index 3a8f1a7..c484f2e 100644 --- a/src/fesa/analysis/hht_reduced_operators.cpp +++ b/src/fesa/analysis/hht_reduced_operators.cpp @@ -6,10 +6,12 @@ #include #include #include +#include #include #include #include "fesa/solvers/linear/tangent_solver.h" +#include "math/dense_blas_internal.h" #include "math/high_low_arithmetic.h" namespace fesa::hht_reduced_internal { @@ -614,10 +616,33 @@ Matrix PhysicalJacobian(const ConstraintOperators& operators, Vector ContractConstraintHessian(const ConstraintOperators& operators, const Vector& left, const Vector& right) { + const dense_blas_internal::ScopedLocalDenseThreads local_threads{}; Vector result{operators.rows.size()}; - for (std::size_t row = 0U; row < operators.rows.size(); ++row) - for (const auto& entry : operators.rows[row].hessian) - result[row] += entry.value * left[entry.row] * right[entry.column]; + for (std::size_t row = 0U; row < operators.rows.size(); ++row) { + const auto& hessian = operators.rows[row].hessian; + Vector coefficients{hessian.size()}; + Vector left_values{hessian.size()}; + Vector right_values{hessian.size()}; + for (std::size_t i = 0U; i < hessian.size(); ++i) { + coefficients[i] = hessian[i].value; + left_values[i] = left[hessian[i].row]; + right_values[i] = right[hessian[i].column]; + } + try { + const Vector first = + dense_blas_internal::MultiplyDiagonal(coefficients, left_values); + const Vector products = + dense_blas_internal::MultiplyDiagonal(first, right_values); + // Retain the original entry order and left-associated products; only + // ordinary products move to MKL. This is failure-only replay evidence. + for (std::size_t i = 0U; i < products.Size(); ++i) + result[row] += products[i]; + } catch (const std::invalid_argument&) { + // A nonfinite input/intermediate cannot yield a finite original sum. + // Keep the existing caller's finite failure without a scalar fallback. + result[row] = std::numeric_limits::quiet_NaN(); + } + } return result; } @@ -1186,13 +1211,19 @@ Result CheckGraphRegularity( Failure("Reduced graph regularity node is invalid.", "dynamic-mass-regularity")); const auto& null = assembly.rotations[node].drill_null; - double metric = 0.0; - for (std::size_t i = 0U; i < 3U; ++i) + Matrix metric_block{3U, 3U}; + Vector direction{3U}; + for (std::size_t i = 0U; i < 3U; ++i) { + direction[i] = null[i]; for (std::size_t j = 0U; j < 3U; ++j) - metric += null[i] * - SparseEntry(assembly.drilling_metric, 6U * node + 3U + i, - 6U * node + 3U + j) * - null[j]; + metric_block(i, j) = SparseEntry( + assembly.drilling_metric, 6U * node + 3U + i, 6U * node + 3U + j); + } + double metric = 0.0; + { + const dense_blas_internal::ScopedLocalDenseThreads local_threads{}; + metric = direction.Dot(metric_block.Multiply(direction)); + } if (!std::isfinite(metric) || !(metric > 0.0)) return Result::Failure(Failure( "Reduced graph requires a positive finite null drilling metric.", diff --git a/src/fesa/elements/mitc4_shell.cpp b/src/fesa/elements/mitc4_shell.cpp index a4f3aa5..129f05b 100644 --- a/src/fesa/elements/mitc4_shell.cpp +++ b/src/fesa/elements/mitc4_shell.cpp @@ -8,6 +8,7 @@ #include #include "fesa/model/model_types.h" +#include "math/dense_blas_internal.h" namespace fesa { @@ -223,28 +224,18 @@ Matrix ScaledMatrix(const Matrix& source, double factor) { return result; } -/// @brief Lifts a local matrix by a deterministic serial congruence product. -/// @note The loop and accumulation order are part of reproducible stiffness. +/// @brief Lifts a local matrix through MKL and mirrors the upper triangle. Matrix Congruence(const Matrix& local, const Matrix& transformation) { if (local.Rows() != local.Columns() || local.Rows() != transformation.Rows()) { throw std::invalid_argument{ "MITC4 congruence dimensions are incompatible."}; } - Matrix result{transformation.Columns(), transformation.Columns()}; + Matrix result = dense_blas_internal::MultiplyTransposeLeft( + transformation, local.Multiply(transformation)); for (std::size_t row = 0U; row < result.Rows(); ++row) { for (std::size_t column = row; column < result.Columns(); ++column) { - double value = 0.0; - for (std::size_t local_row = 0U; local_row < local.Rows(); ++local_row) { - for (std::size_t local_column = 0U; local_column < local.Columns(); - ++local_column) { - value += transformation(local_row, row) * - local(local_row, local_column) * - transformation(local_column, column); - } - } - result(row, column) = value; - result(column, row) = value; + result(column, row) = result(row, column); } } return result; @@ -611,7 +602,8 @@ Matrix Mitc4Shell::TransverseShearSectionMatrix() const { return result; } -Result Mitc4Shell::Stiffness() const { +Result Mitc4Shell::Stiffness() const try { + const dense_blas_internal::ScopedLocalDenseThreads local_threads{}; const Matrix constitutive = MaterialConstitutive5(); const Matrix tying_samples = CovariantTyingShearSamples20(); Matrix physical_local{kPhysicalDofCount, kPhysicalDofCount}; @@ -627,17 +619,12 @@ Result Mitc4Shell::Stiffness() const { const Matrix strain = StrainDisplacement( point.natural_coordinates[0], point.natural_coordinates[1], point.natural_coordinates[2], &tying_samples); + const Matrix integrand = dense_blas_internal::MultiplyTransposeLeft( + strain, constitutive.Multiply(strain)); for (std::size_t row = 0U; row < kPhysicalDofCount; ++row) { for (std::size_t column = row; column < kPhysicalDofCount; ++column) { - double integrand = 0.0; - for (std::size_t first = 0U; first < 5U; ++first) { - for (std::size_t second = 0U; second < 5U; ++second) { - integrand += strain(first, row) * constitutive(first, second) * - strain(second, column); - } - } const double contribution = - integrand * geometry.jacobian * point.weight; + integrand(row, column) * geometry.jacobian * point.weight; physical_local(row, column) += contribution; if (row != column) { physical_local(column, row) += contribution; @@ -714,6 +701,8 @@ Result Mitc4Shell::Stiffness() const { Mitc4Stiffness{std::move(physical_local), std::move(physical_global), std::move(drilling_global), std::move(stabilized_global), drilling_stiffness}); +} catch (const std::exception& error) { + return StiffnessFailure(source_location_, identity_, error.what()); } Result Mitc4Shell::Mass(const double density) const { diff --git a/src/fesa/elements/nonlinear_mitc4_inertia.cpp b/src/fesa/elements/nonlinear_mitc4_inertia.cpp index 1fba4fb..13198fb 100644 --- a/src/fesa/elements/nonlinear_mitc4_inertia.cpp +++ b/src/fesa/elements/nonlinear_mitc4_inertia.cpp @@ -4,11 +4,13 @@ #include #include #include +#include #include #include #include #include "fesa/elements/finite_rotation.h" +#include "math/dense_blas_internal.h" namespace fesa { namespace { @@ -35,6 +37,16 @@ class InertiaSum { double correction_ = 0.0; }; +/// @brief Contracts two nodal 6-by-3 coefficient blocks through dense MKL. +/// @note Nonfinite coefficients remain failures of the final inertia gate. +void DotBlocks(const std::array& left, + const std::array& right, + std::array& result) { + if (!dense_blas_internal::MultiplyTransposeRight( + left.data(), 6U, right.data(), 6U, 3U, result.data())) + result.fill(std::numeric_limits::quiet_NaN()); +} + Vector3 WeightedSum(const std::array& weights, const std::array& values) { Vector3 result{}; @@ -442,12 +454,11 @@ Status CheckPhysicalMass(const Matrix& mass, } double gram_error = 0.0; double orthogonality_error = 0.0; + const Matrix gram = + dense_blas_internal::MultiplyTransposeLeft(complement, complement); for (std::size_t row = 0U; row < kPhysical; ++row) { for (std::size_t column = 0U; column < kPhysical; ++column) { - double entry = row == column ? -1.0 : 0.0; - for (std::size_t inner = 0U; inner < kDofs; ++inner) { - entry += complement(inner, row) * complement(inner, column); - } + const double entry = gram(row, column) - (row == column ? 1.0 : 0.0); gram_error = std::hypot(gram_error, entry); } for (std::size_t node = 0U; node < kNodes; ++node) { @@ -472,55 +483,42 @@ Status CheckPhysicalMass(const Matrix& mass, return fail("invalid-nonlinear-inertia-null", "Scaled mass symmetry and analytic null identities failed."); } - Matrix block{kPhysical, kPhysical}; - Matrix mass_complement{kDofs, kPhysical}; - for (std::size_t row = 0U; row < kDofs; ++row) { - for (std::size_t column = 0U; column < kPhysical; ++column) { - for (std::size_t inner = 0U; inner < kDofs; ++inner) { - mass_complement(row, column) += - scaled(row, inner) * complement(inner, column); - } - } - } + const Matrix mass_complement = scaled.Multiply(complement); + const Matrix block = + dense_blas_internal::MultiplyTransposeLeft(complement, mass_complement); double block_norm = 0.0; for (std::size_t row = 0U; row < kPhysical; ++row) { double row_sum = 0.0; for (std::size_t column = 0U; column < kPhysical; ++column) { - for (std::size_t inner = 0U; inner < kDofs; ++inner) { - block(row, column) += - complement(inner, row) * mass_complement(inner, column); - } row_sum += std::abs(block(row, column)); } block_norm = std::max(block_norm, row_sum); } Matrix lower{kPhysical, kPhysical}; - for (std::size_t row = 0U; row < kPhysical; ++row) { - for (std::size_t column = 0U; column <= row; ++column) { - double entry = block(row, column); - for (std::size_t inner = 0U; inner < column; ++inner) { - entry -= lower(row, inner) * lower(column, inner); - } - if (!std::isfinite(entry) || (row == column && !(entry > 0.0))) { - return fail( - "invalid-nonlinear-inertia-positive-mass", - "Physical mass requires strictly positive Cholesky pivots."); - } - lower(row, column) = - row == column ? std::sqrt(entry) : entry / lower(column, column); - } + Matrix inverse_lower{kPhysical, kPhysical}; + auto factor_failure = dense_blas_internal::CholeskyFailure::kNone; + if (!dense_blas_internal::CholeskyLowerAndInverse(block, lower, inverse_lower, + &factor_failure)) { + if (factor_failure == dense_blas_internal::CholeskyFailure::kInverse) + return fail( + "ill-conditioned-nonlinear-inertia-positive-mass", + "Physical mass reciprocal infinity condition is below 2^-32."); + return fail("invalid-nonlinear-inertia-positive-mass", + "Physical mass requires strictly positive Cholesky pivots."); } std::array inverse_row_sums{}; double reconstruction_error = 0.0; double block_frobenius = 0.0; + Matrix transpose{kPhysical, kPhysical}; + for (std::size_t row = 0U; row < kPhysical; ++row) + for (std::size_t column = 0U; column < kPhysical; ++column) + transpose(row, column) = lower(column, row); + const Matrix reconstructed = lower.Multiply(transpose); for (std::size_t row = 0U; row < kPhysical; ++row) { for (std::size_t column = 0U; column < kPhysical; ++column) { - double reconstructed = 0.0; - for (std::size_t inner = 0U; inner <= std::min(row, column); ++inner) { - reconstructed += lower(row, inner) * lower(column, inner); - } reconstruction_error = - std::hypot(reconstruction_error, reconstructed - block(row, column)); + std::hypot(reconstruction_error, + reconstructed(row, column) - block(row, column)); block_frobenius = std::hypot(block_frobenius, block(row, column)); } } @@ -528,25 +526,11 @@ Status CheckPhysicalMass(const Matrix& mass, return fail("invalid-nonlinear-inertia-positive-mass", "Physical mass Cholesky reconstruction failed."); } - for (std::size_t column = 0U; column < kPhysical; ++column) { - std::array solution{}; - for (std::size_t row = 0U; row < kPhysical; ++row) { - double rhs = row == column ? 1.0 : 0.0; - for (std::size_t inner = 0U; inner < row; ++inner) { - rhs -= lower(row, inner) * solution[inner]; - } - solution[row] = rhs / lower(row, row); - } - for (std::size_t reverse = kPhysical; reverse > 0U; --reverse) { - const std::size_t row = reverse - 1U; - double rhs = solution[row]; - for (std::size_t inner = row + 1U; inner < kPhysical; ++inner) { - rhs -= lower(inner, row) * solution[inner]; - } - solution[row] = rhs / lower(row, row); - inverse_row_sums[row] += std::abs(solution[row]); - } - } + const Matrix inverse = + dense_blas_internal::MultiplyTransposeLeft(inverse_lower, inverse_lower); + for (std::size_t column = 0U; column < kPhysical; ++column) + for (std::size_t row = 0U; row < kPhysical; ++row) + inverse_row_sums[row] += std::abs(inverse(row, column)); const double inverse_norm = *std::max_element(inverse_row_sums.begin(), inverse_row_sums.end()); const double reciprocal_condition = 1.0 / block_norm / inverse_norm; @@ -712,6 +696,7 @@ Result NonlinearMitc4Inertia::Create( Result NonlinearMitc4Inertia::Evaluate( const NonlinearMitc4TrialDofs& trial, const Vector& velocity, const Vector& acceleration) const { + const dense_blas_internal::ScopedLocalDenseThreads local_threads{}; if (velocity.Size() != kDofs || acceleration.Size() != kDofs || !IsFinite(velocity) || !IsFinite(acceleration)) { return Failure( @@ -899,6 +884,42 @@ Result NonlinearMitc4Inertia::Evaluate( result.momentum_rate_scale += scale * acceleration_point.Norm(); result.angular_momentum_rate_scale += scale * position.Cross(acceleration_point).Norm(); + std::array h_block; + std::array velocity_block; + std::array position_block; + for (std::size_t local = 0U; local < 6U; ++local) { + const auto field = 6U * mass_node + local; + for (std::size_t axis = 0U; axis < 3U; ++axis) { + h_block[3U * local + axis] = h[field][axis]; + velocity_block[3U * local + axis] = + velocity_derivative[field][axis]; + position_block[3U * local + axis] = + position_derivative[field][axis]; + } + } + std::array mass_block; + std::array velocity_product; + std::array position_product; + DotBlocks(h_block, h_block, mass_block); + DotBlocks(h_block, velocity_block, velocity_product); + DotBlocks(h_block, position_block, position_product); + std::array, 3U> mass_variation; + for (std::size_t direction = 0U; direction < 3U; ++direction) { + std::array derivative_block; + for (std::size_t local = 0U; local < 6U; ++local) { + const Vector3 derivative = + local < 3U ? Vector3{} + : z * rotations[mass_node] + .director_second[local - 3U][direction]; + for (std::size_t axis = 0U; axis < 3U; ++axis) + derivative_block[3U * local + axis] = derivative[axis]; + } + std::array reverse; + DotBlocks(derivative_block, h_block, mass_variation[direction]); + DotBlocks(h_block, derivative_block, reverse); + for (std::size_t index = 0U; index < 36U; ++index) + mass_variation[direction][index] += reverse[index]; + } for (std::size_t row = 6U * mass_node; row < 6U * (mass_node + 1U); ++row) { // Keep c independent of a; subtracting M*a destroys small c values. @@ -906,10 +927,11 @@ Result NonlinearMitc4Inertia::Evaluate( scale * h[row].Dot(convective_point); for (std::size_t column = 6U * mass_node; column < 6U * (mass_node + 1U); ++column) { - result.mass(row, column) += scale * h[row].Dot(h[column]); + const auto local_index = 6U * (row % 6U) + column % 6U; + result.mass(row, column) += scale * mass_block[local_index]; result.velocity_tangent(row, column) += - scale * h[row].Dot(velocity_derivative[column]); - double position_entry = h[row].Dot(position_derivative[column]); + scale * velocity_product[local_index]; + double position_entry = position_product[local_index]; if (row / 6U == column / 6U && row % 6U >= 3U && column % 6U >= 3U) { const std::size_t node = row / 6U; @@ -922,20 +944,10 @@ Result NonlinearMitc4Inertia::Evaluate( result.position_tangent(row, column) += scale * position_entry; if (column % 6U >= 3U) { const auto direction = column % 6U - 3U; - const Vector3 dh_row = - row % 6U < 3U - ? Vector3{} - : z * rotations[mass_node] - .director_second[row % 6U - 3U][direction]; for (std::size_t local = 0U; local < 6U; ++local) { - const auto field = 6U * mass_node + local; - const Vector3 dh_field = - local < 3U - ? Vector3{} - : z * rotations[mass_node] - .director_second[local - 3U][direction]; mass_derivative[mass_node][direction][row % 6U][local] += - scale * (dh_row.Dot(h[field]) + h[row].Dot(dh_field)); + scale * + mass_variation[direction][6U * (row % 6U) + local]; } } } @@ -1004,8 +1016,15 @@ Result NonlinearMitc4Inertia::Evaluate( "invalid-nonlinear-inertia-evaluation", source_location_, identity_, "Nonlinear MITC4 inertia evaluation must remain finite."); } - const Status physical_mass = CheckPhysicalMass(result.mass, nulls, length, - source_location_, identity_); + Status physical_mass = Status::Ok(); + try { + physical_mass = CheckPhysicalMass(result.mass, nulls, length, + source_location_, identity_); + } catch (const std::exception& error) { + return Failure( + "invalid-nonlinear-inertia-positive-mass", source_location_, identity_, + error.what()); + } if (!physical_mass.IsOk()) { return Result::Failure(physical_mass); } diff --git a/src/fesa/elements/nonlinear_mitc4_shell.cpp b/src/fesa/elements/nonlinear_mitc4_shell.cpp index f61ac61..3dd5baf 100644 --- a/src/fesa/elements/nonlinear_mitc4_shell.cpp +++ b/src/fesa/elements/nonlinear_mitc4_shell.cpp @@ -10,6 +10,7 @@ #include #include "fesa/elements/finite_rotation.h" +#include "math/dense_blas_internal.h" namespace fesa { namespace { @@ -394,22 +395,63 @@ constexpr std::size_t kHessianSize = kGlobalDofCount * kGlobalDofCount; /// @brief Carries exact first and second trial-coordinate variations. class Hessian { public: - Hessian() : values_{std::make_shared>()} {} + Hessian() = default; - double& operator[](const std::size_t index) { return (*values_)[index]; } + double& operator[](const std::size_t index) { return Data()[index]; } const double& operator[](const std::size_t index) const { - return (*values_)[index]; + return Storage()[index]; } - auto begin() const { return values_->begin(); } + auto begin() const { return Storage().begin(); } - auto end() const { return values_->end(); } + auto end() const { return Storage().end(); } + + double* Data() { + if (!values_) { + values_ = std::make_shared(); + values_->entries.fill(0.0); + } + return values_->entries.data(); + } + + const double* ReadData() const { return Storage().data(); } + + /// @brief Gets fresh output storage without redundant zero initialization. + /// @note Only fresh operator results use this pointer; all 576 entries must + /// be written before any read or copy of the result. + double* OverwriteData() { + if (!values_) values_ = std::make_shared(); + return values_->entries.data(); + } private: - std::shared_ptr> values_; + struct Values { + // User-provided construction leaves entries uninitialized for a complete + // output overwrite while make_shared retains the single allocation. + Values() noexcept {} + std::array entries; + }; + + const std::array& Storage() const { + // Constants borrow only this immutable zero. Mutable coefficients retain + // the existing shared ownership; no trial-dependent cache is introduced. + static const std::array kZero{}; + return values_ ? values_->entries : kZero; + } + + std::shared_ptr values_; }; +/// @brief Forms an ordinary derivative outer product without a Matrix copy. +/// @note Rejected nonfinite operands propagate to the existing Jet finite gate. +void Outer(const std::array& left, + const std::array& right, double* result) { + if (!dense_blas_internal::OuterProduct(left.data(), left.size(), right.data(), + right.size(), result)) + std::fill_n(result, kHessianSize, std::numeric_limits::quiet_NaN()); +} + /// @brief Carries exact first and second trial-coordinate variations. struct Jet { double value{}; @@ -455,8 +497,11 @@ Jet operator+(const Jet& left, const Jet& right) { result.mixed_second[index] = left.mixed_second[index] + right.mixed_second[index]; } + const double* left_second = left.second.ReadData(); + const double* right_second = right.second.ReadData(); + double* result_second = result.second.OverwriteData(); for (std::size_t index = 0U; index < kHessianSize; ++index) { - result.second[index] = left.second[index] + right.second[index]; + result_second[index] = left_second[index] + right_second[index]; } return result; } @@ -475,8 +520,11 @@ Jet operator-(const Jet& left, const Jet& right) { result.mixed_second[index] = left.mixed_second[index] - right.mixed_second[index]; } + const double* left_second = left.second.ReadData(); + const double* right_second = right.second.ReadData(); + double* result_second = result.second.OverwriteData(); for (std::size_t index = 0U; index < kHessianSize; ++index) { - result.second[index] = left.second[index] - right.second[index]; + result_second[index] = left_second[index] - right_second[index]; } return result; } @@ -506,13 +554,18 @@ Jet operator*(const Jet& left, const Jet& right) { 2.0 * left.directional_first * right.mixed_first[index] + left.value * right.mixed_second[index]; } + const double* left_second = left.second.ReadData(); + const double* right_second = right.second.ReadData(); + double* result_second = result.second.OverwriteData(); + Outer(left.first, right.first, result_second); + std::array reverse; + Outer(right.first, left.first, reverse.data()); for (std::size_t row = 0U; row < kGlobalDofCount; ++row) { for (std::size_t column = 0U; column < kGlobalDofCount; ++column) { const std::size_t index = row * kGlobalDofCount + column; - result.second[index] = left.second[index] * right.value + - left.first[row] * right.first[column] + - left.first[column] * right.first[row] + - left.value * right.second[index]; + result_second[index] = left_second[index] * right.value + + result_second[index] + reverse[index] + + left.value * right_second[index]; } } return result; @@ -529,8 +582,10 @@ Jet operator*(const double scalar, const Jet& value) { result.mixed_first[index] = scalar * value.mixed_first[index]; result.mixed_second[index] = scalar * value.mixed_second[index]; } + const double* value_second = value.second.ReadData(); + double* result_second = result.second.OverwriteData(); for (std::size_t index = 0U; index < kHessianSize; ++index) { - result.second[index] = scalar * value.second[index]; + result_second[index] = scalar * value_second[index]; } return result; } @@ -559,11 +614,16 @@ Jet Compose(const Jet& input, const double value, const double first, input.first[index] * input.directional_second) + first * input.mixed_second[index]; } + std::array scaled_first{}; + for (std::size_t row = 0U; row < kGlobalDofCount; ++row) + scaled_first[row] = second * input.first[row]; + const double* input_second = input.second.ReadData(); + double* result_second = result.second.OverwriteData(); + Outer(scaled_first, input.first, result_second); for (std::size_t row = 0U; row < kGlobalDofCount; ++row) { for (std::size_t column = 0U; column < kGlobalDofCount; ++column) { const std::size_t index = row * kGlobalDofCount + column; - result.second[index] = first * input.second[index] + - second * input.first[row] * input.first[column]; + result_second[index] = first * input_second[index] + result_second[index]; } } return result; @@ -613,19 +673,32 @@ Jet Atan2(const Jet& y, const Jet& x) { const double derivative_yx = (y.value * y.value - x.value * x.value) / radius_fourth; const double derivative_xx = 2.0 * x.value * y.value / radius_fourth; + std::array scaled_y{}; + std::array scaled_x{}; + for (std::size_t i = 0U; i < kGlobalDofCount; ++i) { + scaled_y[i] = derivative_yy * y.first[i]; + scaled_x[i] = derivative_xx * x.first[i]; + } + std::array, 4U> products; + Outer(scaled_y, y.first, products[0U].data()); + Outer(y.first, x.first, products[1U].data()); + Outer(x.first, y.first, products[2U].data()); + Outer(scaled_x, x.first, products[3U].data()); for (std::size_t index = 0U; index < kGlobalDofCount; ++index) { result.first[index] = derivative_y * y.first[index] + derivative_x * x.first[index]; } + const double* y_second = y.second.ReadData(); + const double* x_second = x.second.ReadData(); + double* result_second = result.second.OverwriteData(); for (std::size_t row = 0U; row < kGlobalDofCount; ++row) { for (std::size_t column = 0U; column < kGlobalDofCount; ++column) { const std::size_t index = row * kGlobalDofCount + column; - result.second[index] = derivative_y * y.second[index] + - derivative_x * x.second[index] + - derivative_yy * y.first[row] * y.first[column] + - derivative_yx * (y.first[row] * x.first[column] + - x.first[row] * y.first[column]) + - derivative_xx * x.first[row] * x.first[column]; + result_second[index] = + derivative_y * y_second[index] + derivative_x * x_second[index] + + products[0U][index] + + derivative_yx * (products[1U][index] + products[2U][index]) + + products[3U][index]; } } // Use a bounded ratio for the directional third derivative only. Retain @@ -1181,10 +1254,14 @@ Jet DrillingEnergy(const std::array& coordinates, DrillSpinRate(trial.rotation[node], initial_directors[node]) * gamma; } energy = energy + 0.5 * drilling_stiffness * gamma * gamma; + std::array scaled_gradient{}; + for (std::size_t row = 0U; row < kGlobalDofCount; ++row) + scaled_gradient[row] = drilling_stiffness * gamma.first[row]; + std::array contribution; + Outer(scaled_gradient, gamma.first, contribution.data()); for (std::size_t row = 0U; row < kGlobalDofCount; ++row) { for (std::size_t column = 0U; column < kGlobalDofCount; ++column) { - (*metric)(row, column) += - drilling_stiffness * gamma.first[row] * gamma.first[column]; + (*metric)(row, column) += contribution[row * kGlobalDofCount + column]; } } } @@ -1505,6 +1582,7 @@ NonlinearMitc4Shell::EvaluateExternalDrillingConstraint( const std::size_t node, const Vector3& phi, const Vector3& phi_dot, const Vector3& force, const Vector3& force_rate, const Vector3& force_acceleration) const { + const dense_blas_internal::ScopedLocalDenseThreads local_threads{}; if (node >= kNodeCount || !phi_dot.IsFinite() || !force.IsFinite() || !force_rate.IsFinite() || !force_acceleration.IsFinite()) { return Failure( @@ -1594,6 +1672,9 @@ NonlinearMitc4Shell::EvaluateExternalDrillingConstraint( Result EvaluateInternalResidualAndTangent( const NonlinearMitc4Shell& element, const NonlinearMitc4TrialDofs& trial, const Vector* coordinate_direction) { + // Element-local MKL work uses one worker only when the caller has no TLS + // override. The guard restores that request before returning to global solve. + const dense_blas_internal::ScopedLocalDenseThreads local_threads{}; if (coordinate_direction != nullptr) { if (coordinate_direction->Size() != kGlobalDofCount) { return Failure( diff --git a/src/fesa/math/dense_blas_internal.cpp b/src/fesa/math/dense_blas_internal.cpp index 413a71c..a544715 100644 --- a/src/fesa/math/dense_blas_internal.cpp +++ b/src/fesa/math/dense_blas_internal.cpp @@ -1,5 +1,6 @@ #include "math/dense_blas_internal.h" +#include #include #include #include @@ -26,8 +27,22 @@ ScopedSingleEquationThreads::~ScopedSingleEquationThreads() { if (restore_default_) static_cast(mkl_set_num_threads_local(0)); } +ScopedLocalDenseThreads::ScopedLocalDenseThreads() noexcept { + const int previous = mkl_set_num_threads_local(1); + if (previous == 0) { + restore_default_ = true; + } else { + static_cast(mkl_set_num_threads_local(previous)); + } +} + +ScopedLocalDenseThreads::~ScopedLocalDenseThreads() { + if (restore_default_) static_cast(mkl_set_num_threads_local(0)); +} + bool CholeskyLowerAndInverse(const Matrix& matrix, Matrix& lower, - Matrix& inverse_lower) { + Matrix& inverse_lower, CholeskyFailure* failure) { + if (failure != nullptr) *failure = CholeskyFailure::kNone; if (matrix.Rows() != matrix.Columns()) throw std::invalid_argument{"Cholesky requires a square matrix."}; const auto size = ToMklSize(matrix.Rows()); @@ -48,17 +63,30 @@ bool CholeskyLowerAndInverse(const Matrix& matrix, Matrix& lower, values[row * matrix.Columns() + column] = matrix(row, column); } if (LAPACKE_dpotrf(LAPACK_ROW_MAJOR, 'L', size.Value(), values.data(), - size.Value()) != 0) return false; + size.Value()) != 0) { + if (failure != nullptr) *failure = CholeskyFailure::kPivot; + return false; + } for (std::size_t row = 0; row < matrix.Rows(); ++row) - for (std::size_t column = 0; column <= row; ++column) + for (std::size_t column = 0; column <= row; ++column) { candidate_lower(row, column) = values[row * matrix.Columns() + column]; + if (!std::isfinite(candidate_lower(row, column))) { + if (failure != nullptr) *failure = CholeskyFailure::kPivot; + return false; + } + } if (LAPACKE_dtrtri(LAPACK_ROW_MAJOR, 'L', 'N', size.Value(), values.data(), - size.Value()) != 0) return false; + size.Value()) != 0) { + if (failure != nullptr) *failure = CholeskyFailure::kInverse; + return false; + } for (std::size_t row = 0; row < matrix.Rows(); ++row) for (std::size_t column = 0; column <= row; ++column) { const auto value = values[row * matrix.Columns() + column]; - if (!std::isfinite(candidate_lower(row, column)) || !std::isfinite(value)) + if (!std::isfinite(value)) { + if (failure != nullptr) *failure = CholeskyFailure::kInverse; return false; + } candidate_inverse(row, column) = value; } lower = std::move(candidate_lower); @@ -96,6 +124,41 @@ Matrix MultiplyTransposeLeft(const Matrix& left, const Matrix& right) { return transposed.Multiply(right); } +bool OuterProduct(const double* left, const std::size_t rows, + const double* right, const std::size_t columns, + double* destination) { + return MultiplyTransposeRight(left, rows, right, columns, 1U, destination); +} + +bool MultiplyTransposeRight(const double* left, const std::size_t rows, + const double* right, const std::size_t columns, + const std::size_t inner, double* destination) { + const auto m = ToMklSize(rows); + const auto n = ToMklSize(columns); + const auto k = ToMklSize(inner); + const auto maximum = std::numeric_limits::max() / sizeof(double); + if (!m.HasValue() || !n.HasValue() || !k.HasValue() || + (columns != 0U && rows > maximum / columns) || + (inner != 0U && (rows > maximum / inner || columns > maximum / inner))) + throw std::length_error{"Dense product dimensions exceed backend limits."}; + if (rows == 0U || columns == 0U) return true; + if (destination == nullptr || + (inner != 0U && (left == nullptr || right == nullptr))) + throw std::invalid_argument{"Dense product requires nonnull storage."}; + if (inner == 0U) { + std::fill_n(destination, rows * columns, 0.0); + return true; + } + for (std::size_t i = 0U; i < rows * inner; ++i) + if (!std::isfinite(left[i])) return false; + for (std::size_t i = 0U; i < columns * inner; ++i) + if (!std::isfinite(right[i])) return false; + cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, m.Value(), n.Value(), + k.Value(), 1.0, left, k.Value(), right, k.Value(), 0.0, + destination, n.Value()); + return true; +} + Result ToMklSize(const std::size_t size) { if (size > static_cast((std::numeric_limits::max)())) { return Result::Failure(Status::Failure( diff --git a/src/fesa/math/dense_blas_internal.h b/src/fesa/math/dense_blas_internal.h index 78f9e77..ea3b107 100644 --- a/src/fesa/math/dense_blas_internal.h +++ b/src/fesa/math/dense_blas_internal.h @@ -26,6 +26,20 @@ class ScopedSingleEquationThreads { bool restore_default_{false}; }; +/// @brief Limits element-local dense work under the caller's default policy. +/// @note Positive caller TLS overrides are retained. Default requests are +/// restored on success and unwind; process-wide settings are never changed. +class ScopedLocalDenseThreads { + public: + ScopedLocalDenseThreads() noexcept; + ScopedLocalDenseThreads(const ScopedLocalDenseThreads&) = delete; + ScopedLocalDenseThreads& operator=(const ScopedLocalDenseThreads&) = delete; + ~ScopedLocalDenseThreads(); + + private: + bool restore_default_{false}; +}; + /// @brief Converts a dense storage length to the private MKL integer type. /// @return The converted length or an uncategorized overflow failure. Result ToMklSize(std::size_t size); @@ -39,13 +53,38 @@ Matrix MultiplyTransposeLeft(const Matrix& left, const Matrix& right); /// @throws std::length_error if the size exceeds the backend integer range. Vector MultiplyDiagonal(const Vector& diagonal, const Vector& rhs); +/// @brief Writes one row-major outer product through MKL without allocation. +/// @note Source arrays and destination must not alias. A nonfinite input +/// returns false before writing; finite product overflow remains for the +/// owner's gate. A zero dimension accepts null pointers and performs no backend +/// call. +/// @throws std::invalid_argument for null storage of a nonempty product. +/// @throws std::length_error for dimensions outside the backend integer range. +bool OuterProduct(const double* left, std::size_t rows, const double* right, + std::size_t columns, double* destination); + +/// @brief Writes left times right-transpose to nonaliased contiguous storage. +/// @note Inputs are row-major rows-by-inner and columns-by-inner. Nonfinite +/// input returns false before writing. Finite output overflow is not hidden. +/// Zero output dimensions accept null pointers; inner zero writes exact zeros. +/// @throws std::invalid_argument for missing nonempty storage. +/// @throws std::length_error for dimension or storage-size overflow. +bool MultiplyTransposeRight(const double* left, std::size_t rows, + const double* right, std::size_t columns, + std::size_t inner, double* destination); + +/// @brief Distinguishes positive-factor failure from inverse failure. +enum class CholeskyFailure { kNone, kPivot, kInverse }; + /// @brief Computes a lower Cholesky factor and its inverse with MKL LAPACK. /// @return False for a nonpositive pivot or nonfinite factor/inverse. /// @note Outputs are replaced only after both finite candidates succeed. +/// The optional failure output is initialized even when validation throws. /// @throws std::invalid_argument for a nonsquare or nonfinite input. /// @throws std::length_error for a backend size overflow. bool CholeskyLowerAndInverse(const Matrix& matrix, Matrix& lower, - Matrix& inverse_lower); + Matrix& inverse_lower, + CholeskyFailure* failure = nullptr); /// @brief Copies contiguous double values through the private BLAS backend. /// @param source Readable storage for size values when size is nonzero. diff --git a/tests/unit/math/dense_blas_internal_test.cpp b/tests/unit/math/dense_blas_internal_test.cpp index eaebe46..d3be2a4 100644 --- a/tests/unit/math/dense_blas_internal_test.cpp +++ b/tests/unit/math/dense_blas_internal_test.cpp @@ -11,6 +11,98 @@ namespace fesa::dense_blas_internal { namespace { +TEST(DenseBlasInternal, + OuterProductUsesRectangularStorageAndPreservesRejection) { + const std::array left{2.0, -3.0}; + const std::array right{4.0, -5.0, 0.25}; + std::array product{}; + ASSERT_TRUE(OuterProduct(left.data(), left.size(), right.data(), right.size(), + product.data())); + const std::array expected{8.0, -10.0, 0.5, -12.0, 15.0, -0.75}; + EXPECT_EQ(product, expected); + const auto accepted = product; + const std::array invalid{ + 1.0, std::numeric_limits::quiet_NaN()}; + EXPECT_FALSE(OuterProduct(invalid.data(), invalid.size(), right.data(), + right.size(), product.data())); + EXPECT_EQ(product, accepted); + EXPECT_TRUE(OuterProduct(nullptr, 0U, nullptr, 3U, nullptr)); + EXPECT_TRUE(OuterProduct(nullptr, 2U, nullptr, 0U, nullptr)); + EXPECT_THROW(OuterProduct(nullptr, 2U, right.data(), 3U, product.data()), + std::invalid_argument); + const auto too_large = + static_cast((std::numeric_limits::max)()) + 1U; + EXPECT_THROW(OuterProduct(nullptr, too_large, nullptr, 1U, nullptr), + std::length_error); + const auto too_many_elements = + std::numeric_limits::max() / sizeof(double) + 1U; + EXPECT_THROW(OuterProduct(nullptr, too_many_elements, nullptr, 1U, nullptr), + std::length_error); +} + +TEST(DenseBlasInternal, LocalDenseScopeRestoresTlsAfterProductsAndUnwind) { + const int previous = mkl_set_num_threads_local(0); + struct Restore { + int previous; + ~Restore() { static_cast(mkl_set_num_threads_local(previous)); } + } restore{previous}; + const int default_threads = mkl_get_max_threads(); + const std::array left{0.0, 2.0}; + const std::array right{3.0, 4.0}; + std::array product{}; + product.fill(std::numeric_limits::quiet_NaN()); + { + const ScopedLocalDenseThreads scope{}; + EXPECT_EQ(mkl_get_max_threads(), 1); + EXPECT_TRUE( + OuterProduct(left.data(), 2U, right.data(), 2U, product.data())); + EXPECT_DOUBLE_EQ(product[0U], 0.0); + EXPECT_DOUBLE_EQ(product[3U], 8.0); + } + EXPECT_EQ(mkl_get_max_threads(), default_threads); + EXPECT_EQ(mkl_set_num_threads_local(0), 0); + try { + const ScopedLocalDenseThreads scope{}; + OuterProduct(nullptr, 2U, right.data(), 2U, product.data()); + } catch (const std::invalid_argument&) { + } + EXPECT_EQ(mkl_set_num_threads_local(0), 0); + static_cast(mkl_set_num_threads_local(3)); + { + const ScopedLocalDenseThreads scope{}; + EXPECT_EQ(mkl_get_max_threads(), 3); + } + EXPECT_EQ(mkl_get_max_threads(), 3); +} + +TEST(DenseBlasInternal, OuterProductLeavesFiniteInputOverflowForOwnerGate) { + const std::array left{std::numeric_limits::max()}; + const std::array right{2.0}; + std::array product{}; + EXPECT_TRUE(OuterProduct(left.data(), 1U, right.data(), 1U, product.data())); + EXPECT_TRUE(std::isinf(product[0U])); +} + +TEST(DenseBlasInternal, TransposeRightProductUsesInnerAndRectangularStrides) { + const std::array left{1, 2, 3, 4, 5, 6}; + const std::array right{7, 8, 9, 10, 11, 12, + 13, 14, 15, 16, 17, 18}; + std::array product; + ASSERT_TRUE(MultiplyTransposeRight(left.data(), 2U, right.data(), 4U, 3U, + product.data())); + const std::array expected{50, 68, 86, 104, 122, 167, 212, 257}; + EXPECT_EQ(product, expected); + EXPECT_TRUE( + MultiplyTransposeRight(nullptr, 2U, nullptr, 4U, 0U, product.data())); + for (const auto value : product) EXPECT_DOUBLE_EQ(value, 0.0); + auto nonfinite = right; + nonfinite[8U] = std::numeric_limits::infinity(); + product.fill(31.0); + EXPECT_FALSE(MultiplyTransposeRight(left.data(), 2U, nonfinite.data(), 4U, 3U, + product.data())); + for (const auto value : product) EXPECT_DOUBLE_EQ(value, 31.0); +} + TEST(DenseBlasInternal, CDup003ConvertsZeroNormalAndOverflowLengths) { const auto zero = ToMklSize(0U); ASSERT_TRUE(zero.HasValue()); @@ -115,6 +207,30 @@ TEST(DenseBlasInternal, CholeskyAndInverseMatchKnownFactorsAndKeepFailedOutput) std::invalid_argument); } +TEST(DenseBlasInternal, + CholeskyFailureReasonPreservesCandidatesAndResetsOnReuse) { + Matrix matrix{2U, 2U}; + matrix(0U, 0U) = 4.0; + matrix(1U, 1U) = 9.0; + Matrix lower{0U, 0U}; + Matrix inverse{0U, 0U}; + auto reason = CholeskyFailure::kInverse; + ASSERT_TRUE(CholeskyLowerAndInverse(matrix, lower, inverse, &reason)); + EXPECT_EQ(reason, CholeskyFailure::kNone); + matrix(1U, 1U) = -1.0; + EXPECT_FALSE(CholeskyLowerAndInverse(matrix, lower, inverse, &reason)); + EXPECT_EQ(reason, CholeskyFailure::kPivot); + EXPECT_DOUBLE_EQ(lower(1U, 1U), 3.0); + EXPECT_DOUBLE_EQ(inverse(1U, 1U), 1.0 / 3.0); + matrix(1U, 1U) = 9.0; + EXPECT_TRUE(CholeskyLowerAndInverse(matrix, lower, inverse, &reason)); + EXPECT_EQ(reason, CholeskyFailure::kNone); + reason = CholeskyFailure::kInverse; + EXPECT_THROW(CholeskyLowerAndInverse(Matrix{2U, 1U}, lower, inverse, &reason), + std::invalid_argument); + EXPECT_EQ(reason, CholeskyFailure::kNone); +} + TEST(DenseBlasInternal, SingleEquationScopeRestoresTlsAndRespectsOverrides) { const int previous = mkl_set_num_threads_local(0); struct Restore {