diff --git a/MathNet.Numerics.sln.DotSettings b/MathNet.Numerics.sln.DotSettings index 569d8949..10562cdc 100644 --- a/MathNet.Numerics.sln.DotSettings +++ b/MathNet.Numerics.sln.DotSettings @@ -54,7 +54,9 @@ 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. </copyright> + BLAS CDF + CUDA DFT FFT ILU diff --git a/src/Numerics.Tests/IntegralTransformsTests/MatchingReferenceTransformTest.cs b/src/Numerics.Tests/IntegralTransformsTests/MatchingReferenceTransformTest.cs index 7d2bca84..f2e35df2 100644 --- a/src/Numerics.Tests/IntegralTransformsTests/MatchingReferenceTransformTest.cs +++ b/src/Numerics.Tests/IntegralTransformsTests/MatchingReferenceTransformTest.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2018 Math.NET +// Copyright (c) 2009-2021 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -321,7 +321,7 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests { // 65536 = 2^16 var samples = Generate.RandomComplex32(65536, GetUniform(1)); - Verify(samples, 5, FourierTransformScaling.NoScaling, FourierTransformControl.CreateManaged().Forward, FourierTransformControl.Provider.Forward); + Verify(samples, 5, FourierTransformScaling.NoScaling, ManagedFourierTransformProvider.Instance.Forward, FourierTransformControl.Provider.Forward); } [Test] @@ -329,7 +329,7 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests { // 65536 = 2^16 var samples = Generate.RandomComplex(65536, GetUniform(1)); - Verify(samples, 10, FourierTransformScaling.NoScaling, FourierTransformControl.CreateManaged().Forward, FourierTransformControl.Provider.Forward); + Verify(samples, 10, FourierTransformScaling.NoScaling, ManagedFourierTransformProvider.Instance.Forward, FourierTransformControl.Provider.Forward); } [Test] @@ -337,7 +337,7 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests { // 30870 = 2*3*3*5*7*7*7 var samples = Generate.RandomComplex32(30870, GetUniform(1)); - Verify(samples, 5, FourierTransformScaling.NoScaling, FourierTransformControl.CreateManaged().Forward, FourierTransformControl.Provider.Forward); + Verify(samples, 5, FourierTransformScaling.NoScaling, ManagedFourierTransformProvider.Instance.Forward, FourierTransformControl.Provider.Forward); } [Test] @@ -345,21 +345,21 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests { // 30870 = 2*3*3*5*7*7*7 var samples = Generate.RandomComplex(30870, GetUniform(1)); - Verify(samples, 10, FourierTransformScaling.NoScaling, FourierTransformControl.CreateManaged().Forward, FourierTransformControl.Provider.Forward); + Verify(samples, 10, FourierTransformScaling.NoScaling, ManagedFourierTransformProvider.Instance.Forward, FourierTransformControl.Provider.Forward); } [Test] public void ProviderMatchesManagedProviderArbitraryLarge32_GH286() { var samples = Generate.RandomComplex32(46500, GetUniform(1)); - Verify(samples, 5, FourierTransformScaling.NoScaling, FourierTransformControl.CreateManaged().Forward, FourierTransformControl.Provider.Forward); + Verify(samples, 5, FourierTransformScaling.NoScaling, ManagedFourierTransformProvider.Instance.Forward, FourierTransformControl.Provider.Forward); } [Test] public void ProviderMatchesManagedProviderArbitraryLarge64_GH286() { var samples = Generate.RandomComplex(46500, GetUniform(1)); - Verify(samples, 10, FourierTransformScaling.NoScaling, FourierTransformControl.CreateManaged().Forward, FourierTransformControl.Provider.Forward); + Verify(samples, 10, FourierTransformScaling.NoScaling, ManagedFourierTransformProvider.Instance.Forward, FourierTransformControl.Provider.Forward); } [Test, Explicit("Long-Running")] diff --git a/src/Numerics/Providers/FourierTransform/FourierTransformControl.cs b/src/Numerics/Providers/FourierTransform/FourierTransformControl.cs index 0fd0828b..43c4c265 100644 --- a/src/Numerics/Providers/FourierTransform/FourierTransformControl.cs +++ b/src/Numerics/Providers/FourierTransform/FourierTransformControl.cs @@ -79,8 +79,7 @@ namespace MathNet.Numerics.Providers.FourierTransform } } - public static IFourierTransformProvider CreateManaged() => new ManagedFourierTransformProvider(); - public static void UseManaged() => Provider = CreateManaged(); + public static void UseManaged() => Provider = ManagedFourierTransformProvider.Instance; public static void UseNativeMKL() => Provider = MklProbe.Create(); public static bool TryUseNativeMKL() => TryUse(MklProbe.TryCreate()); diff --git a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Bluestein.cs b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Bluestein.cs index 17479b63..17efd697 100644 --- a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Bluestein.cs +++ b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Bluestein.cs @@ -2,7 +2,7 @@ // Math.NET Numerics, part of the Math.NET Project // https://numerics.mathdotnet.com // -// Copyright (c) 2009-2018 Math.NET +// Copyright (c) 2009-2021 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -32,7 +32,7 @@ using Complex = System.Numerics.Complex; namespace MathNet.Numerics.Providers.FourierTransform { - internal partial class ManagedFourierTransformProvider + public partial class ManagedFourierTransformProvider { /// /// Sequences with length greater than Math.Sqrt(Int32.MaxValue) + 1 @@ -45,7 +45,7 @@ namespace MathNet.Numerics.Providers.FourierTransform /// /// Number of samples. /// Bluestein sequence exp(I*Pi*k^2/N) - private static Complex32[] BluesteinSequence32(int n) + static Complex32[] BluesteinSequence32(int n) { double s = Constants.Pi / n; var sequence = new Complex32[n]; @@ -77,7 +77,7 @@ namespace MathNet.Numerics.Providers.FourierTransform /// /// Number of samples. /// Bluestein sequence exp(I*Pi*k^2/N) - private static Complex[] BluesteinSequence(int n) + static Complex[] BluesteinSequence(int n) { double s = Constants.Pi / n; var sequence = new Complex[n]; @@ -108,7 +108,7 @@ namespace MathNet.Numerics.Providers.FourierTransform /// Convolution with the bluestein sequence (Parallel Version). /// /// Sample Vector. - private static void BluesteinConvolutionParallel(Complex32[] samples) + static void BluesteinConvolutionParallel(Complex32[] samples) { int n = samples.Length; Complex32[] sequence = BluesteinSequence32(n); @@ -163,7 +163,7 @@ namespace MathNet.Numerics.Providers.FourierTransform /// Convolution with the bluestein sequence (Parallel Version). /// /// Sample Vector. - private static void BluesteinConvolutionParallel(Complex[] samples) + static void BluesteinConvolutionParallel(Complex[] samples) { int n = samples.Length; Complex[] sequence = BluesteinSequence(n); @@ -218,7 +218,7 @@ namespace MathNet.Numerics.Providers.FourierTransform /// Swap the real and imaginary parts of each sample. /// /// Sample Vector. - private static void SwapRealImaginary(Complex32[] samples) + static void SwapRealImaginary(Complex32[] samples) { for (int i = 0; i < samples.Length; i++) { @@ -230,7 +230,7 @@ namespace MathNet.Numerics.Providers.FourierTransform /// Swap the real and imaginary parts of each sample. /// /// Sample Vector. - private static void SwapRealImaginary(Complex[] samples) + static void SwapRealImaginary(Complex[] samples) { for (int i = 0; i < samples.Length; i++) { @@ -241,7 +241,7 @@ namespace MathNet.Numerics.Providers.FourierTransform /// /// Bluestein generic FFT for arbitrary sized sample vectors. /// - private static void BluesteinForward(Complex[] samples) + static void BluesteinForward(Complex[] samples) { BluesteinConvolutionParallel(samples); } @@ -249,7 +249,7 @@ namespace MathNet.Numerics.Providers.FourierTransform /// /// Bluestein generic FFT for arbitrary sized sample vectors. /// - private static void BluesteinInverse(Complex[] spectrum) + static void BluesteinInverse(Complex[] spectrum) { SwapRealImaginary(spectrum); BluesteinConvolutionParallel(spectrum); @@ -259,7 +259,7 @@ namespace MathNet.Numerics.Providers.FourierTransform /// /// Bluestein generic FFT for arbitrary sized sample vectors. /// - private static void BluesteinForward(Complex32[] samples) + static void BluesteinForward(Complex32[] samples) { BluesteinConvolutionParallel(samples); } @@ -267,7 +267,7 @@ namespace MathNet.Numerics.Providers.FourierTransform /// /// Bluestein generic FFT for arbitrary sized sample vectors. /// - private static void BluesteinInverse(Complex32[] spectrum) + static void BluesteinInverse(Complex32[] spectrum) { SwapRealImaginary(spectrum); BluesteinConvolutionParallel(spectrum); diff --git a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Radix2.cs b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Radix2.cs index 69ba4844..fe6c453f 100644 --- a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Radix2.cs +++ b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Radix2.cs @@ -2,7 +2,7 @@ // Math.NET Numerics, part of the Math.NET Project // https://numerics.mathdotnet.com // -// Copyright (c) 2009-2018 Math.NET +// Copyright (c) 2009-2021 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -33,14 +33,14 @@ using Complex = System.Numerics.Complex; namespace MathNet.Numerics.Providers.FourierTransform { - internal partial class ManagedFourierTransformProvider + public partial class ManagedFourierTransformProvider { /// /// Radix-2 Reorder Helper Method /// /// Sample type /// Sample vector - private static void Radix2Reorder(T[] samples) + static void Radix2Reorder(T[] samples) { var j = 0; for (var i = 0; i < samples.Length - 1; i++) @@ -73,7 +73,7 @@ namespace MathNet.Numerics.Providers.FourierTransform #if !NET40 [MethodImpl(MethodImplOptions.AggressiveInlining)] #endif - private static void Radix2Step(Complex32[] samples, int exponentSign, int levelSize, int k) + static void Radix2Step(Complex32[] samples, int exponentSign, int levelSize, int k) { // Twiddle Factor var exponent = (exponentSign * k) * Constants.Pi / levelSize; @@ -99,7 +99,7 @@ namespace MathNet.Numerics.Providers.FourierTransform #if !NET40 [MethodImpl(MethodImplOptions.AggressiveInlining)] #endif - private static void Radix2Step(Complex[] samples, int exponentSign, int levelSize, int k) + static void Radix2Step(Complex[] samples, int exponentSign, int levelSize, int k) { // Twiddle Factor var exponent = (exponentSign * k) * Constants.Pi / levelSize; @@ -118,7 +118,7 @@ namespace MathNet.Numerics.Providers.FourierTransform /// /// Radix-2 generic FFT for power-of-two sized sample vectors. /// - private static void Radix2Forward(Complex32[] data) + static void Radix2Forward(Complex32[] data) { Radix2Reorder(data); for (var levelSize = 1; levelSize < data.Length; levelSize *= 2) @@ -133,7 +133,7 @@ namespace MathNet.Numerics.Providers.FourierTransform /// /// Radix-2 generic FFT for power-of-two sized sample vectors. /// - private static void Radix2Forward(Complex[] data) + static void Radix2Forward(Complex[] data) { Radix2Reorder(data); for (var levelSize = 1; levelSize < data.Length; levelSize *= 2) @@ -148,7 +148,7 @@ namespace MathNet.Numerics.Providers.FourierTransform /// /// Radix-2 generic FFT for power-of-two sized sample vectors. /// - private static void Radix2Inverse(Complex32[] data) + static void Radix2Inverse(Complex32[] data) { Radix2Reorder(data); for (var levelSize = 1; levelSize < data.Length; levelSize *= 2) @@ -163,7 +163,7 @@ namespace MathNet.Numerics.Providers.FourierTransform /// /// Radix-2 generic FFT for power-of-two sized sample vectors. /// - private static void Radix2Inverse(Complex[] data) + static void Radix2Inverse(Complex[] data) { Radix2Reorder(data); for (var levelSize = 1; levelSize < data.Length; levelSize *= 2) @@ -178,7 +178,7 @@ namespace MathNet.Numerics.Providers.FourierTransform /// /// Radix-2 generic FFT for power-of-two sample vectors (Parallel Version). /// - private static void Radix2ForwardParallel(Complex32[] data) + static void Radix2ForwardParallel(Complex32[] data) { Radix2Reorder(data); for (var levelSize = 1; levelSize < data.Length; levelSize *= 2) @@ -198,7 +198,7 @@ namespace MathNet.Numerics.Providers.FourierTransform /// /// Radix-2 generic FFT for power-of-two sample vectors (Parallel Version). /// - private static void Radix2ForwardParallel(Complex[] data) + static void Radix2ForwardParallel(Complex[] data) { Radix2Reorder(data); for (var levelSize = 1; levelSize < data.Length; levelSize *= 2) @@ -218,7 +218,7 @@ namespace MathNet.Numerics.Providers.FourierTransform /// /// Radix-2 generic FFT for power-of-two sample vectors (Parallel Version). /// - private static void Radix2InverseParallel(Complex32[] data) + static void Radix2InverseParallel(Complex32[] data) { Radix2Reorder(data); for (var levelSize = 1; levelSize < data.Length; levelSize *= 2) @@ -238,7 +238,7 @@ namespace MathNet.Numerics.Providers.FourierTransform /// /// Radix-2 generic FFT for power-of-two sample vectors (Parallel Version). /// - private static void Radix2InverseParallel(Complex[] data) + static void Radix2InverseParallel(Complex[] data) { Radix2Reorder(data); for (var levelSize = 1; levelSize < data.Length; levelSize *= 2) diff --git a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Scaling.cs b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Scaling.cs new file mode 100644 index 00000000..f9d0673f --- /dev/null +++ b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Scaling.cs @@ -0,0 +1,90 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// https://numerics.mathdotnet.com +// +// Copyright (c) 2009-2021 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 System; +using Complex = System.Numerics.Complex; + +namespace MathNet.Numerics.Providers.FourierTransform +{ + public partial class ManagedFourierTransformProvider + { + + + /// + /// Fully rescale the FFT result. + /// + /// Sample Vector. + static void FullRescale(Complex32[] samples) + { + var scalingFactor = (float)1.0 / samples.Length; + for (int i = 0; i < samples.Length; i++) + { + samples[i] *= scalingFactor; + } + } + + /// + /// Fully rescale the FFT result. + /// + /// Sample Vector. + static void FullRescale(Complex[] samples) + { + var scalingFactor = 1.0 / samples.Length; + for (int i = 0; i < samples.Length; i++) + { + samples[i] *= scalingFactor; + } + } + + /// + /// Half rescale the FFT result (e.g. for symmetric transforms). + /// + /// Sample Vector. + static void HalfRescale(Complex32[] samples) + { + var scalingFactor = (float)Math.Sqrt(1.0 / samples.Length); + for (int i = 0; i < samples.Length; i++) + { + samples[i] *= scalingFactor; + } + } + + /// + /// Fully rescale the FFT result (e.g. for symmetric transforms). + /// + /// Sample Vector. + static void HalfRescale(Complex[] samples) + { + var scalingFactor = Math.Sqrt(1.0 / samples.Length); + for (int i = 0; i < samples.Length; i++) + { + samples[i] *= scalingFactor; + } + } + } +} diff --git a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs index 8949f544..36c6a4ba 100644 --- a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs @@ -2,7 +2,7 @@ // Math.NET Numerics, part of the Math.NET Project // https://numerics.mathdotnet.com // -// Copyright (c) 2009-2018 Math.NET +// Copyright (c) 2009-2021 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -31,8 +31,10 @@ using Complex = System.Numerics.Complex; namespace MathNet.Numerics.Providers.FourierTransform { - internal partial class ManagedFourierTransformProvider : IFourierTransformProvider + public sealed partial class ManagedFourierTransformProvider : IFourierTransformProvider { + public static ManagedFourierTransformProvider Instance { get; } = new ManagedFourierTransformProvider(); + /// /// Try to find out whether the provider is available, at least in principle. /// Verification may still fail if available, but it will certainly fail if unavailable. @@ -53,7 +55,7 @@ namespace MathNet.Numerics.Providers.FourierTransform /// Frees memory buffers, caches and handles allocated in or to the provider. /// Does not unload the provider itself, it is still usable afterwards. /// - public virtual void FreeResources() + public void FreeResources() { } @@ -335,57 +337,5 @@ namespace MathNet.Numerics.Providers.FourierTransform { throw new NotSupportedException(); } - - /// - /// Fully rescale the FFT result. - /// - /// Sample Vector. - private static void FullRescale(Complex32[] samples) - { - var scalingFactor = (float)1.0 / samples.Length; - for (int i = 0; i < samples.Length; i++) - { - samples[i] *= scalingFactor; - } - } - - /// - /// Fully rescale the FFT result. - /// - /// Sample Vector. - private static void FullRescale(Complex[] samples) - { - var scalingFactor = 1.0 / samples.Length; - for (int i = 0; i < samples.Length; i++) - { - samples[i] *= scalingFactor; - } - } - - /// - /// Half rescale the FFT result (e.g. for symmetric transforms). - /// - /// Sample Vector. - private static void HalfRescale(Complex32[] samples) - { - var scalingFactor = (float)Math.Sqrt(1.0 / samples.Length); - for (int i = 0; i < samples.Length; i++) - { - samples[i] *= scalingFactor; - } - } - - /// - /// Fully rescale the FFT result (e.g. for symmetric transforms). - /// - /// Sample Vector. - private static void HalfRescale(Complex[] samples) - { - var scalingFactor = Math.Sqrt(1.0 / samples.Length); - for (int i = 0; i < samples.Length; i++) - { - samples[i] *= scalingFactor; - } - } } } diff --git a/src/Numerics/Providers/LinearAlgebra/LinearAlgebraControl.cs b/src/Numerics/Providers/LinearAlgebra/LinearAlgebraControl.cs index eb131e0b..7ea86678 100644 --- a/src/Numerics/Providers/LinearAlgebra/LinearAlgebraControl.cs +++ b/src/Numerics/Providers/LinearAlgebra/LinearAlgebraControl.cs @@ -86,8 +86,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra } } - public static ILinearAlgebraProvider CreateManaged() => new ManagedLinearAlgebraProvider(); - public static void UseManaged() => Provider = CreateManaged(); + public static void UseManaged() => Provider = ManagedLinearAlgebraProvider.Instance; public static void UseNativeMKL() => Provider = MklProbe.Create(); public static bool TryUseNativeMKL() => TryUse(MklProbe.TryCreate()); diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs index 7bf0466e..25d47e2a 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2020 Math.NET +// Copyright (c) 2009-2021 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -48,7 +48,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The vector to add to . /// The result of the addition. /// This is similar to the AXPY BLAS routine. - public virtual void AddVectorToScaledVector(Complex[] y, Complex alpha, Complex[] x, Complex[] result) + public void AddVectorToScaledVector(Complex[] y, Complex alpha, Complex[] x, Complex[] result) { if (y == null) { @@ -92,7 +92,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The values to scale. /// This result of the scaling. /// This is similar to the SCAL BLAS routine. - public virtual void ScaleArray(Complex alpha, Complex[] x, Complex[] result) + public void ScaleArray(Complex alpha, Complex[] x, Complex[] result) { if (x == null) { @@ -121,7 +121,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// /// The values to conjugate. /// This result of the conjugation. - public virtual void ConjugateArray(Complex[] x, Complex[] result) + public void ConjugateArray(Complex[] x, Complex[] result) { if (x == null) { @@ -141,7 +141,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The vector y. /// The dot product of x and y. /// This is equivalent to the DOT BLAS routine. - public virtual Complex DotProduct(Complex[] x, Complex[] y) + public Complex DotProduct(Complex[] x, Complex[] y) { if (y == null) { @@ -177,7 +177,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public virtual void AddArrays(Complex[] x, Complex[] y, Complex[] result) + public void AddArrays(Complex[] x, Complex[] y, Complex[] result) { if (y == null) { @@ -215,7 +215,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public virtual void SubtractArrays(Complex[] x, Complex[] y, Complex[] result) + public void SubtractArrays(Complex[] x, Complex[] y, Complex[] result) { if (y == null) { @@ -253,7 +253,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public virtual void PointWiseMultiplyArrays(Complex[] x, Complex[] y, Complex[] result) + public void PointWiseMultiplyArrays(Complex[] x, Complex[] y, Complex[] result) { if (y == null) { @@ -291,7 +291,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public virtual void PointWiseDivideArrays(Complex[] x, Complex[] y, Complex[] result) + public void PointWiseDivideArrays(Complex[] x, Complex[] y, Complex[] result) { if (y == null) { @@ -332,7 +332,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public virtual void PointWisePowerArrays(Complex[] x, Complex[] y, Complex[] result) + public void PointWisePowerArrays(Complex[] x, Complex[] y, Complex[] result) { if (y == null) { @@ -373,7 +373,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// /// The requested of the matrix. /// - public virtual double MatrixNorm(Norm norm, int rows, int columns, Complex[] matrix) + public double MatrixNorm(Norm norm, int rows, int columns, Complex[] matrix) { switch (norm) { @@ -444,7 +444,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// Where to store the result of the multiplication. /// This is a simplified version of the BLAS GEMM routine with alpha /// set to 1.0 and beta set to 0.0, and x and y are not transposed. - public virtual void MatrixMultiply(Complex[] x, int rowsX, int columnsX, Complex[] y, int rowsY, int columnsY, Complex[] result) + public void MatrixMultiply(Complex[] x, int rowsX, int columnsX, Complex[] y, int rowsY, int columnsY, Complex[] result) { if (x == null) { @@ -551,7 +551,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The number of columns in the matrix. /// The value to scale the matrix. /// The c matrix. - public virtual void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] a, int rowsA, int columnsA, Complex[] b, int rowsB, int columnsB, Complex beta, Complex[] c) + public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] a, int rowsA, int columnsA, Complex[] b, int rowsB, int columnsB, Complex beta, Complex[] c) { if (a == null) { @@ -679,7 +679,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The order of the square matrix . /// On exit, it contains the pivot indices. The size of the array must be . /// This is equivalent to the GETRF LAPACK routine. - public virtual void LUFactor(Complex[] data, int order, int[] ipiv) + public void LUFactor(Complex[] data, int order, int[] ipiv) { if (data == null) { @@ -777,7 +777,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The N by N matrix to invert. Contains the inverse On exit. /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. - public virtual void LUInverse(Complex[] a, int order) + public void LUInverse(Complex[] a, int order) { if (a == null) { @@ -801,7 +801,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The order of the square matrix . /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. - public virtual void LUInverseFactored(Complex[] a, int order, int[] ipiv) + public void LUInverseFactored(Complex[] a, int order, int[] ipiv) { if (a == null) { @@ -841,7 +841,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The order of the square matrix . /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. - public virtual void LUSolve(int columnsOfB, Complex[] a, int order, Complex[] b) + public void LUSolve(int columnsOfB, Complex[] a, int order, Complex[] b) { if (a == null) { @@ -884,7 +884,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The pivot indices of . /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRS LAPACK routine. - public virtual void LUSolveFactored(int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b) + public void LUSolveFactored(int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b) { if (a == null) { @@ -983,7 +983,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// the Cholesky factorization. /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. - public virtual void CholeskyFactor(Complex[] a, int order) + public void CholeskyFactor(Complex[] a, int order) { if (a == null) { @@ -1070,7 +1070,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// The number of columns in the B matrix. /// This is equivalent to the POTRF add POTRS LAPACK routines. - public virtual void CholeskySolve(Complex[] a, int orderA, Complex[] b, int columnsB) + public void CholeskySolve(Complex[] a, int orderA, Complex[] b, int columnsB) { if (a == null) { @@ -1106,7 +1106,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The B matrix. /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. - public virtual void CholeskySolveFactored(Complex[] a, int orderA, Complex[] b, int columnsB) + public void CholeskySolveFactored(Complex[] a, int orderA, Complex[] b, int columnsB) { if (a == null) { @@ -1187,7 +1187,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// A min(m,n) vector. On exit, contains additional information /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public virtual void QRFactor(Complex[] r, int rowsR, int columnsR, Complex[] q, Complex[] tau) + public void QRFactor(Complex[] r, int rowsR, int columnsR, Complex[] q, Complex[] tau) { if (r == null) { @@ -1249,7 +1249,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// A min(m,n) vector. On exit, contains additional information /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public virtual void ThinQRFactor(Complex[] a, int rowsA, int columnsA, Complex[] r, Complex[] tau) + public void ThinQRFactor(Complex[] a, int rowsA, int columnsA, Complex[] r, Complex[] tau) { if (r == null) { @@ -1435,7 +1435,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// On exit, the solution matrix. /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. - public virtual void QRSolve(Complex[] a, int rows, int columns, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full) + public void QRSolve(Complex[] a, int rows, int columns, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full) { if (a == null) { @@ -1507,7 +1507,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// On exit, the solution matrix. /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. - public virtual void QRSolveFactored(Complex[] q, Complex[] r, int rowsA, int columnsA, Complex[] tau, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full) + public void QRSolveFactored(Complex[] q, Complex[] r, int rowsA, int columnsA, Complex[] tau, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full) { if (r == null) { @@ -1633,7 +1633,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// If is true, on exit VT contains the transposed /// right singular vectors. /// This is equivalent to the GESVD LAPACK routine. - public virtual void SingularValueDecomposition(bool computeVectors, Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt) + public void SingularValueDecomposition(bool computeVectors, Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt) { if (a == null) { @@ -2294,7 +2294,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public virtual void SvdSolve(Complex[] a, int rowsA, int columnsA, Complex[] b, int columnsB, Complex[] x) + public void SvdSolve(Complex[] a, int rowsA, int columnsA, Complex[] b, int columnsB, Complex[] x) { if (a == null) { @@ -2342,7 +2342,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public virtual void SvdSolveFactored(int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, int columnsB, Complex[] x) + public void SvdSolveFactored(int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, int columnsB, Complex[] x) { if (s == null) { @@ -2437,7 +2437,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// On output, the matrix contains the eigen vectors. The length of the array must be order * order. /// On output, the eigen values (λ) of matrix in ascending value. The length of the array must . /// On output, the block diagonal eigenvalue matrix. The length of the array must be order * order. - public virtual void EigenDecomp(bool isSymmetric, int order, Complex[] matrix, Complex[] matrixEv, Complex[] vectorEv, Complex[] matrixD) + public void EigenDecomp(bool isSymmetric, int order, Complex[] matrix, Complex[] matrixEv, Complex[] vectorEv, Complex[] matrixD) { if (matrix == null) { @@ -3229,7 +3229,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// /// Assumes that and have already been transposed. /// - protected static void GetRow(Transpose transpose, int rowindx, int numRows, int numCols, Complex[] matrix, Complex[] row) + static void GetRow(Transpose transpose, int rowindx, int numRows, int numCols, Complex[] matrix, Complex[] row) { if (transpose == Transpose.DontTranspose) { @@ -3255,7 +3255,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// /// Assumes that and have already been transposed. /// - protected static void GetColumn(Transpose transpose, int colindx, int numRows, int numCols, Complex[] matrix, Complex[] column) + static void GetColumn(Transpose transpose, int colindx, int numRows, int numCols, Complex[] matrix, Complex[] column) { if (transpose == Transpose.DontTranspose) { diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs index f069b7e3..468f4bde 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2020 Math.NET +// Copyright (c) 2009-2021 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -48,7 +48,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The vector to add to . /// The result of the addition. /// This is similar to the AXPY BLAS routine. - public virtual void AddVectorToScaledVector(Complex32[] y, Complex32 alpha, Complex32[] x, Complex32[] result) + public void AddVectorToScaledVector(Complex32[] y, Complex32 alpha, Complex32[] x, Complex32[] result) { if (y == null) { @@ -93,7 +93,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The values to scale. /// This result of the scaling. /// This is similar to the SCAL BLAS routine. - public virtual void ScaleArray(Complex32 alpha, Complex32[] x, Complex32[] result) + public void ScaleArray(Complex32 alpha, Complex32[] x, Complex32[] result) { if (x == null) { @@ -122,7 +122,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// /// The values to conjugate. /// This result of the conjugation. - public virtual void ConjugateArray(Complex32[] x, Complex32[] result) + public void ConjugateArray(Complex32[] x, Complex32[] result) { if (x == null) { @@ -142,7 +142,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The vector y. /// The dot product of x and y. /// This is equivalent to the DOT BLAS routine. - public virtual Complex32 DotProduct(Complex32[] x, Complex32[] y) + public Complex32 DotProduct(Complex32[] x, Complex32[] y) { if (y == null) { @@ -178,7 +178,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public virtual void AddArrays(Complex32[] x, Complex32[] y, Complex32[] result) + public void AddArrays(Complex32[] x, Complex32[] y, Complex32[] result) { if (y == null) { @@ -216,7 +216,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public virtual void SubtractArrays(Complex32[] x, Complex32[] y, Complex32[] result) + public void SubtractArrays(Complex32[] x, Complex32[] y, Complex32[] result) { if (y == null) { @@ -254,7 +254,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public virtual void PointWiseMultiplyArrays(Complex32[] x, Complex32[] y, Complex32[] result) + public void PointWiseMultiplyArrays(Complex32[] x, Complex32[] y, Complex32[] result) { if (y == null) { @@ -292,7 +292,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public virtual void PointWiseDivideArrays(Complex32[] x, Complex32[] y, Complex32[] result) + public void PointWiseDivideArrays(Complex32[] x, Complex32[] y, Complex32[] result) { if (y == null) { @@ -333,7 +333,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public virtual void PointWisePowerArrays(Complex32[] x, Complex32[] y, Complex32[] result) + public void PointWisePowerArrays(Complex32[] x, Complex32[] y, Complex32[] result) { if (y == null) { @@ -372,7 +372,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The number of columns. /// The matrix to compute the norm from. /// The requested of the matrix. - public virtual double MatrixNorm(Norm norm, int rows, int columns, Complex32[] matrix) + public double MatrixNorm(Norm norm, int rows, int columns, Complex32[] matrix) { switch (norm) { @@ -444,7 +444,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// Where to store the result of the multiplication. /// This is a simplified version of the BLAS GEMM routine with alpha /// set to 1.0 and beta set to 0.0, and x and y are not transposed. - public virtual void MatrixMultiply(Complex32[] x, int rowsX, int columnsX, Complex32[] y, int rowsY, int columnsY, Complex32[] result) + public void MatrixMultiply(Complex32[] x, int rowsX, int columnsX, Complex32[] y, int rowsY, int columnsY, Complex32[] result) { if (x == null) { @@ -551,7 +551,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The number of columns in the matrix. /// The value to scale the matrix. /// The c matrix. - public virtual void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex32 alpha, Complex32[] a, int rowsA, int columnsA, Complex32[] b, int rowsB, int columnsB, Complex32 beta, Complex32[] c) + public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex32 alpha, Complex32[] a, int rowsA, int columnsA, Complex32[] b, int rowsB, int columnsB, Complex32 beta, Complex32[] c) { if (a == null) { @@ -679,7 +679,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The order of the square matrix . /// On exit, it contains the pivot indices. The size of the array must be . /// This is equivalent to the GETRF LAPACK routine. - public virtual void LUFactor(Complex32[] data, int order, int[] ipiv) + public void LUFactor(Complex32[] data, int order, int[] ipiv) { if (data == null) { @@ -777,7 +777,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The N by N matrix to invert. Contains the inverse On exit. /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. - public virtual void LUInverse(Complex32[] a, int order) + public void LUInverse(Complex32[] a, int order) { if (a == null) { @@ -801,7 +801,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The order of the square matrix . /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. - public virtual void LUInverseFactored(Complex32[] a, int order, int[] ipiv) + public void LUInverseFactored(Complex32[] a, int order, int[] ipiv) { if (a == null) { @@ -841,7 +841,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The order of the square matrix . /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. - public virtual void LUSolve(int columnsOfB, Complex32[] a, int order, Complex32[] b) + public void LUSolve(int columnsOfB, Complex32[] a, int order, Complex32[] b) { if (a == null) { @@ -884,7 +884,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The pivot indices of . /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRS LAPACK routine. - public virtual void LUSolveFactored(int columnsOfB, Complex32[] a, int order, int[] ipiv, Complex32[] b) + public void LUSolveFactored(int columnsOfB, Complex32[] a, int order, int[] ipiv, Complex32[] b) { if (a == null) { @@ -983,7 +983,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// the Cholesky factorization. /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. - public virtual void CholeskyFactor(Complex32[] a, int order) + public void CholeskyFactor(Complex32[] a, int order) { if (a == null) { @@ -1070,7 +1070,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// The number of columns in the B matrix. /// This is equivalent to the POTRF add POTRS LAPACK routines. - public virtual void CholeskySolve(Complex32[] a, int orderA, Complex32[] b, int columnsB) + public void CholeskySolve(Complex32[] a, int orderA, Complex32[] b, int columnsB) { if (a == null) { @@ -1106,7 +1106,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. - public virtual void CholeskySolveFactored(Complex32[] a, int orderA, Complex32[] b, int columnsB) + public void CholeskySolveFactored(Complex32[] a, int orderA, Complex32[] b, int columnsB) { if (a == null) { @@ -1187,7 +1187,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// A min(m,n) vector. On exit, contains additional information /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public virtual void QRFactor(Complex32[] r, int rowsR, int columnsR, Complex32[] q, Complex32[] tau) + public void QRFactor(Complex32[] r, int rowsR, int columnsR, Complex32[] q, Complex32[] tau) { if (r == null) { @@ -1249,7 +1249,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// A min(m,n) vector. On exit, contains additional information /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public virtual void ThinQRFactor(Complex32[] a, int rowsA, int columnsA, Complex32[] r, Complex32[] tau) + public void ThinQRFactor(Complex32[] a, int rowsA, int columnsA, Complex32[] r, Complex32[] tau) { if (r == null) { @@ -1435,7 +1435,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// On exit, the solution matrix. /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. - public virtual void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full) + public void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full) { if (a == null) { @@ -1505,7 +1505,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// On exit, the solution matrix. /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. - public virtual void QRSolveFactored(Complex32[] q, Complex32[] r, int rowsA, int columnsA, Complex32[] tau, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full) + public void QRSolveFactored(Complex32[] q, Complex32[] r, int rowsA, int columnsA, Complex32[] tau, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full) { if (r == null) { @@ -1631,7 +1631,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// If is true, on exit VT contains the transposed /// right singular vectors. /// This is equivalent to the GESVD LAPACK routine. - public virtual void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt) + public void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt) { if (a == null) { @@ -2292,7 +2292,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public virtual void SvdSolve(Complex32[] a, int rowsA, int columnsA, Complex32[] b, int columnsB, Complex32[] x) + public void SvdSolve(Complex32[] a, int rowsA, int columnsA, Complex32[] b, int columnsB, Complex32[] x) { if (a == null) { @@ -2340,7 +2340,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public virtual void SvdSolveFactored(int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, int columnsB, Complex32[] x) + public void SvdSolveFactored(int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, int columnsB, Complex32[] x) { if (s == null) { @@ -2435,7 +2435,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// On output, the matrix contains the eigen vectors. The length of the array must be order * order. /// On output, the eigen values (λ) of matrix in ascending value. The length of the array must . /// On output, the block diagonal eigenvalue matrix. The length of the array must be order * order. - public virtual void EigenDecomp(bool isSymmetric, int order, Complex32[] matrix, Complex32[] matrixEv, Complex[] vectorEv, Complex32[] matrixD) + public void EigenDecomp(bool isSymmetric, int order, Complex32[] matrix, Complex32[] matrixEv, Complex[] vectorEv, Complex32[] matrixD) { if (matrix == null) { @@ -3231,7 +3231,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// /// Assumes that and have already been transposed. /// - protected static void GetRow(Transpose transpose, int rowindx, int numRows, int numCols, Complex32[] matrix, Complex32[] row) + static void GetRow(Transpose transpose, int rowindx, int numRows, int numCols, Complex32[] matrix, Complex32[] row) { if (transpose == Transpose.DontTranspose) { @@ -3257,7 +3257,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// /// Assumes that and have already been transposed. /// - protected static void GetColumn(Transpose transpose, int colindx, int numRows, int numCols, Complex32[] matrix, Complex32[] column) + static void GetColumn(Transpose transpose, int colindx, int numRows, int numCols, Complex32[] matrix, Complex32[] column) { if (transpose == Transpose.DontTranspose) { diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs index d0f96c2e..1ddd3ff9 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2020 Math.NET +// Copyright (c) 2009-2021 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -48,7 +48,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The vector to add to . /// The result of the addition. /// This is similar to the AXPY BLAS routine. - public virtual void AddVectorToScaledVector(double[] y, double alpha, double[] x, double[] result) + public void AddVectorToScaledVector(double[] y, double alpha, double[] x, double[] result) { if (y == null) { @@ -92,7 +92,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The values to scale. /// This result of the scaling. /// This is similar to the SCAL BLAS routine. - public virtual void ScaleArray(double alpha, double[] x, double[] result) + public void ScaleArray(double alpha, double[] x, double[] result) { if (x == null) { @@ -121,7 +121,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// /// The values to conjugate. /// This result of the conjugation. - public virtual void ConjugateArray(double[] x, double[] result) + public void ConjugateArray(double[] x, double[] result) { if (x == null) { @@ -141,7 +141,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The vector y. /// The dot product of x and y. /// This is equivalent to the DOT BLAS routine. - public virtual double DotProduct(double[] x, double[] y) + public double DotProduct(double[] x, double[] y) { if (y == null) { @@ -177,7 +177,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public virtual void AddArrays(double[] x, double[] y, double[] result) + public void AddArrays(double[] x, double[] y, double[] result) { if (y == null) { @@ -215,7 +215,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public virtual void SubtractArrays(double[] x, double[] y, double[] result) + public void SubtractArrays(double[] x, double[] y, double[] result) { if (y == null) { @@ -253,7 +253,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public virtual void PointWiseMultiplyArrays(double[] x, double[] y, double[] result) + public void PointWiseMultiplyArrays(double[] x, double[] y, double[] result) { if (y == null) { @@ -291,7 +291,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public virtual void PointWiseDivideArrays(double[] x, double[] y, double[] result) + public void PointWiseDivideArrays(double[] x, double[] y, double[] result) { if (y == null) { @@ -332,7 +332,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public virtual void PointWisePowerArrays(double[] x, double[] y, double[] result) + public void PointWisePowerArrays(double[] x, double[] y, double[] result) { if (y == null) { @@ -373,7 +373,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// /// The requested of the matrix. /// - public virtual double MatrixNorm(Norm norm, int rows, int columns, double[] matrix) + public double MatrixNorm(Norm norm, int rows, int columns, double[] matrix) { switch (norm) { @@ -444,7 +444,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// Where to store the result of the multiplication. /// This is a simplified version of the BLAS GEMM routine with alpha /// set to 1.0 and beta set to 0.0, and x and y are not transposed. - public virtual void MatrixMultiply(double[] x, int rowsX, int columnsX, double[] y, int rowsY, int columnsY, double[] result) + public void MatrixMultiply(double[] x, int rowsX, int columnsX, double[] y, int rowsY, int columnsY, double[] result) { if (x == null) { @@ -551,7 +551,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The number of columns in the matrix. /// The value to scale the matrix. /// The c matrix. - public virtual void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, int rowsA, int columnsA, double[] b, int rowsB, int columnsB, double beta, double[] c) + public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, int rowsA, int columnsA, double[] b, int rowsB, int columnsB, double beta, double[] c) { if (a == null) { @@ -679,7 +679,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The order of the square matrix . /// On exit, it contains the pivot indices. The size of the array must be . /// This is equivalent to the GETRF LAPACK routine. - public virtual void LUFactor(double[] data, int order, int[] ipiv) + public void LUFactor(double[] data, int order, int[] ipiv) { if (data == null) { @@ -777,7 +777,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The N by N matrix to invert. Contains the inverse On exit. /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. - public virtual void LUInverse(double[] a, int order) + public void LUInverse(double[] a, int order) { if (a == null) { @@ -801,7 +801,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The order of the square matrix . /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. - public virtual void LUInverseFactored(double[] a, int order, int[] ipiv) + public void LUInverseFactored(double[] a, int order, int[] ipiv) { if (a == null) { @@ -841,7 +841,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The order of the square matrix . /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. - public virtual void LUSolve(int columnsOfB, double[] a, int order, double[] b) + public void LUSolve(int columnsOfB, double[] a, int order, double[] b) { if (a == null) { @@ -884,7 +884,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The pivot indices of . /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRS LAPACK routine. - public virtual void LUSolveFactored(int columnsOfB, double[] a, int order, int[] ipiv, double[] b) + public void LUSolveFactored(int columnsOfB, double[] a, int order, int[] ipiv, double[] b) { if (a == null) { @@ -983,7 +983,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// the Cholesky factorization. /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. - public virtual void CholeskyFactor(double[] a, int order) + public void CholeskyFactor(double[] a, int order) { if (a == null) { @@ -1070,7 +1070,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// The number of columns in the B matrix. /// This is equivalent to the POTRF add POTRS LAPACK routines. - public virtual void CholeskySolve(double[] a, int orderA, double[] b, int columnsB) + public void CholeskySolve(double[] a, int orderA, double[] b, int columnsB) { if (a == null) { @@ -1106,7 +1106,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. - public virtual void CholeskySolveFactored(double[] a, int orderA, double[] b, int columnsB) + public void CholeskySolveFactored(double[] a, int orderA, double[] b, int columnsB) { if (a == null) { @@ -1187,7 +1187,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// A min(m,n) vector. On exit, contains additional information /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public virtual void QRFactor(double[] r, int rowsR, int columnsR, double[] q, double[] tau) + public void QRFactor(double[] r, int rowsR, int columnsR, double[] q, double[] tau) { if (r == null) { @@ -1248,7 +1248,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// A min(m,n) vector. On exit, contains additional information /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public virtual void ThinQRFactor(double[] a, int rowsA, int columnsA, double[] r, double[] tau) + public void ThinQRFactor(double[] a, int rowsA, int columnsA, double[] r, double[] tau) { if (r == null) { @@ -1434,7 +1434,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// On exit, the solution matrix. /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. - public virtual void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full) + public void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full) { if (a == null) { @@ -1504,7 +1504,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// On exit, the solution matrix. /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. - public virtual void QRSolveFactored(double[] q, double[] r, int rowsA, int columnsA, double[] tau, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full) + public void QRSolveFactored(double[] q, double[] r, int rowsA, int columnsA, double[] tau, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full) { if (r == null) { @@ -1630,7 +1630,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// If is true, on exit VT contains the transposed /// right singular vectors. /// This is equivalent to the GESVD LAPACK routine. - public virtual void SingularValueDecomposition(bool computeVectors, double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt) + public void SingularValueDecomposition(bool computeVectors, double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt) { if (a == null) { @@ -2352,7 +2352,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public virtual void SvdSolve(double[] a, int rowsA, int columnsA, double[] b, int columnsB, double[] x) + public void SvdSolve(double[] a, int rowsA, int columnsA, double[] b, int columnsB, double[] x) { if (a == null) { @@ -2399,7 +2399,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public virtual void SvdSolveFactored(int rowsA, int columnsA, double[] s, double[] u, double[] vt, double[] b, int columnsB, double[] x) + public void SvdSolveFactored(int rowsA, int columnsA, double[] s, double[] u, double[] vt, double[] b, int columnsB, double[] x) { if (s == null) { @@ -2494,7 +2494,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// On output, the matrix contains the eigen vectors. The length of the array must be order * order. /// On output, the eigen values (λ) of matrix in ascending value. The length of the array must . /// On output, the block diagonal eigenvalue matrix. The length of the array must be order * order. - public virtual void EigenDecomp(bool isSymmetric, int order, double[] matrix, double[] matrixEv, Complex[] vectorEv, double[] matrixD) + public void EigenDecomp(bool isSymmetric, int order, double[] matrix, double[] matrixEv, Complex[] vectorEv, double[] matrixD) { if (matrix == null) { diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs index 27d21528..1acf5e96 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2020 Math.NET +// Copyright (c) 2009-2021 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -48,7 +48,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The vector to add to . /// The result of the addition. /// This is similar to the AXPY BLAS routine. - public virtual void AddVectorToScaledVector(float[] y, float alpha, float[] x, float[] result) + public void AddVectorToScaledVector(float[] y, float alpha, float[] x, float[] result) { if (y == null) { @@ -92,7 +92,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The values to scale. /// This result of the scaling. /// This is similar to the SCAL BLAS routine. - public virtual void ScaleArray(float alpha, float[] x, float[] result) + public void ScaleArray(float alpha, float[] x, float[] result) { if (x == null) { @@ -121,7 +121,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// /// The values to conjugate. /// This result of the conjugation. - public virtual void ConjugateArray(float[] x, float[] result) + public void ConjugateArray(float[] x, float[] result) { if (x == null) { @@ -141,7 +141,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The vector y. /// The dot product of x and y. /// This is equivalent to the DOT BLAS routine. - public virtual float DotProduct(float[] x, float[] y) + public float DotProduct(float[] x, float[] y) { if (y == null) { @@ -177,7 +177,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public virtual void AddArrays(float[] x, float[] y, float[] result) + public void AddArrays(float[] x, float[] y, float[] result) { if (y == null) { @@ -215,7 +215,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public virtual void SubtractArrays(float[] x, float[] y, float[] result) + public void SubtractArrays(float[] x, float[] y, float[] result) { if (y == null) { @@ -253,7 +253,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public virtual void PointWiseMultiplyArrays(float[] x, float[] y, float[] result) + public void PointWiseMultiplyArrays(float[] x, float[] y, float[] result) { if (y == null) { @@ -291,7 +291,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public virtual void PointWiseDivideArrays(float[] x, float[] y, float[] result) + public void PointWiseDivideArrays(float[] x, float[] y, float[] result) { if (y == null) { @@ -332,7 +332,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public virtual void PointWisePowerArrays(float[] x, float[] y, float[] result) + public void PointWisePowerArrays(float[] x, float[] y, float[] result) { if (y == null) { @@ -373,7 +373,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// /// The requested of the matrix. /// - public virtual double MatrixNorm(Norm norm, int rows, int columns, float[] matrix) + public double MatrixNorm(Norm norm, int rows, int columns, float[] matrix) { switch (norm) { @@ -444,7 +444,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// Where to store the result of the multiplication. /// This is a simplified version of the BLAS GEMM routine with alpha /// set to 1.0 and beta set to 0.0, and x and y are not transposed. - public virtual void MatrixMultiply(float[] x, int rowsX, int columnsX, float[] y, int rowsY, int columnsY, float[] result) + public void MatrixMultiply(float[] x, int rowsX, int columnsX, float[] y, int rowsY, int columnsY, float[] result) { if (x == null) { @@ -551,7 +551,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The number of columns in the matrix. /// The value to scale the matrix. /// The c matrix. - public virtual void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, float alpha, float[] a, int rowsA, int columnsA, float[] b, int rowsB, int columnsB, float beta, float[] c) + public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, float alpha, float[] a, int rowsA, int columnsA, float[] b, int rowsB, int columnsB, float beta, float[] c) { if (a == null) { @@ -679,7 +679,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The order of the square matrix . /// On exit, it contains the pivot indices. The size of the array must be . /// This is equivalent to the GETRF LAPACK routine. - public virtual void LUFactor(float[] data, int order, int[] ipiv) + public void LUFactor(float[] data, int order, int[] ipiv) { if (data == null) { @@ -777,7 +777,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The N by N matrix to invert. Contains the inverse On exit. /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. - public virtual void LUInverse(float[] a, int order) + public void LUInverse(float[] a, int order) { if (a == null) { @@ -801,7 +801,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The order of the square matrix . /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. - public virtual void LUInverseFactored(float[] a, int order, int[] ipiv) + public void LUInverseFactored(float[] a, int order, int[] ipiv) { if (a == null) { @@ -841,7 +841,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The order of the square matrix . /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. - public virtual void LUSolve(int columnsOfB, float[] a, int order, float[] b) + public void LUSolve(int columnsOfB, float[] a, int order, float[] b) { if (a == null) { @@ -884,7 +884,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The pivot indices of . /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRS LAPACK routine. - public virtual void LUSolveFactored(int columnsOfB, float[] a, int order, int[] ipiv, float[] b) + public void LUSolveFactored(int columnsOfB, float[] a, int order, int[] ipiv, float[] b) { if (a == null) { @@ -983,7 +983,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// the Cholesky factorization. /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. - public virtual void CholeskyFactor(float[] a, int order) + public void CholeskyFactor(float[] a, int order) { if (a == null) { @@ -1070,7 +1070,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// The number of columns in the B matrix. /// This is equivalent to the POTRF add POTRS LAPACK routines. - public virtual void CholeskySolve(float[] a, int orderA, float[] b, int columnsB) + public void CholeskySolve(float[] a, int orderA, float[] b, int columnsB) { if (a == null) { @@ -1106,7 +1106,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. - public virtual void CholeskySolveFactored(float[] a, int orderA, float[] b, int columnsB) + public void CholeskySolveFactored(float[] a, int orderA, float[] b, int columnsB) { if (a == null) { @@ -1187,7 +1187,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// A min(m,n) vector. On exit, contains additional information /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public virtual void QRFactor(float[] r, int rowsR, int columnsR, float[] q, float[] tau) + public void QRFactor(float[] r, int rowsR, int columnsR, float[] q, float[] tau) { if (r == null) { @@ -1249,7 +1249,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// A min(m,n) vector. On exit, contains additional information /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public virtual void ThinQRFactor(float[] a, int rowsA, int columnsA, float[] r, float[] tau) + public void ThinQRFactor(float[] a, int rowsA, int columnsA, float[] r, float[] tau) { if (r == null) { @@ -1437,7 +1437,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// On exit, the solution matrix. /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. - public virtual void QRSolve(float[] a, int rows, int columns, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full) + public void QRSolve(float[] a, int rows, int columns, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full) { if (a == null) { @@ -1507,7 +1507,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// On exit, the solution matrix. /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. - public virtual void QRSolveFactored(float[] q, float[] r, int rowsA, int columnsA, float[] tau, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full) + public void QRSolveFactored(float[] q, float[] r, int rowsA, int columnsA, float[] tau, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full) { if (r == null) { @@ -1633,7 +1633,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// If is true, on exit VT contains the transposed /// right singular vectors. /// This is equivalent to the GESVD LAPACK routine. - public virtual void SingularValueDecomposition(bool computeVectors, float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt) + public void SingularValueDecomposition(bool computeVectors, float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt) { if (a == null) { @@ -2357,7 +2357,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public virtual void SvdSolve(float[] a, int rowsA, int columnsA, float[] b, int columnsB, float[] x) + public void SvdSolve(float[] a, int rowsA, int columnsA, float[] b, int columnsB, float[] x) { if (a == null) { @@ -2405,7 +2405,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public virtual void SvdSolveFactored(int rowsA, int columnsA, float[] s, float[] u, float[] vt, float[] b, int columnsB, float[] x) + public void SvdSolveFactored(int rowsA, int columnsA, float[] s, float[] u, float[] vt, float[] b, int columnsB, float[] x) { if (s == null) { @@ -2500,7 +2500,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// On output, the matrix contains the eigen vectors. The length of the array must be order * order. /// On output, the eigen values (λ) of matrix in ascending value. The length of the array must . /// On output, the block diagonal eigenvalue matrix. The length of the array must be order * order. - public virtual void EigenDecomp(bool isSymmetric, int order, float[] matrix, float[] matrixEv, Complex[] vectorEv, float[] matrixD) + public void EigenDecomp(bool isSymmetric, int order, float[] matrix, float[] matrixEv, Complex[] vectorEv, float[] matrixD) { if (matrix == null) { diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.cs index 13cc6b71..810dafc4 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2020 Math.NET +// Copyright (c) 2009-2021 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -34,13 +34,15 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// /// The managed linear algebra provider. /// - public partial class ManagedLinearAlgebraProvider : ILinearAlgebraProvider + public sealed partial class ManagedLinearAlgebraProvider : ILinearAlgebraProvider { + public static ManagedLinearAlgebraProvider Instance { get; } = new ManagedLinearAlgebraProvider(); + /// /// Try to find out whether the provider is available, at least in principle. /// Verification may still fail if available, but it will certainly fail if unavailable. /// - public virtual bool IsAvailable() + public bool IsAvailable() { return true; } @@ -48,7 +50,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// /// Initialize and verify that the provided is indeed available. If not, fall back to alternatives like the managed provider /// - public virtual void InitializeVerify() + public void InitializeVerify() { } @@ -56,7 +58,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// Frees memory buffers, caches and handles allocated in or to the provider. /// Does not unload the provider itself, it is still usable afterwards. /// - public virtual void FreeResources() + public void FreeResources() { } diff --git a/src/Numerics/Providers/SparseSolver/ManagedSparseSolverProvider.cs b/src/Numerics/Providers/SparseSolver/ManagedSparseSolverProvider.cs index 9042b0e1..141412aa 100644 --- a/src/Numerics/Providers/SparseSolver/ManagedSparseSolverProvider.cs +++ b/src/Numerics/Providers/SparseSolver/ManagedSparseSolverProvider.cs @@ -6,8 +6,10 @@ namespace MathNet.Numerics.Providers.SparseSolver /// /// The managed sparse solver provider /// - internal class ManagedSparseSolverProvider : ISparseSolverProvider + public sealed class ManagedSparseSolverProvider : ISparseSolverProvider { + public static ManagedSparseSolverProvider Instance { get; } = new ManagedSparseSolverProvider(); + /// /// Try to find out whether the provider is available, at least in principle. /// Verification may still fail if available, but it will certainly fail if unavailable. diff --git a/src/Numerics/Providers/SparseSolver/SparseSolverControl.cs b/src/Numerics/Providers/SparseSolver/SparseSolverControl.cs index 7b83f3ad..e9fef226 100644 --- a/src/Numerics/Providers/SparseSolver/SparseSolverControl.cs +++ b/src/Numerics/Providers/SparseSolver/SparseSolverControl.cs @@ -79,8 +79,7 @@ namespace MathNet.Numerics.Providers.SparseSolver } } - public static ISparseSolverProvider CreateManaged() => new ManagedSparseSolverProvider(); - public static void UseManaged() => Provider = CreateManaged(); + public static void UseManaged() => Provider = ManagedSparseSolverProvider.Instance; public static void UseNativeMKL() => Provider = MklProbe.Create(); public static bool TryUseNativeMKL() => TryUse(MklProbe.TryCreate()); @@ -117,19 +116,6 @@ namespace MathNet.Numerics.Providers.SparseSolver } } - public static bool TryUse(Lazy> providerCreator) - { - try - { - return TryUse(providerCreator.Value?.CreateProvider()); - } - catch - { - // intentionally swallow exceptions here - use the explicit variants if you're interested in why - return false; - } - } - /// /// Use the best provider available. /// diff --git a/src/Providers.CUDA/LinearAlgebra/CudaLinearAlgebraProvider.Complex.cs b/src/Providers.CUDA/LinearAlgebra/CudaLinearAlgebraProvider.Complex.cs index e1cf9db1..8b67a6a1 100644 --- a/src/Providers.CUDA/LinearAlgebra/CudaLinearAlgebraProvider.Complex.cs +++ b/src/Providers.CUDA/LinearAlgebra/CudaLinearAlgebraProvider.Complex.cs @@ -29,6 +29,7 @@ using System; using System.Security; +using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Providers.LinearAlgebra; using Complex = System.Numerics.Complex; @@ -39,6 +40,21 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// internal partial class CudaLinearAlgebraProvider { + /// + /// Computes the requested of the matrix. + /// + /// The type of norm to compute. + /// The number of rows. + /// The number of columns. + /// The matrix to compute the norm from. + /// + /// The requested of the matrix. + /// + public double MatrixNorm(Norm norm, int rows, int columns, Complex[] matrix) + { + return ManagedLinearAlgebraProvider.Instance.MatrixNorm(norm, rows, columns, matrix); + } + /// /// Computes the dot product of x and y. /// @@ -47,7 +63,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The dot product of x and y. /// This is equivalent to the DOT BLAS routine. [SecuritySafeCritical] - public override Complex DotProduct(Complex[] x, Complex[] y) + public Complex DotProduct(Complex[] x, Complex[] y) { if (y == null) { @@ -67,6 +83,196 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra return SafeNativeMethods.z_dot_product(_blasHandle, x.Length, x, y); } + /// + /// Does a point wise add of two arrays z = x + y. This can be used + /// to add vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the addition. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void AddArrays(Complex[] x, Complex[] y, Complex[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] + y[i]; + } + } + + /// + /// Does a point wise subtraction of two arrays z = x - y. This can be used + /// to subtract vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the subtraction. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void SubtractArrays(Complex[] x, Complex[] y, Complex[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] - y[i]; + } + } + + /// + /// Does a point wise multiplication of two arrays z = x * y. This can be used + /// to multiple elements of vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the point wise multiplication. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void PointWiseMultiplyArrays(Complex[] x, Complex[] y, Complex[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] * y[i]; + } + } + + /// + /// Does a point wise division of two arrays z = x / y. This can be used + /// to divide elements of vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the point wise division. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void PointWiseDivideArrays(Complex[] x, Complex[] y, Complex[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] / y[i]; + } + } + + /// + /// Does a point wise power of two arrays z = x ^ y. This can be used + /// to raise elements of vectors or matrices to the powers of another vector or matrix. + /// + /// The array x. + /// The array y. + /// The result of the point wise power. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void PointWisePowerArrays(Complex[] x, Complex[] y, Complex[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = Complex.Pow(x[i], y[i]); + } + } + /// /// Adds a scaled vector to another: result = y + alpha*x. /// @@ -76,7 +282,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The result of the addition. /// This is similar to the AXPY BLAS routine. [SecuritySafeCritical] - public override void AddVectorToScaledVector(Complex[] y, Complex alpha, Complex[] x, Complex[] result) + public void AddVectorToScaledVector(Complex[] y, Complex alpha, Complex[] x, Complex[] result) { if (y == null) { @@ -114,7 +320,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// This result of the scaling. /// This is similar to the SCAL BLAS routine. [SecuritySafeCritical] - public override void ScaleArray(Complex alpha, Complex[] x, Complex[] result) + public void ScaleArray(Complex alpha, Complex[] x, Complex[] result) { if (x == null) { @@ -134,6 +340,24 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra SafeNativeMethods.z_scale(_blasHandle, x.Length, alpha, result); } + /// + /// Conjugates an array. Can be used to conjugate a vector and a matrix. + /// + /// The values to conjugate. + /// This result of the conjugation. + public void ConjugateArray(Complex[] x, Complex[] result) + { + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i].Conjugate(); + } + } + /// /// Multiples two matrices. result = x * y /// @@ -146,7 +370,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// Where to store the result of the multiplication. /// This is a simplified version of the BLAS GEMM routine with alpha /// set to Complex.One and beta set to Complex.Zero, and x and y are not transposed. - public override void MatrixMultiply(Complex[] x, int rowsX, int columnsX, Complex[] y, int rowsY, int columnsY, Complex[] result) + public void MatrixMultiply(Complex[] x, int rowsX, int columnsX, Complex[] y, int rowsY, int columnsY, Complex[] result) { MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex.One, x, rowsX, columnsX, y, rowsY, columnsY, Complex.Zero, result); } @@ -166,7 +390,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The value to scale the matrix. /// The c matrix. [SecuritySafeCritical] - public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] a, int rowsA, int columnsA, Complex[] b, int rowsB, int columnsB, Complex beta, Complex[] c) + public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] a, int rowsA, int columnsA, Complex[] b, int rowsB, int columnsB, Complex beta, Complex[] c) { if (a == null) { @@ -211,7 +435,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// On exit, it contains the pivot indices. The size of the array must be . /// This is equivalent to the GETRF LAPACK routine. [SecuritySafeCritical] - public override void LUFactor(Complex[] data, int order, int[] ipiv) + public void LUFactor(Complex[] data, int order, int[] ipiv) { if (data == null) { @@ -243,7 +467,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. [SecuritySafeCritical] - public override void LUInverse(Complex[] a, int order) + public void LUInverse(Complex[] a, int order) { if (a == null) { @@ -266,7 +490,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. [SecuritySafeCritical] - public override void LUInverseFactored(Complex[] a, int order, int[] ipiv) + public void LUInverseFactored(Complex[] a, int order, int[] ipiv) { if (a == null) { @@ -300,7 +524,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. [SecuritySafeCritical] - public override void LUSolve(int columnsOfB, Complex[] a, int order, Complex[] b) + public void LUSolve(int columnsOfB, Complex[] a, int order, Complex[] b) { if (a == null) { @@ -335,7 +559,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRS LAPACK routine. [SecuritySafeCritical] - public override void LUSolveFactored(int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b) + public void LUSolveFactored(int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b) { if (a == null) { @@ -378,7 +602,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. [SecuritySafeCritical] - public override void CholeskyFactor(Complex[] a, int order) + public void CholeskyFactor(Complex[] a, int order) { if (a == null) { @@ -408,7 +632,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// This is equivalent to the POTRF add POTRS LAPACK routines. /// [SecuritySafeCritical] - public override void CholeskySolve(Complex[] a, int orderA, Complex[] b, int columnsB) + public void CholeskySolve(Complex[] a, int orderA, Complex[] b, int columnsB) { if (a == null) { @@ -442,7 +666,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. [SecuritySafeCritical] - public override void CholeskySolveFactored(Complex[] a, int orderA, Complex[] b, int columnsB) + public void CholeskySolveFactored(Complex[] a, int orderA, Complex[] b, int columnsB) { if (a == null) { @@ -467,6 +691,75 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra Solver(SafeNativeMethods.z_cholesky_solve_factored(_solverHandle, orderA, columnsB, a, b)); } + /// + /// Computes the QR factorization of A. + /// + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the R matrix of the QR factorization. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// On exit, A M by M matrix that holds the Q matrix of the + /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + public void QRFactor(Complex[] r, int rowsR, int columnsR, Complex[] q, Complex[] tau) + { + ManagedLinearAlgebraProvider.Instance.QRFactor(r, rowsR, columnsR, q, tau); + } + + /// + /// Computes the thin QR factorization of A where M > N. + /// + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the Q matrix of the QR factorization. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// On exit, A N by N matrix that holds the R matrix of the + /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + public void ThinQRFactor(Complex[] q, int rowsA, int columnsA, Complex[] r, Complex[] tau) + { + ManagedLinearAlgebraProvider.Instance.ThinQRFactor(q, rowsA, columnsA, r, tau); + } + + /// + /// Solves A*X=B for X using QR factorization of A. + /// + /// The A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + /// The type of QR factorization to perform. + /// Rows must be greater or equal to columns. + public void QRSolve(Complex[] a, int rows, int columns, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full) + { + ManagedLinearAlgebraProvider.Instance.QRSolve(a, rows, columns, b, columnsB, x, method); + } + + /// + /// Solves A*X=B for X using a previously QR factored matrix. + /// + /// The Q matrix obtained by calling . + /// The R matrix obtained by calling . + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// Contains additional information on Q. Only used for the native solver + /// and can be null for the managed provider. + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + /// The type of QR factorization to perform. + /// Rows must be greater or equal to columns. + public void QRSolveFactored(Complex[] q, Complex[] r, int rowsA, int columnsA, Complex[] tau, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full) + { + ManagedLinearAlgebraProvider.Instance.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, method); + } + /// /// Solves A*X=B for X using the singular value decomposition of A. /// @@ -476,7 +769,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public override void SvdSolve(Complex[] a, int rowsA, int columnsA, Complex[] b, int columnsB, Complex[] x) + public void SvdSolve(Complex[] a, int rowsA, int columnsA, Complex[] b, int columnsB, Complex[] x) { if (a == null) { @@ -513,6 +806,22 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); } + /// + /// Solves A*X=B for X using a previously SVD decomposed matrix. + /// + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The s values returned by . + /// The left singular vectors returned by . + /// The right singular vectors returned by . + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + public void SvdSolveFactored(int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, int columnsB, Complex[] x) + { + ManagedLinearAlgebraProvider.Instance.SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); + } + /// /// Computes the singular value decomposition of A. /// @@ -527,7 +836,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// right singular vectors. /// This is equivalent to the GESVD LAPACK routine. [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt) + public void SingularValueDecomposition(bool computeVectors, Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt) { if (a == null) { @@ -565,8 +874,22 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra } if (columnsA > rowsA || !computeVectors) // see remarks http://docs.nvidia.com/cuda/cusolver/index.html#cuds-lt-t-gt-gesvd - base.SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt); + ManagedLinearAlgebraProvider.Instance.SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt); else Solver(SafeNativeMethods.z_svd_factor(_solverHandle, computeVectors, rowsA, columnsA, a, s, u, vt)); } + + /// + /// Computes the eigenvalues and eigenvectors of a matrix. + /// + /// Whether the matrix is symmetric or not. + /// The order of the matrix. + /// The matrix to decompose. The length of the array must be order * order. + /// On output, the matrix contains the eigen vectors. The length of the array must be order * order. + /// On output, the eigen values (λ) of matrix in ascending value. The length of the array must . + /// On output, the block diagonal eigenvalue matrix. The length of the array must be order * order. + public void EigenDecomp(bool isSymmetric, int order, Complex[] matrix, Complex[] matrixEv, Complex[] vectorEv, Complex[] matrixD) + { + ManagedLinearAlgebraProvider.Instance.EigenDecomp(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD); + } } } diff --git a/src/Providers.CUDA/LinearAlgebra/CudaLinearAlgebraProvider.Complex32.cs b/src/Providers.CUDA/LinearAlgebra/CudaLinearAlgebraProvider.Complex32.cs index 3137a3b7..478caf02 100644 --- a/src/Providers.CUDA/LinearAlgebra/CudaLinearAlgebraProvider.Complex32.cs +++ b/src/Providers.CUDA/LinearAlgebra/CudaLinearAlgebraProvider.Complex32.cs @@ -29,7 +29,9 @@ using System; using System.Security; +using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Providers.LinearAlgebra; +using Complex = System.Numerics.Complex; namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra { @@ -38,6 +40,21 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// internal partial class CudaLinearAlgebraProvider { + /// + /// Computes the requested of the matrix. + /// + /// The type of norm to compute. + /// The number of rows. + /// The number of columns. + /// The matrix to compute the norm from. + /// + /// The requested of the matrix. + /// + public double MatrixNorm(Norm norm, int rows, int columns, Complex32[] matrix) + { + return ManagedLinearAlgebraProvider.Instance.MatrixNorm(norm, rows, columns, matrix); + } + /// /// Computes the dot product of x and y. /// @@ -46,7 +63,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The dot product of x and y. /// This is equivalent to the DOT BLAS routine. [SecuritySafeCritical] - public override Complex32 DotProduct(Complex32[] x, Complex32[] y) + public Complex32 DotProduct(Complex32[] x, Complex32[] y) { if (y == null) { @@ -66,6 +83,196 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra return SafeNativeMethods.c_dot_product(_blasHandle, x.Length, x, y); } + /// + /// Does a point wise add of two arrays z = x + y. This can be used + /// to add vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the addition. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void AddArrays(Complex32[] x, Complex32[] y, Complex32[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] + y[i]; + } + } + + /// + /// Does a point wise subtraction of two arrays z = x - y. This can be used + /// to subtract vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the subtraction. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void SubtractArrays(Complex32[] x, Complex32[] y, Complex32[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] - y[i]; + } + } + + /// + /// Does a point wise multiplication of two arrays z = x * y. This can be used + /// to multiple elements of vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the point wise multiplication. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void PointWiseMultiplyArrays(Complex32[] x, Complex32[] y, Complex32[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] * y[i]; + } + } + + /// + /// Does a point wise division of two arrays z = x / y. This can be used + /// to divide elements of vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the point wise division. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void PointWiseDivideArrays(Complex32[] x, Complex32[] y, Complex32[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < y.Length; i++) + { + result[i] = x[i] / y[i]; + } + } + + /// + /// Does a point wise power of two arrays z = x ^ y. This can be used + /// to raise elements of vectors or matrices to the powers of another vector or matrix. + /// + /// The array x. + /// The array y. + /// The result of the point wise power. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void PointWisePowerArrays(Complex32[] x, Complex32[] y, Complex32[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < y.Length; i++) + { + result[i] = Complex32.Pow(x[i], y[i]); + } + } + /// /// Adds a scaled vector to another: result = y + alpha*x. /// @@ -75,7 +282,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The result of the addition. /// This is similar to the AXPY BLAS routine. [SecuritySafeCritical] - public override void AddVectorToScaledVector(Complex32[] y, Complex32 alpha, Complex32[] x, Complex32[] result) + public void AddVectorToScaledVector(Complex32[] y, Complex32 alpha, Complex32[] x, Complex32[] result) { if (y == null) { @@ -113,7 +320,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// This result of the scaling. /// This is similar to the SCAL BLAS routine. [SecuritySafeCritical] - public override void ScaleArray(Complex32 alpha, Complex32[] x, Complex32[] result) + public void ScaleArray(Complex32 alpha, Complex32[] x, Complex32[] result) { if (x == null) { @@ -133,6 +340,24 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra SafeNativeMethods.c_scale(_blasHandle, x.Length, alpha, result); } + /// + /// Conjugates an array. Can be used to conjugate a vector and a matrix. + /// + /// The values to conjugate. + /// This result of the conjugation. + public void ConjugateArray(Complex32[] x, Complex32[] result) + { + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i].Conjugate(); + } + } + /// /// Multiples two matrices. result = x * y /// @@ -145,7 +370,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// Where to store the result of the multiplication. /// This is a simplified version of the BLAS GEMM routine with alpha /// set to Complex32.One and beta set to Complex32.Zero, and x and y are not transposed. - public override void MatrixMultiply(Complex32[] x, int rowsX, int columnsX, Complex32[] y, int rowsY, int columnsY, Complex32[] result) + public void MatrixMultiply(Complex32[] x, int rowsX, int columnsX, Complex32[] y, int rowsY, int columnsY, Complex32[] result) { MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex32.One, x, rowsX, columnsX, y, rowsY, columnsY, Complex32.Zero, result); } @@ -165,7 +390,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The value to scale the matrix. /// The c matrix. [SecuritySafeCritical] - public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex32 alpha, Complex32[] a, int rowsA, int columnsA, Complex32[] b, int rowsB, int columnsB, Complex32 beta, Complex32[] c) + public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex32 alpha, Complex32[] a, int rowsA, int columnsA, Complex32[] b, int rowsB, int columnsB, Complex32 beta, Complex32[] c) { if (a == null) { @@ -210,7 +435,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// On exit, it contains the pivot indices. The size of the array must be . /// This is equivalent to the GETRF LAPACK routine. [SecuritySafeCritical] - public override void LUFactor(Complex32[] data, int order, int[] ipiv) + public void LUFactor(Complex32[] data, int order, int[] ipiv) { if (data == null) { @@ -242,7 +467,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. [SecuritySafeCritical] - public override void LUInverse(Complex32[] a, int order) + public void LUInverse(Complex32[] a, int order) { if (a == null) { @@ -265,7 +490,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. [SecuritySafeCritical] - public override void LUInverseFactored(Complex32[] a, int order, int[] ipiv) + public void LUInverseFactored(Complex32[] a, int order, int[] ipiv) { if (a == null) { @@ -299,7 +524,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. [SecuritySafeCritical] - public override void LUSolve(int columnsOfB, Complex32[] a, int order, Complex32[] b) + public void LUSolve(int columnsOfB, Complex32[] a, int order, Complex32[] b) { if (a == null) { @@ -334,7 +559,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRS LAPACK routine. [SecuritySafeCritical] - public override void LUSolveFactored(int columnsOfB, Complex32[] a, int order, int[] ipiv, Complex32[] b) + public void LUSolveFactored(int columnsOfB, Complex32[] a, int order, int[] ipiv, Complex32[] b) { if (a == null) { @@ -377,7 +602,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. [SecuritySafeCritical] - public override void CholeskyFactor(Complex32[] a, int order) + public void CholeskyFactor(Complex32[] a, int order) { if (a == null) { @@ -407,7 +632,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// This is equivalent to the POTRF add POTRS LAPACK routines. /// [SecuritySafeCritical] - public override void CholeskySolve(Complex32[] a, int orderA, Complex32[] b, int columnsB) + public void CholeskySolve(Complex32[] a, int orderA, Complex32[] b, int columnsB) { if (a == null) { @@ -441,7 +666,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. [SecuritySafeCritical] - public override void CholeskySolveFactored(Complex32[] a, int orderA, Complex32[] b, int columnsB) + public void CholeskySolveFactored(Complex32[] a, int orderA, Complex32[] b, int columnsB) { if (a == null) { @@ -466,6 +691,75 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra Solver(SafeNativeMethods.c_cholesky_solve_factored(_solverHandle, orderA, columnsB, a, b)); } + /// + /// Computes the QR factorization of A. + /// + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the R matrix of the QR factorization. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// On exit, A M by M matrix that holds the Q matrix of the + /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + public void QRFactor(Complex32[] r, int rowsR, int columnsR, Complex32[] q, Complex32[] tau) + { + ManagedLinearAlgebraProvider.Instance.QRFactor(r, rowsR, columnsR, q, tau); + } + + /// + /// Computes the thin QR factorization of A where M > N. + /// + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the Q matrix of the QR factorization. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// On exit, A N by N matrix that holds the R matrix of the + /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + public void ThinQRFactor(Complex32[] q, int rowsA, int columnsA, Complex32[] r, Complex32[] tau) + { + ManagedLinearAlgebraProvider.Instance.ThinQRFactor(q, rowsA, columnsA, r, tau); + } + + /// + /// Solves A*X=B for X using QR factorization of A. + /// + /// The A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + /// The type of QR factorization to perform. + /// Rows must be greater or equal to columns. + public void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full) + { + ManagedLinearAlgebraProvider.Instance.QRSolve(a, rows, columns, b, columnsB, x, method); + } + + /// + /// Solves A*X=B for X using a previously QR factored matrix. + /// + /// The Q matrix obtained by calling . + /// The R matrix obtained by calling . + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// Contains additional information on Q. Only used for the native solver + /// and can be null for the managed provider. + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + /// The type of QR factorization to perform. + /// Rows must be greater or equal to columns. + public void QRSolveFactored(Complex32[] q, Complex32[] r, int rowsA, int columnsA, Complex32[] tau, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full) + { + ManagedLinearAlgebraProvider.Instance.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, method); + } + /// /// Solves A*X=B for X using the singular value decomposition of A. /// @@ -475,7 +769,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public override void SvdSolve(Complex32[] a, int rowsA, int columnsA, Complex32[] b, int columnsB, Complex32[] x) + public void SvdSolve(Complex32[] a, int rowsA, int columnsA, Complex32[] b, int columnsB, Complex32[] x) { if (a == null) { @@ -512,6 +806,22 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); } + /// + /// Solves A*X=B for X using a previously SVD decomposed matrix. + /// + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The s values returned by . + /// The left singular vectors returned by . + /// The right singular vectors returned by . + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + public void SvdSolveFactored(int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, int columnsB, Complex32[] x) + { + ManagedLinearAlgebraProvider.Instance.SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); + } + /// /// Computes the singular value decomposition of A. /// @@ -526,7 +836,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// right singular vectors. /// This is equivalent to the GESVD LAPACK routine. [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt) + public void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt) { if (a == null) { @@ -564,8 +874,22 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra } if (columnsA > rowsA || !computeVectors) // see remarks http://docs.nvidia.com/cuda/cusolver/index.html#cuds-lt-t-gt-gesvd - base.SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt); + ManagedLinearAlgebraProvider.Instance.SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt); else Solver(SafeNativeMethods.c_svd_factor(_solverHandle, computeVectors, rowsA, columnsA, a, s, u, vt)); } + + /// + /// Computes the eigenvalues and eigenvectors of a matrix. + /// + /// Whether the matrix is symmetric or not. + /// The order of the matrix. + /// The matrix to decompose. The length of the array must be order * order. + /// On output, the matrix contains the eigen vectors. The length of the array must be order * order. + /// On output, the eigen values (λ) of matrix in ascending value. The length of the array must . + /// On output, the block diagonal eigenvalue matrix. The length of the array must be order * order. + public void EigenDecomp(bool isSymmetric, int order, Complex32[] matrix, Complex32[] matrixEv, Complex[] vectorEv, Complex32[] matrixD) + { + ManagedLinearAlgebraProvider.Instance.EigenDecomp(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD); + } } } diff --git a/src/Providers.CUDA/LinearAlgebra/CudaLinearAlgebraProvider.Double.cs b/src/Providers.CUDA/LinearAlgebra/CudaLinearAlgebraProvider.Double.cs index f9b4ebed..44bdeafe 100644 --- a/src/Providers.CUDA/LinearAlgebra/CudaLinearAlgebraProvider.Double.cs +++ b/src/Providers.CUDA/LinearAlgebra/CudaLinearAlgebraProvider.Double.cs @@ -29,7 +29,9 @@ using System; using System.Security; +using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Providers.LinearAlgebra; +using Complex = System.Numerics.Complex; namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra { @@ -38,6 +40,21 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// internal partial class CudaLinearAlgebraProvider { + /// + /// Computes the requested of the matrix. + /// + /// The type of norm to compute. + /// The number of rows. + /// The number of columns. + /// The matrix to compute the norm from. + /// + /// The requested of the matrix. + /// + public double MatrixNorm(Norm norm, int rows, int columns, double[] matrix) + { + return ManagedLinearAlgebraProvider.Instance.MatrixNorm(norm, rows, columns, matrix); + } + /// /// Computes the dot product of x and y. /// @@ -46,7 +63,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The dot product of x and y. /// This is equivalent to the DOT BLAS routine. [SecuritySafeCritical] - public override double DotProduct(double[] x, double[] y) + public double DotProduct(double[] x, double[] y) { if (y == null) { @@ -66,6 +83,196 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra return SafeNativeMethods.d_dot_product(_blasHandle, x.Length, x, y); } + /// + /// Does a point wise add of two arrays z = x + y. This can be used + /// to add vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the addition. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void AddArrays(double[] x, double[] y, double[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] + y[i]; + } + } + + /// + /// Does a point wise subtraction of two arrays z = x - y. This can be used + /// to subtract vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the subtraction. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void SubtractArrays(double[] x, double[] y, double[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] - y[i]; + } + } + + /// + /// Does a point wise multiplication of two arrays z = x * y. This can be used + /// to multiple elements of vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the point wise multiplication. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void PointWiseMultiplyArrays(double[] x, double[] y, double[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] * y[i]; + } + } + + /// + /// Does a point wise division of two arrays z = x / y. This can be used + /// to divide elements of vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the point wise division. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void PointWiseDivideArrays(double[] x, double[] y, double[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] / y[i]; + } + } + + /// + /// Does a point wise power of two arrays z = x ^ y. This can be used + /// to raise elements of vectors or matrices to the powers of another vector or matrix. + /// + /// The array x. + /// The array y. + /// The result of the point wise power. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void PointWisePowerArrays(double[] x, double[] y, double[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = Math.Pow(x[i], y[i]); + } + } + /// /// Adds a scaled vector to another: result = y + alpha*x. /// @@ -75,7 +282,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The result of the addition. /// This is similar to the AXPY BLAS routine. [SecuritySafeCritical] - public override void AddVectorToScaledVector(double[] y, double alpha, double[] x, double[] result) + public void AddVectorToScaledVector(double[] y, double alpha, double[] x, double[] result) { if (y == null) { @@ -113,7 +320,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// This result of the scaling. /// This is similar to the SCAL BLAS routine. [SecuritySafeCritical] - public override void ScaleArray(double alpha, double[] x, double[] result) + public void ScaleArray(double alpha, double[] x, double[] result) { if (x == null) { @@ -133,6 +340,24 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra SafeNativeMethods.d_scale(_blasHandle, x.Length, alpha, result); } + /// + /// Conjugates an array. Can be used to conjugate a vector and a matrix. + /// + /// The values to conjugate. + /// This result of the conjugation. + public void ConjugateArray(double[] x, double[] result) + { + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (!ReferenceEquals(x, result)) + { + x.CopyTo(result, 0); + } + } + /// /// Multiples two matrices. result = x * y /// @@ -145,7 +370,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// Where to store the result of the multiplication. /// This is a simplified version of the BLAS GEMM routine with alpha /// set to 1.0 and beta set to 0.0, and x and y are not transposed. - public override void MatrixMultiply(double[] x, int rowsX, int columnsX, double[] y, int rowsY, int columnsY, double[] result) + public void MatrixMultiply(double[] x, int rowsX, int columnsX, double[] y, int rowsY, int columnsY, double[] result) { MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, 1.0, x, rowsX, columnsX, y, rowsY, columnsY, 0.0, result); } @@ -165,7 +390,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The value to scale the matrix. /// The c matrix. [SecuritySafeCritical] - public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, int rowsA, int columnsA, double[] b, int rowsB, int columnsB, double beta, double[] c) + public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, int rowsA, int columnsA, double[] b, int rowsB, int columnsB, double beta, double[] c) { if (a == null) { @@ -210,7 +435,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// On exit, it contains the pivot indices. The size of the array must be . /// This is equivalent to the GETRF LAPACK routine. [SecuritySafeCritical] - public override void LUFactor(double[] data, int order, int[] ipiv) + public void LUFactor(double[] data, int order, int[] ipiv) { if (data == null) { @@ -242,7 +467,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. [SecuritySafeCritical] - public override void LUInverse(double[] a, int order) + public void LUInverse(double[] a, int order) { if (a == null) { @@ -265,7 +490,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. [SecuritySafeCritical] - public override void LUInverseFactored(double[] a, int order, int[] ipiv) + public void LUInverseFactored(double[] a, int order, int[] ipiv) { if (a == null) { @@ -299,7 +524,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. [SecuritySafeCritical] - public override void LUSolve(int columnsOfB, double[] a, int order, double[] b) + public void LUSolve(int columnsOfB, double[] a, int order, double[] b) { if (a == null) { @@ -334,7 +559,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRS LAPACK routine. [SecuritySafeCritical] - public override void LUSolveFactored(int columnsOfB, double[] a, int order, int[] ipiv, double[] b) + public void LUSolveFactored(int columnsOfB, double[] a, int order, int[] ipiv, double[] b) { if (a == null) { @@ -377,7 +602,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. [SecuritySafeCritical] - public override void CholeskyFactor(double[] a, int order) + public void CholeskyFactor(double[] a, int order) { if (a == null) { @@ -407,7 +632,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// This is equivalent to the POTRF add POTRS LAPACK routines. /// [SecuritySafeCritical] - public override void CholeskySolve(double[] a, int orderA, double[] b, int columnsB) + public void CholeskySolve(double[] a, int orderA, double[] b, int columnsB) { if (a == null) { @@ -441,7 +666,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. [SecuritySafeCritical] - public override void CholeskySolveFactored(double[] a, int orderA, double[] b, int columnsB) + public void CholeskySolveFactored(double[] a, int orderA, double[] b, int columnsB) { if (a == null) { @@ -466,6 +691,75 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra Solver(SafeNativeMethods.d_cholesky_solve_factored(_solverHandle, orderA, columnsB, a, b)); } + /// + /// Computes the QR factorization of A. + /// + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the R matrix of the QR factorization. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// On exit, A M by M matrix that holds the Q matrix of the + /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + public void QRFactor(double[] r, int rowsR, int columnsR, double[] q, double[] tau) + { + ManagedLinearAlgebraProvider.Instance.QRFactor(r, rowsR, columnsR, q, tau); + } + + /// + /// Computes the thin QR factorization of A where M > N. + /// + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the Q matrix of the QR factorization. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// On exit, A N by N matrix that holds the R matrix of the + /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + public void ThinQRFactor(double[] q, int rowsA, int columnsA, double[] r, double[] tau) + { + ManagedLinearAlgebraProvider.Instance.ThinQRFactor(q, rowsA, columnsA, r, tau); + } + + /// + /// Solves A*X=B for X using QR factorization of A. + /// + /// The A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + /// The type of QR factorization to perform. + /// Rows must be greater or equal to columns. + public void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full) + { + ManagedLinearAlgebraProvider.Instance.QRSolve(a, rows, columns, b, columnsB, x, method); + } + + /// + /// Solves A*X=B for X using a previously QR factored matrix. + /// + /// The Q matrix obtained by calling . + /// The R matrix obtained by calling . + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// Contains additional information on Q. Only used for the native solver + /// and can be null for the managed provider. + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + /// The type of QR factorization to perform. + /// Rows must be greater or equal to columns. + public void QRSolveFactored(double[] q, double[] r, int rowsA, int columnsA, double[] tau, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full) + { + ManagedLinearAlgebraProvider.Instance.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, method); + } + /// /// Solves A*X=B for X using the singular value decomposition of A. /// @@ -475,7 +769,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public override void SvdSolve(double[] a, int rowsA, int columnsA, double[] b, int columnsB, double[] x) + public void SvdSolve(double[] a, int rowsA, int columnsA, double[] b, int columnsB, double[] x) { if (a == null) { @@ -512,6 +806,22 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); } + /// + /// Solves A*X=B for X using a previously SVD decomposed matrix. + /// + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The s values returned by . + /// The left singular vectors returned by . + /// The right singular vectors returned by . + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + public void SvdSolveFactored(int rowsA, int columnsA, double[] s, double[] u, double[] vt, double[] b, int columnsB, double[] x) + { + ManagedLinearAlgebraProvider.Instance.SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); + } + /// /// Computes the singular value decomposition of A. /// @@ -526,7 +836,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// right singular vectors. /// This is equivalent to the GESVD LAPACK routine. [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt) + public void SingularValueDecomposition(bool computeVectors, double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt) { if (a == null) { @@ -564,8 +874,22 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra } if (columnsA > rowsA || !computeVectors) // see remarks http://docs.nvidia.com/cuda/cusolver/index.html#cuds-lt-t-gt-gesvd - base.SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt); + ManagedLinearAlgebraProvider.Instance.SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt); else Solver (SafeNativeMethods.d_svd_factor(_solverHandle, computeVectors, rowsA, columnsA, a, s, u, vt)); } + + /// + /// Computes the eigenvalues and eigenvectors of a matrix. + /// + /// Whether the matrix is symmetric or not. + /// The order of the matrix. + /// The matrix to decompose. The length of the array must be order * order. + /// On output, the matrix contains the eigen vectors. The length of the array must be order * order. + /// On output, the eigen values (λ) of matrix in ascending value. The length of the array must . + /// On output, the block diagonal eigenvalue matrix. The length of the array must be order * order. + public void EigenDecomp(bool isSymmetric, int order, double[] matrix, double[] matrixEv, Complex[] vectorEv, double[] matrixD) + { + ManagedLinearAlgebraProvider.Instance.EigenDecomp(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD); + } } } diff --git a/src/Providers.CUDA/LinearAlgebra/CudaLinearAlgebraProvider.Single.cs b/src/Providers.CUDA/LinearAlgebra/CudaLinearAlgebraProvider.Single.cs index 38847a4c..d70c784a 100644 --- a/src/Providers.CUDA/LinearAlgebra/CudaLinearAlgebraProvider.Single.cs +++ b/src/Providers.CUDA/LinearAlgebra/CudaLinearAlgebraProvider.Single.cs @@ -29,7 +29,9 @@ using System; using System.Security; +using MathNet.Numerics.LinearAlgebra.Factorization; using MathNet.Numerics.Providers.LinearAlgebra; +using Complex = System.Numerics.Complex; namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra { @@ -38,6 +40,21 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// internal partial class CudaLinearAlgebraProvider { + /// + /// Computes the requested of the matrix. + /// + /// The type of norm to compute. + /// The number of rows. + /// The number of columns. + /// The matrix to compute the norm from. + /// + /// The requested of the matrix. + /// + public double MatrixNorm(Norm norm, int rows, int columns, float[] matrix) + { + return ManagedLinearAlgebraProvider.Instance.MatrixNorm(norm, rows, columns, matrix); + } + /// /// Computes the dot product of x and y. /// @@ -46,7 +63,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The dot product of x and y. /// This is equivalent to the DOT BLAS routine. [SecuritySafeCritical] - public override float DotProduct(float[] x, float[] y) + public float DotProduct(float[] x, float[] y) { if (y == null) { @@ -66,6 +83,196 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra return SafeNativeMethods.s_dot_product(_blasHandle, x.Length, x, y); } + /// + /// Does a point wise add of two arrays z = x + y. This can be used + /// to add vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the addition. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void AddArrays(float[] x, float[] y, float[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] + y[i]; + } + } + + /// + /// Does a point wise subtraction of two arrays z = x - y. This can be used + /// to subtract vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the subtraction. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void SubtractArrays(float[] x, float[] y, float[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] - y[i]; + } + } + + /// + /// Does a point wise multiplication of two arrays z = x * y. This can be used + /// to multiple elements of vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the point wise multiplication. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void PointWiseMultiplyArrays(float[] x, float[] y, float[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] * y[i]; + } + } + + /// + /// Does a point wise division of two arrays z = x / y. This can be used + /// to divide elements of vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the point wise division. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void PointWiseDivideArrays(float[] x, float[] y, float[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] / y[i]; + } + } + + /// + /// Does a point wise power of two arrays z = x ^ y. This can be used + /// to raise elements of vectors or matrices to the powers of another vector or matrix. + /// + /// The array x. + /// The array y. + /// The result of the point wise power. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void PointWisePowerArrays(float[] x, float[] y, float[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = (float)Math.Pow(x[i], y[i]); + } + } + /// /// Adds a scaled vector to another: result = y + alpha*x. /// @@ -75,7 +282,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The result of the addition. /// This is similar to the AXPY BLAS routine. [SecuritySafeCritical] - public override void AddVectorToScaledVector(float[] y, float alpha, float[] x, float[] result) + public void AddVectorToScaledVector(float[] y, float alpha, float[] x, float[] result) { if (y == null) { @@ -113,7 +320,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// This result of the scaling. /// This is similar to the SCAL BLAS routine. [SecuritySafeCritical] - public override void ScaleArray(float alpha, float[] x, float[] result) + public void ScaleArray(float alpha, float[] x, float[] result) { if (x == null) { @@ -133,6 +340,24 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra SafeNativeMethods.s_scale(_blasHandle, x.Length, alpha, result); } + /// + /// Conjugates an array. Can be used to conjugate a vector and a matrix. + /// + /// The values to conjugate. + /// This result of the conjugation. + public void ConjugateArray(float[] x, float[] result) + { + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (!ReferenceEquals(x, result)) + { + x.CopyTo(result, 0); + } + } + /// /// Multiples two matrices. result = x * y /// @@ -145,7 +370,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// Where to store the result of the multiplication. /// This is a simplified version of the BLAS GEMM routine with alpha /// set to 1.0f and beta set to 0.0f, and x and y are not transposed. - public override void MatrixMultiply(float[] x, int rowsX, int columnsX, float[] y, int rowsY, int columnsY, float[] result) + public void MatrixMultiply(float[] x, int rowsX, int columnsX, float[] y, int rowsY, int columnsY, float[] result) { MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, 1.0f, x, rowsX, columnsX, y, rowsY, columnsY, 0.0f, result); } @@ -165,7 +390,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The value to scale the matrix. /// The c matrix. [SecuritySafeCritical] - public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, float alpha, float[] a, int rowsA, int columnsA, float[] b, int rowsB, int columnsB, float beta, float[] c) + public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, float alpha, float[] a, int rowsA, int columnsA, float[] b, int rowsB, int columnsB, float beta, float[] c) { if (a == null) { @@ -210,7 +435,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// On exit, it contains the pivot indices. The size of the array must be . /// This is equivalent to the GETRF LAPACK routine. [SecuritySafeCritical] - public override void LUFactor(float[] data, int order, int[] ipiv) + public void LUFactor(float[] data, int order, int[] ipiv) { if (data == null) { @@ -242,7 +467,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. [SecuritySafeCritical] - public override void LUInverse(float[] a, int order) + public void LUInverse(float[] a, int order) { if (a == null) { @@ -265,7 +490,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. [SecuritySafeCritical] - public override void LUInverseFactored(float[] a, int order, int[] ipiv) + public void LUInverseFactored(float[] a, int order, int[] ipiv) { if (a == null) { @@ -299,7 +524,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. [SecuritySafeCritical] - public override void LUSolve(int columnsOfB, float[] a, int order, float[] b) + public void LUSolve(int columnsOfB, float[] a, int order, float[] b) { if (a == null) { @@ -334,7 +559,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRS LAPACK routine. [SecuritySafeCritical] - public override void LUSolveFactored(int columnsOfB, float[] a, int order, int[] ipiv, float[] b) + public void LUSolveFactored(int columnsOfB, float[] a, int order, int[] ipiv, float[] b) { if (a == null) { @@ -377,7 +602,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. [SecuritySafeCritical] - public override void CholeskyFactor(float[] a, int order) + public void CholeskyFactor(float[] a, int order) { if (a == null) { @@ -407,7 +632,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// This is equivalent to the POTRF add POTRS LAPACK routines. /// [SecuritySafeCritical] - public override void CholeskySolve(float[] a, int orderA, float[] b, int columnsB) + public void CholeskySolve(float[] a, int orderA, float[] b, int columnsB) { if (a == null) { @@ -441,7 +666,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. [SecuritySafeCritical] - public override void CholeskySolveFactored(float[] a, int orderA, float[] b, int columnsB) + public void CholeskySolveFactored(float[] a, int orderA, float[] b, int columnsB) { if (a == null) { @@ -466,6 +691,75 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra Solver(SafeNativeMethods.s_cholesky_solve_factored(_solverHandle, orderA, columnsB, a, b)); } + /// + /// Computes the QR factorization of A. + /// + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the R matrix of the QR factorization. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// On exit, A M by M matrix that holds the Q matrix of the + /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + public void QRFactor(float[] r, int rowsR, int columnsR, float[] q, float[] tau) + { + ManagedLinearAlgebraProvider.Instance.QRFactor(r, rowsR, columnsR, q, tau); + } + + /// + /// Computes the thin QR factorization of A where M > N. + /// + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the Q matrix of the QR factorization. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// On exit, A N by N matrix that holds the R matrix of the + /// QR factorization. + /// A min(m,n) vector. On exit, contains additional information + /// to be used by the QR solve routine. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + public void ThinQRFactor(float[] q, int rowsA, int columnsA, float[] r, float[] tau) + { + ManagedLinearAlgebraProvider.Instance.ThinQRFactor(q, rowsA, columnsA, r, tau); + } + + /// + /// Solves A*X=B for X using QR factorization of A. + /// + /// The A matrix. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + /// The type of QR factorization to perform. + /// Rows must be greater or equal to columns. + public void QRSolve(float[] a, int rows, int columns, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full) + { + ManagedLinearAlgebraProvider.Instance.QRSolve(a, rows, columns, b, columnsB, x, method); + } + + /// + /// Solves A*X=B for X using a previously QR factored matrix. + /// + /// The Q matrix obtained by calling . + /// The R matrix obtained by calling . + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// Contains additional information on Q. Only used for the native solver + /// and can be null for the managed provider. + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + /// The type of QR factorization to perform. + /// Rows must be greater or equal to columns. + public void QRSolveFactored(float[] q, float[] r, int rowsA, int columnsA, float[] tau, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full) + { + ManagedLinearAlgebraProvider.Instance.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, method); + } + /// /// Solves A*X=B for X using the singular value decomposition of A. /// @@ -475,7 +769,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public override void SvdSolve(float[] a, int rowsA, int columnsA, float[] b, int columnsB, float[] x) + public void SvdSolve(float[] a, int rowsA, int columnsA, float[] b, int columnsB, float[] x) { if (a == null) { @@ -512,6 +806,22 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); } + /// + /// Solves A*X=B for X using a previously SVD decomposed matrix. + /// + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The s values returned by . + /// The left singular vectors returned by . + /// The right singular vectors returned by . + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + public void SvdSolveFactored(int rowsA, int columnsA, float[] s, float[] u, float[] vt, float[] b, int columnsB, float[] x) + { + ManagedLinearAlgebraProvider.Instance.SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); + } + /// /// Computes the singular value decomposition of A. /// @@ -526,7 +836,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// right singular vectors. /// This is equivalent to the GESVD LAPACK routine. [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt) + public void SingularValueDecomposition(bool computeVectors, float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt) { if (a == null) { @@ -564,8 +874,22 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra } if (columnsA > rowsA || !computeVectors) // see remarks http://docs.nvidia.com/cuda/cusolver/index.html#cuds-lt-t-gt-gesvd - base.SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt); + ManagedLinearAlgebraProvider.Instance.SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt); else Solver(SafeNativeMethods.s_svd_factor(_solverHandle, computeVectors, rowsA, columnsA, a, s, u, vt)); } + + /// + /// Computes the eigenvalues and eigenvectors of a matrix. + /// + /// Whether the matrix is symmetric or not. + /// The order of the matrix. + /// The matrix to decompose. The length of the array must be order * order. + /// On output, the matrix contains the eigen vectors. The length of the array must be order * order. + /// On output, the eigen values (λ) of matrix in ascending value. The length of the array must . + /// On output, the block diagonal eigenvalue matrix. The length of the array must be order * order. + public void EigenDecomp(bool isSymmetric, int order, float[] matrix, float[] matrixEv, Complex[] vectorEv, float[] matrixD) + { + ManagedLinearAlgebraProvider.Instance.EigenDecomp(isSymmetric, order, matrix, matrixEv, vectorEv, matrixD); + } } } diff --git a/src/Providers.CUDA/LinearAlgebra/CudaLinearAlgebraProvider.cs b/src/Providers.CUDA/LinearAlgebra/CudaLinearAlgebraProvider.cs index edef95a9..c1fa5e09 100644 --- a/src/Providers.CUDA/LinearAlgebra/CudaLinearAlgebraProvider.cs +++ b/src/Providers.CUDA/LinearAlgebra/CudaLinearAlgebraProvider.cs @@ -35,7 +35,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// /// NVidia's CUDA Toolkit linear algebra provider. /// - internal partial class CudaLinearAlgebraProvider : ManagedLinearAlgebraProvider, IDisposable + internal sealed partial class CudaLinearAlgebraProvider : ILinearAlgebraProvider, IDisposable { const int MinimumCompatibleRevision = 1; @@ -53,7 +53,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// Try to find out whether the provider is available, at least in principle. /// Verification may still fail if available, but it will certainly fail if unavailable. /// - public override bool IsAvailable() + public bool IsAvailable() { return CudaProvider.IsAvailable(hintPath: _hintPath); } @@ -62,7 +62,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// Initialize and verify that the provided is indeed available. /// If calling this method fails, consider to fall back to alternatives like the managed provider. /// - public override void InitializeVerify() + public void InitializeVerify() { int revision = CudaProvider.Load(hintPath: _hintPath); if (revision < MinimumCompatibleRevision) @@ -86,7 +86,7 @@ namespace MathNet.Numerics.Providers.CUDA.LinearAlgebra /// Frees memory buffers, caches and handles allocated in or to the provider. /// Does not unload the provider itself, it is still usable afterwards. /// - public override void FreeResources() + public void FreeResources() { CudaProvider.FreeResources(); } diff --git a/src/Providers.MKL/FourierTransform/MklFourierTransformProvider.cs b/src/Providers.MKL/FourierTransform/MklFourierTransformProvider.cs index a59b28c9..df44b829 100644 --- a/src/Providers.MKL/FourierTransform/MklFourierTransformProvider.cs +++ b/src/Providers.MKL/FourierTransform/MklFourierTransformProvider.cs @@ -34,7 +34,7 @@ using Complex = System.Numerics.Complex; namespace MathNet.Numerics.Providers.MKL.FourierTransform { - internal class MklFourierTransformProvider : IFourierTransformProvider, IDisposable + internal sealed class MklFourierTransformProvider : IFourierTransformProvider, IDisposable { const int MinimumCompatibleRevision = 11; @@ -91,7 +91,7 @@ namespace MathNet.Numerics.Providers.MKL.FourierTransform /// Does not unload the provider itself, it is still usable afterwards. /// [SecuritySafeCritical] - public virtual void FreeResources() + public void FreeResources() { Kernel kernel = Interlocked.Exchange(ref _kernel, null); if (kernel != null) diff --git a/src/Providers.MKL/LinearAlgebra/MklLinearAlgebraProvider.Complex.cs b/src/Providers.MKL/LinearAlgebra/MklLinearAlgebraProvider.Complex.cs index d14c5474..d4a80289 100644 --- a/src/Providers.MKL/LinearAlgebra/MklLinearAlgebraProvider.Complex.cs +++ b/src/Providers.MKL/LinearAlgebra/MklLinearAlgebraProvider.Complex.cs @@ -51,7 +51,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The requested of the matrix. /// [SecuritySafeCritical] - public override double MatrixNorm(Norm norm, int rows, int columns, Complex[] matrix) + public double MatrixNorm(Norm norm, int rows, int columns, Complex[] matrix) { if (matrix == null) { @@ -84,7 +84,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The dot product of x and y. /// This is equivalent to the DOT BLAS routine. [SecuritySafeCritical] - public override Complex DotProduct(Complex[] x, Complex[] y) + public Complex DotProduct(Complex[] x, Complex[] y) { if (y == null) { @@ -113,7 +113,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The result of the addition. /// This is similar to the AXPY BLAS routine. [SecuritySafeCritical] - public override void AddVectorToScaledVector(Complex[] y, Complex alpha, Complex[] x, Complex[] result) + public void AddVectorToScaledVector(Complex[] y, Complex alpha, Complex[] x, Complex[] result) { if (y == null) { @@ -151,7 +151,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// This result of the scaling. /// This is similar to the SCAL BLAS routine. [SecuritySafeCritical] - public override void ScaleArray(Complex alpha, Complex[] x, Complex[] result) + public void ScaleArray(Complex alpha, Complex[] x, Complex[] result) { if (x == null) { @@ -171,6 +171,24 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra SafeNativeMethods.z_scale(x.Length, alpha, result); } + /// + /// Conjugates an array. Can be used to conjugate a vector and a matrix. + /// + /// The values to conjugate. + /// This result of the conjugation. + public void ConjugateArray(Complex[] x, Complex[] result) + { + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i].Conjugate(); + } + } + /// /// Multiples two matrices. result = x * y /// @@ -183,7 +201,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// Where to store the result of the multiplication. /// This is a simplified version of the BLAS GEMM routine with alpha /// set to Complex.One and beta set to Complex.Zero, and x and y are not transposed. - public override void MatrixMultiply(Complex[] x, int rowsX, int columnsX, Complex[] y, int rowsY, int columnsY, Complex[] result) + public void MatrixMultiply(Complex[] x, int rowsX, int columnsX, Complex[] y, int rowsY, int columnsY, Complex[] result) { MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex.One, x, rowsX, columnsX, y, rowsY, columnsY, Complex.Zero, result); } @@ -203,7 +221,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The value to scale the matrix. /// The c matrix. [SecuritySafeCritical] - public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] a, int rowsA, int columnsA, Complex[] b, int rowsB, int columnsB, Complex beta, Complex[] c) + public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] a, int rowsA, int columnsA, Complex[] b, int rowsB, int columnsB, Complex beta, Complex[] c) { if (a == null) { @@ -248,7 +266,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// On exit, it contains the pivot indices. The size of the array must be . /// This is equivalent to the GETRF LAPACK routine. [SecuritySafeCritical] - public override void LUFactor(Complex[] data, int order, int[] ipiv) + public void LUFactor(Complex[] data, int order, int[] ipiv) { if (data == null) { @@ -285,7 +303,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. [SecuritySafeCritical] - public override void LUInverse(Complex[] a, int order) + public void LUInverse(Complex[] a, int order) { if (a == null) { @@ -323,7 +341,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. [SecuritySafeCritical] - public override void LUInverseFactored(Complex[] a, int order, int[] ipiv) + public void LUInverseFactored(Complex[] a, int order, int[] ipiv) { if (a == null) { @@ -367,7 +385,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. [SecuritySafeCritical] - public override void LUSolve(int columnsOfB, Complex[] a, int order, Complex[] b) + public void LUSolve(int columnsOfB, Complex[] a, int order, Complex[] b) { if (a == null) { @@ -412,7 +430,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRS LAPACK routine. [SecuritySafeCritical] - public override void LUSolveFactored(int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b) + public void LUSolveFactored(int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b) { if (a == null) { @@ -465,7 +483,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. [SecuritySafeCritical] - public override void CholeskyFactor(Complex[] a, int order) + public void CholeskyFactor(Complex[] a, int order) { if (a == null) { @@ -505,7 +523,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// This is equivalent to the POTRF add POTRS LAPACK routines. /// [SecuritySafeCritical] - public override void CholeskySolve(Complex[] a, int orderA, Complex[] b, int columnsB) + public void CholeskySolve(Complex[] a, int orderA, Complex[] b, int columnsB) { if (a == null) { @@ -549,7 +567,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. [SecuritySafeCritical] - public override void CholeskySolveFactored(Complex[] a, int orderA, Complex[] b, int columnsB) + public void CholeskySolveFactored(Complex[] a, int orderA, Complex[] b, int columnsB) { if (a == null) { @@ -592,7 +610,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. [SecuritySafeCritical] - public override void QRFactor(Complex[] r, int rowsR, int columnsR, Complex[] q, Complex[] tau) + public void QRFactor(Complex[] r, int rowsR, int columnsR, Complex[] q, Complex[] tau) { if (r == null) { @@ -640,7 +658,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. [SecuritySafeCritical] - public override void ThinQRFactor(Complex[] q, int rowsA, int columnsA, Complex[] r, Complex[] tau) + public void ThinQRFactor(Complex[] q, int rowsA, int columnsA, Complex[] r, Complex[] tau) { if (r == null) { @@ -687,7 +705,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] - public override void QRSolve(Complex[] a, int rows, int columns, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full) + public void QRSolve(Complex[] a, int rows, int columns, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full) { if (a == null) { @@ -757,7 +775,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] - public override void QRSolveFactored(Complex[] q, Complex[] r, int rowsA, int columnsA, Complex[] tau, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full) + public void QRSolveFactored(Complex[] q, Complex[] r, int rowsA, int columnsA, Complex[] tau, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full) { if (r == null) { @@ -829,7 +847,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra { // we don't have access to the raw Q matrix any more(it is stored in R in the full QR), need to think about this. // let just call the managed version in the meantime. The heavy lifting has already been done. -marcus - base.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin); + ManagedLinearAlgebraProvider.Instance.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin); } } @@ -842,7 +860,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public override void SvdSolve(Complex[] a, int rowsA, int columnsA, Complex[] b, int columnsB, Complex[] x) + public void SvdSolve(Complex[] a, int rowsA, int columnsA, Complex[] b, int columnsB, Complex[] x) { if (a == null) { @@ -879,6 +897,22 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); } + /// + /// Solves A*X=B for X using a previously SVD decomposed matrix. + /// + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The s values returned by . + /// The left singular vectors returned by . + /// The right singular vectors returned by . + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + public void SvdSolveFactored(int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, int columnsB, Complex[] x) + { + ManagedLinearAlgebraProvider.Instance.SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); + } + /// /// Computes the singular value decomposition of A. /// @@ -893,7 +927,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// right singular vectors. /// This is equivalent to the GESVD LAPACK routine. [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt) + public void SingularValueDecomposition(bool computeVectors, Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt) { if (a == null) { @@ -958,7 +992,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public override void AddArrays(Complex[] x, Complex[] y, Complex[] result) + public void AddArrays(Complex[] x, Complex[] y, Complex[] result) { if (y == null) { @@ -993,7 +1027,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public override void SubtractArrays(Complex[] x, Complex[] y, Complex[] result) + public void SubtractArrays(Complex[] x, Complex[] y, Complex[] result) { if (y == null) { @@ -1028,7 +1062,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public override void PointWiseMultiplyArrays(Complex[] x, Complex[] y, Complex[] result) + public void PointWiseMultiplyArrays(Complex[] x, Complex[] y, Complex[] result) { if (y == null) { @@ -1063,13 +1097,8 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public override void PointWisePowerArrays(Complex[] x, Complex[] y, Complex[] result) + public void PointWisePowerArrays(Complex[] x, Complex[] y, Complex[] result) { - if (_vectorFunctionsMajor != 0 || _vectorFunctionsMinor < 1) - { - base.PointWisePowerArrays(x, y, result); - } - if (y == null) { throw new ArgumentNullException(nameof(y)); @@ -1090,6 +1119,16 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra throw new ArgumentException("The array arguments must have the same length."); } + if (_vectorFunctionsMajor != 0 || _vectorFunctionsMinor < 1) + { + for (int i = 0; i < y.Length; i++) + { + result[i] = Complex.Pow(x[i], y[i]); + } + + return; + } + SafeNativeMethods.z_vector_power(x.Length, x, y, result); } @@ -1103,7 +1142,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public override void PointWiseDivideArrays(Complex[] x, Complex[] y, Complex[] result) + public void PointWiseDivideArrays(Complex[] x, Complex[] y, Complex[] result) { if (y == null) { @@ -1137,7 +1176,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// On output, the matrix contains the eigen vectors. The length of the array must be order * order. /// On output, the eigen values (λ) of matrix in ascending value. The length of the array must . /// On output, the block diagonal eigenvalue matrix. The length of the array must be order * order. - public override void EigenDecomp(bool isSymmetric, int order, Complex[] matrix, Complex[] matrixEv, Complex[] vectorEv, Complex[] matrixD) + public void EigenDecomp(bool isSymmetric, int order, Complex[] matrix, Complex[] matrixEv, Complex[] vectorEv, Complex[] matrixD) { if (matrix == null) { diff --git a/src/Providers.MKL/LinearAlgebra/MklLinearAlgebraProvider.Complex32.cs b/src/Providers.MKL/LinearAlgebra/MklLinearAlgebraProvider.Complex32.cs index 4d3c8b3a..f0ba2a88 100644 --- a/src/Providers.MKL/LinearAlgebra/MklLinearAlgebraProvider.Complex32.cs +++ b/src/Providers.MKL/LinearAlgebra/MklLinearAlgebraProvider.Complex32.cs @@ -51,7 +51,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The requested of the matrix. /// [SecuritySafeCritical] - public override double MatrixNorm(Norm norm, int rows, int columns, Complex32[] matrix) + public double MatrixNorm(Norm norm, int rows, int columns, Complex32[] matrix) { if (matrix == null) { @@ -84,7 +84,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The dot product of x and y. /// This is equivalent to the DOT BLAS routine. [SecuritySafeCritical] - public override Complex32 DotProduct(Complex32[] x, Complex32[] y) + public Complex32 DotProduct(Complex32[] x, Complex32[] y) { if (y == null) { @@ -113,7 +113,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The result of the addition. /// This is similar to the AXPY BLAS routine. [SecuritySafeCritical] - public override void AddVectorToScaledVector(Complex32[] y, Complex32 alpha, Complex32[] x, Complex32[] result) + public void AddVectorToScaledVector(Complex32[] y, Complex32 alpha, Complex32[] x, Complex32[] result) { if (y == null) { @@ -151,7 +151,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// This result of the scaling. /// This is similar to the SCAL BLAS routine. [SecuritySafeCritical] - public override void ScaleArray(Complex32 alpha, Complex32[] x, Complex32[] result) + public void ScaleArray(Complex32 alpha, Complex32[] x, Complex32[] result) { if (x == null) { @@ -171,6 +171,24 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra SafeNativeMethods.c_scale(x.Length, alpha, result); } + /// + /// Conjugates an array. Can be used to conjugate a vector and a matrix. + /// + /// The values to conjugate. + /// This result of the conjugation. + public void ConjugateArray(Complex32[] x, Complex32[] result) + { + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i].Conjugate(); + } + } + /// /// Multiples two matrices. result = x * y /// @@ -183,7 +201,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// Where to store the result of the multiplication. /// This is a simplified version of the BLAS GEMM routine with alpha /// set to Complex32.One and beta set to Complex32.Zero, and x and y are not transposed. - public override void MatrixMultiply(Complex32[] x, int rowsX, int columnsX, Complex32[] y, int rowsY, int columnsY, Complex32[] result) + public void MatrixMultiply(Complex32[] x, int rowsX, int columnsX, Complex32[] y, int rowsY, int columnsY, Complex32[] result) { MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex32.One, x, rowsX, columnsX, y, rowsY, columnsY, Complex32.Zero, result); } @@ -203,7 +221,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The value to scale the matrix. /// The c matrix. [SecuritySafeCritical] - public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex32 alpha, Complex32[] a, int rowsA, int columnsA, Complex32[] b, int rowsB, int columnsB, Complex32 beta, Complex32[] c) + public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex32 alpha, Complex32[] a, int rowsA, int columnsA, Complex32[] b, int rowsB, int columnsB, Complex32 beta, Complex32[] c) { if (a == null) { @@ -248,7 +266,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// On exit, it contains the pivot indices. The size of the array must be . /// This is equivalent to the GETRF LAPACK routine. [SecuritySafeCritical] - public override void LUFactor(Complex32[] data, int order, int[] ipiv) + public void LUFactor(Complex32[] data, int order, int[] ipiv) { if (data == null) { @@ -285,7 +303,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. [SecuritySafeCritical] - public override void LUInverse(Complex32[] a, int order) + public void LUInverse(Complex32[] a, int order) { if (a == null) { @@ -323,7 +341,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. [SecuritySafeCritical] - public override void LUInverseFactored(Complex32[] a, int order, int[] ipiv) + public void LUInverseFactored(Complex32[] a, int order, int[] ipiv) { if (a == null) { @@ -367,7 +385,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. [SecuritySafeCritical] - public override void LUSolve(int columnsOfB, Complex32[] a, int order, Complex32[] b) + public void LUSolve(int columnsOfB, Complex32[] a, int order, Complex32[] b) { if (a == null) { @@ -412,7 +430,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRS LAPACK routine. [SecuritySafeCritical] - public override void LUSolveFactored(int columnsOfB, Complex32[] a, int order, int[] ipiv, Complex32[] b) + public void LUSolveFactored(int columnsOfB, Complex32[] a, int order, int[] ipiv, Complex32[] b) { if (a == null) { @@ -460,7 +478,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. [SecuritySafeCritical] - public override void CholeskyFactor(Complex32[] a, int order) + public void CholeskyFactor(Complex32[] a, int order) { if (a == null) { @@ -500,7 +518,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// This is equivalent to the POTRF add POTRS LAPACK routines. /// [SecuritySafeCritical] - public override void CholeskySolve(Complex32[] a, int orderA, Complex32[] b, int columnsB) + public void CholeskySolve(Complex32[] a, int orderA, Complex32[] b, int columnsB) { if (a == null) { @@ -544,7 +562,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. [SecuritySafeCritical] - public override void CholeskySolveFactored(Complex32[] a, int orderA, Complex32[] b, int columnsB) + public void CholeskySolveFactored(Complex32[] a, int orderA, Complex32[] b, int columnsB) { if (a == null) { @@ -587,7 +605,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. [SecuritySafeCritical] - public override void QRFactor(Complex32[] r, int rowsR, int columnsR, Complex32[] q, Complex32[] tau) + public void QRFactor(Complex32[] r, int rowsR, int columnsR, Complex32[] q, Complex32[] tau) { if (r == null) { @@ -635,7 +653,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. [SecuritySafeCritical] - public override void ThinQRFactor(Complex32[] q, int rowsA, int columnsA, Complex32[] r, Complex32[] tau) + public void ThinQRFactor(Complex32[] q, int rowsA, int columnsA, Complex32[] r, Complex32[] tau) { if (r == null) { @@ -682,7 +700,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] - public override void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full) + public void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full) { if (a == null) { @@ -752,7 +770,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] - public override void QRSolveFactored(Complex32[] q, Complex32[] r, int rowsA, int columnsA, Complex32[] tau, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full) + public void QRSolveFactored(Complex32[] q, Complex32[] r, int rowsA, int columnsA, Complex32[] tau, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full) { if (r == null) { @@ -824,7 +842,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra { // we don't have access to the raw Q matrix any more(it is stored in R in the full QR), need to think about this. // let just call the managed version in the meantime. The heavy lifting has already been done. -marcus - base.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin); + ManagedLinearAlgebraProvider.Instance.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin); } } @@ -837,7 +855,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public override void SvdSolve(Complex32[] a, int rowsA, int columnsA, Complex32[] b, int columnsB, Complex32[] x) + public void SvdSolve(Complex32[] a, int rowsA, int columnsA, Complex32[] b, int columnsB, Complex32[] x) { if (a == null) { @@ -874,6 +892,22 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); } + /// + /// Solves A*X=B for X using a previously SVD decomposed matrix. + /// + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The s values returned by . + /// The left singular vectors returned by . + /// The right singular vectors returned by . + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + public void SvdSolveFactored(int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, int columnsB, Complex32[] x) + { + ManagedLinearAlgebraProvider.Instance.SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); + } + /// /// Computes the singular value decomposition of A. /// @@ -888,7 +922,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// right singular vectors. /// This is equivalent to the GESVD LAPACK routine. [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt) + public void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt) { if (a == null) { @@ -953,7 +987,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public override void AddArrays(Complex32[] x, Complex32[] y, Complex32[] result) + public void AddArrays(Complex32[] x, Complex32[] y, Complex32[] result) { if (y == null) { @@ -988,7 +1022,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public override void SubtractArrays(Complex32[] x, Complex32[] y, Complex32[] result) + public void SubtractArrays(Complex32[] x, Complex32[] y, Complex32[] result) { if (y == null) { @@ -1023,7 +1057,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public override void PointWiseMultiplyArrays(Complex32[] x, Complex32[] y, Complex32[] result) + public void PointWiseMultiplyArrays(Complex32[] x, Complex32[] y, Complex32[] result) { if (y == null) { @@ -1058,7 +1092,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public override void PointWiseDivideArrays(Complex32[] x, Complex32[] y, Complex32[] result) + public void PointWiseDivideArrays(Complex32[] x, Complex32[] y, Complex32[] result) { if (y == null) { @@ -1093,13 +1127,8 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public override void PointWisePowerArrays(Complex32[] x, Complex32[] y, Complex32[] result) + public void PointWisePowerArrays(Complex32[] x, Complex32[] y, Complex32[] result) { - if (_vectorFunctionsMajor != 0 || _vectorFunctionsMinor < 1) - { - base.PointWisePowerArrays(x, y, result); - } - if (y == null) { throw new ArgumentNullException(nameof(y)); @@ -1120,6 +1149,16 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra throw new ArgumentException("The array arguments must have the same length."); } + if (_vectorFunctionsMajor != 0 || _vectorFunctionsMinor < 1) + { + for (int i = 0; i < y.Length; i++) + { + result[i] = Complex32.Pow(x[i], y[i]); + } + + return; + } + SafeNativeMethods.c_vector_power(x.Length, x, y, result); } @@ -1132,7 +1171,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// On output, the matrix contains the eigen vectors. The length of the array must be order * order. /// On output, the eigen values (λ) of matrix in ascending value. The length of the array must . /// On output, the block diagonal eigenvalue matrix. The length of the array must be order * order. - public override void EigenDecomp(bool isSymmetric, int order, Complex32[] matrix, Complex32[] matrixEv, Complex[] vectorEv, Complex32[] matrixD) + public void EigenDecomp(bool isSymmetric, int order, Complex32[] matrix, Complex32[] matrixEv, Complex[] vectorEv, Complex32[] matrixD) { if (matrix == null) { diff --git a/src/Providers.MKL/LinearAlgebra/MklLinearAlgebraProvider.Double.cs b/src/Providers.MKL/LinearAlgebra/MklLinearAlgebraProvider.Double.cs index 9c97195d..6dbc8ae9 100644 --- a/src/Providers.MKL/LinearAlgebra/MklLinearAlgebraProvider.Double.cs +++ b/src/Providers.MKL/LinearAlgebra/MklLinearAlgebraProvider.Double.cs @@ -51,7 +51,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The requested of the matrix. /// [SecuritySafeCritical] - public override double MatrixNorm(Norm norm, int rows, int columns, double[] matrix) + public double MatrixNorm(Norm norm, int rows, int columns, double[] matrix) { if (matrix == null) { @@ -84,7 +84,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The dot product of x and y. /// This is equivalent to the DOT BLAS routine. [SecuritySafeCritical] - public override double DotProduct(double[] x, double[] y) + public double DotProduct(double[] x, double[] y) { if (y == null) { @@ -113,7 +113,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The result of the addition. /// This is similar to the AXPY BLAS routine. [SecuritySafeCritical] - public override void AddVectorToScaledVector(double[] y, double alpha, double[] x, double[] result) + public void AddVectorToScaledVector(double[] y, double alpha, double[] x, double[] result) { if (y == null) { @@ -151,7 +151,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// This result of the scaling. /// This is similar to the SCAL BLAS routine. [SecuritySafeCritical] - public override void ScaleArray(double alpha, double[] x, double[] result) + public void ScaleArray(double alpha, double[] x, double[] result) { if (x == null) { @@ -171,6 +171,24 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra SafeNativeMethods.d_scale(x.Length, alpha, result); } + /// + /// Conjugates an array. Can be used to conjugate a vector and a matrix. + /// + /// The values to conjugate. + /// This result of the conjugation. + public void ConjugateArray(double[] x, double[] result) + { + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (!ReferenceEquals(x, result)) + { + x.CopyTo(result, 0); + } + } + /// /// Multiples two matrices. result = x * y /// @@ -183,7 +201,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// Where to store the result of the multiplication. /// This is a simplified version of the BLAS GEMM routine with alpha /// set to 1.0 and beta set to 0.0, and x and y are not transposed. - public override void MatrixMultiply(double[] x, int rowsX, int columnsX, double[] y, int rowsY, int columnsY, double[] result) + public void MatrixMultiply(double[] x, int rowsX, int columnsX, double[] y, int rowsY, int columnsY, double[] result) { MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, 1.0, x, rowsX, columnsX, y, rowsY, columnsY, 0.0, result); } @@ -203,7 +221,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The value to scale the matrix. /// The c matrix. [SecuritySafeCritical] - public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, int rowsA, int columnsA, double[] b, int rowsB, int columnsB, double beta, double[] c) + public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, int rowsA, int columnsA, double[] b, int rowsB, int columnsB, double beta, double[] c) { if (a == null) { @@ -248,7 +266,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// On exit, it contains the pivot indices. The size of the array must be . /// This is equivalent to the GETRF LAPACK routine. [SecuritySafeCritical] - public override void LUFactor(double[] data, int order, int[] ipiv) + public void LUFactor(double[] data, int order, int[] ipiv) { if (data == null) { @@ -285,7 +303,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. [SecuritySafeCritical] - public override void LUInverse(double[] a, int order) + public void LUInverse(double[] a, int order) { if (a == null) { @@ -323,7 +341,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. [SecuritySafeCritical] - public override void LUInverseFactored(double[] a, int order, int[] ipiv) + public void LUInverseFactored(double[] a, int order, int[] ipiv) { if (a == null) { @@ -367,7 +385,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. [SecuritySafeCritical] - public override void LUSolve(int columnsOfB, double[] a, int order, double[] b) + public void LUSolve(int columnsOfB, double[] a, int order, double[] b) { if (a == null) { @@ -412,7 +430,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRS LAPACK routine. [SecuritySafeCritical] - public override void LUSolveFactored(int columnsOfB, double[] a, int order, int[] ipiv, double[] b) + public void LUSolveFactored(int columnsOfB, double[] a, int order, int[] ipiv, double[] b) { if (a == null) { @@ -460,7 +478,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. [SecuritySafeCritical] - public override void CholeskyFactor(double[] a, int order) + public void CholeskyFactor(double[] a, int order) { if (a == null) { @@ -505,7 +523,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// This is equivalent to the POTRF add POTRS LAPACK routines. /// [SecuritySafeCritical] - public override void CholeskySolve(double[] a, int orderA, double[] b, int columnsB) + public void CholeskySolve(double[] a, int orderA, double[] b, int columnsB) { if (a == null) { @@ -549,7 +567,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. [SecuritySafeCritical] - public override void CholeskySolveFactored(double[] a, int orderA, double[] b, int columnsB) + public void CholeskySolveFactored(double[] a, int orderA, double[] b, int columnsB) { if (a == null) { @@ -592,7 +610,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. [SecuritySafeCritical] - public override void QRFactor(double[] r, int rowsR, int columnsR, double[] q, double[] tau) + public void QRFactor(double[] r, int rowsR, int columnsR, double[] q, double[] tau) { if (r == null) { @@ -640,7 +658,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. [SecuritySafeCritical] - public override void ThinQRFactor(double[] q, int rowsA, int columnsA, double[] r, double[] tau) + public void ThinQRFactor(double[] q, int rowsA, int columnsA, double[] r, double[] tau) { if (r == null) { @@ -687,7 +705,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] - public override void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full) + public void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full) { if (a == null) { @@ -757,7 +775,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] - public override void QRSolveFactored(double[] q, double[] r, int rowsA, int columnsA, double[] tau, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full) + public void QRSolveFactored(double[] q, double[] r, int rowsA, int columnsA, double[] tau, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full) { if (r == null) { @@ -829,7 +847,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra { // we don't have access to the raw Q matrix any more(it is stored in R in the full QR), need to think about this. // let just call the managed version in the meantime. The heavy lifting has already been done. -marcus - base.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin); + ManagedLinearAlgebraProvider.Instance.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin); } } @@ -842,7 +860,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public override void SvdSolve(double[] a, int rowsA, int columnsA, double[] b, int columnsB, double[] x) + public void SvdSolve(double[] a, int rowsA, int columnsA, double[] b, int columnsB, double[] x) { if (a == null) { @@ -879,6 +897,22 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); } + /// + /// Solves A*X=B for X using a previously SVD decomposed matrix. + /// + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The s values returned by . + /// The left singular vectors returned by . + /// The right singular vectors returned by . + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + public void SvdSolveFactored(int rowsA, int columnsA, double[] s, double[] u, double[] vt, double[] b, int columnsB, double[] x) + { + ManagedLinearAlgebraProvider.Instance.SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); + } + /// /// Computes the singular value decomposition of A. /// @@ -893,7 +927,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// right singular vectors. /// This is equivalent to the GESVD LAPACK routine. [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt) + public void SingularValueDecomposition(bool computeVectors, double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt) { if (a == null) { @@ -958,7 +992,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public override void AddArrays(double[] x, double[] y, double[] result) + public void AddArrays(double[] x, double[] y, double[] result) { if (y == null) { @@ -993,7 +1027,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public override void SubtractArrays(double[] x, double[] y, double[] result) + public void SubtractArrays(double[] x, double[] y, double[] result) { if (y == null) { @@ -1028,7 +1062,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public override void PointWiseMultiplyArrays(double[] x, double[] y, double[] result) + public void PointWiseMultiplyArrays(double[] x, double[] y, double[] result) { if (y == null) { @@ -1063,7 +1097,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public override void PointWiseDivideArrays(double[] x, double[] y, double[] result) + public void PointWiseDivideArrays(double[] x, double[] y, double[] result) { if (y == null) { @@ -1098,13 +1132,8 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public override void PointWisePowerArrays(double[] x, double[] y, double[] result) + public void PointWisePowerArrays(double[] x, double[] y, double[] result) { - if (_vectorFunctionsMajor != 0 || _vectorFunctionsMinor < 1) - { - base.PointWisePowerArrays(x, y, result); - } - if (y == null) { throw new ArgumentNullException(nameof(y)); @@ -1125,6 +1154,16 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra throw new ArgumentException("The array arguments must have the same length."); } + if (_vectorFunctionsMajor != 0 || _vectorFunctionsMinor < 1) + { + for (int i = 0; i < y.Length; i++) + { + result[i] = Math.Pow(x[i], y[i]); + } + + return; + } + SafeNativeMethods.d_vector_power(x.Length, x, y, result); } @@ -1137,7 +1176,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// On output, the matrix contains the eigen vectors. The length of the array must be order * order. /// On output, the eigen values (λ) of matrix in ascending value. The length of the array must . /// On output, the block diagonal eigenvalue matrix. The length of the array must be order * order. - public override void EigenDecomp(bool isSymmetric, int order, double[] matrix, double[] matrixEv, Complex[] vectorEv, double[] matrixD) + public void EigenDecomp(bool isSymmetric, int order, double[] matrix, double[] matrixEv, Complex[] vectorEv, double[] matrixD) { if (matrix == null) { diff --git a/src/Providers.MKL/LinearAlgebra/MklLinearAlgebraProvider.Single.cs b/src/Providers.MKL/LinearAlgebra/MklLinearAlgebraProvider.Single.cs index 3cbea25b..d5a544e5 100644 --- a/src/Providers.MKL/LinearAlgebra/MklLinearAlgebraProvider.Single.cs +++ b/src/Providers.MKL/LinearAlgebra/MklLinearAlgebraProvider.Single.cs @@ -51,7 +51,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The requested of the matrix. /// [SecuritySafeCritical] - public override double MatrixNorm(Norm norm, int rows, int columns, float[] matrix) + public double MatrixNorm(Norm norm, int rows, int columns, float[] matrix) { if (matrix == null) { @@ -84,7 +84,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The dot product of x and y. /// This is equivalent to the DOT BLAS routine. [SecuritySafeCritical] - public override float DotProduct(float[] x, float[] y) + public float DotProduct(float[] x, float[] y) { if (y == null) { @@ -113,7 +113,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The result of the addition. /// This is similar to the AXPY BLAS routine. [SecuritySafeCritical] - public override void AddVectorToScaledVector(float[] y, float alpha, float[] x, float[] result) + public void AddVectorToScaledVector(float[] y, float alpha, float[] x, float[] result) { if (y == null) { @@ -151,7 +151,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// This result of the scaling. /// This is similar to the SCAL BLAS routine. [SecuritySafeCritical] - public override void ScaleArray(float alpha, float[] x, float[] result) + public void ScaleArray(float alpha, float[] x, float[] result) { if (x == null) { @@ -171,6 +171,24 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra SafeNativeMethods.s_scale(x.Length, alpha, result); } + /// + /// Conjugates an array. Can be used to conjugate a vector and a matrix. + /// + /// The values to conjugate. + /// This result of the conjugation. + public void ConjugateArray(float[] x, float[] result) + { + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (!ReferenceEquals(x, result)) + { + x.CopyTo(result, 0); + } + } + /// /// Multiples two matrices. result = x * y /// @@ -183,7 +201,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// Where to store the result of the multiplication. /// This is a simplified version of the BLAS GEMM routine with alpha /// set to 1.0f and beta set to 0.0f, and x and y are not transposed. - public override void MatrixMultiply(float[] x, int rowsX, int columnsX, float[] y, int rowsY, int columnsY, float[] result) + public void MatrixMultiply(float[] x, int rowsX, int columnsX, float[] y, int rowsY, int columnsY, float[] result) { MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, 1.0f, x, rowsX, columnsX, y, rowsY, columnsY, 0.0f, result); } @@ -203,7 +221,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The value to scale the matrix. /// The c matrix. [SecuritySafeCritical] - public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, float alpha, float[] a, int rowsA, int columnsA, float[] b, int rowsB, int columnsB, float beta, float[] c) + public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, float alpha, float[] a, int rowsA, int columnsA, float[] b, int rowsB, int columnsB, float beta, float[] c) { if (a == null) { @@ -248,7 +266,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// On exit, it contains the pivot indices. The size of the array must be . /// This is equivalent to the GETRF LAPACK routine. [SecuritySafeCritical] - public override void LUFactor(float[] data, int order, int[] ipiv) + public void LUFactor(float[] data, int order, int[] ipiv) { if (data == null) { @@ -285,7 +303,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. [SecuritySafeCritical] - public override void LUInverse(float[] a, int order) + public void LUInverse(float[] a, int order) { if (a == null) { @@ -323,7 +341,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. [SecuritySafeCritical] - public override void LUInverseFactored(float[] a, int order, int[] ipiv) + public void LUInverseFactored(float[] a, int order, int[] ipiv) { if (a == null) { @@ -367,7 +385,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. [SecuritySafeCritical] - public override void LUSolve(int columnsOfB, float[] a, int order, float[] b) + public void LUSolve(int columnsOfB, float[] a, int order, float[] b) { if (a == null) { @@ -412,7 +430,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRS LAPACK routine. [SecuritySafeCritical] - public override void LUSolveFactored(int columnsOfB, float[] a, int order, int[] ipiv, float[] b) + public void LUSolveFactored(int columnsOfB, float[] a, int order, int[] ipiv, float[] b) { if (a == null) { @@ -460,7 +478,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. [SecuritySafeCritical] - public override void CholeskyFactor(float[] a, int order) + public void CholeskyFactor(float[] a, int order) { if (a == null) { @@ -500,7 +518,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// This is equivalent to the POTRF add POTRS LAPACK routines. /// [SecuritySafeCritical] - public override void CholeskySolve(float[] a, int orderA, float[] b, int columnsB) + public void CholeskySolve(float[] a, int orderA, float[] b, int columnsB) { if (a == null) { @@ -544,7 +562,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. [SecuritySafeCritical] - public override void CholeskySolveFactored(float[] a, int orderA, float[] b, int columnsB) + public void CholeskySolveFactored(float[] a, int orderA, float[] b, int columnsB) { if (a == null) { @@ -587,7 +605,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. [SecuritySafeCritical] - public override void QRFactor(float[] r, int rowsR, int columnsR, float[] q, float[] tau) + public void QRFactor(float[] r, int rowsR, int columnsR, float[] q, float[] tau) { if (r == null) { @@ -635,7 +653,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. [SecuritySafeCritical] - public override void ThinQRFactor(float[] q, int rowsA, int columnsA, float[] r, float[] tau) + public void ThinQRFactor(float[] q, int rowsA, int columnsA, float[] r, float[] tau) { if (r == null) { @@ -682,7 +700,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] - public override void QRSolve(float[] a, int rows, int columns, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full) + public void QRSolve(float[] a, int rows, int columns, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full) { if (a == null) { @@ -752,7 +770,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] - public override void QRSolveFactored(float[] q, float[] r, int rowsA, int columnsA, float[] tau, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full) + public void QRSolveFactored(float[] q, float[] r, int rowsA, int columnsA, float[] tau, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full) { if (r == null) { @@ -824,7 +842,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra { // we don't have access to the raw Q matrix any more(it is stored in R in the full QR), need to think about this. // let just call the managed version in the meantime. The heavy lifting has already been done. -marcus - base.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin); + ManagedLinearAlgebraProvider.Instance.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin); } } @@ -837,7 +855,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public override void SvdSolve(float[] a, int rowsA, int columnsA, float[] b, int columnsB, float[] x) + public void SvdSolve(float[] a, int rowsA, int columnsA, float[] b, int columnsB, float[] x) { if (a == null) { @@ -874,6 +892,22 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); } + /// + /// Solves A*X=B for X using a previously SVD decomposed matrix. + /// + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The s values returned by . + /// The left singular vectors returned by . + /// The right singular vectors returned by . + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + public void SvdSolveFactored(int rowsA, int columnsA, float[] s, float[] u, float[] vt, float[] b, int columnsB, float[] x) + { + ManagedLinearAlgebraProvider.Instance.SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); + } + /// /// Computes the singular value decomposition of A. /// @@ -888,7 +922,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// right singular vectors. /// This is equivalent to the GESVD LAPACK routine. [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt) + public void SingularValueDecomposition(bool computeVectors, float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt) { if (a == null) { @@ -953,7 +987,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public override void AddArrays(float[] x, float[] y, float[] result) + public void AddArrays(float[] x, float[] y, float[] result) { if (y == null) { @@ -988,7 +1022,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public override void SubtractArrays(float[] x, float[] y, float[] result) + public void SubtractArrays(float[] x, float[] y, float[] result) { if (y == null) { @@ -1023,7 +1057,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public override void PointWiseMultiplyArrays(float[] x, float[] y, float[] result) + public void PointWiseMultiplyArrays(float[] x, float[] y, float[] result) { if (y == null) { @@ -1058,7 +1092,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public override void PointWiseDivideArrays(float[] x, float[] y, float[] result) + public void PointWiseDivideArrays(float[] x, float[] y, float[] result) { if (y == null) { @@ -1093,13 +1127,8 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// There is no equivalent BLAS routine, but many libraries /// provide optimized (parallel and/or vectorized) versions of this /// routine. - public override void PointWisePowerArrays(float[] x, float[] y, float[] result) + public void PointWisePowerArrays(float[] x, float[] y, float[] result) { - if (_vectorFunctionsMajor != 0 || _vectorFunctionsMinor < 1) - { - base.PointWisePowerArrays(x, y, result); - } - if (y == null) { throw new ArgumentNullException(nameof(y)); @@ -1120,6 +1149,16 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra throw new ArgumentException("The array arguments must have the same length."); } + if (_vectorFunctionsMajor != 0 || _vectorFunctionsMinor < 1) + { + for (int i = 0; i < y.Length; i++) + { + result[i] = (float)Math.Pow(x[i], y[i]); + } + + return; + } + SafeNativeMethods.s_vector_power(x.Length, x, y, result); } @@ -1132,7 +1171,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// On output, the matrix contains the eigen vectors. The length of the array must be order * order. /// On output, the eigen values (λ) of matrix in ascending value. The length of the array must . /// On output, the block diagonal eigenvalue matrix. The length of the array must be order * order. - public override void EigenDecomp(bool isSymmetric, int order, float[] matrix, float[] matrixEv, Complex[] vectorEv, float[] matrixD) + public void EigenDecomp(bool isSymmetric, int order, float[] matrix, float[] matrixEv, Complex[] vectorEv, float[] matrixD) { if (matrix == null) { diff --git a/src/Providers.MKL/LinearAlgebra/MklLinearAlgebraProvider.cs b/src/Providers.MKL/LinearAlgebra/MklLinearAlgebraProvider.cs index 857d63fb..e3fa9e25 100644 --- a/src/Providers.MKL/LinearAlgebra/MklLinearAlgebraProvider.cs +++ b/src/Providers.MKL/LinearAlgebra/MklLinearAlgebraProvider.cs @@ -28,6 +28,7 @@ // using System; +using System.Numerics; using System.Security; using MathNet.Numerics.Providers.LinearAlgebra; @@ -47,7 +48,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// /// Intel's Math Kernel Library (MKL) linear algebra provider. /// - internal partial class MklLinearAlgebraProvider : ManagedLinearAlgebraProvider, IDisposable + internal sealed partial class MklLinearAlgebraProvider : ILinearAlgebraProvider, IDisposable { const int MinimumCompatibleRevision = 4; @@ -80,7 +81,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// Try to find out whether the provider is available, at least in principle. /// Verification may still fail if available, but it will certainly fail if unavailable. /// - public override bool IsAvailable() + public bool IsAvailable() { return MklProvider.IsAvailable(hintPath: _hintPath); } @@ -90,7 +91,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// If calling this method fails, consider to fall back to alternatives like the managed provider. /// [SecuritySafeCritical] - public override void InitializeVerify() + public void InitializeVerify() { int revision = MklProvider.Load(_hintPath, _consistency, _precision, _accuracy); if (revision < MinimumCompatibleRevision) @@ -114,7 +115,7 @@ namespace MathNet.Numerics.Providers.MKL.LinearAlgebra /// Frees memory buffers, caches and handles allocated in or to the provider. /// Does not unload the provider itself, it is still usable afterwards. /// - public override void FreeResources() + public void FreeResources() { MklProvider.FreeResources(); } diff --git a/src/Providers.MKL/SparseSolver/MklSparseSolverProvider.cs b/src/Providers.MKL/SparseSolver/MklSparseSolverProvider.cs index e6922904..08dc30e5 100644 --- a/src/Providers.MKL/SparseSolver/MklSparseSolverProvider.cs +++ b/src/Providers.MKL/SparseSolver/MklSparseSolverProvider.cs @@ -6,7 +6,7 @@ namespace MathNet.Numerics.Providers.MKL.SparseSolver /// /// Intel's Math Kernel Library (MKL) sparse solver provider. /// - internal partial class MklSparseSolverProvider : ISparseSolverProvider, IDisposable + internal sealed partial class MklSparseSolverProvider : ISparseSolverProvider, IDisposable { const int MinimumCompatibleRevision = 14; diff --git a/src/Providers.OpenBLAS/LinearAlgebra/OpenBlasLinearAlgebraProvider.Complex.cs b/src/Providers.OpenBLAS/LinearAlgebra/OpenBlasLinearAlgebraProvider.Complex.cs index c54ddee3..415b12aa 100644 --- a/src/Providers.OpenBLAS/LinearAlgebra/OpenBlasLinearAlgebraProvider.Complex.cs +++ b/src/Providers.OpenBLAS/LinearAlgebra/OpenBlasLinearAlgebraProvider.Complex.cs @@ -51,7 +51,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The requested of the matrix. /// [SecuritySafeCritical] - public override double MatrixNorm(Norm norm, int rows, int columns, Complex[] matrix) + public double MatrixNorm(Norm norm, int rows, int columns, Complex[] matrix) { if (matrix == null) { @@ -76,6 +76,24 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra return SafeNativeMethods.z_matrix_norm((byte)norm, rows, columns, matrix); } + /// + /// Conjugates an array. Can be used to conjugate a vector and a matrix. + /// + /// The values to conjugate. + /// This result of the conjugation. + public void ConjugateArray(Complex[] x, Complex[] result) + { + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i].Conjugate(); + } + } + /// /// Computes the dot product of x and y. /// @@ -84,7 +102,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The dot product of x and y. /// This is equivalent to the DOT BLAS routine. [SecuritySafeCritical] - public override Complex DotProduct(Complex[] x, Complex[] y) + public Complex DotProduct(Complex[] x, Complex[] y) { if (y == null) { @@ -104,6 +122,194 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra return SafeNativeMethods.z_dot_product(x.Length, x, y); } + /// + /// Does a point wise add of two arrays z = x + y. This can be used + /// to add vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the addition. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void AddArrays(Complex[] x, Complex[] y, Complex[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] + y[i]; + } + } + + /// + /// Does a point wise subtraction of two arrays z = x - y. This can be used + /// to subtract vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the subtraction. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void SubtractArrays(Complex[] x, Complex[] y, Complex[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] - y[i]; + } + }/// + /// Does a point wise multiplication of two arrays z = x * y. This can be used + /// to multiple elements of vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the point wise multiplication. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void PointWiseMultiplyArrays(Complex[] x, Complex[] y, Complex[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] * y[i]; + } + } + + /// + /// Does a point wise division of two arrays z = x / y. This can be used + /// to divide elements of vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the point wise division. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void PointWiseDivideArrays(Complex[] x, Complex[] y, Complex[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] / y[i]; + } + } + + /// + /// Does a point wise power of two arrays z = x ^ y. This can be used + /// to raise elements of vectors or matrices to the powers of another vector or matrix. + /// + /// The array x. + /// The array y. + /// The result of the point wise power. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void PointWisePowerArrays(Complex[] x, Complex[] y, Complex[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = Complex.Pow(x[i], y[i]); + } + } + /// /// Adds a scaled vector to another: result = y + alpha*x. /// @@ -113,7 +319,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The result of the addition. /// This is similar to the AXPY BLAS routine. [SecuritySafeCritical] - public override void AddVectorToScaledVector(Complex[] y, Complex alpha, Complex[] x, Complex[] result) + public void AddVectorToScaledVector(Complex[] y, Complex alpha, Complex[] x, Complex[] result) { if (y == null) { @@ -151,7 +357,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// This result of the scaling. /// This is similar to the SCAL BLAS routine. [SecuritySafeCritical] - public override void ScaleArray(Complex alpha, Complex[] x, Complex[] result) + public void ScaleArray(Complex alpha, Complex[] x, Complex[] result) { if (x == null) { @@ -183,7 +389,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// Where to store the result of the multiplication. /// This is a simplified version of the BLAS GEMM routine with alpha /// set to Complex.One and beta set to Complex.Zero, and x and y are not transposed. - public override void MatrixMultiply(Complex[] x, int rowsX, int columnsX, Complex[] y, int rowsY, int columnsY, Complex[] result) + public void MatrixMultiply(Complex[] x, int rowsX, int columnsX, Complex[] y, int rowsY, int columnsY, Complex[] result) { MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex.One, x, rowsX, columnsX, y, rowsY, columnsY, Complex.Zero, result); } @@ -203,7 +409,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The value to scale the matrix. /// The c matrix. [SecuritySafeCritical] - public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] a, int rowsA, int columnsA, Complex[] b, int rowsB, int columnsB, Complex beta, Complex[] c) + public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] a, int rowsA, int columnsA, Complex[] b, int rowsB, int columnsB, Complex beta, Complex[] c) { if (a == null) { @@ -248,7 +454,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// On exit, it contains the pivot indices. The size of the array must be . /// This is equivalent to the GETRF LAPACK routine. [SecuritySafeCritical] - public override void LUFactor(Complex[] data, int order, int[] ipiv) + public void LUFactor(Complex[] data, int order, int[] ipiv) { if (data == null) { @@ -285,7 +491,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. [SecuritySafeCritical] - public override void LUInverse(Complex[] a, int order) + public void LUInverse(Complex[] a, int order) { if (a == null) { @@ -323,7 +529,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. [SecuritySafeCritical] - public override void LUInverseFactored(Complex[] a, int order, int[] ipiv) + public void LUInverseFactored(Complex[] a, int order, int[] ipiv) { if (a == null) { @@ -367,7 +573,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. [SecuritySafeCritical] - public override void LUSolve(int columnsOfB, Complex[] a, int order, Complex[] b) + public void LUSolve(int columnsOfB, Complex[] a, int order, Complex[] b) { if (a == null) { @@ -412,7 +618,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRS LAPACK routine. [SecuritySafeCritical] - public override void LUSolveFactored(int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b) + public void LUSolveFactored(int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b) { if (a == null) { @@ -465,7 +671,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. [SecuritySafeCritical] - public override void CholeskyFactor(Complex[] a, int order) + public void CholeskyFactor(Complex[] a, int order) { if (a == null) { @@ -505,7 +711,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// This is equivalent to the POTRF add POTRS LAPACK routines. /// [SecuritySafeCritical] - public override void CholeskySolve(Complex[] a, int orderA, Complex[] b, int columnsB) + public void CholeskySolve(Complex[] a, int orderA, Complex[] b, int columnsB) { if (a == null) { @@ -549,7 +755,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. [SecuritySafeCritical] - public override void CholeskySolveFactored(Complex[] a, int orderA, Complex[] b, int columnsB) + public void CholeskySolveFactored(Complex[] a, int orderA, Complex[] b, int columnsB) { if (a == null) { @@ -592,7 +798,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. [SecuritySafeCritical] - public override void QRFactor(Complex[] r, int rowsR, int columnsR, Complex[] q, Complex[] tau) + public void QRFactor(Complex[] r, int rowsR, int columnsR, Complex[] q, Complex[] tau) { if (r == null) { @@ -640,7 +846,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. [SecuritySafeCritical] - public override void ThinQRFactor(Complex[] q, int rowsA, int columnsA, Complex[] r, Complex[] tau) + public void ThinQRFactor(Complex[] q, int rowsA, int columnsA, Complex[] r, Complex[] tau) { if (r == null) { @@ -687,7 +893,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] - public override void QRSolve(Complex[] a, int rows, int columns, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full) + public void QRSolve(Complex[] a, int rows, int columns, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full) { if (a == null) { @@ -757,7 +963,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] - public override void QRSolveFactored(Complex[] q, Complex[] r, int rowsA, int columnsA, Complex[] tau, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full) + public void QRSolveFactored(Complex[] q, Complex[] r, int rowsA, int columnsA, Complex[] tau, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full) { if (r == null) { @@ -829,7 +1035,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra { // we don't have access to the raw Q matrix any more(it is stored in R in the full QR), need to think about this. // let just call the managed version in the meantime. The heavy lifting has already been done. -marcus - base.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin); + ManagedLinearAlgebraProvider.Instance.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin); } } @@ -842,7 +1048,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public override void SvdSolve(Complex[] a, int rowsA, int columnsA, Complex[] b, int columnsB, Complex[] x) + public void SvdSolve(Complex[] a, int rowsA, int columnsA, Complex[] b, int columnsB, Complex[] x) { if (a == null) { @@ -879,6 +1085,22 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); } + /// + /// Solves A*X=B for X using a previously SVD decomposed matrix. + /// + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The s values returned by . + /// The left singular vectors returned by . + /// The right singular vectors returned by . + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + public void SvdSolveFactored(int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, int columnsB, Complex[] x) + { + ManagedLinearAlgebraProvider.Instance.SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); + } + /// /// Computes the singular value decomposition of A. /// @@ -893,7 +1115,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// right singular vectors. /// This is equivalent to the GESVD LAPACK routine. [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt) + public void SingularValueDecomposition(bool computeVectors, Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt) { if (a == null) { @@ -957,7 +1179,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// On output, the matrix contains the eigen vectors. The length of the array must be order * order. /// On output, the eigen values (λ) of matrix in ascending value. The length of the array must . /// On output, the block diagonal eigenvalue matrix. The length of the array must be order * order. - public override void EigenDecomp(bool isSymmetric, int order, Complex[] matrix, Complex[] matrixEv, Complex[] vectorEv, Complex[] matrixD) + public void EigenDecomp(bool isSymmetric, int order, Complex[] matrix, Complex[] matrixEv, Complex[] vectorEv, Complex[] matrixD) { if (matrix == null) { diff --git a/src/Providers.OpenBLAS/LinearAlgebra/OpenBlasLinearAlgebraProvider.Complex32.cs b/src/Providers.OpenBLAS/LinearAlgebra/OpenBlasLinearAlgebraProvider.Complex32.cs index 0a648a42..d199b7f1 100644 --- a/src/Providers.OpenBLAS/LinearAlgebra/OpenBlasLinearAlgebraProvider.Complex32.cs +++ b/src/Providers.OpenBLAS/LinearAlgebra/OpenBlasLinearAlgebraProvider.Complex32.cs @@ -51,7 +51,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The requested of the matrix. /// [SecuritySafeCritical] - public override double MatrixNorm(Norm norm, int rows, int columns, Complex32[] matrix) + public double MatrixNorm(Norm norm, int rows, int columns, Complex32[] matrix) { if (matrix == null) { @@ -76,6 +76,24 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra return SafeNativeMethods.c_matrix_norm((byte)norm, rows, columns, matrix); } + /// + /// Conjugates an array. Can be used to conjugate a vector and a matrix. + /// + /// The values to conjugate. + /// This result of the conjugation. + public void ConjugateArray(Complex32[] x, Complex32[] result) + { + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i].Conjugate(); + } + } + /// /// Computes the dot product of x and y. /// @@ -84,7 +102,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The dot product of x and y. /// This is equivalent to the DOT BLAS routine. [SecuritySafeCritical] - public override Complex32 DotProduct(Complex32[] x, Complex32[] y) + public Complex32 DotProduct(Complex32[] x, Complex32[] y) { if (y == null) { @@ -104,6 +122,196 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra return SafeNativeMethods.c_dot_product(x.Length, x, y); } + /// + /// Does a point wise add of two arrays z = x + y. This can be used + /// to add vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the addition. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void AddArrays(Complex32[] x, Complex32[] y, Complex32[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] + y[i]; + } + } + + /// + /// Does a point wise subtraction of two arrays z = x - y. This can be used + /// to subtract vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the subtraction. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void SubtractArrays(Complex32[] x, Complex32[] y, Complex32[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] - y[i]; + } + } + + /// + /// Does a point wise multiplication of two arrays z = x * y. This can be used + /// to multiple elements of vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the point wise multiplication. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void PointWiseMultiplyArrays(Complex32[] x, Complex32[] y, Complex32[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] * y[i]; + } + } + + /// + /// Does a point wise division of two arrays z = x / y. This can be used + /// to divide elements of vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the point wise division. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void PointWiseDivideArrays(Complex32[] x, Complex32[] y, Complex32[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] / y[i]; + } + } + + /// + /// Does a point wise power of two arrays z = x ^ y. This can be used + /// to raise elements of vectors or matrices to the powers of another vector or matrix. + /// + /// The array x. + /// The array y. + /// The result of the point wise power. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void PointWisePowerArrays(Complex32[] x, Complex32[] y, Complex32[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = Complex32.Pow(x[i], y[i]); + } + } + /// /// Adds a scaled vector to another: result = y + alpha*x. /// @@ -113,7 +321,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The result of the addition. /// This is similar to the AXPY BLAS routine. [SecuritySafeCritical] - public override void AddVectorToScaledVector(Complex32[] y, Complex32 alpha, Complex32[] x, Complex32[] result) + public void AddVectorToScaledVector(Complex32[] y, Complex32 alpha, Complex32[] x, Complex32[] result) { if (y == null) { @@ -151,7 +359,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// This result of the scaling. /// This is similar to the SCAL BLAS routine. [SecuritySafeCritical] - public override void ScaleArray(Complex32 alpha, Complex32[] x, Complex32[] result) + public void ScaleArray(Complex32 alpha, Complex32[] x, Complex32[] result) { if (x == null) { @@ -183,7 +391,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// Where to store the result of the multiplication. /// This is a simplified version of the BLAS GEMM routine with alpha /// set to Complex32.One and beta set to Complex32.Zero, and x and y are not transposed. - public override void MatrixMultiply(Complex32[] x, int rowsX, int columnsX, Complex32[] y, int rowsY, int columnsY, Complex32[] result) + public void MatrixMultiply(Complex32[] x, int rowsX, int columnsX, Complex32[] y, int rowsY, int columnsY, Complex32[] result) { MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex32.One, x, rowsX, columnsX, y, rowsY, columnsY, Complex32.Zero, result); } @@ -203,7 +411,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The value to scale the matrix. /// The c matrix. [SecuritySafeCritical] - public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex32 alpha, Complex32[] a, int rowsA, int columnsA, Complex32[] b, int rowsB, int columnsB, Complex32 beta, Complex32[] c) + public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex32 alpha, Complex32[] a, int rowsA, int columnsA, Complex32[] b, int rowsB, int columnsB, Complex32 beta, Complex32[] c) { if (a == null) { @@ -248,7 +456,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// On exit, it contains the pivot indices. The size of the array must be . /// This is equivalent to the GETRF LAPACK routine. [SecuritySafeCritical] - public override void LUFactor(Complex32[] data, int order, int[] ipiv) + public void LUFactor(Complex32[] data, int order, int[] ipiv) { if (data == null) { @@ -285,7 +493,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. [SecuritySafeCritical] - public override void LUInverse(Complex32[] a, int order) + public void LUInverse(Complex32[] a, int order) { if (a == null) { @@ -323,7 +531,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. [SecuritySafeCritical] - public override void LUInverseFactored(Complex32[] a, int order, int[] ipiv) + public void LUInverseFactored(Complex32[] a, int order, int[] ipiv) { if (a == null) { @@ -367,7 +575,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. [SecuritySafeCritical] - public override void LUSolve(int columnsOfB, Complex32[] a, int order, Complex32[] b) + public void LUSolve(int columnsOfB, Complex32[] a, int order, Complex32[] b) { if (a == null) { @@ -412,7 +620,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRS LAPACK routine. [SecuritySafeCritical] - public override void LUSolveFactored(int columnsOfB, Complex32[] a, int order, int[] ipiv, Complex32[] b) + public void LUSolveFactored(int columnsOfB, Complex32[] a, int order, int[] ipiv, Complex32[] b) { if (a == null) { @@ -460,7 +668,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. [SecuritySafeCritical] - public override void CholeskyFactor(Complex32[] a, int order) + public void CholeskyFactor(Complex32[] a, int order) { if (a == null) { @@ -500,7 +708,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// This is equivalent to the POTRF add POTRS LAPACK routines. /// [SecuritySafeCritical] - public override void CholeskySolve(Complex32[] a, int orderA, Complex32[] b, int columnsB) + public void CholeskySolve(Complex32[] a, int orderA, Complex32[] b, int columnsB) { if (a == null) { @@ -544,7 +752,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. [SecuritySafeCritical] - public override void CholeskySolveFactored(Complex32[] a, int orderA, Complex32[] b, int columnsB) + public void CholeskySolveFactored(Complex32[] a, int orderA, Complex32[] b, int columnsB) { if (a == null) { @@ -587,7 +795,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. [SecuritySafeCritical] - public override void QRFactor(Complex32[] r, int rowsR, int columnsR, Complex32[] q, Complex32[] tau) + public void QRFactor(Complex32[] r, int rowsR, int columnsR, Complex32[] q, Complex32[] tau) { if (r == null) { @@ -635,7 +843,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. [SecuritySafeCritical] - public override void ThinQRFactor(Complex32[] q, int rowsA, int columnsA, Complex32[] r, Complex32[] tau) + public void ThinQRFactor(Complex32[] q, int rowsA, int columnsA, Complex32[] r, Complex32[] tau) { if (r == null) { @@ -682,7 +890,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] - public override void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full) + public void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full) { if (a == null) { @@ -752,7 +960,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] - public override void QRSolveFactored(Complex32[] q, Complex32[] r, int rowsA, int columnsA, Complex32[] tau, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full) + public void QRSolveFactored(Complex32[] q, Complex32[] r, int rowsA, int columnsA, Complex32[] tau, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full) { if (r == null) { @@ -824,7 +1032,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra { // we don't have access to the raw Q matrix any more(it is stored in R in the full QR), need to think about this. // let just call the managed version in the meantime. The heavy lifting has already been done. -marcus - base.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin); + ManagedLinearAlgebraProvider.Instance.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin); } } @@ -837,7 +1045,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public override void SvdSolve(Complex32[] a, int rowsA, int columnsA, Complex32[] b, int columnsB, Complex32[] x) + public void SvdSolve(Complex32[] a, int rowsA, int columnsA, Complex32[] b, int columnsB, Complex32[] x) { if (a == null) { @@ -874,6 +1082,22 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); } + /// + /// Solves A*X=B for X using a previously SVD decomposed matrix. + /// + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The s values returned by . + /// The left singular vectors returned by . + /// The right singular vectors returned by . + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + public void SvdSolveFactored(int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, int columnsB, Complex32[] x) + { + ManagedLinearAlgebraProvider.Instance.SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); + } + /// /// Computes the singular value decomposition of A. /// @@ -888,7 +1112,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// right singular vectors. /// This is equivalent to the GESVD LAPACK routine. [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt) + public void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt) { if (a == null) { @@ -952,7 +1176,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// On output, the matrix contains the eigen vectors. The length of the array must be order * order. /// On output, the eigen values (λ) of matrix in ascending value. The length of the array must . /// On output, the block diagonal eigenvalue matrix. The length of the array must be order * order. - public override void EigenDecomp(bool isSymmetric, int order, Complex32[] matrix, Complex32[] matrixEv, Complex[] vectorEv, Complex32[] matrixD) + public void EigenDecomp(bool isSymmetric, int order, Complex32[] matrix, Complex32[] matrixEv, Complex[] vectorEv, Complex32[] matrixD) { if (matrix == null) { diff --git a/src/Providers.OpenBLAS/LinearAlgebra/OpenBlasLinearAlgebraProvider.Double.cs b/src/Providers.OpenBLAS/LinearAlgebra/OpenBlasLinearAlgebraProvider.Double.cs index cc8d8b4e..18c99a83 100644 --- a/src/Providers.OpenBLAS/LinearAlgebra/OpenBlasLinearAlgebraProvider.Double.cs +++ b/src/Providers.OpenBLAS/LinearAlgebra/OpenBlasLinearAlgebraProvider.Double.cs @@ -51,7 +51,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The requested of the matrix. /// [SecuritySafeCritical] - public override double MatrixNorm(Norm norm, int rows, int columns, double[] matrix) + public double MatrixNorm(Norm norm, int rows, int columns, double[] matrix) { if (matrix == null) { @@ -76,6 +76,24 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra return SafeNativeMethods.d_matrix_norm((byte)norm, rows, columns, matrix); } + /// + /// Conjugates an array. Can be used to conjugate a vector and a matrix. + /// + /// The values to conjugate. + /// This result of the conjugation. + public void ConjugateArray(double[] x, double[] result) + { + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (!ReferenceEquals(x, result)) + { + x.CopyTo(result, 0); + } + } + /// /// Computes the dot product of x and y. /// @@ -84,7 +102,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The dot product of x and y. /// This is equivalent to the DOT BLAS routine. [SecuritySafeCritical] - public override double DotProduct(double[] x, double[] y) + public double DotProduct(double[] x, double[] y) { if (y == null) { @@ -104,6 +122,196 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra return SafeNativeMethods.d_dot_product(x.Length, x, y); } + /// + /// Does a point wise add of two arrays z = x + y. This can be used + /// to add vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the addition. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void AddArrays(double[] x, double[] y, double[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] + y[i]; + } + } + + /// + /// Does a point wise subtraction of two arrays z = x - y. This can be used + /// to subtract vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the subtraction. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void SubtractArrays(double[] x, double[] y, double[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] - y[i]; + } + } + + /// + /// Does a point wise multiplication of two arrays z = x * y. This can be used + /// to multiple elements of vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the point wise multiplication. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void PointWiseMultiplyArrays(double[] x, double[] y, double[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] * y[i]; + } + } + + /// + /// Does a point wise division of two arrays z = x / y. This can be used + /// to divide elements of vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the point wise division. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void PointWiseDivideArrays(double[] x, double[] y, double[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] / y[i]; + } + } + + /// + /// Does a point wise power of two arrays z = x ^ y. This can be used + /// to raise elements of vectors or matrices to the powers of another vector or matrix. + /// + /// The array x. + /// The array y. + /// The result of the point wise power. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void PointWisePowerArrays(double[] x, double[] y, double[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = Math.Pow(x[i], y[i]); + } + } + /// /// Adds a scaled vector to another: result = y + alpha*x. /// @@ -113,7 +321,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The result of the addition. /// This is similar to the AXPY BLAS routine. [SecuritySafeCritical] - public override void AddVectorToScaledVector(double[] y, double alpha, double[] x, double[] result) + public void AddVectorToScaledVector(double[] y, double alpha, double[] x, double[] result) { if (y == null) { @@ -151,7 +359,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// This result of the scaling. /// This is similar to the SCAL BLAS routine. [SecuritySafeCritical] - public override void ScaleArray(double alpha, double[] x, double[] result) + public void ScaleArray(double alpha, double[] x, double[] result) { if (x == null) { @@ -183,7 +391,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// Where to store the result of the multiplication. /// This is a simplified version of the BLAS GEMM routine with alpha /// set to 1.0 and beta set to 0.0, and x and y are not transposed. - public override void MatrixMultiply(double[] x, int rowsX, int columnsX, double[] y, int rowsY, int columnsY, double[] result) + public void MatrixMultiply(double[] x, int rowsX, int columnsX, double[] y, int rowsY, int columnsY, double[] result) { MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, 1.0, x, rowsX, columnsX, y, rowsY, columnsY, 0.0, result); } @@ -203,7 +411,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The value to scale the matrix. /// The c matrix. [SecuritySafeCritical] - public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, int rowsA, int columnsA, double[] b, int rowsB, int columnsB, double beta, double[] c) + public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, int rowsA, int columnsA, double[] b, int rowsB, int columnsB, double beta, double[] c) { if (a == null) { @@ -248,7 +456,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// On exit, it contains the pivot indices. The size of the array must be . /// This is equivalent to the GETRF LAPACK routine. [SecuritySafeCritical] - public override void LUFactor(double[] data, int order, int[] ipiv) + public void LUFactor(double[] data, int order, int[] ipiv) { if (data == null) { @@ -285,7 +493,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. [SecuritySafeCritical] - public override void LUInverse(double[] a, int order) + public void LUInverse(double[] a, int order) { if (a == null) { @@ -323,7 +531,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. [SecuritySafeCritical] - public override void LUInverseFactored(double[] a, int order, int[] ipiv) + public void LUInverseFactored(double[] a, int order, int[] ipiv) { if (a == null) { @@ -367,7 +575,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. [SecuritySafeCritical] - public override void LUSolve(int columnsOfB, double[] a, int order, double[] b) + public void LUSolve(int columnsOfB, double[] a, int order, double[] b) { if (a == null) { @@ -412,7 +620,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRS LAPACK routine. [SecuritySafeCritical] - public override void LUSolveFactored(int columnsOfB, double[] a, int order, int[] ipiv, double[] b) + public void LUSolveFactored(int columnsOfB, double[] a, int order, int[] ipiv, double[] b) { if (a == null) { @@ -460,7 +668,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. [SecuritySafeCritical] - public override void CholeskyFactor(double[] a, int order) + public void CholeskyFactor(double[] a, int order) { if (a == null) { @@ -505,7 +713,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// This is equivalent to the POTRF add POTRS LAPACK routines. /// [SecuritySafeCritical] - public override void CholeskySolve(double[] a, int orderA, double[] b, int columnsB) + public void CholeskySolve(double[] a, int orderA, double[] b, int columnsB) { if (a == null) { @@ -549,7 +757,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. [SecuritySafeCritical] - public override void CholeskySolveFactored(double[] a, int orderA, double[] b, int columnsB) + public void CholeskySolveFactored(double[] a, int orderA, double[] b, int columnsB) { if (a == null) { @@ -592,7 +800,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. [SecuritySafeCritical] - public override void QRFactor(double[] r, int rowsR, int columnsR, double[] q, double[] tau) + public void QRFactor(double[] r, int rowsR, int columnsR, double[] q, double[] tau) { if (r == null) { @@ -640,7 +848,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. [SecuritySafeCritical] - public override void ThinQRFactor(double[] q, int rowsA, int columnsA, double[] r, double[] tau) + public void ThinQRFactor(double[] q, int rowsA, int columnsA, double[] r, double[] tau) { if (r == null) { @@ -687,7 +895,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] - public override void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full) + public void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full) { if (a == null) { @@ -757,7 +965,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] - public override void QRSolveFactored(double[] q, double[] r, int rowsA, int columnsA, double[] tau, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full) + public void QRSolveFactored(double[] q, double[] r, int rowsA, int columnsA, double[] tau, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full) { if (r == null) { @@ -829,7 +1037,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra { // we don't have access to the raw Q matrix any more(it is stored in R in the full QR), need to think about this. // let just call the managed version in the meantime. The heavy lifting has already been done. -marcus - base.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin); + ManagedLinearAlgebraProvider.Instance.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin); } } @@ -842,7 +1050,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public override void SvdSolve(double[] a, int rowsA, int columnsA, double[] b, int columnsB, double[] x) + public void SvdSolve(double[] a, int rowsA, int columnsA, double[] b, int columnsB, double[] x) { if (a == null) { @@ -879,6 +1087,22 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); } + /// + /// Solves A*X=B for X using a previously SVD decomposed matrix. + /// + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The s values returned by . + /// The left singular vectors returned by . + /// The right singular vectors returned by . + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + public void SvdSolveFactored(int rowsA, int columnsA, double[] s, double[] u, double[] vt, double[] b, int columnsB, double[] x) + { + ManagedLinearAlgebraProvider.Instance.SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); + } + /// /// Computes the singular value decomposition of A. /// @@ -893,7 +1117,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// right singular vectors. /// This is equivalent to the GESVD LAPACK routine. [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt) + public void SingularValueDecomposition(bool computeVectors, double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt) { if (a == null) { @@ -957,7 +1181,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// On output, the matrix contains the eigen vectors. The length of the array must be order * order. /// On output, the eigen values (λ) of matrix in ascending value. The length of the array must . /// On output, the block diagonal eigenvalue matrix. The length of the array must be order * order. - public override void EigenDecomp(bool isSymmetric, int order, double[] matrix, double[] matrixEv, Complex[] vectorEv, double[] matrixD) + public void EigenDecomp(bool isSymmetric, int order, double[] matrix, double[] matrixEv, Complex[] vectorEv, double[] matrixD) { if (matrix == null) { diff --git a/src/Providers.OpenBLAS/LinearAlgebra/OpenBlasLinearAlgebraProvider.Single.cs b/src/Providers.OpenBLAS/LinearAlgebra/OpenBlasLinearAlgebraProvider.Single.cs index ee0a7b62..793f4752 100644 --- a/src/Providers.OpenBLAS/LinearAlgebra/OpenBlasLinearAlgebraProvider.Single.cs +++ b/src/Providers.OpenBLAS/LinearAlgebra/OpenBlasLinearAlgebraProvider.Single.cs @@ -51,7 +51,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The requested of the matrix. /// [SecuritySafeCritical] - public override double MatrixNorm(Norm norm, int rows, int columns, float[] matrix) + public double MatrixNorm(Norm norm, int rows, int columns, float[] matrix) { if (matrix == null) { @@ -76,6 +76,24 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra return SafeNativeMethods.s_matrix_norm((byte)norm, rows, columns, matrix); } + /// + /// Conjugates an array. Can be used to conjugate a vector and a matrix. + /// + /// The values to conjugate. + /// This result of the conjugation. + public void ConjugateArray(float[] x, float[] result) + { + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (!ReferenceEquals(x, result)) + { + x.CopyTo(result, 0); + } + } + /// /// Computes the dot product of x and y. /// @@ -84,7 +102,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The dot product of x and y. /// This is equivalent to the DOT BLAS routine. [SecuritySafeCritical] - public override float DotProduct(float[] x, float[] y) + public float DotProduct(float[] x, float[] y) { if (y == null) { @@ -104,6 +122,196 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra return SafeNativeMethods.s_dot_product(x.Length, x, y); } + /// + /// Does a point wise add of two arrays z = x + y. This can be used + /// to add vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the addition. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void AddArrays(float[] x, float[] y, float[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] + y[i]; + } + } + + /// + /// Does a point wise subtraction of two arrays z = x - y. This can be used + /// to subtract vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the subtraction. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void SubtractArrays(float[] x, float[] y, float[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] - y[i]; + } + } + + /// + /// Does a point wise multiplication of two arrays z = x * y. This can be used + /// to multiple elements of vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the point wise multiplication. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void PointWiseMultiplyArrays(float[] x, float[] y, float[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] * y[i]; + } + } + + /// + /// Does a point wise division of two arrays z = x / y. This can be used + /// to divide elements of vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the point wise division. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void PointWiseDivideArrays(float[] x, float[] y, float[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = x[i] / y[i]; + } + } + + /// + /// Does a point wise power of two arrays z = x ^ y. This can be used + /// to raise elements of vectors or matrices to the powers of another vector or matrix. + /// + /// The array x. + /// The array y. + /// The result of the point wise power. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + public void PointWisePowerArrays(float[] x, float[] y, float[] result) + { + if (y == null) + { + throw new ArgumentNullException(nameof(y)); + } + + if (x == null) + { + throw new ArgumentNullException(nameof(x)); + } + + if (result == null) + { + throw new ArgumentNullException(nameof(result)); + } + + if (y.Length != x.Length || y.Length != result.Length) + { + throw new ArgumentException("All vectors must have the same dimensionality."); + } + + for (int i = 0; i < result.Length; i++) + { + result[i] = (float)Math.Pow(x[i], y[i]); + } + } + /// /// Adds a scaled vector to another: result = y + alpha*x. /// @@ -113,7 +321,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The result of the addition. /// This is similar to the AXPY BLAS routine. [SecuritySafeCritical] - public override void AddVectorToScaledVector(float[] y, float alpha, float[] x, float[] result) + public void AddVectorToScaledVector(float[] y, float alpha, float[] x, float[] result) { if (y == null) { @@ -151,7 +359,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// This result of the scaling. /// This is similar to the SCAL BLAS routine. [SecuritySafeCritical] - public override void ScaleArray(float alpha, float[] x, float[] result) + public void ScaleArray(float alpha, float[] x, float[] result) { if (x == null) { @@ -183,7 +391,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// Where to store the result of the multiplication. /// This is a simplified version of the BLAS GEMM routine with alpha /// set to 1.0f and beta set to 0.0f, and x and y are not transposed. - public override void MatrixMultiply(float[] x, int rowsX, int columnsX, float[] y, int rowsY, int columnsY, float[] result) + public void MatrixMultiply(float[] x, int rowsX, int columnsX, float[] y, int rowsY, int columnsY, float[] result) { MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, 1.0f, x, rowsX, columnsX, y, rowsY, columnsY, 0.0f, result); } @@ -203,7 +411,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The value to scale the matrix. /// The c matrix. [SecuritySafeCritical] - public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, float alpha, float[] a, int rowsA, int columnsA, float[] b, int rowsB, int columnsB, float beta, float[] c) + public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, float alpha, float[] a, int rowsA, int columnsA, float[] b, int rowsB, int columnsB, float beta, float[] c) { if (a == null) { @@ -248,7 +456,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// On exit, it contains the pivot indices. The size of the array must be . /// This is equivalent to the GETRF LAPACK routine. [SecuritySafeCritical] - public override void LUFactor(float[] data, int order, int[] ipiv) + public void LUFactor(float[] data, int order, int[] ipiv) { if (data == null) { @@ -285,7 +493,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. [SecuritySafeCritical] - public override void LUInverse(float[] a, int order) + public void LUInverse(float[] a, int order) { if (a == null) { @@ -323,7 +531,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. [SecuritySafeCritical] - public override void LUInverseFactored(float[] a, int order, int[] ipiv) + public void LUInverseFactored(float[] a, int order, int[] ipiv) { if (a == null) { @@ -367,7 +575,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. [SecuritySafeCritical] - public override void LUSolve(int columnsOfB, float[] a, int order, float[] b) + public void LUSolve(int columnsOfB, float[] a, int order, float[] b) { if (a == null) { @@ -412,7 +620,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// On entry the B matrix; on exit the X matrix. /// This is equivalent to the GETRS LAPACK routine. [SecuritySafeCritical] - public override void LUSolveFactored(int columnsOfB, float[] a, int order, int[] ipiv, float[] b) + public void LUSolveFactored(int columnsOfB, float[] a, int order, int[] ipiv, float[] b) { if (a == null) { @@ -460,7 +668,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. [SecuritySafeCritical] - public override void CholeskyFactor(float[] a, int order) + public void CholeskyFactor(float[] a, int order) { if (a == null) { @@ -500,7 +708,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// This is equivalent to the POTRF add POTRS LAPACK routines. /// [SecuritySafeCritical] - public override void CholeskySolve(float[] a, int orderA, float[] b, int columnsB) + public void CholeskySolve(float[] a, int orderA, float[] b, int columnsB) { if (a == null) { @@ -544,7 +752,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. [SecuritySafeCritical] - public override void CholeskySolveFactored(float[] a, int orderA, float[] b, int columnsB) + public void CholeskySolveFactored(float[] a, int orderA, float[] b, int columnsB) { if (a == null) { @@ -587,7 +795,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. [SecuritySafeCritical] - public override void QRFactor(float[] r, int rowsR, int columnsR, float[] q, float[] tau) + public void QRFactor(float[] r, int rowsR, int columnsR, float[] q, float[] tau) { if (r == null) { @@ -635,7 +843,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// to be used by the QR solve routine. /// This is similar to the GEQRF and ORGQR LAPACK routines. [SecuritySafeCritical] - public override void ThinQRFactor(float[] q, int rowsA, int columnsA, float[] r, float[] tau) + public void ThinQRFactor(float[] q, int rowsA, int columnsA, float[] r, float[] tau) { if (r == null) { @@ -682,7 +890,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] - public override void QRSolve(float[] a, int rows, int columns, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full) + public void QRSolve(float[] a, int rows, int columns, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full) { if (a == null) { @@ -752,7 +960,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. [SecuritySafeCritical] - public override void QRSolveFactored(float[] q, float[] r, int rowsA, int columnsA, float[] tau, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full) + public void QRSolveFactored(float[] q, float[] r, int rowsA, int columnsA, float[] tau, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full) { if (r == null) { @@ -824,7 +1032,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra { // we don't have access to the raw Q matrix any more(it is stored in R in the full QR), need to think about this. // let just call the managed version in the meantime. The heavy lifting has already been done. -marcus - base.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin); + ManagedLinearAlgebraProvider.Instance.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin); } } @@ -837,7 +1045,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - public override void SvdSolve(float[] a, int rowsA, int columnsA, float[] b, int columnsB, float[] x) + public void SvdSolve(float[] a, int rowsA, int columnsA, float[] b, int columnsB, float[] x) { if (a == null) { @@ -874,6 +1082,22 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); } + /// + /// Solves A*X=B for X using a previously SVD decomposed matrix. + /// + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The s values returned by . + /// The left singular vectors returned by . + /// The right singular vectors returned by . + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + public void SvdSolveFactored(int rowsA, int columnsA, float[] s, float[] u, float[] vt, float[] b, int columnsB, float[] x) + { + ManagedLinearAlgebraProvider.Instance.SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); + } + /// /// Computes the singular value decomposition of A. /// @@ -888,7 +1112,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// right singular vectors. /// This is equivalent to the GESVD LAPACK routine. [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt) + public void SingularValueDecomposition(bool computeVectors, float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt) { if (a == null) { @@ -952,7 +1176,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// On output, the matrix contains the eigen vectors. The length of the array must be order * order. /// On output, the eigen values (λ) of matrix in ascending value. The length of the array must . /// On output, the block diagonal eigenvalue matrix. The length of the array must be order * order. - public override void EigenDecomp(bool isSymmetric, int order, float[] matrix, float[] matrixEv, Complex[] vectorEv, float[] matrixD) + public void EigenDecomp(bool isSymmetric, int order, float[] matrix, float[] matrixEv, Complex[] vectorEv, float[] matrixD) { if (matrix == null) { diff --git a/src/Providers.OpenBLAS/LinearAlgebra/OpenBlasLinearAlgebraProvider.cs b/src/Providers.OpenBLAS/LinearAlgebra/OpenBlasLinearAlgebraProvider.cs index 6555f7f1..6803566b 100644 --- a/src/Providers.OpenBLAS/LinearAlgebra/OpenBlasLinearAlgebraProvider.cs +++ b/src/Providers.OpenBLAS/LinearAlgebra/OpenBlasLinearAlgebraProvider.cs @@ -53,7 +53,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// /// OpenBLAS linear algebra provider. /// - internal partial class OpenBlasLinearAlgebraProvider : ManagedLinearAlgebraProvider, IDisposable + internal sealed partial class OpenBlasLinearAlgebraProvider : ILinearAlgebraProvider, IDisposable { const int MinimumCompatibleRevision = 1; @@ -69,7 +69,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// Try to find out whether the provider is available, at least in principle. /// Verification may still fail if available, but it will certainly fail if unavailable. /// - public override bool IsAvailable() + public bool IsAvailable() { return OpenBlasProvider.IsAvailable(hintPath: _hintPath); } @@ -78,7 +78,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// Initialize and verify that the provided is indeed available. /// If not, fall back to alternatives like the managed provider /// - public override void InitializeVerify() + public void InitializeVerify() { int revision = OpenBlasProvider.Load(hintPath: _hintPath); if (revision < MinimumCompatibleRevision) @@ -99,7 +99,7 @@ namespace MathNet.Numerics.Providers.OpenBLAS.LinearAlgebra /// Frees memory buffers, caches and handles allocated in or to the provider. /// Does not unload the provider itself, it is still usable afterwards. /// - public override void FreeResources() + public void FreeResources() { OpenBlasProvider.FreeResources(); }