Browse Source

sparse vector optimizations

la-knuth
Marcus Cuda 16 years ago
parent
commit
9b72b71ebf
  1. 13
      src/MathNet.Numerics.5.1.ReSharper
  2. 16
      src/Numerics/LinearAlgebra/Complex/DenseVector.cs
  3. 2
      src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
  4. 576
      src/Numerics/LinearAlgebra/Complex/SparseVector.cs
  5. 101
      src/Numerics/LinearAlgebra/Complex/Vector.cs
  6. 16
      src/Numerics/LinearAlgebra/Complex32/DenseVector.cs
  7. 2
      src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
  8. 548
      src/Numerics/LinearAlgebra/Complex32/SparseVector.cs
  9. 99
      src/Numerics/LinearAlgebra/Complex32/Vector.cs
  10. 16
      src/Numerics/LinearAlgebra/Double/DenseVector.cs
  11. 24
      src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
  12. 566
      src/Numerics/LinearAlgebra/Double/SparseVector.cs
  13. 96
      src/Numerics/LinearAlgebra/Double/Vector.cs
  14. 36
      src/Numerics/LinearAlgebra/Generic/Vector.cs
  15. 16
      src/Numerics/LinearAlgebra/Single/DenseVector.cs
  16. 2
      src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
  17. 539
      src/Numerics/LinearAlgebra/Single/SparseVector.cs
  18. 98
      src/Numerics/LinearAlgebra/Single/Vector.cs
  19. 22
      src/UnitTests/LinearAlgebraTests/Complex/UserDefinedVectorTests.cs
  20. 23
      src/UnitTests/LinearAlgebraTests/Complex32/UserDefinedVectorTests.cs
  21. 23
      src/UnitTests/LinearAlgebraTests/Double/UserDefinedVectorTests.cs
  22. 23
      src/UnitTests/LinearAlgebraTests/Single/UserDefinedVectorTests.cs

13
src/MathNet.Numerics.5.1.ReSharper

@ -15,7 +15,10 @@ cancelled
pre
preconditioner
Preconditioners
indetermined</UserWords>
indetermined
indices
Pointwise
Frobenius</UserWords>
</CustomDictionary>
</Dictionaries>
</CustomDictionaries>
@ -561,7 +564,7 @@ indetermined</UserWords>
<MustNotHaveSuffixes />
</NamingConventionRule>
<NamingConventionRule>
<IsDisabled>false</IsDisabled>
<IsDisabled>true</IsDisabled>
<Matches>
<Match>
<AccessLevel>Any</AccessLevel>
@ -593,7 +596,7 @@ indetermined</UserWords>
<MustNotHaveSuffixes />
</NamingConventionRule>
<NamingConventionRule>
<IsDisabled>false</IsDisabled>
<IsDisabled>true</IsDisabled>
<Matches>
<Match>
<AccessLevel>Any</AccessLevel>
@ -657,7 +660,7 @@ indetermined</UserWords>
<MustNotHaveSuffixes />
</NamingConventionRule>
<NamingConventionRule>
<IsDisabled>false</IsDisabled>
<IsDisabled>true</IsDisabled>
<Matches>
<Match>
<AccessLevel>Any</AccessLevel>
@ -686,7 +689,7 @@ indetermined</UserWords>
<MustNotHaveSuffixes />
</NamingConventionRule>
<NamingConventionRule>
<IsDisabled>false</IsDisabled>
<IsDisabled>true</IsDisabled>
<Matches>
<Match>
<AccessLevel>Any</AccessLevel>

16
src/Numerics/LinearAlgebra/Complex/DenseVector.cs

@ -1334,5 +1334,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
index => otherVector.Data[index] = Data[index].Conjugate());
}
}
/// <summary>Gets the value at the given <paramref name="index"/>.</summary>
/// <param name="index">The index of the value to get or set.</param>
/// <returns>The value of the vector at the given <paramref name="index"/>.</returns>
internal protected override Complex At(int index)
{
return Data[index];
}
/// <summary>Sets the <paramref name="value"/> at the given <paramref name="index"/>.</summary>
/// <param name="index">The index of the value to get or set.</param>
/// <param name="value">The value to set.</param>
internal protected override void At(int index, Complex value)
{
Data[index] = value;
}
}
}

2
src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs

