Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

Sparse Linear Algebra: fix subtraction bug gh-20 gh-18

  • Loading branch information...
commit e0c69362cb04531fe66baae4740f16fd17304802 1 parent be4d769
Christoph Ruegg cdrnet authored
69 src/Numerics/LinearAlgebra/Complex/SparseVector.cs
View
@@ -550,9 +550,74 @@ protected override void DoSubtract(Vector<Complex> other, Vector<Complex> result
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++;
+ }
+ }
}
}
69 src/Numerics/LinearAlgebra/Complex32/SparseVector.cs
View
@@ -580,9 +580,74 @@ protected override void DoSubtract(Vector<Complex32> other, Vector<Complex32> re
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++;
+ }
+ }
}
}
69 src/Numerics/LinearAlgebra/Double/SparseVector.cs
View
@@ -498,9 +498,74 @@ protected override void DoSubtract(Vector<double> other, Vector<double> result)
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++;
+ }
+ }
}
}
69 src/Numerics/LinearAlgebra/Single/SparseVector.cs
View
@@ -528,9 +528,74 @@ protected override void DoSubtract(Vector<float> other, Vector<float> result)
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++;
+ }
+ }
}
}
Please sign in to comment.
Something went wrong with that request. Please try again.