diff --git a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib.cxx b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib.cxx index b3de14e41a..e9fa1d8228 100644 --- a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib.cxx +++ b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib.cxx @@ -361,17 +361,16 @@ int BSplCLib::MinKnotMult(const Array1OfInteger& Mults, const int FromK1, const int BSplCLib::NbPoles(const int Degree, const bool Periodic, const NCollection_Array1& Mults) { - int i, sigma = 0; - int f = Mults.Lower(); - int l = Mults.Upper(); + const int f = Mults.Lower(); + const int l = Mults.Upper(); const int* pmu = &Mults(f); - pmu -= f; - int Mf = pmu[f]; - int Ml = pmu[l]; + const int Mf = pmu[0]; + const int Ml = pmu[l - f]; if (Mf <= 0) return 0; if (Ml <= 0) return 0; + int sigma; if (Periodic) { if (Mf > Degree) @@ -384,7 +383,7 @@ int BSplCLib::NbPoles(const int Degree, const bool Periodic, const NCollection_A } else { - int Deg1 = Degree + 1; + const int Deg1 = Degree + 1; if (Mf > Deg1) return 0; if (Ml > Deg1) @@ -392,7 +391,7 @@ int BSplCLib::NbPoles(const int Degree, const bool Periodic, const NCollection_A sigma = Mf + Ml - Deg1; } - for (i = f + 1; i < l; i++) + for (int i = 1; i < l - f; i++) { if (pmu[i] <= 0) return 0; @@ -3008,39 +3007,45 @@ int BSplCLib::SolveBandedSystem(const math_Matrix& Matrix, return 1; } - double* PolesArray = &Array; + double* PolesArray = &Array; + const int aLRow = Matrix.LowerRow(); + const int aNCols = Matrix.ColNumber(); + const double* aMatData = &Matrix(aLRow, 1); - for (int ii = Matrix.LowerRow() + 1; ii <= Matrix.UpperRow(); ii++) + for (int ii = aLRow + 1; ii <= Matrix.UpperRow(); ii++) { - const int aMinIndex = - (ii - LowerBandWidth >= Matrix.LowerRow() ? ii - LowerBandWidth : Matrix.LowerRow()); + const int aMinIndex = (ii - LowerBandWidth >= aLRow ? ii - LowerBandWidth : aLRow); + const double* aRowII = aMatData + (ii - aLRow) * aNCols; for (int jj = aMinIndex; jj < ii; jj++) { + const double aCoeff = aRowII[jj - ii + LowerBandWidth]; for (int kk = 0; kk < ArrayDimension; kk++) { PolesArray[(ii - 1) * ArrayDimension + kk] += - PolesArray[(jj - 1) * ArrayDimension + kk] * Matrix(ii, jj - ii + LowerBandWidth + 1); + PolesArray[(jj - 1) * ArrayDimension + kk] * aCoeff; } } } - for (int ii = Matrix.UpperRow(); ii >= Matrix.LowerRow(); ii--) + for (int ii = Matrix.UpperRow(); ii >= aLRow; ii--) { const int aMaxIndex = (ii + UpperBandWidth <= Matrix.UpperRow() ? ii + UpperBandWidth : Matrix.UpperRow()); + const double* aRowII = aMatData + (ii - aLRow) * aNCols; for (int jj = aMaxIndex; jj > ii; jj--) { + const double aCoeff = aRowII[jj - ii + LowerBandWidth]; for (int kk = 0; kk < ArrayDimension; kk++) { PolesArray[(ii - 1) * ArrayDimension + kk] -= - PolesArray[(jj - 1) * ArrayDimension + kk] * Matrix(ii, jj - ii + LowerBandWidth + 1); + PolesArray[(jj - 1) * ArrayDimension + kk] * aCoeff; } } // Fixing a bug PRO18577 to avoid division by zero - const double aDivisor = Matrix(ii, LowerBandWidth + 1); + const double aDivisor = aRowII[LowerBandWidth]; constexpr double THE_TOLERANCE = 1.0e-16; if (std::abs(aDivisor) <= THE_TOLERANCE) { @@ -3147,9 +3152,14 @@ void BSplCLib::Interpolate(const int Degree, double& Poles, int& InversionProblem) { - int ErrorCode, UpperBandWidth, LowerBandWidth; - // double *PolesArray = &Poles ; - math_Matrix InterpolationMatrix(1, Parameters.Length(), 1, 2 * Degree + 1); + int ErrorCode, UpperBandWidth, LowerBandWidth; + constexpr int THE_STACK_LIMIT = 2048; + const int aMatSize = Parameters.Length() * (2 * Degree + 1); + double aStackBuf[THE_STACK_LIMIT]; + math_Matrix InterpolationMatrix = + (aMatSize <= THE_STACK_LIMIT) + ? math_Matrix(aStackBuf, 1, Parameters.Length(), 1, 2 * Degree + 1) + : math_Matrix(1, Parameters.Length(), 1, 2 * Degree + 1); ErrorCode = BSplCLib::BuildBSpMatrix(Parameters, ContactOrderArray, FlatKnots, @@ -3187,9 +3197,14 @@ void BSplCLib::Interpolate(const int Degree, double& Weights, int& InversionProblem) { - int ErrorCode, UpperBandWidth, LowerBandWidth; - - math_Matrix InterpolationMatrix(1, Parameters.Length(), 1, 2 * Degree + 1); + int ErrorCode, UpperBandWidth, LowerBandWidth; + constexpr int THE_STACK_LIMIT = 2048; + const int aMatSize = Parameters.Length() * (2 * Degree + 1); + double aStackBuf[THE_STACK_LIMIT]; + math_Matrix InterpolationMatrix = + (aMatSize <= THE_STACK_LIMIT) + ? math_Matrix(aStackBuf, 1, Parameters.Length(), 1, 2 * Degree + 1) + : math_Matrix(1, Parameters.Length(), 1, 2 * Degree + 1); ErrorCode = BSplCLib::BuildBSpMatrix(Parameters, ContactOrderArray, FlatKnots, diff --git a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_2.cxx b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_2.cxx index 7827f07c52..16ba056a1d 100644 --- a/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_2.cxx +++ b/src/FoundationClasses/TKMath/BSplCLib/BSplCLib_2.cxx @@ -316,12 +316,20 @@ int BSplCLib::BuildBSpMatrix(const NCollection_Array1& Parameters, return 1; } - math_Matrix aBSplineBasis(1, aMaxOrder, 1, aMaxOrder); + double aBasisBuf[aMaxOrder * aMaxOrder]; + math_Matrix aBSplineBasis(aBasisBuf, 1, aMaxOrder, 1, aMaxOrder); + + // Zero the entire matrix once instead of per-row zero-fill loops. + Matrix.Init(0.0); + const int aMatLRow = Matrix.LowerRow(); + const int aMatNCols = Matrix.ColNumber(); + double* aMatData = &Matrix(aMatLRow, 1); for (int i = Parameters.Lower(); i <= Parameters.Upper(); i++) { int aFirstNonZeroIndex = 0; - const int anErrorCode = BSplCLib::EvalBsplineBasis(ContactOrderArray(i), + const int aContactOrder = ContactOrderArray(i); + const int anErrorCode = BSplCLib::EvalBsplineBasis(aContactOrder, anOrder, FlatKnots, Parameters(i), @@ -332,19 +340,12 @@ int BSplCLib::BuildBSpMatrix(const NCollection_Array1& Parameters, return 2; } - int anIndex = LowerBandWidth + 1 + aFirstNonZeroIndex - i; - for (int j = 1; j < anIndex; j++) + int anIndex = LowerBandWidth + 1 + aFirstNonZeroIndex - i; + double* aRowData = aMatData + (i - aMatLRow) * aMatNCols; + const double* aBasisSrc = aBasisBuf + aContactOrder * aMaxOrder; + for (int j = 0; j < anOrder; j++) { - Matrix.Value(i, j) = 0.0; - } - for (int j = 1; j <= anOrder; j++) - { - Matrix.Value(i, anIndex) = aBSplineBasis(ContactOrderArray(i) + 1, j); - anIndex += 1; - } - for (int j = anIndex; j <= aBandWidth; j++) - { - Matrix.Value(i, j) = 0.0; + aRowData[anIndex + j - 1] = aBasisSrc[j]; } } @@ -361,14 +362,20 @@ int BSplCLib::FactorBandedMatrix(math_Matrix& Matrix, const int aBandWidth = UpperBandWidth + LowerBandWidth + 1; PivotIndexProblem = 0; - for (int i = Matrix.LowerRow() + 1; i <= Matrix.UpperRow(); i++) + const int aLRow = Matrix.LowerRow(); + const int aNCols = Matrix.ColNumber(); + double* aData = &Matrix(aLRow, 1); + + for (int i = aLRow + 1; i <= Matrix.UpperRow(); i++) { const int aMinIndex = (LowerBandWidth - i + 2 >= 1 ? LowerBandWidth - i + 2 : 1); + double* aRowI = aData + (i - aLRow) * aNCols; for (int j = aMinIndex; j <= LowerBandWidth; j++) { - const int anIndex = i - LowerBandWidth + j - 1; - const double aPivot = Matrix(anIndex, LowerBandWidth + 1); + const int anIndex = i - LowerBandWidth + j - 1; + const double* aRowIdx = aData + (anIndex - aLRow) * aNCols; + const double aPivot = aRowIdx[LowerBandWidth]; if (std::abs(aPivot) <= RealSmall()) { PivotIndexProblem = anIndex; @@ -376,12 +383,12 @@ int BSplCLib::FactorBandedMatrix(math_Matrix& Matrix, } const double anInverse = -1.0 / aPivot; - Matrix(i, j) = Matrix(i, j) * anInverse; + aRowI[j - 1] = aRowI[j - 1] * anInverse; const int aMaxIndex = aBandWidth + anIndex - i; for (int k = j + 1; k <= aMaxIndex; k++) { - Matrix(i, k) += Matrix(i, j) * Matrix(anIndex, k + i - anIndex); + aRowI[k - 1] += aRowI[j - 1] * aRowIdx[k + i - anIndex - 1]; } } } @@ -448,8 +455,17 @@ int BSplCLib::EvalBsplineBasis(const int DerivativeReque FirstNonZeroBsplineIndex = ii - Order + 1; - BsplineBasis(1, 1) = 1.0; - aLocalRequest = DerivativeRequest; + // Use raw pointers for BsplineBasis and FlatKnots to bypass accessor overhead + // (math_Matrix::Value -> math_DoubleTab::Value -> NCollection_Array2::Value -> + // NCollection_Array1::at with bounds checks and DYLD stubs in tight loops). + const int aBasisNCols = BsplineBasis.ColNumber(); + double* aBasisData = &BsplineBasis(1, 1); + const double* aKnotsData = &FlatKnots(FlatKnots.Lower()); + constexpr double aResolution = gp::Resolution(); + ii -= FlatKnots.Lower(); // rebase to zero-based indexing into aKnotsData + + aBasisData[0] = 1.0; + aLocalRequest = DerivativeRequest; if (DerivativeRequest >= Order) { aLocalRequest = Order - 1; @@ -457,21 +473,21 @@ int BSplCLib::EvalBsplineBasis(const int DerivativeReque for (int qq = 2; qq <= Order - aLocalRequest; qq++) { - BsplineBasis(1, qq) = 0.0; + aBasisData[qq - 1] = 0.0; for (int pp = 1; pp <= qq - 1; pp++) { - const double aScale = FlatKnots(ii + pp) - FlatKnots(ii - qq + pp + 1); - if (std::abs(aScale) < gp::Resolution()) + const double aScale = aKnotsData[ii + pp] - aKnotsData[ii - qq + pp + 1]; + if (std::abs(aScale) < aResolution) { return 2; } - const double aFactor = (Parameter - FlatKnots(ii - qq + pp + 1)) / aScale; - const double aSaved = aFactor * BsplineBasis(1, pp); - BsplineBasis(1, pp) *= (1.0 - aFactor); - BsplineBasis(1, pp) += BsplineBasis(1, qq); - BsplineBasis(1, qq) = aSaved; + const double aFactor = (Parameter - aKnotsData[ii - qq + pp + 1]) / aScale; + const double aSaved = aFactor * aBasisData[pp - 1]; + aBasisData[pp - 1] *= (1.0 - aFactor); + aBasisData[pp - 1] += aBasisData[qq - 1]; + aBasisData[qq - 1] = aSaved; } } @@ -479,37 +495,38 @@ int BSplCLib::EvalBsplineBasis(const int DerivativeReque { for (int pp = 1; pp <= qq - 1; pp++) { - BsplineBasis(Order - qq + 2, pp) = BsplineBasis(1, pp); + aBasisData[(Order - qq + 1) * aBasisNCols + (pp - 1)] = aBasisData[pp - 1]; } - BsplineBasis(1, qq) = 0.0; + aBasisData[qq - 1] = 0.0; for (int ss = Order - aLocalRequest + 1; ss <= qq; ss++) { - BsplineBasis(Order - ss + 2, qq) = 0.0; + aBasisData[(Order - ss + 1) * aBasisNCols + (qq - 1)] = 0.0; } for (int pp = 1; pp <= qq - 1; pp++) { - const double aScale = FlatKnots(ii + pp) - FlatKnots(ii - qq + pp + 1); - if (std::abs(aScale) < gp::Resolution()) + const double aScale = aKnotsData[ii + pp] - aKnotsData[ii - qq + pp + 1]; + if (std::abs(aScale) < aResolution) { return 2; } const double anInverse = 1.0 / aScale; - const double aFactor = (Parameter - FlatKnots(ii - qq + pp + 1)) * anInverse; - double aSaved = aFactor * BsplineBasis(1, pp); - BsplineBasis(1, pp) *= (1.0 - aFactor); - BsplineBasis(1, pp) += BsplineBasis(1, qq); - BsplineBasis(1, qq) = aSaved; + const double aFactor = (Parameter - aKnotsData[ii - qq + pp + 1]) * anInverse; + double aSaved = aFactor * aBasisData[pp - 1]; + aBasisData[pp - 1] *= (1.0 - aFactor); + aBasisData[pp - 1] += aBasisData[qq - 1]; + aBasisData[qq - 1] = aSaved; const double aLocalInverse = static_cast(qq - 1) * anInverse; for (int ss = Order - aLocalRequest + 1; ss <= qq; ss++) { - aSaved = aLocalInverse * BsplineBasis(Order - ss + 2, pp); - BsplineBasis(Order - ss + 2, pp) *= -aLocalInverse; - BsplineBasis(Order - ss + 2, pp) += BsplineBasis(Order - ss + 2, qq); - BsplineBasis(Order - ss + 2, qq) = aSaved; + double* aRowS = aBasisData + (Order - ss + 1) * aBasisNCols; + aSaved = aLocalInverse * aRowS[pp - 1]; + aRowS[pp - 1] *= -aLocalInverse; + aRowS[pp - 1] += aRowS[qq - 1]; + aRowS[qq - 1] = aSaved; } } } diff --git a/src/FoundationClasses/TKMath/Convert/Convert_CompPolynomialToPoles.cxx b/src/FoundationClasses/TKMath/Convert/Convert_CompPolynomialToPoles.cxx index 57842c0885..b8b9c9f0b3 100644 --- a/src/FoundationClasses/TKMath/Convert/Convert_CompPolynomialToPoles.cxx +++ b/src/FoundationClasses/TKMath/Convert/Convert_CompPolynomialToPoles.cxx @@ -197,7 +197,10 @@ void Convert_CompPolynomialToPoles::Perform(const int Nu myFlatKnots = NCollection_Array1(1, num_flat_knots); BSplCLib::KnotSequence(myKnots, myMults, myDegree, false, myFlatKnots); - NCollection_Array1 parameters(1, num_poles); + constexpr int THE_MAX_POLES = 128; + double aParamBuf[THE_MAX_POLES]; + int aContactBuf[THE_MAX_POLES]; + NCollection_Array1 parameters(aParamBuf[0], 1, num_poles, num_poles <= THE_MAX_POLES); BSplCLib::BuildSchoenbergPoints(myDegree, myFlatKnots, parameters); myPoles = NCollection_Array2(1, num_poles, 1, Dimension); index = 2; @@ -205,7 +208,7 @@ void Convert_CompPolynomialToPoles::Perform(const int Nu Pindex = PolynomialIntervals.LowerRow(); poles_array = (double*)&myPoles.ChangeValue(1, 1); - NCollection_Array1 contact_array(1, num_poles); + NCollection_Array1 contact_array(aContactBuf[0], 1, num_poles, num_poles <= THE_MAX_POLES); poles_index = 0; for (ii = 1; ii <= num_poles; ii++, poles_index += Dimension) diff --git a/src/FoundationClasses/TKMath/math/math_VectorBase.lxx b/src/FoundationClasses/TKMath/math/math_VectorBase.lxx index be2d52556a..c73afc4dd6 100644 --- a/src/FoundationClasses/TKMath/math/math_VectorBase.lxx +++ b/src/FoundationClasses/TKMath/math/math_VectorBase.lxx @@ -101,66 +101,60 @@ void math_VectorBase::SetLower(const int theLower) template double math_VectorBase::Norm() const { - // 4-way unrolled accumulation for better vectorization - double aSum1 = 0.0, aSum2 = 0.0, aSum3 = 0.0, aSum4 = 0.0; - int anIndex = Lower(); - int anUpper = Upper(); - int anUpper4 = anUpper - 3; + const int aLen = Length(); + const TheItemType* aPtr = &Array(Lower()); + double aSum1 = 0.0, aSum2 = 0.0, aSum3 = 0.0, aSum4 = 0.0; + int i = 0; + int aLen4 = aLen - 3; - // Process 4 elements at a time - for (; anIndex <= anUpper4; anIndex += 4) + for (; i < aLen4; i += 4) { - const double aVal0 = static_cast(Array(anIndex)); - const double aVal1 = static_cast(Array(anIndex + 1)); - const double aVal2 = static_cast(Array(anIndex + 2)); - const double aVal3 = static_cast(Array(anIndex + 3)); + const double aVal0 = static_cast(aPtr[i]); + const double aVal1 = static_cast(aPtr[i + 1]); + const double aVal2 = static_cast(aPtr[i + 2]); + const double aVal3 = static_cast(aPtr[i + 3]); aSum1 += aVal0 * aVal0; aSum2 += aVal1 * aVal1; aSum3 += aVal2 * aVal2; aSum4 += aVal3 * aVal3; } - // Process remaining elements - for (; anIndex <= anUpper; ++anIndex) + for (; i < aLen; ++i) { - const double aVal = static_cast(Array(anIndex)); + const double aVal = static_cast(aPtr[i]); aSum1 += aVal * aVal; } - // Combine partial sums (pairwise for better numerical stability) return std::sqrt((aSum1 + aSum2) + (aSum3 + aSum4)); } template double math_VectorBase::Norm2() const { - // 4-way unrolled accumulation for better vectorization - double aSum1 = 0.0, aSum2 = 0.0, aSum3 = 0.0, aSum4 = 0.0; - int anIndex = Lower(); - int anUpper = Upper(); - int anUpper4 = anUpper - 3; + const int aLen = Length(); + const TheItemType* aPtr = &Array(Lower()); + double aSum1 = 0.0, aSum2 = 0.0, aSum3 = 0.0, aSum4 = 0.0; + int i = 0; + int aLen4 = aLen - 3; - // Process 4 elements at a time - for (; anIndex <= anUpper4; anIndex += 4) + for (; i < aLen4; i += 4) { - const double aVal0 = static_cast(Array(anIndex)); - const double aVal1 = static_cast(Array(anIndex + 1)); - const double aVal2 = static_cast(Array(anIndex + 2)); - const double aVal3 = static_cast(Array(anIndex + 3)); + const double aVal0 = static_cast(aPtr[i]); + const double aVal1 = static_cast(aPtr[i + 1]); + const double aVal2 = static_cast(aPtr[i + 2]); + const double aVal3 = static_cast(aPtr[i + 3]); aSum1 += aVal0 * aVal0; aSum2 += aVal1 * aVal1; aSum3 += aVal2 * aVal2; aSum4 += aVal3 * aVal3; } - // Process remaining elements - for (; anIndex <= anUpper; ++anIndex) + for (; i < aLen; ++i) { - const double aVal = static_cast(Array(anIndex)); + const double aVal = static_cast(aPtr[i]); aSum1 += aVal * aVal; } - // Combine partial sums (pairwise for better numerical stability) return (aSum1 + aSum2) + (aSum3 + aSum4); } @@ -446,17 +440,22 @@ void math_VectorBase::Multiply(const math_Matrix& (Length() != theLeft.RowNumber()) || (theLeft.ColNumber() != theRight.Length()), "math_VectorBase::Multiply() - input matrix and /or vector have wrong dimensions"); - int Index = Lower(); - for (int I = theLeft.LowerRow(); I <= theLeft.UpperRow(); I++) + // result[r] = sum_c theLeft(r, c) * theRight[c] + const int aNRows = theLeft.RowNumber(); + const int aNCols = theLeft.ColNumber(); + const double* aMatData = &theLeft(theLeft.LowerRow(), theLeft.LowerCol()); + const TheItemType* aRightPtr = &theRight.Array(theRight.Lower()); + TheItemType* aResPtr = &Array(Lower()); + + for (int r = 0; r < aNRows; r++) { - Array(Index) = 0.0; - int K = theRight.Lower(); - for (int J = theLeft.LowerCol(); J <= theLeft.UpperCol(); J++) + TheItemType aSum = 0.0; + const double* aRowData = aMatData + r * aNCols; + for (int c = 0; c < aNCols; c++) { - Array(Index) = Array(Index) + theLeft(I, J) * theRight.Array(K); - K++; + aSum += aRowData[c] * aRightPtr[c]; } - Index++; + aResPtr[r] = aSum; } } @@ -468,17 +467,21 @@ void math_VectorBase::Multiply(const math_VectorBase& (Length() != theRight.ColNumber()) || (theLeft.Length() != theRight.RowNumber()), "math_VectorBase::Multiply() - input matrix and /or vector have wrong dimensions"); - int Index = Lower(); - for (int J = theRight.LowerCol(); J <= theRight.UpperCol(); J++) + // result[c] = sum_r theLeft[r] * theRight(r, c) + const int aNCols = theRight.ColNumber(); + const int aNRows = theRight.RowNumber(); + const double* aMatData = &theRight(theRight.LowerRow(), theRight.LowerCol()); + const TheItemType* aLeftPtr = &theLeft.Array(theLeft.Lower()); + TheItemType* aResPtr = &Array(Lower()); + + for (int c = 0; c < aNCols; c++) { - Array(Index) = 0.0; - int K = theLeft.Lower(); - for (int I = theRight.LowerRow(); I <= theRight.UpperRow(); I++) + TheItemType aSum = 0.0; + for (int r = 0; r < aNRows; r++) { - Array(Index) = Array(Index) + theLeft.Array(K) * theRight(I, J); - K++; + aSum += aLeftPtr[r] * aMatData[r * aNCols + c]; } - Index++; + aResPtr[c] = aSum; } } @@ -490,17 +493,21 @@ void math_VectorBase::TMultiply(const math_Matrix& (Length() != theTLeft.ColNumber()) || (theTLeft.RowNumber() != theRight.Length()), "math_VectorBase::TMultiply() - input matrix and /or vector have wrong dimensions"); - int Index = Lower(); - for (int I = theTLeft.LowerCol(); I <= theTLeft.UpperCol(); I++) + // result[c] = sum_r theTLeft(r, c) * theRight[r] (transpose-multiply) + const int aNCols = theTLeft.ColNumber(); + const int aNRows = theTLeft.RowNumber(); + const double* aMatData = &theTLeft(theTLeft.LowerRow(), theTLeft.LowerCol()); + const TheItemType* aRightPtr = &theRight.Array(theRight.Lower()); + TheItemType* aResPtr = &Array(Lower()); + + for (int c = 0; c < aNCols; c++) { - Array(Index) = 0.0; - int K = theRight.Lower(); - for (int J = theTLeft.LowerRow(); J <= theTLeft.UpperRow(); J++) + TheItemType aSum = 0.0; + for (int r = 0; r < aNRows; r++) { - Array(Index) = Array(Index) + theTLeft(J, I) * theRight.Array(K); - K++; + aSum += aMatData[r * aNCols + c] * aRightPtr[r]; } - Index++; + aResPtr[c] = aSum; } } @@ -512,17 +519,22 @@ void math_VectorBase::TMultiply(const math_VectorBase& (Length() != theTRight.RowNumber()) || (theLeft.Length() != theTRight.ColNumber()), "math_VectorBase::TMultiply() - input matrix and /or vector have wrong dimensions"); - int Index = Lower(); - for (int J = theTRight.LowerRow(); J <= theTRight.UpperRow(); J++) + // result[r] = sum_c theLeft[c] * theTRight(r, c) + const int aNCols = theTRight.ColNumber(); + const int aNRows = theTRight.RowNumber(); + const double* aMatData = &theTRight(theTRight.LowerRow(), theTRight.LowerCol()); + const TheItemType* aLeftPtr = &theLeft.Array(theLeft.Lower()); + TheItemType* aResPtr = &Array(Lower()); + + for (int r = 0; r < aNRows; r++) { - Array(Index) = 0.0; - int K = theLeft.Lower(); - for (int I = theTRight.LowerCol(); I <= theTRight.UpperCol(); I++) + TheItemType aSum = 0.0; + const double* aRowData = aMatData + r * aNCols; + for (int c = 0; c < aNCols; c++) { - Array(Index) = Array(Index) + theLeft.Array(K) * theTRight(J, I); - K++; + aSum += aLeftPtr[c] * aRowData[c]; } - Index++; + aResPtr[r] = aSum; } } @@ -530,19 +542,19 @@ template TheItemType math_VectorBase::Multiplied( const math_VectorBase& theRight) const { - double Result = 0; - Standard_DimensionError_Raise_if( Length() != theRight.Length(), "math_VectorBase::Multiplied() - input vector has wrong dimensions"); - int I = theRight.Lower(); - for (int Index = Lower(); Index <= Upper(); Index++) + const int aLen = Length(); + const TheItemType* aLeftPtr = &Array(Lower()); + const TheItemType* aRightPtr = &theRight.Array(theRight.Lower()); + TheItemType aResult = 0; + for (int i = 0; i < aLen; i++) { - Result = Result + Array(Index) * theRight.Array(I); - I++; + aResult += aLeftPtr[i] * aRightPtr[i]; } - return Result; + return aResult; } template diff --git a/src/ModelingAlgorithms/TKFillet/BRepBlend/BRepBlend_AppFuncRoot.cxx b/src/ModelingAlgorithms/TKFillet/BRepBlend/BRepBlend_AppFuncRoot.cxx index 67bb9bb169..abae5f11ce 100644 --- a/src/ModelingAlgorithms/TKFillet/BRepBlend/BRepBlend_AppFuncRoot.cxx +++ b/src/ModelingAlgorithms/TKFillet/BRepBlend/BRepBlend_AppFuncRoot.cxx @@ -75,8 +75,14 @@ BRepBlend_AppFuncRoot::BRepBlend_AppFuncRoot(occ::handle& Line, { myBary.SetCoord(0, 0, 0); } + + mySolver = std::make_unique(Func, myTolerance, 30); } +//================================================================================================== + +BRepBlend_AppFuncRoot::~BRepBlend_AppFuncRoot() = default; + //================================================================================ // Function: D0 // Purpose : Calculation of section for v = Param, if calculation fails @@ -337,28 +343,26 @@ bool BRepBlend_AppFuncRoot::SearchPoint(Blend_AppFunction& Func, // (2) Calculation of the solution ------------------------ Func.Set(Param); Func.GetBounds(X1, X2); - math_FunctionSetRoot rsnld(Func, myTolerance, 30); + mySolver->Perform(Func, XInit, X1, X2); - rsnld.Perform(Func, XInit, X1, X2); - - if (!rsnld.IsDone()) + if (!mySolver->IsDone()) { #ifdef BREPBLEND_DEB std::cout << "AppFunc : RNLD Not done en t = " << Param << std::endl; #endif return false; } - rsnld.Root(Sol); + mySolver->Root(Sol); // (3) Storage of the point Point(Func, Param, Sol, Pnt); // (4) Insertion of the point if the calculation seems long. - if ((!Trouve) && (rsnld.NbIterations() > 3)) + if ((!Trouve) && (mySolver->NbIterations() > 3)) { #ifdef OCCT_DEBUG std::cout << "Evaluation in t = " << Param << "given" << std::endl; - rsnld.Dump(std::cout); + mySolver->Dump(std::cout); #endif myLine->InsertBefore(Index + 1, Pnt); } diff --git a/src/ModelingAlgorithms/TKFillet/BRepBlend/BRepBlend_AppFuncRoot.hxx b/src/ModelingAlgorithms/TKFillet/BRepBlend/BRepBlend_AppFuncRoot.hxx index 07e8904e61..1897f9f96f 100644 --- a/src/ModelingAlgorithms/TKFillet/BRepBlend/BRepBlend_AppFuncRoot.hxx +++ b/src/ModelingAlgorithms/TKFillet/BRepBlend/BRepBlend_AppFuncRoot.hxx @@ -30,8 +30,12 @@ #include #include #include + +#include + class BRepBlend_Line; class Blend_AppFunction; +class math_FunctionSetRoot; //! Function to approximate by AppSurface class BRepBlend_AppFuncRoot : public Approx_SweepFunction @@ -152,6 +156,8 @@ public: DEFINE_STANDARD_RTTIEXT(BRepBlend_AppFuncRoot, Approx_SweepFunction) + Standard_EXPORT ~BRepBlend_AppFuncRoot(); + protected: Standard_EXPORT BRepBlend_AppFuncRoot(occ::handle& Line, Blend_AppFunction& Func, @@ -166,15 +172,16 @@ private: const int LastIndex, int& ParamIndex) const; - occ::handle myLine; - void* myFunc; - math_Vector myTolerance; - Blend_Point myPnt; - gp_Pnt myBary; - math_Vector X1; - math_Vector X2; - math_Vector XInit; - math_Vector Sol; + occ::handle myLine; + void* myFunc; + math_Vector myTolerance; + Blend_Point myPnt; + gp_Pnt myBary; + math_Vector X1; + math_Vector X2; + math_Vector XInit; + math_Vector Sol; + std::unique_ptr mySolver; }; #endif // _BRepBlend_AppFuncRoot_HeaderFile diff --git a/src/ModelingAlgorithms/TKGeomAlgo/GeomFill/GeomFill_PolynomialConvertor.cxx b/src/ModelingAlgorithms/TKGeomAlgo/GeomFill/GeomFill_PolynomialConvertor.cxx index 10dfc1e13d..b460155a80 100644 --- a/src/ModelingAlgorithms/TKGeomAlgo/GeomFill/GeomFill_PolynomialConvertor.cxx +++ b/src/ModelingAlgorithms/TKGeomAlgo/GeomFill/GeomFill_PolynomialConvertor.cxx @@ -23,9 +23,7 @@ #include #include #include -#include #include -#include #include GeomFill_PolynomialConvertor::GeomFill_PolynomialConvertor() @@ -43,63 +41,49 @@ bool GeomFill_PolynomialConvertor::Initialized() const void GeomFill_PolynomialConvertor::Init() { if (myinit) - return; // On n'initialise qu'une fois - int ii, jj; - double terme; - math_Matrix H(1, Ordre, 1, Ordre), B(1, Ordre, 1, Ordre); - occ::handle> Coeffs = - new (NCollection_HArray1)(1, Ordre * Ordre), - TrueInter = new (NCollection_HArray1)(1, 2); - - occ::handle> Poles1d = - new (NCollection_HArray2)(1, Ordre, 1, Ordre), - Inter = new (NCollection_HArray2)(1, 1, 1, 2); - - // Calcul de B - Inter->SetValue(1, 1, -1); - Inter->SetValue(1, 2, 1); - TrueInter->SetValue(1, -1); - TrueInter->SetValue(2, 1); - - Coeffs->Init(0); - for (ii = 1; ii <= Ordre; ii++) - { - Coeffs->SetValue(ii + (ii - 1) * Ordre, 1); - } - - // Convertion ancienne formules - occ::handle> Ncf = new (NCollection_HArray1)(1, 1); - Ncf->Init(Ordre); - - Convert_CompPolynomialToPoles AConverter(1, 1, 8, 8, Ncf, Coeffs, Inter, TrueInter); - /* Convert_CompPolynomialToPoles - AConverter(8, Ordre-1, Ordre-1, - Coeffs, - Inter, - TrueInter); En attente du bon Geomlite*/ - Poles1d = new NCollection_HArray2(AConverter.Poles()); - - for (jj = 1; jj <= Ordre; jj++) - { - for (ii = 1; ii <= Ordre; ii++) - { - terme = Poles1d->Value(ii, jj); - if (std::abs(terme - 1) < 1.e-9) - terme = 1; // petite retouche - if (std::abs(terme + 1) < 1.e-9) - terme = -1; - B(ii, jj) = terme; - } - } - // Calcul de H - myinit = PLib::HermiteCoefficients(-1, 1, Ordre / 2 - 1, Ordre / 2 - 1, H); - H.Transpose(); - - if (!myinit) return; - // reste l'essentiel - BH = B * H; + // BH = B * H where B is the monomial-to-BSpline conversion matrix on [-1,1], degree 7, + // and H is the Hermite coefficients matrix. Both are mathematical constants computed once. + static const math_Matrix THE_BH_MATRIX = []() { + constexpr int anOrdre = 8; + NCollection_Array1 aCoeffs(1, anOrdre * anOrdre); + NCollection_Array1 anInter(1, 2), aTrueInter(1, 2); + anInter.SetValue(1, -1); + anInter.SetValue(2, 1); + aTrueInter.SetValue(1, -1); + aTrueInter.SetValue(2, 1); + aCoeffs.Init(0); + for (int ii = 1; ii <= anOrdre; ii++) + aCoeffs.SetValue(ii + (ii - 1) * anOrdre, 1); + + Convert_CompPolynomialToPoles aConverter(anOrdre, + anOrdre - 1, + anOrdre - 1, + aCoeffs, + anInter, + aTrueInter); + const NCollection_Array2& aPoles = aConverter.Poles(); + math_Matrix aB(1, anOrdre, 1, anOrdre); + for (int jj = 1; jj <= anOrdre; jj++) + for (int ii = 1; ii <= anOrdre; ii++) + { + double aTerm = aPoles.Value(ii, jj); + if (std::abs(aTerm - 1) < 1.e-9) + aTerm = 1; + if (std::abs(aTerm + 1) < 1.e-9) + aTerm = -1; + aB(ii, jj) = aTerm; + } + + math_Matrix aH(1, anOrdre, 1, anOrdre); + PLib::HermiteCoefficients(-1, 1, anOrdre / 2 - 1, anOrdre / 2 - 1, aH); + aH.Transpose(); + return math_Matrix(aB * aH); + }(); + + BH = THE_BH_MATRIX; + myinit = true; } void GeomFill_PolynomialConvertor::Section(const gp_Pnt& FirstPnt, diff --git a/src/ModelingAlgorithms/TKGeomAlgo/GeomFill/GeomFill_QuasiAngularConvertor.cxx b/src/ModelingAlgorithms/TKGeomAlgo/GeomFill/GeomFill_QuasiAngularConvertor.cxx index be6466d898..538980fcaa 100644 --- a/src/ModelingAlgorithms/TKGeomAlgo/GeomFill/GeomFill_QuasiAngularConvertor.cxx +++ b/src/ModelingAlgorithms/TKGeomAlgo/GeomFill/GeomFill_QuasiAngularConvertor.cxx @@ -23,7 +23,6 @@ #include #include #include -#include #define NullAngle 1.e-6 @@ -72,43 +71,45 @@ bool GeomFill_QuasiAngularConvertor::Initialized() const void GeomFill_QuasiAngularConvertor::Init() { if (myinit) - return; // On n'initialise qu'une fois - int ii, jj, Ordre = 7; - double terme; - NCollection_Array1 Coeffs(1, Ordre * Ordre), TrueInter(1, 2), Inter(1, 2); - occ::handle> Poles1d = - new (NCollection_HArray2)(1, Ordre, 1, Ordre); + return; - // Calcul de B - Inter.SetValue(1, -1); - Inter.SetValue(2, 1); - TrueInter.SetValue(1, -1); - TrueInter.SetValue(2, 1); + // B is the monomial-to-BSpline conversion matrix on [-1,1], degree 6. + // It is a mathematical constant computed once. + static const math_Matrix THE_BASIS = []() { + constexpr int anOrdre = 7; + NCollection_Array1 aCoeffs(1, anOrdre * anOrdre); + NCollection_Array1 anInter(1, 2), aTrueInter(1, 2); + anInter.SetValue(1, -1); + anInter.SetValue(2, 1); + aTrueInter.SetValue(1, -1); + aTrueInter.SetValue(2, 1); + aCoeffs.Init(0); + for (int ii = 1; ii <= anOrdre; ii++) + aCoeffs.SetValue(ii + (ii - 1) * anOrdre, 1); - Coeffs.Init(0); - for (ii = 1; ii <= Ordre; ii++) - { - Coeffs.SetValue(ii + (ii - 1) * Ordre, 1); - } + Convert_CompPolynomialToPoles aConverter(anOrdre, + anOrdre - 1, + anOrdre - 1, + aCoeffs, + anInter, + aTrueInter); + const NCollection_Array2& aPoles = aConverter.Poles(); + math_Matrix aResult(1, anOrdre, 1, anOrdre); + for (int jj = 1; jj <= anOrdre; jj++) + for (int ii = 1; ii <= anOrdre; ii++) + { + double aTerm = aPoles.Value(ii, jj); + if (std::abs(aTerm - 1) < 1.e-9) + aTerm = 1; + if (std::abs(aTerm + 1) < 1.e-9) + aTerm = -1; + aResult(ii, jj) = aTerm; + } + return aResult; + }(); - // Convertion - Convert_CompPolynomialToPoles AConverter(Ordre, Ordre - 1, Ordre - 1, Coeffs, Inter, TrueInter); - Poles1d = new NCollection_HArray2(AConverter.Poles()); + B = THE_BASIS; - for (jj = 1; jj <= Ordre; jj++) - { - for (ii = 1; ii <= Ordre; ii++) - { - terme = Poles1d->Value(ii, jj); - if (std::abs(terme - 1) < 1.e-9) - terme = 1; // petite retouche - if (std::abs(terme + 1) < 1.e-9) - terme = -1; - B(ii, jj) = terme; - } - } - - // Init des polynomes Vx.Init(0); Vx(1) = 1; Vy.Init(0); diff --git a/src/ModelingData/TKGeomBase/Approx/Approx_Curve2d.cxx b/src/ModelingData/TKGeomBase/Approx/Approx_Curve2d.cxx index c1550bc2c0..88e94fa65b 100644 --- a/src/ModelingData/TKGeomBase/Approx/Approx_Curve2d.cxx +++ b/src/ModelingData/TKGeomBase/Approx/Approx_Curve2d.cxx @@ -149,12 +149,13 @@ Approx_Curve2d::Approx_Curve2d(const occ::handle& C2D, if (myHasResult) { - NCollection_Array1 Poles2d(1, aApprox.NbPoles()); - NCollection_Array1 Poles1dU(1, aApprox.NbPoles()); + const int aNbPoles = aApprox.NbPoles(); + NCollection_Array1 Poles2d(1, aNbPoles); + NCollection_Array1 Poles1dU(1, aNbPoles); aApprox.Poles1d(1, Poles1dU); - NCollection_Array1 Poles1dV(1, aApprox.NbPoles()); + NCollection_Array1 Poles1dV(1, aNbPoles); aApprox.Poles1d(2, Poles1dV); - for (int i = 1; i <= aApprox.NbPoles(); i++) + for (int i = 1; i <= aNbPoles; i++) Poles2d.SetValue(i, gp_Pnt2d(Poles1dU.Value(i), Poles1dV.Value(i))); occ::handle> Knots = aApprox.Knots(); diff --git a/src/ModelingData/TKGeomBase/Approx/Approx_CurveOnSurface.cxx b/src/ModelingData/TKGeomBase/Approx/Approx_CurveOnSurface.cxx index d86c7e8ca5..c6736467c2 100644 --- a/src/ModelingData/TKGeomBase/Approx/Approx_CurveOnSurface.cxx +++ b/src/ModelingData/TKGeomBase/Approx/Approx_CurveOnSurface.cxx @@ -502,21 +502,22 @@ void Approx_CurveOnSurface::Perform(const int theMaxSegments, occ::handle> Mults = aApprox.Multiplicities(); int Degree = aApprox.Degree(); + const int aNbPoles = aApprox.NbPoles(); if (!theOnly2d) { - NCollection_Array1 Poles(1, aApprox.NbPoles()); + NCollection_Array1 Poles(1, aNbPoles); aApprox.Poles(1, Poles); myCurve3d = new Geom_BSplineCurve(Poles, Knots->Array1(), Mults->Array1(), Degree); myError3d = aApprox.MaxError(3, 1); } if (!theOnly3d) { - NCollection_Array1 Poles2d(1, aApprox.NbPoles()); - NCollection_Array1 Poles1dU(1, aApprox.NbPoles()); + NCollection_Array1 Poles2d(1, aNbPoles); + NCollection_Array1 Poles1dU(1, aNbPoles); aApprox.Poles1d(1, Poles1dU); - NCollection_Array1 Poles1dV(1, aApprox.NbPoles()); + NCollection_Array1 Poles1dV(1, aNbPoles); aApprox.Poles1d(2, Poles1dV); - for (int i = 1; i <= aApprox.NbPoles(); i++) + for (int i = 1; i <= aNbPoles; i++) Poles2d.SetValue(i, gp_Pnt2d(Poles1dU.Value(i), Poles1dV.Value(i))); myCurve2d = new Geom2d_BSplineCurve(Poles2d, Knots->Array1(), Mults->Array1(), Degree); diff --git a/src/ModelingData/TKGeomBase/Approx/Approx_SweepApproximation.cxx b/src/ModelingData/TKGeomBase/Approx/Approx_SweepApproximation.cxx index e29a70c8f7..5cac28d076 100644 --- a/src/ModelingData/TKGeomBase/Approx/Approx_SweepApproximation.cxx +++ b/src/ModelingData/TKGeomBase/Approx/Approx_SweepApproximation.cxx @@ -303,12 +303,13 @@ void Approx_SweepApproximation::Approximation( // --> Fill Champs of the surface ---- int ii, jj; - vdeg = Approx.Degree(); + vdeg = Approx.Degree(); + const int aNbPoles = Approx.NbPoles(); // Unfortunately Adv_Approx stores the transposition of the required // so, writing tabPoles = Approx.Poles() will give an erroneous result // It is only possible to allocate and recopy term by term... - tabPoles = new (NCollection_HArray2)(1, Num3DSS, 1, Approx.NbPoles()); - tabWeights = new (NCollection_HArray2)(1, Num3DSS, 1, Approx.NbPoles()); + tabPoles = new (NCollection_HArray2)(1, Num3DSS, 1, aNbPoles); + tabWeights = new (NCollection_HArray2)(1, Num3DSS, 1, aNbPoles); if (Num1DSS == Num3DSS) { @@ -316,7 +317,7 @@ void Approx_SweepApproximation::Approximation( gp_Pnt P; for (ii = 1; ii <= Num3DSS; ii++) { - for (jj = 1; jj <= Approx.NbPoles(); jj++) + for (jj = 1; jj <= aNbPoles; jj++) { P = Approx.Poles()->Value(jj, ii); wpoid = Approx.Poles1d()->Value(jj, ii); @@ -332,7 +333,7 @@ void Approx_SweepApproximation::Approximation( tabWeights->Init(1); for (ii = 1; ii <= Num3DSS; ii++) { - for (jj = 1; jj <= Approx.NbPoles(); jj++) + for (jj = 1; jj <= aNbPoles; jj++) { tabPoles->SetValue(ii, jj, Approx.Poles()->Value(jj, ii)); } @@ -355,10 +356,10 @@ void Approx_SweepApproximation::Approximation( { TrsfInv = AAffin->Value(ii).Inverted(); occ::handle> P2d = - new (NCollection_HArray1)(1, Approx.NbPoles()); + new (NCollection_HArray1)(1, aNbPoles); Approx.Poles2d(ii, P2d->ChangeArray1()); // do not forget to apply inverted homothety. - for (jj = 1; jj <= Approx.NbPoles(); jj++) + for (jj = 1; jj <= aNbPoles; jj++) { TrsfInv.Transforms(P2d->ChangeValue(jj).ChangeCoord()); }