@ -194,7 +194,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// </summary>
/// <param name="matrix">The matrix to copy.</param>
public SparseMatrix(Matrix<Complex> matrix)
: base(matrix.RowCount, matrix.ColumnCount)
: this(matrix.RowCount, matrix.ColumnCount)
{
var sparseMatrix = matrix as SparseMatrix;

576
src/Numerics/LinearAlgebra/Complex/SparseVector.cs

@ -38,13 +38,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <summary>
/// A vector with sparse storage.
/// </summary>
/// <remarks>The sparse vector is not thread safe.</remarks>
public class SparseVector : Vector
{
/// <summary>
/// Lock object for the indexer.
/// </summary>
private readonly object _lockObject = new object();
/// <summary>
/// Gets the vector's internal data. The array containing the actual values; only the non-zero values are stored.
/// </summary>
@ -130,7 +126,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
{
for (var i = 0; i < other.Count; i++)
{
this[i] = other[i];
this[i] = other.At(i);
}
}
else
@ -174,9 +170,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// </summary>
/// <param name="array">The array to create this vector from.</param>
/// <remarks>The vector copy the array. Any changes to the vector will NOT change the array.</remarks>
public SparseVector(Complex[] array) : this(array.Length)
public SparseVector(IList<Complex> array) : this(array.Count)
{
for (var i = 0; i < array.Length; i++)
for (var i = 0; i < array.Count; i++)
{
this[i] = array[i];
}
@ -193,7 +189,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
var matrix = new SparseMatrix(Count, 1);
for (var i = 0; i < NonZerosCount; i++)
{
matrix[_nonZeroIndices[i], 0] = _nonZeroValues[i];
matrix.At(_nonZeroIndices[i], 0, _nonZeroValues[i]);
}
return matrix;
@ -208,55 +204,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
var matrix = new SparseMatrix(1, Count);
for (var i = 0; i < NonZerosCount; i++)
{
matrix[0, _nonZeroIndices[i]] = _nonZeroValues[i];
matrix.At(0, _nonZeroIndices[i], _nonZeroValues[i]);
}
return matrix;
}
/// <summary>Gets or sets the value at the given <paramref name="index"/>.</summary>
/// <param name="index">The index of the value to get or set.</param>
/// <returns>The value of the vector at the given <paramref name="index"/>.</returns>
/// <exception cref="IndexOutOfRangeException">If <paramref name="index"/> is negative or
/// greater than the size of the vector.</exception>
public override Complex this[int index]
{
get
{
// If index is out of bounds
if ((index < 0) || (index >= Count))
{
throw new IndexOutOfRangeException();
}
lock (_lockObject)
{
// Search if item idex exists in NonZeroIndices array in range "0 - complex nonzero values count"
var itemIndex = Array.BinarySearch(_nonZeroIndices, 0, NonZerosCount, index);
if (itemIndex >= 0)
{
return _nonZeroValues[itemIndex];
}
}
return Complex.Zero;
}
set
{
// If index is out of bounds
if ((index < 0) || (index >= Count))
{
throw new IndexOutOfRangeException();
}
lock (_lockObject)
{
SetValue(index, value);
}
}
}
/// <summary>
/// Creates a matrix with the given dimensions using the same storage type
/// as this vector.
@ -330,10 +283,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
var otherVector = target as SparseVector;
if (otherVector == null)
{
CommonParallel.For(
0,
Count,
index => target[index] = this[index]);
target.Clear();
for (var index = 0; index < NonZerosCount; index++)
{
target.At(_nonZeroIndices[index], _nonZeroValues[index]);
}
}
else
{
@ -392,147 +347,42 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
}
}
}
#region Operators and supplementary functions
/// <summary>
/// Adds a complex to each element of the vector.
/// </summary>
/// <param name="complex">The complex to add.</param>
/// <returns>A copy of the vector with the complex added.</returns>
public override Vector<Complex> Add(Complex complex)
{
if (complex == Complex.Zero)
{
return Clone();
}
var copy = (SparseVector)Clone();
for (var i = 0; i < Count; i++)
{
copy[i] += complex;
}
return copy;
}
/// <summary>
/// Adds a complex to each element of the vector and stores the result in the result vector.
/// </summary>
/// <param name="complex">The complex to add.</param>
/// <param name="result">The vector to store the result of the addition.</param>
/// <exception cref="ArgumentNullException">If the result vector is <see langword="null" />.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
public override void Add(Complex complex, Vector<Complex> result)
{
if (result == null)
{
throw new ArgumentNullException("result");
}
if (Count != result.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
}
base.Add(complex, result);
}
/// <summary>
/// Adds another vector to this vector.
/// Adds a scalar to each element of the vector and stores the result in the result vector.
/// </summary>
/// <param name="other">The vector to add to this one.</param>
/// <returns>A new vector containing the sum of both vectors.</returns>
/// <exception cref="ArgumentNullException">If the other vector is <see langword="null" />.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
public override Vector<Complex> Add(Vector<Complex> other)
{
if (other == null)
{
throw new ArgumentNullException("other");
}
if (Count != other.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
}
var sparseVector = other as SparseVector;
if (sparseVector == null)
{
return base.Add(other);
}
var copy = (SparseVector)Clone();
copy.AddScaledSparseVector(Complex.One, sparseVector);
return copy;
}
/// <summary>
/// Adds the scaled sparse vector.
/// </summary>
/// <param name="alpha">The alpha.</param>
/// <param name="other">The other.</param>
private void AddScaledSparseVector(Complex alpha, SparseVector other)
/// <param name="scalar">
/// The scalar to add.
/// </param>
/// <param name="result">
/// The vector to store the result of the addition.
/// </param>
protected override void DoAdd(Complex scalar, Vector<Complex> result)
{
if (other == null)
{
throw new ArgumentNullException("other");
}
if (Count != other.Count)
if (scalar == Complex.Zero)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
}
if (!ReferenceEquals(this, result))
{
CopyTo(result);
}
if (alpha == Complex.Zero)
{
return;
}
// I don't use ILinearAlgebraProvider because we will get no benefit due to "lock" in this[index]
// Possible fucniton in ILinearAlgebraProvider may be AddSparseVectorToScaledSparseVector(T[] y, int[] yIndices, T alpha, T[] x, int[] xIndices);
// But it require to develop value setting algorithm and due to "lock" it will be even more greedy then implemented below
if (ReferenceEquals(this, other))
if (ReferenceEquals(this, result))
{
// Adding the same instance of sparse vector. That means if we modify "this" then "other" will be modified too.
// To avoid such problem lets change values in internal storage of "this"
if (alpha == Complex.One)
{
for (var i = 0; i < NonZerosCount; i++)
{
_nonZeroValues[i] += _nonZeroValues[i];
}
}
else if (alpha == -Complex.One)
{
Clear(); // Vector is subtracted from itself
return;
}
else
{
for (var i = 0; i < NonZerosCount; i++)
{
_nonZeroValues[i] += alpha * _nonZeroValues[i];
}
}
CommonParallel.For(
0,
NonZerosCount,
index => _nonZeroValues[index] += scalar);
}
else
{
// "this" and "other" are different objects, so by modifying "this" the "other" object will not be changed
if (alpha == Complex.One)
for (var index = 0; index < Count; index++)
{
for (var i = 0; i < other.NonZerosCount; i++)
{
this[other._nonZeroIndices[i]] += other._nonZeroValues[i];
}
}
else
{
for (var i = 0; i < other.NonZerosCount; i++)
{
this[other._nonZeroIndices[i]] += alpha * other._nonZeroValues[i];
}
result.At(index, At(index) + scalar);
}
}
}
@ -540,53 +390,26 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <summary>
/// Adds another vector to this vector and stores the result into the result vector.
/// </summary>
/// <param name="other">The vector to add to this one.</param>
/// <param name="result">The vector to store the result of the addition.</param>
/// <exception cref="ArgumentNullException">If the other vector is <see langword="null" />.</exception>
/// <exception cref="ArgumentNullException">If the result vector is <see langword="null" />.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
public override void Add(Vector<Complex> other, Vector<Complex> result)
/// <param name="other">
/// The vector to add to this one.
/// </param>
/// <param name="result">
/// The vector to store the result of the addition.
/// </param>
protected override void DoAdd(Vector<Complex> other, Vector<Complex> result)
{
if (result == null)
{
throw new ArgumentNullException("result");
}
if (Count != other.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
}
if (Count != result.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
}
if (ReferenceEquals(this, result) || ReferenceEquals(other, result))
if (ReferenceEquals(this, result))
{
var tmp = Add(other);
tmp.CopyTo(result);
CommonParallel.For(
0,
NonZerosCount,
index => _nonZeroValues[index] += _nonZeroValues[index]);
}
else
{
var sparse = result as SparseVector;
if (sparse == null)
for (var index = 0; index < Count; index++)
{
base.Add(other, result);
}
else
{
var sparseother = other as SparseVector;
if (sparseother == null)
{
base.Add(other, result);
}
else
{
CopyTo(result);
sparse.AddScaledSparseVector(Complex.One, sparseother);
}
result.At(index, At(index) + other.At(index));
}
}
}
@ -637,130 +460,39 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
}
/// <summary>
/// Subtracts a complex from each element of the vector.
/// </summary>
/// <param name="complex">The complex to subtract.</param>
/// <returns>A new vector containing the subtraction of this vector and the complex.</returns>
public override Vector<Complex> Subtract(Complex complex)
{
if (complex == Complex.Zero)
{
return Clone();
}
var copy = (SparseVector)Clone();
for (var i = 0; i < Count; i++)
{
copy[i] -= complex;
}
return copy;
}
/// <summary>
/// Subtracts a complex from each element of the vector and stores the result in the result vector.
/// </summary>
/// <param name="complex">The complex to subtract.</param>
/// <param name="result">The vector to store the result of the subtraction.</param>
/// <exception cref="ArgumentNullException">If the result vector is <see langword="null"/>.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
public override void Subtract(Complex complex, Vector<Complex> result)
{
if (result == null)
{
throw new ArgumentNullException("result");
}
if (Count != result.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
}
base.Subtract(complex, result);
}
/// <summary>
/// Subtracts another vector from this vector.
/// Subtracts a scalar from each element of the vector and stores the result in the result vector.
/// </summary>
/// <param name="other">The vector to subtract from this one.</param>
/// <returns>A new vector containing the subtraction of the the two vectors.</returns>
/// <exception cref="ArgumentNullException">If the other vector is <see langword="null"/>.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
public override Vector<Complex> Subtract(Vector<Complex> other)
/// <param name="scalar">
/// The scalar to subtract.
/// </param>
/// <param name="result">
/// The vector to store the result of the subtraction.
/// </param>
protected override void DoSubtract(Complex scalar, Vector<Complex> result)
{
if (other == null)
{
throw new ArgumentNullException("other");
}
if (Count != other.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
}
var sparseVector = other as SparseVector;
if (sparseVector == null)
{
return base.Subtract(other);
}
var copy = (SparseVector)Clone();
copy.AddScaledSparseVector(-Complex.One, sparseVector);
return copy;
DoAdd(-scalar, result);
}
/// <summary>
/// Subtracts another vector to this vector and stores the result into the result vector.
/// </summary>
/// <param name="other">The vector to subtract from this one.</param>
/// <param name="result">The vector to store the result of the subtraction.</param>
/// <exception cref="ArgumentNullException">If the other vector is <see langword="null"/>.</exception>
/// <exception cref="ArgumentNullException">If the result vector is <see langword="null"/>.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
public override void Subtract(Vector<Complex> other, Vector<Complex> result)
/// <param name="other">
/// The vector to subtract from this one.
/// </param>
/// <param name="result">
/// The vector to store the result of the subtraction.
/// </param>
protected override void DoSubtract(Vector<Complex> other, Vector<Complex> result)
{
if (result == null)
{
throw new ArgumentNullException("result");
}
if (Count != other.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
}
if (Count != result.Count)
if (ReferenceEquals(this, other))
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
result.Clear();
return;
}
if (ReferenceEquals(this, result) || ReferenceEquals(other, result))
{
var tmp = Subtract(other);
tmp.CopyTo(result);
}
else
for (var index = 0; index < Count; index++)
{
var sparse = result as SparseVector;
if (sparse == null)
{
base.Subtract(other, result);
}
else
{
var sparseother = other as SparseVector;
if (sparseother == null)
{
base.Subtract(other, result);
}
else
{
CopyTo(result);
sparse.AddScaledSparseVector(-Complex.One, sparseother);
}
}
result.At(index, At(index) - other.At(index));
}
}
@ -835,52 +567,81 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
}
/// <summary>
/// Multiplies a complex to each element of the vector.
/// Multiplies a scalar to each element of the vector and stores the result in the result vector.
/// </summary>
/// <param name="complex">The complex to multiply.</param>
/// <returns>A new vector that is the multiplication of the vector and the complex.</returns>
public override Vector<Complex> Multiply(Complex complex)
/// <param name="scalar">
/// The scalar to multiply.
/// </param>
/// <param name="result">
/// The vector to store the result of the multiplication.
/// </param>
protected override void DoMultiply(Complex scalar, Vector<Complex> result)
{
if (complex == Complex.One)
if (scalar == Complex.One)
{
return Clone();
if (!ReferenceEquals(this, result))
{
CopyTo(result);
}
return;
}
if (complex == Complex.Zero)
if (scalar == Complex.Zero)
{
return new SparseVector(Count);
result.Clear();
return;
}
var copy = new SparseVector(this);
Control.LinearAlgebraProvider.ScaleArray(complex, copy._nonZeroValues, copy._nonZeroValues);
return copy;
var sparseResult = result as SparseVector;
if (sparseResult == null)
{
result.Clear();
for (var index = 0; index < NonZerosCount; index++)
{
result.At(_nonZeroIndices[index], scalar * _nonZeroValues[index]);
}
}
else
{
if (!ReferenceEquals(this, result))
{
sparseResult.NonZerosCount = NonZerosCount;
sparseResult._nonZeroIndices = new int[NonZerosCount];
Buffer.BlockCopy(_nonZeroIndices, 0, sparseResult._nonZeroIndices, 0, _nonZeroIndices.Length * Constants.SizeOfInt);
sparseResult._nonZeroValues = new Complex[_nonZeroValues.Length];
}
Control.LinearAlgebraProvider.ScaleArray(scalar, _nonZeroValues, sparseResult._nonZeroValues);
}
}
/// <summary>
/// Computes the dot product between this vector and another vector.
/// </summary>
/// <param name="other">The other vector to add.</param>
/// <returns>The result of the addition.</returns>
/// <exception cref="ArgumentException">If <paramref name="other"/> is not of the same size.</exception>
/// <exception cref="ArgumentNullException">If <paramref name="other"/> is <see langword="null" />.</exception>
public override Complex DotProduct(Vector<Complex> other)
/// <param name="other">
/// The other vector to add.
/// </param>
/// <returns>s
/// The result of the addition.
/// </returns>
protected override Complex DoDotProduct(Vector<Complex> other)
{
if (other == null)
{
throw new ArgumentNullException("other");
}
var result = Complex.Zero;
if (Count != other.Count)
if (ReferenceEquals(this, other))
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
for (var i = 0; i < NonZerosCount; i++)
{
result += _nonZeroValues[i] * _nonZeroValues[i];
}
}
var result = Complex.Zero;
// base implementation iterates though all elements, but we need only take non-zeros
for (var i = 0; i < NonZerosCount; i++)
else
{
result += _nonZeroValues[i] * other[_nonZeroIndices[i]];
for (var i = 0; i < NonZerosCount; i++)
{
result += _nonZeroValues[i] * other.At(_nonZeroIndices[i]);
}
}
return result;
@ -1023,7 +784,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
var result = new SparseVector(length);
for (var i = index; i < index + length; i++)
{
result[i - index] = this[i];
result.At(i - index, At(i));
}
return result;
@ -1049,7 +810,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
for (var i = 0; i < values.Length; i++)
{
this[i] = values[i];
At(i, values[i]);
}
}
@ -1084,35 +845,51 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
}
/// <summary>
/// Pointwise multiplies this vector with another vector.
/// Pointwise multiplies this vector with another vector and stores the result into the result vector.
/// </summary>
/// <param name="other">The vector to pointwise multiply with this one.</param>
/// <returns>A new vector which is the pointwise multiplication of the two vectors.</returns>
/// <exception cref="ArgumentNullException">If the other vector is <see langword="null" />.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
public override Vector<Complex> PointwiseMultiply(Vector<Complex> other)
/// <param name="result">The vector to store the result of the pointwise multiplication.</param>
protected override void DoPointwiseMultiply(Vector<Complex> other, Vector<Complex> result)
{
if (other == null)
if (ReferenceEquals(this, other))
{
throw new ArgumentNullException("other");
for (var i = 0; i < NonZerosCount; i++)
{
_nonZeroValues[i] *= _nonZeroValues[i];
}
}
if (Count != other.Count)
else
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
for (var i = 0; i < NonZerosCount; i++)
{
var index = _nonZeroIndices[i];
result.At(index, other.At(index) * _nonZeroValues[i]);
}
}
}
var copy = new SparseVector(Count);
for (var i = 0; i < _nonZeroIndices.Length; i++)
/// <summary>
/// Pointwise multiplies this vector with another vector and stores the result into the result vector.
/// </summary>
/// <param name="other">The vector to pointwise multiply with this one.</param>
/// <param name="result">The vector to store the result of the pointwise multiplication.</param>
protected override void DoPointwiseDivide(Vector<Complex> other, Vector<Complex> result)
{
if (ReferenceEquals(this, other))
{
var d = _nonZeroValues[i] * other[_nonZeroIndices[i]];
if (d != Complex.Zero)
for (var i = 0; i < NonZerosCount; i++)
{
copy[_nonZeroIndices[i]] = d;
_nonZeroValues[i] /= _nonZeroValues[i];
}
}
else
{
for (var i = 0; i < NonZerosCount; i++)
{
var index = _nonZeroIndices[i];
result.At(index, _nonZeroValues[i] / other.At(index));
}
}
return copy;
}
/// <summary>
@ -1123,7 +900,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <returns>Matrix M[i,j] = u[i]*v[j] </returns>
/// <exception cref="ArgumentNullException">If the u vector is <see langword="null" />.</exception>
/// <exception cref="ArgumentNullException">If the v vector is <see langword="null" />.</exception>
public static Matrix<Complex> /*SparseMatrix*/ OuterProduct(SparseVector u, SparseVector v)
public static Matrix<Complex> OuterProduct(SparseVector u, SparseVector v)
{
if (u == null)
{
@ -1361,13 +1138,25 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
}
#endregion
/// <summary>
/// Gets the value at the given index.
/// </summary>
/// <param name="index">Value real index in array</param>
/// <returns>The value at the given index.</returns>
internal protected override Complex At(int index)
{
// Search if item idex exists in NonZeroIndices array in range "0 - real nonzero values count"
var itemIndex = Array.BinarySearch(_nonZeroIndices, 0, NonZerosCount, index);
return itemIndex >= 0 ? _nonZeroValues[itemIndex] : Complex.Zero;
}
/// <summary>
/// Delete, Add or Update the value in NonZeroValues and NonZeroIndices
/// </summary>
/// <param name="index">Value real index in array</param>
/// <param name="value">The value to set.</param>
/// <remarks>This method assume that index is between 0 and Array Size</remarks>
private void SetValue(int index, Complex value)
internal protected override void At(int index, Complex value)
{
// Search if "index" already exists in range "0 - complex nonzero values count"
var itemIndex = Array.BinarySearch(_nonZeroIndices, 0, NonZerosCount, index);
@ -1549,5 +1338,22 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
yield return new KeyValuePair<int, Complex>(_nonZeroIndices[i], _nonZeroValues[i]);
}
}
/// <summary>
/// Returns the data contained in the vector as an array.
/// </summary>
/// <returns>
/// The vector's data as an array.
/// </returns>
public override Complex[] ToArray()
{
var ret = new Complex[Count];
for (var i = 0; i < NonZerosCount; i++)
{
ret[_nonZeroIndices[i]] = _nonZeroValues[i];
}
return ret;
}
}
}

101
src/Numerics/LinearAlgebra/Complex/Vector.cs

@ -63,10 +63,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// </param>
protected override void DoAdd(Complex scalar, Vector<Complex> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] + scalar);
for (var index = 0; index < Count; index++)
{
result.At(index, At(index) + scalar);
}
}
/// <summary>
@ -80,10 +80,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// </param>
protected override void DoAdd(Vector<Complex> other, Vector<Complex> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] + other[index]);
for (var index = 0; index < Count; index++)
{
result.At(index, At(index) + other.At(index));
}
}
/// <summary>
@ -111,10 +111,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// </param>
protected override void DoSubtract(Vector<Complex> other, Vector<Complex> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] - other[index]);
for (var index = 0; index < Count; index++)
{
result.At(index, At(index) - other.At(index));
}
}
/// <summary>
@ -128,10 +128,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// </param>
protected override void DoMultiply(Complex scalar, Vector<Complex> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] * scalar);
for (var index = 0; index < Count; index++)
{
result.At(index, At(index) * scalar);
}
}
/// <summary>
@ -145,10 +145,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// </param>
protected override void DoDivide(Complex scalar, Vector<Complex> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] / scalar);
DoMultiply(1 / scalar, result);
}
/// <summary>
@ -158,10 +155,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <param name="result">The vector to store the result of the pointwise multiplication.</param>
protected override void DoPointwiseMultiply(Vector<Complex> other, Vector<Complex> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] * other[index]);
for (var index = 0; index < Count; index++)
{
result.At(index, At(index) * other.At(index));
}
}
/// <summary>
@ -171,10 +168,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <param name="result">The vector to store the result of the pointwise division.</param>
protected override void DoPointwiseDivide(Vector<Complex> other, Vector<Complex> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] / other[index]);
for (var index = 0; index < Count; index++)
{
result.At(index, At(index) / other.At(index));
}
}
/// <summary>
@ -183,7 +180,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <param name="other">
/// The other vector to add.
/// </param>
/// <returns>s
/// <returns>
/// The result of the addition.
/// </returns>
protected override Complex DoDotProduct(Vector<Complex> other)
@ -191,7 +188,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
return CommonParallel.Aggregate(
0,
Count,
i => this[i] * other[i]);
i => At(i) * other.At(i));
}
/// <summary>
@ -200,7 +197,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <returns>The value of the absolute minimum element.</returns>
public override Complex AbsoluteMinimum()
{
return this[AbsoluteMinimumIndex()].Magnitude;
return At(AbsoluteMinimumIndex()).Magnitude;
}
/// <summary>
@ -210,10 +207,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
public override int AbsoluteMinimumIndex()
{
var index = 0;
var min = this[index].Magnitude;
var min = At(index).Magnitude;
for (var i = 1; i < Count; i++)
{
var test = this[i].Magnitude;
var test = At(i).Magnitude;
if (test < min)
{
index = i;
@ -230,7 +227,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <returns>The value of the absolute maximum element.</returns>
public override Complex AbsoluteMaximum()
{
return this[AbsoluteMaximumIndex()].Magnitude;
return At(AbsoluteMaximumIndex()).Magnitude;
}
/// <summary>
@ -240,10 +237,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
public override int AbsoluteMaximumIndex()
{
var index = 0;
var max = this[index].Magnitude;
var max = At(index).Magnitude;
for (var i = 1; i < Count; i++)
{
var test = this[i].Magnitude;
var test = At(i).Magnitude;
if (test > max)
{
index = i;
@ -263,7 +260,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
return CommonParallel.Aggregate(
0,
Count,
i => this[i]);
i => At(i));
}
/// <summary>
@ -275,7 +272,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
return CommonParallel.Aggregate(
0,
Count,
i => this[i].Magnitude);
i => At(i).Magnitude);
}
/// <summary>
@ -285,7 +282,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// The p value.
/// </param>
/// <returns>
/// <c>Scalar ret = (sum(abs(this[i])^p))^(1/p)</c>
/// <c>Scalar ret = (sum(abs(At(i))^p))^(1/p)</c>
/// </returns>
public override Complex Norm(double p)
{
@ -299,14 +296,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
return CommonParallel.Select(
0,
Count,
(index, localData) => Math.Max(localData, this[index].Magnitude),
(index, localData) => Math.Max(localData, At(index).Magnitude),
Math.Max);
}
var sum = CommonParallel.Aggregate(
0,
Count,
index => Math.Pow(this[index].Magnitude, p));
index => Math.Pow(At(index).Magnitude, p));
return Math.Pow(sum, 1.0 / p);
}
@ -317,11 +314,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <param name="target">Target vector</param>
protected override void DoConjugate(Vector<Complex> target)
{
CopyTo(target);
CommonParallel.For(
0,
Count,
index => target[index] = this[index].Conjugate());
for (var index = 0; index < Count; index++)
{
target.At(index, At(index).Conjugate());
}
}
/// <summary>
@ -336,10 +332,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
public override Vector<Complex> Negate()
{
var result = CreateVector(Count);
CommonParallel.For(
0,
Count,
index => result[index] = -this[index]);
for (var index = 0; index < Count; index++)
{
result.At(index, -At(index));
}
return result;
}
@ -364,7 +361,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
var vector = CreateVector(length);
for (var index = 0; index < length; index++)
{
vector[index] = new Complex(randomDistribution.Sample(), randomDistribution.Sample());
vector.At(index, new Complex(randomDistribution.Sample(), randomDistribution.Sample()));
}
return vector;
@ -390,7 +387,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
var vector = CreateVector(length);
for (var index = 0; index < length; index++)
{
vector[index] = new Complex(randomDistribution.Sample(), randomDistribution.Sample());
vector.At(index, new Complex(randomDistribution.Sample(), randomDistribution.Sample()));
}
return vector;

