From ae7e259e17ba2cedc19cda06a8d16b3f199620d7 Mon Sep 17 00:00:00 2001 From: Pasukhin Dmitry Date: Tue, 3 Mar 2026 14:18:18 +0000 Subject: [PATCH] Foundation Classes - align modern Math* APIs with legacy math_* behavior (#1134) MathLin: - Return full matrix solutions for multi-RHS APIs. - Add LinearMultipleResult for matrix RHS solve results. MathSys: - Fix Newton2D/3D/4D tiny-step exit logic: re-check residual at updated point and return OK when converged. MathUtils / MathInteg: - Add modern Gauss points/weights implementation in MathUtils_Gauss.cxx. - Keep legacy-table parity for orders 1..61 and compute fallback for higher orders. - Make GaussAdaptive use IntegConfig InitialOrder/MaxOrder with bounds validation. - Propagate ordered Gauss points/weights retrieval failures in set/multiple integration. - Extend BracketMinimum API with bounded/options-based behavior. Tests: - Extend MathLin, MathSys and MathInteg tests for new behavior and regressions. - Add MathUtils bracketing tests. - Add MathLin_EigenSearch parity test coverage against legacy solver. Documentation: - Update MathLin/MathInteg/MathUtils READMEs to match current APIs and behavior. --- .../TKMath/GTests/FILES.cmake | 2 + .../GTests/MathInteg_Comparison_Test.cxx | 13 + .../TKMath/GTests/MathInteg_Test.cxx | 52 +++- .../GTests/MathLin_EigenSearch_Test.cxx | 284 ++++++++++++++++++ .../TKMath/GTests/MathLin_Test.cxx | 115 +++++++ .../TKMath/GTests/MathSys_Newton2D_Test.cxx | 90 +++++- .../TKMath/GTests/MathSys_Newton3D_Test.cxx | 75 +++++ .../TKMath/GTests/MathSys_Newton4D_Test.cxx | 77 +++++ .../TKMath/GTests/MathUtils_Bracket_Test.cxx | 94 ++++++ .../TKMath/MathInteg/MathInteg_Gauss.hxx | 49 ++- .../TKMath/MathInteg/MathInteg_Multiple.hxx | 6 +- .../TKMath/MathInteg/MathInteg_Set.hxx | 6 +- .../TKMath/MathInteg/README.md | 4 +- .../TKMath/MathLin/MathLin_Gauss.hxx | 18 +- .../TKMath/MathLin/MathLin_Householder.hxx | 61 ++-- .../TKMath/MathLin/README.md | 7 +- .../TKMath/MathSys/MathSys_Newton2D.hxx | 42 ++- .../TKMath/MathSys/MathSys_Newton3D.hxx | 13 +- .../TKMath/MathSys/MathSys_Newton4D.hxx | 13 +- .../TKMath/MathUtils/FILES.cmake | 1 + .../TKMath/MathUtils/MathUtils_Bracket.hxx | 212 ++++++++++++- .../TKMath/MathUtils/MathUtils_Gauss.cxx | 115 +++++++ .../TKMath/MathUtils/MathUtils_Gauss.hxx | 269 +---------------- .../MathUtils_GaussKronrodWeights.cxx | 7 +- .../TKMath/MathUtils/MathUtils_Types.hxx | 15 + .../TKMath/MathUtils/README.md | 5 +- 26 files changed, 1278 insertions(+), 367 deletions(-) create mode 100644 src/FoundationClasses/TKMath/GTests/MathLin_EigenSearch_Test.cxx create mode 100644 src/FoundationClasses/TKMath/GTests/MathUtils_Bracket_Test.cxx create mode 100644 src/FoundationClasses/TKMath/MathUtils/MathUtils_Gauss.cxx diff --git a/src/FoundationClasses/TKMath/GTests/FILES.cmake b/src/FoundationClasses/TKMath/GTests/FILES.cmake index e645860cb1..1dc191190c 100644 --- a/src/FoundationClasses/TKMath/GTests/FILES.cmake +++ b/src/FoundationClasses/TKMath/GTests/FILES.cmake @@ -86,6 +86,7 @@ set(OCCT_TKMath_GTests_FILES math_Uzawa_Test.cxx math_Vector_Test.cxx # MathUtils tests + MathUtils_Bracket_Test.cxx MathUtils_Functor_Test.cxx # MathPoly tests MathPoly_Test.cxx @@ -93,6 +94,7 @@ set(OCCT_TKMath_GTests_FILES MathPoly_Laguerre_Test.cxx # MathLin tests MathLin_Test.cxx + MathLin_EigenSearch_Test.cxx MathLin_Comparison_Test.cxx # MathOpt tests MathOpt_1D_Test.cxx diff --git a/src/FoundationClasses/TKMath/GTests/MathInteg_Comparison_Test.cxx b/src/FoundationClasses/TKMath/GTests/MathInteg_Comparison_Test.cxx index f8ef2cc64e..efb4c5c582 100644 --- a/src/FoundationClasses/TKMath/GTests/MathInteg_Comparison_Test.cxx +++ b/src/FoundationClasses/TKMath/GTests/MathInteg_Comparison_Test.cxx @@ -443,6 +443,19 @@ TEST(MathInteg_ComparisonTest, Order21_Comparison) EXPECT_NEAR(anOldInteg.Value(), *aNewResult.Value, THE_TOLERANCE); } +TEST(MathInteg_ComparisonTest, Order41_Comparison) +{ + SinFuncOld anOldFunc; + SinFuncNew aNewFunc; + + math_GaussSingleIntegration anOldInteg(anOldFunc, 0.0, THE_PI, 41); + MathInteg::IntegResult aNewResult = MathInteg::Gauss(aNewFunc, 0.0, THE_PI, 41); + + ASSERT_TRUE(anOldInteg.IsDone()); + ASSERT_TRUE(aNewResult.IsDone()); + EXPECT_NEAR(anOldInteg.Value(), *aNewResult.Value, THE_TOLERANCE); +} + // ============================================================================ // Higher order accuracy comparison // ============================================================================ diff --git a/src/FoundationClasses/TKMath/GTests/MathInteg_Test.cxx b/src/FoundationClasses/TKMath/GTests/MathInteg_Test.cxx index 182f104ecd..1d1c1c968a 100644 --- a/src/FoundationClasses/TKMath/GTests/MathInteg_Test.cxx +++ b/src/FoundationClasses/TKMath/GTests/MathInteg_Test.cxx @@ -275,11 +275,28 @@ TEST(MathInteg_GaussTest, Order21) EXPECT_NEAR(*aResult.Value, 2.0, THE_TOLERANCE); } -TEST(MathInteg_GaussTest, InvalidOrder) +TEST(MathInteg_GaussTest, Order9) { - SinFunc aFunc; - // Order 9 is not supported (supported orders: 3, 4, 5, 6, 7, 8, 10, 15, 21, 31) + SinFunc aFunc; MathInteg::IntegResult aResult = MathInteg::Gauss(aFunc, 0.0, THE_PI, 9); + ASSERT_TRUE(aResult.IsDone()); + EXPECT_EQ(aResult.NbPoints, 9); + EXPECT_NEAR(*aResult.Value, 2.0, THE_TOLERANCE); +} + +TEST(MathInteg_GaussTest, Order61) +{ + SinFunc aFunc; + MathInteg::IntegResult aResult = MathInteg::Gauss(aFunc, 0.0, THE_PI, 61); + ASSERT_TRUE(aResult.IsDone()); + EXPECT_EQ(aResult.NbPoints, 61); + EXPECT_NEAR(*aResult.Value, 2.0, THE_TOLERANCE); +} + +TEST(MathInteg_GaussTest, InvalidOrderNonPositive) +{ + SinFunc aFunc; + MathInteg::IntegResult aResult = MathInteg::Gauss(aFunc, 0.0, THE_PI, 0); EXPECT_EQ(aResult.Status, MathInteg::Status::InvalidInput); } @@ -353,8 +370,28 @@ TEST(MathInteg_GaussAdaptiveTest, ProvidesErrorEstimate) MathInteg::IntegResult aResult = MathInteg::GaussAdaptive(aFunc, 0.0, THE_PI, aConfig); ASSERT_TRUE(aResult.IsDone()); - EXPECT_GT(aResult.AbsoluteError, 0.0); - EXPECT_LT(aResult.AbsoluteError, 1.0e-6); + ASSERT_TRUE(aResult.AbsoluteError.has_value()); + ASSERT_TRUE(aResult.RelativeError.has_value()); + EXPECT_TRUE(std::isfinite(*aResult.AbsoluteError)); + EXPECT_TRUE(std::isfinite(*aResult.RelativeError)); + EXPECT_GE(*aResult.AbsoluteError, 0.0); + EXPECT_GE(*aResult.RelativeError, 0.0); + EXPECT_LT(*aResult.AbsoluteError, 1.0e-6); +} + +TEST(MathInteg_GaussAdaptiveTest, UsesConfiguredOrders) +{ + QuadraticFunc aFunc; + MathInteg::IntegConfig aConfig; + aConfig.InitialOrder = 9; + aConfig.MaxOrder = 18; + aConfig.Tolerance = 1.0e-12; + aConfig.MaxIterations = 2; + + MathInteg::IntegResult aResult = MathInteg::GaussAdaptive(aFunc, 0.0, 1.0, aConfig); + ASSERT_TRUE(aResult.IsDone()); + EXPECT_EQ(aResult.NbPoints, 18u); + EXPECT_NEAR(*aResult.Value, 1.0 / 3.0, THE_TOLERANCE); } // ============================================================================ @@ -479,9 +516,8 @@ TEST(MathInteg_BoolConversionTest, SuccessfulResultIsTrue) TEST(MathInteg_BoolConversionTest, InvalidInputIsFalse) { - SinFunc aFunc; - // Order 9 is not supported - MathInteg::IntegResult aResult = MathInteg::Gauss(aFunc, 0.0, THE_PI, 9); + SinFunc aFunc; + MathInteg::IntegResult aResult = MathInteg::Gauss(aFunc, 0.0, THE_PI, 0); EXPECT_FALSE(static_cast(aResult)); } diff --git a/src/FoundationClasses/TKMath/GTests/MathLin_EigenSearch_Test.cxx b/src/FoundationClasses/TKMath/GTests/MathLin_EigenSearch_Test.cxx new file mode 100644 index 0000000000..ff0a11f147 --- /dev/null +++ b/src/FoundationClasses/TKMath/GTests/MathLin_EigenSearch_Test.cxx @@ -0,0 +1,284 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +namespace +{ +constexpr double THE_EIGEN_TOL = 1.0e-10; +constexpr double THE_RESIDUAL_TOL = 1.0e-10; +constexpr double THE_ORTHOGONAL_TOL = 1.0e-10; +constexpr double THE_NORMALIZED_TOL = 1.0e-10; +constexpr int THE_RANDOM_NB_CASES = 120; + +math_Matrix BuildSymmetricTridiagonal(const math_Vector& theDiag, const math_Vector& theSubdiag) +{ + const int aN = theDiag.Length(); + math_Matrix aM(1, aN, 1, aN, 0.0); + for (int i = 1; i <= aN; ++i) + { + aM(i, i) = theDiag(theDiag.Lower() + i - 1); + } + for (int i = 2; i <= aN; ++i) + { + const double aE = theSubdiag(theSubdiag.Lower() + i - 1); + aM(i, i - 1) = aE; + aM(i - 1, i) = aE; + } + return aM; +} + +void BuildLegacyArrays(const math_Vector& theDiag, + const math_Vector& theSubdiag, + NCollection_Array1& theDiagLegacy, + NCollection_Array1& theSubdiagLegacy) +{ + const int aN = theDiag.Length(); + for (int i = 1; i <= aN; ++i) + { + theDiagLegacy(i) = theDiag(theDiag.Lower() + i - 1); + theSubdiagLegacy(i) = theSubdiag(theSubdiag.Lower() + i - 1); + } +} + +std::vector SortedEigenValuesFromLegacy(const math_EigenValuesSearcher& theLegacy) +{ + const int aN = theLegacy.Dimension(); + std::vector aVals; + aVals.reserve(static_cast(aN)); + for (int i = 1; i <= aN; ++i) + { + aVals.push_back(theLegacy.EigenValue(i)); + } + std::sort(aVals.begin(), aVals.end()); + return aVals; +} + +std::vector SortedEigenValuesFromModern(const MathLin::EigenResult& theModern) +{ + std::vector aVals; + if (!theModern.EigenValues.has_value()) + { + return aVals; + } + + const math_Vector& aEig = *theModern.EigenValues; + aVals.reserve(static_cast(aEig.Length())); + for (int i = aEig.Lower(); i <= aEig.Upper(); ++i) + { + aVals.push_back(aEig(i)); + } + std::sort(aVals.begin(), aVals.end()); + return aVals; +} + +double VectorNorm2(const math_Vector& theVec) +{ + double aNorm2 = 0.0; + for (int i = theVec.Lower(); i <= theVec.Upper(); ++i) + { + aNorm2 += theVec(i) * theVec(i); + } + return std::sqrt(aNorm2); +} + +double PairResidualInfinity(const math_Matrix& theMatrix, + double theLambda, + const math_Vector& theVector) +{ + const int aN = theMatrix.RowNumber(); + double aMax = 0.0; + for (int i = 1; i <= aN; ++i) + { + double aAx = 0.0; + for (int j = 1; j <= aN; ++j) + { + aAx += theMatrix(i, j) * theVector(j); + } + const double aRes = std::abs(aAx - theLambda * theVector(i)); + if (aRes > aMax) + { + aMax = aRes; + } + } + return aMax; +} + +double DotProduct(const math_Vector& theV1, const math_Vector& theV2) +{ + const int aN = theV1.Length(); + double aDot = 0.0; + for (int i = 1; i <= aN; ++i) + { + aDot += theV1(i) * theV2(i); + } + return aDot; +} + +} // namespace + +TEST(MathLin_EigenSearch_Test, BasicParityWithLegacy_3x3) +{ + math_Vector aDiag(1, 3); + math_Vector aSubdiag(1, 3); + aDiag(1) = 4.0; + aDiag(2) = 4.0; + aDiag(3) = 4.0; + aSubdiag(1) = 0.0; + aSubdiag(2) = 1.0; + aSubdiag(3) = 1.0; + + NCollection_Array1 aDiagLegacy(1, 3); + NCollection_Array1 aSubdiagLegacy(1, 3); + BuildLegacyArrays(aDiag, aSubdiag, aDiagLegacy, aSubdiagLegacy); + + math_EigenValuesSearcher aLegacy(aDiagLegacy, aSubdiagLegacy); + MathLin::EigenResult aModern = MathLin::EigenTridiagonal(aDiag, aSubdiag); + + ASSERT_TRUE(aLegacy.IsDone()); + ASSERT_TRUE(aModern.IsDone()); + ASSERT_TRUE(aModern.EigenValues.has_value()); + ASSERT_TRUE(aModern.EigenVectors.has_value()); + + const std::vector aLegacySorted = SortedEigenValuesFromLegacy(aLegacy); + const std::vector aModernSorted = SortedEigenValuesFromModern(aModern); + ASSERT_EQ(aLegacySorted.size(), aModernSorted.size()); + for (size_t i = 0; i < aLegacySorted.size(); ++i) + { + EXPECT_NEAR(aLegacySorted[i], aModernSorted[i], THE_EIGEN_TOL); + } + + const math_Matrix aA = BuildSymmetricTridiagonal(aDiag, aSubdiag); + const math_Vector& aEigVals = *aModern.EigenValues; + for (int i = 1; i <= 3; ++i) + { + const math_Vector aVec = MathLin::GetEigenVector(aModern, i); + EXPECT_NEAR(VectorNorm2(aVec), 1.0, THE_NORMALIZED_TOL); + EXPECT_NEAR(PairResidualInfinity(aA, aEigVals(i), aVec), 0.0, THE_RESIDUAL_TOL); + } +} + +TEST(MathLin_EigenSearch_Test, HandlesNonOneLowerBounds) +{ + math_Vector aDiag(-2, 2); + math_Vector aSubdiag(-2, 2); + aDiag(-2) = 1.0; + aDiag(-1) = 3.0; + aDiag(0) = -2.0; + aDiag(1) = 5.0; + aDiag(2) = 4.0; + aSubdiag(-2) = 0.0; + aSubdiag(-1) = 0.3; + aSubdiag(0) = -0.2; + aSubdiag(1) = 0.7; + aSubdiag(2) = -1.1; + + NCollection_Array1 aDiagLegacy(-2, 2); + NCollection_Array1 aSubdiagLegacy(-2, 2); + for (int i = -2; i <= 2; ++i) + { + aDiagLegacy(i) = aDiag(i); + aSubdiagLegacy(i) = aSubdiag(i); + } + + math_EigenValuesSearcher aLegacy(aDiagLegacy, aSubdiagLegacy); + MathLin::EigenResult aModern = MathLin::EigenTridiagonal(aDiag, aSubdiag); + + ASSERT_EQ(aLegacy.IsDone(), aModern.IsDone()); + ASSERT_TRUE(aLegacy.IsDone()); + + const std::vector aLegacySorted = SortedEigenValuesFromLegacy(aLegacy); + const std::vector aModernSorted = SortedEigenValuesFromModern(aModern); + ASSERT_EQ(aLegacySorted.size(), aModernSorted.size()); + for (size_t i = 0; i < aLegacySorted.size(); ++i) + { + EXPECT_NEAR(aLegacySorted[i], aModernSorted[i], THE_EIGEN_TOL); + } +} + +TEST(MathLin_EigenSearch_Test, RandomParityAndOrthogonality) +{ + std::mt19937 aGen(123456u); + std::uniform_int_distribution aDimDist(2, 32); + std::uniform_real_distribution aValDist(-100.0, 100.0); + + for (int aCase = 0; aCase < THE_RANDOM_NB_CASES; ++aCase) + { + const int aN = aDimDist(aGen); + + math_Vector aDiag(1, aN); + math_Vector aSubdiag(1, aN); + for (int i = 1; i <= aN; ++i) + { + aDiag(i) = aValDist(aGen); + aSubdiag(i) = (i == 1) ? 0.0 : aValDist(aGen); + } + + NCollection_Array1 aDiagLegacy(1, aN); + NCollection_Array1 aSubdiagLegacy(1, aN); + BuildLegacyArrays(aDiag, aSubdiag, aDiagLegacy, aSubdiagLegacy); + + math_EigenValuesSearcher aLegacy(aDiagLegacy, aSubdiagLegacy); + MathLin::EigenResult aModern = MathLin::EigenTridiagonal(aDiag, aSubdiag); + + ASSERT_EQ(aLegacy.IsDone(), aModern.IsDone()) << "case=" << aCase; + if (!aLegacy.IsDone()) + { + continue; + } + + ASSERT_TRUE(aModern.EigenValues.has_value()) << "case=" << aCase; + ASSERT_TRUE(aModern.EigenVectors.has_value()) << "case=" << aCase; + + const std::vector aLegacySorted = SortedEigenValuesFromLegacy(aLegacy); + const std::vector aModernSorted = SortedEigenValuesFromModern(aModern); + ASSERT_EQ(aLegacySorted.size(), aModernSorted.size()) << "case=" << aCase; + for (size_t i = 0; i < aLegacySorted.size(); ++i) + { + EXPECT_NEAR(aLegacySorted[i], aModernSorted[i], THE_EIGEN_TOL) << "case=" << aCase; + } + + const math_Matrix aA = BuildSymmetricTridiagonal(aDiag, aSubdiag); + const math_Vector& aEigVals = *aModern.EigenValues; + + for (int i = 1; i <= aN; ++i) + { + const math_Vector aVec = MathLin::GetEigenVector(aModern, i); + EXPECT_NEAR(VectorNorm2(aVec), 1.0, THE_NORMALIZED_TOL) << "case=" << aCase; + EXPECT_NEAR(PairResidualInfinity(aA, aEigVals(i), aVec), 0.0, THE_RESIDUAL_TOL) + << "case=" << aCase; + } + + for (int i = 1; i <= aN; ++i) + { + const math_Vector aVecI = MathLin::GetEigenVector(aModern, i); + for (int j = i + 1; j <= aN; ++j) + { + const math_Vector aVecJ = MathLin::GetEigenVector(aModern, j); + EXPECT_NEAR(DotProduct(aVecI, aVecJ), 0.0, THE_ORTHOGONAL_TOL) << "case=" << aCase; + } + } + } +} diff --git a/src/FoundationClasses/TKMath/GTests/MathLin_Test.cxx b/src/FoundationClasses/TKMath/GTests/MathLin_Test.cxx index e51dc5adb3..877d18b30b 100644 --- a/src/FoundationClasses/TKMath/GTests/MathLin_Test.cxx +++ b/src/FoundationClasses/TKMath/GTests/MathLin_Test.cxx @@ -262,6 +262,121 @@ TEST(MathLin_SVD_Test, ConditionNumber) EXPECT_GT(aCondH, 100.0); // Hilbert matrices are ill-conditioned } +// ============================================================================ +// Multi-RHS linear solve tests +// ============================================================================ + +TEST(MathLin_Gauss_Test, SolveMultiple_ReturnsFullMatrix) +{ + math_Matrix aA(1, 3, 1, 3); + aA(1, 1) = 4.0; + aA(1, 2) = 1.0; + aA(1, 3) = 2.0; + aA(2, 1) = 0.0; + aA(2, 2) = 3.0; + aA(2, 3) = 1.0; + aA(3, 1) = 2.0; + aA(3, 2) = 1.0; + aA(3, 3) = 5.0; + + math_Matrix aXExpected(1, 3, 1, 2); + aXExpected(1, 1) = 1.0; + aXExpected(2, 1) = 2.0; + aXExpected(3, 1) = -1.0; + aXExpected(1, 2) = -2.0; + aXExpected(2, 2) = 0.5; + aXExpected(3, 2) = 3.0; + + const math_Matrix aB = MatMul(aA, aXExpected); + + auto aResult = MathLin::SolveMultiple(aA, aB); + ASSERT_TRUE(aResult.IsDone()); + ASSERT_TRUE(aResult.Solutions.has_value()); + + const math_Matrix& aX = *aResult.Solutions; + for (int i = 1; i <= 3; ++i) + { + for (int j = 1; j <= 2; ++j) + { + EXPECT_NEAR(aX(i, j), aXExpected(i, j), THE_TOLERANCE); + } + } + + const math_Matrix aCheckB = MatMul(aA, aX); + for (int i = 1; i <= 3; ++i) + { + for (int j = 1; j <= 2; ++j) + { + EXPECT_NEAR(aCheckB(i, j), aB(i, j), THE_TOLERANCE); + } + } +} + +TEST(MathLin_Gauss_Test, SolveMultiple_DimensionMismatch) +{ + math_Matrix aA(1, 3, 1, 3, 0.0); + for (int i = 1; i <= 3; ++i) + { + aA(i, i) = 1.0; + } + + math_Matrix aBWrong(1, 2, 1, 1, 0.0); + auto aResult = MathLin::SolveMultiple(aA, aBWrong); + EXPECT_EQ(aResult.Status, MathUtils::Status::InvalidInput); +} + +TEST(MathLin_Householder_Test, SolveQRMultiple_ReturnsFullMatrix) +{ + math_Matrix aA(1, 3, 1, 2); + aA(1, 1) = 1.0; + aA(1, 2) = 2.0; + aA(2, 1) = 3.0; + aA(2, 2) = 1.0; + aA(3, 1) = -1.0; + aA(3, 2) = 1.0; + + math_Matrix aXExpected(1, 2, 1, 2); + aXExpected(1, 1) = 2.0; + aXExpected(2, 1) = -1.0; + aXExpected(1, 2) = -0.5; + aXExpected(2, 2) = 3.0; + + const math_Matrix aB = MatMul(aA, aXExpected); + + auto aResult = MathLin::SolveQRMultiple(aA, aB); + ASSERT_TRUE(aResult.IsDone()); + ASSERT_TRUE(aResult.Solutions.has_value()); + + const math_Matrix& aX = *aResult.Solutions; + for (int i = 1; i <= 2; ++i) + { + for (int j = 1; j <= 2; ++j) + { + EXPECT_NEAR(aX(i, j), aXExpected(i, j), THE_TOLERANCE); + } + } + + const math_Matrix aCheckB = MatMul(aA, aX); + for (int i = 1; i <= 3; ++i) + { + for (int j = 1; j <= 2; ++j) + { + EXPECT_NEAR(aCheckB(i, j), aB(i, j), THE_TOLERANCE); + } + } +} + +TEST(MathLin_Householder_Test, SolveQRMultiple_DimensionMismatch) +{ + math_Matrix aA(1, 3, 1, 2, 0.0); + aA(1, 1) = 1.0; + aA(2, 2) = 1.0; + + math_Matrix aBWrong(1, 2, 1, 2, 0.0); + auto aResult = MathLin::SolveQRMultiple(aA, aBWrong); + EXPECT_EQ(aResult.Status, MathUtils::Status::InvalidInput); +} + // ============================================================================ // Householder QR tests // ============================================================================ diff --git a/src/FoundationClasses/TKMath/GTests/MathSys_Newton2D_Test.cxx b/src/FoundationClasses/TKMath/GTests/MathSys_Newton2D_Test.cxx index ebe5f21f3d..f5d0a19441 100644 --- a/src/FoundationClasses/TKMath/GTests/MathSys_Newton2D_Test.cxx +++ b/src/FoundationClasses/TKMath/GTests/MathSys_Newton2D_Test.cxx @@ -24,21 +24,44 @@ namespace class QuadraticFunc { public: - bool ValueAndJacobian(double theU, - double theV, - double& theF1, - double& theF2, - double& theJ11, - double& theJ12, - double& theJ21, - double& theJ22) const + bool operator()(double theU, double theV, double theF[2], double theJ[2][2]) const { - theF1 = 2.0 * theU; - theF2 = 2.0 * theV; - theJ11 = 2.0; - theJ12 = 0.0; - theJ21 = 0.0; - theJ22 = 2.0; + theF[0] = 2.0 * theU; + theF[1] = 2.0 * theV; + theJ[0][0] = 2.0; + theJ[0][1] = 0.0; + theJ[1][0] = 0.0; + theJ[1][1] = 2.0; + return true; + } +}; + +class GenericLinearExactStep +{ +public: + bool operator()(double theU, double theV, double theF[2], double theJ[2][2]) const + { + theF[0] = theU - 1.0; + theF[1] = theV - 2.0; + theJ[0][0] = 1.0; + theJ[0][1] = 0.0; + theJ[1][0] = 0.0; + theJ[1][1] = 1.0; + return true; + } +}; + +class GenericHugeJacobianConstantResidual +{ +public: + bool operator()(double /*theU*/, double /*theV*/, double theF[2], double theJ[2][2]) const + { + theF[0] = 1.0; + theF[1] = 1.0; + theJ[0][0] = 1.0e20; + theJ[0][1] = 0.0; + theJ[1][0] = 0.0; + theJ[1][1] = 1.0e20; return true; } }; @@ -162,6 +185,45 @@ TEST(MathSys_Newton2DTest, Solve2D_Quadratic_Converges) EXPECT_LT(aResult.ResidualNorm, 1.0e-10); } +TEST(MathSys_Newton2DTest, Solve2D_SmallStepAtRoot_ReturnsOK) +{ + GenericLinearExactStep aFunc; + + MathSys::NewtonBoundsN<2> aBounds; + aBounds.Min = {-10.0, -10.0}; + aBounds.Max = {10.0, 10.0}; + + MathSys::NewtonOptions aOptions; + aOptions.FTolerance = 1.0e-12; + aOptions.XTolerance = 100.0; + aOptions.MaxIterations = 5; + + const MathSys::NewtonResultN<2> aResult = MathSys::Solve2D(aFunc, {0.0, 0.0}, aBounds, aOptions); + EXPECT_TRUE(aResult.IsDone()); + EXPECT_EQ(aResult.Status, MathUtils::Status::OK); + EXPECT_NEAR(aResult.X[0], 1.0, 1.0e-14); + EXPECT_NEAR(aResult.X[1], 2.0, 1.0e-14); +} + +TEST(MathSys_Newton2DTest, Solve2D_TinyStepLargeResidual_ReturnsMaxIterations) +{ + GenericHugeJacobianConstantResidual aFunc; + + MathSys::NewtonBoundsN<2> aBounds; + aBounds.Min = {-1.0, -1.0}; + aBounds.Max = {1.0, 1.0}; + + MathSys::NewtonOptions aOptions; + aOptions.FTolerance = 1.0e-8; + aOptions.XTolerance = 1.0e-16; + aOptions.MaxIterations = 10; + + const MathSys::NewtonResultN<2> aResult = MathSys::Solve2D(aFunc, {0.0, 0.0}, aBounds, aOptions); + EXPECT_FALSE(aResult.IsDone()); + EXPECT_EQ(aResult.Status, MathUtils::Status::MaxIterations); + EXPECT_GT(aResult.ResidualNorm, 1.0e-2); +} + TEST(MathSys_Newton2DTest, Solve2DSymmetric_Target_Converges) { SymmetricDistFunc aFunc(3.5, 7.2); diff --git a/src/FoundationClasses/TKMath/GTests/MathSys_Newton3D_Test.cxx b/src/FoundationClasses/TKMath/GTests/MathSys_Newton3D_Test.cxx index 6a9d7651f7..f47601366b 100644 --- a/src/FoundationClasses/TKMath/GTests/MathSys_Newton3D_Test.cxx +++ b/src/FoundationClasses/TKMath/GTests/MathSys_Newton3D_Test.cxx @@ -94,6 +94,81 @@ TEST_F(MathSys_Newton3DTest, Solve3D_NonlinearSystem) EXPECT_NEAR(aResult.X[2], 1.0, 1.0e-5); } +TEST_F(MathSys_Newton3DTest, Solve3D_SmallStepAtRoot_ReturnsOK) +{ + auto aFunc = + [](double theX1, double theX2, double theX3, double theF[3], double theJ[3][3]) -> bool { + theF[0] = theX1 - 1.0; + theF[1] = theX2 - 2.0; + theF[2] = theX3 - 3.0; + + theJ[0][0] = 1.0; + theJ[0][1] = 0.0; + theJ[0][2] = 0.0; + theJ[1][0] = 0.0; + theJ[1][1] = 1.0; + theJ[1][2] = 0.0; + theJ[2][0] = 0.0; + theJ[2][1] = 0.0; + theJ[2][2] = 1.0; + return true; + }; + + MathSys::NewtonBoundsN<3> aBounds; + aBounds.HasBounds = false; + + MathSys::NewtonOptions aOptions; + aOptions.FTolerance = 1.0e-12; + aOptions.XTolerance = 100.0; + aOptions.MaxIterations = 5; + aOptions.MaxStepRatio = 100.0; + + const MathSys::NewtonResultN<3> aResult = + MathSys::Solve3D(aFunc, {0.0, 0.0, 0.0}, aBounds, aOptions); + EXPECT_TRUE(aResult.IsDone()); + EXPECT_EQ(aResult.Status, MathUtils::Status::OK); + EXPECT_NEAR(aResult.X[0], 1.0, 1.0e-14); + EXPECT_NEAR(aResult.X[1], 2.0, 1.0e-14); + EXPECT_NEAR(aResult.X[2], 3.0, 1.0e-14); +} + +TEST_F(MathSys_Newton3DTest, Solve3D_TinyStepLargeResidual_ReturnsMaxIterations) +{ + auto aFunc = [](double /*theX1*/, + double /*theX2*/, + double /*theX3*/, + double theF[3], + double theJ[3][3]) -> bool { + theF[0] = 1.0; + theF[1] = 1.0; + theF[2] = 1.0; + theJ[0][0] = 1.0e20; + theJ[0][1] = 0.0; + theJ[0][2] = 0.0; + theJ[1][0] = 0.0; + theJ[1][1] = 1.0e20; + theJ[1][2] = 0.0; + theJ[2][0] = 0.0; + theJ[2][1] = 0.0; + theJ[2][2] = 1.0e20; + return true; + }; + + MathSys::NewtonBoundsN<3> aBounds; + aBounds.HasBounds = false; + + MathSys::NewtonOptions aOptions; + aOptions.FTolerance = 1.0e-8; + aOptions.XTolerance = 1.0e-16; + aOptions.MaxIterations = 10; + + const MathSys::NewtonResultN<3> aResult = + MathSys::Solve3D(aFunc, {0.0, 0.0, 0.0}, aBounds, aOptions); + EXPECT_FALSE(aResult.IsDone()); + EXPECT_EQ(aResult.Status, MathUtils::Status::MaxIterations); + EXPECT_GT(aResult.ResidualNorm, 1.0e-2); +} + TEST_F(MathSys_Newton3DTest, Solve3D_Bounded) { auto aFunc = diff --git a/src/FoundationClasses/TKMath/GTests/MathSys_Newton4D_Test.cxx b/src/FoundationClasses/TKMath/GTests/MathSys_Newton4D_Test.cxx index e6016e8f95..5071d9d685 100644 --- a/src/FoundationClasses/TKMath/GTests/MathSys_Newton4D_Test.cxx +++ b/src/FoundationClasses/TKMath/GTests/MathSys_Newton4D_Test.cxx @@ -108,6 +108,83 @@ TEST_F(MathSys_Newton4DTest, Solve4D_Bounded) EXPECT_NEAR(aResult.X[3], 4.0, 1.0e-12); } +TEST_F(MathSys_Newton4DTest, Solve4D_SmallStepAtRoot_ReturnsOK) +{ + auto aFunc = + [](double theX1, double theX2, double theX3, double theX4, double theF[4], double theJ[4][4]) + -> bool { + theF[0] = theX1 - 1.0; + theF[1] = theX2 - 2.0; + theF[2] = theX3 - 3.0; + theF[3] = theX4 - 4.0; + + for (int r = 0; r < 4; ++r) + { + for (int c = 0; c < 4; ++c) + { + theJ[r][c] = (r == c) ? 1.0 : 0.0; + } + } + return true; + }; + + MathSys::NewtonBoundsN<4> aBounds; + aBounds.HasBounds = false; + + MathSys::NewtonOptions aOptions; + aOptions.FTolerance = 1.0e-12; + aOptions.XTolerance = 100.0; + aOptions.MaxIterations = 5; + aOptions.MaxStepRatio = 100.0; + + const MathSys::NewtonResultN<4> aResult = + MathSys::Solve4D(aFunc, {0.0, 0.0, 0.0, 0.0}, aBounds, aOptions); + EXPECT_TRUE(aResult.IsDone()); + EXPECT_EQ(aResult.Status, MathUtils::Status::OK); + EXPECT_NEAR(aResult.X[0], 1.0, 1.0e-14); + EXPECT_NEAR(aResult.X[1], 2.0, 1.0e-14); + EXPECT_NEAR(aResult.X[2], 3.0, 1.0e-14); + EXPECT_NEAR(aResult.X[3], 4.0, 1.0e-14); +} + +TEST_F(MathSys_Newton4DTest, Solve4D_TinyStepLargeResidual_ReturnsMaxIterations) +{ + auto aFunc = [](double /*theX1*/, + double /*theX2*/, + double /*theX3*/, + double /*theX4*/, + double theF[4], + double theJ[4][4]) -> bool { + theF[0] = 1.0; + theF[1] = 1.0; + theF[2] = 1.0; + theF[3] = 1.0; + + for (int r = 0; r < 4; ++r) + { + for (int c = 0; c < 4; ++c) + { + theJ[r][c] = (r == c) ? 1.0e20 : 0.0; + } + } + return true; + }; + + MathSys::NewtonBoundsN<4> aBounds; + aBounds.HasBounds = false; + + MathSys::NewtonOptions aOptions; + aOptions.FTolerance = 1.0e-8; + aOptions.XTolerance = 1.0e-16; + aOptions.MaxIterations = 10; + + const MathSys::NewtonResultN<4> aResult = + MathSys::Solve4D(aFunc, {0.0, 0.0, 0.0, 0.0}, aBounds, aOptions); + EXPECT_FALSE(aResult.IsDone()); + EXPECT_EQ(aResult.Status, MathUtils::Status::MaxIterations); + EXPECT_GT(aResult.ResidualNorm, 1.0e-2); +} + TEST_F(MathSys_Newton4DTest, Solve4D_InvalidInput) { auto aFunc = [](double /*theX1*/, diff --git a/src/FoundationClasses/TKMath/GTests/MathUtils_Bracket_Test.cxx b/src/FoundationClasses/TKMath/GTests/MathUtils_Bracket_Test.cxx new file mode 100644 index 0000000000..c3ecd49236 --- /dev/null +++ b/src/FoundationClasses/TKMath/GTests/MathUtils_Bracket_Test.cxx @@ -0,0 +1,94 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include + +#include + +#include + +namespace +{ +class QuadraticMinimum +{ +public: + bool Value(double theX, double& theF) const + { + theF = (theX - 2.0) * (theX - 2.0); + return true; + } +}; + +class QuadraticWithPrecomputedEndpoints +{ +public: + bool Value(double theX, double& theF) const + { + if (std::abs(theX - 0.0) < 1.0e-15 || std::abs(theX - 1.0) < 1.0e-15) + { + return false; + } + theF = (theX - 2.0) * (theX - 2.0); + return true; + } +}; +} // namespace + +TEST(MathUtils_BracketTest, BracketMinimum_WithLimits_Succeeds) +{ + QuadraticMinimum aFunc; + + MathUtils::MinBracketOptions anOptions; + anOptions.MaxIterations = 50; + anOptions.UseLimits = true; + anOptions.LeftLimit = 0.0; + anOptions.RightLimit = 5.0; + + const MathUtils::MinBracketResult aResult = MathUtils::BracketMinimum(aFunc, 0.0, 1.0, anOptions); + ASSERT_TRUE(aResult.IsValid); + EXPECT_GE(aResult.A, anOptions.LeftLimit); + EXPECT_LE(aResult.C, anOptions.RightLimit); + EXPECT_LT(aResult.Fb, aResult.Fa); + EXPECT_LT(aResult.Fb, aResult.Fc); +} + +TEST(MathUtils_BracketTest, BracketMinimum_WithRestrictiveLimits_Fails) +{ + QuadraticMinimum aFunc; + + MathUtils::MinBracketOptions anOptions; + anOptions.MaxIterations = 50; + anOptions.UseLimits = true; + anOptions.LeftLimit = 0.0; + anOptions.RightLimit = 1.0; + + const MathUtils::MinBracketResult aResult = MathUtils::BracketMinimum(aFunc, 0.0, 0.5, anOptions); + EXPECT_FALSE(aResult.IsValid); +} + +TEST(MathUtils_BracketTest, BracketMinimum_UsesPrecomputedEndpointValues) +{ + QuadraticWithPrecomputedEndpoints aFunc; + + MathUtils::MinBracketOptions anOptions; + anOptions.MaxIterations = 50; + anOptions.HasFA = true; + anOptions.HasFB = true; + anOptions.FA = 4.0; + anOptions.FB = 1.0; + + const MathUtils::MinBracketResult aResult = MathUtils::BracketMinimum(aFunc, 0.0, 1.0, anOptions); + ASSERT_TRUE(aResult.IsValid); + EXPECT_LT(aResult.Fb, aResult.Fa); + EXPECT_LT(aResult.Fb, aResult.Fc); +} diff --git a/src/FoundationClasses/TKMath/MathInteg/MathInteg_Gauss.hxx b/src/FoundationClasses/TKMath/MathInteg/MathInteg_Gauss.hxx index d96f649df9..947d73ed08 100644 --- a/src/FoundationClasses/TKMath/MathInteg/MathInteg_Gauss.hxx +++ b/src/FoundationClasses/TKMath/MathInteg/MathInteg_Gauss.hxx @@ -19,6 +19,7 @@ #include #include +#include #include //! Numerical integration algorithms. @@ -40,16 +41,21 @@ using namespace MathUtils; //! @param theFunc function to integrate //! @param theLower lower integration bound //! @param theUpper upper integration bound -//! @param theNbPoints number of quadrature points (3, 4, 5, 6, 7, 8, 10, 15, 21, or 31) +//! @param theNbPoints number of quadrature points (>= 1) //! @return result containing integral value template IntegResult Gauss(Function& theFunc, double theLower, double theUpper, int theNbPoints = 15) { IntegResult aResult; + if (theNbPoints < 1) + { + aResult.Status = Status::InvalidInput; + return aResult; + } // Get quadrature points and weights - const double* aPoints = nullptr; - const double* aWeights = nullptr; + math_Vector aPoints(1, theNbPoints); + math_Vector aWeights(1, theNbPoints); if (!MathUtils::GetGaussPointsAndWeights(theNbPoints, aPoints, aWeights)) { @@ -62,16 +68,16 @@ IntegResult Gauss(Function& theFunc, double theLower, double theUpper, int theNb const double aMid = 0.5 * (theUpper + theLower); double aSum = 0.0; - for (int i = 0; i < theNbPoints; ++i) + for (int i = 1; i <= theNbPoints; ++i) { - const double aX = aMid + aHalfLen * aPoints[i]; + const double aX = aMid + aHalfLen * aPoints(i); double aF = 0.0; if (!theFunc.Value(aX, aF)) { aResult.Status = Status::NumericalError; return aResult; } - aSum += aWeights[i] * aF; + aSum += aWeights(i) * aF; } aResult.Status = Status::OK; @@ -103,14 +109,35 @@ IntegResult GaussAdaptive(Function& theFunc, { IntegResult aResult; + if (theConfig.InitialOrder < 1 || theConfig.MaxOrder < theConfig.InitialOrder + || theConfig.MaxOrder > 61 || theConfig.MaxIterations < 1) + { + aResult.Status = Status::InvalidInput; + return aResult; + } + + int aCoarseOrder = theConfig.InitialOrder; + int aFineOrder = std::min(theConfig.MaxOrder, std::min(61, 2 * aCoarseOrder)); + if (aFineOrder == aCoarseOrder) + { + if (aCoarseOrder > 1) + { + aCoarseOrder -= 1; + } + else if (theConfig.MaxOrder > 1) + { + aFineOrder = 2; + } + } + // Compute with coarse and fine grids - IntegResult aCoarse = Gauss(theFunc, theLower, theUpper, 7); + IntegResult aCoarse = Gauss(theFunc, theLower, theUpper, aCoarseOrder); if (!aCoarse.IsDone()) { return aCoarse; } - IntegResult aFine = Gauss(theFunc, theLower, theUpper, 15); + IntegResult aFine = Gauss(theFunc, theLower, theUpper, aFineOrder); if (!aFine.IsDone()) { return aFine; @@ -126,7 +153,7 @@ IntegResult GaussAdaptive(Function& theFunc, aResult.Value = *aFine.Value; aResult.AbsoluteError = aError; aResult.RelativeError = aError / aScale; - aResult.NbPoints = 15; + aResult.NbPoints = static_cast(aFineOrder); aResult.NbIterations = 1; return aResult; } @@ -138,7 +165,7 @@ IntegResult GaussAdaptive(Function& theFunc, aResult.Value = *aFine.Value; aResult.AbsoluteError = aError; aResult.RelativeError = aError / aScale; - aResult.NbPoints = 15; + aResult.NbPoints = static_cast(aFineOrder); aResult.NbIterations = 1; return aResult; } @@ -183,7 +210,7 @@ IntegResult GaussAdaptive(Function& theFunc, //! @param theLower lower integration bound //! @param theUpper upper integration bound //! @param theNbIntervals number of subintervals -//! @param theNbPoints Gauss points per interval (3, 4, 5, 6, 7, 8, 10, 15, 21, or 31) +//! @param theNbPoints Gauss points per interval (>= 1) //! @return result containing integral value template IntegResult GaussComposite(Function& theFunc, diff --git a/src/FoundationClasses/TKMath/MathInteg/MathInteg_Multiple.hxx b/src/FoundationClasses/TKMath/MathInteg/MathInteg_Multiple.hxx index a8909e4788..f03b3fb965 100644 --- a/src/FoundationClasses/TKMath/MathInteg/MathInteg_Multiple.hxx +++ b/src/FoundationClasses/TKMath/MathInteg/MathInteg_Multiple.hxx @@ -108,7 +108,11 @@ IntegResult GaussMultiple(Func& theFunc, math_Vector aGP(1, aOrd(i)); math_Vector aGW(1, aOrd(i)); - GetOrderedGaussPointsAndWeights(aOrd(i), aGP, aGW); + if (!GetOrderedGaussPointsAndWeights(aOrd(i), aGP, aGW)) + { + aResult.Status = Status::InvalidInput; + return aResult; + } for (int k = 0; k < aOrd(i); ++k) { diff --git a/src/FoundationClasses/TKMath/MathInteg/MathInteg_Set.hxx b/src/FoundationClasses/TKMath/MathInteg/MathInteg_Set.hxx index 62a346dca9..ed5ecb4078 100644 --- a/src/FoundationClasses/TKMath/MathInteg/MathInteg_Set.hxx +++ b/src/FoundationClasses/TKMath/MathInteg/MathInteg_Set.hxx @@ -77,7 +77,11 @@ SetResult GaussSet(Func& theFunc, double theLower, double theUpper, int theOrder // Get Gauss points and weights math_Vector aGP(1, aOrder); math_Vector aGW(1, aOrder); - GetOrderedGaussPointsAndWeights(aOrder, aGP, aGW); + if (!GetOrderedGaussPointsAndWeights(aOrder, aGP, aGW)) + { + aResult.Status = Status::InvalidInput; + return aResult; + } math_Vector aPoints(0, aOrder - 1); math_Vector aWeights(0, aOrder - 1); diff --git a/src/FoundationClasses/TKMath/MathInteg/README.md b/src/FoundationClasses/TKMath/MathInteg/README.md index 402220e82e..1b94fa17e3 100644 --- a/src/FoundationClasses/TKMath/MathInteg/README.md +++ b/src/FoundationClasses/TKMath/MathInteg/README.md @@ -10,8 +10,8 @@ The MathInteg package provides a collection of numerical integration methods for ### MathInteg_Gauss.hxx Gauss-Legendre quadrature methods: -- `Gauss` - Fixed-order Gauss-Legendre integration -- `GaussAdaptive` - Adaptive subdivision with error control +- `Gauss` - Fixed-order Gauss-Legendre integration (orders >= 1) +- `GaussAdaptive` - Adaptive subdivision with error control using `IntegConfig.InitialOrder/MaxOrder` - `GaussComposite` - Composite rule over multiple subintervals ### MathInteg_Kronrod.hxx diff --git a/src/FoundationClasses/TKMath/MathLin/MathLin_Gauss.hxx b/src/FoundationClasses/TKMath/MathLin/MathLin_Gauss.hxx index 85314b7f27..6f915f1c0a 100644 --- a/src/FoundationClasses/TKMath/MathLin/MathLin_Gauss.hxx +++ b/src/FoundationClasses/TKMath/MathLin/MathLin_Gauss.hxx @@ -236,11 +236,11 @@ inline LinearResult Solve(const math_Matrix& theA, //! @param theB right-hand side matrix //! @param theMinPivot minimum pivot value //! @return result containing solution matrix -inline LinearResult SolveMultiple(const math_Matrix& theA, - const math_Matrix& theB, - double theMinPivot = 1.0e-20) +inline LinearMultipleResult SolveMultiple(const math_Matrix& theA, + const math_Matrix& theB, + double theMinPivot = 1.0e-20) { - LinearResult aResult; + LinearMultipleResult aResult; // Perform LU decomposition LUResult aLURes = LU(theA, theMinPivot); @@ -264,7 +264,7 @@ inline LinearResult SolveMultiple(const math_Matrix& theA, const math_IntegerVector& aPivot = *aLURes.Pivot; // Solve for each column of B - math_Matrix aX(theB.LowerRow(), theB.UpperRow(), theB.LowerCol(), theB.UpperCol()); + math_Matrix aX(aRowLower, aRowUpper, theB.LowerCol(), theB.UpperCol()); for (int col = theB.LowerCol(); col <= theB.UpperCol(); ++col) { @@ -317,13 +317,7 @@ inline LinearResult SolveMultiple(const math_Matrix& theA, aResult.Status = Status::OK; aResult.Determinant = aLURes.Determinant; - // Store first column as solution vector for compatibility - math_Vector aSol(aRowLower, aRowUpper); - for (int i = aRowLower; i <= aRowUpper; ++i) - { - aSol(i) = aX(i, theB.LowerCol()); - } - aResult.Solution = aSol; + aResult.Solutions = aX; return aResult; } diff --git a/src/FoundationClasses/TKMath/MathLin/MathLin_Householder.hxx b/src/FoundationClasses/TKMath/MathLin/MathLin_Householder.hxx index c17950e9c6..c60b7aab9c 100644 --- a/src/FoundationClasses/TKMath/MathLin/MathLin_Householder.hxx +++ b/src/FoundationClasses/TKMath/MathLin/MathLin_Householder.hxx @@ -281,12 +281,14 @@ inline LinearResult SolveQR(const math_Matrix& theA, //! @param theB right-hand side matrix (m x p) //! @param theTolerance for singularity detection //! @return result containing solution matrix (n x p) -inline LinearResult SolveQRMultiple(const math_Matrix& theA, - const math_Matrix& theB, - double theTolerance = 1.0e-20) +inline LinearMultipleResult SolveQRMultiple(const math_Matrix& theA, + const math_Matrix& theB, + double theTolerance = 1.0e-20) { - LinearResult aResult; + LinearMultipleResult aResult; + const int aRowLower = theA.LowerRow(); + const int aRowUpper = theA.UpperRow(); const int aColLower = theA.LowerCol(); const int aColUpper = theA.UpperCol(); @@ -305,36 +307,51 @@ inline LinearResult SolveQRMultiple(const math_Matrix& theA, return aResult; } + const math_Matrix& aQ = *aQR.Q; + const math_Matrix& aR = *aQR.R; + // Solve for each column of B - math_Vector aFirstSol(aColLower, aColUpper); - bool aFirstDone = false; + math_Matrix aX(aColLower, aColUpper, theB.LowerCol(), theB.UpperCol(), 0.0); for (int j = theB.LowerCol(); j <= theB.UpperCol(); ++j) { - // Extract column j of B - math_Vector aBj(theB.LowerRow(), theB.UpperRow()); - for (int i = theB.LowerRow(); i <= theB.UpperRow(); ++i) + // Compute c = Q^T * b_j + math_Vector aC(aRowLower, aRowUpper, 0.0); + for (int i = aRowLower; i <= aRowUpper; ++i) { - aBj(i) = theB(i, j); + double aSum = 0.0; + for (int k = aRowLower; k <= aRowUpper; ++k) + { + const int aBRow = theB.LowerRow() + (k - aRowLower); + aSum += aQ(k, i) * theB(aBRow, j); + } + aC(i) = aSum; } - // Solve with QR (we should refactor to avoid re-decomposition) - LinearResult aColResult = SolveQR(theA, aBj, theTolerance); - if (!aColResult.IsDone()) + // Back substitution: R[1:n,1:n] * x_j = c[1:n] + for (int i = aColUpper; i >= aColLower; --i) { - aResult.Status = aColResult.Status; - return aResult; - } + const int aIOffset = i - aColLower; + const int aRRow = aRowLower + aIOffset; + double aDiag = aR(aRRow, i); - if (!aFirstDone) - { - aFirstSol = *aColResult.Solution; - aFirstDone = true; + if (std::abs(aDiag) < theTolerance) + { + aResult.Status = Status::Singular; + return aResult; + } + + double aSum = aC(aRRow); + for (int k = i + 1; k <= aColUpper; ++k) + { + aSum -= aR(aRRow, k) * aX(k, j); + } + aX(i, j) = aSum / aDiag; } } - aResult.Solution = aFirstSol; - aResult.Status = Status::OK; + aResult.Solutions = aX; + aResult.Status = Status::OK; return aResult; } diff --git a/src/FoundationClasses/TKMath/MathLin/README.md b/src/FoundationClasses/TKMath/MathLin/README.md index 8aaf4a2dba..eae79b6b2a 100644 --- a/src/FoundationClasses/TKMath/MathLin/README.md +++ b/src/FoundationClasses/TKMath/MathLin/README.md @@ -8,7 +8,7 @@ The MathLin package provides modern C++ implementations of linear algebra solver LU decomposition with partial pivoting for solving general linear systems Ax = b. - `LU()` - LU decomposition of a square matrix - `Solve()` - Solve linear system using LU decomposition -- `SolveMultiple()` - Solve multiple right-hand sides +- `SolveMultiple()` - Solve multiple right-hand sides (returns full solution matrix) - `Determinant()` - Compute matrix determinant - `Invert()` - Compute matrix inverse @@ -30,7 +30,7 @@ Singular Value Decomposition for general and ill-conditioned matrices. QR decomposition using Householder reflections. - `QR()` - QR decomposition of a matrix - `SolveQR()` - Solve overdetermined system (least squares) -- `SolveQRMultiple()` - Solve multiple right-hand sides +- `SolveQRMultiple()` - Solve multiple right-hand sides (returns full solution matrix) ### MathLin_Jacobi.hxx Jacobi method for eigenvalue decomposition of symmetric matrices. @@ -74,6 +74,9 @@ if (result.IsDone()) } ``` +For matrix right-hand sides (`SolveMultiple`, `SolveQRMultiple`), APIs return +`LinearMultipleResult` with `Solutions` (`math_Matrix`) instead of a single vector. + ## Dependencies The MathLin package depends on: diff --git a/src/FoundationClasses/TKMath/MathSys/MathSys_Newton2D.hxx b/src/FoundationClasses/TKMath/MathSys/MathSys_Newton2D.hxx index 3588c250cd..1f5238d9aa 100644 --- a/src/FoundationClasses/TKMath/MathSys/MathSys_Newton2D.hxx +++ b/src/FoundationClasses/TKMath/MathSys/MathSys_Newton2D.hxx @@ -165,10 +165,8 @@ inline bool SolveSymmetric2x2SVD(double theJ11, //! Solve a general 2x2 nonlinear system by Newton iteration. //! Function contract: -//! bool ValueAndJacobian(double u, double v, -//! double& f1, double& f2, -//! double& j11, double& j12, -//! double& j21, double& j22) const; +//! bool operator()(double u, double v, +//! double f[2], double j[2][2]) const; template NewtonResultN<2> Solve2D(const Function& theFunc, const std::array& theX0, @@ -193,14 +191,15 @@ NewtonResultN<2> Solve2D(const Function& theFunc, { aRes.NbIterations = static_cast(anIter + 1); - double aF1, aF2, aJ11, aJ12, aJ21, aJ22; - if (!theFunc.ValueAndJacobian(aRes.X[0], aRes.X[1], aF1, aF2, aJ11, aJ12, aJ21, aJ22)) + double aF[2]; + double aJ[2][2]; + if (!theFunc(aRes.X[0], aRes.X[1], aF, aJ)) { aRes.Status = MathUtils::Status::NumericalError; return aRes; } - const double aFNormSq = aF1 * aF1 + aF2 * aF2; + const double aFNormSq = aF[0] * aF[0] + aF[1] * aF[1]; aRes.ResidualNorm = std::sqrt(aFNormSq); if (aFNormSq <= aTolSq) @@ -212,11 +211,11 @@ NewtonResultN<2> Solve2D(const Function& theFunc, double aDU = 0.0; double aDV = 0.0; - const double aDet = aJ11 * aJ22 - aJ12 * aJ21; + const double aDet = aJ[0][0] * aJ[1][1] - aJ[0][1] * aJ[1][0]; if (std::abs(aDet) < detail::THE_SINGULAR_DET_TOL) { - const double aGradU = aJ11 * aF1 + aJ21 * aF2; - const double aGradV = aJ12 * aF1 + aJ22 * aF2; + const double aGradU = aJ[0][0] * aF[0] + aJ[1][0] * aF[1]; + const double aGradV = aJ[0][1] * aF[0] + aJ[1][1] * aF[1]; const double aGradSq = aGradU * aGradU + aGradV * aGradV; if (aGradSq < detail::THE_CRITICAL_GRAD_SQ) { @@ -231,8 +230,8 @@ NewtonResultN<2> Solve2D(const Function& theFunc, else { const double aInvDet = 1.0 / aDet; - aDU = (-aF1 * aJ22 + aF2 * aJ12) * aInvDet; - aDV = (-aF2 * aJ11 + aF1 * aJ21) * aInvDet; + aDU = (-aF[0] * aJ[1][1] + aF[1] * aJ[0][1]) * aInvDet; + aDV = (-aF[1] * aJ[0][0] + aF[0] * aJ[1][0]) * aInvDet; } const double aStepNormSq = aDU * aDU + aDV * aDV; @@ -254,19 +253,30 @@ NewtonResultN<2> Solve2D(const Function& theFunc, const double aScaleRef = std::max(1.0, std::max(std::abs(aRes.X[0]), std::abs(aRes.X[1]))); if (aRes.StepNorm <= theOptions.XTolerance * aScaleRef) { - aRes.Status = MathUtils::Status::MaxIterations; + double aCheckF[2]; + double aCheckJ[2][2]; + if (!theFunc(aRes.X[0], aRes.X[1], aCheckF, aCheckJ)) + { + aRes.Status = MathUtils::Status::NumericalError; + return aRes; + } + + aRes.ResidualNorm = std::sqrt(aCheckF[0] * aCheckF[0] + aCheckF[1] * aCheckF[1]); + aRes.Status = (aRes.ResidualNorm <= theOptions.FTolerance) ? MathUtils::Status::OK + : MathUtils::Status::MaxIterations; return aRes; } } - double aF1, aF2, aJ11, aJ12, aJ21, aJ22; - if (!theFunc.ValueAndJacobian(aRes.X[0], aRes.X[1], aF1, aF2, aJ11, aJ12, aJ21, aJ22)) + double aF[2]; + double aJ[2][2]; + if (!theFunc(aRes.X[0], aRes.X[1], aF, aJ)) { aRes.Status = MathUtils::Status::NumericalError; return aRes; } - aRes.ResidualNorm = std::sqrt(aF1 * aF1 + aF2 * aF2); + aRes.ResidualNorm = std::sqrt(aF[0] * aF[0] + aF[1] * aF[1]); aRes.Status = (aRes.ResidualNorm <= theOptions.FTolerance) ? MathUtils::Status::OK : MathUtils::Status::MaxIterations; return aRes; diff --git a/src/FoundationClasses/TKMath/MathSys/MathSys_Newton3D.hxx b/src/FoundationClasses/TKMath/MathSys/MathSys_Newton3D.hxx index b538e0d3de..fb09d492c6 100644 --- a/src/FoundationClasses/TKMath/MathSys/MathSys_Newton3D.hxx +++ b/src/FoundationClasses/TKMath/MathSys/MathSys_Newton3D.hxx @@ -247,7 +247,18 @@ NewtonResultN<3> Solve3D(const Function& theFunc, std::max(std::abs(aRes.X[0]), std::max(std::abs(aRes.X[1]), std::abs(aRes.X[2])))); if (aRes.StepNorm <= theOptions.XTolerance * aScaleRef) { - aRes.Status = MathUtils::Status::MaxIterations; + double aCheckF[3]; + double aCheckJ[3][3]; + if (!theFunc(aRes.X[0], aRes.X[1], aRes.X[2], aCheckF, aCheckJ)) + { + aRes.Status = MathUtils::Status::NumericalError; + return aRes; + } + + aRes.ResidualNorm = + std::sqrt(aCheckF[0] * aCheckF[0] + aCheckF[1] * aCheckF[1] + aCheckF[2] * aCheckF[2]); + aRes.Status = (aRes.ResidualNorm <= theOptions.FTolerance) ? MathUtils::Status::OK + : MathUtils::Status::MaxIterations; return aRes; } } diff --git a/src/FoundationClasses/TKMath/MathSys/MathSys_Newton4D.hxx b/src/FoundationClasses/TKMath/MathSys/MathSys_Newton4D.hxx index ba19ae3725..e5e3898a1b 100644 --- a/src/FoundationClasses/TKMath/MathSys/MathSys_Newton4D.hxx +++ b/src/FoundationClasses/TKMath/MathSys/MathSys_Newton4D.hxx @@ -294,7 +294,18 @@ NewtonResultN<4> Solve4D(const Function& theFunc, std::max(std::abs(aRes.X[1]), std::max(std::abs(aRes.X[2]), std::abs(aRes.X[3]))))); if (aRes.StepNorm <= theOptions.XTolerance * aScaleRef) { - aRes.Status = MathUtils::Status::MaxIterations; + double aCheckF[4]; + double aCheckJ[4][4]; + if (!theFunc(aRes.X[0], aRes.X[1], aRes.X[2], aRes.X[3], aCheckF, aCheckJ)) + { + aRes.Status = MathUtils::Status::NumericalError; + return aRes; + } + + aRes.ResidualNorm = std::sqrt(aCheckF[0] * aCheckF[0] + aCheckF[1] * aCheckF[1] + + aCheckF[2] * aCheckF[2] + aCheckF[3] * aCheckF[3]); + aRes.Status = (aRes.ResidualNorm <= theOptions.FTolerance) ? MathUtils::Status::OK + : MathUtils::Status::MaxIterations; return aRes; } } diff --git a/src/FoundationClasses/TKMath/MathUtils/FILES.cmake b/src/FoundationClasses/TKMath/MathUtils/FILES.cmake index 1dc319d352..7a72561ef2 100644 --- a/src/FoundationClasses/TKMath/MathUtils/FILES.cmake +++ b/src/FoundationClasses/TKMath/MathUtils/FILES.cmake @@ -13,6 +13,7 @@ set(OCCT_MathUtils_FILES MathUtils_Random.hxx MathUtils_Bracket.hxx MathUtils_Gauss.hxx + MathUtils_Gauss.cxx MathUtils_Deriv.hxx MathUtils_LineSearch.hxx MathUtils_GaussKronrodWeights.hxx diff --git a/src/FoundationClasses/TKMath/MathUtils/MathUtils_Bracket.hxx b/src/FoundationClasses/TKMath/MathUtils/MathUtils_Bracket.hxx index fe5a46c417..a5e293fb91 100644 --- a/src/FoundationClasses/TKMath/MathUtils/MathUtils_Bracket.hxx +++ b/src/FoundationClasses/TKMath/MathUtils/MathUtils_Bracket.hxx @@ -16,6 +16,7 @@ #include +#include #include #include @@ -105,26 +106,109 @@ struct MinBracketResult double Fc = 0.0; //!< Function value at C }; +//! Options for minimum bracketing. +struct MinBracketOptions +{ + int MaxIterations = 50; //!< Maximum iterations + bool UseLimits = false; //!< Enable hard limits for parameter + double LeftLimit = 0.0; //!< Left hard limit (inclusive) + double RightLimit = 0.0; //!< Right hard limit (inclusive) + bool HasFA = false; //!< True if FA is precomputed + bool HasFB = false; //!< True if FB is precomputed + double FA = 0.0; //!< Precomputed f(A) + double FB = 0.0; //!< Precomputed f(B) +}; + +namespace detail +{ +inline double Limited(double theValue, const MinBracketOptions& theOptions) +{ + if (!theOptions.UseLimits) + { + return theValue; + } + return std::max(theOptions.LeftLimit, std::min(theOptions.RightLimit, theValue)); +} + +template +bool LimitAndMayBeSwap(Function& theFunc, + const MinBracketOptions& theOptions, + const double theA, + double& theB, + double& theFB, + double& theC, + double& theFC) +{ + theC = Limited(theC, theOptions); + if (std::abs(theB - theC) < THE_ZERO_TOL) + { + return false; + } + if (!theFunc.Value(theC, theFC)) + { + return false; + } + + // Keep B between A and C + if ((theA - theB) * (theB - theC) < 0.0) + { + std::swap(theB, theC); + std::swap(theFB, theFC); + } + return true; +} +} // namespace detail + //! Bracket a minimum by finding three points a < b < c with f(b) < f(a) and f(b) < f(c). //! Uses golden section expansion with parabolic interpolation. //! @tparam Function type with Value(double theX, double& theF) method //! @param theFunc function to bracket //! @param theA initial point A //! @param theB initial point B (should be to the right of A in descent direction) -//! @param theMaxIter maximum iterations +//! @param theOptions bracketing options //! @return bracketing result template -MinBracketResult BracketMinimum(Function& theFunc, double theA, double theB, int theMaxIter = 50) +MinBracketResult BracketMinimum(Function& theFunc, + double theA, + double theB, + const MinBracketOptions& theOptions = MinBracketOptions()) { MinBracketResult aResult; - aResult.A = theA; - aResult.B = theB; - - if (!theFunc.Value(aResult.A, aResult.Fa)) + if (theOptions.MaxIterations < 1) { return aResult; } - if (!theFunc.Value(aResult.B, aResult.Fb)) + if (theOptions.UseLimits && theOptions.LeftLimit > theOptions.RightLimit) + { + return aResult; + } + + aResult.A = detail::Limited(theA, theOptions); + aResult.B = detail::Limited(theB, theOptions); + if (std::abs(aResult.A - aResult.B) < THE_ZERO_TOL) + { + return aResult; + } + + const bool isUseFA = + theOptions.HasFA && (!theOptions.UseLimits || std::abs(aResult.A - theA) < THE_ZERO_TOL); + const bool isUseFB = + theOptions.HasFB && (!theOptions.UseLimits || std::abs(aResult.B - theB) < THE_ZERO_TOL); + + if (isUseFA) + { + aResult.Fa = theOptions.FA; + } + else if (!theFunc.Value(aResult.A, aResult.Fa)) + { + return aResult; + } + + if (isUseFB) + { + aResult.Fb = theOptions.FB; + } + else if (!theFunc.Value(aResult.B, aResult.Fb)) { return aResult; } @@ -138,13 +222,26 @@ MinBracketResult BracketMinimum(Function& theFunc, double theA, double theB, int // Initial guess for C using golden ratio aResult.C = aResult.B + THE_GOLDEN_RATIO * (aResult.B - aResult.A); - if (!theFunc.Value(aResult.C, aResult.Fc)) + if (theOptions.UseLimits) + { + if (!detail::LimitAndMayBeSwap(theFunc, + theOptions, + aResult.A, + aResult.B, + aResult.Fb, + aResult.C, + aResult.Fc)) + { + return aResult; + } + } + else if (!theFunc.Value(aResult.C, aResult.Fc)) { return aResult; } // Keep expanding until we bracket a minimum - for (int anIter = 0; anIter < theMaxIter && aResult.Fb >= aResult.Fc; ++anIter) + for (int anIter = 0; anIter < theOptions.MaxIterations && aResult.Fb >= aResult.Fc; ++anIter) { // Parabolic extrapolation const double aR = (aResult.B - aResult.A) * (aResult.Fb - aResult.Fc); @@ -153,8 +250,12 @@ MinBracketResult BracketMinimum(Function& theFunc, double theA, double theB, int double aU = aResult.B - ((aResult.B - aResult.C) * aQ - (aResult.B - aResult.A) * aR) / aDenom; - const double aULim = aResult.B + 100.0 * (aResult.C - aResult.B); - double aFu = 0.0; + double aULim = aResult.B + 100.0 * (aResult.C - aResult.B); + if (theOptions.UseLimits) + { + aULim = detail::Limited(aULim, theOptions); + } + double aFu = 0.0; if ((aResult.B - aU) * (aU - aResult.C) > 0.0) { @@ -183,7 +284,20 @@ MinBracketResult BracketMinimum(Function& theFunc, double theA, double theB, int // Parabolic step didn't help, use golden section aU = aResult.C + THE_GOLDEN_RATIO * (aResult.C - aResult.B); - if (!theFunc.Value(aU, aFu)) + if (theOptions.UseLimits) + { + if (!detail::LimitAndMayBeSwap(theFunc, + theOptions, + aResult.B, + aResult.C, + aResult.Fc, + aU, + aFu)) + { + return aResult; + } + } + else if (!theFunc.Value(aU, aFu)) { return aResult; } @@ -191,7 +305,20 @@ MinBracketResult BracketMinimum(Function& theFunc, double theA, double theB, int else if ((aResult.C - aU) * (aU - aULim) > 0.0) { // U is between C and limit - if (!theFunc.Value(aU, aFu)) + if (theOptions.UseLimits) + { + if (!detail::LimitAndMayBeSwap(theFunc, + theOptions, + aResult.B, + aResult.C, + aResult.Fc, + aU, + aFu)) + { + return aResult; + } + } + else if (!theFunc.Value(aU, aFu)) { return aResult; } @@ -203,7 +330,20 @@ MinBracketResult BracketMinimum(Function& theFunc, double theA, double theB, int aU = aResult.C + THE_GOLDEN_RATIO * (aResult.C - aResult.B); aResult.Fb = aResult.Fc; aResult.Fc = aFu; - if (!theFunc.Value(aU, aFu)) + if (theOptions.UseLimits) + { + if (!detail::LimitAndMayBeSwap(theFunc, + theOptions, + aResult.B, + aResult.C, + aResult.Fc, + aU, + aFu)) + { + return aResult; + } + } + else if (!theFunc.Value(aU, aFu)) { return aResult; } @@ -213,7 +353,20 @@ MinBracketResult BracketMinimum(Function& theFunc, double theA, double theB, int { // U is beyond limit aU = aULim; - if (!theFunc.Value(aU, aFu)) + if (theOptions.UseLimits) + { + if (!detail::LimitAndMayBeSwap(theFunc, + theOptions, + aResult.B, + aResult.C, + aResult.Fc, + aU, + aFu)) + { + return aResult; + } + } + else if (!theFunc.Value(aU, aFu)) { return aResult; } @@ -222,7 +375,20 @@ MinBracketResult BracketMinimum(Function& theFunc, double theA, double theB, int { // Default golden section step aU = aResult.C + THE_GOLDEN_RATIO * (aResult.C - aResult.B); - if (!theFunc.Value(aU, aFu)) + if (theOptions.UseLimits) + { + if (!detail::LimitAndMayBeSwap(theFunc, + theOptions, + aResult.B, + aResult.C, + aResult.Fc, + aU, + aFu)) + { + return aResult; + } + } + else if (!theFunc.Value(aU, aFu)) { return aResult; } @@ -246,9 +412,23 @@ MinBracketResult BracketMinimum(Function& theFunc, double theA, double theB, int std::swap(aResult.Fa, aResult.Fc); } + if (aResult.IsValid && !(aResult.A < aResult.B && aResult.B < aResult.C)) + { + aResult.IsValid = false; + } + return aResult; } +//! Backward-compatible convenience overload with only max-iterations argument. +template +MinBracketResult BracketMinimum(Function& theFunc, double theA, double theB, int theMaxIter) +{ + MinBracketOptions anOptions; + anOptions.MaxIterations = theMaxIter; + return BracketMinimum(theFunc, theA, theB, anOptions); +} + } // namespace MathUtils #endif // _MathUtils_Bracket_HeaderFile diff --git a/src/FoundationClasses/TKMath/MathUtils/MathUtils_Gauss.cxx b/src/FoundationClasses/TKMath/MathUtils/MathUtils_Gauss.cxx new file mode 100644 index 0000000000..6c6e3d1aa6 --- /dev/null +++ b/src/FoundationClasses/TKMath/MathUtils/MathUtils_Gauss.cxx @@ -0,0 +1,115 @@ +// Copyright (c) 2025 OPEN CASCADE SAS +// +// This file is part of Open CASCADE Technology software library. +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License version 2.1 as published +// by the Free Software Foundation, with special exception defined in the file +// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT +// distribution for complete text of the license and disclaimer of any warranty. +// +// Alternatively, this file may be used under the terms of Open CASCADE +// commercial license or contractual agreement. + +#include + +#include + +#include +#include +#include + +#include +#include + +namespace +{ +struct ValueAndWeight +{ + double Value = 0.0; + double Weight = 0.0; + + bool operator<(const ValueAndWeight& theOther) const { return Value < theOther.Value; } +}; + +bool ComputeGaussLegendre(const int theOrder, math_Vector& thePoints, math_Vector& theWeights) +{ + if (theOrder < 1 || thePoints.Length() != theOrder || theWeights.Length() != theOrder) + { + return false; + } + + try + { + math_Vector aDiag(1, theOrder); + math_Vector aSubDiag(1, theOrder); + + for (int i = 1; i <= theOrder; ++i) + { + aDiag(i) = 0.0; + aSubDiag(i) = 0.0; + if (i > 1) + { + const int aSqrIm1 = (i - 1) * (i - 1); + aSubDiag(i) = std::sqrt(static_cast(aSqrIm1) / (4.0 * aSqrIm1 - 1.0)); + } + } + + const MathLin::EigenResult anEigen = MathLin::EigenTridiagonal(aDiag, aSubDiag); + if (!anEigen.IsDone() || !anEigen.EigenValues.has_value() || !anEigen.EigenVectors.has_value()) + { + return false; + } + + const math_Vector& aEigenValues = *anEigen.EigenValues; + const math_Matrix& aEigenVecs = *anEigen.EigenVectors; + + NCollection_Array1 aValuesAndWeights(1, theOrder); + const int aVecLowerRow = aEigenVecs.LowerRow(); + const int aVecLowerCol = aEigenVecs.LowerCol(); + const int aValLower = aEigenValues.Lower(); + for (int i = 1; i <= theOrder; ++i) + { + const double aWeight = 2.0 * aEigenVecs(aVecLowerRow, aVecLowerCol + i - 1) + * aEigenVecs(aVecLowerRow, aVecLowerCol + i - 1); + aValuesAndWeights(i) = {aEigenValues(aValLower + i - 1), aWeight}; + } + + std::sort(aValuesAndWeights.begin(), aValuesAndWeights.end()); + + const int aPointLower = thePoints.Lower(); + const int aWeightLower = theWeights.Lower(); + for (int i = 1; i <= theOrder; ++i) + { + thePoints(aPointLower + i - 1) = aValuesAndWeights(i).Value; + theWeights(aWeightLower + i - 1) = aValuesAndWeights(i).Weight; + } + + return true; + } + catch (Standard_Failure const&) + { + return false; + } +} + +} // namespace + +//================================================================================================== + +bool MathUtils::GetGaussPointsAndWeights(int theOrder, + math_Vector& thePoints, + math_Vector& theWeights) +{ + if (theOrder < 1 || thePoints.Length() != theOrder || theWeights.Length() != theOrder) + { + return false; + } + + if (theOrder <= math::GaussPointsMax()) + { + return math::OrderedGaussPointsAndWeights(theOrder, thePoints, theWeights); + } + + return ComputeGaussLegendre(theOrder, thePoints, theWeights); +} diff --git a/src/FoundationClasses/TKMath/MathUtils/MathUtils_Gauss.hxx b/src/FoundationClasses/TKMath/MathUtils/MathUtils_Gauss.hxx index b8f26c3cc1..49151a988b 100644 --- a/src/FoundationClasses/TKMath/MathUtils/MathUtils_Gauss.hxx +++ b/src/FoundationClasses/TKMath/MathUtils/MathUtils_Gauss.hxx @@ -14,267 +14,22 @@ #ifndef _MathUtils_Gauss_HeaderFile #define _MathUtils_Gauss_HeaderFile +#include +#include + //! Modern math solver utilities. namespace MathUtils { -//! Gauss-Legendre points for n=3. -inline constexpr double THE_GAUSS_POINTS_3[] = {-0.7745966692414834, 0.0, 0.7745966692414834}; - -//! Gauss-Legendre weights for n=3. -inline constexpr double THE_GAUSS_WEIGHTS_3[] = {0.5555555555555556, - 0.8888888888888888, - 0.5555555555555556}; - -//! Gauss-Legendre points for n=4. -inline constexpr double THE_GAUSS_POINTS_4[] = {-0.8611363115940526, - -0.3399810435848563, - 0.3399810435848563, - 0.8611363115940526}; - -//! Gauss-Legendre weights for n=4. -inline constexpr double THE_GAUSS_WEIGHTS_4[] = {0.3478548451374538, - 0.6521451548625461, - 0.6521451548625461, - 0.3478548451374538}; - -//! Gauss-Legendre points for n=5. -inline constexpr double THE_GAUSS_POINTS_5[] = {-0.9061798459386640, - -0.5384693101056831, - 0.0, - 0.5384693101056831, - 0.9061798459386640}; - -//! Gauss-Legendre weights for n=5. -inline constexpr double THE_GAUSS_WEIGHTS_5[] = {0.2369268850561891, - 0.4786286704993665, - 0.5688888888888889, - 0.4786286704993665, - 0.2369268850561891}; - -//! Gauss-Legendre points for n=6. -inline constexpr double THE_GAUSS_POINTS_6[] = {-0.9324695142031521, - -0.6612093864662645, - -0.2386191860831969, - 0.2386191860831969, - 0.6612093864662645, - 0.9324695142031521}; - -//! Gauss-Legendre weights for n=6. -inline constexpr double THE_GAUSS_WEIGHTS_6[] = {0.1713244923791704, - 0.3607615730481386, - 0.4679139345726910, - 0.4679139345726910, - 0.3607615730481386, - 0.1713244923791704}; - -//! Gauss-Legendre points for n=7. -inline constexpr double THE_GAUSS_POINTS_7[] = {-0.9491079123427585, - -0.7415311855993945, - -0.4058451513773972, - 0.0, - 0.4058451513773972, - 0.7415311855993945, - 0.9491079123427585}; - -//! Gauss-Legendre weights for n=7. -inline constexpr double THE_GAUSS_WEIGHTS_7[] = {0.1294849661688697, - 0.2797053914892766, - 0.3818300505051189, - 0.4179591836734694, - 0.3818300505051189, - 0.2797053914892766, - 0.1294849661688697}; - -//! Gauss-Legendre points for n=8. -inline constexpr double THE_GAUSS_POINTS_8[] = {-0.9602898564975363, - -0.7966664774136268, - -0.5255324099163290, - -0.1834346424956498, - 0.1834346424956498, - 0.5255324099163290, - 0.7966664774136268, - 0.9602898564975363}; - -//! Gauss-Legendre weights for n=8. -inline constexpr double THE_GAUSS_WEIGHTS_8[] = {0.1012285362903763, - 0.2223810344533745, - 0.3137066458778873, - 0.3626837833783620, - 0.3626837833783620, - 0.3137066458778873, - 0.2223810344533745, - 0.1012285362903763}; - -//! Gauss-Legendre points for n=10. -inline constexpr double THE_GAUSS_POINTS_10[] = {-0.9739065285171717, - -0.8650633666889845, - -0.6794095682990244, - -0.4333953941292472, - -0.1488743389816312, - 0.1488743389816312, - 0.4333953941292472, - 0.6794095682990244, - 0.8650633666889845, - 0.9739065285171717}; - -//! Gauss-Legendre weights for n=10. -inline constexpr double THE_GAUSS_WEIGHTS_10[] = {0.0666713443086881, - 0.1494513491505806, - 0.2190863625159820, - 0.2692667193099963, - 0.2955242247147529, - 0.2955242247147529, - 0.2692667193099963, - 0.2190863625159820, - 0.1494513491505806, - 0.0666713443086881}; - -//! Gauss-Legendre points for n=15. -inline constexpr double THE_GAUSS_POINTS_15[] = {-0.9879925180204854, - -0.9372733924007060, - -0.8482065834104272, - -0.7244177313601701, - -0.5709721726085388, - -0.3941513470775634, - -0.2011940939974345, - 0.0, - 0.2011940939974345, - 0.3941513470775634, - 0.5709721726085388, - 0.7244177313601701, - 0.8482065834104272, - 0.9372733924007060, - 0.9879925180204854}; - -//! Gauss-Legendre weights for n=15. -inline constexpr double THE_GAUSS_WEIGHTS_15[] = {0.0307532419961173, - 0.0703660474881081, - 0.1071592204671719, - 0.1395706779261543, - 0.1662692058169939, - 0.1861610000155622, - 0.1984314853271116, - 0.2025782419255613, - 0.1984314853271116, - 0.1861610000155622, - 0.1662692058169939, - 0.1395706779261543, - 0.1071592204671719, - 0.0703660474881081, - 0.0307532419961173}; - -//! Gauss-Legendre points for n=21. -inline constexpr double THE_GAUSS_POINTS_21[] = {-0.9937521706203895, - -0.9672268385663063, - -0.9200993341504008, - -0.8533633645833173, - -0.7684399634756779, - -0.6671388041974123, - -0.5516188358872198, - -0.4243421202074388, - -0.2880213168024011, - -0.1455618541608951, - 0.0, - 0.1455618541608951, - 0.2880213168024011, - 0.4243421202074388, - 0.5516188358872198, - 0.6671388041974123, - 0.7684399634756779, - 0.8533633645833173, - 0.9200993341504008, - 0.9672268385663063, - 0.9937521706203895}; - -//! Gauss-Legendre weights for n=21. -inline constexpr double THE_GAUSS_WEIGHTS_21[] = { - 0.0160172282577743, 0.0369537897708525, 0.0571344254268572, 0.0761001136283793, - 0.0934444234560339, 0.1087972991671484, 0.1218314160537285, 0.1322689386333375, - 0.1398873947910731, 0.1445244039899700, 0.1460811336496904, 0.1445244039899700, - 0.1398873947910731, 0.1322689386333375, 0.1218314160537285, 0.1087972991671484, - 0.0934444234560339, 0.0761001136283793, 0.0571344254268572, 0.0369537897708525, - 0.0160172282577743}; - -//! Gauss-Legendre points for n=31 (high precision). -inline constexpr double THE_GAUSS_POINTS_31[] = { - -0.9970874818194770, -0.9846859096651652, -0.9625039250929496, -0.9307569978966481, - -0.8897600299482696, -0.8399203201462673, -0.7817331484166244, -0.7157767845868534, - -0.6427067229242604, -0.5632491614071489, -0.4781937820449025, -0.3883859016082329, - -0.2947180699817016, -0.1981211993355706, -0.0995553121523415, 0.0, - 0.0995553121523415, 0.1981211993355706, 0.2947180699817016, 0.3883859016082329, - 0.4781937820449025, 0.5632491614071489, 0.6427067229242604, 0.7157767845868534, - 0.7817331484166244, 0.8399203201462673, 0.8897600299482696, 0.9307569978966481, - 0.9625039250929496, 0.9846859096651652, 0.9970874818194770}; - -//! Gauss-Legendre weights for n=31. -inline constexpr double THE_GAUSS_WEIGHTS_31[] = { - 0.0074708315792487, 0.0172953547354097, 0.0269785893254440, 0.0364259099519139, - 0.0455433538665749, 0.0542378613250555, 0.0624191330972525, 0.0700003462636801, - 0.0768994045904914, 0.0830398923041908, 0.0883519271671607, 0.0927724753653041, - 0.0962462948268430, 0.0987263019095116, 0.1001737388011984, 0.1005588858060619, - 0.1001737388011984, 0.0987263019095116, 0.0962462948268430, 0.0927724753653041, - 0.0883519271671607, 0.0830398923041908, 0.0768994045904914, 0.0700003462636801, - 0.0624191330972525, 0.0542378613250555, 0.0455433538665749, 0.0364259099519139, - 0.0269785893254440, 0.0172953547354097, 0.0074708315792487}; - -//! Get Gauss-Legendre points and weights for given order. -//! @param theOrder number of quadrature points (3, 4, 5, 6, 7, 8, 10, 15, 21, or 31) -//! @param[out] thePoints pointer to points array -//! @param[out] theWeights pointer to weights array -//! @return true if order is supported -inline bool GetGaussPointsAndWeights(int theOrder, - const double*& thePoints, - const double*& theWeights) -{ - switch (theOrder) - { - case 3: - thePoints = THE_GAUSS_POINTS_3; - theWeights = THE_GAUSS_WEIGHTS_3; - return true; - case 4: - thePoints = THE_GAUSS_POINTS_4; - theWeights = THE_GAUSS_WEIGHTS_4; - return true; - case 5: - thePoints = THE_GAUSS_POINTS_5; - theWeights = THE_GAUSS_WEIGHTS_5; - return true; - case 6: - thePoints = THE_GAUSS_POINTS_6; - theWeights = THE_GAUSS_WEIGHTS_6; - return true; - case 7: - thePoints = THE_GAUSS_POINTS_7; - theWeights = THE_GAUSS_WEIGHTS_7; - return true; - case 8: - thePoints = THE_GAUSS_POINTS_8; - theWeights = THE_GAUSS_WEIGHTS_8; - return true; - case 10: - thePoints = THE_GAUSS_POINTS_10; - theWeights = THE_GAUSS_WEIGHTS_10; - return true; - case 15: - thePoints = THE_GAUSS_POINTS_15; - theWeights = THE_GAUSS_WEIGHTS_15; - return true; - case 21: - thePoints = THE_GAUSS_POINTS_21; - theWeights = THE_GAUSS_WEIGHTS_21; - return true; - case 31: - thePoints = THE_GAUSS_POINTS_31; - theWeights = THE_GAUSS_WEIGHTS_31; - return true; - default: - thePoints = nullptr; - theWeights = nullptr; - return false; - } -} +//! Get ordered Gauss-Legendre points and weights for given order. +//! Points are returned in ascending order on [-1, 1]. +//! @param theOrder number of quadrature points (>= 1) +//! @param[out] thePoints points array +//! @param[out] theWeights weights array +//! @return true if points/weights are available +Standard_EXPORT bool GetGaussPointsAndWeights(int theOrder, + math_Vector& thePoints, + math_Vector& theWeights); } // namespace MathUtils diff --git a/src/FoundationClasses/TKMath/MathUtils/MathUtils_GaussKronrodWeights.cxx b/src/FoundationClasses/TKMath/MathUtils/MathUtils_GaussKronrodWeights.cxx index 6ac1a02394..17edabe9cb 100644 --- a/src/FoundationClasses/TKMath/MathUtils/MathUtils_GaussKronrodWeights.cxx +++ b/src/FoundationClasses/TKMath/MathUtils/MathUtils_GaussKronrodWeights.cxx @@ -12,6 +12,7 @@ // commercial license or contractual agreement. #include "MathUtils_GaussKronrodWeights.hxx" +#include #include //================================================================================================= @@ -29,5 +30,9 @@ bool MathUtils::GetOrderedGaussPointsAndWeights(int theNbGauss, math_Vector& thePoints, math_Vector& theWeights) { - return math::OrderedGaussPointsAndWeights(theNbGauss, thePoints, theWeights); + if (theNbGauss < 1 || thePoints.Length() != theNbGauss || theWeights.Length() != theNbGauss) + { + return false; + } + return MathUtils::GetGaussPointsAndWeights(theNbGauss, thePoints, theWeights); } diff --git a/src/FoundationClasses/TKMath/MathUtils/MathUtils_Types.hxx b/src/FoundationClasses/TKMath/MathUtils/MathUtils_Types.hxx index b95d541f08..323ed7ed94 100644 --- a/src/FoundationClasses/TKMath/MathUtils/MathUtils_Types.hxx +++ b/src/FoundationClasses/TKMath/MathUtils/MathUtils_Types.hxx @@ -111,6 +111,21 @@ struct LinearResult explicit operator bool() const { return IsDone(); } }; +//! Result for multiple linear systems solving (AX = B with matrix RHS). +//! Contains the full solution matrix and determinant if computed. +struct LinearMultipleResult +{ + MathUtils::Status Status = MathUtils::Status::NotConverged; //!< Computation status + std::optional Solutions; //!< Solution matrix X in AX = B (set by solver) + std::optional Determinant; //!< Determinant of matrix (if computed) + + //! Returns true if computation succeeded. + bool IsDone() const { return Status == MathUtils::Status::OK; } + + //! Conversion to bool for convenient checking. + explicit operator bool() const { return IsDone(); } +}; + //! Result for eigenvalue/eigenvector computation. //! Contains eigenvalues and optionally eigenvectors. struct EigenResult diff --git a/src/FoundationClasses/TKMath/MathUtils/README.md b/src/FoundationClasses/TKMath/MathUtils/README.md index 569c440e7c..5135489a9f 100644 --- a/src/FoundationClasses/TKMath/MathUtils/README.md +++ b/src/FoundationClasses/TKMath/MathUtils/README.md @@ -23,8 +23,8 @@ The MathUtils package provides foundational utilities used by all other modern m - `MathUtils_Convergence.hxx` - Convergence testing utilities - `MathUtils_Poly.hxx` - Polynomial evaluation and manipulation - `MathUtils_Domain.hxx` - 1D/2D parameter domain helpers (contains/clamp/normalize/equality checks) -- `MathUtils_Bracket.hxx` - Root and minimum bracketing algorithms -- `MathUtils_Gauss.hxx` - Gauss-Legendre quadrature points and weights +- `MathUtils_Bracket.hxx` - Root and minimum bracketing algorithms (including bounded options for minima) +- `MathUtils_Gauss.hxx` - Gauss-Legendre quadrature points and weights (orders >= 1) - `MathUtils_Deriv.hxx` - Numerical differentiation utilities - `MathUtils_LineSearch.hxx` - Line search algorithms for optimization @@ -93,6 +93,7 @@ via `MathSys_NewtonTypes.hxx`. - `ScalarResult` - For 1D root finding results - `PolyResult` - For polynomial root results (up to 4 roots) - `VectorResult` - For N-D optimization results +- `LinearMultipleResult` - For linear systems with matrix right-hand side (`AX=B`) - `IntegResult` - For integration results with error estimate ### Configuration