From 851733707559819f90efdd2ef07e4f649a9158c2 Mon Sep 17 00:00:00 2001 From: Francesco Lannutti Date: Mon, 20 Jun 2016 13:05:36 +0200 Subject: [PATCH] Fixed the sign inversion calculation for the determinant of KLU --- src/maths/KLU/klusmp.c | 42 ++++++++++++++++++++++++++++++++++-------- 1 file changed, 34 insertions(+), 8 deletions(-) diff --git a/src/maths/KLU/klusmp.c b/src/maths/KLU/klusmp.c index 82d25e3a5..b6648ab18 100644 --- a/src/maths/KLU/klusmp.c +++ b/src/maths/KLU/klusmp.c @@ -557,7 +557,7 @@ spDeterminant_KLU (SMPmatrix *Matrix, int *pExponent, RealNumber *pDeterminant, int *P, *Q ; double *Rs, *Ux, *Uz ; - unsigned int nSwap ; + unsigned int nSwap, nSwapP, nSwapQ ; #define NORM(a) (nr = ABS((a).Real), ni = ABS((a).Imag), MAX (nr,ni)) @@ -631,14 +631,27 @@ spDeterminant_KLU (SMPmatrix *Matrix, int *pExponent, RealNumber *pDeterminant, printf ("U - Value: %-.9g\t%-.9g\n", Ux [I], Uz [I]) ; } */ - nSwap = 0 ; + nSwapP = 0 ; for (I = 0 ; I < Matrix->CKTkluN ; I++) { - if (P [I] != I || Q [I] != I) + if (P [I] != I) { - nSwap++ ; + nSwapP++ ; } } + nSwapP /= 2 ; + + nSwapQ = 0 ; + for (I = 0 ; I < Matrix->CKTkluN ; I++) + { + if (Q [I] != I) + { + nSwapQ++ ; + } + } + nSwapQ /= 2 ; + + nSwap = nSwapP + nSwapQ ; /* free (Lp) ; free (Li) ; @@ -722,18 +735,31 @@ spDeterminant_KLU (SMPmatrix *Matrix, int *pExponent, RealNumber *pDeterminant, klu_extract_Udiag (Matrix->CKTkluNumeric, Matrix->CKTkluSymbolic, Ux, P, Q, Rs, Matrix->CKTkluCommon) ; - nSwap = 0 ; + nSwapP = 0 ; for (I = 0 ; I < Matrix->CKTkluN ; I++) { - if (P [I] != I || Q [I] != I) + if (P [I] != I) { - nSwap++ ; + nSwapP++ ; } } + nSwapP /= 2 ; + + nSwapQ = 0 ; + for (I = 0 ; I < Matrix->CKTkluN ; I++) + { + if (Q [I] != I) + { + nSwapQ++ ; + } + } + nSwapQ /= 2 ; + + nSwap = nSwapP + nSwapQ ; while (I < Size) { - *pDeterminant /= Ux [I] ; + *pDeterminant /= (Ux [I] * Rs [I]) ; /* Scale Determinant. */ if (*pDeterminant != 0.0)