16
src/Numerics/LinearAlgebra/Complex32/DenseVector.cs

@ -1387,5 +1387,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
index => otherVector.Data[index] = Data[index].Conjugate());
}
}
/// <summary>Gets the value at the given <paramref name="index"/>.</summary>
/// <param name="index">The index of the value to get or set.</param>
/// <returns>The value of the vector at the given <paramref name="index"/>.</returns>
internal protected override Complex32 At(int index)
{
return Data[index];
}
/// <summary>Sets the <paramref name="value"/> at the given <paramref name="index"/>.</summary>
/// <param name="index">The index of the value to get or set.</param>
/// <param name="value">The value to set.</param>
internal protected override void At(int index, Complex32 value)
{
Data[index] = value;
}
}
}

2
src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs

@ -189,7 +189,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// </summary>
/// <param name="matrix">The matrix to copy.</param>
public SparseMatrix(Matrix<Complex32> matrix)
: base(matrix.RowCount, matrix.ColumnCount)
: this(matrix.RowCount, matrix.ColumnCount)
{
var sparseMatrix = matrix as SparseMatrix;

548
src/Numerics/LinearAlgebra/Complex32/SparseVector.cs

@ -38,13 +38,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <summary>
/// A vector with sparse storage.
/// </summary>
/// <remarks>The sparse vector is not thread safe.</remarks>
public class SparseVector : Vector
{
/// <summary>
/// Lock object for the indexer.
/// </summary>
private readonly object _lockObject = new object();
/// <summary>
/// Gets the vector's internal data. The array containing the actual values; only the non-zero values are stored.
/// </summary>
@ -130,7 +126,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
{
for (var i = 0; i < other.Count; i++)
{
this[i] = other[i];
this[i] = other.At(i);
}
}
else
@ -174,9 +170,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// </summary>
/// <param name="array">The array to create this vector from.</param>
/// <remarks>The vector copy the array. Any changes to the vector will NOT change the array.</remarks>
public SparseVector(Complex32[] array) : this(array.Length)
public SparseVector(IList<Complex32> array) : this(array.Count)
{
for (var i = 0; i < array.Length; i++)
for (var i = 0; i < array.Count; i++)
{
this[i] = array[i];
}
@ -193,7 +189,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
var matrix = new SparseMatrix(Count, 1);
for (var i = 0; i < NonZerosCount; i++)
{
matrix[_nonZeroIndices[i], 0] = _nonZeroValues[i];
matrix.At(_nonZeroIndices[i], 0, _nonZeroValues[i]);
}
return matrix;
@ -208,7 +204,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
var matrix = new SparseMatrix(1, Count);
for (var i = 0; i < NonZerosCount; i++)
{
matrix[0, _nonZeroIndices[i]] = _nonZeroValues[i];
matrix.At(0, _nonZeroIndices[i], _nonZeroValues[i]);
}
return matrix;
@ -229,17 +225,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
throw new IndexOutOfRangeException();
}
lock (_lockObject)
{
// Search if item idex exists in NonZeroIndices array in range "0 - complex nonzero values count"
var itemIndex = Array.BinarySearch(_nonZeroIndices, 0, NonZerosCount, index);
if (itemIndex >= 0)
{
return _nonZeroValues[itemIndex];
}
}
return Complex32.Zero;
return At(index);
}
set
@ -250,10 +236,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
throw new IndexOutOfRangeException();
}
lock (_lockObject)
{
SetValue(index, value);
}
At(index, value);
}
}
@ -330,10 +313,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
var otherVector = target as SparseVector;
if (otherVector == null)
{
CommonParallel.For(
0,
Count,
index => target[index] = this[index]);
target.Clear();
for (var index = 0; index < NonZerosCount; index++)
{
target.At(_nonZeroIndices[index], _nonZeroValues[index]);
}
}
else
{
@ -392,147 +377,42 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
}
}
}
#region Operators and supplementary functions
/// <summary>
/// Adds a complex to each element of the vector.
/// Adds a scalar to each element of the vector and stores the result in the result vector.
/// </summary>
/// <param name="complex">The complex to add.</param>
/// <returns>A copy of the vector with the complex added.</returns>
public override Vector<Complex32> Add(Complex32 complex)
{
if (complex == Complex32.Zero)
{
return Clone();
}
var copy = (SparseVector)Clone();
for (var i = 0; i < Count; i++)
{
copy[i] += complex;
}
return copy;
}
/// <summary>
/// Adds a complex to each element of the vector and stores the result in the result vector.
/// </summary>
/// <param name="complex">The complex to add.</param>
/// <param name="result">The vector to store the result of the addition.</param>
/// <exception cref="ArgumentNullException">If the result vector is <see langword="null" />.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
public override void Add(Complex32 complex, Vector<Complex32> result)
{
if (result == null)
{
throw new ArgumentNullException("result");
}
if (Count != result.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
}
base.Add(complex, result);
}
/// <summary>
/// Adds another vector to this vector.
/// </summary>
/// <param name="other">The vector to add to this one.</param>
/// <returns>A new vector containing the sum of both vectors.</returns>
/// <exception cref="ArgumentNullException">If the other vector is <see langword="null" />.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
public override Vector<Complex32> Add(Vector<Complex32> other)
{
if (other == null)
{
throw new ArgumentNullException("other");
}
if (Count != other.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
}
var sparseVector = other as SparseVector;
if (sparseVector == null)
{
return base.Add(other);
}
var copy = (SparseVector)Clone();
copy.AddScaledSparseVector(Complex32.One, sparseVector);
return copy;
}
/// <summary>
/// Adds the scaled sparse vector.
/// </summary>
/// <param name="alpha">The alpha.</param>
/// <param name="other">The other.</param>
private void AddScaledSparseVector(Complex32 alpha, SparseVector other)
/// <param name="scalar">
/// The scalar to add.
/// </param>
/// <param name="result">
/// The vector to store the result of the addition.
/// </param>
protected override void DoAdd(Complex32 scalar, Vector<Complex32> result)
{
if (other == null)
if (scalar == Complex32.Zero)
{
throw new ArgumentNullException("other");
}
if (Count != other.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
}
if (!ReferenceEquals(this, result))
{
CopyTo(result);
}
if (alpha == Complex32.Zero)
{
return;
}
// I don't use ILinearAlgebraProvider because we will get no benefit due to "lock" in this[index]
// Possible fucniton in ILinearAlgebraProvider may be AddSparseVectorToScaledSparseVector(T[] y, int[] yIndices, T alpha, T[] x, int[] xIndices);
// But it require to develop value setting algorithm and due to "lock" it will be even more greedy then implemented below
if (ReferenceEquals(this, other))
if (ReferenceEquals(this, result))
{
// Adding the same instance of sparse vector. That means if we modify "this" then "other" will be modified too.
// To avoid such problem lets change values in internal storage of "this"
if (alpha == Complex32.One)
{
for (var i = 0; i < NonZerosCount; i++)
{
_nonZeroValues[i] += _nonZeroValues[i];
}
}
else if (alpha == -Complex32.One)
{
Clear(); // Vector is subtracted from itself
return;
}
else
{
for (var i = 0; i < NonZerosCount; i++)
{
_nonZeroValues[i] += alpha * _nonZeroValues[i];
}
}
CommonParallel.For(
0,
NonZerosCount,
index => _nonZeroValues[index] += scalar);
}
else
{
// "this" and "other" are different objects, so by modifying "this" the "other" object will not be changed
if (alpha == Complex32.One)
for (var index = 0; index < Count; index++)
{
for (var i = 0; i < other.NonZerosCount; i++)
{
this[other._nonZeroIndices[i]] += other._nonZeroValues[i];
}
}
else
{
for (var i = 0; i < other.NonZerosCount; i++)
{
this[other._nonZeroIndices[i]] += alpha * other._nonZeroValues[i];
}
result.At(index, At(index) + scalar);
}
}
}
@ -540,53 +420,26 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <summary>
/// Adds another vector to this vector and stores the result into the result vector.
/// </summary>
/// <param name="other">The vector to add to this one.</param>
/// <param name="result">The vector to store the result of the addition.</param>
/// <exception cref="ArgumentNullException">If the other vector is <see langword="null" />.</exception>
/// <exception cref="ArgumentNullException">If the result vector is <see langword="null" />.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
public override void Add(Vector<Complex32> other, Vector<Complex32> result)
/// <param name="other">
/// The vector to add to this one.
/// </param>
/// <param name="result">
/// The vector to store the result of the addition.
/// </param>
protected override void DoAdd(Vector<Complex32> other, Vector<Complex32> result)
{
if (result == null)
{
throw new ArgumentNullException("result");
}
if (Count != other.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
}
if (Count != result.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
}
if (ReferenceEquals(this, result) || ReferenceEquals(other, result))
if (ReferenceEquals(this, result))
{
var tmp = Add(other);
tmp.CopyTo(result);
CommonParallel.For(
0,
NonZerosCount,
index => _nonZeroValues[index] += _nonZeroValues[index]);
}
else
{
var sparse = result as SparseVector;
if (sparse == null)
for (var index = 0; index < Count; index++)
{
base.Add(other, result);
}
else
{
var sparseother = other as SparseVector;
if (sparseother == null)
{
base.Add(other, result);
}
else
{
CopyTo(result);
sparse.AddScaledSparseVector(Complex32.One, sparseother);
}
result.At(index, At(index) + other.At(index));
}
}
}
@ -637,130 +490,39 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
}
/// <summary>
/// Subtracts a complex from each element of the vector.
/// </summary>
/// <param name="complex">The complex to subtract.</param>
/// <returns>A new vector containing the subtraction of this vector and the complex.</returns>
public override Vector<Complex32> Subtract(Complex32 complex)
{
if (complex == Complex32.Zero)
{
return Clone();
}
var copy = (SparseVector)Clone();
for (var i = 0; i < Count; i++)
{
copy[i] -= complex;
}
return copy;
}
/// <summary>
/// Subtracts a complex from each element of the vector and stores the result in the result vector.
/// </summary>
/// <param name="complex">The complex to subtract.</param>
/// <param name="result">The vector to store the result of the subtraction.</param>
/// <exception cref="ArgumentNullException">If the result vector is <see langword="null"/>.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
public override void Subtract(Complex32 complex, Vector<Complex32> result)
{
if (result == null)
{
throw new ArgumentNullException("result");
}
if (Count != result.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
}
base.Subtract(complex, result);
}
/// <summary>
/// Subtracts another vector from this vector.
/// Subtracts a scalar from each element of the vector and stores the result in the result vector.
/// </summary>
/// <param name="other">The vector to subtract from this one.</param>
/// <returns>A new vector containing the subtraction of the the two vectors.</returns>
/// <exception cref="ArgumentNullException">If the other vector is <see langword="null"/>.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
public override Vector<Complex32> Subtract(Vector<Complex32> other)
/// <param name="scalar">
/// The scalar to subtract.
/// </param>
/// <param name="result">
/// The vector to store the result of the subtraction.
/// </param>
protected override void DoSubtract(Complex32 scalar, Vector<Complex32> result)
{
if (other == null)
{
throw new ArgumentNullException("other");
}
if (Count != other.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
}
var sparseVector = other as SparseVector;
if (sparseVector == null)
{
return base.Subtract(other);
}
var copy = (SparseVector)Clone();
copy.AddScaledSparseVector(-Complex32.One, sparseVector);
return copy;
DoAdd(-scalar, result);
}
/// <summary>
/// Subtracts another vector to this vector and stores the result into the result vector.
/// </summary>
/// <param name="other">The vector to subtract from this one.</param>
/// <param name="result">The vector to store the result of the subtraction.</param>
/// <exception cref="ArgumentNullException">If the other vector is <see langword="null"/>.</exception>
/// <exception cref="ArgumentNullException">If the result vector is <see langword="null"/>.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
public override void Subtract(Vector<Complex32> other, Vector<Complex32> result)
/// <param name="other">
/// The vector to subtract from this one.
/// </param>
/// <param name="result">
/// The vector to store the result of the subtraction.
/// </param>
protected override void DoSubtract(Vector<Complex32> other, Vector<Complex32> result)
{
if (result == null)
{
throw new ArgumentNullException("result");
}
if (Count != other.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
}
if (Count != result.Count)
if (ReferenceEquals(this, other))
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
result.Clear();
return;
}
if (ReferenceEquals(this, result) || ReferenceEquals(other, result))
{
var tmp = Subtract(other);
tmp.CopyTo(result);
}
else
for (var index = 0; index < Count; index++)
{
var sparse = result as SparseVector;
if (sparse == null)
{
base.Subtract(other, result);
}
else
{
var sparseother = other as SparseVector;
if (sparseother == null)
{
base.Subtract(other, result);
}
else
{
CopyTo(result);
sparse.AddScaledSparseVector(-Complex32.One, sparseother);
}
}
result.At(index, At(index) - other.At(index));
}
}
@ -835,52 +597,81 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
}
/// <summary>
/// Multiplies a complex to each element of the vector.
/// Multiplies a scalar to each element of the vector and stores the result in the result vector.
/// </summary>
/// <param name="complex">The complex to multiply.</param>
/// <returns>A new vector that is the multiplication of the vector and the complex.</returns>
public override Vector<Complex32> Multiply(Complex32 complex)
/// <param name="scalar">
/// The scalar to multiply.
/// </param>
/// <param name="result">
/// The vector to store the result of the multiplication.
/// </param>
protected override void DoMultiply(Complex32 scalar, Vector<Complex32> result)
{
if (complex == Complex32.One)
if (scalar == Complex32.One)
{
return Clone();
if (!ReferenceEquals(this, result))
{
CopyTo(result);
}
return;
}
if (complex == Complex32.Zero)
if (scalar == Complex32.Zero)
{
return new SparseVector(Count);
result.Clear();
return;
}
var copy = new SparseVector(this);
Control.LinearAlgebraProvider.ScaleArray(complex, copy._nonZeroValues, copy._nonZeroValues);
return copy;
var sparseResult = result as SparseVector;
if (sparseResult == null)
{
result.Clear();
for (var index = 0; index < NonZerosCount; index++)
{
result.At(_nonZeroIndices[index], scalar * _nonZeroValues[index]);
}
}
else
{
if (!ReferenceEquals(this, result))
{
sparseResult.NonZerosCount = NonZerosCount;
sparseResult._nonZeroIndices = new int[NonZerosCount];
Buffer.BlockCopy(_nonZeroIndices, 0, sparseResult._nonZeroIndices, 0, _nonZeroIndices.Length * Constants.SizeOfInt);
sparseResult._nonZeroValues = new Complex32[_nonZeroValues.Length];
}
Control.LinearAlgebraProvider.ScaleArray(scalar, _nonZeroValues, sparseResult._nonZeroValues);
}
}
/// <summary>
/// Computes the dot product between this vector and another vector.
/// </summary>
/// <param name="other">The other vector to add.</param>
/// <returns>The result of the addition.</returns>
/// <exception cref="ArgumentException">If <paramref name="other"/> is not of the same size.</exception>
/// <exception cref="ArgumentNullException">If <paramref name="other"/> is <see langword="null" />.</exception>
public override Complex32 DotProduct(Vector<Complex32> other)
/// <param name="other">
/// The other vector to add.
/// </param>
/// <returns>s
/// The result of the addition.
/// </returns>
protected override Complex32 DoDotProduct(Vector<Complex32> other)
{
if (other == null)
{
throw new ArgumentNullException("other");
}
var result = Complex32.Zero;
if (Count != other.Count)
if (ReferenceEquals(this, other))
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
for (var i = 0; i < NonZerosCount; i++)
{
result += _nonZeroValues[i] * _nonZeroValues[i];
}
}
var result = Complex32.Zero;
// base implementation iterates though all elements, but we need only take non-zeros
for (var i = 0; i < NonZerosCount; i++)
else
{
result += _nonZeroValues[i] * other[_nonZeroIndices[i]];
for (var i = 0; i < NonZerosCount; i++)
{
result += _nonZeroValues[i] * other.At(_nonZeroIndices[i]);
}
}
return result;
@ -1023,7 +814,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
var result = new SparseVector(length);
for (var i = index; i < index + length; i++)
{
result[i - index] = this[i];
result.At(i - index, At(i));
}
return result;
@ -1049,7 +840,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
for (var i = 0; i < values.Length; i++)
{
this[i] = values[i];
At(i, values[i]);
}
}
@ -1084,35 +875,51 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
}
/// <summary>
/// Pointwise multiplies this vector with another vector.
/// Pointwise multiplies this vector with another vector and stores the result into the result vector.
/// </summary>
/// <param name="other">The vector to pointwise multiply with this one.</param>
/// <returns>A new vector which is the pointwise multiplication of the two vectors.</returns>
/// <exception cref="ArgumentNullException">If the other vector is <see langword="null" />.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
public override Vector<Complex32> PointwiseMultiply(Vector<Complex32> other)
/// <param name="result">The vector to store the result of the pointwise multiplication.</param>
protected override void DoPointwiseMultiply(Vector<Complex32> other, Vector<Complex32> result)
{
if (other == null)
if (ReferenceEquals(this, other))
{
throw new ArgumentNullException("other");
for (var i = 0; i < NonZerosCount; i++)
{
_nonZeroValues[i] *= _nonZeroValues[i];
}
}
if (Count != other.Count)
else
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
for (var i = 0; i < NonZerosCount; i++)
{
var index = _nonZeroIndices[i];
result.At(index, other.At(index) * _nonZeroValues[i]);
}
}
}
var copy = new SparseVector(Count);
for (var i = 0; i < _nonZeroIndices.Length; i++)
/// <summary>
/// Pointwise multiplies this vector with another vector and stores the result into the result vector.
/// </summary>
/// <param name="other">The vector to pointwise multiply with this one.</param>
/// <param name="result">The vector to store the result of the pointwise multiplication.</param>
protected override void DoPointwiseDivide(Vector<Complex32> other, Vector<Complex32> result)
{
if (ReferenceEquals(this, other))
{
var d = _nonZeroValues[i] * other[_nonZeroIndices[i]];
if (d != Complex32.Zero)
for (var i = 0; i < NonZerosCount; i++)
{
copy[_nonZeroIndices[i]] = d;
_nonZeroValues[i] /= _nonZeroValues[i];
}
}
else
{
for (var i = 0; i < NonZerosCount; i++)
{
var index = _nonZeroIndices[i];
result.At(index, _nonZeroValues[i] / other.At(index));
}
}
return copy;
}
/// <summary>
@ -1361,13 +1168,25 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
}
#endregion
/// <summary>
/// Gets the value at the given index.
/// </summary>
/// <param name="index">Value real index in array</param>
/// <returns>The value at the given index.</returns>
internal protected override Complex32 At(int index)
{
// Search if item idex exists in NonZeroIndices array in range "0 - real nonzero values count"
var itemIndex = Array.BinarySearch(_nonZeroIndices, 0, NonZerosCount, index);
return itemIndex >= 0 ? _nonZeroValues[itemIndex] : Complex32.Zero;
}
/// <summary>
/// Delete, Add or Update the value in NonZeroValues and NonZeroIndices
/// </summary>
/// <param name="index">Value real index in array</param>
/// <param name="value">The value to set.</param>
/// <remarks>This method assume that index is between 0 and Array Size</remarks>
private void SetValue(int index, Complex32 value)
internal protected override void At(int index, Complex32 value)
{
// Search if "index" already exists in range "0 - complex nonzero values count"
var itemIndex = Array.BinarySearch(_nonZeroIndices, 0, NonZerosCount, index);
@ -1549,5 +1368,22 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
yield return new KeyValuePair<int, Complex32>(_nonZeroIndices[i], _nonZeroValues[i]);
}
}
/// <summary>
/// Returns the data contained in the vector as an array.
/// </summary>
/// <returns>
/// The vector's data as an array.
/// </returns>
public override Complex32[] ToArray()
{
var ret = new Complex32[Count];
for (var i = 0; i < NonZerosCount; i++)
{
ret[_nonZeroIndices[i]] = _nonZeroValues[i];
}
return ret;
}
}
}

