diff --git a/src/Benchmark/LinearAlgebra/SparseVector.cs b/src/Benchmark/LinearAlgebra/SparseVector.cs new file mode 100644 index 000000000..5e1f0bb8f --- /dev/null +++ b/src/Benchmark/LinearAlgebra/SparseVector.cs @@ -0,0 +1,64 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// https://numerics.mathdotnet.com +// https://github.com/mathnet/mathnet-numerics +// +// Copyright (c) 2009-$CURRENT_YEAR$ Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using BenchmarkDotNet.Attributes; +using MathNet.Numerics.Distributions; +using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.Random; +using System.Linq; + +namespace Benchmark.LinearAlgebra +{ + public class SparseVector + { + private Matrix matrix; + private Vector vector; + + [Params(100, 1000, 10000)] + public int N { get; set; } + + [Params(0.01, 0.1, 0.5)] + public double SparseRatio { get; set; } + + [GlobalSetup] + public void GlobalSetup() + { + const int Seed = 42; + + var rnd = new SystemRandomSource(Seed, true); + + matrix = Matrix.Build.Random(N, N, new Normal(rnd)); + + vector = Vector.Build.SparseOfEnumerable(rnd.NextDoubleSequence().Take(N).Select(x => x < SparseRatio ? x : 0)); + } + + [Benchmark] + public Vector DenseMatrixSparseVector() => matrix.Multiply(vector); + } +} diff --git a/src/Benchmark/Program.cs b/src/Benchmark/Program.cs index 25cd5e411..88814f9cb 100644 --- a/src/Benchmark/Program.cs +++ b/src/Benchmark/Program.cs @@ -16,6 +16,7 @@ public static void Main(string[] args) typeof(Transforms.FFT), typeof(LinearAlgebra.DenseMatrixProduct), typeof(LinearAlgebra.DenseVector), + typeof(LinearAlgebra.SparseVector), }); switcher.Run(args); diff --git a/src/Numerics.Tests/LinearAlgebraTests/Complex/MatrixTests.Arithmetic.cs b/src/Numerics.Tests/LinearAlgebraTests/Complex/MatrixTests.Arithmetic.cs index 6cb5cc3f8..aa5f78237 100644 --- a/src/Numerics.Tests/LinearAlgebraTests/Complex/MatrixTests.Arithmetic.cs +++ b/src/Numerics.Tests/LinearAlgebraTests/Complex/MatrixTests.Arithmetic.cs @@ -41,6 +41,28 @@ namespace MathNet.Numerics.Tests.LinearAlgebraTests.Complex /// public abstract partial class MatrixTests { + private static Vector GetVector(string name, Complex[] data) + { + switch (name) + { + case "Dense": return Vector.Build.Dense(data); + case "Sparse": return Vector.Build.SparseOfArray(data); + case "User": return new UserDefinedVector(data); + default: throw new NotImplementedException($"{nameof(GetVector)}(string, Complex[]) for {nameof(name)}=\"{name}\""); + } + } + + private static Vector GetVector(string name, int size) + { + switch (name) + { + case "Dense": return Vector.Build.Dense(size); + case "Sparse": return Vector.Build.Sparse(size); + case "User": return new UserDefinedVector(size); + default: throw new NotImplementedException($"{nameof(GetVector)}(string, int) for {nameof(name)}=\"{name}\""); + } + } + /// /// Can multiply with a complex number. /// @@ -68,10 +90,10 @@ public void CanMultiplyWithComplex(double real) /// Can multiply with a vector. /// [Test] - public void CanMultiplyWithVector() + public void CanMultiplyWithVector([Values("Dense", "Sparse", "User")] string vec) { var matrix = TestMatrices["Singular3x3"]; - var x = new DenseVector(new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); + var x = GetVector(vec, new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); var y = matrix * x; Assert.AreEqual(matrix.RowCount, y.Count); @@ -88,11 +110,13 @@ public void CanMultiplyWithVector() /// Can multiply with a vector into a result. /// [Test] - public void CanMultiplyWithVectorIntoResult() + public void CanMultiplyWithVectorIntoResult( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) { var matrix = TestMatrices["Singular3x3"]; - var x = new DenseVector(new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); - var y = new DenseVector(3); + var x = GetVector(vecX, new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); + var y = GetVector(vecY, 3); matrix.Multiply(x, y); for (var i = 0; i < matrix.RowCount; i++) @@ -107,16 +131,18 @@ public void CanMultiplyWithVectorIntoResult() /// Can multiply with a vector into result when updating input argument. /// [Test] - public void CanMultiplyWithVectorIntoResultWhenUpdatingInputArgument() + public void CanMultiplyWithVectorIntoResultWhenUpdatingInputArgument( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) { var matrix = TestMatrices["Singular3x3"]; - var x = new DenseVector(new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); + var x = GetVector(vecX, new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); var y = x; matrix.Multiply(x, x); Assert.AreSame(y, x); - y = new DenseVector(new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); + y = GetVector(vecY, new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); for (var i = 0; i < matrix.RowCount; i++) { var ar = matrix.Row(i); @@ -129,11 +155,13 @@ public void CanMultiplyWithVectorIntoResultWhenUpdatingInputArgument() /// Multiply with a vector into too large result throws ArgumentException. /// [Test] - public void MultiplyWithVectorIntoLargerResultThrowsArgumentException() + public void MultiplyWithVectorIntoLargerResultThrowsArgumentException( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) { var matrix = TestMatrices["Singular3x3"]; - var x = new DenseVector(new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); - Vector y = new DenseVector(4); + var x = GetVector(vecX, new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); + var y = GetVector(vecY, 4); Assert.That(() => matrix.Multiply(x, y), Throws.ArgumentException); } @@ -614,10 +642,10 @@ public virtual void CanMultiplyMatrixWithMatrixIntoResult(string nameA, string n /// Can multiply transposed matrix with a vector. /// [Test] - public void CanTransposeThisAndMultiplyWithVector() + public void CanTransposeThisAndMultiplyWithVector([Values("Dense", "Sparse", "User")] string vec) { var matrix = TestMatrices["Singular3x3"]; - var x = new DenseVector(new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); + var x = GetVector(vec, new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); var y = matrix.TransposeThisAndMultiply(x); Assert.AreEqual(matrix.ColumnCount, y.Count); @@ -634,11 +662,13 @@ public void CanTransposeThisAndMultiplyWithVector() /// Can multiply transposed matrix with a vector into a result. /// [Test] - public void CanTransposeThisAndMultiplyWithVectorIntoResult() + public void CanTransposeThisAndMultiplyWithVectorIntoResult( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) { var matrix = TestMatrices["Singular3x3"]; - var x = new DenseVector(new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); - var y = new DenseVector(3); + var x = GetVector(vecX, new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); + var y = GetVector(vecY, 3); matrix.TransposeThisAndMultiply(x, y); for (var j = 0; j < matrix.ColumnCount; j++) @@ -653,16 +683,18 @@ public void CanTransposeThisAndMultiplyWithVectorIntoResult() /// Can multiply transposed matrix with a vector into result when updating input argument. /// [Test] - public void CanTransposeThisAndMultiplyWithVectorIntoResultWhenUpdatingInputArgument() + public void CanTransposeThisAndMultiplyWithVectorIntoResultWhenUpdatingInputArgument( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) { var matrix = TestMatrices["Singular3x3"]; - var x = new DenseVector(new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); + var x = GetVector(vecX, new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); var y = x; matrix.TransposeThisAndMultiply(x, x); Assert.AreSame(y, x); - y = new DenseVector(new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); + y = GetVector(vecY, new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); for (var j = 0; j < matrix.ColumnCount; j++) { var ar = matrix.Column(j); @@ -675,14 +707,95 @@ public void CanTransposeThisAndMultiplyWithVectorIntoResultWhenUpdatingInputArgu /// Multiply transposed matrix with a vector into too large result throws ArgumentException. /// [Test] - public void TransposeThisAndMultiplyWithVectorIntoLargerResultThrowsArgumentException() + public void TransposeThisAndMultiplyWithVectorIntoLargerResultThrowsArgumentException( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) { var matrix = TestMatrices["Singular3x3"]; - var x = new DenseVector(new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); - Vector y = new DenseVector(4); + var x = GetVector(vecX, new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); + var y = GetVector(vecY, 4); Assert.That(() => matrix.TransposeThisAndMultiply(x, y), Throws.ArgumentException); } + /// + /// Can multiply conjugate transposed matrix with a vector. + /// + [Test] + public void CanConjugateTransposeThisAndMultiplyWithVector([Values("Dense", "Sparse", "User")] string vec) + { + var matrix = TestMatrices["Singular3x3"]; + var x = GetVector(vec, new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); + var y = matrix.ConjugateTransposeThisAndMultiply(x); + + Assert.AreEqual(matrix.ColumnCount, y.Count); + + for (var j = 0; j < matrix.ColumnCount; j++) + { + var ar = matrix.Column(j); + var dot = ar.ConjugateDotProduct(x); + Assert.AreEqual(dot, y[j]); + } + } + + /// + /// Can multiply conjugate transposed matrix with a vector into a result. + /// + [Test] + public void CanConjugateTransposeThisAndMultiplyWithVectorIntoResult( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) + { + var matrix = TestMatrices["Singular3x3"]; + var x = GetVector(vecX, new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); + var y = GetVector(vecY, 3); + matrix.ConjugateTransposeThisAndMultiply(x, y); + + for (var j = 0; j < matrix.ColumnCount; j++) + { + var ar = matrix.Column(j); + var dot = ar.ConjugateDotProduct(x); + Assert.AreEqual(dot, y[j]); + } + } + + /// + /// Can multiply conjugate transposed matrix with a vector into result when updating input argument. + /// + [Test] + public void CanConjugateTransposeThisAndMultiplyWithVectorIntoResultWhenUpdatingInputArgument( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) + { + var matrix = TestMatrices["Singular3x3"]; + var x = GetVector(vecX, new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); + var y = x; + matrix.ConjugateTransposeThisAndMultiply(x, x); + + Assert.AreSame(y, x); + + y = GetVector(vecY, new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); + for (var j = 0; j < matrix.ColumnCount; j++) + { + var ar = matrix.Column(j); + var dot = ar.ConjugateDotProduct(y); + Assert.AreEqual(dot, x[j]); + } + } + + /// + /// Multiply conjugate transposed matrix with a vector into too large result throws ArgumentException. + /// + [Test] + public void ConjugateTransposeThisAndMultiplyWithVectorIntoLargerResultThrowsArgumentException( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) + { + var matrix = TestMatrices["Singular3x3"]; + var x = GetVector(vecX, new[] { new Complex(1, 1), new Complex(2, 1), new Complex(3, 1) }); + var y = GetVector(vecY, 4); + Assert.That(() => matrix.ConjugateTransposeThisAndMultiply(x, y), Throws.ArgumentException); + } + /// /// Can multiply transposed matrix with another matrix. /// diff --git a/src/Numerics.Tests/LinearAlgebraTests/Complex32/MatrixTests.Arithmetic.cs b/src/Numerics.Tests/LinearAlgebraTests/Complex32/MatrixTests.Arithmetic.cs index 45d51fa32..a5faaa75a 100644 --- a/src/Numerics.Tests/LinearAlgebraTests/Complex32/MatrixTests.Arithmetic.cs +++ b/src/Numerics.Tests/LinearAlgebraTests/Complex32/MatrixTests.Arithmetic.cs @@ -42,6 +42,28 @@ namespace MathNet.Numerics.Tests.LinearAlgebraTests.Complex32 /// public abstract partial class MatrixTests { + private static Vector GetVector(string name, Complex32[] data) + { + switch (name) + { + case "Dense": return Vector.Build.Dense(data); + case "Sparse": return Vector.Build.SparseOfArray(data); + case "User": return new UserDefinedVector(data); + default: throw new NotImplementedException($"{nameof(GetVector)}(string, Complex32[]) for {nameof(name)}=\"{name}\""); + } + } + + private static Vector GetVector(string name, int size) + { + switch (name) + { + case "Dense": return Vector.Build.Dense(size); + case "Sparse": return Vector.Build.Sparse(size); + case "User": return new UserDefinedVector(size); + default: throw new NotImplementedException($"{nameof(GetVector)}(string, int) for {nameof(name)}=\"{name}\""); + } + } + /// /// Can multiply with a complex number. /// @@ -69,10 +91,10 @@ public void CanMultiplyWithComplex(float real) /// Can multiply with a vector. /// [Test] - public void CanMultiplyWithVector() + public void CanMultiplyWithVector([Values("Dense", "Sparse", "User")] string vec) { var matrix = TestMatrices["Singular3x3"]; - var x = new DenseVector(new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); + var x = GetVector(vec, new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); var y = matrix * x; Assert.AreEqual(matrix.RowCount, y.Count); @@ -89,11 +111,13 @@ public void CanMultiplyWithVector() /// Can multiply with a vector into a result. /// [Test] - public void CanMultiplyWithVectorIntoResult() + public void CanMultiplyWithVectorIntoResult( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) { var matrix = TestMatrices["Singular3x3"]; - var x = new DenseVector(new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); - var y = new DenseVector(3); + var x = GetVector(vecX, new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); + var y = GetVector(vecY, 3); matrix.Multiply(x, y); for (var i = 0; i < matrix.RowCount; i++) @@ -108,16 +132,18 @@ public void CanMultiplyWithVectorIntoResult() /// Can multiply with a vector into result when updating input argument. /// [Test] - public void CanMultiplyWithVectorIntoResultWhenUpdatingInputArgument() + public void CanMultiplyWithVectorIntoResultWhenUpdatingInputArgument( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) { var matrix = TestMatrices["Singular3x3"]; - var x = new DenseVector(new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); + var x = GetVector(vecX, new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); var y = x; matrix.Multiply(x, x); Assert.AreSame(y, x); - y = new DenseVector(new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); + y = GetVector(vecY, new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); for (var i = 0; i < matrix.RowCount; i++) { var ar = matrix.Row(i); @@ -130,11 +156,13 @@ public void CanMultiplyWithVectorIntoResultWhenUpdatingInputArgument() /// Multiply with a vector into too large result throws ArgumentException. /// [Test] - public void MultiplyWithVectorIntoLargerResultThrowsArgumentException() + public void MultiplyWithVectorIntoLargerResultThrowsArgumentException( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) { var matrix = TestMatrices["Singular3x3"]; - var x = new DenseVector(new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); - Vector y = new DenseVector(4); + var x = GetVector(vecX, new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); + var y = GetVector(vecY, 4); Assert.That(() => matrix.Multiply(x, y), Throws.ArgumentException); } @@ -615,10 +643,10 @@ public virtual void CanMultiplyMatrixWithMatrixIntoResult(string nameA, string n /// Can multiply transposed matrix with a vector. /// [Test] - public void CanTransposeThisAndMultiplyWithVector() + public void CanTransposeThisAndMultiplyWithVector([Values("Dense", "Sparse", "User")] string vec) { var matrix = TestMatrices["Singular3x3"]; - var x = new DenseVector(new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); + var x = GetVector(vec, new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); var y = matrix.TransposeThisAndMultiply(x); Assert.AreEqual(matrix.ColumnCount, y.Count); @@ -635,11 +663,13 @@ public void CanTransposeThisAndMultiplyWithVector() /// Can multiply transposed matrix with a vector into a result. /// [Test] - public void CanTransposeThisAndMultiplyWithVectorIntoResult() + public void CanTransposeThisAndMultiplyWithVectorIntoResult( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) { var matrix = TestMatrices["Singular3x3"]; - var x = new DenseVector(new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); - var y = new DenseVector(3); + var x = GetVector(vecX, new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); + var y = GetVector(vecY, 3); matrix.TransposeThisAndMultiply(x, y); for (var j = 0; j < matrix.ColumnCount; j++) @@ -654,16 +684,18 @@ public void CanTransposeThisAndMultiplyWithVectorIntoResult() /// Can multiply transposed matrix with a vector into result when updating input argument. /// [Test] - public void CanTransposeThisAndMultiplyWithVectorIntoResultWhenUpdatingInputArgument() + public void CanTransposeThisAndMultiplyWithVectorIntoResultWhenUpdatingInputArgument( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) { var matrix = TestMatrices["Singular3x3"]; - var x = new DenseVector(new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); + var x = GetVector(vecX, new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); var y = x; matrix.TransposeThisAndMultiply(x, x); Assert.AreSame(y, x); - y = new DenseVector(new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); + y = GetVector(vecY, new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); for (var j = 0; j < matrix.ColumnCount; j++) { var ar = matrix.Column(j); @@ -676,14 +708,95 @@ public void CanTransposeThisAndMultiplyWithVectorIntoResultWhenUpdatingInputArgu /// Multiply transposed matrix with a vector into too large result throws ArgumentException. /// [Test] - public void TransposeThisAndMultiplyWithVectorIntoLargerResultThrowsArgumentException() + public void TransposeThisAndMultiplyWithVectorIntoLargerResultThrowsArgumentException( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) { var matrix = TestMatrices["Singular3x3"]; - var x = new DenseVector(new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); - Vector y = new DenseVector(4); + var x = GetVector(vecX, new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); + var y = GetVector(vecY, 4); Assert.That(() => matrix.TransposeThisAndMultiply(x, y), Throws.ArgumentException); } + /// + /// Can multiply conjugate transposed matrix with a vector. + /// + [Test] + public void CanConjugateTransposeThisAndMultiplyWithVector([Values("Dense", "Sparse", "User")] string vec) + { + var matrix = TestMatrices["Singular3x3"]; + var x = GetVector(vec, new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); + var y = matrix.ConjugateTransposeThisAndMultiply(x); + + Assert.AreEqual(matrix.ColumnCount, y.Count); + + for (var j = 0; j < matrix.ColumnCount; j++) + { + var ar = matrix.Column(j); + var dot = ar.ConjugateDotProduct(x); + AssertHelpers.AlmostEqual(dot, y[j], 5); + } + } + + /// + /// Can multiply conjugate transposed matrix with a vector into a result. + /// + [Test] + public void CanConjugateTransposeThisAndMultiplyWithVectorIntoResult( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) + { + var matrix = TestMatrices["Singular3x3"]; + var x = GetVector(vecX, new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); + var y = GetVector(vecY, 3); + matrix.ConjugateTransposeThisAndMultiply(x, y); + + for (var j = 0; j < matrix.ColumnCount; j++) + { + var ar = matrix.Column(j); + var dot = ar.ConjugateDotProduct(x); + AssertHelpers.AlmostEqual(dot, y[j], 5); + } + } + + /// + /// Can multiply conjugate transposed matrix with a vector into result when updating input argument. + /// + [Test] + public void CanConjugateTransposeThisAndMultiplyWithVectorIntoResultWhenUpdatingInputArgument( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) + { + var matrix = TestMatrices["Singular3x3"]; + var x = GetVector(vecX, new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); + var y = x; + matrix.ConjugateTransposeThisAndMultiply(x, x); + + Assert.AreSame(y, x); + + y = GetVector(vecY, new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); + for (var j = 0; j < matrix.ColumnCount; j++) + { + var ar = matrix.Column(j); + var dot = ar.ConjugateDotProduct(y); + AssertHelpers.AlmostEqual(dot, x[j], 5); + } + } + + /// + /// Multiply conjugate transposed matrix with a vector into too large result throws ArgumentException. + /// + [Test] + public void ConjugateTransposeThisAndMultiplyWithVectorIntoLargerResultThrowsArgumentException( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) + { + var matrix = TestMatrices["Singular3x3"]; + var x = GetVector(vecX, new[] { new Complex32(1, 1), new Complex32(2, 1), new Complex32(3, 1) }); + var y = GetVector(vecY, 4); + Assert.That(() => matrix.ConjugateTransposeThisAndMultiply(x, y), Throws.ArgumentException); + } + /// /// Can multiply transposed matrix with another matrix. /// diff --git a/src/Numerics.Tests/LinearAlgebraTests/Double/MatrixTests.Arithmetic.cs b/src/Numerics.Tests/LinearAlgebraTests/Double/MatrixTests.Arithmetic.cs index 33fff02d5..8a38d2706 100644 --- a/src/Numerics.Tests/LinearAlgebraTests/Double/MatrixTests.Arithmetic.cs +++ b/src/Numerics.Tests/LinearAlgebraTests/Double/MatrixTests.Arithmetic.cs @@ -39,6 +39,28 @@ namespace MathNet.Numerics.Tests.LinearAlgebraTests.Double /// public abstract partial class MatrixTests { + private static Vector GetVector(string name, double[] data) + { + switch (name) + { + case "Dense": return Vector.Build.Dense(data); + case "Sparse": return Vector.Build.SparseOfArray(data); + case "User": return new UserDefinedVector(data); + default: throw new NotImplementedException($"{nameof(GetVector)}(string, double[]) for {nameof(name)}=\"{name}\""); + } + } + + private static Vector GetVector(string name, int size) + { + switch (name) + { + case "Dense": return Vector.Build.Dense(size); + case "Sparse": return Vector.Build.Sparse(size); + case "User": return new UserDefinedVector(size); + default: throw new NotImplementedException($"{nameof(GetVector)}(string, int) for {nameof(name)}=\"{name}\""); + } + } + /// /// Can multiply with a scalar. /// @@ -65,10 +87,10 @@ public void CanMultiplyWithScalar(double scalar) /// Can multiply with a vector. /// [Test] - public void CanMultiplyWithVector() + public void CanMultiplyWithVector([Values("Dense", "Sparse", "User")] string vec) { var matrix = TestMatrices["Singular3x3"]; - var x = Vector.Build.Dense(new[] { 1.0, 2.0, 3.0 }); + var x = GetVector(vec, new[] { 1.0, 2.0, 3.0 }); var y = matrix*x; Assert.AreEqual(matrix.RowCount, y.Count); @@ -85,11 +107,13 @@ public void CanMultiplyWithVector() /// Can multiply with a vector into a result. /// [Test] - public void CanMultiplyWithVectorIntoResult() + public void CanMultiplyWithVectorIntoResult( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) { var matrix = TestMatrices["Singular3x3"]; - var x = Vector.Build.Dense(new[] { 1.0, 2.0, 3.0 }); - var y = Vector.Build.Dense(3); + var x = GetVector(vecX, new[] { 1.0, 2.0, 3.0 }); + var y = GetVector(vecY, 3); matrix.Multiply(x, y); for (var i = 0; i < matrix.RowCount; i++) @@ -104,16 +128,18 @@ public void CanMultiplyWithVectorIntoResult() /// Can multiply with a vector into result when updating input argument. /// [Test] - public void CanMultiplyWithVectorIntoResultWhenUpdatingInputArgument() + public void CanMultiplyWithVectorIntoResultWhenUpdatingInputArgument( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) { var matrix = TestMatrices["Singular3x3"]; - var x = Vector.Build.Dense(new[] { 1.0, 2.0, 3.0 }); + var x = GetVector(vecX, new[] { 1.0, 2.0, 3.0 }); var y = x; matrix.Multiply(x, x); Assert.AreSame(y, x); - y = Vector.Build.Dense(new[] { 1.0, 2.0, 3.0 }); + y = GetVector(vecY, new[] { 1.0, 2.0, 3.0 }); for (var i = 0; i < matrix.RowCount; i++) { var ar = matrix.Row(i); @@ -126,11 +152,13 @@ public void CanMultiplyWithVectorIntoResultWhenUpdatingInputArgument() /// Multiply with a vector into too large result throws ArgumentException. /// [Test] - public void MultiplyWithVectorIntoLargerResultThrowsArgumentException() + public void MultiplyWithVectorIntoLargerResultThrowsArgumentException( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) { var matrix = TestMatrices["Singular3x3"]; - var x = Vector.Build.Dense(new[] { 1.0, 2.0, 3.0 }); - Vector y = Vector.Build.Dense(4); + var x = GetVector(vecX, new[] { 1.0, 2.0, 3.0 }); + var y = GetVector(vecY, 4); Assert.That(() => matrix.Multiply(x, y), Throws.ArgumentException); } @@ -608,10 +636,10 @@ public virtual void CanMultiplyMatrixWithMatrixIntoResult(string nameA, string n /// Can multiply transposed matrix with a vector. /// [Test] - public void CanTransposeThisAndMultiplyWithVector() + public void CanTransposeThisAndMultiplyWithVector([Values("Dense", "Sparse", "User")] string vec) { var matrix = TestMatrices["Singular3x3"]; - var x = Vector.Build.Dense(new[] { 1.0, 2.0, 3.0 }); + var x = GetVector(vec, new[] { 1.0, 2.0, 3.0 }); var y = matrix.TransposeThisAndMultiply(x); Assert.AreEqual(matrix.ColumnCount, y.Count); @@ -628,11 +656,13 @@ public void CanTransposeThisAndMultiplyWithVector() /// Can multiply transposed matrix with a vector into a result. /// [Test] - public void CanTransposeThisAndMultiplyWithVectorIntoResult() + public void CanTransposeThisAndMultiplyWithVectorIntoResult( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) { var matrix = TestMatrices["Singular3x3"]; - var x = Vector.Build.Dense(new[] { 1.0, 2.0, 3.0 }); - var y = Vector.Build.Dense(3); + var x = GetVector(vecX, new[] { 1.0, 2.0, 3.0 }); + var y = GetVector(vecY, 3); matrix.TransposeThisAndMultiply(x, y); for (var j = 0; j < matrix.ColumnCount; j++) @@ -647,16 +677,18 @@ public void CanTransposeThisAndMultiplyWithVectorIntoResult() /// Can multiply transposed matrix with a vector into result when updating input argument. /// [Test] - public void CanTransposeThisAndMultiplyWithVectorIntoResultWhenUpdatingInputArgument() + public void CanTransposeThisAndMultiplyWithVectorIntoResultWhenUpdatingInputArgument( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) { var matrix = TestMatrices["Singular3x3"]; - var x = Vector.Build.Dense(new[] { 1.0, 2.0, 3.0 }); + var x = GetVector(vecX, new[] { 1.0, 2.0, 3.0 }); var y = x; matrix.TransposeThisAndMultiply(x, x); Assert.AreSame(y, x); - y = Vector.Build.Dense(new[] { 1.0, 2.0, 3.0 }); + y = GetVector(vecY, new[] { 1.0, 2.0, 3.0 }); for (var j = 0; j < matrix.ColumnCount; j++) { var ar = matrix.Column(j); @@ -669,11 +701,13 @@ public void CanTransposeThisAndMultiplyWithVectorIntoResultWhenUpdatingInputArgu /// Multiply transposed matrix with a vector into too large result throws ArgumentException. /// [Test] - public void TransposeThisAndMultiplyWithVectorIntoLargerResultThrowsArgumentException() + public void TransposeThisAndMultiplyWithVectorIntoLargerResultThrowsArgumentException( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) { var matrix = TestMatrices["Singular3x3"]; - var x = Vector.Build.Dense(new[] { 1.0, 2.0, 3.0 }); - Vector y = Vector.Build.Dense(4); + var x = GetVector(vecX, new[] { 1.0, 2.0, 3.0 }); + var y = GetVector(vecY, 4); Assert.That(() => matrix.TransposeThisAndMultiply(x, y), Throws.ArgumentException); } diff --git a/src/Numerics.Tests/LinearAlgebraTests/MatrixTests.cs b/src/Numerics.Tests/LinearAlgebraTests/MatrixTests.cs index d461d9c35..838226b33 100644 --- a/src/Numerics.Tests/LinearAlgebraTests/MatrixTests.cs +++ b/src/Numerics.Tests/LinearAlgebraTests/MatrixTests.cs @@ -116,5 +116,52 @@ public void PointwiseDivision_SparseReturnsSameResultAsDenseTest() var result = SparseMatrix.OfDiagonalArray(new double[] { 10, 5, 2 }); Assert.AreEqual(result, x); } + + [Test] + public void MatrixVectorMultiplicationSameResult( + [Values("Dense", "Sparse")] string mtx, + [Values("Dense", "Sparse")] string vec, + [Values(0, 1, -1, double.Epsilon, double.MaxValue, double.MinValue, double.NaN, double.NegativeInfinity, double.PositiveInfinity)] double m_01, + [Values(0, 1, -1, double.Epsilon, double.MaxValue, double.MinValue, double.NaN, double.NegativeInfinity, double.PositiveInfinity)] double v_1 + ) + { + double m_00 = 1; + double m_10 = 0; + double m_11 = -2; + + double v_0 = 0; + + var matrixValues = new double[,] { { m_00, m_01 }, { m_10, m_11 } }; + var vectorValues = new double[] { v_0, v_1 }; + + var m = GetMatrix(mtx, matrixValues); + var v = GetVector(vec, vectorValues); + + var result = m * v; + + Assert.AreEqual(2, result.Count, "Bad result dimension"); + Assert.AreEqual((m_00 * v_0) + (m_01 * v_1), result[0], "Bad result index 0"); + Assert.AreEqual((m_10 * v_0) + (m_11 * v_1), result[1], "Bad result index 1"); + } + + private static Matrix GetMatrix(string name, double[,] data) + { + switch (name) + { + case "Dense": return Matrix.Build.DenseOfArray(data); + case "Sparse": return Matrix.Build.SparseOfArray(data); + default: throw new NotImplementedException($"{nameof(GetMatrix)}(string, double[,]) for {nameof(name)}=\"{name}\""); + } + } + + private static Vector GetVector(string name, double[] data) + { + switch (name) + { + case "Dense": return Vector.Build.Dense(data); + case "Sparse": return Vector.Build.SparseOfArray(data); + default: throw new NotImplementedException($"{nameof(GetVector)}(string, double[]) for {nameof(name)}=\"{name}\""); + } + } } } diff --git a/src/Numerics.Tests/LinearAlgebraTests/Single/MatrixTests.Arithmetic.cs b/src/Numerics.Tests/LinearAlgebraTests/Single/MatrixTests.Arithmetic.cs index 494a05d95..7953a22c4 100644 --- a/src/Numerics.Tests/LinearAlgebraTests/Single/MatrixTests.Arithmetic.cs +++ b/src/Numerics.Tests/LinearAlgebraTests/Single/MatrixTests.Arithmetic.cs @@ -40,6 +40,28 @@ namespace MathNet.Numerics.Tests.LinearAlgebraTests.Single /// public abstract partial class MatrixTests { + private static Vector GetVector(string name, float[] data) + { + switch (name) + { + case "Dense": return Vector.Build.Dense(data); + case "Sparse": return Vector.Build.SparseOfArray(data); + case "User": return new UserDefinedVector(data); + default: throw new NotImplementedException($"{nameof(GetVector)}(string, float[]) for {nameof(name)}=\"{name}\""); + } + } + + private static Vector GetVector(string name, int size) + { + switch (name) + { + case "Dense": return Vector.Build.Dense(size); + case "Sparse": return Vector.Build.Sparse(size); + case "User": return new UserDefinedVector(size); + default: throw new NotImplementedException($"{nameof(GetVector)}(string, int) for {nameof(name)}=\"{name}\""); + } + } + /// /// Can multiply with a scalar. /// @@ -66,10 +88,10 @@ public void CanMultiplyWithScalar(float scalar) /// Can multiply with a vector. /// [Test] - public void CanMultiplyWithVector() + public void CanMultiplyWithVector([Values("Dense", "Sparse", "User")] string vec) { var matrix = TestMatrices["Singular3x3"]; - var x = new DenseVector(new[] { 1.0f, 2.0f, 3.0f }); + var x = GetVector(vec, new[] { 1.0f, 2.0f, 3.0f }); var y = matrix * x; Assert.AreEqual(matrix.RowCount, y.Count); @@ -86,11 +108,13 @@ public void CanMultiplyWithVector() /// Can multiply with a vector into a result. /// [Test] - public void CanMultiplyWithVectorIntoResult() + public void CanMultiplyWithVectorIntoResult( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) { var matrix = TestMatrices["Singular3x3"]; - var x = new DenseVector(new[] { 1.0f, 2.0f, 3.0f }); - var y = new DenseVector(3); + var x = GetVector(vecX, new[] { 1.0f, 2.0f, 3.0f }); + var y = GetVector(vecY, 3); matrix.Multiply(x, y); for (var i = 0; i < matrix.RowCount; i++) @@ -105,16 +129,18 @@ public void CanMultiplyWithVectorIntoResult() /// Can multiply with a vector into result when updating input argument. /// [Test] - public void CanMultiplyWithVectorIntoResultWhenUpdatingInputArgument() + public void CanMultiplyWithVectorIntoResultWhenUpdatingInputArgument( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) { var matrix = TestMatrices["Singular3x3"]; - var x = new DenseVector(new[] { 1.0f, 2.0f, 3.0f }); + var x = GetVector(vecX, new[] { 1.0f, 2.0f, 3.0f }); var y = x; matrix.Multiply(x, x); Assert.AreSame(y, x); - y = new DenseVector(new[] { 1.0f, 2.0f, 3.0f }); + y = GetVector(vecY, new[] { 1.0f, 2.0f, 3.0f }); for (var i = 0; i < matrix.RowCount; i++) { var ar = matrix.Row(i); @@ -127,11 +153,13 @@ public void CanMultiplyWithVectorIntoResultWhenUpdatingInputArgument() /// Multiply with a vector into too large result throws ArgumentException. /// [Test] - public void MultiplyWithVectorIntoLargerResultThrowsArgumentException() + public void MultiplyWithVectorIntoLargerResultThrowsArgumentException( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) { var matrix = TestMatrices["Singular3x3"]; - var x = new DenseVector(new[] { 1.0f, 2.0f, 3.0f }); - Vector y = new DenseVector(4); + var x = GetVector(vecX, new[] { 1.0f, 2.0f, 3.0f }); + var y = GetVector(vecY, 4); Assert.That(() => matrix.Multiply(x, y), Throws.ArgumentException); } @@ -609,10 +637,10 @@ public virtual void CanMultiplyMatrixWithMatrixIntoResult(string nameA, string n /// Can multiply transposed matrix with a vector. /// [Test] - public void CanTransposeThisAndMultiplyWithVector() + public void CanTransposeThisAndMultiplyWithVector([Values("Dense", "Sparse", "User")] string vec) { var matrix = TestMatrices["Singular3x3"]; - var x = new DenseVector(new[] { 1.0f, 2.0f, 3.0f }); + var x = GetVector(vec, new[] { 1.0f, 2.0f, 3.0f }); var y = matrix.TransposeThisAndMultiply(x); Assert.AreEqual(matrix.ColumnCount, y.Count); @@ -627,11 +655,13 @@ public void CanTransposeThisAndMultiplyWithVector() /// Can multiply transposed matrix with a vector into a result. /// [Test] - public void CanTransposeThisAndMultiplyWithVectorIntoResult() + public void CanTransposeThisAndMultiplyWithVectorIntoResult( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) { var matrix = TestMatrices["Singular3x3"]; - var x = new DenseVector(new[] { 1.0f, 2.0f, 3.0f }); - var y = new DenseVector(3); + var x = GetVector(vecX, new[] { 1.0f, 2.0f, 3.0f }); + var y = GetVector(vecY, 3); matrix.TransposeThisAndMultiply(x, y); for (var j = 0; j < matrix.ColumnCount; j++) @@ -644,16 +674,18 @@ public void CanTransposeThisAndMultiplyWithVectorIntoResult() /// Can multiply transposed matrix with a vector into result when updating input argument. /// [Test] - public void CanTransposeThisAndMultiplyWithVectorIntoResultWhenUpdatingInputArgument() + public void CanTransposeThisAndMultiplyWithVectorIntoResultWhenUpdatingInputArgument( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) { var matrix = TestMatrices["Singular3x3"]; - var x = new DenseVector(new[] { 1.0f, 2.0f, 3.0f }); + var x = GetVector(vecX, new[] { 1.0f, 2.0f, 3.0f }); var y = x; matrix.TransposeThisAndMultiply(x, x); Assert.AreSame(y, x); - y = new DenseVector(new[] { 1.0f, 2.0f, 3.0f }); + y = GetVector(vecY, new[] { 1.0f, 2.0f, 3.0f }); for (var j = 0; j < matrix.ColumnCount; j++) { AssertHelpers.AlmostEqual(matrix.Column(j) * y, x[j], 6); @@ -664,11 +696,13 @@ public void CanTransposeThisAndMultiplyWithVectorIntoResultWhenUpdatingInputArgu /// Multiply transposed matrix with a vector into too large result throws ArgumentException. /// [Test] - public void TransposeThisAndMultiplyWithVectorIntoLargerResultThrowsArgumentException() + public void TransposeThisAndMultiplyWithVectorIntoLargerResultThrowsArgumentException( + [Values("Dense", "Sparse", "User")] string vecX, + [Values("Dense", "Sparse", "User")] string vecY) { var matrix = TestMatrices["Singular3x3"]; - var x = new DenseVector(new[] { 1.0f, 2.0f, 3.0f }); - Vector y = new DenseVector(4); + var x = GetVector(vecX, new[] { 1.0f, 2.0f, 3.0f }); + var y = GetVector(vecY, 4); Assert.That(() => matrix.TransposeThisAndMultiply(x, y), Throws.ArgumentException); } diff --git a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs index 6dca6191e..cffa14974 100644 --- a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs @@ -598,11 +598,24 @@ protected override void DoMultiply(Vector rightSide, Vector re denseRight.Count, 1, denseResult.Values); + return; } - else + + if (rightSide.Storage is SparseVectorStorage sparseRight) { - base.DoMultiply(rightSide, result); + for (var i = 0; i < RowCount; i++) + { + var s = Complex.Zero; + for (var j = 0; j < sparseRight.ValueCount; j++) + { + s += At(i, sparseRight.Indices[j])*sparseRight.Values[j]; + } + result[i] = s; + } + return; } + + base.DoMultiply(rightSide, result); } /// @@ -772,6 +785,20 @@ protected override void DoTransposeThisAndMultiply(Vector rightSide, Ve return; } + if (rightSide.Storage is SparseVectorStorage sparseRight) + { + for (var j = 0; j < ColumnCount; j++) + { + var s = Complex.Zero; + for (var i = 0; i < sparseRight.ValueCount; i++) + { + s += At(sparseRight.Indices[i], j)*sparseRight.Values[i]; + } + result[j] = s; + } + return; + } + base.DoTransposeThisAndMultiply(rightSide, result); } @@ -796,6 +823,19 @@ protected override void DoConjugateTransposeThisAndMultiply(Vector righ 1, 0.0, denseResult.Values); + } + + if (rightSide.Storage is SparseVectorStorage sparseRight) + { + for (var j = 0; j < ColumnCount; j++) + { + var s = Complex.Zero; + for (var i = 0; i < sparseRight.ValueCount; i++) + { + s += At(sparseRight.Indices[i], j).Conjugate()*sparseRight.Values[i]; + } + result[j] = s; + } return; } diff --git a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs index ed0d8ca96..8b8099b4a 100644 --- a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs @@ -598,11 +598,24 @@ protected override void DoMultiply(Vector rightSide, Vector sparseRight) { - base.DoMultiply(rightSide, result); + for (var i = 0; i < RowCount; i++) + { + var s = Complex32.Zero; + for (var j = 0; j < sparseRight.ValueCount; j++) + { + s += At(i, sparseRight.Indices[j])*sparseRight.Values[j]; + } + result[i] = s; + } + return; } + + base.DoMultiply(rightSide, result); } /// @@ -769,11 +782,24 @@ protected override void DoTransposeThisAndMultiply(Vector rightSide, 1, 0.0f, denseResult.Values); + return; } - else + + if (rightSide.Storage is SparseVectorStorage sparseRight) { - base.DoTransposeThisAndMultiply(rightSide, result); + for (var j = 0; j < ColumnCount; j++) + { + var s = Complex32.Zero; + for (var i = 0; i < sparseRight.ValueCount; i++) + { + s += At(sparseRight.Indices[i], j)*sparseRight.Values[i]; + } + result[j] = s; + } + return; } + + base.DoTransposeThisAndMultiply(rightSide, result); } /// @@ -800,6 +826,20 @@ protected override void DoConjugateTransposeThisAndMultiply(Vector ri return; } + if (rightSide.Storage is SparseVectorStorage sparseRight) + { + for (var j = 0; j < ColumnCount; j++) + { + var s = Complex32.Zero; + for (var i = 0; i < sparseRight.ValueCount; i++) + { + s += At(sparseRight.Indices[i], j).Conjugate()*sparseRight.Values[i]; + } + result[j] = s; + } + return; + } + base.DoConjugateTransposeThisAndMultiply(rightSide, result); } diff --git a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs index c5bd91e64..a4c795af8 100644 --- a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs @@ -581,11 +581,24 @@ protected override void DoMultiply(Vector rightSide, Vector resu denseRight.Count, 1, denseResult.Values); + return; } - else + + if (rightSide.Storage is SparseVectorStorage sparseRight) { - base.DoMultiply(rightSide, result); + for (var i = 0; i < RowCount; i++) + { + var s = 0.0; + for (var j = 0; j < sparseRight.ValueCount; j++) + { + s += At(i, sparseRight.Indices[j])*sparseRight.Values[j]; + } + result[i] = s; + } + return; } + + base.DoMultiply(rightSide, result); } /// @@ -699,11 +712,24 @@ protected override void DoTransposeThisAndMultiply(Vector rightSide, Vec 1, 0.0, denseResult.Values); + return; } - else + + if (rightSide.Storage is SparseVectorStorage sparseRight) { - base.DoTransposeThisAndMultiply(rightSide, result); + for (var j = 0; j < ColumnCount; j++) + { + var s = 0.0; + for (var i = 0; i < sparseRight.ValueCount; i++) + { + s += At(sparseRight.Indices[i], j)*sparseRight.Values[i]; + } + result[j] = s; + } + return; } + + base.DoTransposeThisAndMultiply(rightSide, result); } /// diff --git a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs index 9f45cf212..2327299cd 100644 --- a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs @@ -581,11 +581,24 @@ protected override void DoMultiply(Vector rightSide, Vector result denseRight.Count, 1, denseResult.Values); + return; } - else + + if (rightSide.Storage is SparseVectorStorage sparseRight) { - base.DoMultiply(rightSide, result); + for (var i = 0; i < RowCount; i++) + { + var s = 0.0f; + for (var j = 0; j < sparseRight.ValueCount; j++) + { + s += At(i, sparseRight.Indices[j])*sparseRight.Values[j]; + } + result[i] = s; + } + return; } + + base.DoMultiply(rightSide, result); } /// @@ -699,11 +712,24 @@ protected override void DoTransposeThisAndMultiply(Vector rightSide, Vect 1, 0.0f, denseResult.Values); + return; } - else + + if (rightSide.Storage is SparseVectorStorage sparseRight) { - base.DoTransposeThisAndMultiply(rightSide, result); + for (var j = 0; j < ColumnCount; j++) + { + var s = 0.0f; + for (var i = 0; i < sparseRight.ValueCount; i++) + { + s += At(sparseRight.Indices[i], j)*sparseRight.Values[i]; + } + result[j] = s; + } + return; } + + base.DoTransposeThisAndMultiply(rightSide, result); } ///