diff options
| author | Dennis Francis <dennisfrancis.in@gmail.com> | 2015-10-23 03:28:37 +0530 |
|---|---|---|
| committer | Eike Rathke <erack@redhat.com> | 2015-10-24 14:33:25 +0000 |
| commit | e1e73a97b20af862f2fb914cb16a4f74c3ad31cb (patch) | |
| tree | aabadeb65c3b5196f6b82fc6dc1207a680c32269 /sc | |
| parent | 5f64aba978b117eb4372fef913823243a800fbea (diff) | |
tdf#32834 : improve the precision of MDETERM calculation
MDETERM uses lcl_LUP_decompose() which as the name suggests does LUP decomposition.
This patch allows additive cancellations to occur in a cleaner way while doing the *only*
additive operation in LUP decomposition, that is while computing Shur complement.
This patch does not change the high level semantics of the algorithm.
Also note that this change makes improvement only in the case where matrix elements entered
by the user are *integers*. The change allows MDETERM to evaluate to exact 0.0 for
singular integer matrices.
The steps to calculate Shur complement are :
for i = k+1 to n
aik = aik / akk;
for j = k+1 to n
aij = aij - akj*aik
This is now modified as :
for i = k+1 to n
for j = k+1 to n
aij = ( aij*akk - akj*aik ) / akk
Without this change MDETERM() for certain singular matrices used to evaluate to a tiny non zero value,
which also caused MINVERSE() to generate a wrong output.
An example of such a matrix is :
1 2 3
4 5 6
7 8 9
Change-Id: Idd4211ddceab1b758fd05bfd57f7eecd5d4fd1a0
Reviewed-on: https://gerrit.libreoffice.org/19534
Reviewed-by: Eike Rathke <erack@redhat.com>
Tested-by: Eike Rathke <erack@redhat.com>
Diffstat (limited to 'sc')
| -rw-r--r-- | sc/source/core/tool/interpr5.cxx | 9 |
1 files changed, 5 insertions, 4 deletions
diff --git a/sc/source/core/tool/interpr5.cxx b/sc/source/core/tool/interpr5.cxx index 88cd12dd9a5c..eac7cd46b2ab 100644 --- a/sc/source/core/tool/interpr5.cxx +++ b/sc/source/core/tool/interpr5.cxx @@ -711,11 +711,12 @@ static int lcl_LUP_decompose( ScMatrix* mA, const SCSIZE n, // Compute Schur complement. for (SCSIZE i = k+1; i < n; ++i) { - double fTmp = mA->GetDouble( k, i) / mA->GetDouble( k, k); - mA->PutDouble( fTmp, k, i); + double fNum = mA->GetDouble( k, i); + double fDen = mA->GetDouble( k, k); + mA->PutDouble( fNum/fDen, k, i); for (SCSIZE j = k+1; j < n; ++j) - mA->PutDouble( mA->GetDouble( j, i) - fTmp * mA->GetDouble( j, - k), j, i); + mA->PutDouble( ( mA->GetDouble( j, i) * fDen - + fNum * mA->GetDouble( j, k) ) / fDen, j, i); } } #if OSL_DEBUG_LEVEL > 1 |