99
src/Numerics/LinearAlgebra/Complex32/Vector.cs

@ -63,10 +63,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// </param>
protected override void DoAdd(Complex32 scalar, Vector<Complex32> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] + scalar);
for (var index = 0; index < Count; index++)
{
result.At(index, At(index) + scalar);
}
}
/// <summary>
@ -80,10 +80,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// </param>
protected override void DoAdd(Vector<Complex32> other, Vector<Complex32> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] + other[index]);
for (var index = 0; index < Count; index++)
{
result.At(index, At(index) + other.At(index));
}
}
/// <summary>
@ -111,10 +111,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// </param>
protected override void DoSubtract(Vector<Complex32> other, Vector<Complex32> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] - other[index]);
for (var index = 0; index < Count; index++)
{
result.At(index, At(index) - other.At(index));
}
}
/// <summary>
@ -128,10 +128,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// </param>
protected override void DoMultiply(Complex32 scalar, Vector<Complex32> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] * scalar);
for (var index = 0; index < Count; index++)
{
result.At(index, At(index) * scalar);
}
}
/// <summary>
@ -145,10 +145,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// </param>
protected override void DoDivide(Complex32 scalar, Vector<Complex32> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] / scalar);
DoMultiply(1 / scalar, result);
}
/// <summary>
@ -158,10 +155,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <param name="result">The vector to store the result of the pointwise multiplication.</param>
protected override void DoPointwiseMultiply(Vector<Complex32> other, Vector<Complex32> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] * other[index]);
for (var index = 0; index < Count; index++)
{
result.At(index, At(index) * other.At(index));
}
}
/// <summary>
@ -171,10 +168,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <param name="result">The vector to store the result of the pointwise division.</param>
protected override void DoPointwiseDivide(Vector<Complex32> other, Vector<Complex32> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] / other[index]);
for (var index = 0; index < Count; index++)
{
result.At(index, At(index) / other.At(index));
}
}
/// <summary>
@ -191,7 +188,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
return CommonParallel.Aggregate(
0,
Count,
i => this[i] * other[i]);
i => At(i) * other.At(i));
}
/// <summary>
@ -200,7 +197,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <returns>The value of the absolute minimum element.</returns>
public override Complex32 AbsoluteMinimum()
{
return this[AbsoluteMinimumIndex()].Magnitude;
return At(AbsoluteMinimumIndex()).Magnitude;
}
/// <summary>
@ -210,10 +207,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
public override int AbsoluteMinimumIndex()
{
var index = 0;
var min = this[index].Magnitude;
var min = At(index).Magnitude;
for (var i = 1; i < Count; i++)
{
var test = this[i].Magnitude;
var test = At(i).Magnitude;
if (test < min)
{
index = i;
@ -230,7 +227,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <returns>The value of the absolute maximum element.</returns>
public override Complex32 AbsoluteMaximum()
{
return this[AbsoluteMaximumIndex()].Magnitude;
return At(AbsoluteMaximumIndex()).Magnitude;
}
/// <summary>
@ -240,10 +237,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
public override int AbsoluteMaximumIndex()
{
var index = 0;
var max = this[index].Magnitude;
var max = At(index).Magnitude;
for (var i = 1; i < Count; i++)
{
var test = this[i].Magnitude;
var test = At(i).Magnitude;
if (test > max)
{
index = i;
@ -263,7 +260,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
return CommonParallel.Aggregate(
0,
Count,
i => this[i]);
i => At(i));
}
/// <summary>
@ -275,7 +272,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
return CommonParallel.Aggregate(
0,
Count,
i => this[i].Magnitude);
i => At(i).Magnitude);
}
/// <summary>
@ -285,7 +282,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// The p value.
/// </param>
/// <returns>
/// <c>Scalar ret = (sum(abs(this[i])^p))^(1/p)</c>
/// <c>Scalar ret = (sum(abs(At(i))^p))^(1/p)</c>
/// </returns>
public override Complex32 Norm(double p)
{
@ -299,14 +296,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
return CommonParallel.Select(
0,
Count,
(index, localData) => Math.Max(localData, this[index].Magnitude),
(index, localData) => Math.Max(localData, At(index).Magnitude),
Common.Max);
}
var sum = CommonParallel.Aggregate(
0,
Count,
index => Math.Pow(this[index].Magnitude, p));
index => Math.Pow(At(index).Magnitude, p));
return (float)Math.Pow(sum, 1.0 / p);
}
@ -317,11 +314,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <param name="target">Target vector</param>
protected override void DoConjugate(Vector<Complex32> target)
{
CopyTo(target);
CommonParallel.For(
0,
Count,
index => target[index] = this[index].Conjugate());
for (var index = 0; index < Count; index++)
{
target.At(index, At(index).Conjugate());
}
}
/// <summary>
@ -336,10 +332,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
public override Vector<Complex32> Negate()
{
var result = CreateVector(Count);
CommonParallel.For(
0,
Count,
index => result[index] = -this[index]);
for (var index = 0; index < Count; index++)
{
result.At(index, -At(index));
}
return result;
}
@ -410,7 +407,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
var vector = CreateVector(length);
for (var index = 0; index < length; index++)
{
vector[index] = new Complex32(Convert.ToSingle(randomDistribution.Sample()), Convert.ToSingle(randomDistribution.Sample()));
vector.At(index, new Complex32(Convert.ToSingle(randomDistribution.Sample()), Convert.ToSingle(randomDistribution.Sample())));
}
return vector;
@ -436,7 +433,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
var vector = CreateVector(length);
for (var index = 0; index < length; index++)
{
vector[index] = new Complex32(Convert.ToSingle(randomDistribution.Sample()), Convert.ToSingle(randomDistribution.Sample()));
vector.At(index, new Complex32(Convert.ToSingle(randomDistribution.Sample()), Convert.ToSingle(randomDistribution.Sample())));
}
return vector;

16
src/Numerics/LinearAlgebra/Double/DenseVector.cs

@ -1393,5 +1393,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double
{
Array.Clear(Data, 0, Data.Length);
}
/// <summary>Gets the value at the given <paramref name="index"/>.</summary>
/// <param name="index">The index of the value to get or set.</param>
/// <returns>The value of the vector at the given <paramref name="index"/>.</returns>
internal protected override double At(int index)
{
return Data[index];
}
/// <summary>Sets the <paramref name="value"/> at the given <paramref name="index"/>.</summary>
/// <param name="index">The index of the value to get or set.</param>
/// <param name="value">The value to set.</param>
internal protected override void At(int index, double value)
{
Data[index] = value;
}
}
}

