sc/source/core/tool/interpr5.cxx | 50 +++++++++++++++++++--------------------
1 file changed, 25 insertions(+), 25 deletions(-)
New commits:
commit 2b90807fd474cb4e26c80ba5407aca4b66f6d3a2
Author: dante <dante19031...@gmail.com>
AuthorDate: Wed May 5 11:06:54 2021 +0200
Commit: Mike Kaganski <mike.kagan...@collabora.com>
CommitDate: Fri May 7 00:22:18 2021 +0200
tdf#137679 Use Kahan summation for interpr5.cxx
Change-Id: Iefd09f2af4c77a464e47367c75b9e7e303f2c2a9
Reviewed-on: https://gerrit.libreoffice.org/c/core/+/115128
Tested-by: Jenkins
Reviewed-by: Mike Kaganski <mike.kagan...@collabora.com>
diff --git a/sc/source/core/tool/interpr5.cxx b/sc/source/core/tool/interpr5.cxx
index 4f10bb8961d2..f7dd527780ea 100644
--- a/sc/source/core/tool/interpr5.cxx
+++ b/sc/source/core/tool/interpr5.cxx
@@ -90,15 +90,14 @@ struct MatrixPow
void lcl_MFastMult(const ScMatrixRef& pA, const ScMatrixRef& pB, const
ScMatrixRef& pR,
SCSIZE n, SCSIZE m, SCSIZE l)
{
- double sum;
for (SCSIZE row = 0; row < n; row++)
{
for (SCSIZE col = 0; col < l; col++)
{ // result element(col, row) =sum[ (row of A) * (column of B)]
- sum = 0.0;
+ KahanSum fSum = 0.0;
for (SCSIZE k = 0; k < m; k++)
- sum += pA->GetDouble(k,row) * pB->GetDouble(col,k);
- pR->PutDouble(sum, col, row);
+ fSum += pA->GetDouble(k,row) * pB->GetDouble(col,k);
+ pR->PutDouble(fSum.get(), col, row);
}
}
}
@@ -886,7 +885,7 @@ static void lcl_LUP_solve( const ScMatrix* mLU, const
SCSIZE n,
// Define y=Ux and solve for y in Ly=Pb using forward substitution.
for (SCSIZE i=0; i < n; ++i)
{
- double fSum = B[P[i]];
+ KahanSum fSum = B[P[i]];
// Matrix inversion comes with a lot of zeros in the B vectors, we
// don't have to do all the computing with results multiplied by zero.
// Until then, simply lookout for the position of the first nonzero
@@ -894,19 +893,19 @@ static void lcl_LUP_solve( const ScMatrix* mLU, const
SCSIZE n,
if (nFirst != SCSIZE_MAX)
{
for (SCSIZE j = nFirst; j < i; ++j)
- fSum -= mLU->GetDouble( j, i) * X[j]; // X[j] === y[j]
+ fSum -= mLU->GetDouble( j, i) * X[j]; // X[j] === y[j]
}
- else if (fSum)
+ else if (fSum != 0)
nFirst = i;
- X[i] = fSum; // X[i] === y[i]
+ X[i] = fSum.get(); // X[i] === y[i]
}
// Solve for x in Ux=y using back substitution.
for (SCSIZE i = n; i--; )
{
- double fSum = X[i]; // X[i] === y[i]
+ KahanSum fSum = X[i]; // X[i] === y[i]
for (SCSIZE j = i+1; j < n; ++j)
- fSum -= mLU->GetDouble( j, i) * X[j]; // X[j] === x[j]
- X[i] = fSum / mLU->GetDouble( i, i); // X[i] === x[i]
+ fSum -= mLU->GetDouble( j, i) * X[j]; // X[j] === x[j]
+ X[i] = fSum.get() / mLU->GetDouble( i, i); // X[i] === x[i]
}
#ifdef DEBUG_SC_LUP_DECOMPOSITION
fprintf( stderr, "\n%s\n", "lcl_LUP_solve():");
@@ -1098,17 +1097,17 @@ void ScInterpreter::ScMatMult()
pRMat = GetNewMat(nC2, nR1, /*bEmpty*/true);
if (pRMat)
{
- double sum;
+ KahanSum fSum;
for (SCSIZE i = 0; i < nR1; i++)
{
for (SCSIZE j = 0; j < nC2; j++)
{
- sum = 0.0;
+ fSum = 0.0;
for (SCSIZE k = 0; k < nC1; k++)
{
- sum +=
pMat1->GetDouble(k,i)*pMat2->GetDouble(j,k);
+ fSum +=
pMat1->GetDouble(k,i)*pMat2->GetDouble(j,k);
}
- pRMat->PutDouble(sum, j, i);
+ pRMat->PutDouble(fSum.get(), j, i);
}
}
PushMatrix(pRMat);
@@ -1804,7 +1803,8 @@ void ScInterpreter::CalculateSumX2MY2SumX2DY2(bool
_bSumX2DY2)
PushNoValue();
return;
}
- double fVal, fSum = 0.0;
+ double fVal;
+ KahanSum fSum = 0.0;
for (i = 0; i < nC1; i++)
for (j = 0; j < nR1; j++)
if (!pMat1->IsStringOrEmpty(i,j) && !pMat2->IsStringOrEmpty(i,j))
@@ -1817,7 +1817,7 @@ void ScInterpreter::CalculateSumX2MY2SumX2DY2(bool
_bSumX2DY2)
else
fSum -= fVal * fVal;
}
- PushDouble(fSum);
+ PushDouble(fSum.get());
}
void ScInterpreter::ScSumX2DY2()
@@ -2701,7 +2701,7 @@ void ScInterpreter::CalculateRGPRKP(bool _bRKP)
// = RMSE * sqrt( Xmean * (R' R)^(-1) * Xmean' + 1/N)
// (R' R)^(-1) = R^(-1) * (R')^(-1). Do not calculate it as
// a whole matrix, but iterate over unit vectors.
- double fSigmaIntercept = 0.0;
+ KahanSum aSigmaIntercept = 0.0;
double fPart; // for Xmean * single column of (R' R)^(-1)
for (SCSIZE col = 0; col < K; col++)
{
@@ -2719,13 +2719,13 @@ void ScInterpreter::CalculateRGPRKP(bool _bRKP)
if (bConstant)
{
fPart = lcl_GetSumProduct(pMeans, pMatZ, K);
- fSigmaIntercept += fPart * pMeans->GetDouble(col);
+ aSigmaIntercept += fPart * pMeans->GetDouble(col);
}
}
if (bConstant)
{
- fSigmaIntercept = fRMSE
- * sqrt(fSigmaIntercept + 1.0 /
static_cast<double>(N));
+ double fSigmaIntercept = fRMSE
+ * sqrt( (aSigmaIntercept + 1.0 /
static_cast<double>(N) ).get() );
pResMat->PutDouble(fSigmaIntercept, K, 1);
}
else
@@ -2859,7 +2859,7 @@ void ScInterpreter::CalculateRGPRKP(bool _bRKP)
// (R' R)^(-1) = R^(-1) * (R')^(-1). Do not calculate it as
// a whole matrix, but iterate over unit vectors.
// (R' R) ^(-1) is symmetric
- double fSigmaIntercept = 0.0;
+ KahanSum aSigmaIntercept = 0.0;
double fPart; // for Xmean * single col of (R' R)^(-1)
for (SCSIZE row = 0; row < K; row++)
{
@@ -2876,13 +2876,13 @@ void ScInterpreter::CalculateRGPRKP(bool _bRKP)
if (bConstant)
{
fPart = lcl_GetSumProduct(pMeans, pMatZ, K);
- fSigmaIntercept += fPart * pMeans->GetDouble(row);
+ aSigmaIntercept += fPart * pMeans->GetDouble(row);
}
}
if (bConstant)
{
- fSigmaIntercept = fRMSE
- * sqrt(fSigmaIntercept + 1.0 /
static_cast<double>(N));
+ double fSigmaIntercept = fRMSE
+ * sqrt( (aSigmaIntercept + 1.0 /
static_cast<double>(N) ).get() );
pResMat->PutDouble(fSigmaIntercept, K, 1);
}
else
_______________________________________________
Libreoffice-commits mailing list
libreoffice-comm...@lists.freedesktop.org
https://lists.freedesktop.org/mailman/listinfo/libreoffice-commits
