diff --git a/docs/calculation-performance/mkl-coverage.md b/docs/calculation-performance/mkl-coverage.md new file mode 100644 index 0000000..3e0baf7 --- /dev/null +++ b/docs/calculation-performance/mkl-coverage.md @@ -0,0 +1,250 @@ +# MKL 적용 범위와 수명 근거 + +FEM-REQ-PERF-004/005/007/009의 T2 구현 기록이다. 공통 요구사항·물리·수치 +기준과 reference는 변경하지 않았다. T3/T4 인계 목록은 해당 owner가 완료한 +실제 diff와 검사를 통합한 뒤 갱신한다. 이 기록만으로 전체 production 적용이나 +세 모델의 공식 시간 목표 통과를 주장하지 않는다. + +## 공통 연산 경계 + +| 위치 | 연산 / 실제 private 호출 | 유지하는 의미 | +| --- | --- | --- | +| `math/matrix.cpp::Multiply` | dense gemv/gemm, `cblas_dgemv/dgemm` | 기존 shape/empty/checked MKL_INT 및 새 actual operand finite preflight | +| `math/vector.cpp` | Dot/Scale/Axpy와 copy, dense BLAS | 기존 공개 값·실패 의미 유지 | +| `math/dense_blas_internal.cpp::MultiplyTransposeLeft` | transpose storage copy 뒤 `Matrix::Multiply`; diagonal DSBMV(k0); LAPACK dpotrf/dtrtri | transpose-left shape 확인; ordinary product는 MKL | +| `math/sparse_matrix.cpp::Multiply` | `mkl_sparse_d_create_csr` → `mkl_sparse_d_mv` → destroy | canonical CSR와 structural zeros, general full matrix, zero-based index | +| `solvers/linear/mkl_pardiso_solver.cpp` | PARDISO SPD factorization/solve | 기존 경로와 scoped n=1 policy | +| `solvers/linear/mkl_pardiso_tangent_solver.cpp` | fresh phase 11/22, solve 33 | 새 분해 전 release; 실패 뒤 stale 사용 금지; 동일 성공 분해의 여러 RHS 허용 | +| `results/global_equilibrium_recovery.cpp` | full K*d, free/constrained masked RHS의 각 K*d | Kff*df/Kfc*dc scale를 별도로 유지; 원래 free rows만 norm 판정 | +| `analysis/linear_frequency_procedure.cpp` / `solid_frequency_procedure.cpp` | original/whitened pencil sparse action, dense Dot | original-pencil finite/independent residual/energy/Gram 검사 유지 | +| `solvers/eigen/mkl_feast_eigen_solver.cpp` | independent scaled full/mass/mirrored sparse actions; QR ddot/daxpy | FEAST actual 호출·count/solve budget·seed/Householder 절차 유지; backend residual로 독립 residual 대체하지 않음 | +| `analysis/eigen_reduction.cpp` | nodal/block congruence, mode restoration, nodal mass action | nodal analytic null split과 scaled original pencil 유지 | +| `analysis/thermal_transient_algebra.cpp` | BE RHS entry의 compact diagonal DSBMV, reduced/full/dissipation CSR actions | raw old constrained state·signed temperature·actual delta·independent interval companion 보존 | +| `elements/euler_beam_3d.cpp` | B*d, BᵀDB, TᵀKT, interpolationᵀq | 2 Gauss points·closed stiffness invariant·local force/section identity 유지 | + +### Sparse 단발 자원 + +입력 의존성은 현재 matrix의 전체 CSR shape/index/value와 현재 RHS이다. +`Validate()`가 canonical CSR와 finite coefficient를 확인한 뒤 actual `MKL_INT` +폭에 맞춰 rows/columns/nnz/각 index를 checked 변환한다. 공개 `size_t` storage는 +변경하지 않고 기존 private `ToMklSize`를 재사용한다. 이 설치의 ILP64를 public +API에 하드코딩하지 않는다. + +호출 stack이 변환 index 배열·복사 values와 handle을 소유한다. handle은 배열 뒤에 +생성되며 RAII destroy가 배열 소멸보다 먼저 실행된다. MKL이 CSR 배열을 참조한다는 +[Intel create CSR 계약](https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2026-0/mkl-sparse-create-csr.html)을 따른다. +직접 private values 복사는 const-cast 없이 const 동시 호출을 분리한다. +호출 종료가 invalidation이며 지속 handle/cache가 없다. 저장량은 단발 +`(rows+1+nnz)*sizeof(MKL_INT)+nnz*sizeof(double)` 및 output이다. +측정 근거 없는 symbolic/persistent cache는 채택하지 않았다. 아래 n=1 PARDISO의 제한된 TLS 정책만 실제 baseline/보조 근거로 채택했다. + +0 rows/0 columns/no entries는 검증과 index range 확인 뒤 zero vector를 반환한다. +nonempty matrix의 구조상 빈 row는 MKL 결과의 zero row다. Copy/move/assignment는 +공개 CSR 값 storage 수명만 갖고 backend pointer를 보관하지 않는다. Moved-from +불완전 CSR은 사용 전에 validation에서 거부되므로 indexing UB로 진행하지 않는다. + +backend failure는 scalar fallback이나 성공 candidate 없이 `runtime_error`로, +integer overflow는 `length_error`로 전달한다. vendor enum/원시 오류 번호는 공개하지 +않는다. `GlobalEquilibriumRecovery`와 frequency/eigen/thermal consumer는 해당 +새 실패를 `Status/Result`로 변환하며 이전 state/history/final 결과를 commit하지 않는다. +실제 allocation/backend fault를 강제한 native 시험은 실행하지 않았다. 해당 +실패 branch의 fallback 부재·exception 변환·RAII 해제는 정적 검토 범위이며, +일반 실패/반복·copy/move·동시 const 호출과 bad RHS 복구는 실제 검사 범위다. + +### Finite의 owner + +CSR coefficient와 모든 RHS는 MKL 호출 전 math에서 finite 검증한다. 비유한 RHS는 +기존 값 반환/exception 형식의 `invalid_argument`로 거부하며, public signature는 +변경하지 않는다. Finite-input product overflow는 consumer의 기존 output finite +검증으로 실패한다. Empty/zero shape도 입력 finite/shape/range 검증을 유지한다. + +Dense gemv/gemm도양쪽actualownedstorage/VectorRHS를중앙Matrix::Multiply에서finite-first로검사한다. PrivateTransposeLeft는이중앙dispatch를사용하며 nonfinite입력에서invalid_argument를전달한다. BuildEigenReduction/linearfrequency/FEAST owner는해당typedfailure를기존Result taxonomy로변환한다. MakeLift actual basis/scaled operands와FEAST normalizedphi/mass_modes producer도dispatch전에검사한다. + +Production caller finite 입력은 별도로 확인한다. Recovery는 기존 +`IsFiniteVector`→Multiply→finite output 순서를 유지한다. Linear frequency의 physical +및 whitened x, solid frequency의 raw/normalized mode는 finite-first이며, FEAST가 +반환한 벡터의 finite gate도 기존대로 유지한다. Thermal은 기존 Validate·solver +candidate finite와 interval enclosures를 사용한다. Thermal RHS의 각 entry는 +`(c/delta)*(old-g)-k*g`라는 원래 association을 먼저 계산한 뒤 원래 row 순서로 +누적한다. 두 행곱을 별도로 합산한 뒤 빼는 후보는 stable b=4를 잃어 기존 +Cancellation 검사를 실패했으므로 폐기했다. 현재 private diagonal MKL은 O(nnz) +compact band k=0 storage로 entry별 곱만 계산한다. Exact companion과 construction +error·gamma·underflow certification·free/reaction 판정은 그대로다. +비유한 product를 정상값으로 clamp하거나 tolerance를 완화하지 않는다. + +### PARDISO 수명 + +현재 구현 그대로 `release → CSR/shape/width 확인 → private arrays → phase11 → phase22` +를 수행한다. 값/pattern이 달라지면 `FactorizeFresh`를 다시 호출한다. 새 singular +분해 실패 뒤 `SolveCorrection`은 not-factorized이며 caller correction은 보존된다. +잘못된 RHS의 실패는 output 보존과 유효한 기존 matrix 분해의 수명을 구별한다. +여러 RHS는 같은 성공 분해의 phase33을 사용할 수 있다. 새로운 symbolic 재사용이나 +다중 RHS framework를 도입하지 않았다. + +### One-equation PARDISO의 scoped 자원 정책 + +[T1 실제 VTune](verification.md)의 OMP worker/spin 비용과 아래 toy n=1 lifecycle +근거로, SPD/Tangent의 `CallPardiso`에서 **현재 retained equation_count=1**인 경우만 +호출 기본 thread 정책을1로 제한한다. Phase11/22/33/−1 각각의 호출 범위다. +`mkl_set_num_threads_local(1)`의 이전 TLS override가0인 경우만 scope1을 유지하고 +끝에0을 복원한다. 이전값이비영이면 즉시 그값을복원하여 caller override대로 +backend를실행한다. n=0/2이상·전역 MKL/OpenMP 요청·dispatch·iparm·pivot/scaling/ +regularization·acceptance는변경하지않는다. 넓은size threshold는추정하지않았다. + +Release의oldfactor는old retained count, newfactor는newcount를각call에서사용하므로 +작은oldmatrix release 뒤 큰newmatrix가 single-thread로 오염되지 않는다. +현재원시기록의 외부 MKL/OMP 요청은8/dynamicFALSE로baseline과같고 실제 source +정책은default1×1 calls에서달라진다. 이 차이를환경동일이라는표현으로숨기지않는다. +Phase 성공/오류/exception·잘못된 RHS 및 concurrent worker의 TLS 복원과 candidate +보존을실제검사한다. Native process crash는C++ RAII unwind를보장하지않으며, 해당 +진단의scope설정은종료한process의TLS여서다음실행/전역환경에남지않는다. + +### Baseline PARDISO heap failure와 제한 + +기존staticpartialfixture는full24×24/free1×1이다. `backend-isolation`의SpMV-only +1000회create/mv/destroy는native0이며, PARDISO-only100회예정실행은native +`−1073740940 (0xc0000374)`로중단됐다. 별도반복staticfixture도iteration3에서같은 +heap-corruption을재현했다. 이를새SpMV수명회귀나단순testPASS로분류하지않는다. + +같은freshauxTU를T1의기존Release library에`/O2 /Ob2 /MD /DNDEBUG /Zi`로link해 +실제baseline을대조했다. Baseline default8은61회완료뒤62에서같은nativecrash, +scoped1은100회native0와각scope끝effective8복원을기록했다. CurrentDebug +`/Od /Ob0 /MDd /RTC1 /Zi` recordedprobe의8/1은각100회완료했으나이전Debug +probe의2/3회crash를지우지않는다. 실패가간헐적이므로100회통과를vendor문제 +해결확정으로쓰지않는다. Debug/Release차이도그대로기록한다. + +Ordered full/free CSR offsets/columns/rawdoublebits와RHS inputSHA는양쪽 +`CC7E72000783E9694BAF76627D95AB7B420A31198B210B5A60C7ACE93E5E2258`로같다. +ILP64·mtype2·actual pardisoinit defaultcontrols(26 checker/34 zerobased 포함), 동일 +MKL/runtimeDLL hash는[recorded-probes metadata](../../.harness/performance/t2/recorded-probes-metadata.json)와 +[native 기록](../../.harness/performance/t2/recorded-probes-native.json)에있다. +기존PARDISO/backend영역실패를실제baseline에서재현했으며vendor내부원인은미확정이다. +[Intel release notes](https://www.intel.com/content/www/us/en/developer/articles/release-notes/onemkl/2026.html)의 +non-AVX2 PARDISO sporadic pivoting issue는조사근거이고현재원인확정이아니다. +설치·dispatch·전역registry/pageheap변경은하지않았다. + +같은toy n=1 보조비용실행에서default8은11회뒤crash/0.24s, +scope1은100회완료/0.15s였다. 이는동일입력준비/로그를포함한auxiliary이며 +실패/표본수차이가있어공식whole-model speedup이나시간목표판정에사용하지않는다. +[Tiny native 비용기록](../../.harness/performance/t2/tiny-pardiso-cost-native.json)에실제범위를기록했다. + +## 특수 정밀도 개선과 예외 + +T1의 실제 Hex8 VTune stack은 `TwoSum → Add → Multiply → hex8_high_low::operator*` +비용 및 frexp/ldexp 분류 비용을 확인했다. T2는 arithmetic association과 underflow +정책을 바꾸지 않고 두 helper만 개선했다. + +- `high_low_arithmetic.h::IsNormalized`: finite와 canonical exact-zero guard 뒤 + `RN(high+low)==high`. 이 경우 기존 TwoSum의 bridge=0/restored=high/error=low이며, + 반대 방향도 기존 predicate의 high equality가 같은 조건이다. overflow/nonfinite와 + noncanonical signed zero를 먼저 거부한다. +- `LowestBitExponent`: IEEE binary64를 memcpy로 읽어 정상은 + `(encoded exponent−1075)+trailing(implicit53bits)`, subnormal은 + `−1074+trailing(fraction)`으로 계산한다. 부호를 제외하고 actual least bit를 보존한다. + TwoProduct의 least-bit-sum < −1074 거부와 FMA residual·TwoSum 정규화는 그대로다. +- zero/low-zero shortcut, nonzero multiplication association 변경은 채택하지 않았다. + +[독립 수치 검토](../../.harness/performance/t2/arithmetic-numerical-review.md)는 +위 두 helper proof/source 범위 PASS이며 실제 검사/속도와 전체 MKL diff 검토를 +대체하지 않는다. Frozen original oracle는 pinned `b4fc1bb`의 원본이며 include guard와 +namespace만 바꾼 test 전용 `high_low_arithmetic_original.h`이다. 실제 모든 정상 +exponent×53 trailing-bit×양 부호, subnormal boundaries, signed zero/ties/overflow/ +invalid/raw random 입력의 bit 대조 및 기존 analytical checks를 실행했다. + +| 파일 / 함수 | 연산·크기·정밀도 | 예외 이유 / 검사·owner | +| --- | --- | --- | +| `math/high_low_arithmetic.h` 전체, `elements/hex8_high_low.cpp::EvaluateHighLow` | normalized binary64 high+low, 3×3 tensor 및 24×24 pair tangent | ordinary Float64 BLAS로 low/tail rejection을 잃을 수 없음; frozen oracle/Hex8HighLow·T3 | +| `elements/hex8_geometry.cpp`의 interval/relative coordinate 인증 | 3×3 inverse/Jacobian 및 outward interval | whole-cell certificate·직접 geometry predicate 보존; Hex8 geometry 검사·T3 | +| `elements/hex8_element.cpp` gradient/inverse·F/S/P/FFT/stress의 3×3 tensor | fixed3-coordinate contraction | 3-coordinate geometry/constitutive 의미; 24×24 ordinary tangent는 이 예외로 일괄 분류하지 않음·T3 | +| `elements/finite_rotation.cpp` Rodrigues/skew/invariant | fixed3×3 axis coordinate | 좌표 rotation/invariant; finite rotation 검사 | +| `analysis/eigen_reduction.cpp::MakeLift` analytic director/frame basis와 null-coordinate checks | fixed6 nodal coordinates, 최대6 positive basis | 고정 director/frame의 analytic 좌표 선택만 예외. 일반 mass whitening/Cholesky/inverse/congruence/reconstruction은 MKL 전환·EigenReduction | +| `elements/thermal_shell.cpp` shape/metric contraction | fixed4 shape·fixed3 gradient·2D metric, separable integration coefficient | local coordinate/shape quadrature; global capacity/conduction actions는 MKL·thermal local oracle | +| `analysis/thermal_transient_algebra.cpp`의 `I`/Add/Mul/Sub/Div/Norm/Error/Roundoff/Tiny | outward interval·construction/underflow enclosure | floating product와 분리된 정확대수 인증, gamma/epsilon/limit 변경없음·TransientHeatProcedure/ThermalPhysics | +| `analysis/hht_reduced_operators.cpp` Projection/Embed/Sum, BasisValue, reduced blocks | compensated sums 및 outward basis radius | low/error compensation 유지, ordinary Hessian/product는 T4 판별/전환 | +| `analysis/hht_dynamic_algebra.cpp` AddCompensated·constraint radius 및 3-vector director dot | compensated error·fixed coordinate | 승인 N3/N7·scale 인증; 일반 CSR contraction과 구별·T4 | +| `elements/nonlinear_mitc4_shell.cpp` / `nonlinear_mitc4_inertia.cpp` Jet 산술 | analytic derivative jets, fixed3 rotation | AD derivative association/coordinate 의미; ordinary24×24 Hessian/dense congruence blanket 예외 금지·T4 | +| `io/hdf5/nonlinear_hdf5_results_writer.cpp` spatial moment transform | fixed3×3 rotation coordinate | 공개 generalized/spatial quantity 변환 보존·T4 | + +## T3/T4 인계 + +전수 검색의 raw 위치 목록은 [ordinary-product-audit](../../.harness/performance/t2/ordinary-product-audit.txt)다. +검색은 +=*만으로 완료하지 않고 Multiply·중첩 루프와 각 함수의 precision/owner를 +함께 확인했다. 이 목록은 후속 diff에서 갱신해야 한다. + +- T3: Hex8 pair tangent/relative geometry/mean-volume certification은 pair/interval + 예외로 유지한다. Float64 24×24 response/tangent의 일반 contractions와 procedure + 직접 products는 ordinary MKL 적용 여부를 확인한다. +- T4: `mitc4_shell.cpp::` transformation congruence(기존241 주변), strainᵀCstrain + stiffness(기존635 주변), `nonlinear_mitc4_inertia.cpp` complement Gram/projection + (기존449/456), `hht_reduced_operators.cpp` bilinear Hessian action(기존620), + reduced block/gradient products(기존970 이후), `hht_dynamic_algebra.cpp` 일반 + stiffness projection(기존442)은 MKL 또는 구체적인 compensated 예외 근거가 필요하다. + `Jet*`의 derivative array 크기가24라는 이유만으로 모든24×24 Float64 products를 + 특수 정밀도 예외로 분류하지 않는다. +- 공통 recovery/eigen/thermal/beam은 T2 확장 owner로 직접 변경했다. Shared + assembly tuple reduction/source identity와 task별 CMake 등록은 root owner이며 변경하지 않았다. + +## 실제 검증과 남은 작업 + +Parent artifact `.harness/performance/t2/`는 worktree cleanup 뒤에도 남는다. +각 native JSON은 argv/cwd/실제 exit와 source 또는 snapshot 구분을 보존한다. +아래 링크는 통합한 저장소의 feature 문서 위치를 기준으로 한다. + +| 실행 범위 / source | 결과 | 실제 기록 | +| --- | --- | --- | +| 최초 Debug build | include 누락 exit2; 수정 뒤 build0 | `build-debug-initial-failed.txt`, `build-native-initial-failed.json` | +| 초기 개발 snapshot | core43/43·consumer71/71 PASS; whole-source hash 없는 precommit 자료 | `core-native.json`, `consumer-native.json` | +| Extended66 개발 snapshot | 64 PASS·1 thermal FAIL·1 large-mesh interrupted; whole native−1 | `extended-native.json`, `extended-checks.txt`, `extended-interruption.json`, `extended-cleanup.json` | +| 59c8e93 후보100개 | 99 PASS·열 RHS 상쇄1 FAIL; native8 | `final-native.json`, `final-checks.txt` | +| Thermal/finite/Gram/whitening 수정 snapshot82개 | 82/82 PASS; native0 | `corrected-native.json`, `corrected-checks.txt` | +| 후속 eigen20개 | 19 PASS·기존 PARDISO heap fixture1 native crash; native8 | `null-products-native.json`, `null-products-checks.txt` | +| 770fa2a scoped-policy 후보 | build0·99/99 PASS; native0 | `policy-build-native.json`, `policy-native.json` | +| 9c4e685 중앙 dense finite 후보 | build0·120/120 PASS; native0 | [build](../../.harness/performance/t2/finite-build-native.json), [native](../../.harness/performance/t2/finite-native.json), [log](../../.harness/performance/t2/finite-checks.txt) | +| 6fc7011 tests-only delta | build0·1↔2 retained-size 전환2/2 PASS; native0 | [build](../../.harness/performance/t2/transition-build-native.json), [native](../../.harness/performance/t2/transition-native.json), [log](../../.harness/performance/t2/transition-checks.txt) | + +`MklPardisoSolver.SingleEquationScopeTransitionsBetweenRetainedSizes`와 +`TangentSolver.SingleEquationScopeTransitionsBetweenRetainedSizes`는 각각 +1→2→1 factorize/solve와 old/new retained count별 TLS 복원을 직접 확인한다. +동시 default TLS의 성공→다른값→bad RHS→새 분해 실패→stale solve 거부는 +`SingleEquationPolicyPreservesLifecycleAndConcurrentTls` 및 +`SingleEquationScopeRestoresTlsAcrossFreshFailureAndReuse`에서 확인한다. +Caller의 비영 TLS override·exception unwind와 큰 크기 정책은 +`DenseBlasInternal.SingleEquationScopeRestoresTlsAndRespectsOverrides`의 직접 관찰이다. + +SpMV destroy/PARDISO release-error를 live 주입한 검사는 미실행이다. 유효 handle에서 +이를 안전하게 발생시킬 근거가 없고 production mutable fault hook은 추가하지 않았다. +정상 release·분해 오류 후 release·반복 및 크기 전환은 실제 검사이며, 원시 release +오류 분기의 candidate 실패 전달과 자원 수명은 독립 정적 검토 범위다. + +9c4e685 production source의 [독립 conformance 재검토](../../.harness/performance/t2/source-conformance-review.md)는 PASS다. 이후 6fc7011은 +검토자가 요청한 transition tests-only delta이며 production source는 같다. +High+low 두 helper의 독립 proof는 앞 절의 실제 보고서와 연결했다. +검토 PASS와 actual execution PASS는 서로 대체하지 않는다. + +최초 실패·상쇄 회귀·간헐적 heap 실패·중단 자료를 삭제하거나 현재 PASS로 바꾸지 않았다. +전체 MSVC Debug/CTest, 원본3+Hex8 필수5 비교, VTune 최종 raw 및 공식 +no-warmup/one-measurement 평가는 T5에서 미실행 항목을 완료해야 한다. +Task source를 fresh build했으며 다른 source의 object를 수동 복사하지 않았다. +Baseline-linked auxiliary는 root가 승인한 보존된 T1 Release library를 직접 link하고 +그 flags/runtime/input을 별도로 기록했으므로 현재 Debug evidence로 대체하지 않는다. + +High+low standalone 보조 microbenchmark는 같은 compiler/한 exe에서 +original/optimized 각400000 products를3회 대조했다. Original은 +0.129315/0.129951/0.129560s, optimized는0.069739/0.070177/0.068700s이며 +모든 checksum이 같고 native0이다. +[보조기록](../../.harness/performance/t2/microbench.jsonl)과 +[native](../../.harness/performance/t2/microbench-native.json)는 helper 비용 근거이며 +3모델의 공식 wall 평가를 대체하지 않는다. + +CSR consumer 전수 보강 검색은 +[csr-consumer-audit](../../.harness/performance/t2/csr-consumer-audit.txt)에 있다. +`assembly/load_assembler.cpp::EffectiveFreeRhs`의 Kfc*dc는 root owner가 T2 통합 뒤 +entrywise MKL diagonal products와 기존 개별 product/ordered row-sum finite gate로 +연결한다. 이 공유 파일은 T2가 수정하지 않았다. 최종 production 전수 완료는 +이 인계와 T3/T4 일반 연산 적용까지 확인해야 한다. + +최종 source/test hash·production 검토 이후 source 변경 없음·874개 검사 discovery native0은 +[final-source-manifest](../../.harness/performance/t2/final-source-manifest.json)에 기록했다. +[final CMakeCache](../../.harness/performance/t2/final-CMakeCache.txt)는 actual dependency와 Debug flags를 보존한다. +Whole changed-file clang-format dry-run은 기존 및 후보 formatting 위반으로 exit1이었다. +[실제 기록](../../.harness/performance/t2/clang-format-native.json)을 보존하며 formatter PASS로 보고하지 않는다. diff --git a/include/fesa/math/matrix.h b/include/fesa/math/matrix.h index eb700a8..af51673 100644 --- a/include/fesa/math/matrix.h +++ b/include/fesa/math/matrix.h @@ -41,11 +41,11 @@ class Matrix { const double& operator()(std::size_t row, std::size_t column) const; /// @brief Multiplies this row-major matrix by a dense vector. - /// @throws std::invalid_argument if the dimensions are incompatible. + /// @throws std::invalid_argument for incompatible shapes or nonfinite input. Vector Multiply(const Vector& rhs) const; /// @brief Multiplies this row-major matrix by another dense matrix. - /// @throws std::invalid_argument if the dimensions are incompatible. + /// @throws std::invalid_argument for incompatible shapes or nonfinite input. Matrix Multiply(const Matrix& rhs) const; private: diff --git a/include/fesa/math/sparse_matrix.h b/include/fesa/math/sparse_matrix.h index dedf73e..1aabf21 100644 --- a/include/fesa/math/sparse_matrix.h +++ b/include/fesa/math/sparse_matrix.h @@ -46,8 +46,12 @@ class SparseMatrix { /// @brief Returns CSR values including preserved structural zeros. const std::vector& Values() const noexcept; - /// @brief Multiplies this matrix by a dense vector in stable CSR order. - /// @throws std::invalid_argument if the dimensions are incompatible. + /// @brief Multiplies canonical CSR by a dense vector through private MKL. + /// @throws std::invalid_argument if dimensions, CSR, or RHS are invalid. + /// @throws std::length_error if a size exceeds the backend integer range. + /// @throws std::runtime_error if the sparse backend fails. + /// @note Calls retain no handle; concurrent const calls own their resources. + /// Finite-input overflow is retained for consumer output finite validation. Vector Multiply(const Vector& rhs) const; /// @brief Validates shape, indices, ordering, and finite CSR values. diff --git a/src/fesa/analysis/eigen_reduction.cpp b/src/fesa/analysis/eigen_reduction.cpp index c02514e..2b172e4 100644 --- a/src/fesa/analysis/eigen_reduction.cpp +++ b/src/fesa/analysis/eigen_reduction.cpp @@ -5,11 +5,13 @@ #include #include #include +#include #include #include #include #include "fesa/fem/dof_manager.h" +#include "math/dense_blas_internal.h" namespace fesa::eigen_internal { namespace { @@ -117,18 +119,31 @@ Result MakeLift(const Matrix& mass, if (positive == 0U) return Result::Success(std::move(result)); const Matrix scaled = ScaledMass(mass, length); - Matrix mp(positive, positive); + // Raw frame tangents can enter this internal boundary independently of + // Domain validation. Check the actual operands before any dense dispatch. + for (std::size_t row = 0; row < 6U; ++row) { + for (std::size_t j = 0; j < basis.Columns(); ++j) { + if (!std::isfinite(basis(row, j))) + return Result::Failure( + Failure("invalid-eigen-mass-space", node, + "Invalid positive mass basis.")); + } + for (std::size_t j = 0; j < 6U; ++j) { + if (!std::isfinite(scaled(row, j))) + return Result::Failure( + Failure("invalid-eigen-mass-space", node, + "Scaled mass is nonfinite.")); + } + } + Matrix mp = dense_blas_internal::MultiplyTransposeLeft(basis, scaled) + .Multiply(basis); double scale = 0.0; double basis_error = 0.0; + const auto basis_gram = + dense_blas_internal::MultiplyTransposeLeft(basis, basis); for (std::size_t i = 0; i < positive; ++i) { for (std::size_t j = 0; j < positive; ++j) { - double gram = 0.0; - for (std::size_t a = 0; a < 6U; ++a) { - gram += basis(a, i) * basis(a, j); - for (std::size_t b = 0; b < 6U; ++b) { - mp(i, j) += basis(a, i) * scaled(a, b) * basis(b, j); - } - } + const double gram = basis_gram(i, j); basis_error = std::hypot(basis_error, gram - (i == j ? 1.0 : 0.0)); scale = (std::max)(scale, std::abs(mp(i, j))); } @@ -138,32 +153,34 @@ Result MakeLift(const Matrix& mass, return Result::Failure(Failure("invalid-eigen-mass-space", node, "Invalid positive mass basis.")); } + Matrix normalized_mp(positive, positive); + for (std::size_t i = 0; i < positive; ++i) + for (std::size_t j = 0; j < positive; ++j) { + normalized_mp(i, j) = mp(i, j) / scale; + if (!std::isfinite(normalized_mp(i, j))) + return Result::Failure( + Failure("eigen-mass-conditioning", node, + "A positive-space mass pivot is not resolvable.")); + } Matrix lower(positive, positive); - for (std::size_t i = 0; i < positive; ++i) { - for (std::size_t j = 0; j <= i; ++j) { - double value = mp(i, j) / scale; - for (std::size_t k = 0; k < j; ++k) value -= lower(i, k) * lower(j, k); - if (i == j) { - if (!(value > 0.0) || !std::isfinite(value)) { - return Result::Failure( - Failure("eigen-mass-conditioning", node, - "A positive-space mass pivot is not resolvable.")); - } - lower(i, j) = std::sqrt(value); - } else { - lower(i, j) = value / lower(j, j); - } - } - } Matrix inverse_lower(positive, positive); - for (std::size_t j = 0; j < positive; ++j) { - for (std::size_t i = 0; i < positive; ++i) { - double value = i == j ? 1.0 : 0.0; - for (std::size_t k = 0; k < i; ++k) - value -= lower(i, k) * inverse_lower(k, j); - inverse_lower(i, j) = value / lower(i, i); - } + if (!dense_blas_internal::CholeskyLowerAndInverse(normalized_mp, lower, + inverse_lower)) { + return Result::Failure( + Failure("eigen-mass-conditioning", node, + "A positive-space mass pivot is not resolvable.")); } + Matrix transposed_lower(positive, positive); + Matrix scaled_inverse_transpose(positive, positive); + const double root_scale = std::sqrt(scale); + for (std::size_t i = 0; i < positive; ++i) + for (std::size_t j = 0; j < positive; ++j) { + transposed_lower(i, j) = lower(j, i); + scaled_inverse_transpose(i, j) = inverse_lower(j, i) / root_scale; + } + const auto inverse_mass = dense_blas_internal::MultiplyTransposeLeft( + inverse_lower, inverse_lower); + const auto rebuilt_mass = lower.Multiply(transposed_lower); double norm_mp = 0.0; double norm_inverse = 0.0; double reconstruction = 0.0; @@ -173,12 +190,8 @@ Result MakeLift(const Matrix& mass, double row_mp = 0.0; double row_inverse = 0.0; for (std::size_t j = 0; j < positive; ++j) { - double inverse = 0.0; - double rebuilt = 0.0; - for (std::size_t k = 0; k < positive; ++k) { - inverse += inverse_lower(k, i) * inverse_lower(k, j); - rebuilt += lower(i, k) * lower(j, k); - } + const double inverse = inverse_mass(i, j); + const double rebuilt = rebuilt_mass(i, j); row_inverse += std::abs(inverse); row_mp += std::abs(mp(i, j) / scale); reconstruction = std::hypot(reconstruction, rebuilt - mp(i, j) / scale); @@ -196,13 +209,9 @@ Result MakeLift(const Matrix& mass, "eigen-mass-conditioning", node, "Nodal mass fails rcond_inf >= 2^-32 or Cholesky reconstruction.")); } - const double root_scale = std::sqrt(scale); + result.scaled_lift = basis.Multiply(scaled_inverse_transpose); for (std::size_t c = 0; c < 6U; ++c) { for (std::size_t j = 0; j < positive; ++j) { - for (std::size_t k = 0; k < positive; ++k) { - result.scaled_lift(c, j) += - basis(c, k) * inverse_lower(j, k) / root_scale; - } if (!std::isfinite(result.scaled_lift(c, j))) { return Result::Failure(Failure( "eigen-mass-conditioning", node, "Nodal whitening is nonfinite.")); @@ -218,7 +227,7 @@ Result BuildEigenReduction( const SparseMatrix& stiffness, const std::vector& mass, const std::vector& frames, const std::vector& mask, const std::vector& drilling, - const double length) { + const double length) try { const std::size_t count = mass.size(); const Status valid = stiffness.Validate(); if (!valid.IsOk()) return Result::Failure(valid); @@ -281,19 +290,23 @@ Result BuildEigenReduction( free_count += lift.Value().record.free_count; null_count += lift.Value().record.null_count; const auto& z = lift.Value().null_vector; + Vector null_coordinate{6U}; + Matrix null_column{6U, 1U}; + for (std::size_t a = 0; a < 6U; ++a) { + null_coordinate[a] = z[a]; + null_column(a, 0U) = z[a]; + } + const auto mass_null_action = scaled.Multiply(null_coordinate); + const auto positive_null = dense_blas_internal::MultiplyTransposeLeft( + lift.Value().positive_basis, null_column); for (std::size_t a = 0; a < 6U; ++a) { null_norm = std::hypot(null_norm, z[a]); expected_kz_norm = std::hypot(expected_kz_norm, drilling[node] * z[a]); if (mask[6U * node + a] != 0U) continue; - double action = 0.0; - for (std::size_t b = 0; b < 6U; ++b) action += scaled(a, b) * z[b]; - mass_null_error = std::hypot(mass_null_error, action); + mass_null_error = std::hypot(mass_null_error, mass_null_action[a]); } for (std::size_t j = 0; j < lift.Value().positive_basis.Columns(); ++j) { - double dot = 0.0; - for (std::size_t a = 0; a < 6U; ++a) - dot += lift.Value().positive_basis(a, j) * z[a]; - pz_norm = std::hypot(pz_norm, dot); + pz_norm = std::hypot(pz_norm, positive_null(j, 0U)); } lifts.push_back(std::move(lift.Value())); } @@ -319,10 +332,6 @@ Result BuildEigenReduction( for (std::size_t a = 0; a < 6U; ++a) { const auto row = node * 6U + a; if (mask[row] != 0U) continue; - std::map null_action; - if (lifts[node].record.null_count != 0U) { - null_action[node] = -drilling[node] * lifts[node].null_vector[a]; - } for (std::size_t pos = offsets[row]; pos < offsets[row + 1U]; ++pos) { const auto col = columns[pos]; if (mask[col] != 0U) continue; @@ -335,25 +344,40 @@ Result BuildEigenReduction( "Scaled stiffness is nonfinite.")); } stiffness_norm = std::hypot(stiffness_norm, value); - null_action[other] += value * lifts[other].null_vector[b]; auto block = blocks.try_emplace(other, 6U, 6U).first; block->second(a, b) = value; } - for (const auto& item : null_action) - kz_error = std::hypot(kz_error, item.second); } const auto& left = lifts[node]; + std::map null_actions; + for (const auto& entry : blocks) { + Vector null_coordinate{6U}; + for (std::size_t b = 0; b < 6U; ++b) + null_coordinate[b] = lifts[entry.first].null_vector[b]; + null_actions.emplace(entry.first, entry.second.Multiply(null_coordinate)); + } + if (left.record.null_count != 0U) { + auto& local = null_actions.try_emplace(node, 6U).first->second; + for (std::size_t a = 0; a < 6U; ++a) + local[a] -= drilling[node] * left.null_vector[a]; + } + for (std::size_t a = 0; a < 6U; ++a) { + if (mask[node * 6U + a] != 0U) continue; + for (const auto& action : null_actions) + kz_error = std::hypot(kz_error, action.second[a]); + } + std::map projected_blocks; + for (const auto& entry : blocks) { + projected_blocks.emplace( + entry.first, dense_blas_internal::MultiplyTransposeLeft( + left.scaled_lift, entry.second) + .Multiply(lifts[entry.first].scaled_lift)); + } for (std::size_t i = 0; i < left.scaled_lift.Columns(); ++i) { - for (const auto& entry : blocks) { + for (const auto& entry : projected_blocks) { const auto& right = lifts[entry.first]; for (std::size_t j = 0; j < right.scaled_lift.Columns(); ++j) { - double value = 0.0; - for (std::size_t a = 0; a < 6U; ++a) { - for (std::size_t b = 0; b < 6U; ++b) { - value += left.scaled_lift(a, i) * entry.second(a, b) * - right.scaled_lift(b, j); - } - } + const double value = entry.second(i, j); pattern.column_indices.push_back(right.offset + j); contributions.push_back({left.offset + i, right.offset + j, value, node, pattern.column_indices.size() - 1U}); @@ -378,6 +402,10 @@ Result BuildEigenReduction( return Result::Success( {std::move(reduced.Value()), std::move(lifts), free_count, null_count, length, eta_mz, eta_kz, pz_norm}); +} catch (const std::invalid_argument&) { + return Result::Failure( + Failure("invalid-eigen-mass-space", 0U, + "Eigen product operands must be finite and shape compatible.")); } Vector RestoreScaledMode(const EigenReduction& reduction, const Matrix& vectors, @@ -385,11 +413,12 @@ Vector RestoreScaledMode(const EigenReduction& reduction, const Matrix& vectors, Vector result(reduction.nodes.size() * 6U); for (std::size_t node = 0; node < reduction.nodes.size(); ++node) { const auto& lift = reduction.nodes[node]; + Vector reduced{lift.scaled_lift.Columns()}; + for (std::size_t j = 0; j < reduced.Size(); ++j) + reduced[j] = vectors(lift.offset + j, column); + const auto restored = lift.scaled_lift.Multiply(reduced); for (std::size_t a = 0; a < 6U; ++a) { - for (std::size_t j = 0; j < lift.scaled_lift.Columns(); ++j) { - result[6U * node + a] += - lift.scaled_lift(a, j) * vectors(lift.offset + j, column); - } + result[6U * node + a] = restored[a]; } } return result; @@ -399,10 +428,11 @@ Vector MultiplyNodalMass(const std::vector& mass, const Vector& vector) { Vector result(vector.Size()); for (std::size_t node = 0; node < mass.size(); ++node) { + Vector local{6U}; + for (std::size_t b = 0; b < 6U; ++b) local[b] = vector[6U * node + b]; + const auto product = mass[node].Multiply(local); for (std::size_t a = 0; a < 6U; ++a) { - for (std::size_t b = 0; b < 6U; ++b) { - result[6U * node + a] += mass[node](a, b) * vector[6U * node + b]; - } + result[6U * node + a] = product[a]; } } return result; diff --git a/src/fesa/analysis/linear_frequency_procedure.cpp b/src/fesa/analysis/linear_frequency_procedure.cpp index 1064c4e..fe3110c 100644 --- a/src/fesa/analysis/linear_frequency_procedure.cpp +++ b/src/fesa/analysis/linear_frequency_procedure.cpp @@ -32,11 +32,7 @@ Status Failure(const std::string& code, const std::string& identity, } double Dot(const Vector& left, const Vector& right) { - double result = 0.0; - for (std::size_t index = 0U; index < left.Size(); ++index) { - result += left[index] * right[index]; - } - return result; + return left.Dot(right); } double Norm(const Vector& vector) { @@ -279,6 +275,10 @@ Result ExecuteCandidate(const Domain& domain) { std::size_t pivot = 0U; for (std::size_t dof = 0; dof < y.Size(); ++dof) { physical[dof] = D(dof, length) * y[dof]; + if (!std::isfinite(physical[dof])) + return Result::Failure( + MetricFailure("initial_modal_mass", physical[dof], + state.backend.column_ids[column])); if (std::abs(y[dof]) > std::abs(y[pivot])) pivot = dof; } Vector mass_times = eigen_internal::MultiplyNodalMass(nodal_mass, physical); @@ -318,8 +318,13 @@ Result ExecuteCandidate(const Domain& domain) { (pencil_scale * y_norm) * y_norm); const double rayleigh = energy / modal_mass; Vector x(a.Rows()); - for (std::size_t row = 0; row < x.Size(); ++row) + for (std::size_t row = 0; row < x.Size(); ++row) { x[row] = state.backend.eigenvectors(row, column) * normalization; + if (!std::isfinite(x[row])) + return Result::Failure( + Failure("invalid-eigenvalue-result", std::to_string(column), + "Recovered mode is nonfinite or violates constraints.")); + } Vector ax = a.Multiply(x); Vector whitened_residual = ax; whitened_residual.Axpy(-lambda, x); @@ -412,6 +417,10 @@ Result LinearFrequencyProcedureExecutor::Execute( std::move(diagnostics))); } return result; + } catch (const std::invalid_argument&) { + return Result::Failure( + Failure("invalid-eigenvalue-result", "modal-workspace", + "Modal product inputs must be finite and shape compatible.")); } catch (const std::bad_alloc&) { return Result::Failure( Failure("eigen-allocation-failed", "modal-workspace", @@ -421,6 +430,10 @@ Result LinearFrequencyProcedureExecutor::Execute( return Result::Failure( Failure("eigen-index-overflow", "modal-workspace", "Requested workspace exceeds container capacity.")); + } catch (const std::runtime_error&) { + return Result::Failure( + Failure("eigen-multiplication-failed", "modal-workspace", + "Modal sparse multiplication failed.")); } } diff --git a/src/fesa/analysis/solid_frequency_procedure.cpp b/src/fesa/analysis/solid_frequency_procedure.cpp index 0c92757..f3a33dc 100644 --- a/src/fesa/analysis/solid_frequency_procedure.cpp +++ b/src/fesa/analysis/solid_frequency_procedure.cpp @@ -248,8 +248,13 @@ Result ExecuteCandidate(const Domain& domain) { mass_phi.reserve(count); for (std::size_t c = 0; c < count; ++c) { Vector vector(k.Rows()); - for (std::size_t row = 0; row < k.Rows(); ++row) + for (std::size_t row = 0; row < k.Rows(); ++row) { vector[row] = state.backend.eigenvectors(row, c); + if (!std::isfinite(vector[row])) + return Result::Failure( + Failure("solid-eigenpair-validation-failed", + "Modal mass is not finite positive.")); + } Vector mv = m.Multiply(vector); const double raw_mass = vector.Dot(mv); if (!(raw_mass > 0.0) || !std::isfinite(raw_mass)) @@ -257,6 +262,12 @@ Result ExecuteCandidate(const Domain& domain) { Failure("solid-eigenpair-validation-failed", "Modal mass is not finite positive.")); vector.Scale(1.0 / std::sqrt(raw_mass)); + for (std::size_t row = 0; row < vector.Size(); ++row) { + if (!std::isfinite(vector[row])) + return Result::Failure( + Failure("solid-eigenpair-validation-failed", + "Normalized modal vector is not finite.")); + } mv = m.Multiply(vector); const Vector kv = k.Multiply(vector); const double lambda = state.backend.eigenvalues[c]; @@ -362,6 +373,10 @@ Result SolidFrequencyProcedureExecutor::Execute( return Result::Failure( Failure("eigen-allocation-failed", "Solid modal workspace exceeds addressable storage.")); + } catch (const std::runtime_error&) { + return Result::Failure( + Failure("eigen-multiplication-failed", + "Solid modal sparse multiplication failed.")); } } } // namespace fesa diff --git a/src/fesa/analysis/thermal_transient_algebra.cpp b/src/fesa/analysis/thermal_transient_algebra.cpp index 7552e83..2784d43 100644 --- a/src/fesa/analysis/thermal_transient_algebra.cpp +++ b/src/fesa/analysis/thermal_transient_algebra.cpp @@ -3,10 +3,12 @@ #include #include #include +#include #include #include #include "thermal_transient_bounds.h" +#include "math/dense_blas_internal.h" namespace fesa::thermal_transient_internal { namespace { @@ -116,6 +118,9 @@ Result Build(const ThermalSystem& s, const ThermalPartition& p, constrained[p.constrained[i]] = i; SparsePattern pattern{{0}, {}}; std::vector coo; + const auto entries = s.c.Values().size(); + Vector capacity_diagonal{entries}, conductivity_diagonal{entries}; + Vector capacity_input{entries}, conductivity_input{entries}; Vector b{n}; std::vector exact_a, exact_b(n); for (std::size_t f = 0; f < n; ++f) { @@ -126,22 +131,39 @@ Result Build(const ThermalSystem& s, const ThermalPartition& p, const auto col = s.c.ColumnIndices()[j]; const double c = s.c.Values()[j], k = s.k.Values()[j]; const I cd = Div(I{c}, I{delta}); + capacity_diagonal[j] = c / delta; if (free[col] != n) { + capacity_input[j] = old[col]; coo.push_back({f, free[col], c / delta + k, 0, coo.size()}); pattern.column_indices.push_back(free[col]); exact_a.push_back(Add(cd, I{k})); - b[f] += (c / delta) * old[col]; exact_b[f] = Add(exact_b[f], Mul(cd, I{old[col]})); } else { const double g = p.prescribed[constrained[col]]; + capacity_input[j] = old[col] - g; + conductivity_diagonal[j] = k; + conductivity_input[j] = g; // Preserve raw old constrained capacity; never project frame zero. - b[f] += (c / delta) * (old[col] - g) - k * g; exact_b[f] = Add(exact_b[f], Sub(Mul(cd, Sub(I{old[col]}, I{g})), Mul(I{k}, I{g}))); } } pattern.row_offsets.push_back(coo.size()); } + const auto capacity_action = dense_blas_internal::MultiplyDiagonal( + capacity_diagonal, capacity_input); + const auto conductivity_action = dense_blas_internal::MultiplyDiagonal( + conductivity_diagonal, conductivity_input); + for (std::size_t f = 0; f < n; ++f) { + const auto row = p.free[f]; + for (auto j = s.c.RowOffsets()[row]; j < s.c.RowOffsets()[row + 1]; ++j) { + const auto col = s.c.ColumnIndices()[j]; + // Preserve the original per-entry subtraction before the row sum. + // Separately summing C*delta and K*g loses small stable RHS terms. + b[f] += free[col] != n ? capacity_action[j] + : capacity_action[j] - conductivity_action[j]; + } + } auto a = SparseMatrix::FromCoo(n, n, std::move(coo), pattern); if (!a.HasValue()) return Result::Failure(a.GetStatus()); if (!Finite(b)) @@ -178,13 +200,15 @@ Status Check(const ThermalSystem& s, const ThermalPartition& p, I a_norm, t_norm, b_norm = Norm(eq.b); Vector red{n}, red_error{n}; std::vector construction(n); + Vector free_temperature{n}; + for (std::size_t f = 0; f < n; ++f) + free_temperature[f] = candidate.temperature[p.free[f]]; + const auto reduced_action = eq.a.Multiply(free_temperature); for (std::size_t f = 0; f < n; ++f) { I row_norm, scale, ea, identity; - double dot = 0; for (auto j = eq.a.RowOffsets()[f]; j < eq.a.RowOffsets()[f + 1]; ++j) { const double a = eq.a.Values()[j], t = candidate.temperature[p.free[eq.a.ColumnIndices()[j]]]; - dot += a * t; row_norm = Add(row_norm, Magnitude(a)); scale = Add(scale, Mul(Magnitude(a), Magnitude(t))); const I error{Error(eq.exact_a[j], a)}; @@ -193,7 +217,7 @@ Status Check(const ThermalSystem& s, const ThermalPartition& p, } a_norm = Max(a_norm, row_norm); t_norm = Max(t_norm, Magnitude(candidate.temperature[p.free[f]])); - red[f] = dot - eq.b[f]; + red[f] = reduced_action[f] - eq.b[f]; scale = Add(scale, Magnitude(eq.b[f])); const double nu = static_cast(eq.a.RowOffsets()[f + 1] - eq.a.RowOffsets()[f]); @@ -216,18 +240,19 @@ Status Check(const ThermalSystem& s, const ThermalPartition& p, for (std::size_t i = 0; i < m; ++i) change[i] = candidate.temperature[i] - old[i]; if (!Finite(change)) return Failure("Nonfinite temperature change."); + const auto capacity_change = s.c.Multiply(change); + const auto conductivity_new = s.k.Multiply(candidate.temperature); I c_change, k_new, q_abs, c_temperatures, underflow_sum; std::vector k_columns(m); double energy = 0; for (std::size_t i = 0; i < m; ++i) { - double cd = 0, kt = 0, ksum = 0; + const double cd = capacity_change[i], kt = conductivity_new[i]; + double ksum = 0; I cs, ks, ts, cabs, nullscale, nullerror, coefficient; for (auto j = s.c.RowOffsets()[i]; j < s.c.RowOffsets()[i + 1]; ++j) { const auto col = s.c.ColumnIndices()[j]; const double c = s.c.Values()[j], k = s.k.Values()[j], t = candidate.temperature[col]; - cd += c * change[col]; - kt += k * t; ksum += k; cs = Add(cs, Mul(Magnitude(c), Magnitude(change[col]))); ks = Add(ks, Mul(Magnitude(k), Magnitude(t))); @@ -334,10 +359,13 @@ Status Check(const ThermalSystem& s, const ThermalPartition& p, if (p.constrained.empty() && Zero(q_abs)) { ev.dissipation_checked = true; double cn = 0, co = 0, cd = 0, kn = 0; + const auto capacity_new = s.c.Multiply(candidate.temperature); + const auto capacity_old = s.c.Multiply(old); I sc, subscale, t_one, d_square, o_square, skew_square, quadratic_underflow; for (std::size_t i = 0; i < m; ++i) { const double t = candidate.temperature[i]; - double row_cn = 0, row_co = 0, row_cd = 0, row_kn = 0; + const double row_cn = capacity_new[i], row_co = capacity_old[i], + row_cd = capacity_change[i], row_kn = conductivity_new[i]; I scale_cn, scale_co, scale_cd, scale_kn; t_one = Add(t_one, Magnitude(t)); d_square = Add(d_square, Mul(Magnitude(change[i]), Magnitude(change[i]))); @@ -345,10 +373,6 @@ Status Check(const ThermalSystem& s, const ThermalPartition& p, for (auto j = s.c.RowOffsets()[i]; j < s.c.RowOffsets()[i + 1]; ++j) { const auto col = s.c.ColumnIndices()[j]; const double c = s.c.Values()[j], k = s.k.Values()[j]; - row_cn += c * candidate.temperature[col]; - row_co += c * old[col]; - row_cd += c * change[col]; - row_kn += k * candidate.temperature[col]; scale_cn = Add( scale_cn, Mul(Magnitude(c), Magnitude(candidate.temperature[col]))); scale_co = Add(scale_co, Mul(Magnitude(c), Magnitude(old[col]))); @@ -421,7 +445,7 @@ Status Check(const ThermalSystem& s, const ThermalPartition& p, Result SolveThermalIncrement( const ThermalSystem& s, const ThermalPartition& p, const Vector& old, - double delta, LinearSolver& solver, ThermalFactorization& cache) { + double delta, LinearSolver& solver, ThermalFactorization& cache) try { auto status = Validate(s, p, old, delta); if (!status.IsOk()) return Result::Failure(status); auto equations = Build(s, p, old, delta); @@ -461,6 +485,15 @@ Result SolveThermalIncrement( status = Check(s, p, old, delta, equations.Value(), candidate); if (!status.IsOk()) return Result::Failure(status); return Result::Success(std::move(candidate)); +} catch (const std::invalid_argument&) { + return Result::Failure( + Failure("Nonfinite backward Euler RHS.")); +} catch (const std::length_error&) { + return Result::Failure( + Failure("Thermal sparse action exceeds the backend integer range.")); +} catch (const std::runtime_error&) { + return Result::Failure( + Failure("Thermal sparse action failed.")); } Status RunThermalHistory(const ThermalAssembly& assembly, diff --git a/src/fesa/elements/euler_beam_3d.cpp b/src/fesa/elements/euler_beam_3d.cpp index ea43ced..82ca250 100644 --- a/src/fesa/elements/euler_beam_3d.cpp +++ b/src/fesa/elements/euler_beam_3d.cpp @@ -12,6 +12,7 @@ #include "fesa/math/vector3.h" #include "fesa/model/model_types.h" +#include "math/dense_blas_internal.h" namespace fesa { @@ -228,10 +229,9 @@ double NormalizedMatrixError(const Matrix& lhs, const Matrix& rhs) { std::array GeneralizedStrain( const Matrix& b, const Vector& local_displacement) { std::array strain{}; + const auto product = b.Multiply(local_displacement); for (std::size_t component = 0; component < strain.size(); ++component) { - for (std::size_t dof = 0; dof < local_displacement.Size(); ++dof) { - strain[component] += b(component, dof) * local_displacement[dof]; - } + strain[component] = product[component]; } return strain; } @@ -395,13 +395,15 @@ Matrix EulerBeam3D::LocalStiffness() const { // Both Gauss points are required: a one-point bending rule loses two ranks. for (const double xi : gauss_points) { const Matrix b = StrainDisplacement(xi, length_); + Matrix weighted_b = b; + for (std::size_t component = 0; component < diagonal.size(); ++component) + for (std::size_t column = 0; column < kElementDofCount; ++column) + weighted_b(component, column) *= diagonal[component]; + const auto contribution = + dense_blas_internal::MultiplyTransposeLeft(b, weighted_b); for (std::size_t row = 0; row < kElementDofCount; ++row) { for (std::size_t column = 0; column < kElementDofCount; ++column) { - for (std::size_t component = 0; component < diagonal.size(); - ++component) { - stiffness(row, column) += b(component, row) * diagonal[component] * - b(component, column) * jacobian; - } + stiffness(row, column) += contribution(row, column) * jacobian; } } } @@ -418,22 +420,17 @@ Matrix EulerBeam3D::GlobalStiffness() const { const Matrix local = LocalStiffness(); const Matrix transform = Transformation(rotation_); const Matrix local_times_transform = local.Multiply(transform); - Matrix global{kElementDofCount, kElementDofCount}; // Kg=T^T*Kl*T while dl=T*dg. - for (std::size_t row = 0; row < kElementDofCount; ++row) { - for (std::size_t column = 0; column < kElementDofCount; ++column) { - for (std::size_t inner = 0; inner < kElementDofCount; ++inner) { - global(row, column) += - transform(inner, row) * local_times_transform(inner, column); - } - } - } - return global; + return dense_blas_internal::MultiplyTransposeLeft(transform, + local_times_transform); } Vector EulerBeam3D::LocalEquivalentLoad( const ConstantLocalLineLoad& load) const { const std::array components = {load.px, load.py, load.pz, load.mx}; + Matrix component_column{components.size(), 1U}; + for (std::size_t component = 0; component < components.size(); ++component) + component_column(component, 0U) = components[component]; Vector equivalent{kElementDofCount}; const double inverse_sqrt_three = 1.0 / std::sqrt(3.0); const std::array gauss_points = {-inverse_sqrt_three, @@ -441,12 +438,10 @@ Vector EulerBeam3D::LocalEquivalentLoad( const double jacobian = 0.5 * length_; for (const double xi : gauss_points) { const Matrix interpolation = KinematicInterpolation(xi, length_); + const auto contribution = dense_blas_internal::MultiplyTransposeLeft( + interpolation, component_column); for (std::size_t dof = 0; dof < equivalent.Size(); ++dof) { - for (std::size_t component = 0; component < components.size(); - ++component) { - equivalent[dof] += - interpolation(component, dof) * components[component] * jacobian; - } + equivalent[dof] += contribution(dof, 0U) * jacobian; } } return equivalent; diff --git a/src/fesa/math/dense_blas_internal.cpp b/src/fesa/math/dense_blas_internal.cpp index e19948c..413a71c 100644 --- a/src/fesa/math/dense_blas_internal.cpp +++ b/src/fesa/math/dense_blas_internal.cpp @@ -1,10 +1,101 @@ #include "math/dense_blas_internal.h" +#include #include #include +#include +#include namespace fesa::dense_blas_internal { +ScopedSingleEquationThreads::ScopedSingleEquationThreads( + const std::size_t equation_count) noexcept { + // The measured 1x1 PARDISO case has backend/thread overhead and an + // intermittently failing default-thread lifecycle. Broader sizes are not + // inferred from this auxiliary evidence. + if (equation_count != 1U) return; + const int previous = mkl_set_num_threads_local(1); + if (previous == 0) { + restore_default_ = true; + } else { + static_cast(mkl_set_num_threads_local(previous)); + } +} + +ScopedSingleEquationThreads::~ScopedSingleEquationThreads() { + if (restore_default_) static_cast(mkl_set_num_threads_local(0)); +} + +bool CholeskyLowerAndInverse(const Matrix& matrix, Matrix& lower, + Matrix& inverse_lower) { + if (matrix.Rows() != matrix.Columns()) + throw std::invalid_argument{"Cholesky requires a square matrix."}; + const auto size = ToMklSize(matrix.Rows()); + if (!size.HasValue()) + throw std::length_error{"Cholesky size exceeds the MKL integer range."}; + Matrix candidate_lower{matrix.Rows(), matrix.Columns()}; + Matrix candidate_inverse{matrix.Rows(), matrix.Columns()}; + if (matrix.Rows() == 0U) { + lower = std::move(candidate_lower); + inverse_lower = std::move(candidate_inverse); + return true; + } + std::vector values(matrix.Rows() * matrix.Columns()); + for (std::size_t row = 0; row < matrix.Rows(); ++row) + for (std::size_t column = 0; column < matrix.Columns(); ++column) { + if (!std::isfinite(matrix(row, column))) + throw std::invalid_argument{"Cholesky inputs must be finite."}; + values[row * matrix.Columns() + column] = matrix(row, column); + } + if (LAPACKE_dpotrf(LAPACK_ROW_MAJOR, 'L', size.Value(), values.data(), + size.Value()) != 0) return false; + for (std::size_t row = 0; row < matrix.Rows(); ++row) + for (std::size_t column = 0; column <= row; ++column) + candidate_lower(row, column) = values[row * matrix.Columns() + column]; + if (LAPACKE_dtrtri(LAPACK_ROW_MAJOR, 'L', 'N', size.Value(), values.data(), + size.Value()) != 0) 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)) + return false; + candidate_inverse(row, column) = value; + } + lower = std::move(candidate_lower); + inverse_lower = std::move(candidate_inverse); + return true; +} + +Vector MultiplyDiagonal(const Vector& diagonal, const Vector& rhs) { + if (diagonal.Size() != rhs.Size()) + throw std::invalid_argument{"Diagonal product dimensions differ."}; + const auto size = ToMklSize(diagonal.Size()); + if (!size.HasValue()) + throw std::length_error{"Diagonal size exceeds the MKL integer range."}; + for (std::size_t i = 0; i < rhs.Size(); ++i) { + if (!std::isfinite(diagonal[i]) || !std::isfinite(rhs[i])) + throw std::invalid_argument{"Diagonal product inputs must be finite."}; + } + Vector candidate{rhs.Size()}; + if (rhs.Size() == 0U) return candidate; + // k=0 stores exactly one coefficient per independent row. MKL performs + // each product without regrouping a caller's subsequent certified sums. + cblas_dsbmv(CblasColMajor, CblasUpper, size.Value(), 0, 1.0, diagonal.Data(), + 1, rhs.Data(), 1, 0.0, candidate.Data(), 1); + return candidate; +} + +Matrix MultiplyTransposeLeft(const Matrix& left, const Matrix& right) { + if (left.Rows() != right.Rows()) { + throw std::invalid_argument{"Transpose product dimensions differ."}; + } + Matrix transposed{left.Columns(), left.Rows()}; + for (std::size_t row = 0; row < left.Rows(); ++row) + for (std::size_t column = 0; column < left.Columns(); ++column) + transposed(column, row) = left(row, column); + return transposed.Multiply(right); +} + 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 056a20a..78f9e77 100644 --- a/src/fesa/math/dense_blas_internal.h +++ b/src/fesa/math/dense_blas_internal.h @@ -6,13 +6,47 @@ #include #include "fesa/core/status.h" +#include "fesa/math/matrix.h" +#include "fesa/math/vector.h" namespace fesa::dense_blas_internal { +/// @brief Limits only default-policy one-equation backend calls to one thread. +/// @note Existing caller TLS overrides are retained. Destruction restores the +/// default request; process-wide MKL/OpenMP settings are never changed. +class ScopedSingleEquationThreads { + public: + explicit ScopedSingleEquationThreads(std::size_t equation_count) noexcept; + ScopedSingleEquationThreads(const ScopedSingleEquationThreads&) = delete; + ScopedSingleEquationThreads& operator=(const ScopedSingleEquationThreads&) = + delete; + ~ScopedSingleEquationThreads(); + + 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); +/// @brief Computes left-transpose times right through ordinary dense MKL. +/// @throws std::invalid_argument for mismatched rows or nonfinite operands. +Matrix MultiplyTransposeLeft(const Matrix& left, const Matrix& right); + +/// @brief Multiplies a compact diagonal matrix by a vector through MKL. +/// @throws std::invalid_argument if shapes differ or inputs are nonfinite. +/// @throws std::length_error if the size exceeds the backend integer range. +Vector MultiplyDiagonal(const Vector& diagonal, const Vector& rhs); + +/// @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. +/// @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); + /// @brief Copies contiguous double values through the private BLAS backend. /// @param source Readable storage for size values when size is nonzero. /// @param size Number of contiguous values to copy. diff --git a/src/fesa/math/high_low_arithmetic.h b/src/fesa/math/high_low_arithmetic.h index c724d1e..05ccad6 100644 --- a/src/fesa/math/high_low_arithmetic.h +++ b/src/fesa/math/high_low_arithmetic.h @@ -4,6 +4,8 @@ #include #include #include +#include +#include #include #if defined(_M_X64) || defined(__SSE2__) #include @@ -56,9 +58,9 @@ inline bool IsNormalized(HighLow value) { if (!IsFinite(value)) return false; if (value.high == 0 && value.low == 0) return !std::signbit(value.high) && !std::signbit(value.low); - const auto normalized = TwoSum(value.high, value.low); - return IsFinite(normalized) && normalized.high == value.high && - normalized.low == value.low; + // If RN(high + low) == high, TwoSum's bridge is zero and its error is + // exactly low. The finite and canonical exact-zero guards remain above. + return value.high + value.low == value.high; } /// @brief Rounds the full physical value to one finite output value. inline double Rounded(HighLow value) { @@ -92,15 +94,27 @@ inline HighLow Subtract(HighLow left, HighLow right) { /// @brief Returns the exponent of the least nonzero binary64 significand bit. /// @pre value is finite and nonzero; gradual underflow is enabled. inline int LowestBitExponent(double value) { - int exponent = 0; - const double fraction = std::frexp(std::abs(value), &exponent); - auto significand = static_cast(std::ldexp(fraction, 53)); + static_assert(sizeof(double) == sizeof(std::uint64_t) && + std::numeric_limits::is_iec559 && + std::numeric_limits::digits == 53); + std::uint64_t bits = 0; + std::memcpy(&bits, &value, sizeof(bits)); + const auto encoded_exponent = (bits >> 52U) & 0x7ffU; + auto significand = bits & 0xfffffffffffffULL; + const int exponent = encoded_exponent == 0U + ? -1074 + : static_cast(encoded_exponent) - 1075; + if (encoded_exponent != 0U) significand |= 1ULL << 52U; int trailing = 0; - while ((significand & 1U) == 0U) { - significand >>= 1U; - ++trailing; + // Strip large zero groups first without vendor/compiler intrinsics. + for (const int shift : {32, 16, 8, 4, 2, 1}) { + const auto mask = (1ULL << shift) - 1U; + if ((significand & mask) == 0U) { + significand >>= shift; + trailing += shift; + } } - return exponent - 53 + trailing; + return exponent + trailing; } /// @brief Gives a product and its FMA residual when the dyadic tail fits. /// @note Rejects underflow-risk exact tails, even if the rounded high is finite. diff --git a/src/fesa/math/matrix.cpp b/src/fesa/math/matrix.cpp index fa9ce89..a28e167 100644 --- a/src/fesa/math/matrix.cpp +++ b/src/fesa/math/matrix.cpp @@ -1,5 +1,6 @@ #include "fesa/math/matrix.h" +#include #include #include @@ -107,6 +108,12 @@ Vector Matrix::Multiply(const Vector& rhs) const { throw std::invalid_argument{ "Matrix-vector multiplication has incompatible dimensions."}; } + for (const double value : values_) + if (!std::isfinite(value)) + throw std::invalid_argument{"Matrix product inputs must be finite."}; + for (std::size_t i = 0; i < rhs.Size(); ++i) + if (!std::isfinite(rhs.Data()[i])) + throw std::invalid_argument{"Matrix product inputs must be finite."}; Vector result{rows_}; if (rows_ == 0 || columns_ == 0) { @@ -126,6 +133,12 @@ Matrix Matrix::Multiply(const Matrix& rhs) const { throw std::invalid_argument{ "Matrix multiplication has incompatible dimensions."}; } + for (const double value : values_) + if (!std::isfinite(value)) + throw std::invalid_argument{"Matrix product inputs must be finite."}; + for (const double value : rhs.values_) + if (!std::isfinite(value)) + throw std::invalid_argument{"Matrix product inputs must be finite."}; Matrix result{rows_, rhs.columns_}; if (rows_ == 0 || columns_ == 0 || rhs.columns_ == 0) { diff --git a/src/fesa/math/sparse_matrix.cpp b/src/fesa/math/sparse_matrix.cpp index f0ff8c8..988964e 100644 --- a/src/fesa/math/sparse_matrix.cpp +++ b/src/fesa/math/sparse_matrix.cpp @@ -8,11 +8,49 @@ #include #include +#include + #include "fesa/fem/dof_manager.h" +#include "math/dense_blas_internal.h" namespace fesa { namespace { +/// @brief Checks the actual private backend integer width before conversion. +MKL_INT SparseSize(const std::size_t size) { + const auto converted = dense_blas_internal::ToMklSize(size); + if (!converted.HasValue()) { + throw std::length_error{"Sparse size exceeds the MKL integer range."}; + } + return converted.Value(); +} + +/// @brief Releases a call-local sparse handle before its borrowed arrays die. +class SparseHandle { + public: + SparseHandle() = default; + SparseHandle(const SparseHandle&) = delete; + SparseHandle& operator=(const SparseHandle&) = delete; + ~SparseHandle() { + if (handle != nullptr) static_cast(mkl_sparse_destroy(handle)); + } + /// @brief Reports normal release failure before publishing the candidate. + sparse_status_t Destroy() { + const auto owned = handle; + handle = nullptr; + return mkl_sparse_destroy(owned); + } + sparse_matrix_t handle{nullptr}; +}; + +/// @brief Fails explicitly instead of substituting a scalar backend. +void CheckSparseStatus(const sparse_status_t status, const char* operation) { + if (status != SPARSE_STATUS_SUCCESS) { + throw std::runtime_error{std::string{"Sparse backend "} + operation + + " failed."}; + } +} + /// @brief Builds a structured sparse-matrix model failure. Status SparseFailure(const std::string& code, const std::string& identity, const std::string& message) { @@ -170,16 +208,47 @@ Vector SparseMatrix::Multiply(const Vector& rhs) const { throw std::invalid_argument{ "Sparse matrix-vector multiplication has incompatible dimensions."}; } - - Vector result{rows_}; - for (std::size_t row = 0U; row < rows_; ++row) { - double value = 0.0; - for (std::size_t position = row_offsets_[row]; - position < row_offsets_[row + 1U]; ++position) { - value += values_[position] * rhs[column_indices_[position]]; - } - result[row] = value; + if (!Validate().IsOk()) { + throw std::invalid_argument{"Sparse multiplication requires valid CSR."}; } + for (std::size_t i = 0; i < rhs.Size(); ++i) { + if (!std::isfinite(rhs[i])) + throw std::invalid_argument{"Sparse multiplication RHS must be finite."}; + } + const MKL_INT rows = SparseSize(rows_); + const MKL_INT columns = SparseSize(columns_); + static_cast(SparseSize(values_.size())); + Vector result{rows_}; + if (rows == 0 || columns == 0 || values_.empty()) return result; + + // Each const call owns its converted arrays and handle. Neither copy/move + // nor concurrent calls retain pointers into another call's storage. + std::vector row_offsets; + std::vector column_indices; + row_offsets.reserve(row_offsets_.size()); + column_indices.reserve(column_indices_.size()); + for (const auto offset : row_offsets_) { + row_offsets.push_back(SparseSize(offset)); + } + for (const auto column : column_indices_) { + column_indices.push_back(SparseSize(column)); + } + auto values = values_; + SparseHandle sparse; + CheckSparseStatus( + mkl_sparse_d_create_csr(&sparse.handle, SPARSE_INDEX_BASE_ZERO, rows, + columns, row_offsets.data(), + row_offsets.data() + 1, column_indices.data(), + values.data()), + "create CSR"); + const matrix_descr descriptor{SPARSE_MATRIX_TYPE_GENERAL, + SPARSE_FILL_MODE_FULL, + SPARSE_DIAG_NON_UNIT}; + CheckSparseStatus( + mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, 1.0, sparse.handle, + descriptor, rhs.Data(), 0.0, result.Data()), + "multiply"); + CheckSparseStatus(sparse.Destroy(), "destroy CSR"); return result; } diff --git a/src/fesa/results/global_equilibrium_recovery.cpp b/src/fesa/results/global_equilibrium_recovery.cpp index 1fa096b..3f36bca 100644 --- a/src/fesa/results/global_equilibrium_recovery.cpp +++ b/src/fesa/results/global_equilibrium_recovery.cpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include @@ -26,33 +27,25 @@ double IndexedNorm(const Vector& values, return result; } -/// @brief Computes separate Kff*df and Kfc*dc norms in stable CSR order. +/// @brief Computes separate Kff*df and Kfc*dc norms through sparse MKL. std::array FreeEquationInternalTermNorms( const SparseMatrix& stiffness, const Vector& displacement, const DofManager& dofs) { - std::vector free_column(dofs.FullDofCount(), 0U); + Vector free_displacement{displacement.Size()}; + Vector constrained_displacement = displacement; for (const std::size_t full_dof : dofs.FreeDofs()) { - free_column[full_dof] = 1U; + free_displacement[full_dof] = displacement[full_dof]; + constrained_displacement[full_dof] = 0.0; } - + const auto free_term = stiffness.Multiply(free_displacement); + const auto constrained_term = stiffness.Multiply(constrained_displacement); double free_term_norm = 0.0; double constrained_term_norm = 0.0; for (const std::size_t row : dofs.FreeDofs()) { - double free_term = 0.0; - double constrained_term = 0.0; - for (std::size_t position = stiffness.RowOffsets()[row]; - position < stiffness.RowOffsets()[row + 1U]; ++position) { - const std::size_t column = stiffness.ColumnIndices()[position]; - const double contribution = - stiffness.Values()[position] * displacement[column]; - if (free_column[column] != 0U) { - free_term += contribution; - } else { - constrained_term += contribution; - } - } - free_term_norm = std::hypot(free_term_norm, free_term); - constrained_term_norm = std::hypot(constrained_term_norm, constrained_term); + // Only the original free-row terms enter these acceptance scales. + free_term_norm = std::hypot(free_term_norm, free_term[row]); + constrained_term_norm = + std::hypot(constrained_term_norm, constrained_term[row]); } return {free_term_norm, constrained_term_norm}; } @@ -96,7 +89,7 @@ double NormalizedBalance(const std::array& balance, Result GlobalEquilibriumRecovery( const Domain& domain, const DofManager& dofs, const SparseMatrix& full_stiffness, const Vector& full_displacement, - const Vector& full_external_force) { + const Vector& full_external_force) try { const std::size_t full_count = dofs.FullDofCount(); if (full_stiffness.Rows() != full_count || full_stiffness.Columns() != full_count || @@ -171,6 +164,15 @@ Result GlobalEquilibriumRecovery( // shell global equilibrium below separately filters constrained reactions. result.reaction = result.residual; return Result::Success(std::move(result)); +} catch (const std::length_error&) { + return RecoveryResultFailure( + "recovery-index-overflow", {domain.SourcePath(), 0U}, + domain.SourceContentIdentity(), + "Recovery multiplication exceeds the backend integer range."); +} catch (const std::runtime_error&) { + return RecoveryResultFailure( + "recovery-multiplication-failed", {domain.SourcePath(), 0U}, + domain.SourceContentIdentity(), "Recovery sparse multiplication failed."); } Status RecoverShellGlobalEvidence(const Domain& domain, const DofManager& dofs, diff --git a/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp b/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp index af2179f..7e3a237 100644 --- a/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp +++ b/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp @@ -12,14 +12,17 @@ #include #include #include +#include #include #include #include #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 { @@ -186,6 +189,37 @@ Result CopyUpperOneBased(const SparseMatrix& matrix, const double scale) { return Result::Success(std::move(copy)); } +/// @brief Keeps the independently scaled full or upper-mirrored operator. +Result ScaledActionMatrix(const SparseMatrix& matrix, + const double scale, + const bool mirror_upper = false) { + std::vector contributions; + std::vector> 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 IdentityOneBased(const std::size_t size) { if (!FitsMkl(size) || !FitsMkl(size + 1U)) { return Result::Failure( @@ -321,13 +355,11 @@ bool FillInitialBasis(std::vector& vectors, const std::size_t n, } for (double& value : reflector) value /= reflector_norm; for (std::size_t target = column; target < columns; ++target) { - double dot = 0.0; - for (std::size_t row = column; row < n; ++row) { - dot += reflector[row - column] * vectors[row + target * n]; - } - for (std::size_t row = column; row < n; ++row) { - vectors[row + target * n] -= 2.0 * reflector[row - column] * dot; - } + const auto length = static_cast(n - column); + auto* target_column = vectors.data() + column + target * n; + const double dot = cblas_ddot(length, reflector.data(), 1, + target_column, 1); + cblas_daxpy(length, -2.0 * dot, reflector.data(), 1, target_column, 1); } reflectors.push_back(std::move(reflector)); } @@ -338,13 +370,11 @@ bool FillInitialBasis(std::vector& vectors, const std::size_t n, const std::size_t column = index - 1U; const auto& reflector = reflectors[column]; for (std::size_t target = 0U; target < columns; ++target) { - double dot = 0.0; - for (std::size_t row = column; row < n; ++row) { - dot += reflector[row - column] * vectors[row + target * n]; - } - for (std::size_t row = column; row < n; ++row) { - vectors[row + target * n] -= 2.0 * reflector[row - column] * dot; - } + const auto length = static_cast(n - column); + auto* target_column = vectors.data() + column + target * n; + const double dot = cblas_ddot(length, reflector.data(), 1, + target_column, 1); + cblas_daxpy(length, -2.0 * dot, reflector.data(), 1, target_column, 1); } } for (std::size_t column = 0U; column < columns; ++column) @@ -353,10 +383,9 @@ bool FillInitialBasis(std::vector& vectors, const std::size_t n, double gram_error = 0.0; for (std::size_t left = 0U; left < columns; ++left) { for (std::size_t right = 0U; right < columns; ++right) { - double dot = 0.0; - for (std::size_t row = 0U; row < n; ++row) { - dot += vectors[row + left * n] * vectors[row + right * n]; - } + const double dot = cblas_ddot(static_cast(n), + vectors.data() + left * n, 1, + vectors.data() + right * n, 1); gram_error = std::hypot(gram_error, dot - (left == right ? 1.0 : 0.0)); } } @@ -390,7 +419,7 @@ Result MklFeastEigenSolver::Solve( Result 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 { + 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(); @@ -472,6 +501,9 @@ Result MklFeastEigenSolver::SolveImpl( auto b = mass != nullptr ? CopyUpperOneBased(*mass, 1.0) : IdentityOneBased(n); if (!b.HasValue()) return Result::Failure(b.GetStatus()); + auto action_matrix = ScaledActionMatrix(whitened_operator, operator_scale); + if (!action_matrix.HasValue()) + return Result::Failure(action_matrix.GetStatus()); double lower = -1.0; double upper = 1.0; @@ -1045,13 +1077,11 @@ Result MklFeastEigenSolver::SolveImpl( "eigen-endpoint-separation", std::to_string(input), "A returned FEAST eigenvalue is too close to a solve endpoint.")); } - double norm = 0.0; + Vector raw(n); for (std::size_t row = 0U; row < n; ++row) - norm += vectors[row + input * n] * vectors[row + input * n]; + raw[row] = vectors[row + input * n]; + double norm = raw.Dot(raw); if (mass != nullptr) { - Vector raw(n); - for (std::size_t row = 0U; row < n; ++row) - raw[row] = vectors[row + input * n]; norm = raw.Dot(mass->Multiply(raw)); } norm = std::sqrt(norm); @@ -1060,23 +1090,16 @@ Result MklFeastEigenSolver::SolveImpl( 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) { - double ax = 0.0; - for (std::size_t p = offsets[row]; p < offsets[row + 1U]; ++p) { - ax += values[p] / operator_scale * - vectors[whitened_operator.ColumnIndices()[p] + input * n]; - } - double mx = vectors[row + input * n]; - if (mass != nullptr) { - mx = 0.0; - for (std::size_t p = mass->RowOffsets()[row]; - p < mass->RowOffsets()[row + 1U]; ++p) - mx += - mass->Values()[p] * vectors[mass->ColumnIndices()[p] + input * n]; - } - const double difference = ax - mu * mx; + 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::Failure( + TraceFailure("eigenpair-validation-failed", std::to_string(input), + "FEAST returned a non-normalizable eigenvector.")); } double residual_denominator = norm; if (mass != nullptr) { @@ -1109,20 +1132,11 @@ Result MklFeastEigenSolver::SolveImpl( FormatDouble(rho_b) + ", mu=" + FormatDouble(mu) + "; calls=" + CallHistory(calls))); const FeastBackendCall& last_call = calls.back(); - std::vector mirrored_action(n, 0.0); - std::vector full_action(n, 0.0); - for (std::size_t row = 0U; row < n; ++row) { - for (std::size_t p = offsets[row]; p < offsets[row + 1U]; ++p) { - const std::size_t column = whitened_operator.ColumnIndices()[p]; - const double entry = values[p] / operator_scale; - full_action[row] += entry * vectors[column + input * n]; - if (column < row) continue; - mirrored_action[row] += entry * vectors[column + input * n]; - if (column != row) { - mirrored_action[column] += entry * vectors[row + input * n]; - } - } - } + auto mirrored_matrix = + ScaledActionMatrix(whitened_operator, operator_scale, true); + if (!mirrored_matrix.HasValue()) + return Result::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) { @@ -1130,7 +1144,7 @@ Result MklFeastEigenSolver::SolveImpl( mirrored_action[row] - mu * vectors[row + input * n]; mirrored_residual_squared += mirrored_residual * mirrored_residual; const double action_difference = - full_action[row] - mirrored_action[row]; + action[row] - mirrored_action[row]; action_difference_squared += action_difference * action_difference; } const double rho_b_upper = std::sqrt(mirrored_residual_squared) / norm; @@ -1174,22 +1188,25 @@ Result MklFeastEigenSolver::SolveImpl( } 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::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) { - double dot = 0.0; - for (std::size_t row = 0U; row < n; ++row) - dot += ordered_vectors(row, left) * ordered_vectors(row, right); - if (mass != nullptr) { - dot = 0.0; - for (std::size_t row = 0U; row < n; ++row) { - double action = 0.0; - for (std::size_t p = mass->RowOffsets()[row]; - p < mass->RowOffsets()[row + 1U]; ++p) - action += mass->Values()[p] * - ordered_vectors(mass->ColumnIndices()[p], right); - dot += ordered_vectors(row, left) * action; - } - } + 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; @@ -1233,6 +1250,18 @@ Result MklFeastEigenSolver::SolveImpl( solve_upper, upper_search_termination}; return Result::Success(std::move(result)); +} catch (const std::invalid_argument&) { + return Result::Failure( + Failure("eigenpair-validation-failed", "independent-action", + "Modal product inputs must be finite and shape compatible.")); +} catch (const std::length_error&) { + return Result::Failure( + Failure("feast-integer-overflow", "independent-action", + "Independent sparse action exceeds the backend range.")); +} catch (const std::runtime_error&) { + return Result::Failure( + Failure("eigen-multiplication-failed", "independent-action", + "Independent sparse action failed.")); } } // namespace fesa diff --git a/src/fesa/solvers/linear/mkl_pardiso_solver.cpp b/src/fesa/solvers/linear/mkl_pardiso_solver.cpp index add8672..0cc9adc 100644 --- a/src/fesa/solvers/linear/mkl_pardiso_solver.cpp +++ b/src/fesa/solvers/linear/mkl_pardiso_solver.cpp @@ -14,6 +14,7 @@ #include "fesa/math/sparse_matrix.h" #include "fesa/math/vector.h" +#include "math/dense_blas_internal.h" namespace fesa { namespace { @@ -286,6 +287,8 @@ class MklPardisoSolver::Impl { /// @brief Calls PARDISO with retained private arrays and phase state. void CallPardiso(const MKL_INT phase, double* rhs, double* solution, MKL_INT& error) { + const dense_blas_internal::ScopedSingleEquationThreads thread_scope{ + static_cast(equation_count_)}; pardiso(pt_.data(), &max_factorizations_, &matrix_number_, &mtype_, &phase, &equation_count_, values_.data(), row_offsets_.data(), column_indices_.data(), permutation_.data(), &rhs_count_, diff --git a/src/fesa/solvers/linear/mkl_pardiso_tangent_solver.cpp b/src/fesa/solvers/linear/mkl_pardiso_tangent_solver.cpp index 1554b67..b544312 100644 --- a/src/fesa/solvers/linear/mkl_pardiso_tangent_solver.cpp +++ b/src/fesa/solvers/linear/mkl_pardiso_tangent_solver.cpp @@ -13,6 +13,7 @@ #include "fesa/math/sparse_matrix.h" #include "fesa/math/vector.h" +#include "math/dense_blas_internal.h" namespace fesa { namespace { @@ -230,6 +231,8 @@ class MklPardisoTangentSolver::Impl { /// @brief Calls PARDISO with retained private arrays and phase state. void CallPardiso(const MKL_INT phase, double* rhs, double* correction, MKL_INT& error) { + const dense_blas_internal::ScopedSingleEquationThreads thread_scope{ + static_cast(equation_count_)}; pardiso(pt_.data(), &max_factorizations_, &matrix_number_, &mtype_, &phase, &equation_count_, values_.data(), row_offsets_.data(), column_indices_.data(), permutation_.data(), &rhs_count_, diff --git a/tests/unit/analysis/eigen_reduction_test.cpp b/tests/unit/analysis/eigen_reduction_test.cpp index bc31dd3..a7f9108 100644 --- a/tests/unit/analysis/eigen_reduction_test.cpp +++ b/tests/unit/analysis/eigen_reduction_test.cpp @@ -6,6 +6,7 @@ #include #include #include +#include #include #include @@ -92,6 +93,38 @@ TEST(EigenReduction, RemovesOnlyActualFreeDirectorNull) { .HasValue()); } +TEST(EigenReduction, RejectsNonfiniteUsedBasisAndPreservesPriorCandidate) { + const fesa::ShellNodeInitialFrame frame{0U, {0, 0, 1}, {1, 0, 0}, {0, 1, 0}}; + const auto mass = NodalMass(frame.director); + fesa::Matrix stiffness(6U, 6U); + for (std::size_t i = 0; i < 5U; ++i) stiffness(i, i) = mass(i, i) * 7.0; + stiffness(5U, 5U) = .25; + const auto matrix = Sparse(stiffness); + const std::vector mask(6U, 0U); + const auto prior = fesa::eigen_internal::BuildEigenReduction( + matrix, {mass}, {frame}, mask, {.25}, 2.0); + ASSERT_TRUE(prior.HasValue()); + const auto prior_values = prior.Value().whitened_operator.Values(); + for (const double invalid : {std::numeric_limits::quiet_NaN(), + std::numeric_limits::infinity()}) { + for (const bool first_tangent : {false, true}) { + auto bad_frame = frame; + (first_tangent ? bad_frame.tangent_a : bad_frame.tangent_b)[0] = invalid; + const auto rejected = fesa::eigen_internal::BuildEigenReduction( + matrix, {mass}, {bad_frame}, mask, {.25}, 2.0); + ASSERT_FALSE(rejected.HasValue()); + EXPECT_EQ(rejected.GetStatus().Category(), fesa::FailureCategory::kSolver); + EXPECT_EQ(rejected.GetStatus().Diagnostics()[0].code, + "invalid-eigen-mass-space"); + EXPECT_EQ(prior.Value().whitened_operator.Values(), prior_values); + } + } + const auto rerun = fesa::eigen_internal::BuildEigenReduction( + matrix, {mass}, {frame}, mask, {.25}, 2.0); + ASSERT_TRUE(rerun.HasValue()); + EXPECT_EQ(rerun.Value().whitened_operator.Values(), prior_values); +} + TEST(EigenReduction, RejectsPositiveButUnresolvableMassWithoutDeletingDofs) { const double p = 1e-8; const std::array d{2 * p / (1 + p * p), 0, diff --git a/tests/unit/math/dense_blas_internal_test.cpp b/tests/unit/math/dense_blas_internal_test.cpp index d036800..eaebe46 100644 --- a/tests/unit/math/dense_blas_internal_test.cpp +++ b/tests/unit/math/dense_blas_internal_test.cpp @@ -3,8 +3,10 @@ #include #include +#include #include #include +#include namespace fesa::dense_blas_internal { namespace { @@ -35,5 +37,114 @@ TEST(DenseBlasInternal, CDup003CopiesZeroAndNonzeroContiguousValues) { EXPECT_EQ(destination, source); } +TEST(DenseBlasInternal, TransposeProductChecksShapeAndKnownRectangularAction) { + Matrix left{3U, 2U}; + Matrix right{3U, 1U}; + for (std::size_t row = 0; row < 3U; ++row) { + left(row, 0U) = static_cast(row + 1U); + left(row, 1U) = 2.0; + right(row, 0U) = static_cast(row); + } + const auto product = MultiplyTransposeLeft(left, right); + ASSERT_EQ(product.Rows(), 2U); + ASSERT_EQ(product.Columns(), 1U); + EXPECT_DOUBLE_EQ(product(0U, 0U), 8.0); + EXPECT_DOUBLE_EQ(product(1U, 0U), 6.0); + EXPECT_THROW(static_cast(MultiplyTransposeLeft(left, Matrix{2U, 1U})), + std::invalid_argument); + const auto empty = MultiplyTransposeLeft(Matrix{0U, 2U}, Matrix{0U, 3U}); + EXPECT_EQ(empty.Rows(), 2U); + EXPECT_EQ(empty.Columns(), 3U); + EXPECT_DOUBLE_EQ(empty(1U, 2U), 0.0); + left(0U, 0U) = std::numeric_limits::quiet_NaN(); + EXPECT_THROW(static_cast(MultiplyTransposeLeft(left, right)), + std::invalid_argument); + left(0U, 0U) = 1.0; + right(0U, 0U) = std::numeric_limits::infinity(); + EXPECT_THROW(static_cast(MultiplyTransposeLeft(left, right)), + std::invalid_argument); +} + +TEST(DenseBlasInternal, DiagonalProductKeepsEntrywiseCancellationAndChecks) { + Vector diagonal{3U}; + diagonal[0U] = 4; + diagonal[1U] = 2; + diagonal[2U] = -1; + Vector input{3U}; + const double large = std::ldexp(1.0, 56); + input[0U] = 1; + input[1U] = -large; + input[2U] = 2 * large; + const auto product = MultiplyDiagonal(diagonal, input); + EXPECT_DOUBLE_EQ(product[0U] + (product[1U] - product[2U]), 4); + EXPECT_EQ(MultiplyDiagonal(Vector{0U}, Vector{0U}).Size(), 0U); + EXPECT_THROW(static_cast(MultiplyDiagonal(diagonal, Vector{2U})), + std::invalid_argument); + input[0U] = std::numeric_limits::infinity(); + EXPECT_THROW(static_cast(MultiplyDiagonal(diagonal, input)), + std::invalid_argument); + input[0U] = 1; + diagonal[0U] = std::numeric_limits::quiet_NaN(); + EXPECT_THROW(static_cast(MultiplyDiagonal(diagonal, input)), + std::invalid_argument); +} + +TEST(DenseBlasInternal, CholeskyAndInverseMatchKnownFactorsAndKeepFailedOutput) { + Matrix matrix{2U, 2U}; + matrix(0U, 0U) = 4; + matrix(0U, 1U) = 2; + matrix(1U, 0U) = 2; + matrix(1U, 1U) = 5; + Matrix lower{0U, 0U}, inverse{0U, 0U}; + ASSERT_TRUE(CholeskyLowerAndInverse(matrix, lower, inverse)); + EXPECT_DOUBLE_EQ(lower(0U, 0U), 2); + EXPECT_DOUBLE_EQ(lower(1U, 0U), 1); + EXPECT_DOUBLE_EQ(lower(1U, 1U), 2); + EXPECT_DOUBLE_EQ(lower(0U, 1U), 0); + EXPECT_DOUBLE_EQ(inverse(0U, 0U), .5); + EXPECT_DOUBLE_EQ(inverse(1U, 0U), -.25); + EXPECT_DOUBLE_EQ(inverse(1U, 1U), .5); + matrix(1U, 1U) = 0; + EXPECT_FALSE(CholeskyLowerAndInverse(matrix, lower, inverse)); + EXPECT_DOUBLE_EQ(lower(1U, 1U), 2); + EXPECT_DOUBLE_EQ(inverse(1U, 1U), .5); + EXPECT_THROW(static_cast(CholeskyLowerAndInverse( + Matrix{2U, 1U}, lower, inverse)), std::invalid_argument); + matrix(0U, 0U) = std::numeric_limits::quiet_NaN(); + EXPECT_THROW(static_cast(CholeskyLowerAndInverse(matrix, lower, inverse)), + std::invalid_argument); +} + +TEST(DenseBlasInternal, SingleEquationScopeRestoresTlsAndRespectsOverrides) { + 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 ScopedSingleEquationThreads scope{1U}; + EXPECT_EQ(mkl_get_max_threads(), 1); + } + EXPECT_EQ(mkl_get_max_threads(), default_threads); + EXPECT_EQ(mkl_set_num_threads_local(0), 0); + { + const ScopedSingleEquationThreads scope{2U}; + EXPECT_EQ(mkl_get_max_threads(), default_threads); + } + try { + const ScopedSingleEquationThreads scope{1U}; + throw std::runtime_error{"scope failure"}; + } catch (const std::runtime_error&) { + } + EXPECT_EQ(mkl_set_num_threads_local(0), 0); + static_cast(mkl_set_num_threads_local(3)); + { + const ScopedSingleEquationThreads scope{1U}; + EXPECT_EQ(mkl_get_max_threads(), 3); + } + EXPECT_EQ(mkl_get_max_threads(), 3); +} + } // namespace } // namespace fesa::dense_blas_internal diff --git a/tests/unit/math/high_low_arithmetic_original.h b/tests/unit/math/high_low_arithmetic_original.h new file mode 100644 index 0000000..072fe6b --- /dev/null +++ b/tests/unit/math/high_low_arithmetic_original.h @@ -0,0 +1,162 @@ +#ifndef FESA_TEST_HIGH_LOW_ARITHMETIC_ORIGINAL_H_ +#define FESA_TEST_HIGH_LOW_ARITHMETIC_ORIGINAL_H_ + +#include +#include +#include +#include +#if defined(_M_X64) || defined(__SSE2__) +#include +#endif + +namespace fesa::high_low_original { +/// @brief Checks the rounding and gradual-underflow premises without mutation. +inline bool ArithmeticEnvironmentSupported() { + if (std::fegetround() != FE_TONEAREST) return false; +#if defined(_M_X64) || defined(__SSE2__) + // MXCSR FTZ, DAZ and rounding-control bits must all be clear. + if ((_mm_getcsr() & 0xE040U) != 0U) return false; +#endif + return true; +} +/// @brief Two binary64 components of one physical value; raw parts are kept. +/// @note Arithmetic requires normalized finite parts and nearest rounding. +/// This is an expansion with bounded rounding, not an IEEE quad type. +struct HighLow { + double high{0}; + double low{0}; + HighLow() = default; + // Scalar lifting canonicalizes exact signed zero; raw pair inputs do not. + HighLow(double value) : high(value == 0 ? 0.0 : value) {} + HighLow(double high_part, double low_part) + : high(high_part), low(low_part) {} +}; + +/// @brief Gives an invalid arithmetic result for a checked caller to reject. +inline HighLow InvalidHighLow() { + return {std::numeric_limits::quiet_NaN(), 0}; +} +/// @brief Checks both physical components, before comparisons or publication. +inline bool IsFinite(HighLow value) { + return std::isfinite(value.high) && std::isfinite(value.low); +} +/// @brief Recovers the exact error of supported finite binary64 addition. +inline HighLow TwoSum(double left, double right) { + const double sum = left + right; + const double bridge = sum - left; + const double restored_left = sum - bridge; + const double error = (left - restored_left) + (right - bridge); + if (!std::isfinite(sum) || !std::isfinite(bridge) || + !std::isfinite(restored_left) || !std::isfinite(error)) + return InvalidHighLow(); + return sum == 0 && error == 0 ? HighLow{} : HighLow{sum, error}; +} +/// @brief Requires the raw pair to equal its nearest-rounding normalization. +inline bool IsNormalized(HighLow value) { + if (!IsFinite(value)) return false; + if (value.high == 0 && value.low == 0) + return !std::signbit(value.high) && !std::signbit(value.low); + const auto normalized = TwoSum(value.high, value.low); + return IsFinite(normalized) && normalized.high == value.high && + normalized.low == value.low; +} +/// @brief Rounds the full physical value to one finite output value. +inline double Rounded(HighLow value) { + return IsNormalized(value) ? value.high + value.low + : std::numeric_limits::quiet_NaN(); +} +/// @brief Adds normalized expansions; rounding of discarded tail is bounded. +inline HighLow Add(HighLow left, HighLow right) { + if (!IsNormalized(left) || !IsNormalized(right)) return InvalidHighLow(); + const auto head = TwoSum(left.high, right.high); + const auto tail = TwoSum(left.low, right.low); + const auto middle = TwoSum(head.low, tail.high); + const auto joined = TwoSum(head.high, middle.high); + const double remainder = (joined.low + middle.low) + tail.low; + if (!IsFinite(head) || !IsFinite(tail) || !IsFinite(middle) || + !IsFinite(joined) || !std::isfinite(remainder)) + return InvalidHighLow(); + return TwoSum(joined.high, remainder); +} +/// @brief Negates both components without collapsing the physical value. +inline HighLow Negate(HighLow value) { + if (!IsNormalized(value)) return InvalidHighLow(); + return value.high == 0 && value.low == 0 + ? HighLow{} + : HighLow{-value.high, -value.low}; +} +/// @brief Subtracts full expansions, including cancellations between highs. +inline HighLow Subtract(HighLow left, HighLow right) { + return Add(left, Negate(right)); +} +/// @brief Returns the exponent of the least nonzero binary64 significand bit. +/// @pre value is finite and nonzero; gradual underflow is enabled. +inline int LowestBitExponent(double value) { + int exponent = 0; + const double fraction = std::frexp(std::abs(value), &exponent); + auto significand = static_cast(std::ldexp(fraction, 53)); + int trailing = 0; + while ((significand & 1U) == 0U) { + significand >>= 1U; + ++trailing; + } + return exponent - 53 + trailing; +} +/// @brief Gives a product and its FMA residual when the dyadic tail fits. +/// @note Rejects underflow-risk exact tails, even if the rounded high is finite. +inline HighLow TwoProduct(double left, double right) { + if (!std::isfinite(left) || !std::isfinite(right)) return InvalidHighLow(); + if (left == 0 || right == 0) return {}; + if (LowestBitExponent(left) + LowestBitExponent(right) < -1074) + return InvalidHighLow(); + const double product = left * right; + if (!std::isfinite(product)) return InvalidHighLow(); + const double error = std::fma(left, right, -product); + if (!std::isfinite(error)) return InvalidHighLow(); + return TwoSum(product, error); +} +/// @brief Retains all four component products in deterministic order. +inline HighLow Multiply(HighLow left, HighLow right) { + if (!IsNormalized(left) || !IsNormalized(right)) return InvalidHighLow(); + auto product = TwoProduct(left.high, right.high); + product = Add(product, TwoProduct(left.high, right.low)); + product = Add(product, TwoProduct(left.low, right.high)); + return Add(product, TwoProduct(left.low, right.low)); +} +/// @brief Orders canonical values; 2 means unordered/unsupported input. +inline int Compare(HighLow left, HighLow right) { + if (!IsNormalized(left) || !IsNormalized(right)) return 2; + if (left.high < right.high) return -1; + if (left.high > right.high) return 1; + if (left.low < right.low) return -1; + if (left.low > right.low) return 1; + return 0; +} +/// @brief Takes the absolute value using the full physical sign. +inline HighLow Abs(HighLow value) { + const int sign = Compare(value, HighLow{}); + return sign == 2 ? InvalidHighLow() : (sign < 0 ? Negate(value) : value); +} +/// @brief Refines a binary64 quotient twice using the full residual. +/// @note Unsupported products/quotient underflow propagate an invalid result. +inline HighLow Divide(HighLow numerator, HighLow denominator) { + if (!IsNormalized(numerator) || !IsNormalized(denominator) || + Compare(denominator, HighLow{}) == 0) + return InvalidHighLow(); + if (Compare(numerator, HighLow{}) == 0) return {}; + const double first = numerator.high / denominator.high; + if (!std::isfinite(first) || first == 0) return InvalidHighLow(); + auto quotient = HighLow{first}; + for (int iteration = 0; iteration < 2; ++iteration) { + const auto residual = Subtract(numerator, Multiply(denominator, quotient)); + if (!IsNormalized(residual)) return InvalidHighLow(); + const double correction = residual.high / denominator.high; + if (!std::isfinite(correction) || (residual.high != 0 && correction == 0)) + return InvalidHighLow(); + quotient = Add(quotient, HighLow{correction}); + } + return quotient; +} +} // namespace fesa::high_low_original + +#endif // FESA_TEST_HIGH_LOW_ARITHMETIC_ORIGINAL_H_ diff --git a/tests/unit/math/high_low_arithmetic_test.cpp b/tests/unit/math/high_low_arithmetic_test.cpp index 092410d..03a94ec 100644 --- a/tests/unit/math/high_low_arithmetic_test.cpp +++ b/tests/unit/math/high_low_arithmetic_test.cpp @@ -3,7 +3,13 @@ #include #include +#include +#include #include +#include +#include + +#include "high_low_arithmetic_original.h" namespace { using namespace fesa::math_internal; @@ -125,4 +131,92 @@ TEST(HighLowArithmetic, RejectsChangedRoundingAndFlushToZeroEnvironment) { EXPECT_FALSE(ArithmeticEnvironmentSupported()); #endif } + +std::uint64_t Bits(const double value) { + std::uint64_t result = 0; + std::memcpy(&result, &value, sizeof(result)); + return result; +} + +double FromBits(const std::uint64_t bits) { + double result = 0; + std::memcpy(&result, &bits, sizeof(result)); + return result; +} + +void ExpectSamePair(const HighLow actual, + const fesa::high_low_original::HighLow expected) { + EXPECT_EQ(Bits(actual.high), Bits(expected.high)); + EXPECT_EQ(Bits(actual.low), Bits(expected.low)); +} + +TEST(HighLowArithmetic, BitExponentMatchesFrozenOracleAcrossEveryExponent) { + ASSERT_TRUE(ArithmeticEnvironmentSupported()); + // Every normal exponent, every possible trailing significand-bit count, + // both signs; then every subnormal single-bit and mixed-bit boundary. + for (std::uint64_t exponent = 1; exponent < 2047; ++exponent) { + for (unsigned trailing = 0; trailing <= 52U; ++trailing) { + const auto fraction = trailing == 52U ? 0ULL : 1ULL << trailing; + for (const auto sign : {0ULL, 1ULL << 63U}) { + const double value = FromBits(sign | (exponent << 52U) | fraction); + ASSERT_EQ(LowestBitExponent(value), + fesa::high_low_original::LowestBitExponent(value)); + } + } + } + for (unsigned bit = 0; bit < 52U; ++bit) { + for (const auto significand : {1ULL << bit, (1ULL << 52U) - (1ULL << bit)}) { + const double value = FromBits(significand); + EXPECT_EQ(LowestBitExponent(value), + fesa::high_low_original::LowestBitExponent(value)); + } + } +} + +TEST(HighLowArithmetic, OptimizedArithmeticMatchesFrozenOriginalBits) { + ASSERT_TRUE(ArithmeticEnvironmentSupported()); + std::vector values{0.0, -0.0, + std::numeric_limits::denorm_min(), + -std::numeric_limits::denorm_min(), + std::numeric_limits::min(), + std::nextafter(std::numeric_limits::min(), 0.0), + std::numeric_limits::max(), + -std::numeric_limits::max(), 1.0, -1.0, + std::nextafter(1.0, 2.0), std::ldexp(1.0, -53), + std::ldexp(1.0, -54), std::numeric_limits::infinity(), + -std::numeric_limits::infinity(), + std::numeric_limits::quiet_NaN()}; + std::mt19937_64 random{0x58a0fU}; + for (int sample = 0; sample < 4000; ++sample) { + values.push_back(FromBits(random())); + } + for (std::size_t index = 0; index < values.size(); ++index) { + const double high = values[index]; + const double low = values[(index * 13U + 7U) % values.size()]; + const fesa::high_low_original::HighLow original{high, low}; + EXPECT_EQ(IsNormalized({high, low}), + fesa::high_low_original::IsNormalized(original)); + // Normalize a raw pair with the frozen operation, never the new predicate. + const auto a = fesa::high_low_original::TwoSum(high, low); + const auto b = fesa::high_low_original::TwoSum( + values[(index + 1U) % values.size()], + values[(index * 3U + 2U) % values.size()]); + ExpectSamePair(Add({a.high, a.low}, {b.high, b.low}), + fesa::high_low_original::Add(a, b)); + ExpectSamePair(Multiply({a.high, a.low}, {b.high, b.low}), + fesa::high_low_original::Multiply(a, b)); + ExpectSamePair(Divide({a.high, a.low}, {b.high, b.low}), + fesa::high_low_original::Divide(a, b)); + ExpectSamePair(TwoProduct(high, low), + fesa::high_low_original::TwoProduct(high, low)); + } + for (double high : {1.0, std::nextafter(1.0, 2.0), -1.0, + std::numeric_limits::max()}) { + for (double low : {0.0, -0.0, std::ldexp(1.0, -53), + -std::ldexp(1.0, -53)}) { + EXPECT_EQ(IsNormalized({high, low}), + fesa::high_low_original::IsNormalized({high, low})); + } + } +} } // namespace diff --git a/tests/unit/math/matrix_test.cpp b/tests/unit/math/matrix_test.cpp index 575b848..b315b44 100644 --- a/tests/unit/math/matrix_test.cpp +++ b/tests/unit/math/matrix_test.cpp @@ -109,5 +109,23 @@ TEST(DenseMath, RowMajorMatrixMatchesKnownGemvGemm) { std::invalid_argument); } +TEST(DenseMath, MatrixProductsRejectNonfiniteOperandsAndPreserveValidReuse) { + Matrix left{2U, 2U, 1.0}; + Matrix right{2U, 1U, 2.0}; + Vector rhs{2U, 2.0}; + rhs[1U] = std::numeric_limits::infinity(); + EXPECT_THROW(static_cast(left.Multiply(rhs)), std::invalid_argument); + rhs[1U] = 2.0; + left(0U, 1U) = std::numeric_limits::quiet_NaN(); + EXPECT_THROW(static_cast(left.Multiply(rhs)), std::invalid_argument); + EXPECT_THROW(static_cast(left.Multiply(right)), std::invalid_argument); + left(0U, 1U) = 1.0; + right(1U, 0U) = std::numeric_limits::infinity(); + EXPECT_THROW(static_cast(left.Multiply(right)), std::invalid_argument); + right(1U, 0U) = 2.0; + EXPECT_DOUBLE_EQ(left.Multiply(rhs)[0U], 4.0); + EXPECT_DOUBLE_EQ(left.Multiply(right)(1U, 0U), 4.0); +} + } // namespace } // namespace fesa diff --git a/tests/unit/math/sparse_matrix_test.cpp b/tests/unit/math/sparse_matrix_test.cpp index 7754a2f..7185f91 100644 --- a/tests/unit/math/sparse_matrix_test.cpp +++ b/tests/unit/math/sparse_matrix_test.cpp @@ -3,6 +3,8 @@ #include #include +#include +#include #include #include #include @@ -131,4 +133,104 @@ TEST(SparseAssembly, PreservesExpectedStructuralZeros) { 4); } +TEST(SparseAssembly, MultipliesRectangularCopyMoveAndConcurrentConstCalls) { + auto assembled = SparseMatrix::FromCoo( + 3U, 5U, {{0U, 4U, 3.0, 0U, 0U}, {2U, 1U, -2.0, 1U, 0U}}, + {{0U, 1U, 1U, 2U}, {4U, 1U}}); + ASSERT_TRUE(assembled.HasValue()); + const auto original = assembled.Value(); + auto copied = original; + auto moved = std::move(copied); + Vector rhs{5U, 2.0}; + const auto evaluate = [&rhs](const SparseMatrix& matrix) { + for (int repeat = 0; repeat < 20; ++repeat) { + const auto product = matrix.Multiply(rhs); + if (product.Size() != 3U || product[0U] != 6.0 || + product[1U] != 0.0 || product[2U] != -4.0) return false; + } + return true; + }; + auto first = std::async(std::launch::async, [&] { return evaluate(moved); }); + auto second = std::async(std::launch::async, [&] { return evaluate(moved); }); + EXPECT_TRUE(first.get()); + EXPECT_TRUE(second.get()); + EXPECT_TRUE(evaluate(original)); + auto replacement = SparseMatrix::FromCoo( + 1U, 1U, {{0U, 0U, 7.0, 0U, 0U}}, {{0U, 1U}, {0U}}); + ASSERT_TRUE(replacement.HasValue()); + moved = replacement.Value(); + EXPECT_DOUBLE_EQ(moved.Multiply(Vector{1U, 3.0})[0U], 21.0); + EXPECT_TRUE(evaluate(original)); +} + +TEST(SparseAssembly, PreservesZeroShapesAndEmptyRows) { + auto no_rows = SparseMatrix::FromCoo(0U, 3U, {}, {{0U}, {}}); + auto no_columns = SparseMatrix::FromCoo(2U, 0U, {}, {{0U, 0U, 0U}, {}}); + auto no_entries = SparseMatrix::FromCoo(2U, 3U, {}, {{0U, 0U, 0U}, {}}); + ASSERT_TRUE(no_rows.HasValue()); + ASSERT_TRUE(no_columns.HasValue()); + ASSERT_TRUE(no_entries.HasValue()); + EXPECT_EQ(no_rows.Value().Multiply(Vector{3U}).Size(), 0U); + const auto empty_columns = no_columns.Value().Multiply(Vector{0U}); + const auto empty_entries = no_entries.Value().Multiply(Vector{3U}); + ASSERT_EQ(empty_columns.Size(), 2U); + ASSERT_EQ(empty_entries.Size(), 2U); + for (std::size_t row = 0; row < 2U; ++row) { + EXPECT_DOUBLE_EQ(empty_columns[row], 0.0); + EXPECT_DOUBLE_EQ(empty_entries[row], 0.0); + } +} + +TEST(SparseAssembly, RejectsNonfiniteRhsAndRecoversAfterFiniteOverflow) { + auto matrix = SparseMatrix::FromCoo( + 2U, 2U, {{0U, 0U, 2.0, 0U, 0U}}, {{0U, 1U, 2U}, {0U, 1U}}); + ASSERT_TRUE(matrix.HasValue()); + Vector rhs{2U, (std::numeric_limits::infinity)()}; + EXPECT_THROW(static_cast(matrix.Value().Multiply(rhs)), + std::invalid_argument); + rhs[1U] = 4.0; + rhs[0U] = (std::numeric_limits::max)(); + EXPECT_TRUE(std::isinf(matrix.Value().Multiply(rhs)[0U])); + rhs[0U] = (std::numeric_limits::quiet_NaN)(); + EXPECT_THROW(static_cast(matrix.Value().Multiply(rhs)), + std::invalid_argument); + EXPECT_THROW(static_cast(matrix.Value().Multiply(Vector{1U})), + std::invalid_argument); + rhs[0U] = 3.0; + rhs[1U] = 4.0; + const auto finite = matrix.Value().Multiply(rhs); + EXPECT_DOUBLE_EQ(finite[0U], 6.0); + EXPECT_DOUBLE_EQ(finite[1U], 0.0); +} + +TEST(SparseAssembly, MatchesIndependentLongDoubleResidualOnUnevenRows) { + const std::vector entries{1.25, -0.375, 4.5, -3.0, 0.125, 2.0}; + const SparsePattern pattern{{0U, 3U, 3U, 5U, 6U}, {0U, 2U, 4U, 1U, 3U, 2U}}; + std::vector contributions; + for (std::size_t row = 0; row < 4U; ++row) { + for (auto position = pattern.row_offsets[row]; + position < pattern.row_offsets[row + 1U]; ++position) { + contributions.push_back( + {row, pattern.column_indices[position], entries[position], row, + position}); + } + } + auto matrix = SparseMatrix::FromCoo(4U, 5U, contributions, pattern); + ASSERT_TRUE(matrix.HasValue()); + Vector rhs{5U}; + for (std::size_t column = 0; column < 5U; ++column) { + rhs[column] = 0.3 + static_cast(column) * 0.7; + } + const auto product = matrix.Value().Multiply(rhs); + for (std::size_t row = 0; row < 4U; ++row) { + long double expected = 0; + for (auto position = pattern.row_offsets[row]; + position < pattern.row_offsets[row + 1U]; ++position) { + expected += static_cast(entries[position]) * + rhs[pattern.column_indices[position]]; + } + EXPECT_NEAR(product[row], static_cast(expected), 2.0e-14); + } +} + } // namespace diff --git a/tests/unit/results/result_recovery_components_test.cpp b/tests/unit/results/result_recovery_components_test.cpp index 91bcc79..1f9a057 100644 --- a/tests/unit/results/result_recovery_components_test.cpp +++ b/tests/unit/results/result_recovery_components_test.cpp @@ -6,6 +6,7 @@ #include #include #include +#include #include #include #include @@ -353,6 +354,41 @@ TEST(ResultRecoveryComponents, } } +TEST(ResultRecoveryComponents, PartitionedScaleKeepsCancelledInternalTerms) { + auto definition = MakeShellDefinition(); + definition.node_sets = {{"Fixed", std::optional{"Shell-1"}, + {0U}, {definition.source_path, 60U}}}; + definition.steps[0].boundaries = { + {"Fixed", 1, 1, 1.0, {definition.source_path, 61U}}}; + const auto fixture = MakeShellFixture(std::move(definition)); + const auto count = fixture.dofs->FullDofCount(); + fesa::SparsePattern pattern; + pattern.row_offsets = {0U, 1U, 3U}; + pattern.row_offsets.resize(count + 1U, 3U); + pattern.column_indices = {0U, 0U, 1U}; + auto matrix = fesa::SparseMatrix::FromCoo( + count, count, + {{0U, 0U, 1.0, 0U, 0U}, {1U, 0U, -1.0, 0U, 1U}, + {1U, 1U, 1.0, 0U, 2U}}, pattern); + ASSERT_TRUE(matrix.HasValue()); + fesa::Vector displacement{count}; + displacement[0U] = 1.0; + displacement[1U] = 1.0; + fesa::Vector external{count}; + external[1U] = 1.0e-12; + const auto result = fesa::results_internal::GlobalEquilibriumRecovery( + *fixture.domain, *fixture.dofs, matrix.Value(), displacement, external); + ASSERT_TRUE(result.HasValue()); + EXPECT_DOUBLE_EQ(result.Value().internal_force[1U], 0.0); + EXPECT_DOUBLE_EQ(result.Value().normalized_free_residual, 1.0e-12); + displacement[1U] = std::numeric_limits::quiet_NaN(); + const auto rejected = fesa::results_internal::GlobalEquilibriumRecovery( + *fixture.domain, *fixture.dofs, matrix.Value(), displacement, external); + ASSERT_FALSE(rejected.HasValue()); + ExpectStatusCode(rejected.GetStatus(), "nonfinite-recovery-value"); + EXPECT_DOUBLE_EQ(result.Value().reaction[1U], -1.0e-12); +} + // C-MODULE-002 TEST(ResultRecoveryComponents, BeamResultRecoveryKeepsEndpointGaussAndStressIdentities) { diff --git a/tests/unit/solvers/linear/mkl_pardiso_solver_test.cpp b/tests/unit/solvers/linear/mkl_pardiso_solver_test.cpp index 1a6953e..4758e70 100644 --- a/tests/unit/solvers/linear/mkl_pardiso_solver_test.cpp +++ b/tests/unit/solvers/linear/mkl_pardiso_solver_test.cpp @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -11,6 +12,7 @@ #include "fesa/fem/dof_manager.h" #include "fesa/math/sparse_matrix.h" +#include "math/dense_blas_internal.h" namespace { @@ -123,6 +125,42 @@ TEST(MklPardisoSolver, ReusesOneFactorizationForRepeatedRhs) { EXPECT_LE(NormalizedResidual(matrix, second, rhs_second), 1.0e-10); } +TEST(MklPardisoSolver, SingleEquationPolicyPreservesLifecycleAndConcurrentTls) { + const auto worker = [] { + 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(); + fesa::MklPardisoSolver solver; + for (int repeat = 0; repeat < 10; ++repeat) { + if (!solver.Factorize(MakeDenseCsr(1U, {2.0})).IsOk()) return false; + auto solution = MakeVector({9.0}); + if (!solver.Solve(MakeVector({4.0}), solution).IsOk() || + solution[0U] != 2.0) return false; + if (!solver.Factorize(MakeDenseCsr(1U, {4.0})).IsOk()) return false; + if (!solver.Solve(MakeVector({4.0}), solution).IsOk() || + solution[0U] != 1.0) return false; + auto invalid_rhs = MakeVector({std::numeric_limits::quiet_NaN()}); + if (solver.Solve(invalid_rhs, solution).IsOk() || solution[0U] != 1.0) + return false; + if (!solver.Solve(MakeVector({8.0}), solution).IsOk() || + solution[0U] != 2.0) return false; + if (solver.Factorize(MakeDenseCsr(1U, {-1.0})).IsOk()) return false; + if (solver.Solve(MakeVector({8.0}), solution).IsOk() || + solution[0U] != 2.0) return false; + if (mkl_get_max_threads() != default_threads || + mkl_set_num_threads_local(0) != 0) return false; + } + return true; + }; + auto first = std::async(std::launch::async, worker); + auto second = std::async(std::launch::async, worker); + EXPECT_TRUE(first.get()); + EXPECT_TRUE(second.get()); +} + TEST(MklPardisoSolver, StillRejectsNonsymmetricPublicCsr) { const auto nonsymmetric = MakeDenseCsr(2U, {2.0, 1.0, 0.0, 3.0}); fesa::MklPardisoSolver solver; @@ -155,6 +193,32 @@ TEST(MklPardisoSolver, RefactorizesWithoutLeakingState) { 1.0e-10); } +TEST(MklPardisoSolver, SingleEquationScopeTransitionsBetweenRetainedSizes) { + 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(); + fesa::MklPardisoSolver solver; + ASSERT_TRUE(solver.Factorize(MakeDenseCsr(1U, {2.0})).IsOk()); + auto single = MakeVector({0.0}); + ASSERT_TRUE(solver.Solve(MakeVector({4.0}), single).IsOk()); + EXPECT_DOUBLE_EQ(single[0U], 2.0); + ASSERT_TRUE(solver.Factorize(MakeDenseCsr(2U, {2, 0, 0, 3})).IsOk()); + fesa::Vector pair{2U}; + ASSERT_TRUE(solver.Solve(MakeVector({4, 6}), pair).IsOk()); + EXPECT_NEAR(pair[0U], 2, 1e-14); + EXPECT_NEAR(pair[1U], 2, 1e-14); + EXPECT_EQ(mkl_get_max_threads(), default_threads); + EXPECT_EQ(mkl_set_num_threads_local(0), 0); + ASSERT_TRUE(solver.Factorize(MakeDenseCsr(1U, {4.0})).IsOk()); + ASSERT_TRUE(solver.Solve(MakeVector({8.0}), single).IsOk()); + EXPECT_DOUBLE_EQ(single[0U], 2.0); + EXPECT_EQ(mkl_get_max_threads(), default_threads); + EXPECT_EQ(mkl_set_num_threads_local(0), 0); +} + TEST(MklPardisoSolver, ClassifiesSingularIndefiniteAndNonfiniteFailures) { fesa::MklPardisoSolver solver; diff --git a/tests/unit/solvers/linear/tangent_solver_test.cpp b/tests/unit/solvers/linear/tangent_solver_test.cpp index 0dab893..bb7096f 100644 --- a/tests/unit/solvers/linear/tangent_solver_test.cpp +++ b/tests/unit/solvers/linear/tangent_solver_test.cpp @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -13,6 +14,7 @@ #include "fesa/math/sparse_matrix.h" #include "fesa/math/vector.h" #include "fesa/solvers/linear/tangent_solver.h" +#include "math/dense_blas_internal.h" namespace { @@ -129,4 +131,104 @@ TEST(TangentSolver, RejectsSolveBeforeSuccessfulFactorization) { EXPECT_DOUBLE_EQ(correction[1U], -9.0); } +TEST(TangentSolver, ReusesRhsButFreshValuesAndFailureCannotReuseStaleState) { + fesa::MklPardisoTangentSolver solver; + ASSERT_TRUE(solver.FactorizeFresh(MakeDenseCsr(2U, {2.0, 1.0, 0.0, 3.0})) + .IsOk()); + fesa::Vector correction{2U}; + for (double scale : {1.0, -2.0, 0.5}) { + ASSERT_TRUE(solver.SolveCorrection(MakeVector({4.0 * scale, 6.0 * scale}), + correction).IsOk()); + EXPECT_NEAR(correction[0U], scale, 1.0e-14); + EXPECT_NEAR(correction[1U], 2.0 * scale, 1.0e-14); + } + const auto bad_rhs = MakeVector({std::numeric_limits::quiet_NaN(), 0}); + correction = MakeVector({17, -4}); + EXPECT_FALSE(solver.SolveCorrection(bad_rhs, correction).IsOk()); + EXPECT_DOUBLE_EQ(correction[0U], 17); + EXPECT_DOUBLE_EQ(correction[1U], -4); + ASSERT_TRUE(solver.SolveCorrection(MakeVector({4, 6}), correction).IsOk()); + ASSERT_TRUE(solver.FactorizeFresh(MakeDenseCsr(2U, {4, 0, 1, 2})).IsOk()); + ASSERT_TRUE(solver.SolveCorrection(MakeVector({4, 5}), correction).IsOk()); + EXPECT_NEAR(correction[0U], 1, 1.0e-14); + EXPECT_NEAR(correction[1U], 2, 1.0e-14); + auto diagonal = fesa::SparseMatrix::FromCoo( + 2U, 2U, {{0U, 0U, 5, 0U, 0U}, {1U, 1U, 7, 1U, 0U}}, + {{0U, 1U, 2U}, {0U, 1U}}); + ASSERT_TRUE(diagonal.HasValue()); + ASSERT_TRUE(solver.FactorizeFresh(diagonal.Value()).IsOk()); + ASSERT_TRUE(solver.SolveCorrection(MakeVector({10, -7}), correction).IsOk()); + EXPECT_NEAR(correction[0U], 2, 1.0e-14); + EXPECT_NEAR(correction[1U], -1, 1.0e-14); + EXPECT_FALSE(solver.FactorizeFresh(MakeDenseCsr(2U, {0, 0, 0, 0})).IsOk()); + correction = MakeVector({19, -8}); + const auto stale = solver.SolveCorrection(MakeVector({10, -7}), correction); + ASSERT_FALSE(stale.IsOk()); + EXPECT_EQ(stale.Diagnostics().front().code, "solver-not-factorized"); + EXPECT_DOUBLE_EQ(correction[0U], 19); + EXPECT_DOUBLE_EQ(correction[1U], -8); + ASSERT_TRUE(solver.FactorizeFresh(diagonal.Value()).IsOk()); + EXPECT_TRUE(solver.SolveCorrection(MakeVector({10, -7}), correction).IsOk()); +} + +TEST(TangentSolver, SingleEquationScopeRestoresTlsAcrossFreshFailureAndReuse) { + const auto worker = [] { + 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(); + fesa::MklPardisoTangentSolver solver; + for (int repeat = 0; repeat < 10; ++repeat) { + auto correction = MakeVector({17.0}); + if (!solver.FactorizeFresh(MakeDenseCsr(1U, {2.0})).IsOk()) return false; + if (!solver.SolveCorrection(MakeVector({4.0}), correction).IsOk() || + correction[0U] != 2.0) return false; + if (!solver.FactorizeFresh(MakeDenseCsr(1U, {4.0})).IsOk()) return false; + if (!solver.SolveCorrection(MakeVector({4.0}), correction).IsOk() || + correction[0U] != 1.0) return false; + if (solver.SolveCorrection(MakeVector({ + std::numeric_limits::infinity()}), correction).IsOk() || + correction[0U] != 1.0) return false; + if (solver.FactorizeFresh(MakeDenseCsr(1U, {0.0})).IsOk()) return false; + if (solver.SolveCorrection(MakeVector({4.0}), correction).IsOk() || + correction[0U] != 1.0) return false; + if (mkl_get_max_threads() != default_threads || + mkl_set_num_threads_local(0) != 0) return false; + } + return true; + }; + auto first = std::async(std::launch::async, worker); + auto second = std::async(std::launch::async, worker); + EXPECT_TRUE(first.get()); + EXPECT_TRUE(second.get()); +} + +TEST(TangentSolver, SingleEquationScopeTransitionsBetweenRetainedSizes) { + 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(); + fesa::MklPardisoTangentSolver solver; + ASSERT_TRUE(solver.FactorizeFresh(MakeDenseCsr(1U, {2.0})).IsOk()); + auto single = MakeVector({0.0}); + ASSERT_TRUE(solver.SolveCorrection(MakeVector({4.0}), single).IsOk()); + EXPECT_DOUBLE_EQ(single[0U], 2.0); + ASSERT_TRUE(solver.FactorizeFresh(MakeDenseCsr(2U, {2, 1, 0, 3})).IsOk()); + fesa::Vector pair{2U}; + ASSERT_TRUE(solver.SolveCorrection(MakeVector({6, 6}), pair).IsOk()); + EXPECT_NEAR(pair[0U], 2, 1e-14); + EXPECT_NEAR(pair[1U], 2, 1e-14); + EXPECT_EQ(mkl_get_max_threads(), default_threads); + EXPECT_EQ(mkl_set_num_threads_local(0), 0); + ASSERT_TRUE(solver.FactorizeFresh(MakeDenseCsr(1U, {4.0})).IsOk()); + ASSERT_TRUE(solver.SolveCorrection(MakeVector({8.0}), single).IsOk()); + EXPECT_DOUBLE_EQ(single[0U], 2.0); + EXPECT_EQ(mkl_get_max_threads(), default_threads); + EXPECT_EQ(mkl_set_num_threads_local(0), 0); +} + } // namespace