24
src/Numerics/LinearAlgebra/Double/SparseMatrix.cs

@ -187,7 +187,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// the values from the given matrix.
/// </summary>
/// <param name="matrix">The matrix to copy.</param>
public SparseMatrix(Matrix<double> matrix) : base(matrix.RowCount, matrix.ColumnCount)
public SparseMatrix(Matrix<double> matrix) : this(matrix.RowCount, matrix.ColumnCount)
{
var sparseMatrix = matrix as SparseMatrix;
@ -608,8 +608,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
{
lock (_lockObject)
{
var index = FindItem(row, column);
return index >= 0 ? _nonZeroValues[index] : 0.0;
return GetValueAt(row, column);
}
}
@ -634,6 +633,25 @@ namespace MathNet.Numerics.LinearAlgebra.Double
}
#region Internal methods - CRS storage implementation
/// <summary>
/// Retrieves the requested element without range checking.
/// </summary>
/// <param name="row">
/// The row of the element.
/// </param>
/// <param name="column">
/// The column of the element.
/// </param>
/// <returns>
/// The requested element.
/// </returns>
private double GetValueAt(int row, int column)
{
var index = FindItem(row, column);
return index >= 0 ? _nonZeroValues[index] : 0.0;
}
/// <summary>
/// Created this method because we cannot call "virtual At" in constructor of the class, but we need to do it
/// </summary>

566
src/Numerics/LinearAlgebra/Double/SparseVector.cs

@ -38,13 +38,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <summary>
/// A vector with sparse storage.
/// </summary>
/// <remarks>The sparse vector is not thread safe.</remarks>
public class SparseVector : Vector
{
/// <summary>
/// Lock object for the indexer.
/// </summary>
private readonly object _lockObject = new object();
/// <summary>
/// Gets the vector's internal data. The array containing the actual values; only the non-zero values are stored.
/// </summary>
@ -130,7 +126,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
{
for (var i = 0; i < other.Count; i++)
{
this[i] = other[i];
this[i] = other.At(i);
}
}
else
@ -168,9 +164,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// </summary>
/// <param name="array">The array to create this vector from.</param>
/// <remarks>The vector copy the array. Any changes to the vector will NOT change the array.</remarks>
public SparseVector(double[] array) : this(array.Length)
public SparseVector(IList<double> array) : this(array.Count)
{
for (var i = 0; i < array.Length; i++)
for (var i = 0; i < array.Count; i++)
{
this[i] = array[i];
}
@ -187,7 +183,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
var matrix = new SparseMatrix(Count, 1);
for (var i = 0; i < NonZerosCount; i++)
{
matrix[_nonZeroIndices[i], 0] = _nonZeroValues[i];
matrix.At(_nonZeroIndices[i], 0, _nonZeroValues[i]);
}
return matrix;
@ -202,55 +198,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double
var matrix = new SparseMatrix(1, Count);
for (var i = 0; i < NonZerosCount; i++)
{
matrix[0, _nonZeroIndices[i]] = _nonZeroValues[i];
matrix.At(0, _nonZeroIndices[i], _nonZeroValues[i]);
}
return matrix;
}
/// <summary>Gets or sets the value at the given <paramref name="index"/>.</summary>
/// <param name="index">The index of the value to get or set.</param>
/// <returns>The value of the vector at the given <paramref name="index"/>.</returns>
/// <exception cref="IndexOutOfRangeException">If <paramref name="index"/> is negative or
/// greater than the size of the vector.</exception>
public override double this[int index]
{
get
{
// If index is out of bounds
if ((index < 0) || (index >= Count))
{
throw new IndexOutOfRangeException();
}
lock (_lockObject)
{
// Search if item idex exists in NonZeroIndices array in range "0 - real nonzero values count"
var itemIndex = Array.BinarySearch(_nonZeroIndices, 0, NonZerosCount, index);
if (itemIndex >= 0)
{
return _nonZeroValues[itemIndex];
}
}
return 0.0;
}
set
{
// If index is out of bounds
if ((index < 0) || (index >= Count))
{
throw new IndexOutOfRangeException();
}
lock (_lockObject)
{
SetValue(index, value);
}
}
}
/// <summary>
/// Creates a matrix with the given dimensions using the same storage type
/// as this vector.
@ -324,10 +277,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double
var otherVector = target as SparseVector;
if (otherVector == null)
{
CommonParallel.For(
0,
Count,
index => target[index] = this[index]);
target.Clear();
for (var index = 0; index < NonZerosCount; index++)
{
target.At(_nonZeroIndices[index], _nonZeroValues[index]);
}
}
else
{
@ -343,145 +298,39 @@ namespace MathNet.Numerics.LinearAlgebra.Double
#region Operators and supplementary functions
/// <summary>
/// Adds a scalar to each element of the vector.
/// </summary>
/// <param name="scalar">The scalar to add.</param>
/// <returns>A copy of the vector with the scalar added.</returns>
public override Vector<double> Add(double scalar)
{
if (scalar == 0.0)
{
return Clone();
}
var copy = (SparseVector)Clone();
for (var i = 0; i < Count; i++)
{
copy[i] += scalar;
}
return copy;
}
/// <summary>
/// Adds a scalar to each element of the vector and stores the result in the result vector.
/// </summary>
/// <param name="scalar">The scalar to add.</param>
/// <param name="result">The vector to store the result of the addition.</param>
/// <exception cref="ArgumentNullException">If the result vector is <see langword="null" />.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
public override void Add(double scalar, Vector<double> result)
{
if (result == null)
{
throw new ArgumentNullException("result");
}
if (Count != result.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
}
base.Add(scalar, result);
}
/// <summary>
/// Adds another vector to this vector.
/// </summary>
/// <param name="other">The vector to add to this one.</param>
/// <returns>A new vector containing the sum of both vectors.</returns>
/// <exception cref="ArgumentNullException">If the other vector is <see langword="null" />.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
public override Vector<double> Add(Vector<double> other)
{
if (other == null)
{
throw new ArgumentNullException("other");
}
if (Count != other.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
}
var sparseVector = other as SparseVector;
if (sparseVector == null)
{
return base.Add(other);
}
var copy = (SparseVector)Clone();
copy.AddScaledSparseVector(1.0, sparseVector);
return copy;
}
/// <summary>
/// Adds the scaled sparse vector.
/// </summary>
/// <param name="alpha">The alpha.</param>
/// <param name="other">The other.</param>
private void AddScaledSparseVector(double alpha, SparseVector other)
/// <param name="scalar">
/// The scalar to add.
/// </param>
/// <param name="result">
/// The vector to store the result of the addition.
/// </param>
protected override void DoAdd(double scalar, Vector<double> result)
{
if (other == null)
{
throw new ArgumentNullException("other");
}
if (Count != other.Count)
if (scalar == 0.0)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
}
if (!ReferenceEquals(this, result))
{
CopyTo(result);
}
if (alpha == 0.0)
{
return;
}
// I don't use ILinearAlgebraProvider because we will get no benefit due to "lock" in this[index]
// Possible fucniton in ILinearAlgebraProvider may be AddSparseVectorToScaledSparseVector(T[] y, int[] yIndices, T alpha, T[] x, int[] xIndices);
// But it require to develop value setting algorithm and due to "lock" it will be even more greedy then implemented below
if (ReferenceEquals(this, other))
if (ReferenceEquals(this, result))
{
// Adding the same instance of sparse vector. That means if we modify "this" then "other" will be modified too.
// To avoid such problem lets change values in internal storage of "this"
if (alpha == 1.0)
{
for (var i = 0; i < NonZerosCount; i++)
{
_nonZeroValues[i] += _nonZeroValues[i];
}
}
else if (alpha == -1.0)
{
Clear(); // Vector is subtracted from itself
return;
}
else
{
for (var i = 0; i < NonZerosCount; i++)
{
_nonZeroValues[i] += alpha * _nonZeroValues[i];
}
}
CommonParallel.For(
0,
NonZerosCount,
index => _nonZeroValues[index] += scalar);
}
else
{
// "this" and "other" are different objects, so by modifying "this" the "other" object will not be changed
if (alpha == 1.0)
for (var index = 0; index < Count; index++)
{
for (var i = 0; i < other.NonZerosCount; i++)
{
this[other._nonZeroIndices[i]] += other._nonZeroValues[i];
}
}
else
{
for (var i = 0; i < other.NonZerosCount; i++)
{
this[other._nonZeroIndices[i]] += alpha * other._nonZeroValues[i];
}
result.At(index, At(index) + scalar);
}
}
}
@ -489,53 +338,26 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <summary>
/// Adds another vector to this vector and stores the result into the result vector.
/// </summary>
/// <param name="other">The vector to add to this one.</param>
/// <param name="result">The vector to store the result of the addition.</param>
/// <exception cref="ArgumentNullException">If the other vector is <see langword="null" />.</exception>
/// <exception cref="ArgumentNullException">If the result vector is <see langword="null" />.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
public override void Add(Vector<double> other, Vector<double> result)
/// <param name="other">
/// The vector to add to this one.
/// </param>
/// <param name="result">
/// The vector to store the result of the addition.
/// </param>
protected override void DoAdd(Vector<double> other, Vector<double> result)
{
if (result == null)
{
throw new ArgumentNullException("result");
}
if (Count != other.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
}
if (Count != result.Count)
if (ReferenceEquals(this, result))
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
}
if (ReferenceEquals(this, result) || ReferenceEquals(other, result))
{
var tmp = Add(other);
tmp.CopyTo(result);
CommonParallel.For(
0,
NonZerosCount,
index => _nonZeroValues[index] += _nonZeroValues[index]);
}
else
{
var sparse = result as SparseVector;
if (sparse == null)
{
base.Add(other, result);
}
else
for (var index = 0; index < Count; index++)
{
var sparseother = other as SparseVector;
if (sparseother == null)
{
base.Add(other, result);
}
else
{
CopyTo(result);
sparse.AddScaledSparseVector(1.0, sparseother);
}
result.At(index, At(index) + other.At(index));
}
}
}
@ -585,131 +407,40 @@ namespace MathNet.Numerics.LinearAlgebra.Double
return leftSide.Add(rightSide);
}
/// <summary>
/// Subtracts a scalar from each element of the vector.
/// </summary>
/// <param name="scalar">The scalar to subtract.</param>
/// <returns>A new vector containing the subtraction of this vector and the scalar.</returns>
public override Vector<double> Subtract(double scalar)
{
if (scalar == 0.0)
{
return Clone();
}
var copy = (SparseVector)Clone();
for (var i = 0; i < Count; i++)
{
copy[i] -= scalar;
}
return copy;
}
/// <summary>
/// Subtracts a scalar from each element of the vector and stores the result in the result vector.
/// </summary>
/// <param name="scalar">The scalar to subtract.</param>
/// <param name="result">The vector to store the result of the subtraction.</param>
/// <exception cref="ArgumentNullException">If the result vector is <see langword="null"/>.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
public override void Subtract(double scalar, Vector<double> result)
{
if (result == null)
{
throw new ArgumentNullException("result");
}
if (Count != result.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
}
base.Subtract(scalar, result);
}
/// <summary>
/// Subtracts another vector from this vector.
/// </summary>
/// <param name="other">The vector to subtract from this one.</param>
/// <returns>A new vector containing the subtraction of the the two vectors.</returns>
/// <exception cref="ArgumentNullException">If the other vector is <see langword="null"/>.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
public override Vector<double> Subtract(Vector<double> other)
/// <param name="scalar">
/// The scalar to subtract.
/// </param>
/// <param name="result">
/// The vector to store the result of the subtraction.
/// </param>
protected override void DoSubtract(double scalar, Vector<double> result)
{
if (other == null)
{
throw new ArgumentNullException("other");
}
if (Count != other.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
}
var sparseVector = other as SparseVector;
if (sparseVector == null)
{
return base.Subtract(other);
}
var copy = (SparseVector)Clone();
copy.AddScaledSparseVector(-1.0, sparseVector);
return copy;
DoAdd(-scalar, result);
}
/// <summary>
/// Subtracts another vector to this vector and stores the result into the result vector.
/// </summary>
/// <param name="other">The vector to subtract from this one.</param>
/// <param name="result">The vector to store the result of the subtraction.</param>
/// <exception cref="ArgumentNullException">If the other vector is <see langword="null"/>.</exception>
/// <exception cref="ArgumentNullException">If the result vector is <see langword="null"/>.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
public override void Subtract(Vector<double> other, Vector<double> result)
/// <param name="other">
/// The vector to subtract from this one.
/// </param>
/// <param name="result">
/// The vector to store the result of the subtraction.
/// </param>
protected override void DoSubtract(Vector<double> other, Vector<double> result)
{
if (result == null)
{
throw new ArgumentNullException("result");
}
if (Count != other.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
}
if (Count != result.Count)
if (ReferenceEquals(this, other))
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
result.Clear();
return;
}
if (ReferenceEquals(this, result) || ReferenceEquals(other, result))
for (var index = 0; index < Count; index++)
{
var tmp = Subtract(other);
tmp.CopyTo(result);
}
else
{
var sparse = result as SparseVector;
if (sparse == null)
{
base.Subtract(other, result);
}
else
{
var sparseother = other as SparseVector;
if (sparseother == null)
{
base.Subtract(other, result);
}
else
{
CopyTo(result);
sparse.AddScaledSparseVector(-1.0, sparseother);
}
}
result.At(index, At(index) - other.At(index));
}
}
@ -777,6 +508,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
0,
NonZerosCount,
index => result._nonZeroValues[index] = -_nonZeroValues[index]);
Buffer.BlockCopy(_nonZeroIndices, 0, result._nonZeroIndices, 0, NonZerosCount * Constants.SizeOfInt);
}
@ -806,30 +538,81 @@ namespace MathNet.Numerics.LinearAlgebra.Double
}
/// <summary>
/// Computes the dot product between this vector and another vector.
/// Multiplies a scalar to each element of the vector and stores the result in the result vector.
/// </summary>
/// <param name="other">The other vector to add.</param>
/// <returns>The result of the addition.</returns>
/// <exception cref="ArgumentException">If <paramref name="other"/> is not of the same size.</exception>
/// <exception cref="ArgumentNullException">If <paramref name="other"/> is <see langword="null" />.</exception>
public override double DotProduct(Vector<double> other)
/// <param name="scalar">
/// The scalar to multiply.
/// </param>
/// <param name="result">
/// The vector to store the result of the multiplication.
/// </param>
protected override void DoMultiply(double scalar, Vector<double> result)
{
if (other == null)
if (scalar == 1.0)
{
throw new ArgumentNullException("other");
if (!ReferenceEquals(this, result))
{
CopyTo(result);
}
return;
}
if (Count != other.Count)
if (scalar == 0)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
result.Clear();
return;
}
double result = 0;
var sparseResult = result as SparseVector;
if (sparseResult == null)
{
result.Clear();
for (var index = 0; index < NonZerosCount; index++)
{
result.At(_nonZeroIndices[index], scalar * _nonZeroValues[index]);
}
}
else
{
if (!ReferenceEquals(this, result))
{
sparseResult.NonZerosCount = NonZerosCount;
sparseResult._nonZeroIndices = new int[NonZerosCount];
Buffer.BlockCopy(_nonZeroIndices, 0, sparseResult._nonZeroIndices, 0, _nonZeroIndices.Length * Constants.SizeOfInt);
sparseResult._nonZeroValues = new double[_nonZeroValues.Length];
}
// base implementation iterates though all elements, but we need only take non-zeros
for (var i = 0; i < NonZerosCount; i++)
Control.LinearAlgebraProvider.ScaleArray(scalar, _nonZeroValues, sparseResult._nonZeroValues);
}
}
/// <summary>
/// Computes the dot product between this vector and another vector.
/// </summary>
/// <param name="other">
/// The other vector to add.
/// </param>
/// <returns>s
/// The result of the addition.
/// </returns>
protected override double DoDotProduct(Vector<double> other)
{
var result = 0.0;
if (ReferenceEquals(this, other))
{
for (var i = 0; i < NonZerosCount; i++)
{
result += _nonZeroValues[i] * _nonZeroValues[i];
}
}
else
{
result += _nonZeroValues[i] * other[_nonZeroIndices[i]];
for (var i = 0; i < NonZerosCount; i++)
{
result += _nonZeroValues[i] * other.At(_nonZeroIndices[i]);
}
}
return result;
@ -972,7 +755,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
var result = new SparseVector(length);
for (var i = index; i < index + length; i++)
{
result[i - index] = this[i];
result.At(i - index, At(i));
}
return result;
@ -998,7 +781,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
for (var i = 0; i < values.Length; i++)
{
this[i] = values[i];
At(i, values[i]);
}
}
@ -1083,35 +866,51 @@ namespace MathNet.Numerics.LinearAlgebra.Double
}
/// <summary>
/// Pointwise multiplies this vector with another vector.
/// Pointwise multiplies this vector with another vector and stores the result into the result vector.
/// </summary>
/// <param name="other">The vector to pointwise multiply with this one.</param>
/// <returns>A new vector which is the pointwise multiplication of the two vectors.</returns>
/// <exception cref="ArgumentNullException">If the other vector is <see langword="null" />.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
public override Vector<double> PointwiseMultiply(Vector<double> other)
/// <param name="result">The vector to store the result of the pointwise multiplication.</param>
protected override void DoPointwiseMultiply(Vector<double> other, Vector<double> result)
{
if (other == null)
if (ReferenceEquals(this, other))
{
throw new ArgumentNullException("other");
for (var i = 0; i < NonZerosCount; i++)
{
_nonZeroValues[i] *= _nonZeroValues[i];
}
}
if (Count != other.Count)
else
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
for (var i = 0; i < NonZerosCount; i++)
{
var index = _nonZeroIndices[i];
result.At(index, other.At(index) * _nonZeroValues[i]);
}
}
}
var copy = new SparseVector(Count);
for (var i = 0; i < _nonZeroIndices.Length; i++)
/// <summary>
/// Pointwise multiplies this vector with another vector and stores the result into the result vector.
/// </summary>
/// <param name="other">The vector to pointwise multiply with this one.</param>
/// <param name="result">The vector to store the result of the pointwise multiplication.</param>
protected override void DoPointwiseDivide(Vector<double> other, Vector<double> result)
{
if (ReferenceEquals(this, other))
{
var d = _nonZeroValues[i] * other[_nonZeroIndices[i]];
if (d != 0.0)
for (var i = 0; i < NonZerosCount; i++)
{
copy[_nonZeroIndices[i]] = d;
_nonZeroValues[i] /= _nonZeroValues[i];
}
}
else
{
for (var i = 0; i < NonZerosCount; i++)
{
var index = _nonZeroIndices[i];
result.At(index, _nonZeroValues[i] / other.At(index));
}
}
return copy;
}
/// <summary>
@ -1122,7 +921,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <returns>Matrix M[i,j] = u[i]*v[j] </returns>
/// <exception cref="ArgumentNullException">If the u vector is <see langword="null" />.</exception>
/// <exception cref="ArgumentNullException">If the v vector is <see langword="null" />.</exception>
public static Matrix<double> /*SparseMatrix*/ OuterProduct(SparseVector u, SparseVector v)
public static Matrix<double> OuterProduct(SparseVector u, SparseVector v)
{
if (u == null)
{
@ -1361,13 +1160,25 @@ namespace MathNet.Numerics.LinearAlgebra.Double
#endregion
/// <summary>
/// Gets the value at the given index.
/// </summary>
/// <param name="index">Value real index in array</param>
/// <returns>The value at the given index.</returns>
internal protected override double At(int index)
{
// Search if item idex exists in NonZeroIndices array in range "0 - real nonzero values count"
var itemIndex = Array.BinarySearch(_nonZeroIndices, 0, NonZerosCount, index);
return itemIndex >= 0 ? _nonZeroValues[itemIndex] : 0.0;
}
/// <summary>
/// Delete, Add or Update the value in NonZeroValues and NonZeroIndices
/// </summary>
/// <param name="index">Value real index in array</param>
/// <param name="value">The value to set.</param>
/// <remarks>This method assume that index is between 0 and Array Size</remarks>
private void SetValue(int index, double value)
internal protected override void At(int index, double value)
{
// Search if "index" already exists in range "0 - real nonzero values count"
var itemIndex = Array.BinarySearch(_nonZeroIndices, 0, NonZerosCount, index);
@ -1562,5 +1373,22 @@ namespace MathNet.Numerics.LinearAlgebra.Double
yield return new KeyValuePair<int, double>(_nonZeroIndices[i], _nonZeroValues[i]);
}
}
/// <summary>
/// Returns the data contained in the vector as an array.
/// </summary>
/// <returns>
/// The vector's data as an array.
/// </returns>
public override double[] ToArray()
{
var ret = new double[Count];
for (var i = 0; i < NonZerosCount; i++)
{
ret[_nonZeroIndices[i]] = _nonZeroValues[i];
}
return ret;
}
}
}

