lili/OCCT
1
0
Fork 0
mirror of https://github.com/Open-Cascade-SAS/OCCT.git synced 2026-07-06 08:20:46 +08:00
Code Issues Packages Projects Releases Wiki Activity

Foundation Classes - align modern Math* APIs with legacy math_* behavior (#1134)

Browse Source 
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.
This commit is contained in:
Pasukhin Dmitry
2026-03-03 14:18:18 +00:00
committed by GitHub
parent 33688d1049
commit ae7e259e17
26 changed files with 1278 additions and 367 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/FoundationClasses/TKMath/GTests/FILES.cmake
View File

@@ -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

13
src/FoundationClasses/TKMath/GTests/MathInteg_Comparison_Test.cxx
View File

@@ -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
// ============================================================================

52
src/FoundationClasses/TKMath/GTests/MathInteg_Test.cxx
View File

@@ -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<bool>(aResult));
}

284
src/FoundationClasses/TKMath/GTests/MathLin_EigenSearch_Test.cxx Normal file
View File

@@ -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 <gtest/gtest.h>
#include <MathLin_EigenSearch.hxx>
#include <math_EigenValuesSearcher.hxx>
#include <NCollection_Array1.hxx>
#include <math_Matrix.hxx>
#include <math_Vector.hxx>
#include <algorithm>
#include <cmath>
#include <random>
#include <vector>
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<double>& theDiagLegacy,
NCollection_Array1<double>& 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<double> SortedEigenValuesFromLegacy(const math_EigenValuesSearcher& theLegacy)
{
const int aN = theLegacy.Dimension();
std::vector<double> aVals;
aVals.reserve(static_cast<size_t>(aN));
for (int i = 1; i <= aN; ++i)
{
aVals.push_back(theLegacy.EigenValue(i));
}
std::sort(aVals.begin(), aVals.end());
return aVals;
}
std::vector<double> SortedEigenValuesFromModern(const MathLin::EigenResult& theModern)
{
std::vector<double> aVals;
if (!theModern.EigenValues.has_value())
{
return aVals;
}
const math_Vector& aEig = *theModern.EigenValues;
aVals.reserve(static_cast<size_t>(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<double> aDiagLegacy(1, 3);
NCollection_Array1<double> 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<double> aLegacySorted = SortedEigenValuesFromLegacy(aLegacy);
const std::vector<double> 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<double> aDiagLegacy(-2, 2);
NCollection_Array1<double> 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<double> aLegacySorted = SortedEigenValuesFromLegacy(aLegacy);
const std::vector<double> 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<int> aDimDist(2, 32);
std::uniform_real_distribution<double> 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<double> aDiagLegacy(1, aN);
NCollection_Array1<double> 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<double> aLegacySorted = SortedEigenValuesFromLegacy(aLegacy);
const std::vector<double> 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;
}
}
}
}

115
src/FoundationClasses/TKMath/GTests/MathLin_Test.cxx
View File

@@ -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
// ============================================================================

90
src/FoundationClasses/TKMath/GTests/MathSys_Newton2D_Test.cxx
View File

@@ -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);

75
src/FoundationClasses/TKMath/GTests/MathSys_Newton3D_Test.cxx
View File

@@ -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 =

77
src/FoundationClasses/TKMath/GTests/MathSys_Newton4D_Test.cxx
View File

@@ -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*/,

94
src/FoundationClasses/TKMath/GTests/MathUtils_Bracket_Test.cxx Normal file
View File

@@ -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 <gtest/gtest.h>
#include <MathUtils_Bracket.hxx>
#include <cmath>
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);
}

49
src/FoundationClasses/TKMath/MathInteg/MathInteg_Gauss.hxx
View File

@@ -19,6 +19,7 @@
#include <MathUtils_Core.hxx>
#include <MathUtils_Gauss.hxx>
#include <algorithm>
#include <cmath>
//! 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 <typename Function>
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<size_t>(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<size_t>(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 <typename Function>
IntegResult GaussComposite(Function& theFunc,

6
src/FoundationClasses/TKMath/MathInteg/MathInteg_Multiple.hxx
View File

@@ -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)
{

6
src/FoundationClasses/TKMath/MathInteg/MathInteg_Set.hxx
View File

@@ -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);

4
src/FoundationClasses/TKMath/MathInteg/README.md
View File

@@ -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

18
src/FoundationClasses/TKMath/MathLin/MathLin_Gauss.hxx
View File

@@ -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;
}

61
src/FoundationClasses/TKMath/MathLin/MathLin_Householder.hxx
View File

@@ -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;
}

7
src/FoundationClasses/TKMath/MathLin/README.md
View File

@@ -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:

42
src/FoundationClasses/TKMath/MathSys/MathSys_Newton2D.hxx
View File

@@ -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 <typename Function>
NewtonResultN<2> Solve2D(const Function& theFunc,
const std::array<double, 2>& theX0,
@@ -193,14 +191,15 @@ NewtonResultN<2> Solve2D(const Function& theFunc,
{
aRes.NbIterations = static_cast<size_t>(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;

13
src/FoundationClasses/TKMath/MathSys/MathSys_Newton3D.hxx
View File

@@ -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;
}
}

13
src/FoundationClasses/TKMath/MathSys/MathSys_Newton4D.hxx
View File

@@ -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;
}
}

1
src/FoundationClasses/TKMath/MathUtils/FILES.cmake
View File

@@ -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

212
src/FoundationClasses/TKMath/MathUtils/MathUtils_Bracket.hxx
View File

@@ -16,6 +16,7 @@
#include <MathUtils_Core.hxx>
#include <algorithm>
#include <cmath>
#include <utility>
@@ -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 <typename Function>
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 <typename Function>
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 <typename Function>
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

115
src/FoundationClasses/TKMath/MathUtils/MathUtils_Gauss.cxx Normal file
View File

@@ -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 <MathUtils_Gauss.hxx>
#include <MathLin_EigenSearch.hxx>
#include <math.hxx>
#include <NCollection_Array1.hxx>
#include <Standard_Failure.hxx>
#include <algorithm>
#include <cmath>
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<double>(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<ValueAndWeight> 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);
}

269
src/FoundationClasses/TKMath/MathUtils/MathUtils_Gauss.hxx
View File

@@ -14,267 +14,22 @@
#ifndef _MathUtils_Gauss_HeaderFile
#define _MathUtils_Gauss_HeaderFile
#include <Standard_Macro.hxx>
#include <math_Vector.hxx>
//! 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

7
src/FoundationClasses/TKMath/MathUtils/MathUtils_GaussKronrodWeights.cxx
View File

@@ -12,6 +12,7 @@
// commercial license or contractual agreement.
#include "MathUtils_GaussKronrodWeights.hxx"
#include <MathUtils_Gauss.hxx>
#include <math.hxx>
//=================================================================================================
@@ -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);
}

15
src/FoundationClasses/TKMath/MathUtils/MathUtils_Types.hxx
View File

@@ -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<math_Matrix> Solutions; //!< Solution matrix X in AX = B (set by solver)
std::optional<double> 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

5
src/FoundationClasses/TKMath/MathUtils/README.md
View File

@@ -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