From ea4e6a6d806f3ecf47839774b08143a958c871c5 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Tue, 11 Oct 2011 21:55:18 +0200 Subject: [PATCH] Sparse Linear Algebra: fix subtraction bug gh-20 gh-18 --- .../LinearAlgebra/Complex/SparseVector.cs | 69 ++++++++++++++++++- .../LinearAlgebra/Complex32/SparseVector.cs | 69 ++++++++++++++++++- .../LinearAlgebra/Double/SparseVector.cs | 69 ++++++++++++++++++- .../LinearAlgebra/Single/SparseVector.cs | 69 ++++++++++++++++++- 4 files changed, 268 insertions(+), 8 deletions(-) diff --git a/src/Numerics/LinearAlgebra/Complex/SparseVector.cs b/src/Numerics/LinearAlgebra/Complex/SparseVector.cs index 38edfa32..4e33e2a0 100644 --- a/src/Numerics/LinearAlgebra/Complex/SparseVector.cs +++ b/src/Numerics/LinearAlgebra/Complex/SparseVector.cs @@ -550,9 +550,74 @@ namespace MathNet.Numerics.LinearAlgebra.Complex return; } - for (var index = 0; index < Count; index++) + var otherSparse = other as SparseVector; + if (otherSparse == null) + { + base.DoSubtract(other, result); + return; + } + + var resultSparse = result as SparseVector; + if (resultSparse == null) + { + base.DoSubtract(other, result); + return; + } + + // TODO (ruegg, 2011-10-11): Options to optimize? + + if (ReferenceEquals(this, resultSparse)) + { + int i = 0, j = 0; + while (i < NonZerosCount || j < otherSparse.NonZerosCount) + { + if (i < NonZerosCount && j < otherSparse.NonZerosCount && _nonZeroIndices[i] == otherSparse._nonZeroIndices[j]) + { + _nonZeroValues[i++] -= otherSparse._nonZeroValues[j++]; + } + else if (j >= otherSparse.NonZerosCount || i < NonZerosCount && _nonZeroIndices[i] < otherSparse._nonZeroIndices[j]) + { + _nonZeroValues[i] -= otherSparse.At(_nonZeroIndices[i]); + i++; + } + else + { + var otherValue = otherSparse._nonZeroValues[j]; + if (otherValue != Complex.Zero) + { + InsertAtUnchecked(i++, otherSparse._nonZeroIndices[j], -otherValue); + } + j++; + } + } + } + else { - result.At(index, At(index) - other.At(index)); + result.Clear(); + int i = 0, j = 0, last = -1; + while (i < NonZerosCount || j < otherSparse.NonZerosCount) + { + if (j >= otherSparse.NonZerosCount || i < NonZerosCount && _nonZeroIndices[i] <= otherSparse._nonZeroIndices[j]) + { + var next = _nonZeroIndices[i]; + if (next != last) + { + last = next; + result.At(next, _nonZeroValues[i] - otherSparse.At(next)); + } + i++; + } + else + { + var next = otherSparse._nonZeroIndices[j]; + if (next != last) + { + last = next; + result.At(next, At(next) - otherSparse._nonZeroValues[j]); + } + j++; + } + } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/SparseVector.cs b/src/Numerics/LinearAlgebra/Complex32/SparseVector.cs index 3b82a5a1..a58040fe 100644 --- a/src/Numerics/LinearAlgebra/Complex32/SparseVector.cs +++ b/src/Numerics/LinearAlgebra/Complex32/SparseVector.cs @@ -580,9 +580,74 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 return; } - for (var index = 0; index < Count; index++) + var otherSparse = other as SparseVector; + if (otherSparse == null) + { + base.DoSubtract(other, result); + return; + } + + var resultSparse = result as SparseVector; + if (resultSparse == null) + { + base.DoSubtract(other, result); + return; + } + + // TODO (ruegg, 2011-10-11): Options to optimize? + + if (ReferenceEquals(this, resultSparse)) + { + int i = 0, j = 0; + while (i < NonZerosCount || j < otherSparse.NonZerosCount) + { + if (i < NonZerosCount && j < otherSparse.NonZerosCount && _nonZeroIndices[i] == otherSparse._nonZeroIndices[j]) + { + _nonZeroValues[i++] -= otherSparse._nonZeroValues[j++]; + } + else if (j >= otherSparse.NonZerosCount || i < NonZerosCount && _nonZeroIndices[i] < otherSparse._nonZeroIndices[j]) + { + _nonZeroValues[i] -= otherSparse.At(_nonZeroIndices[i]); + i++; + } + else + { + var otherValue = otherSparse._nonZeroValues[j]; + if (otherValue != Complex32.Zero) + { + InsertAtUnchecked(i++, otherSparse._nonZeroIndices[j], -otherValue); + } + j++; + } + } + } + else { - result.At(index, At(index) - other.At(index)); + result.Clear(); + int i = 0, j = 0, last = -1; + while (i < NonZerosCount || j < otherSparse.NonZerosCount) + { + if (j >= otherSparse.NonZerosCount || i < NonZerosCount && _nonZeroIndices[i] <= otherSparse._nonZeroIndices[j]) + { + var next = _nonZeroIndices[i]; + if (next != last) + { + last = next; + result.At(next, _nonZeroValues[i] - otherSparse.At(next)); + } + i++; + } + else + { + var next = otherSparse._nonZeroIndices[j]; + if (next != last) + { + last = next; + result.At(next, At(next) - otherSparse._nonZeroValues[j]); + } + j++; + } + } } } diff --git a/src/Numerics/LinearAlgebra/Double/SparseVector.cs b/src/Numerics/LinearAlgebra/Double/SparseVector.cs index 4d516bfe..cffed1c5 100644 --- a/src/Numerics/LinearAlgebra/Double/SparseVector.cs +++ b/src/Numerics/LinearAlgebra/Double/SparseVector.cs @@ -498,9 +498,74 @@ namespace MathNet.Numerics.LinearAlgebra.Double return; } - for (var index = 0; index < Count; index++) + var otherSparse = other as SparseVector; + if (otherSparse == null) + { + base.DoSubtract(other, result); + return; + } + + var resultSparse = result as SparseVector; + if (resultSparse == null) + { + base.DoSubtract(other, result); + return; + } + + // TODO (ruegg, 2011-10-11): Options to optimize? + + if (ReferenceEquals(this, resultSparse)) + { + int i = 0, j = 0; + while (i < NonZerosCount || j < otherSparse.NonZerosCount) + { + if (i < NonZerosCount && j < otherSparse.NonZerosCount && _nonZeroIndices[i] == otherSparse._nonZeroIndices[j]) + { + _nonZeroValues[i++] -= otherSparse._nonZeroValues[j++]; + } + else if (j >= otherSparse.NonZerosCount || i < NonZerosCount && _nonZeroIndices[i] < otherSparse._nonZeroIndices[j]) + { + _nonZeroValues[i] -= otherSparse.At(_nonZeroIndices[i]); + i++; + } + else + { + var otherValue = otherSparse._nonZeroValues[j]; + if (otherValue != 0.0) + { + InsertAtUnchecked(i++, otherSparse._nonZeroIndices[j], -otherValue); + } + j++; + } + } + } + else { - result.At(index, At(index) - other.At(index)); + result.Clear(); + int i = 0, j = 0, last = -1; + while (i < NonZerosCount || j < otherSparse.NonZerosCount) + { + if (j >= otherSparse.NonZerosCount || i < NonZerosCount && _nonZeroIndices[i] <= otherSparse._nonZeroIndices[j]) + { + var next = _nonZeroIndices[i]; + if (next != last) + { + last = next; + result.At(next, _nonZeroValues[i] - otherSparse.At(next)); + } + i++; + } + else + { + var next = otherSparse._nonZeroIndices[j]; + if (next != last) + { + last = next; + result.At(next, At(next) - otherSparse._nonZeroValues[j]); + } + j++; + } + } } } diff --git a/src/Numerics/LinearAlgebra/Single/SparseVector.cs b/src/Numerics/LinearAlgebra/Single/SparseVector.cs index 1c8ce6f8..b1c48feb 100644 --- a/src/Numerics/LinearAlgebra/Single/SparseVector.cs +++ b/src/Numerics/LinearAlgebra/Single/SparseVector.cs @@ -528,9 +528,74 @@ namespace MathNet.Numerics.LinearAlgebra.Single return; } - for (var index = 0; index < Count; index++) + var otherSparse = other as SparseVector; + if (otherSparse == null) + { + base.DoSubtract(other, result); + return; + } + + var resultSparse = result as SparseVector; + if (resultSparse == null) + { + base.DoSubtract(other, result); + return; + } + + // TODO (ruegg, 2011-10-11): Options to optimize? + + if (ReferenceEquals(this, resultSparse)) + { + int i = 0, j = 0; + while (i < NonZerosCount || j < otherSparse.NonZerosCount) + { + if (i < NonZerosCount && j < otherSparse.NonZerosCount && _nonZeroIndices[i] == otherSparse._nonZeroIndices[j]) + { + _nonZeroValues[i++] -= otherSparse._nonZeroValues[j++]; + } + else if (j >= otherSparse.NonZerosCount || i < NonZerosCount && _nonZeroIndices[i] < otherSparse._nonZeroIndices[j]) + { + _nonZeroValues[i] -= otherSparse.At(_nonZeroIndices[i]); + i++; + } + else + { + var otherValue = otherSparse._nonZeroValues[j]; + if (otherValue != 0.0) + { + InsertAtUnchecked(i++, otherSparse._nonZeroIndices[j], -otherValue); + } + j++; + } + } + } + else { - result.At(index, At(index) - other.At(index)); + result.Clear(); + int i = 0, j = 0, last = -1; + while (i < NonZerosCount || j < otherSparse.NonZerosCount) + { + if (j >= otherSparse.NonZerosCount || i < NonZerosCount && _nonZeroIndices[i] <= otherSparse._nonZeroIndices[j]) + { + var next = _nonZeroIndices[i]; + if (next != last) + { + last = next; + result.At(next, _nonZeroValues[i] - otherSparse.At(next)); + } + i++; + } + else + { + var next = otherSparse._nonZeroIndices[j]; + if (next != last) + { + last = next; + result.At(next, At(next) - otherSparse._nonZeroValues[j]); + } + j++; + } + } } }