96
src/Numerics/LinearAlgebra/Double/Vector.cs

@ -62,10 +62,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// </param>
protected override void DoAdd(double scalar, Vector<double> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] + scalar);
for (var index = 0; index < Count; index++)
{
result.At(index, At(index) + scalar);
}
}
/// <summary>
@ -79,10 +79,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// </param>
protected override void DoAdd(Vector<double> other, Vector<double> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] + other[index]);
for (var index = 0; index < Count; index++)
{
result.At(index, At(index) + other.At(index));
}
}
/// <summary>
@ -110,10 +110,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// </param>
protected override void DoSubtract(Vector<double> other, Vector<double> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] - other[index]);
for (var index = 0; index < Count; index++)
{
result.At(index, At(index) - other.At(index));
}
}
/// <summary>
@ -127,10 +127,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// </param>
protected override void DoMultiply(double scalar, Vector<double> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] * scalar);
for (var index = 0; index < Count; index++)
{
result.At(index, At(index) * scalar);
}
}
/// <summary>
@ -144,10 +144,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// </param>
protected override void DoDivide(double scalar, Vector<double> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] / scalar);
DoMultiply(1 / scalar, result);
}
/// <summary>
@ -157,10 +154,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">The vector to store the result of the pointwise multiplication.</param>
protected override void DoPointwiseMultiply(Vector<double> other, Vector<double> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] * other[index]);
for (var index = 0; index < Count; index++)
{
result.At(index, At(index) * other.At(index));
}
}
/// <summary>
@ -170,10 +167,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">The vector to store the result of the pointwise division.</param>
protected override void DoPointwiseDivide(Vector<double> other, Vector<double> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] / other[index]);
for (var index = 0; index < Count; index++)
{
result.At(index, At(index) / other.At(index));
}
}
/// <summary>
@ -190,7 +187,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
return CommonParallel.Aggregate(
0,
Count,
i => this[i] * other[i]);
i => At(i) * other.At(i));
}
/// <summary>
@ -199,7 +196,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <returns>The value of the absolute minimum element.</returns>
public override double AbsoluteMinimum()
{
return Math.Abs(this[AbsoluteMinimumIndex()]);
return Math.Abs(At(AbsoluteMinimumIndex()));
}
/// <summary>
@ -209,10 +206,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
public override int AbsoluteMinimumIndex()
{
var index = 0;
var min = Math.Abs(this[index]);
var min = Math.Abs(At(index));
for (var i = 1; i < Count; i++)
{
var test = Math.Abs(this[i]);
var test = Math.Abs(At(i));
if (test < min)
{
index = i;
@ -229,7 +226,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <returns>The value of the absolute maximum element.</returns>
public override double AbsoluteMaximum()
{
return Math.Abs(this[AbsoluteMaximumIndex()]);
return Math.Abs(At(AbsoluteMaximumIndex()));
}
/// <summary>
@ -239,10 +236,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
public override int AbsoluteMaximumIndex()
{
var index = 0;
var max = Math.Abs(this[index]);
var max = Math.Abs(At(index));
for (var i = 1; i < Count; i++)
{
var test = Math.Abs(this[i]);
var test = Math.Abs(At(i));
if (test > max)
{
index = i;
@ -262,7 +259,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
return CommonParallel.Aggregate(
0,
Count,
i => this[i]);
i => At(i));
}
/// <summary>
@ -274,7 +271,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
return CommonParallel.Aggregate(
0,
Count,
i => Math.Abs(this[i]));
i => Math.Abs(At(i)));
}
/// <summary>
@ -298,14 +295,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double
return CommonParallel.Select(
0,
Count,
(index, localData) => Math.Max(localData, Math.Abs(this[index])),
(index, localData) => Math.Max(localData, Math.Abs(At(index))),
Math.Max);
}
var sum = CommonParallel.Aggregate(
0,
Count,
index => Math.Pow(Math.Abs(this[index]), p));
index => Math.Pow(Math.Abs(At(index)), p));
return Math.Pow(sum, 1.0 / p);
}
@ -336,10 +333,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double
public override Vector<double> Negate()
{
var result = CreateVector(Count);
CommonParallel.For(
0,
Count,
index => result[index] = -this[index]);
for (var index = 0; index < Count; index++)
{
result.At(index, -At(index));
}
return result;
}
@ -351,10 +349,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
public override int MaximumIndex()
{
var index = 0;
var max = this[index];
var max = At(index);
for (var i = 1; i < Count; i++)
{
var test = this[i];
var test = At(i);
if (test > max)
{
index = i;
@ -372,10 +370,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
public override int MinimumIndex()
{
var index = 0;
var min = this[index];
var min = At(index);
for (var i = 1; i < Count; i++)
{
var test = this[i];
var test = At(i);
if (test < min)
{
index = i;
@ -434,7 +432,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
var v = CreateVector(length);
for (var index = 0; index < Count; index++)
{
v[index] = randomDistribution.Sample();
v.At(index, randomDistribution.Sample());
}
return v;
@ -460,7 +458,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
var v = CreateVector(length);
for (var index = 0; index < Count; index++)
{
this[index] = randomDistribution.Sample();
At(index, randomDistribution.Sample());
}
return v;

36
src/Numerics/LinearAlgebra/Generic/Vector.cs

@ -93,10 +93,27 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
/// <returns>The value of the vector at the given <paramref name="index"/>.</returns>
/// <exception cref="IndexOutOfRangeException">If <paramref name="index"/> is negative or
/// greater than the size of the vector.</exception>
public abstract T this[int index]
public virtual T this[int index]
{
get;
set;
get
{
if (index < 0 || index >= Count)
{
throw new ArgumentOutOfRangeException("index");
}
return At(index);
}
set
{
if (index < 0 || index >= Count)
{
throw new ArgumentOutOfRangeException("index");
}
At(index, value);
}
}
/// <summary>
@ -715,7 +732,6 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
/// Pointwise divide this vector with another vector and stores the result into the result vector.
/// </summary>
/// <param name="other">The vector to pointwise divide this one by.</param>
/// <param name="result">The vector to store the result of the pointwise division.</param>
protected abstract void DoPointwiseDivide(Vector<T> other, Vector<T> result);
/// <summary>
@ -1207,7 +1223,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
var ret = new T[Count];
for (var i = 0; i < ret.Length; i++)
{
ret[i] = this[i];
ret[i] = At(i);
}
return ret;
@ -1534,5 +1550,15 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
{
CommonParallel.For(0, Count, index => this[index] = default(T));
}
/// <summary>Gets the value at the given <paramref name="index"/>.</summary>
/// <param name="index">The index of the value to get or set.</param>
/// <returns>The value of the vector at the given <paramref name="index"/>.</returns>
internal protected abstract T At(int index);
/// <summary>Sets the <paramref name="value"/> at the given <paramref name="index"/>.</summary>
/// <param name="index">The index of the value to get or set.</param>
/// <param name="value">The value to set.</param>
internal protected abstract void At(int index, T value);
}
}

16
src/Numerics/LinearAlgebra/Single/DenseVector.cs

@ -1393,5 +1393,21 @@ namespace MathNet.Numerics.LinearAlgebra.Single
{
Array.Clear(Data, 0, Data.Length);
}
/// <summary>Gets the value at the given <paramref name="index"/>.</summary>
/// <param name="index">The index of the value to get or set.</param>
/// <returns>The value of the vector at the given <paramref name="index"/>.</returns>
internal protected override float At(int index)
{
return Data[index];
}
/// <summary>Sets the <paramref name="value"/> at the given <paramref name="index"/>.</summary>
/// <param name="index">The index of the value to get or set.</param>
/// <param name="value">The value to set.</param>
internal protected override void At(int index, float value)
{
Data[index] = value;
}
}
}

2
src/Numerics/LinearAlgebra/Single/SparseMatrix.cs

@ -188,7 +188,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// </summary>
/// <param name="matrix">The matrix to copy.</param>
public SparseMatrix(Matrix<float> matrix)
: base(matrix.RowCount, matrix.ColumnCount)
: this(matrix.RowCount, matrix.ColumnCount)
{
var sparseMatrix = matrix as SparseMatrix;

539
src/Numerics/LinearAlgebra/Single/SparseVector.cs

@ -38,13 +38,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <summary>
/// A vector with sparse storage.
/// </summary>
/// <remarks>The sparse vector is not thread safe.</remarks>
public class SparseVector : Vector
{
/// <summary>
/// Lock object for the indexer.
/// </summary>
private readonly object _lockObject = new object();
/// <summary>
/// Gets the vector's internal data. The array containing the actual values; only the non-zero values are stored.
/// </summary>
@ -130,7 +126,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
{
for (var i = 0; i < other.Count; i++)
{
this[i] = other[i];
this[i] = other.At(i);
}
}
else
@ -168,9 +164,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// </summary>
/// <param name="array">The array to create this vector from.</param>
/// <remarks>The vector copy the array. Any changes to the vector will NOT change the array.</remarks>
public SparseVector(float[] array) : this(array.Length)
public SparseVector(IList<float> array) : this(array.Count)
{
for (var i = 0; i < array.Length; i++)
for (var i = 0; i < array.Count; i++)
{
this[i] = array[i];
}
@ -187,7 +183,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
var matrix = new SparseMatrix(Count, 1);
for (var i = 0; i < NonZerosCount; i++)
{
matrix[_nonZeroIndices[i], 0] = _nonZeroValues[i];
matrix.At(_nonZeroIndices[i], 0, _nonZeroValues[i]);
}
return matrix;
@ -202,7 +198,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
var matrix = new SparseMatrix(1, Count);
for (var i = 0; i < NonZerosCount; i++)
{
matrix[0, _nonZeroIndices[i]] = _nonZeroValues[i];
matrix.At(0, _nonZeroIndices[i], _nonZeroValues[i]);
}
return matrix;
@ -223,17 +219,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
throw new IndexOutOfRangeException();
}
lock (_lockObject)
{
// Search if item idex exists in NonZeroIndices array in range "0 - real nonzero values count"
var itemIndex = Array.BinarySearch(_nonZeroIndices, 0, NonZerosCount, index);
if (itemIndex >= 0)
{
return _nonZeroValues[itemIndex];
}
}
return 0.0f;
return At(index);
}
set
@ -244,10 +230,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
throw new IndexOutOfRangeException();
}
lock (_lockObject)
{
SetValue(index, value);
}
At(index, value);
}
}
@ -324,10 +307,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single
var otherVector = target as SparseVector;
if (otherVector == null)
{
CommonParallel.For(
0,
Count,
index => target[index] = this[index]);
target.Clear();
for (var index = 0; index < NonZerosCount; index++)
{
target.At(_nonZeroIndices[index], _nonZeroValues[index]);
}
}
else
{
@ -343,145 +328,39 @@ namespace MathNet.Numerics.LinearAlgebra.Single
#region Operators and supplementary functions
/// <summary>
/// Adds a scalar to each element of the vector.
/// </summary>
/// <param name="scalar">The scalar to add.</param>
/// <returns>A copy of the vector with the scalar added.</returns>
public override Vector<float> Add(float scalar)
{
if (scalar == 0.0)
{
return Clone();
}
var copy = (SparseVector)Clone();
for (var i = 0; i < Count; i++)
{
copy[i] += scalar;
}
return copy;
}
/// <summary>
/// Adds a scalar to each element of the vector and stores the result in the result vector.
/// </summary>
/// <param name="scalar">The scalar to add.</param>
/// <param name="result">The vector to store the result of the addition.</param>
/// <exception cref="ArgumentNullException">If the result vector is <see langword="null" />.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
public override void Add(float scalar, Vector<float> result)
{
if (result == null)
{
throw new ArgumentNullException("result");
}
if (Count != result.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
}
base.Add(scalar, result);
}
/// <summary>
/// Adds another vector to this vector.
/// </summary>
/// <param name="other">The vector to add to this one.</param>
/// <returns>A new vector containing the sum of both vectors.</returns>
/// <exception cref="ArgumentNullException">If the other vector is <see langword="null" />.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
public override Vector<float> Add(Vector<float> other)
{
if (other == null)
{
throw new ArgumentNullException("other");
}
if (Count != other.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
}
var sparseVector = other as SparseVector;
if (sparseVector == null)
{
return base.Add(other);
}
var copy = (SparseVector)Clone();
copy.AddScaledSparseVector(1.0f, sparseVector);
return copy;
}
/// <summary>
/// Adds the scaled sparse vector.
/// </summary>
/// <param name="alpha">The alpha.</param>
/// <param name="other">The other.</param>
private void AddScaledSparseVector(float alpha, SparseVector other)
/// <param name="scalar">
/// The scalar to add.
/// </param>
/// <param name="result">
/// The vector to store the result of the addition.
/// </param>
protected override void DoAdd(float scalar, Vector<float> result)
{
if (other == null)
{
throw new ArgumentNullException("other");
}
if (Count != other.Count)
if (scalar == 0.0f)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
}
if (!ReferenceEquals(this, result))
{
CopyTo(result);
}
if (alpha == 0.0)
{
return;
}
// I don't use ILinearAlgebraProvider because we will get no benefit due to "lock" in this[index]
// Possible fucniton in ILinearAlgebraProvider may be AddSparseVectorToScaledSparseVector(T[] y, int[] yIndices, T alpha, T[] x, int[] xIndices);
// But it require to develop value setting algorithm and due to "lock" it will be even more greedy then implemented below
if (ReferenceEquals(this, other))
if (ReferenceEquals(this, result))
{
// Adding the same instance of sparse vector. That means if we modify "this" then "other" will be modified too.
// To avoid such problem lets change values in internal storage of "this"
if (alpha == 1.0)
{
for (var i = 0; i < NonZerosCount; i++)
{
_nonZeroValues[i] += _nonZeroValues[i];
}
}
else if (alpha == -1.0)
{
Clear(); // Vector is subtracted from itself
return;
}
else
{
for (var i = 0; i < NonZerosCount; i++)
{
_nonZeroValues[i] += alpha * _nonZeroValues[i];
}
}
CommonParallel.For(
0,
NonZerosCount,
index => _nonZeroValues[index] += scalar);
}
else
{
// "this" and "other" are different objects, so by modifying "this" the "other" object will not be changed
if (alpha == 1.0)
{
for (var i = 0; i < other.NonZerosCount; i++)
{
this[other._nonZeroIndices[i]] += other._nonZeroValues[i];
}
}
else
for (var index = 0; index < Count; index++)
{
for (var i = 0; i < other.NonZerosCount; i++)
{
this[other._nonZeroIndices[i]] += alpha * other._nonZeroValues[i];
}
result.At(index, At(index) + scalar);
}
}
}
@ -489,53 +368,26 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <summary>
/// Adds another vector to this vector and stores the result into the result vector.
/// </summary>
/// <param name="other">The vector to add to this one.</param>
/// <param name="result">The vector to store the result of the addition.</param>
/// <exception cref="ArgumentNullException">If the other vector is <see langword="null" />.</exception>
/// <exception cref="ArgumentNullException">If the result vector is <see langword="null" />.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
public override void Add(Vector<float> other, Vector<float> result)
/// <param name="other">
/// The vector to add to this one.
/// </param>
/// <param name="result">
/// The vector to store the result of the addition.
/// </param>
protected override void DoAdd(Vector<float> other, Vector<float> result)
{
if (result == null)
{
throw new ArgumentNullException("result");
}
if (Count != other.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
}
if (Count != result.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
}
if (ReferenceEquals(this, result) || ReferenceEquals(other, result))
if (ReferenceEquals(this, result))
{
var tmp = Add(other);
tmp.CopyTo(result);
CommonParallel.For(
0,
NonZerosCount,
index => _nonZeroValues[index] += _nonZeroValues[index]);
}
else
{
var sparse = result as SparseVector;
if (sparse == null)
for (var index = 0; index < Count; index++)
{
base.Add(other, result);
}
else
{
var sparseother = other as SparseVector;
if (sparseother == null)
{
base.Add(other, result);
}
else
{
CopyTo(result);
sparse.AddScaledSparseVector(1.0f, sparseother);
}
result.At(index, At(index) + other.At(index));
}
}
}
@ -585,131 +437,40 @@ namespace MathNet.Numerics.LinearAlgebra.Single
return leftSide.Add(rightSide);
}
/// <summary>
/// Subtracts a scalar from each element of the vector.
/// </summary>
/// <param name="scalar">The scalar to subtract.</param>
/// <returns>A new vector containing the subtraction of this vector and the scalar.</returns>
public override Vector<float> Subtract(float scalar)
{
if (scalar == 0.0)
{
return Clone();
}
var copy = (SparseVector)Clone();
for (var i = 0; i < Count; i++)
{
copy[i] -= scalar;
}
return copy;
}
/// <summary>
/// Subtracts a scalar from each element of the vector and stores the result in the result vector.
/// </summary>
/// <param name="scalar">The scalar to subtract.</param>
/// <param name="result">The vector to store the result of the subtraction.</param>
/// <exception cref="ArgumentNullException">If the result vector is <see langword="null"/>.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
public override void Subtract(float scalar, Vector<float> result)
{
if (result == null)
{
throw new ArgumentNullException("result");
}
if (Count != result.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
}
base.Subtract(scalar, result);
}
/// <summary>
/// Subtracts another vector from this vector.
/// </summary>
/// <param name="other">The vector to subtract from this one.</param>
/// <returns>A new vector containing the subtraction of the the two vectors.</returns>
/// <exception cref="ArgumentNullException">If the other vector is <see langword="null"/>.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
public override Vector<float> Subtract(Vector<float> other)
/// <param name="scalar">
/// The scalar to subtract.
/// </param>
/// <param name="result">
/// The vector to store the result of the subtraction.
/// </param>
protected override void DoSubtract(float scalar, Vector<float> result)
{
if (other == null)
{
throw new ArgumentNullException("other");
}
if (Count != other.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
}
var sparseVector = other as SparseVector;
if (sparseVector == null)
{
return base.Subtract(other);
}
var copy = (SparseVector)Clone();
copy.AddScaledSparseVector(-1.0f, sparseVector);
return copy;
DoAdd(-scalar, result);
}
/// <summary>
/// Subtracts another vector to this vector and stores the result into the result vector.
/// </summary>
/// <param name="other">The vector to subtract from this one.</param>
/// <param name="result">The vector to store the result of the subtraction.</param>
/// <exception cref="ArgumentNullException">If the other vector is <see langword="null"/>.</exception>
/// <exception cref="ArgumentNullException">If the result vector is <see langword="null"/>.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="result"/> are not the same size.</exception>
public override void Subtract(Vector<float> other, Vector<float> result)
/// <param name="other">
/// The vector to subtract from this one.
/// </param>
/// <param name="result">
/// The vector to store the result of the subtraction.
/// </param>
protected override void DoSubtract(Vector<float> other, Vector<float> result)
{
if (result == null)
{
throw new ArgumentNullException("result");
}
if (Count != other.Count)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
}
if (Count != result.Count)
if (ReferenceEquals(this, other))
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "result");
result.Clear();
return;
}
if (ReferenceEquals(this, result) || ReferenceEquals(other, result))
{
var tmp = Subtract(other);
tmp.CopyTo(result);
}
else
for (var index = 0; index < Count; index++)
{
var sparse = result as SparseVector;
if (sparse == null)
{
base.Subtract(other, result);
}
else
{
var sparseother = other as SparseVector;
if (sparseother == null)
{
base.Subtract(other, result);
}
else
{
CopyTo(result);
sparse.AddScaledSparseVector(-1.0f, sparseother);
}
}
result.At(index, At(index) - other.At(index));
}
}
@ -784,52 +545,81 @@ namespace MathNet.Numerics.LinearAlgebra.Single
}
/// <summary>
/// Multiplies a scalar to each element of the vector.
/// Multiplies a scalar to each element of the vector and stores the result in the result vector.
/// </summary>
/// <param name="scalar">The scalar to multiply.</param>
/// <returns>A new vector that is the multiplication of the vector and the scalar.</returns>
public override Vector<float> Multiply(float scalar)
/// <param name="scalar">
/// The scalar to multiply.
/// </param>
/// <param name="result">
/// The vector to store the result of the multiplication.
/// </param>
protected override void DoMultiply(float scalar, Vector<float> result)
{
if (scalar == 1.0)
{
return Clone();
if (!ReferenceEquals(this, result))
{
CopyTo(result);
}
return;
}
if (scalar == 0)
{
return new SparseVector(Count);
result.Clear();
return;
}
var copy = new SparseVector(this);
Control.LinearAlgebraProvider.ScaleArray(scalar, copy._nonZeroValues, copy._nonZeroValues);
return copy;
var sparseResult = result as SparseVector;
if (sparseResult == null)
{
result.Clear();
for (var index = 0; index < NonZerosCount; index++)
{
result.At(_nonZeroIndices[index], scalar * _nonZeroValues[index]);
}
}
else
{
if (!ReferenceEquals(this, result))
{
sparseResult.NonZerosCount = NonZerosCount;
sparseResult._nonZeroIndices = new int[NonZerosCount];
Buffer.BlockCopy(_nonZeroIndices, 0, sparseResult._nonZeroIndices, 0, _nonZeroIndices.Length * Constants.SizeOfInt);
sparseResult._nonZeroValues = new float[_nonZeroValues.Length];
}
Control.LinearAlgebraProvider.ScaleArray(scalar, _nonZeroValues, sparseResult._nonZeroValues);
}
}
/// <summary>
/// Computes the dot product between this vector and another vector.
/// </summary>
/// <param name="other">The other vector to add.</param>
/// <returns>The result of the addition.</returns>
/// <exception cref="ArgumentException">If <paramref name="other"/> is not of the same size.</exception>
/// <exception cref="ArgumentNullException">If <paramref name="other"/> is <see langword="null" />.</exception>
public override float DotProduct(Vector<float> other)
/// <param name="other">
/// The other vector to add.
/// </param>
/// <returns>s
/// The result of the addition.
/// </returns>
protected override float DoDotProduct(Vector<float> other)
{
if (other == null)
{
throw new ArgumentNullException("other");
}
var result = 0.0f;
if (Count != other.Count)
if (ReferenceEquals(this, other))
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
for (var i = 0; i < NonZerosCount; i++)
{
result += _nonZeroValues[i] * _nonZeroValues[i];
}
}
float result = 0;
// base implementation iterates though all elements, but we need only take non-zeros
for (var i = 0; i < NonZerosCount; i++)
else
{
result += _nonZeroValues[i] * other[_nonZeroIndices[i]];
for (var i = 0; i < NonZerosCount; i++)
{
result += _nonZeroValues[i] * other.At(_nonZeroIndices[i]);
}
}
return result;
@ -998,7 +788,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
for (var i = 0; i < values.Length; i++)
{
this[i] = values[i];
At(i, values[i]);
}
}
@ -1083,35 +873,51 @@ namespace MathNet.Numerics.LinearAlgebra.Single
}
/// <summary>
/// Pointwise multiplies this vector with another vector.
/// Pointwise multiplies this vector with another vector and stores the result into the result vector.
/// </summary>
/// <param name="other">The vector to pointwise multiply with this one.</param>
/// <returns>A new vector which is the pointwise multiplication of the two vectors.</returns>
/// <exception cref="ArgumentNullException">If the other vector is <see langword="null" />.</exception>
/// <exception cref="ArgumentException">If this vector and <paramref name="other"/> are not the same size.</exception>
public override Vector<float> PointwiseMultiply(Vector<float> other)
/// <param name="result">The vector to store the result of the pointwise multiplication.</param>
protected override void DoPointwiseMultiply(Vector<float> other, Vector<float> result)
{
if (other == null)
if (ReferenceEquals(this, other))
{
throw new ArgumentNullException("other");
for (var i = 0; i < NonZerosCount; i++)
{
_nonZeroValues[i] *= _nonZeroValues[i];
}
}
if (Count != other.Count)
else
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");
for (var i = 0; i < NonZerosCount; i++)
{
var index = _nonZeroIndices[i];
result.At(index, other.At(index) * _nonZeroValues[i]);
}
}
}
var copy = new SparseVector(Count);
for (var i = 0; i < _nonZeroIndices.Length; i++)
/// <summary>
/// Pointwise multiplies this vector with another vector and stores the result into the result vector.
/// </summary>
/// <param name="other">The vector to pointwise multiply with this one.</param>
/// <param name="result">The vector to store the result of the pointwise multiplication.</param>
protected override void DoPointwiseDivide(Vector<float> other, Vector<float> result)
{
if (ReferenceEquals(this, other))
{
var d = _nonZeroValues[i] * other[_nonZeroIndices[i]];
if (d != 0.0)
for (var i = 0; i < NonZerosCount; i++)
{
copy[_nonZeroIndices[i]] = d;
_nonZeroValues[i] /= _nonZeroValues[i];
}
}
else
{
for (var i = 0; i < NonZerosCount; i++)
{
var index = _nonZeroIndices[i];
result.At(index, _nonZeroValues[i] / other.At(index));
}
}
return copy;
}
/// <summary>
@ -1122,7 +928,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <returns>Matrix M[i,j] = u[i]*v[j] </returns>
/// <exception cref="ArgumentNullException">If the u vector is <see langword="null" />.</exception>
/// <exception cref="ArgumentNullException">If the v vector is <see langword="null" />.</exception>
public static Matrix<float> /*SparseMatrix*/ OuterProduct(SparseVector u, SparseVector v)
public static Matrix<float> OuterProduct(SparseVector u, SparseVector v)
{
if (u == null)
{
@ -1361,13 +1167,25 @@ namespace MathNet.Numerics.LinearAlgebra.Single
#endregion
/// <summary>
/// Gets the value at the given index.
/// </summary>
/// <param name="index">Value real index in array</param>
/// <returns>The value at the given index.</returns>
internal protected override float At(int index)
{
// Search if item idex exists in NonZeroIndices array in range "0 - real nonzero values count"
var itemIndex = Array.BinarySearch(_nonZeroIndices, 0, NonZerosCount, index);
return itemIndex >= 0 ? _nonZeroValues[itemIndex] : 0.0f;
}
/// <summary>
/// Delete, Add or Update the value in NonZeroValues and NonZeroIndices
/// </summary>
/// <param name="index">Value real index in array</param>
/// <param name="value">The value to set.</param>
/// <remarks>This method assume that index is between 0 and Array Size</remarks>
private void SetValue(int index, float value)
internal protected override void At(int index, float value)
{
// Search if "index" already exists in range "0 - real nonzero values count"
var itemIndex = Array.BinarySearch(_nonZeroIndices, 0, NonZerosCount, index);
@ -1562,5 +1380,22 @@ namespace MathNet.Numerics.LinearAlgebra.Single
yield return new KeyValuePair<int, float>(_nonZeroIndices[i], _nonZeroValues[i]);
}
}
/// <summary>
/// Returns the data contained in the vector as an array.
/// </summary>
/// <returns>
/// The vector's data as an array.
/// </returns>
public override float[] ToArray()
{
var ret = new float[Count];
for (var i = 0; i < NonZerosCount; i++)
{
ret[_nonZeroIndices[i]] = _nonZeroValues[i];
}
return ret;
}
}
}

98
src/Numerics/LinearAlgebra/Single/Vector.cs

@ -63,10 +63,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// </param>
protected override void DoAdd(float scalar, Vector<float> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] + scalar);
for (var index = 0; index < Count; index++)
{
result.At(index, At(index) + scalar);
}
}
/// <summary>
@ -80,10 +80,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// </param>
protected override void DoAdd(Vector<float> other, Vector<float> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] + other[index]);
for (var index = 0; index < Count; index++)
{
result.At(index, At(index) + other.At(index));
}
}
/// <summary>
@ -111,10 +111,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// </param>
protected override void DoSubtract(Vector<float> other, Vector<float> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] - other[index]);
for (var index = 0; index < Count; index++)
{
result.At(index, At(index) - other.At(index));
}
}
/// <summary>
@ -128,10 +128,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// </param>
protected override void DoMultiply(float scalar, Vector<float> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] * scalar);
for (var index = 0; index < Count; index++)
{
result.At(index, At(index) * scalar);
}
}
/// <summary>
@ -145,10 +145,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// </param>
protected override void DoDivide(float scalar, Vector<float> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] / scalar);
DoMultiply(1 / scalar, result);
}
/// <summary>
@ -158,10 +155,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">The vector to store the result of the pointwise multiplication.</param>
protected override void DoPointwiseMultiply(Vector<float> other, Vector<float> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] * other[index]);
for (var index = 0; index < Count; index++)
{
result.At(index, At(index) * other.At(index));
}
}
/// <summary>
@ -171,10 +168,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">The vector to store the result of the pointwise division.</param>
protected override void DoPointwiseDivide(Vector<float> other, Vector<float> result)
{
CommonParallel.For(
0,
Count,
index => result[index] = this[index] / other[index]);
for (var index = 0; index < Count; index++)
{
result.At(index, At(index) / other.At(index));
}
}
/// <summary>
@ -191,7 +188,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
return CommonParallel.Aggregate(
0,
Count,
i => this[i] * other[i]);
i => At(i) * other.At(i));
}
/// <summary>
@ -200,7 +197,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <returns>The value of the absolute minimum element.</returns>
public override float AbsoluteMinimum()
{
return Math.Abs(this[AbsoluteMinimumIndex()]);
return Math.Abs(At(AbsoluteMinimumIndex()));
}
/// <summary>
@ -210,10 +207,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single
public override int AbsoluteMinimumIndex()
{
var index = 0;
var min = Math.Abs(this[index]);
var min = Math.Abs(At(index));
for (var i = 1; i < Count; i++)
{
var test = Math.Abs(this[i]);
var test = Math.Abs(At(i));
if (test < min)
{
index = i;
@ -230,7 +227,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <returns>The value of the absolute maximum element.</returns>
public override float AbsoluteMaximum()
{
return Math.Abs(this[AbsoluteMaximumIndex()]);
return Math.Abs(At(AbsoluteMaximumIndex()));
}
/// <summary>
@ -240,10 +237,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single
public override int AbsoluteMaximumIndex()
{
var index = 0;
var max = Math.Abs(this[index]);
var max = Math.Abs(At(index));
for (var i = 1; i < Count; i++)
{
var test = Math.Abs(this[i]);
var test = Math.Abs(At(i));
if (test > max)
{
index = i;
@ -263,7 +260,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
return CommonParallel.Aggregate(
0,
Count,
i => this[i]);
i => At(i));
}
/// <summary>
@ -275,7 +272,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
return CommonParallel.Aggregate(
0,
Count,
i => Math.Abs(this[i]));
i => Math.Abs(At(i)));
}
/// <summary>
@ -285,7 +282,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// The p value.
/// </param>
/// <returns>
/// <c>Scalar ret = (sum(abs(this[i])^p))^(1/p)</c>
/// <c>Scalar ret = (sum(abs(At(i))^p))^(1/p)</c>
/// </returns>
public override float Norm(double p)
{
@ -299,14 +296,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single
return CommonParallel.Select(
0,
Count,
(index, localData) => Math.Max(localData, Math.Abs(this[index])),
(index, localData) => Math.Max(localData, Math.Abs(At(index))),
Common.Max);
}
var sum = CommonParallel.Aggregate(
0,
Count,
index => Math.Pow(Math.Abs(this[index]), p));
index => Math.Pow(Math.Abs(At(index)), p));
return (float)Math.Pow(sum, 1.0 / p);
}
@ -337,10 +334,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single
public override Vector<float> Negate()
{
var result = CreateVector(Count);
CommonParallel.For(
0,
Count,
index => result[index] = -this[index]);
for (var index = 0; index < Count; index++)
{
result.At(index, -At(index));
}
return result;
}
@ -352,10 +350,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single
public override int MaximumIndex()
{
var index = 0;
var max = this[index];
var max = At(index);
for (var i = 1; i < Count; i++)
{
var test = this[i];
var test = At(i);
if (test > max)
{
index = i;
@ -373,10 +371,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single
public override int MinimumIndex()
{
var index = 0;
var min = this[index];
var min = At(index);
for (var i = 1; i < Count; i++)
{
var test = this[i];
var test = At(i);
if (test < min)
{
index = i;
@ -435,7 +433,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
var v = CreateVector(length);
for (var index = 0; index < Count; index++)
{
v[index] = Convert.ToSingle(randomDistribution.Sample());
v.At(index, Convert.ToSingle(randomDistribution.Sample()));
}
return v;
@ -461,7 +459,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
var v = CreateVector(length);
for (var index = 0; index < Count; index++)
{
v[index] = Convert.ToSingle(randomDistribution.Sample());
v.At(index, Convert.ToSingle(randomDistribution.Sample()));
}
return v;

22
src/UnitTests/LinearAlgebraTests/Complex/UserDefinedVectorTests.cs

@ -54,19 +54,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex
_data = (Complex[])data.Clone();
}
public override Complex this[int index]
{
get
{
return _data[index];
}
set
{
_data[index] = value;
}
}
public override Matrix<Complex> CreateMatrix(int rows, int columns)
{
return new UserDefinedMatrix(rows, columns);
@ -77,6 +64,15 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex
return new UserDefinedVector(size);
}
protected internal override Complex At(int index)
{
return _data[index];
}
protected internal override void At(int index, Complex value)
{
_data[index] = value;
}
}
public class UserDefinedVectorTests : VectorTests

23
src/UnitTests/LinearAlgebraTests/Complex32/UserDefinedVectorTests.cs

@ -47,19 +47,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32
_data = (Complex32[])data.Clone();
}
public override Complex32 this[int index]
{
get
{
return _data[index];
}
set
{
_data[index] = value;
}
}
public override Matrix<Complex32> CreateMatrix(int rows, int columns)
{
return new UserDefinedMatrix(rows, columns);
@ -69,6 +56,16 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32
{
return new UserDefinedVector(size);
}
protected internal override Complex32 At(int index)
{
return _data[index];
}
protected internal override void At(int index, Complex32 value)
{
_data[index] = value;
}
}
public class UserDefinedVectorTests : VectorTests

23
src/UnitTests/LinearAlgebraTests/Double/UserDefinedVectorTests.cs

@ -46,19 +46,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
_data = (double[])data.Clone();
}
public override double this[int index]
{
get
{
return _data[index];
}
set
{
_data[index] = value;
}
}
public override Matrix<double> CreateMatrix(int rows, int columns)
{
return new UserDefinedMatrix(rows, columns);
@ -68,6 +55,16 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
{
return new UserDefinedVector(size);
}
protected internal override double At(int index)
{
return _data[index];
}
protected internal override void At(int index, double value)
{
_data[index] = value;
}
}
public class UserDefinedVectorTests : VectorTests

23
src/UnitTests/LinearAlgebraTests/Single/UserDefinedVectorTests.cs

@ -55,19 +55,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single
_data = (float[])data.Clone();
}
public override float this[int index]
{
get
{
return _data[index];
}
set
{
_data[index] = value;
}
}
public override Matrix<float> CreateMatrix(int rows, int columns)
{
return new UserDefinedMatrix(rows, columns);
@ -77,6 +64,16 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single
{
return new UserDefinedVector(size);
}
protected internal override float At(int index)
{
return _data[index];
}
protected internal override void At(int index, float value)
{
_data[index] = value;
}
}
public class UserDefinedVectorTests : VectorTests

Loading…
Cancel
Save