Browse Source

removed work arrays from the LA provider interface

pull/338/head
Marcus Cuda 11 years ago
parent
commit
e1e2bb294e
  1. 44
      MathNet.Numerics.NativeProviders.sln
  2. 136
      src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex.cs
  3. 137
      src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex32.cs
  4. 135
      src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Double.cs
  5. 137
      src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Single.cs
  6. 112
      src/Numerics/Providers/LinearAlgebra/ILinearAlgebraProvider.cs
  7. 291
      src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs
  8. 293
      src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs
  9. 298
      src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs
  10. 291
      src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs
  11. 602
      src/UnitTests/LinearAlgebraProviderTests/Complex/LinearAlgebraProviderTests.cs
  12. 558
      src/UnitTests/LinearAlgebraProviderTests/Complex32/LinearAlgebraProviderTests.cs
  13. 558
      src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs
  14. 558
      src/UnitTests/LinearAlgebraProviderTests/Single/LinearAlgebraProviderTests.cs
  15. 4
      src/UnitTests/UnitTests-CUDA.csproj
  16. 15
      src/UnitTests/UnitTests-MKL.csproj
  17. 4
      src/UnitTests/UnitTests-OpenBLAS.csproj

44
MathNet.Numerics.NativeProviders.sln

@ -26,6 +26,8 @@ Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "OpenBLAS", "src\NativeProvi
EndProject
Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "UnitTests-OpenBLAS", "src\UnitTests\UnitTests-OpenBLAS.csproj", "{96B903EF-3EE1-4569-803C-0482D2F5ED37}"
EndProject
Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "TestData", "src\TestData\TestData.csproj", "{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}"
EndProject
Global
GlobalSection(SolutionConfigurationPlatforms) = preSolution
Debug|Any CPU = Debug|Any CPU
@ -344,6 +346,48 @@ Global
{96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-Signed|Mixed Platforms.Build.0 = Release|Any CPU
{96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-Signed|Win32.ActiveCfg = Release|Any CPU
{96B903EF-3EE1-4569-803C-0482D2F5ED37}.Release-Signed|x64.ActiveCfg = Release|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Debug|Any CPU.Build.0 = Debug|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Debug|Mixed Platforms.ActiveCfg = Debug|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Debug|Mixed Platforms.Build.0 = Debug|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Debug|Win32.ActiveCfg = Debug|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Debug|x64.ActiveCfg = Debug|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release|Any CPU.ActiveCfg = Release|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release|Any CPU.Build.0 = Release|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release|Mixed Platforms.ActiveCfg = Release|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release|Mixed Platforms.Build.0 = Release|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release|Win32.ActiveCfg = Release|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release|x64.ActiveCfg = Release|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-CUDA|Any CPU.ActiveCfg = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-CUDA|Any CPU.Build.0 = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-CUDA|Mixed Platforms.ActiveCfg = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-CUDA|Mixed Platforms.Build.0 = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-CUDA|Win32.ActiveCfg = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-CUDA|Win32.Build.0 = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-CUDA|x64.ActiveCfg = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-CUDA|x64.Build.0 = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-MKL|Any CPU.ActiveCfg = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-MKL|Any CPU.Build.0 = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-MKL|Mixed Platforms.ActiveCfg = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-MKL|Mixed Platforms.Build.0 = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-MKL|Win32.ActiveCfg = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-MKL|Win32.Build.0 = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-MKL|x64.ActiveCfg = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-MKL|x64.Build.0 = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-OpenBLAS|Any CPU.ActiveCfg = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-OpenBLAS|Any CPU.Build.0 = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-OpenBLAS|Mixed Platforms.ActiveCfg = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-OpenBLAS|Mixed Platforms.Build.0 = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-OpenBLAS|Win32.ActiveCfg = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-OpenBLAS|Win32.Build.0 = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-OpenBLAS|x64.ActiveCfg = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-OpenBLAS|x64.Build.0 = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-Signed|Any CPU.ActiveCfg = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-Signed|Any CPU.Build.0 = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-Signed|Mixed Platforms.ActiveCfg = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-Signed|Mixed Platforms.Build.0 = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-Signed|Win32.ActiveCfg = Release-Signed|Any CPU
{AF3253C9-4DB5-45A0-98CF-C105FDA9DA47}.Release-Signed|x64.ActiveCfg = Release-Signed|Any CPU
EndGlobalSection
GlobalSection(SolutionProperties) = preSolution
HideSolutionNode = FALSE

136
src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex.cs

@ -4,7 +4,7 @@
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 Math.NET
// Copyright (c) 2009-2015 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
@ -33,7 +33,6 @@
using System;
using System.Numerics;
using System.Security;
using MathNet.Numerics.LinearAlgebra.Factorization;
using MathNet.Numerics.Properties;
namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
@ -295,73 +294,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
BLAS(SafeNativeMethods.z_lu_inverse_factored(_blasHandle, order, a, ipiv));
}
/// <summary>
/// Computes the inverse of matrix using LU factorization.
/// </summary>
/// <param name="a">The N by N matrix to invert. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="work">Not supported. Should be left null.</param>
/// <remarks>This is equivalent to the GETRF and GETRI LAPACK routines.</remarks>
[SecuritySafeCritical]
public override void LUInverse(Complex[] a, int order, Complex[] work)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (a.Length != order*order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
}
if (work != null)
{
throw new ArgumentException(Resources.UserWorkBufferNotSupported);
}
Solver(SafeNativeMethods.z_lu_inverse(_solverHandle, _blasHandle, order, a));
}
/// <summary>
/// Computes the inverse of a previously factored matrix.
/// </summary>
/// <param name="a">The LU factored N by N matrix. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="work">Not supported. Should be left null.</param>
/// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
[SecuritySafeCritical]
public override void LUInverseFactored(Complex[] a, int order, int[] ipiv, Complex[] work)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (ipiv == null)
{
throw new ArgumentNullException("ipiv");
}
if (a.Length != order*order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
}
if (ipiv.Length != order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv");
}
if (work != null)
{
throw new ArgumentException(Resources.UserWorkBufferNotSupported);
}
BLAS(SafeNativeMethods.z_lu_inverse_factored(_blasHandle, order, a, ipiv));
}
/// <summary>
/// Solves A*X=B for X using LU factorization.
/// </summary>
@ -538,60 +470,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
Solver(SafeNativeMethods.z_cholesky_solve_factored(_solverHandle, orderA, columnsB, a, b));
}
/// <summary>
/// Computes the singular value decomposition of A.
/// </summary>
/// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">If <paramref name="computeVectors"/> is <c>true</c>, on exit U contains the left
/// singular vectors.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
[SecuritySafeCritical]
public override void SingularValueDecomposition(bool computeVectors, Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (s == null)
{
throw new ArgumentNullException("s");
}
if (u == null)
{
throw new ArgumentNullException("u");
}
if (vt == null)
{
throw new ArgumentNullException("vt");
}
if (u.Length != rowsA*rowsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
}
if (vt.Length != columnsA*columnsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt");
}
if (s.Length != Math.Min(rowsA, columnsA))
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
}
SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, null);
}
/// <summary>
/// Solves A*X=B for X using the singular value decomposition of A.
/// </summary>
@ -634,7 +512,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
var clone = new Complex[a.Length];
a.Copy(clone);
SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt, null);
SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt);
SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x);
}
@ -650,10 +528,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
/// singular vectors.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <param name="work">User work buffers are not supported. Should be null.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
[SecuritySafeCritical]
public override void SingularValueDecomposition(bool computeVectors, Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt, Complex[] work)
public override void SingularValueDecomposition(bool computeVectors, Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt)
{
if (a == null)
{
@ -675,11 +552,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
throw new ArgumentNullException("vt");
}
if (work != null)
{
throw new ArgumentException(Resources.UserWorkBufferNotSupported);
}
if (u.Length != rowsA*rowsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
@ -696,7 +568,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
}
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, new Complex[rowsA]);
base.SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt);
else Solver(SafeNativeMethods.z_svd_factor(_solverHandle, computeVectors, rowsA, columnsA, a, s, u, vt));
}
}

137
src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex32.cs

@ -4,7 +4,7 @@
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 Math.NET
// Copyright (c) 2009-2015 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
@ -31,9 +31,7 @@
#if NATIVE
using System;
using System.Numerics;
using System.Security;
using MathNet.Numerics.LinearAlgebra.Factorization;
using MathNet.Numerics.Properties;
namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
@ -295,73 +293,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
BLAS(SafeNativeMethods.c_lu_inverse_factored(_blasHandle, order, a, ipiv));
}
/// <summary>
/// Computes the inverse of matrix using LU factorization.
/// </summary>
/// <param name="a">The N by N matrix to invert. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="work">Not supported. Should be left null.</param>
/// <remarks>This is equivalent to the GETRF and GETRI LAPACK routines.</remarks>
[SecuritySafeCritical]
public override void LUInverse(Complex32[] a, int order, Complex32[] work)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (a.Length != order*order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
}
if (work != null)
{
throw new ArgumentException(Resources.UserWorkBufferNotSupported);
}
Solver(SafeNativeMethods.c_lu_inverse(_solverHandle, _blasHandle, order, a));
}
/// <summary>
/// Computes the inverse of a previously factored matrix.
/// </summary>
/// <param name="a">The LU factored N by N matrix. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="work">Not supported. Should be left null.</param>
/// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
[SecuritySafeCritical]
public override void LUInverseFactored(Complex32[] a, int order, int[] ipiv, Complex32[] work)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (ipiv == null)
{
throw new ArgumentNullException("ipiv");
}
if (a.Length != order*order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
}
if (ipiv.Length != order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv");
}
if (work != null)
{
throw new ArgumentException(Resources.UserWorkBufferNotSupported);
}
BLAS(SafeNativeMethods.c_lu_inverse_factored(_blasHandle, order, a, ipiv));
}
/// <summary>
/// Solves A*X=B for X using LU factorization.
/// </summary>
@ -538,60 +469,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
Solver(SafeNativeMethods.c_cholesky_solve_factored(_solverHandle, orderA, columnsB, a, b));
}
/// <summary>
/// Computes the singular value decomposition of A.
/// </summary>
/// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">If <paramref name="computeVectors"/> is <c>true</c>, on exit U contains the left
/// singular vectors.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
[SecuritySafeCritical]
public override void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (s == null)
{
throw new ArgumentNullException("s");
}
if (u == null)
{
throw new ArgumentNullException("u");
}
if (vt == null)
{
throw new ArgumentNullException("vt");
}
if (u.Length != rowsA*rowsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
}
if (vt.Length != columnsA*columnsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt");
}
if (s.Length != Math.Min(rowsA, columnsA))
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
}
SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, null);
}
/// <summary>
/// Solves A*X=B for X using the singular value decomposition of A.
/// </summary>
@ -634,7 +511,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
var clone = new Complex32[a.Length];
a.Copy(clone);
SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt, null);
SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt);
SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x);
}
@ -650,10 +527,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
/// singular vectors.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <param name="work">Not supported. Should be left null.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
[SecuritySafeCritical]
public override void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] work)
public override void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt)
{
if (a == null)
{
@ -675,11 +551,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
throw new ArgumentNullException("vt");
}
if (work != null)
{
throw new ArgumentException(Resources.UserWorkBufferNotSupported);
}
if (u.Length != rowsA*rowsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
@ -696,7 +567,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
}
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, new Complex32[rowsA]);
base.SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt);
else Solver(SafeNativeMethods.c_svd_factor(_solverHandle, computeVectors, rowsA, columnsA, a, s, u, vt));
}
}

135
src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Double.cs

@ -4,7 +4,7 @@
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 Math.NET
// Copyright (c) 2009-2015 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
@ -31,9 +31,7 @@
#if NATIVE
using System;
using System.Numerics;
using System.Security;
using MathNet.Numerics.LinearAlgebra.Factorization;
using MathNet.Numerics.Properties;
namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
@ -295,73 +293,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
BLAS(SafeNativeMethods.d_lu_inverse_factored(_blasHandle, order, a, ipiv));
}
/// <summary>
/// Computes the inverse of matrix using LU factorization.
/// </summary>
/// <param name="a">The N by N matrix to invert. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="work">Not supported. Should be left null.</param>
/// <remarks>This is equivalent to the GETRF and GETRI LAPACK routines.</remarks>
[SecuritySafeCritical]
public override void LUInverse(double[] a, int order, double[] work)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (a.Length != order*order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
}
if (work != null)
{
throw new ArgumentException(Resources.UserWorkBufferNotSupported);
}
Solver(SafeNativeMethods.d_lu_inverse(_solverHandle, _blasHandle, order, a));
}
/// <summary>
/// Computes the inverse of a previously factored matrix.
/// </summary>
/// <param name="a">The LU factored N by N matrix. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="work">Not supported. Should be left null.</param>
/// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
[SecuritySafeCritical]
public override void LUInverseFactored(double[] a, int order, int[] ipiv, double[] work)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (ipiv == null)
{
throw new ArgumentNullException("ipiv");
}
if (a.Length != order*order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
}
if (ipiv.Length != order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv");
}
if (work != null)
{
throw new ArgumentException(Resources.UserWorkBufferNotSupported);
}
BLAS(SafeNativeMethods.d_lu_inverse_factored(_blasHandle, order, a, ipiv));
}
/// <summary>
/// Solves A*X=B for X using LU factorization.
/// </summary>
@ -538,60 +469,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
Solver(SafeNativeMethods.d_cholesky_solve_factored(_solverHandle, orderA, columnsB, a, b));
}
/// <summary>
/// Computes the singular value decomposition of A.
/// </summary>
/// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">If <paramref name="computeVectors"/> is <c>true</c>, on exit U contains the left
/// singular vectors.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
[SecuritySafeCritical]
public override void SingularValueDecomposition(bool computeVectors, double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (s == null)
{
throw new ArgumentNullException("s");
}
if (u == null)
{
throw new ArgumentNullException("u");
}
if (vt == null)
{
throw new ArgumentNullException("vt");
}
if (u.Length != rowsA*rowsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
}
if (vt.Length != columnsA*columnsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt");
}
if (s.Length != Math.Min(rowsA, columnsA))
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
}
SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, null);
}
/// <summary>
/// Solves A*X=B for X using the singular value decomposition of A.
/// </summary>
@ -650,10 +527,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
/// singular vectors.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <param name="work">Not supported. Should be left null.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
[SecuritySafeCritical]
public override void SingularValueDecomposition(bool computeVectors, double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt, double[] work)
public override void SingularValueDecomposition(bool computeVectors, double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt)
{
if (a == null)
{
@ -675,11 +551,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
throw new ArgumentNullException("vt");
}
if (work != null)
{
throw new ArgumentException(Resources.UserWorkBufferNotSupported);
}
if (u.Length != rowsA*rowsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
@ -696,7 +567,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
}
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, new double[rowsA]);
base.SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt);
else Solver (SafeNativeMethods.d_svd_factor(_solverHandle, computeVectors, rowsA, columnsA, a, s, u, vt));
}
}

137
src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Single.cs

@ -4,7 +4,7 @@
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 Math.NET
// Copyright (c) 2009-2015 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
@ -31,9 +31,7 @@
#if NATIVE
using System;
using System.Numerics;
using System.Security;
using MathNet.Numerics.LinearAlgebra.Factorization;
using MathNet.Numerics.Properties;
namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
@ -295,73 +293,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
BLAS(SafeNativeMethods.s_lu_inverse_factored(_blasHandle, order, a, ipiv));
}
/// <summary>
/// Computes the inverse of matrix using LU factorization.
/// </summary>
/// <param name="a">The N by N matrix to invert. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="work">Not supported. Should be left null.</param>
/// <remarks>This is equivalent to the GETRF and GETRI LAPACK routines.</remarks>
[SecuritySafeCritical]
public override void LUInverse(float[] a, int order, float[] work)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (a.Length != order*order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
}
if (work != null)
{
throw new ArgumentException(Resources.UserWorkBufferNotSupported);
}
Solver(SafeNativeMethods.s_lu_inverse(_solverHandle, _blasHandle, order, a));
}
/// <summary>
/// Computes the inverse of a previously factored matrix.
/// </summary>
/// <param name="a">The LU factored N by N matrix. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="work">Not supported. This should be left null.</param>
/// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
[SecuritySafeCritical]
public override void LUInverseFactored(float[] a, int order, int[] ipiv, float[] work)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (ipiv == null)
{
throw new ArgumentNullException("ipiv");
}
if (a.Length != order*order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
}
if (ipiv.Length != order)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv");
}
if (work != null)
{
throw new ArgumentException(Resources.UserWorkBufferNotSupported);
}
BLAS(SafeNativeMethods.s_lu_inverse_factored(_blasHandle, order, a, ipiv));
}
/// <summary>
/// Solves A*X=B for X using LU factorization.
/// </summary>
@ -538,60 +469,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
Solver(SafeNativeMethods.s_cholesky_solve_factored(_solverHandle, orderA, columnsB, a, b));
}
/// <summary>
/// Computes the singular value decomposition of A.
/// </summary>
/// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">If <paramref name="computeVectors"/> is <c>true</c>, on exit U contains the left
/// singular vectors.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
[SecuritySafeCritical]
public override void SingularValueDecomposition(bool computeVectors, float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (s == null)
{
throw new ArgumentNullException("s");
}
if (u == null)
{
throw new ArgumentNullException("u");
}
if (vt == null)
{
throw new ArgumentNullException("vt");
}
if (u.Length != rowsA*rowsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
}
if (vt.Length != columnsA*columnsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt");
}
if (s.Length != Math.Min(rowsA, columnsA))
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
}
SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, null);
}
/// <summary>
/// Solves A*X=B for X using the singular value decomposition of A.
/// </summary>
@ -634,7 +511,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
var clone = new float[a.Length];
a.Copy(clone);
SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt, null);
SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt);
SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x);
}
@ -650,10 +527,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
/// singular vectors.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <param name="work">Not supported. Should be left null.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
[SecuritySafeCritical]
public override void SingularValueDecomposition(bool computeVectors, float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt, float[] work)
public override void SingularValueDecomposition(bool computeVectors, float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt)
{
if (a == null)
{
@ -675,11 +551,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
throw new ArgumentNullException("vt");
}
if (work != null)
{
throw new ArgumentException(Resources.UserWorkBufferNotSupported);
}
if (u.Length != rowsA*rowsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
@ -696,7 +567,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda
}
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, new float[rowsA]);
base.SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt);
else Solver(SafeNativeMethods.s_svd_factor(_solverHandle, computeVectors, rowsA, columnsA, a, s, u, vt));
}
}

112
src/Numerics/Providers/LinearAlgebra/ILinearAlgebraProvider.cs

@ -269,29 +269,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
/// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
void LUInverseFactored(T[] a, int order, int[] ipiv);
/// <summary>
/// Computes the inverse of matrix using LU factorization.
/// </summary>
/// <param name="a">The N by N matrix to invert. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <remarks>This is equivalent to the GETRF and GETRI LAPACK routines.</remarks>
void LUInverse(T[] a, int order, T[] work);
/// <summary>
/// Computes the inverse of a previously factored matrix.
/// </summary>
/// <param name="a">The LU factored N by N matrix. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
void LUInverseFactored(T[] a, int order, int[] ipiv, T[] work);
/// <summary>
/// Solves A*X=B for X using LU factorization.
/// </summary>
@ -356,23 +333,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
/// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
void QRFactor(T[] a, int rowsA, int columnsA, T[] q, T[] tau);
/// <summary>
/// Computes the full QR factorization of A.
/// </summary>
/// <param name="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.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
/// QR factorization.</param>
/// <param name="tau">A min(m,n) vector. On exit, contains additional information
/// to be used by the QR solve routine.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
void QRFactor(T[] a, int rowsA, int columnsA, T[] q, T[] tau, T[] work);
/// <summary>
/// Computes the thin QR factorization of A where M &gt; N.
/// </summary>
@ -387,23 +347,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
/// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
void ThinQRFactor(T[] a, int rowsA, int columnsA, T[] r, T[] tau);
/// <summary>
/// Computes the thin QR factorization of A where M &gt; N.
/// </summary>
/// <param name="a">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.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="r">On exit, A N by N matrix that holds the R matrix of the
/// QR factorization.</param>
/// <param name="tau">A min(m,n) vector. On exit, contains additional information
/// to be used by the QR solve routine.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
void ThinQRFactor(T[] a, int rowsA, int columnsA, T[] r, T[] tau, T[] work);
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// </summary>
@ -417,22 +360,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
/// <remarks>Rows must be greater or equal to columns.</remarks>
void QRSolve(T[] a, int rows, int columns, T[] b, int columnsB, T[] x, QRMethod method = QRMethod.Full);
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// </summary>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
void QRSolve(T[] a, int rows, int columns, T[] b, int columnsB, T[] x, T[] work, QRMethod method = QRMethod.Full);
/// <summary>
/// Solves A*X=B for X using a previously QR factored matrix.
/// </summary>
@ -450,26 +377,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
void QRSolveFactored(T[] q, T[] r, int rowsA, int columnsA, T[] tau, T[] b, int columnsB, T[] x, QRMethod method = QRMethod.Full);
/// <summary>
/// Solves A*X=B for X using a previously QR factored matrix.
/// </summary>
/// <param name="q">The Q matrix obtained by QR factor. This is only used for the managed provider and can be
/// <c>null</c> for the native provider. The native provider uses the Q portion stored in the R matrix.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(T[],int,int,T[],T[])"/>. </param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
/// and can be <c>null</c> for the managed provider.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">The work array - only used in the native provider. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
void QRSolveFactored(T[] q, T[] r, int rowsA, int columnsA, T[] tau, T[] b, int columnsB, T[] x, T[] work, QRMethod method = QRMethod.Full);
/// <summary>
/// Computes the singular value decomposition of A.
/// </summary>
@ -485,23 +392,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
void SingularValueDecomposition(bool computeVectors, T[] a, int rowsA, int columnsA, T[] s, T[] u, T[] vt);
/// <summary>
/// Computes the singular value decomposition of A.
/// </summary>
/// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value. </param>
/// <param name="u">If <paramref name="computeVectors"/> is <c>true</c>, on exit U contains the left
/// singular vectors.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <param name="work">The work array. On exit, work[0] contains the optimal work size value.
/// </param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
void SingularValueDecomposition(bool computeVectors, T[] a, int rowsA, int columnsA, T[] s, T[] u, T[] vt, T[] work);
/// <summary>
/// Solves A*X=B for X using the singular value decomposition of A.
/// </summary>
@ -520,7 +410,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The s values returned by <see cref="SingularValueDecomposition(bool,T[],int,int,T[],T[],T[])"/>.</param>
/// <param name="u">The left singular vectors returned by <see cref="SingularValueDecomposition(bool,T[],int,int, T[],T[],T[])"/>.</param>
/// <param name="vt">The right singular vectors returned by <see cref="SingularValueDecomposition(bool,T[],int,int,T[],T[],T[],T[])"/>.</param>
/// <param name="vt">The right singular vectors returned by <see cref="SingularValueDecomposition(bool,T[],int,int,T[],T[],T[])"/>.</param>
/// <param name="b">The B matrix</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>

291
src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs

@ -4,7 +4,7 @@
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 Math.NET
// Copyright (c) 2009-2015 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
@ -1058,35 +1058,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
inverse.Copy(a);
}
/// <summary>
/// Computes the inverse of matrix using LU factorization.
/// </summary>
/// <param name="a">The N by N matrix to invert. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <remarks>This is equivalent to the GETRF and GETRI LAPACK routines.</remarks>
public virtual void LUInverse(Complex[] a, int order, Complex[] work)
{
LUInverse(a, order);
}
/// <summary>
/// Computes the inverse of a previously factored matrix.
/// </summary>
/// <param name="a">The LU factored N by N matrix. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
public virtual void LUInverseFactored(Complex[] a, int order, int[] ipiv, Complex[] work)
{
LUInverseFactored(a, order, ipiv);
}
/// <summary>
/// Solves A*X=B for X using LU factorization.
/// </summary>
@ -1469,72 +1440,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
}
var work = columnsR > rowsR ? new Complex[rowsR*rowsR] : new Complex[rowsR*columnsR];
QRFactor(r, rowsR, columnsR, q, tau, work);
}
/// <summary>
/// Computes the QR factorization of A.
/// </summary>
/// <param name="r">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. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
/// QR factorization.</param>
/// <param name="tau">A min(m,n) vector. On exit, contains additional information
/// to be used by the QR solve routine.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
public virtual void QRFactor(Complex[] r, int rowsR, int columnsR, Complex[] q, Complex[] tau, Complex[] work)
{
if (r == null)
{
throw new ArgumentNullException("r");
}
if (q == null)
{
throw new ArgumentNullException("q");
}
if (work == null)
{
throw new ArgumentNullException("q");
}
if (r.Length != rowsR*columnsR)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r");
}
if (tau.Length < Math.Min(rowsR, columnsR))
{
throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau");
}
if (q.Length != rowsR*rowsR)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q");
}
if (columnsR > rowsR)
{
if (work.Length < rowsR*rowsR)
{
work[0] = rowsR*rowsR;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
}
else
{
if (work.Length < rowsR*columnsR)
{
work[0] = rowsR*columnsR;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
}
CommonParallel.For(0, rowsR, (a, b) =>
{
@ -1555,8 +1460,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
{
ComputeQR(work, i, q, i, rowsR, i, rowsR, Control.MaxDegreeOfParallelism);
}
work[0] = columnsR > rowsR ? rowsR*rowsR : rowsR*columnsR;
}
/// <summary>
@ -1599,61 +1502,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
}
var work = new Complex[rowsA*columnsA];
ThinQRFactor(a, rowsA, columnsA, r, tau, work);
}
/// <summary>
/// Computes the QR factorization of A where M &gt; N.
/// </summary>
/// <param name="a">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.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="r">On exit, A N by N matrix that holds the R matrix of the
/// QR factorization.</param>
/// <param name="tau">A min(m,n) vector. On exit, contains additional information
/// to be used by the QR solve routine.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
public virtual void ThinQRFactor(Complex[] a, int rowsA, int columnsA, Complex[] r, Complex[] tau, Complex[] work)
{
if (r == null)
{
throw new ArgumentNullException("r");
}
if (a == null)
{
throw new ArgumentNullException("a");
}
if (work == null)
{
throw new ArgumentNullException("work");
}
if (a.Length != rowsA*columnsA)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "a");
}
if (tau.Length < Math.Min(rowsA, columnsA))
{
throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau");
}
if (r.Length != columnsA*columnsA)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r");
}
if (work.Length < rowsA*columnsA)
{
work[0] = rowsA*columnsA;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
var minmn = Math.Min(rowsA, columnsA);
for (var i = 0; i < minmn; i++)
@ -1684,8 +1532,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
{
ComputeQR(work, i, a, i, rowsA, i, columnsA, Control.MaxDegreeOfParallelism);
}
work[0] = rowsA*columnsA;
}
@ -1815,26 +1661,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
public virtual void QRSolve(Complex[] a, int rows, int columns, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full)
{
var work = new Complex[rows*columns];
QRSolve(a, rows, columns, b, columnsB, x, work, method);
}
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// </summary>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
public virtual void QRSolve(Complex[] a, int rows, int columns, Complex[] b, int columnsB, Complex[] x, Complex[] work, QRMethod method = QRMethod.Full)
{
if (a == null)
{
@ -1851,10 +1677,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
throw new ArgumentNullException("x");
}
if (work == null)
{
throw new ArgumentNullException("work");
}
if (a.Length != rows*columns)
{
@ -1876,11 +1699,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
throw new ArgumentException(Resources.ArgumentArraysSameLength, "x");
}
if (work.Length < rows*columns)
{
work[0] = rows*columns;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
var work = new Complex[rows * columns];
var clone = new Complex[a.Length];
a.Copy(clone);
@ -1897,31 +1716,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
ThinQRFactor(clone, rows, columns, r, work);
QRSolveFactored(clone, r, rows, columns, null, b, columnsB, x, method);
}
work[0] = rows*columns;
}
/// <summary>
/// Solves A*X=B for X using a previously QR factored matrix.
/// </summary>
/// <param name="q">The Q matrix obtained by QR factor. This is only used for the managed provider and can be
/// <c>null</c> for the native provider. The native provider uses the Q portion stored in the R matrix.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(Complex[],int,int,Complex[],Complex[])"/>. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
/// and can be <c>null</c> for the managed provider.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">The work array - only used in the native provider. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
public virtual void QRSolveFactored(Complex[] q, Complex[] r, int rowsR, int columnsR, Complex[] tau, Complex[] b, int columnsB, Complex[] x, Complex[] work, QRMethod method = QRMethod.Full)
{
QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, method);
}
/// <summary>
@ -2102,77 +1896,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
}
var work = new Complex[rowsA];
SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work);
}
/// <summary>
/// Computes the singular value decomposition of A.
/// </summary>
/// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">If <paramref name="computeVectors"/> is <c>true</c>, on exit U contains the left
/// singular vectors.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <param name="work">The work array. Length should be at least <paramref name="rowsA"/>.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
/// <exception cref="NonConvergenceException"></exception>
public virtual void SingularValueDecomposition(bool computeVectors, Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt, Complex[] work)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (s == null)
{
throw new ArgumentNullException("s");
}
if (u == null)
{
throw new ArgumentNullException("u");
}
if (vt == null)
{
throw new ArgumentNullException("vt");
}
if (work == null)
{
throw new ArgumentNullException("work");
}
if (u.Length != rowsA*rowsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
}
if (vt.Length != columnsA*columnsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt");
}
if (s.Length != Math.Min(rowsA, columnsA))
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
}
if (work.Length == 0)
{
throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work");
}
if (work.Length < rowsA)
{
work[0] = rowsA;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
const int maxiter = 1000;
var e = new Complex[columnsA];
@ -2784,10 +2508,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
// a singular vector of length rows+1 when rows < columns. The last element is not used and needs to be removed.
// We should port lapack's svd routine to remove this problem.
Array.Copy(stemp, s, Math.Min(rowsA, columnsA));
// On return the first element of the work array stores the min size of the work array could have been
// work[0] = Math.Max(3 * Math.Min(aRows, aColumns) + Math.Max(aRows, aColumns), 5 * Math.Min(aRows, aColumns));
work[0] = rowsA;
}
/// <summary>
@ -2826,14 +2546,13 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
var work = new Complex[rowsA];
var s = new Complex[Math.Min(rowsA, columnsA)];
var u = new Complex[rowsA*rowsA];
var vt = new Complex[columnsA*columnsA];
var clone = new Complex[a.Length];
a.Copy(clone);
SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt, work);
SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt);
SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x);
}

293
src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs

@ -4,7 +4,7 @@
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 Math.NET
// Copyright (c) 2009-2015 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
@ -1055,35 +1055,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
inverse.Copy(a);
}
/// <summary>
/// Computes the inverse of matrix using LU factorization.
/// </summary>
/// <param name="a">The N by N matrix to invert. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <remarks>This is equivalent to the GETRF and GETRI LAPACK routines.</remarks>
public virtual void LUInverse(Complex32[] a, int order, Complex32[] work)
{
LUInverse(a, order);
}
/// <summary>
/// Computes the inverse of a previously factored matrix.
/// </summary>
/// <param name="a">The LU factored N by N matrix. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
public virtual void LUInverseFactored(Complex32[] a, int order, int[] ipiv, Complex32[] work)
{
LUInverseFactored(a, order, ipiv);
}
/// <summary>
/// Solves A*X=B for X using LU factorization.
/// </summary>
@ -1466,72 +1437,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
}
var work = columnsR > rowsR ? new Complex32[rowsR*rowsR] : new Complex32[rowsR*columnsR];
QRFactor(r, rowsR, columnsR, q, tau, work);
}
/// <summary>
/// Computes the QR factorization of A.
/// </summary>
/// <param name="r">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. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
/// QR factorization.</param>
/// <param name="tau">A min(m,n) vector. On exit, contains additional information
/// to be used by the QR solve routine.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
public virtual void QRFactor(Complex32[] r, int rowsR, int columnsR, Complex32[] q, Complex32[] tau, Complex32[] work)
{
if (r == null)
{
throw new ArgumentNullException("r");
}
if (q == null)
{
throw new ArgumentNullException("q");
}
if (work == null)
{
throw new ArgumentNullException("q");
}
if (r.Length != rowsR*columnsR)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r");
}
if (tau.Length < Math.Min(rowsR, columnsR))
{
throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau");
}
if (q.Length != rowsR*rowsR)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q");
}
if (columnsR > rowsR)
{
if (work.Length < rowsR*rowsR)
{
work[0] = rowsR*rowsR;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
}
else
{
if (work.Length < rowsR*columnsR)
{
work[0] = rowsR*columnsR;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
}
CommonParallel.For(0, rowsR, (a, b) =>
{
@ -1552,8 +1457,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
{
ComputeQR(work, i, q, i, rowsR, i, rowsR, Control.MaxDegreeOfParallelism);
}
work[0] = columnsR > rowsR ? rowsR*rowsR : rowsR*columnsR;
}
/// <summary>
@ -1596,62 +1499,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
}
var work = new Complex32[rowsA*columnsA];
ThinQRFactor(a, rowsA, columnsA, r, tau, work);
}
/// <summary>
/// Computes the QR factorization of A where M &gt; N.
/// </summary>
/// <param name="a">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.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="r">On exit, A N by N matrix that holds the R matrix of the
/// QR factorization.</param>
/// <param name="tau">A min(m,n) vector. On exit, contains additional information
/// to be used by the QR solve routine.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
public virtual void ThinQRFactor(Complex32[] a, int rowsA, int columnsA, Complex32[] r, Complex32[] tau, Complex32[] work)
{
if (r == null)
{
throw new ArgumentNullException("r");
}
if (a == null)
{
throw new ArgumentNullException("a");
}
if (work == null)
{
throw new ArgumentNullException("work");
}
if (a.Length != rowsA*columnsA)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "a");
}
if (tau.Length < Math.Min(rowsA, columnsA))
{
throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau");
}
if (r.Length != columnsA*columnsA)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r");
}
if (work.Length < rowsA*columnsA)
{
work[0] = rowsA*columnsA;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
var minmn = Math.Min(rowsA, columnsA);
for (var i = 0; i < minmn; i++)
{
@ -1681,8 +1529,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
{
ComputeQR(work, i, a, i, rowsA, i, columnsA, Control.MaxDegreeOfParallelism);
}
work[0] = rowsA*columnsA;
}
@ -1812,26 +1658,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
public virtual void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full)
{
var work = new Complex32[rows*columns];
QRSolve(a, rows, columns, b, columnsB, x, work, method);
}
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// </summary>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
public virtual void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x, Complex32[] work, QRMethod method = QRMethod.Full)
{
if (a == null)
{
@ -1848,11 +1674,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
throw new ArgumentNullException("x");
}
if (work == null)
{
throw new ArgumentNullException("work");
}
if (a.Length != rows*columns)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
@ -1873,11 +1694,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
throw new ArgumentException(Resources.RowsLessThanColumns);
}
if (work.Length < rows*columns)
{
work[0] = rows*columns;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
var work = new Complex32[rows * columns];
var clone = new Complex32[a.Length];
a.Copy(clone);
@ -1894,31 +1711,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
ThinQRFactor(clone, rows, columns, r, work);
QRSolveFactored(clone, r, rows, columns, null, b, columnsB, x, method);
}
work[0] = rows*columns;
}
/// <summary>
/// Solves A*X=B for X using a previously QR factored matrix.
/// </summary>
/// <param name="q">The Q matrix obtained by QR factor. This is only used for the managed provider and can be
/// <c>null</c> for the native provider. The native provider uses the Q portion stored in the R matrix.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(Complex32[],int,int,Complex32[],Complex32[])"/>. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
/// and can be <c>null</c> for the managed provider.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">The work array - only used in the native provider. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
public virtual void QRSolveFactored(Complex32[] q, Complex32[] r, int rowsR, int columnsR, Complex32[] tau, Complex32[] b, int columnsB, Complex32[] x, Complex32[] work, QRMethod method = QRMethod.Full)
{
QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, method);
}
/// <summary>
@ -2099,77 +1891,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
}
var work = new Complex32[rowsA];
SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work);
}
/// <summary>
/// Computes the singular value decomposition of A.
/// </summary>
/// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">If <paramref name="computeVectors"/> is <c>true</c>, on exit U contains the left
/// singular vectors.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <param name="work">The work array. Length should be at least <paramref name="rowsA"/>.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
/// <exception cref="NonConvergenceException"></exception>
public virtual void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] work)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (s == null)
{
throw new ArgumentNullException("s");
}
if (u == null)
{
throw new ArgumentNullException("u");
}
if (vt == null)
{
throw new ArgumentNullException("vt");
}
if (work == null)
{
throw new ArgumentNullException("work");
}
if (u.Length != rowsA*rowsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
}
if (vt.Length != columnsA*columnsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt");
}
if (s.Length != Math.Min(rowsA, columnsA))
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
}
if (work.Length == 0)
{
throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work");
}
if (work.Length < rowsA)
{
work[0] = rowsA;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
const int maxiter = 1000;
var e = new Complex32[columnsA];
@ -2781,10 +2503,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
// a singular vector of length rows+1 when rows < columns. The last element is not used and needs to be removed.
// We should port lapack's svd routine to remove this problem.
Array.Copy(stemp, s, Math.Min(rowsA, columnsA));
// On return the first element of the work array stores the min size of the work array could have been
// work[0] = Math.Max(3 * Math.Min(aRows, aColumns) + Math.Max(aRows, aColumns), 5 * Math.Min(aRows, aColumns));
work[0] = rowsA;
}
/// <summary>
@ -2823,14 +2541,13 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
var work = new Complex32[rowsA];
var s = new Complex32[Math.Min(rowsA, columnsA)];
var u = new Complex32[rowsA*rowsA];
var vt = new Complex32[columnsA*columnsA];
var clone = new Complex32[a.Length];
a.Copy(clone);
SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt, work);
SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt);
SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x);
}

298
src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs

@ -4,7 +4,7 @@
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 Math.NET
// Copyright (c) 2009-2015 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
@ -941,35 +941,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
inverse.Copy(a);
}
/// <summary>
/// Computes the inverse of matrix using LU factorization.
/// </summary>
/// <param name="a">The N by N matrix to invert. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <remarks>This is equivalent to the GETRF and GETRI LAPACK routines.</remarks>
public virtual void LUInverse(double[] a, int order, double[] work)
{
LUInverse(a, order);
}
/// <summary>
/// Computes the inverse of a previously factored matrix.
/// </summary>
/// <param name="a">The LU factored N by N matrix. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
public virtual void LUInverseFactored(double[] a, int order, int[] ipiv, double[] work)
{
LUInverseFactored(a, order, ipiv);
}
/// <summary>
/// Solves A*X=B for X using LU factorization.
/// </summary>
@ -1351,83 +1322,15 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q");
}
var work = columnsR > rowsR ? new double[rowsR*rowsR] : new double[rowsR*columnsR];
QRFactor(r, rowsR, columnsR, q, tau, work);
}
/// <summary>
/// Computes the QR factorization of A.
/// </summary>
/// <param name="r">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. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
/// QR factorization.</param>
/// <param name="tau">A min(m,n) vector. On exit, contains additional information
/// to be used by the QR solve routine.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
public virtual void QRFactor(double[] r, int rowsR, int columnsR, double[] q, double[] tau, double[] work)
{
if (r == null)
{
throw new ArgumentNullException("r");
}
if (q == null)
{
throw new ArgumentNullException("q");
}
if (work == null)
{
throw new ArgumentNullException("q");
}
if (r.Length != rowsR*columnsR)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r");
}
if (tau.Length < Math.Min(rowsR, columnsR))
{
throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau");
}
if (q.Length != rowsR*rowsR)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q");
}
if (columnsR > rowsR)
{
if (work.Length < rowsR*rowsR)
{
work[0] = rowsR*rowsR;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
}
else
{
if (work.Length < rowsR*columnsR)
{
work[0] = rowsR*columnsR;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
}
CommonParallel.For(0, rowsR, (a, b) =>
{
for (int i = a; i < b; i++)
for (var i = a; i < b; i++)
{
q[(i*rowsR) + i] = 1.0;
}
});
var work = columnsR > rowsR ? new double[rowsR * rowsR] : new double[rowsR * columnsR];
var minmn = Math.Min(rowsR, columnsR);
for (var i = 0; i < minmn; i++)
{
@ -1439,8 +1342,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
{
ComputeQR(work, i, q, i, rowsR, i, rowsR, Control.MaxDegreeOfParallelism);
}
work[0] = columnsR > rowsR ? rowsR*rowsR : rowsR*columnsR;
}
/// <summary>
@ -1483,62 +1384,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
}
var work = new double[rowsA*columnsA];
ThinQRFactor(a, rowsA, columnsA, r, tau, work);
}
/// <summary>
/// Computes the thin QR factorization of A where M &gt; N.
/// </summary>
/// <param name="a">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.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="r">On exit, A N by N matrix that holds the R matrix of the
/// QR factorization.</param>
/// <param name="tau">A min(m,n) vector. On exit, contains additional information
/// to be used by the QR solve routine.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
public virtual void ThinQRFactor(double[] a, int rowsA, int columnsA, double[] r, double[] tau, double[] work)
{
if (r == null)
{
throw new ArgumentNullException("r");
}
if (a == null)
{
throw new ArgumentNullException("a");
}
if (work == null)
{
throw new ArgumentNullException("work");
}
if (a.Length != rowsA*columnsA)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "a");
}
if (tau.Length < Math.Min(rowsA, columnsA))
{
throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau");
}
if (r.Length != columnsA*columnsA)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r");
}
if (work.Length < rowsA*columnsA)
{
work[0] = rowsA*columnsA;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
var minmn = Math.Min(rowsA, columnsA);
for (var i = 0; i < minmn; i++)
{
@ -1568,8 +1414,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
{
ComputeQR(work, i, a, i, rowsA, i, columnsA, Control.MaxDegreeOfParallelism);
}
work[0] = rowsA*columnsA;
}
#region QR Factor Helper functions
@ -1699,26 +1543,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
public virtual void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full)
{
var work = new double[rows*columns];
QRSolve(a, rows, columns, b, columnsB, x, work, method);
}
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// </summary>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
public virtual void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x, double[] work, QRMethod method = QRMethod.Full)
{
if (a == null)
{
@ -1735,11 +1559,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
throw new ArgumentNullException("x");
}
if (work == null)
{
throw new ArgumentNullException("work");
}
if (a.Length != rows*columns)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
@ -1760,11 +1579,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
throw new ArgumentException(Resources.RowsLessThanColumns);
}
if (work.Length < rows*columns)
{
work[0] = rows*columns;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
var work = new double[rows * columns];
var clone = new double[a.Length];
a.Copy(clone);
@ -1781,31 +1596,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
ThinQRFactor(clone, rows, columns, r, work);
QRSolveFactored(clone, r, rows, columns, null, b, columnsB, x, method);
}
work[0] = rows*columns;
}
/// <summary>
/// Solves A*X=B for X using a previously QR factored matrix.
/// </summary>
/// <param name="q">The Q matrix obtained by QR factor. This is only used for the managed provider and can be
/// <c>null</c> for the native provider. The native provider uses the Q portion stored in the R matrix.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(double[],int,int,double[],double[])"/>. </param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
/// and can be <c>null</c> for the managed provider.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">The work array - only used in the native provider. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
public virtual void QRSolveFactored(double[] q, double[] r, int rowsA, int columnsA, double[] tau, double[] b, int columnsB, double[] x, double[] work, QRMethod method = QRMethod.Full)
{
QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, method);
}
/// <summary>
@ -1986,76 +1776,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
}
var work = new double[rowsA];
SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work);
}
/// <summary>
/// Computes the singular value decomposition of A.
/// </summary>
/// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">If <paramref name="computeVectors"/> is <c>true</c>, on exit U contains the left
/// singular vectors.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <param name="work">The work array. Length should be at least <paramref name="rowsA"/>.</param>
/// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
/// <exception cref="NonConvergenceException"></exception>
public virtual void SingularValueDecomposition(bool computeVectors, double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt, double[] work)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (s == null)
{
throw new ArgumentNullException("s");
}
if (u == null)
{
throw new ArgumentNullException("u");
}
if (vt == null)
{
throw new ArgumentNullException("vt");
}
if (work == null)
{
throw new ArgumentNullException("work");
}
if (u.Length != rowsA*rowsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
}
if (vt.Length != columnsA*columnsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt");
}
if (s.Length != Math.Min(rowsA, columnsA))
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
}
if (work.Length == 0)
{
throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work");
}
if (work.Length < rowsA)
{
work[0] = rowsA;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
const int maxiter = 1000;
@ -2672,10 +2392,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
// a singular vector of length rows+1 when rows < columns. The last element is not used and needs to be removed.
// We should port lapack's svd routine to remove this problem.
Buffer.BlockCopy(stemp, 0, s, 0, Math.Min(rowsA, columnsA)*Constants.SizeOfDouble);
// On return the first element of the work array stores the min size of the work array could have been
// work[0] = Math.Max(3 * Math.Min(aRows, aColumns) + Math.Max(aRows, aColumns), 5 * Math.Min(aRows, aColumns));
work[0] = rowsA;
}
/// <summary>
@ -2770,15 +2486,13 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
var work = new double[rowsA];
var s = new double[Math.Min(rowsA, columnsA)];
var u = new double[rowsA*rowsA];
var vt = new double[columnsA*columnsA];
var clone = new double[a.Length];
Buffer.BlockCopy(a, 0, clone, 0, a.Length*Constants.SizeOfDouble);
SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt, work);
SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt);
SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x);
}

291
src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs

@ -4,7 +4,7 @@
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 Math.NET
// Copyright (c) 2009-2015 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
@ -941,35 +941,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
inverse.Copy(a);
}
/// <summary>
/// Computes the inverse of matrix using LU factorization.
/// </summary>
/// <param name="a">The N by N matrix to invert. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <remarks>This is equivalent to the GETRF and GETRI LAPACK routines.</remarks>
public virtual void LUInverse(float[] a, int order, float[] work)
{
LUInverse(a, order);
}
/// <summary>
/// Computes the inverse of a previously factored matrix.
/// </summary>
/// <param name="a">The LU factored N by N matrix. Contains the inverse On exit.</param>
/// <param name="order">The order of the square matrix <paramref name="a"/>.</param>
/// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
public virtual void LUInverseFactored(float[] a, int order, int[] ipiv, float[] work)
{
LUInverseFactored(a, order, ipiv);
}
/// <summary>
/// Solves A*X=B for X using LU factorization.
/// </summary>
@ -1352,72 +1323,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
}
var work = columnsR > rowsR ? new float[rowsR*rowsR] : new float[rowsR*columnsR];
QRFactor(r, rowsR, columnsR, q, tau, work);
}
/// <summary>
/// Computes the QR factorization of A.
/// </summary>
/// <param name="r">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. </param>
/// <param name="rowsR">The number of rows in the A matrix.</param>
/// <param name="columnsR">The number of columns in the A matrix.</param>
/// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
/// QR factorization.</param>
/// <param name="tau">A min(m,n) vector. On exit, contains additional information
/// to be used by the QR solve routine.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
public virtual void QRFactor(float[] r, int rowsR, int columnsR, float[] q, float[] tau, float[] work)
{
if (r == null)
{
throw new ArgumentNullException("r");
}
if (q == null)
{
throw new ArgumentNullException("q");
}
if (work == null)
{
throw new ArgumentNullException("q");
}
if (r.Length != rowsR*columnsR)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r");
}
if (tau.Length < Math.Min(rowsR, columnsR))
{
throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau");
}
if (q.Length != rowsR*rowsR)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q");
}
if (columnsR > rowsR)
{
if (work.Length < rowsR*rowsR)
{
work[0] = rowsR*rowsR;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
}
else
{
if (work.Length < rowsR*columnsR)
{
work[0] = rowsR*columnsR;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
}
CommonParallel.For(0, rowsR, (a, b) =>
{
@ -1438,8 +1343,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
{
ComputeQR(work, i, q, i, rowsR, i, rowsR, Control.MaxDegreeOfParallelism);
}
work[0] = columnsR > rowsR ? rowsR*rowsR : rowsR*columnsR;
}
/// <summary>
@ -1482,61 +1385,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
}
var work = new float[rowsA*columnsA];
ThinQRFactor(a, rowsA, columnsA, r, tau, work);
}
/// <summary>
/// Computes the QR factorization of A where M &gt; N.
/// </summary>
/// <param name="a">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.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="r">On exit, A N by N matrix that holds the R matrix of the
/// QR factorization.</param>
/// <param name="tau">A min(m,n) vector. On exit, contains additional information
/// to be used by the QR solve routine.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
public virtual void ThinQRFactor(float[] a, int rowsA, int columnsA, float[] r, float[] tau, float[] work)
{
if (r == null)
{
throw new ArgumentNullException("r");
}
if (a == null)
{
throw new ArgumentNullException("a");
}
if (work == null)
{
throw new ArgumentNullException("work");
}
if (a.Length != rowsA*columnsA)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "a");
}
if (tau.Length < Math.Min(rowsA, columnsA))
{
throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau");
}
if (r.Length != columnsA*columnsA)
{
throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "columnsA * columnsA"), "r");
}
if (work.Length < rowsA*columnsA)
{
work[0] = rowsA*columnsA;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
var minmn = Math.Min(rowsA, columnsA);
for (var i = 0; i < minmn; i++)
@ -1567,8 +1415,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
{
ComputeQR(work, i, a, i, rowsA, i, columnsA, Control.MaxDegreeOfParallelism);
}
work[0] = rowsA*columnsA;
}
@ -1687,22 +1533,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
#endregion
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
/// </summary>
/// <param name="a">The A matrix.</param>
/// <param name="rows">The number of rows in the A matrix.</param>
/// <param name="columns">The number of columns in the A matrix.</param>
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
public virtual void QRSolve(float[] a, int rows, int columns, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full)
{
var work = new float[rows*columns];
QRSolve(a, rows, columns, b, columnsB, x, work, method);
}
/// <summary>
/// Solves A*X=B for X using QR factorization of A.
@ -1713,12 +1543,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
/// <param name="b">The B matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">The work array. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
public virtual void QRSolve(float[] a, int rows, int columns, float[] b, int columnsB, float[] x, float[] work, QRMethod method = QRMethod.Full)
public virtual void QRSolve(float[] a, int rows, int columns, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full)
{
if (a == null)
{
@ -1735,11 +1562,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
throw new ArgumentNullException("x");
}
if (work == null)
{
throw new ArgumentNullException("work");
}
if (a.Length != rows*columns)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
@ -1760,11 +1582,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
throw new ArgumentException(Resources.RowsLessThanColumns);
}
if (work.Length < rows*columns)
{
work[0] = rows*columns;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
var work = new float[rows * columns];
var clone = new float[a.Length];
a.Copy(clone);
@ -1781,31 +1599,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
ThinQRFactor(clone, rows, columns, r, work);
QRSolveFactored(clone, r, rows, columns, null, b, columnsB, x, method);
}
work[0] = rows*columns;
}
/// <summary>
/// Solves A*X=B for X using a previously QR factored matrix.
/// </summary>
/// <param name="q">The Q matrix obtained by QR factor. This is only used for the managed provider and can be
/// <c>null</c> for the native provider. The native provider uses the Q portion stored in the R matrix.</param>
/// <param name="r">The R matrix obtained by calling <see cref="QRFactor(float[],int,int,float[],float[])"/>. </param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="tau">Contains additional information on Q. Only used for the native solver
/// and can be <c>null</c> for the managed provider.</param>
/// <param name="b">On entry the B matrix; on exit the X matrix.</param>
/// <param name="columnsB">The number of columns of B.</param>
/// <param name="x">On exit, the solution matrix.</param>
/// <param name="work">The work array - only used in the native provider. The array must have a length of at least N,
/// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
/// work size value.</param>
/// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
/// <remarks>Rows must be greater or equal to columns.</remarks>
public virtual void QRSolveFactored(float[] q, float[] r, int rowsA, int columnsA, float[] tau, float[] b, int columnsB, float[] x, float[] work, QRMethod method = QRMethod.Full)
{
QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, method);
}
/// <summary>
@ -1986,76 +1779,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
}
var work = new float[rowsA];
SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work);
}
/// <summary>
/// Computes the singular value decomposition of A.
/// </summary>
/// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
/// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
/// <param name="rowsA">The number of rows in the A matrix.</param>
/// <param name="columnsA">The number of columns in the A matrix.</param>
/// <param name="s">The singular values of A in ascending value.</param>
/// <param name="u">If <paramref name="computeVectors"/> is <c>true</c>, on exit U contains the left
/// singular vectors.</param>
/// <param name="vt">If <paramref name="computeVectors"/> is <c>true</c>, on exit VT contains the transposed
/// right singular vectors.</param>
/// <param name="work">The work array. Length should be at least <paramref name="rowsA"/>.</param>
/// <exception cref="NonConvergenceException"></exception>
public virtual void SingularValueDecomposition(bool computeVectors, float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt, float[] work)
{
if (a == null)
{
throw new ArgumentNullException("a");
}
if (s == null)
{
throw new ArgumentNullException("s");
}
if (u == null)
{
throw new ArgumentNullException("u");
}
if (vt == null)
{
throw new ArgumentNullException("vt");
}
if (work == null)
{
throw new ArgumentNullException("work");
}
if (u.Length != rowsA*rowsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "u");
}
if (vt.Length != columnsA*columnsA)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt");
}
if (s.Length != Math.Min(rowsA, columnsA))
{
throw new ArgumentException(Resources.ArgumentArraysSameLength, "s");
}
if (work.Length == 0)
{
throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work");
}
if (work.Length < rowsA)
{
work[0] = rowsA;
throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
}
const int maxiter = 1000;
var e = new float[columnsA];
@ -2673,10 +2397,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
// a singular vector of length rows+1 when rows < columns. The last element is not used and needs to be removed.
// We should port lapack's svd routine to remove this problem.
Buffer.BlockCopy(stemp, 0, s, 0, Math.Min(rowsA, columnsA)*Constants.SizeOfFloat);
// On return the first element of the work array stores the min size of the work array could have been
// work[0] = Math.Max(3 * Math.Min(aRows, aColumns) + Math.Max(aRows, aColumns), 5 * Math.Min(aRows, aColumns));
work[0] = rowsA;
}
/// <summary>
@ -2772,14 +2492,13 @@ namespace MathNet.Numerics.Providers.LinearAlgebra
throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
}
var work = new float[rowsA];
var s = new float[Math.Min(rowsA, columnsA)];
var u = new float[rowsA*rowsA];
var vt = new float[columnsA*columnsA];
var clone = new float[a.Length];
Buffer.BlockCopy(a, 0, clone, 0, a.Length*Constants.SizeOfFloat);
SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt, work);
SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt);
SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x);
}

602
src/UnitTests/LinearAlgebraProviderTests/Complex/LinearAlgebraProviderTests.cs

@ -4,7 +4,7 @@
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 Math.NET
// Copyright (c) 2009-2015 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
@ -433,72 +433,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex
AssertHelpers.AlmostEqualRelative(a[8], -0.113636363636364, 13);
}
#if ! MKL
/// <summary>
/// Can compute the inverse of a matrix using LU factorization
/// with a work array.
/// </summary>
[Test]
public void CanComputeLuInverseWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new Complex[matrix.RowCount*matrix.RowCount];
Array.Copy(matrix.Values, a, a.Length);
#if CUDA
Complex[] work = null;
#else
var work = new Complex[matrix.RowCount];
#endif
Control.LinearAlgebraProvider.LUInverse(a, matrix.RowCount, work);
AssertHelpers.AlmostEqualRelative(a[0], -0.454545454545454, 13);
AssertHelpers.AlmostEqualRelative(a[1], -0.909090909090908, 13);
AssertHelpers.AlmostEqualRelative(a[2], 0.454545454545454, 13);
AssertHelpers.AlmostEqualRelative(a[3], -0.340909090909090, 13);
AssertHelpers.AlmostEqualRelative(a[4], -2.045454545454543, 13);
AssertHelpers.AlmostEqualRelative(a[5], 1.477272727272726, 13);
AssertHelpers.AlmostEqualRelative(a[6], -0.113636363636364, 13);
AssertHelpers.AlmostEqualRelative(a[7], 0.227272727272727, 13);
AssertHelpers.AlmostEqualRelative(a[8], -0.113636363636364, 13);
}
#endif
#if ! MKL
/// <summary>
/// Can compute the inverse of a matrix using LU factorization
/// using a previously factored matrix with a work array.
/// </summary>
[Test]
public void CanComputeLuInverseOnFactoredMatrixWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new Complex[matrix.RowCount*matrix.RowCount];
Array.Copy(matrix.Values, a, a.Length);
var ipiv = new int[matrix.RowCount];
Control.LinearAlgebraProvider.LUFactor(a, matrix.RowCount, ipiv);
#if CUDA
Complex[] work = null;
#else
var work = new Complex[matrix.RowCount];
#endif
Control.LinearAlgebraProvider.LUInverseFactored(a, matrix.RowCount, ipiv, work);
AssertHelpers.AlmostEqualRelative(a[0], -0.454545454545454, 13);
AssertHelpers.AlmostEqualRelative(a[1], -0.909090909090908, 13);
AssertHelpers.AlmostEqualRelative(a[2], 0.454545454545454, 13);
AssertHelpers.AlmostEqualRelative(a[3], -0.340909090909090, 13);
AssertHelpers.AlmostEqualRelative(a[4], -2.045454545454543, 13);
AssertHelpers.AlmostEqualRelative(a[5], 1.477272727272726, 13);
AssertHelpers.AlmostEqualRelative(a[6], -0.113636363636364, 13);
AssertHelpers.AlmostEqualRelative(a[7], 0.227272727272727, 13);
AssertHelpers.AlmostEqualRelative(a[8], -0.113636363636364, 13);
}
#endif
/// <summary>
/// Can solve Ax=b using LU factorization.
/// </summary>
@ -696,92 +630,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex
}
}
#if ! MKL
/// <summary>
/// Can compute QR factorization of a square matrix using a work array.
/// </summary>
[Test]
public void CanComputeQRFactorSquareMatrixWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var r = new Complex[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, r, r.Length);
var tau = new Complex[3];
var q = new Complex[matrix.RowCount*matrix.RowCount];
var work = new Complex[matrix.ColumnCount*Control.BlockSize];
Control.LinearAlgebraProvider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau, work);
var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q);
var mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r).UpperTriangle();
var a = mq*mr;
for (var row = 0; row < matrix.RowCount; row++)
{
for (var col = 0; col < matrix.ColumnCount; col++)
{
AssertHelpers.AlmostEqualRelative(matrix[row, col], a[row, col], 14);
}
}
}
/// <summary>
/// Can compute QR factorization of a tall matrix using a work matrix.
/// </summary>
[Test]
public void CanComputeQRFactorTallMatrixWithWorkArray()
{
var matrix = _matrices["Tall3x2"];
var r = new Complex[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, r, r.Length);
var tau = new Complex[3];
var q = new Complex[matrix.RowCount*matrix.RowCount];
var work = new Complex[matrix.ColumnCount*Control.BlockSize];
Control.LinearAlgebraProvider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau, work);
var mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r).UpperTriangle();
var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q);
var a = mq*mr;
for (var row = 0; row < matrix.RowCount; row++)
{
for (var col = 0; col < matrix.ColumnCount; col++)
{
AssertHelpers.AlmostEqualRelative(matrix[row, col], a[row, col], 14);
}
}
}
/// <summary>
/// Can compute QR factorization of a wide matrix using a work matrix.
/// </summary>
[Test]
public void CanComputeQRFactorWideMatrixWithWorkArray()
{
var matrix = _matrices["Wide2x3"];
var r = new Complex[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, r, r.Length);
var tau = new Complex[3];
var q = new Complex[matrix.RowCount*matrix.RowCount];
var work = new Complex[matrix.ColumnCount*Control.BlockSize];
Control.LinearAlgebraProvider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau, work);
var mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r).UpperTriangle();
var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q);
var a = mq*mr;
for (var row = 0; row < matrix.RowCount; row++)
{
for (var col = 0; col < matrix.ColumnCount; col++)
{
AssertHelpers.AlmostEqualRelative(matrix[row, col], a[row, col], 14);
}
}
}
#endif
/// <summary>
/// Can compute thin QR factorization of a square matrix.
/// </summary>
@ -836,63 +684,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex
}
}
#if ! MKL
/// <summary>
/// Can compute thin QR factorization of a square matrix using a work array.
/// </summary>
[Test]
public void CanComputeThinQRFactorSquareMatrixWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var r = new Complex[matrix.ColumnCount*matrix.ColumnCount];
var tau = new Complex[3];
var q = new Complex[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, q, q.Length);
var work = new Complex[matrix.ColumnCount*Control.BlockSize];
Control.LinearAlgebraProvider.ThinQRFactor(q, matrix.RowCount, matrix.ColumnCount, r, tau, work);
var mq = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, q);
var mr = new DenseMatrix(matrix.ColumnCount, matrix.ColumnCount, r);
var a = mq*mr;
for (var row = 0; row < matrix.RowCount; row++)
{
for (var col = 0; col < matrix.ColumnCount; col++)
{
AssertHelpers.AlmostEqualRelative(matrix[row, col], a[row, col], 14);
}
}
}
/// <summary>
/// Can compute thin QR factorization of a tall matrix using a work matrix.
/// </summary>
[Test]
public void CanComputeThinQRFactorTallMatrixWithWorkArray()
{
var matrix = _matrices["Tall3x2"];
var r = new Complex[matrix.ColumnCount*matrix.ColumnCount];
var tau = new Complex[3];
var q = new Complex[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, q, q.Length);
var work = new Complex[matrix.ColumnCount*Control.BlockSize];
Control.LinearAlgebraProvider.ThinQRFactor(q, matrix.RowCount, matrix.ColumnCount, r, tau, work);
var mq = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, q);
var mr = new DenseMatrix(matrix.ColumnCount, matrix.ColumnCount, r);
var a = mq*mr;
for (var row = 0; row < matrix.RowCount; row++)
{
for (var col = 0; col < matrix.ColumnCount; col++)
{
AssertHelpers.AlmostEqualRelative(matrix[row, col], a[row, col], 14);
}
}
}
#endif
/// <summary>
/// Can solve Ax=b using QR factorization with a square A matrix.
/// </summary>
@ -944,94 +735,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 13);
}
#if ! MKL
/// <summary>
/// Can solve Ax=b using QR factorization with a square A matrix
/// using a work array.
/// </summary>
[Test]
public void CanSolveUsingQRSquareMatrixUsingWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new Complex[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var b = new[] {new Complex(1.0, 0), 2.0, 3.0, 4.0, 5.0, 6.0};
var x = new Complex[matrix.ColumnCount*2];
var work = new Complex[matrix.RowCount*matrix.RowCount];
Control.LinearAlgebraProvider.QRSolve(a, matrix.RowCount, matrix.ColumnCount, b, 2, x, work);
NotModified(3, 3, a, matrix);
var mx = new DenseMatrix(matrix.ColumnCount, 2, x);
var mb = matrix*mx;
AssertHelpers.AlmostEqualRelative(mb[0, 0], b[0], 13);
AssertHelpers.AlmostEqualRelative(mb[1, 0], b[1], 13);
AssertHelpers.AlmostEqualRelative(mb[2, 0], b[2], 13);
AssertHelpers.AlmostEqualRelative(mb[0, 1], b[3], 13);
AssertHelpers.AlmostEqualRelative(mb[1, 1], b[4], 13);
AssertHelpers.AlmostEqualRelative(mb[2, 1], b[5], 13);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a tall A matrix
/// using a work array.
/// </summary>
[Test]
public void CanSolveUsingQRTallMatrixUsingWorkArray()
{
var matrix = _matrices["Tall3x2"];
var a = new Complex[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var b = new[] {new Complex(1.0, 0), 2.0, 3.0, 4.0, 5.0, 6.0};
var x = new Complex[matrix.ColumnCount*2];
var work = new Complex[matrix.RowCount*matrix.RowCount];
Control.LinearAlgebraProvider.QRSolve(a, matrix.RowCount, matrix.ColumnCount, b, 2, x, work);
NotModified(3, 2, a, matrix);
var mb = new DenseMatrix(matrix.RowCount, 2, b);
var test = (matrix.Transpose()*matrix).Inverse()*matrix.Transpose()*mb;
AssertHelpers.AlmostEqualRelative(test[0, 0], x[0], 13);
AssertHelpers.AlmostEqualRelative(test[1, 0], x[1], 13);
AssertHelpers.AlmostEqualRelative(test[0, 1], x[2], 13);
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 13);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a square A matrix
/// using a factored A matrix.
/// </summary>
[Test]
public void CanSolveUsingQRSquareMatrixOnFactoredMatrix()
{
var matrix = _matrices["Square3x3"];
var a = new Complex[matrix.RowCount*matrix.RowCount];
Array.Copy(matrix.Values, a, a.Length);
var tau = new Complex[matrix.ColumnCount];
var q = new Complex[matrix.ColumnCount*matrix.ColumnCount];
Control.LinearAlgebraProvider.QRFactor(a, matrix.RowCount, matrix.ColumnCount, q, tau);
var b = new[] {new Complex(1.0, 0), 2.0, 3.0, 4.0, 5.0, 6.0};
var x = new Complex[matrix.ColumnCount*2];
Control.LinearAlgebraProvider.QRSolveFactored(q, a, matrix.RowCount, matrix.ColumnCount, tau, b, 2, x);
var mx = new DenseMatrix(matrix.ColumnCount, 2, x);
var mb = matrix*mx;
AssertHelpers.AlmostEqualRelative(mb[0, 0], b[0], 13);
AssertHelpers.AlmostEqualRelative(mb[1, 0], b[1], 13);
AssertHelpers.AlmostEqualRelative(mb[2, 0], b[2], 13);
AssertHelpers.AlmostEqualRelative(mb[0, 1], b[3], 13);
AssertHelpers.AlmostEqualRelative(mb[1, 1], b[4], 13);
AssertHelpers.AlmostEqualRelative(mb[2, 1], b[5], 13);
}
#endif
/// <summary>
/// Can solve Ax=b using QR factorization with a tall A matrix
/// using a factored A matrix.
@ -1060,68 +763,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 13);
}
#if ! MKL
/// <summary>
/// Can solve Ax=b using QR factorization with a square A matrix
/// using a factored A matrix with a work array.
/// </summary>
[Test]
public void CanSolveUsingQRSquareMatrixOnFactoredMatrixWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new Complex[matrix.RowCount*matrix.RowCount];
Array.Copy(matrix.Values, a, a.Length);
var tau = new Complex[matrix.ColumnCount];
var q = new Complex[matrix.ColumnCount*matrix.ColumnCount];
var work = new Complex[2048];
Control.LinearAlgebraProvider.QRFactor(a, matrix.RowCount, matrix.ColumnCount, q, tau, work);
var b = new[] {new Complex(1.0, 0), 2.0, 3.0, 4.0, 5.0, 6.0};
var x = new Complex[matrix.ColumnCount*2];
Control.LinearAlgebraProvider.QRSolveFactored(q, a, matrix.RowCount, matrix.ColumnCount, tau, b, 2, x, work);
var mx = new DenseMatrix(matrix.ColumnCount, 2, x);
var mb = matrix*mx;
AssertHelpers.AlmostEqualRelative(mb[0, 0], b[0], 13);
AssertHelpers.AlmostEqualRelative(mb[1, 0], b[1], 13);
AssertHelpers.AlmostEqualRelative(mb[2, 0], b[2], 13);
AssertHelpers.AlmostEqualRelative(mb[0, 1], b[3], 13);
AssertHelpers.AlmostEqualRelative(mb[1, 1], b[4], 13);
AssertHelpers.AlmostEqualRelative(mb[2, 1], b[5], 13);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a tall A matrix
/// using a factored A matrix with a work array.
/// </summary>
[Test]
public void CanSolveUsingQRTallMatrixOnFactoredMatrixWithWorkArray()
{
var matrix = _matrices["Tall3x2"];
var a = new Complex[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var tau = new Complex[matrix.ColumnCount];
var q = new Complex[matrix.RowCount*matrix.RowCount];
var work = new Complex[2048];
Control.LinearAlgebraProvider.QRFactor(a, matrix.RowCount, matrix.ColumnCount, q, tau, work);
var b = new[] {new Complex(1.0, 0), 2.0, 3.0, 4.0, 5.0, 6.0};
var x = new Complex[matrix.ColumnCount*2];
Control.LinearAlgebraProvider.QRSolveFactored(q, a, matrix.RowCount, matrix.ColumnCount, tau, b, 2, x, work);
var mb = new DenseMatrix(matrix.RowCount, 2, b);
var test = (matrix.Transpose()*matrix).Inverse()*matrix.Transpose()*mb;
AssertHelpers.AlmostEqualRelative(test[0, 0], x[0], 13);
AssertHelpers.AlmostEqualRelative(test[1, 0], x[1], 13);
AssertHelpers.AlmostEqualRelative(test[0, 1], x[2], 13);
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 13);
}
#endif
/// <summary>
/// Can solve Ax=b using thin QR factorization with a square A matrix.
/// </summary>
@ -1174,64 +815,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 13);
}
#if ! MKL
/// <summary>
/// Can solve Ax=b using thin QR factorization with a square A matrix
/// using a work array.
/// </summary>
[Test]
public void CanSolveUsingThinQRSquareMatrixUsingWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new Complex[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var b = new[] {new Complex(1.0, 0), 2.0, 3.0, 4.0, 5.0, 6.0};
var x = new Complex[matrix.ColumnCount*2];
var work = new Complex[matrix.RowCount*matrix.ColumnCount];
Control.LinearAlgebraProvider.QRSolve(a, matrix.RowCount, matrix.ColumnCount, b, 2, x, work, QRMethod.Thin);
NotModified(3, 3, a, matrix);
var mx = new DenseMatrix(matrix.ColumnCount, 2, x);
var mb = matrix*mx;
AssertHelpers.AlmostEqualRelative(mb[0, 0], b[0], 13);
AssertHelpers.AlmostEqualRelative(mb[1, 0], b[1], 13);
AssertHelpers.AlmostEqualRelative(mb[2, 0], b[2], 13);
AssertHelpers.AlmostEqualRelative(mb[0, 1], b[3], 13);
AssertHelpers.AlmostEqualRelative(mb[1, 1], b[4], 13);
AssertHelpers.AlmostEqualRelative(mb[2, 1], b[5], 13);
}
/// <summary>
/// Can solve Ax=b using thin QR factorization with a tall A matrix
/// using a work array.
/// </summary>
[Test]
public void CanSolveUsingThinQRTallMatrixUsingWorkArray()
{
var matrix = _matrices["Tall3x2"];
var a = new Complex[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var b = new[] {new Complex(1.0, 0), 2.0, 3.0, 4.0, 5.0, 6.0};
var x = new Complex[matrix.ColumnCount*2];
var work = new Complex[matrix.RowCount*matrix.ColumnCount];
Control.LinearAlgebraProvider.QRSolve(a, matrix.RowCount, matrix.ColumnCount, b, 2, x, work, QRMethod.Thin);
NotModified(3, 2, a, matrix);
var mb = new DenseMatrix(matrix.RowCount, 2, b);
var test = (matrix.Transpose()*matrix).Inverse()*matrix.Transpose()*mb;
AssertHelpers.AlmostEqualRelative(test[0, 0], x[0], 13);
AssertHelpers.AlmostEqualRelative(test[1, 0], x[1], 13);
AssertHelpers.AlmostEqualRelative(test[0, 1], x[2], 13);
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 13);
}
#endif
/// <summary>
/// Can solve Ax=b using thin QR factorization with a square A matrix
/// using a factored A matrix.
@ -1290,66 +873,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 13);
}
/// <summary>
/// Can solve Ax=b using thin QR factorization with a square A matrix
/// using a factored A matrix with a work array.
/// </summary>
[Test]
public void CanSolveUsingThinQRSquareMatrixOnFactoredMatrixWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new Complex[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var tau = new Complex[matrix.ColumnCount];
var r = new Complex[matrix.ColumnCount*matrix.ColumnCount];
var work = new Complex[2048];
Control.LinearAlgebraProvider.ThinQRFactor(a, matrix.RowCount, matrix.ColumnCount, r, tau, work);
var b = new[] {new Complex(1.0, 0), 2.0, 3.0, 4.0, 5.0, 6.0};
var x = new Complex[matrix.ColumnCount*2];
Control.LinearAlgebraProvider.QRSolveFactored(a, r, matrix.RowCount, matrix.ColumnCount, tau, b, 2, x, work, QRMethod.Thin);
var mx = new DenseMatrix(matrix.ColumnCount, 2, x);
var mb = matrix*mx;
AssertHelpers.AlmostEqualRelative(mb[0, 0], b[0], 13);
AssertHelpers.AlmostEqualRelative(mb[1, 0], b[1], 13);
AssertHelpers.AlmostEqualRelative(mb[2, 0], b[2], 13);
AssertHelpers.AlmostEqualRelative(mb[0, 1], b[3], 13);
AssertHelpers.AlmostEqualRelative(mb[1, 1], b[4], 13);
AssertHelpers.AlmostEqualRelative(mb[2, 1], b[5], 13);
}
/// <summary>
/// Can solve Ax=b using thin QR factorization with a tall A matrix
/// using a factored A matrix with a work array.
/// </summary>
[Test]
public void CanSolveUsingThinQRTallMatrixOnFactoredMatrixWithWorkArray()
{
var matrix = _matrices["Tall3x2"];
var a = new Complex[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var tau = new Complex[matrix.ColumnCount];
var r = new Complex[matrix.ColumnCount*matrix.ColumnCount];
var work = new Complex[2048];
Control.LinearAlgebraProvider.ThinQRFactor(a, matrix.RowCount, matrix.ColumnCount, r, tau, work);
var b = new[] {new Complex(1.0, 0), 2.0, 3.0, 4.0, 5.0, 6.0};
var x = new Complex[matrix.ColumnCount*2];
Control.LinearAlgebraProvider.QRSolveFactored(a, r, matrix.RowCount, matrix.ColumnCount, tau, b, 2, x, work, QRMethod.Thin);
var mb = new DenseMatrix(matrix.RowCount, 2, b);
var test = (matrix.Transpose()*matrix).Inverse()*matrix.Transpose()*mb;
AssertHelpers.AlmostEqualRelative(test[0, 0], x[0], 13);
AssertHelpers.AlmostEqualRelative(test[1, 0], x[1], 13);
AssertHelpers.AlmostEqualRelative(test[0, 1], x[2], 13);
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 13);
}
/// <summary>
/// Can compute the SVD factorization of a square matrix.
/// </summary>
@ -1455,129 +978,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex
AssertHelpers.AlmostEqualRelative(matrix[1, 2], result[1, 2], 14);
}
/// <summary>
/// Can compute the SVD factorization of a square matrix using
/// a work array.
/// </summary>
[Test]
public void CanComputeSVDFactorizationOfSquareMatrixWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new Complex[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var s = new Complex[matrix.RowCount];
var u = new Complex[matrix.RowCount*matrix.RowCount];
var vt = new Complex[matrix.ColumnCount*matrix.ColumnCount];
#if CUDA
Complex[] work = null;
#else
var work = new Complex[100];
#endif
Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work);
var w = new DenseMatrix(matrix.RowCount, matrix.ColumnCount);
for (var index = 0; index < s.Length; index++)
{
w[index, index] = s[index];
}
var mU = new DenseMatrix(matrix.RowCount, matrix.RowCount, u);
var mV = new DenseMatrix(matrix.ColumnCount, matrix.ColumnCount, vt);
var result = mU*w*mV;
AssertHelpers.AlmostEqualRelative(matrix[0, 0], result[0, 0], 13);
AssertHelpers.AlmostEqualRelative(matrix[1, 0], result[1, 0], 13);
AssertHelpers.AlmostEqualRelative(matrix[2, 0], result[2, 0], 13);
AssertHelpers.AlmostEqualRelative(matrix[0, 1], result[0, 1], 13);
AssertHelpers.AlmostEqualRelative(matrix[1, 1], result[1, 1], 13);
AssertHelpers.AlmostEqualRelative(matrix[2, 1], result[2, 1], 13);
AssertHelpers.AlmostEqualRelative(matrix[0, 2], result[0, 2], 13);
AssertHelpers.AlmostEqualRelative(matrix[1, 2], result[1, 2], 13);
AssertHelpers.AlmostEqualRelative(matrix[2, 2], result[2, 2], 13);
}
/// <summary>
/// Can compute the SVD factorization of a tall matrix using
/// a work array.
/// </summary>
[Test]
public void CanComputeSVDFactorizationOfTallMatrixWithWorkArray()
{
var matrix = _matrices["Tall3x2"];
var a = new Complex[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var s = new Complex[matrix.ColumnCount];
var u = new Complex[matrix.RowCount*matrix.RowCount];
var vt = new Complex[matrix.ColumnCount*matrix.ColumnCount];
#if CUDA
Complex[] work = null;
#else
var work = new Complex[100];
#endif
Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work);
var w = new DenseMatrix(matrix.RowCount, matrix.ColumnCount);
for (var index = 0; index < s.Length; index++)
{
w[index, index] = s[index];
}
var mU = new DenseMatrix(matrix.RowCount, matrix.RowCount, u);
var mV = new DenseMatrix(matrix.ColumnCount, matrix.ColumnCount, vt);
var result = mU*w*mV;
AssertHelpers.AlmostEqualRelative(matrix[0, 0], result[0, 0], 14);
AssertHelpers.AlmostEqualRelative(matrix[1, 0], result[1, 0], 14);
AssertHelpers.AlmostEqualRelative(matrix[2, 0], result[2, 0], 14);
AssertHelpers.AlmostEqualRelative(matrix[0, 1], result[0, 1], 14);
AssertHelpers.AlmostEqualRelative(matrix[1, 1], result[1, 1], 14);
AssertHelpers.AlmostEqualRelative(matrix[2, 1], result[2, 1], 14);
}
/// <summary>
/// Can compute the SVD factorization of a wide matrix using
/// a work array.
/// </summary>
[Test]
public void CanComputeSVDFactorizationOfWideMatrixWithWorkArray()
{
var matrix = _matrices["Wide2x3"];
var a = new Complex[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var s = new Complex[matrix.RowCount];
var u = new Complex[matrix.RowCount*matrix.RowCount];
var vt = new Complex[matrix.ColumnCount*matrix.ColumnCount];
#if CUDA
Complex[] work = null;
#else
var work = new Complex[100];
#endif
Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work);
var w = new DenseMatrix(matrix.RowCount, matrix.ColumnCount);
for (var index = 0; index < s.Length; index++)
{
w[index, index] = s[index];
}
var mU = new DenseMatrix(matrix.RowCount, matrix.RowCount, u);
var mV = new DenseMatrix(matrix.ColumnCount, matrix.ColumnCount, vt);
var result = mU*w*mV;
AssertHelpers.AlmostEqualRelative(matrix[0, 0], result[0, 0], 14);
AssertHelpers.AlmostEqualRelative(matrix[1, 0], result[1, 0], 14);
AssertHelpers.AlmostEqualRelative(matrix[0, 1], result[0, 1], 14);
AssertHelpers.AlmostEqualRelative(matrix[1, 1], result[1, 1], 14);
AssertHelpers.AlmostEqualRelative(matrix[0, 2], result[0, 2], 14);
AssertHelpers.AlmostEqualRelative(matrix[1, 2], result[1, 2], 14);
}
/// <summary>
/// Can solve Ax=b using SVD factorization with a square A matrix.
/// </summary>

558
src/UnitTests/LinearAlgebraProviderTests/Complex32/LinearAlgebraProviderTests.cs

@ -4,7 +4,7 @@
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 Math.NET
// Copyright (c) 2009-2015 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
@ -437,68 +437,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex32
AssertHelpers.AlmostEqualRelative(a[8], -0.113636363636364f, 5);
}
/// <summary>
/// Can compute the inverse of a matrix using LU factorization
/// with a work array.
/// </summary>
[Test]
public void CanComputeLuInverseWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new Complex32[matrix.RowCount*matrix.RowCount];
Array.Copy(matrix.Values, a, a.Length);
#if CUDA
Complex32[] work = null;
#else
var work = new Complex32[matrix.RowCount];
#endif
Control.LinearAlgebraProvider.LUInverse(a, matrix.RowCount, work);
AssertHelpers.AlmostEqualRelative(a[0], -0.454545454545454f, 5);
AssertHelpers.AlmostEqualRelative(a[1], -0.909090909090908f, 5);
AssertHelpers.AlmostEqualRelative(a[2], 0.454545454545454f, 5);
AssertHelpers.AlmostEqualRelative(a[3], -0.340909090909090f, 5);
AssertHelpers.AlmostEqualRelative(a[4], -2.045454545454543f, 5);
AssertHelpers.AlmostEqualRelative(a[5], 1.477272727272726f, 5);
AssertHelpers.AlmostEqualRelative(a[6], -0.113636363636364f, 5);
AssertHelpers.AlmostEqualRelative(a[7], 0.227272727272727f, 5);
AssertHelpers.AlmostEqualRelative(a[8], -0.113636363636364f, 5);
}
/// <summary>
/// Can compute the inverse of a matrix using LU factorization
/// using a previously factored matrix with a work array.
/// </summary>
[Test]
public void CanComputeLuInverseOnFactoredMatrixWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new Complex32[matrix.RowCount*matrix.RowCount];
Array.Copy(matrix.Values, a, a.Length);
var ipiv = new int[matrix.RowCount];
Control.LinearAlgebraProvider.LUFactor(a, matrix.RowCount, ipiv);
#if CUDA
Complex32[] work = null;
#else
var work = new Complex32[matrix.RowCount];
#endif
Control.LinearAlgebraProvider.LUInverseFactored(a, matrix.RowCount, ipiv, work);
AssertHelpers.AlmostEqualRelative(a[0], -0.454545454545454f, 5);
AssertHelpers.AlmostEqualRelative(a[1], -0.909090909090908f, 5);
AssertHelpers.AlmostEqualRelative(a[2], 0.454545454545454f, 5);
AssertHelpers.AlmostEqualRelative(a[3], -0.340909090909090f, 5);
AssertHelpers.AlmostEqualRelative(a[4], -2.045454545454543f, 5);
AssertHelpers.AlmostEqualRelative(a[5], 1.477272727272726f, 5);
AssertHelpers.AlmostEqualRelative(a[6], -0.113636363636364f, 5);
AssertHelpers.AlmostEqualRelative(a[7], 0.227272727272727f, 5);
AssertHelpers.AlmostEqualRelative(a[8], -0.113636363636364f, 5);
}
/// <summary>
/// Can solve Ax=b using LU factorization.
/// </summary>
@ -696,90 +634,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex32
}
}
/// <summary>
/// Can compute QR factorization of a square matrix using a work array.
/// </summary>
[Test]
public void CanComputeQRFactorSquareMatrixWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var r = new Complex32[matrix.RowCount*matrix.ColumnCount];
var tau = new Complex32[3];
var q = new Complex32[matrix.RowCount*matrix.RowCount];
Array.Copy(matrix.Values, r, r.Length);
var work = new Complex32[matrix.ColumnCount*Control.BlockSize];
Control.LinearAlgebraProvider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau, work);
var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q);
var mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r).UpperTriangle();
var a = mq*mr;
for (var row = 0; row < matrix.RowCount; row++)
{
for (var col = 0; col < matrix.ColumnCount; col++)
{
AssertHelpers.AlmostEqualRelative(matrix[row, col], a[row, col], 5);
}
}
}
/// <summary>
/// Can compute QR factorization of a tall matrix using a work matrix.
/// </summary>
[Test]
public void CanComputeQRFactorTallMatrixWithWorkArray()
{
var matrix = _matrices["Tall3x2"];
var r = new Complex32[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, r, r.Length);
var tau = new Complex32[3];
var q = new Complex32[matrix.RowCount*matrix.RowCount];
var work = new Complex32[matrix.ColumnCount*Control.BlockSize];
Control.LinearAlgebraProvider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau, work);
var mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r).UpperTriangle();
var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q);
var a = mq*mr;
for (var row = 0; row < matrix.RowCount; row++)
{
for (var col = 0; col < matrix.ColumnCount; col++)
{
AssertHelpers.AlmostEqualRelative(matrix[row, col], a[row, col], 5);
}
}
}
/// <summary>
/// Can compute QR factorization of a wide matrix using a work matrix.
/// </summary>
[Test]
public void CanComputeQRFactorWideMatrixWithWorkArray()
{
var matrix = _matrices["Wide2x3"];
var r = new Complex32[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, r, r.Length);
var tau = new Complex32[3];
var q = new Complex32[matrix.RowCount*matrix.RowCount];
var work = new Complex32[matrix.ColumnCount*Control.BlockSize];
Control.LinearAlgebraProvider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau, work);
var mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r).UpperTriangle();
var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q);
var a = mq*mr;
for (var row = 0; row < matrix.RowCount; row++)
{
for (var col = 0; col < matrix.ColumnCount; col++)
{
AssertHelpers.AlmostEqualRelative(matrix[row, col], a[row, col], 5);
}
}
}
/// <summary>
/// Can compute thin QR factorization of a square matrix.
/// </summary>
@ -834,61 +688,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex32
}
}
/// <summary>
/// Can compute thin QR factorization of a square matrix using a work array.
/// </summary>
[Test]
public void CanComputeThinQRFactorSquareMatrixWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var r = new Complex32[matrix.ColumnCount*matrix.ColumnCount];
var tau = new Complex32[3];
var q = new Complex32[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, q, q.Length);
var work = new Complex32[matrix.ColumnCount*Control.BlockSize];
Control.LinearAlgebraProvider.ThinQRFactor(q, matrix.RowCount, matrix.ColumnCount, r, tau, work);
var mq = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, q);
var mr = new DenseMatrix(matrix.ColumnCount, matrix.ColumnCount, r);
var a = mq*mr;
for (var row = 0; row < matrix.RowCount; row++)
{
for (var col = 0; col < matrix.ColumnCount; col++)
{
AssertHelpers.AlmostEqualRelative(matrix[row, col], a[row, col], 5);
}
}
}
/// <summary>
/// Can compute thin QR factorization of a tall matrix using a work matrix.
/// </summary>
[Test]
public void CanComputeThinQRFactorTallMatrixWithWorkArray()
{
var matrix = _matrices["Tall3x2"];
var r = new Complex32[matrix.ColumnCount*matrix.ColumnCount];
var tau = new Complex32[3];
var q = new Complex32[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, q, q.Length);
var work = new Complex32[matrix.ColumnCount*Control.BlockSize];
Control.LinearAlgebraProvider.ThinQRFactor(q, matrix.RowCount, matrix.ColumnCount, r, tau, work);
var mq = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, q);
var mr = new DenseMatrix(matrix.ColumnCount, matrix.ColumnCount, r);
var a = mq*mr;
for (var row = 0; row < matrix.RowCount; row++)
{
for (var col = 0; col < matrix.ColumnCount; col++)
{
AssertHelpers.AlmostEqualRelative(matrix[row, col], a[row, col], 5);
}
}
}
/// <summary>
/// Can solve Ax=b using QR factorization with a square A matrix.
/// </summary>
@ -941,62 +740,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex32
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 5);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a square A matrix
/// using a work array.
/// </summary>
[Test]
public void CanSolveUsingQRSquareMatrixUsingWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new Complex32[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var b = new[] {new Complex32(1.0f, 0.0f), 2.0f, 3.0f, 4.0f, 5.0f, 6.0f};
var x = new Complex32[matrix.ColumnCount*2];
var work = new Complex32[matrix.RowCount*matrix.RowCount];
Control.LinearAlgebraProvider.QRSolve(a, matrix.RowCount, matrix.ColumnCount, b, 2, x, work);
NotModified(3, 3, a, matrix);
var mx = new DenseMatrix(matrix.ColumnCount, 2, x);
var mb = matrix*mx;
AssertHelpers.AlmostEqualRelative(mb[0, 0], b[0], 5);
AssertHelpers.AlmostEqualRelative(mb[1, 0], b[1], 5);
AssertHelpers.AlmostEqualRelative(mb[2, 0], b[2], 5);
AssertHelpers.AlmostEqualRelative(mb[0, 1], b[3], 5);
AssertHelpers.AlmostEqualRelative(mb[1, 1], b[4], 4);
AssertHelpers.AlmostEqualRelative(mb[2, 1], b[5], 4);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a tall A matrix
/// using a work array.
/// </summary>
[Test]
public void CanSolveUsingQRTallMatrixUsingWorkArray()
{
var matrix = _matrices["Tall3x2"];
var a = new Complex32[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var b = new[] {new Complex32(1.0f, 0.0f), 2.0f, 3.0f, 4.0f, 5.0f, 6.0f};
var x = new Complex32[matrix.ColumnCount*2];
var work = new Complex32[matrix.RowCount*matrix.RowCount];
Control.LinearAlgebraProvider.QRSolve(a, matrix.RowCount, matrix.ColumnCount, b, 2, x, work);
NotModified(3, 2, a, matrix);
var mb = new DenseMatrix(matrix.RowCount, 2, b);
var test = (matrix.Transpose()*matrix).Inverse()*matrix.Transpose()*mb;
AssertHelpers.AlmostEqualRelative(test[0, 0], x[0], 5);
AssertHelpers.AlmostEqualRelative(test[1, 0], x[1], 5);
AssertHelpers.AlmostEqualRelative(test[0, 1], x[2], 5);
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 5);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a square A matrix
/// using a factored A matrix.
@ -1055,66 +798,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex32
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 5);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a square A matrix
/// using a factored A matrix with a work array.
/// </summary>
[Test]
public void CanSolveUsingQRSquareMatrixOnFactoredMatrixWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new Complex32[matrix.RowCount*matrix.RowCount];
Array.Copy(matrix.Values, a, a.Length);
var tau = new Complex32[matrix.ColumnCount];
var q = new Complex32[matrix.ColumnCount*matrix.ColumnCount];
var work = new Complex32[2048];
Control.LinearAlgebraProvider.QRFactor(a, matrix.RowCount, matrix.ColumnCount, q, tau, work);
var b = new[] {new Complex32(1.0f, 0.0f), 2.0f, 3.0f, 4.0f, 5.0f, 6.0f};
var x = new Complex32[matrix.ColumnCount*2];
Control.LinearAlgebraProvider.QRSolveFactored(q, a, matrix.RowCount, matrix.ColumnCount, tau, b, 2, x, work);
var mx = new DenseMatrix(matrix.ColumnCount, 2, x);
var mb = matrix*mx;
AssertHelpers.AlmostEqualRelative(mb[0, 0], b[0], 5);
AssertHelpers.AlmostEqualRelative(mb[1, 0], b[1], 5);
AssertHelpers.AlmostEqualRelative(mb[2, 0], b[2], 5);
AssertHelpers.AlmostEqualRelative(mb[0, 1], b[3], 5);
AssertHelpers.AlmostEqualRelative(mb[1, 1], b[4], 4);
AssertHelpers.AlmostEqualRelative(mb[2, 1], b[5], 4);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a tall A matrix
/// using a factored A matrix with a work array.
/// </summary>
[Test]
public void CanSolveUsingQRTallMatrixOnFactoredMatrixWithWorkArray()
{
var matrix = _matrices["Tall3x2"];
var a = new Complex32[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var tau = new Complex32[matrix.ColumnCount];
var q = new Complex32[matrix.RowCount*matrix.RowCount];
var work = new Complex32[2048];
Control.LinearAlgebraProvider.QRFactor(a, matrix.RowCount, matrix.ColumnCount, q, tau, work);
var b = new[] {new Complex32(1.0f, 0.0f), 2.0f, 3.0f, 4.0f, 5.0f, 6.0f};
var x = new Complex32[matrix.ColumnCount*2];
Control.LinearAlgebraProvider.QRSolveFactored(q, a, matrix.RowCount, matrix.ColumnCount, tau, b, 2, x, work);
var mb = new DenseMatrix(matrix.RowCount, 2, b);
var test = (matrix.Transpose()*matrix).Inverse()*matrix.Transpose()*mb;
AssertHelpers.AlmostEqualRelative(test[0, 0], x[0], 5);
AssertHelpers.AlmostEqualRelative(test[1, 0], x[1], 5);
AssertHelpers.AlmostEqualRelative(test[0, 1], x[2], 5);
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 5);
}
/// <summary>
/// Can solve Ax=b using thin QR factorization with a square A matrix.
/// </summary>
@ -1167,62 +850,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex32
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 5);
}
/// <summary>
/// Can solve Ax=b using thin QR factorization with a square A matrix
/// using a work array.
/// </summary>
[Test]
public void CanSolveUsingThinQRSquareMatrixUsingWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new Complex32[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var b = new[] {new Complex32(1.0f, 0.0f), 2.0f, 3.0f, 4.0f, 5.0f, 6.0f};
var x = new Complex32[matrix.ColumnCount*2];
var work = new Complex32[matrix.RowCount*matrix.ColumnCount];
Control.LinearAlgebraProvider.QRSolve(a, matrix.RowCount, matrix.ColumnCount, b, 2, x, work, QRMethod.Thin);
NotModified(3, 3, a, matrix);
var mx = new DenseMatrix(matrix.ColumnCount, 2, x);
var mb = matrix*mx;
AssertHelpers.AlmostEqualRelative(mb[0, 0], b[0], 5);
AssertHelpers.AlmostEqualRelative(mb[1, 0], b[1], 5);
AssertHelpers.AlmostEqualRelative(mb[2, 0], b[2], 5);
AssertHelpers.AlmostEqualRelative(mb[0, 1], b[3], 5);
AssertHelpers.AlmostEqualRelative(mb[1, 1], b[4], 4);
AssertHelpers.AlmostEqualRelative(mb[2, 1], b[5], 4);
}
/// <summary>
/// Can solve Ax=b using thin QR factorization with a tall A matrix
/// using a work array.
/// </summary>
[Test]
public void CanSolveUsingThinQRTallMatrixUsingWorkArray()
{
var matrix = _matrices["Tall3x2"];
var a = new Complex32[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var b = new[] {new Complex32(1.0f, 0.0f), 2.0f, 3.0f, 4.0f, 5.0f, 6.0f};
var x = new Complex32[matrix.ColumnCount*2];
var work = new Complex32[matrix.RowCount*matrix.ColumnCount];
Control.LinearAlgebraProvider.QRSolve(a, matrix.RowCount, matrix.ColumnCount, b, 2, x, work, QRMethod.Thin);
NotModified(3, 2, a, matrix);
var mb = new DenseMatrix(matrix.RowCount, 2, b);
var test = (matrix.Transpose()*matrix).Inverse()*matrix.Transpose()*mb;
AssertHelpers.AlmostEqualRelative(test[0, 0], x[0], 5);
AssertHelpers.AlmostEqualRelative(test[1, 0], x[1], 5);
AssertHelpers.AlmostEqualRelative(test[0, 1], x[2], 5);
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 5);
}
/// <summary>
/// Can solve Ax=b using thin QR factorization with a square A matrix
/// using a factored A matrix.
@ -1281,66 +908,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex32
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 5);
}
/// <summary>
/// Can solve Ax=b using thin QR factorization with a square A matrix
/// using a factored A matrix with a work array.
/// </summary>
[Test]
public void CanSolveUsingThinQRSquareMatrixOnFactoredMatrixWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new Complex32[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var tau = new Complex32[matrix.ColumnCount];
var r = new Complex32[matrix.ColumnCount*matrix.ColumnCount];
var work = new Complex32[2048];
Control.LinearAlgebraProvider.ThinQRFactor(a, matrix.RowCount, matrix.ColumnCount, r, tau, work);
var b = new[] {new Complex32(1.0f, 0.0f), 2.0f, 3.0f, 4.0f, 5.0f, 6.0f};
var x = new Complex32[matrix.ColumnCount*2];
Control.LinearAlgebraProvider.QRSolveFactored(a, r, matrix.RowCount, matrix.ColumnCount, tau, b, 2, x, work, QRMethod.Thin);
var mx = new DenseMatrix(matrix.ColumnCount, 2, x);
var mb = matrix*mx;
AssertHelpers.AlmostEqualRelative(mb[0, 0], b[0], 5);
AssertHelpers.AlmostEqualRelative(mb[1, 0], b[1], 5);
AssertHelpers.AlmostEqualRelative(mb[2, 0], b[2], 5);
AssertHelpers.AlmostEqualRelative(mb[0, 1], b[3], 5);
AssertHelpers.AlmostEqualRelative(mb[1, 1], b[4], 4);
AssertHelpers.AlmostEqualRelative(mb[2, 1], b[5], 4);
}
/// <summary>
/// Can solve Ax=b using thin QR factorization with a tall A matrix
/// using a factored A matrix with a work array.
/// </summary>
[Test]
public void CanSolveUsingThinQRTallMatrixOnFactoredMatrixWithWorkArray()
{
var matrix = _matrices["Tall3x2"];
var a = new Complex32[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var tau = new Complex32[matrix.ColumnCount];
var r = new Complex32[matrix.ColumnCount*matrix.ColumnCount];
var work = new Complex32[2048];
Control.LinearAlgebraProvider.ThinQRFactor(a, matrix.RowCount, matrix.ColumnCount, r, tau, work);
var b = new[] {new Complex32(1.0f, 0.0f), 2.0f, 3.0f, 4.0f, 5.0f, 6.0f};
var x = new Complex32[matrix.ColumnCount*2];
Control.LinearAlgebraProvider.QRSolveFactored(a, r, matrix.RowCount, matrix.ColumnCount, tau, b, 2, x, work, QRMethod.Thin);
var mb = new DenseMatrix(matrix.RowCount, 2, b);
var test = (matrix.Transpose()*matrix).Inverse()*matrix.Transpose()*mb;
AssertHelpers.AlmostEqualRelative(test[0, 0], x[0], 5);
AssertHelpers.AlmostEqualRelative(test[1, 0], x[1], 5);
AssertHelpers.AlmostEqualRelative(test[0, 1], x[2], 5);
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 5);
}
/// <summary>
/// Can compute the SVD factorization of a square matrix.
/// </summary>
@ -1446,129 +1013,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex32
AssertHelpers.AlmostEqualRelative(matrix[1, 2], result[1, 2], 4);
}
/// <summary>
/// Can compute the SVD factorization of a square matrix using
/// a work array.
/// </summary>
[Test]
public void CanComputeSVDFactorizationOfSquareMatrixWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new Complex32[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var s = new Complex32[matrix.RowCount];
var u = new Complex32[matrix.RowCount*matrix.RowCount];
var vt = new Complex32[matrix.ColumnCount*matrix.ColumnCount];
#if CUDA
Complex32[] work = null;
#else
var work = new Complex32[100];
#endif
Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work);
var w = new DenseMatrix(matrix.RowCount, matrix.ColumnCount);
for (var index = 0; index < s.Length; index++)
{
w[index, index] = s[index];
}
var mU = new DenseMatrix(matrix.RowCount, matrix.RowCount, u);
var mV = new DenseMatrix(matrix.ColumnCount, matrix.ColumnCount, vt);
var result = mU*w*mV;
AssertHelpers.AlmostEqualRelative(matrix[0, 0], result[0, 0], 5);
AssertHelpers.AlmostEqualRelative(matrix[1, 0], result[1, 0], 5);
AssertHelpers.AlmostEqualRelative(matrix[2, 0], result[2, 0], 5);
AssertHelpers.AlmostEqualRelative(matrix[0, 1], result[0, 1], 5);
AssertHelpers.AlmostEqualRelative(matrix[1, 1], result[1, 1], 5);
AssertHelpers.AlmostEqualRelative(matrix[2, 1], result[2, 1], 5);
AssertHelpers.AlmostEqualRelative(matrix[0, 2], result[0, 2], 5);
AssertHelpers.AlmostEqualRelative(matrix[1, 2], result[1, 2], 5);
AssertHelpers.AlmostEqualRelative(matrix[2, 2], result[2, 2], 5);
}
/// <summary>
/// Can compute the SVD factorization of a tall matrix using
/// a work array.
/// </summary>
[Test]
public void CanComputeSVDFactorizationOfTallMatrixWithWorkArray()
{
var matrix = _matrices["Tall3x2"];
var a = new Complex32[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var s = new Complex32[matrix.ColumnCount];
var u = new Complex32[matrix.RowCount*matrix.RowCount];
var vt = new Complex32[matrix.ColumnCount*matrix.ColumnCount];
#if CUDA
Complex32[] work = null;
#else
var work = new Complex32[100];
#endif
Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work);
var w = new DenseMatrix(matrix.RowCount, matrix.ColumnCount);
for (var index = 0; index < s.Length; index++)
{
w[index, index] = s[index];
}
var mU = new DenseMatrix(matrix.RowCount, matrix.RowCount, u);
var mV = new DenseMatrix(matrix.ColumnCount, matrix.ColumnCount, vt);
var result = mU*w*mV;
AssertHelpers.AlmostEqualRelative(matrix[0, 0], result[0, 0], 5);
AssertHelpers.AlmostEqualRelative(matrix[1, 0], result[1, 0], 5);
AssertHelpers.AlmostEqualRelative(matrix[2, 0], result[2, 0], 5);
AssertHelpers.AlmostEqualRelative(matrix[0, 1], result[0, 1], 5);
AssertHelpers.AlmostEqualRelative(matrix[1, 1], result[1, 1], 5);
AssertHelpers.AlmostEqualRelative(matrix[2, 1], result[2, 1], 4);
}
/// <summary>
/// Can compute the SVD factorization of a wide matrix using
/// a work array.
/// </summary>
[Test]
public void CanComputeSVDFactorizationOfWideMatrixWithWorkArray()
{
var matrix = _matrices["Wide2x3"];
var a = new Complex32[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var s = new Complex32[matrix.RowCount];
var u = new Complex32[matrix.RowCount*matrix.RowCount];
var vt = new Complex32[matrix.ColumnCount*matrix.ColumnCount];
#if CUDA
Complex32[] work = null;
#else
var work = new Complex32[100];
#endif
Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work);
var w = new DenseMatrix(matrix.RowCount, matrix.ColumnCount);
for (var index = 0; index < s.Length; index++)
{
w[index, index] = s[index];
}
var mU = new DenseMatrix(matrix.RowCount, matrix.RowCount, u);
var mV = new DenseMatrix(matrix.ColumnCount, matrix.ColumnCount, vt);
var result = mU*w*mV;
AssertHelpers.AlmostEqualRelative(matrix[0, 0], result[0, 0], 5);
AssertHelpers.AlmostEqualRelative(matrix[1, 0], result[1, 0], 5);
AssertHelpers.AlmostEqualRelative(matrix[0, 1], result[0, 1], 5);
AssertHelpers.AlmostEqualRelative(matrix[1, 1], result[1, 1], 5);
AssertHelpers.AlmostEqualRelative(matrix[0, 2], result[0, 2], 5);
AssertHelpers.AlmostEqualRelative(matrix[1, 2], result[1, 2], 4);
}
/// <summary>
/// Can solve Ax=b using SVD factorization with a square A matrix.
/// </summary>

558
src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs

@ -4,7 +4,7 @@
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 Math.NET
// Copyright (c) 2009-2015 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
@ -427,68 +427,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double
AssertHelpers.AlmostEqualRelative(a[8], -0.113636363636364, 13);
}
/// <summary>
/// Can compute the inverse of a matrix using LU factorization
/// with a work array.
/// </summary>
[Test]
public void CanComputeLuInverseWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new double[matrix.RowCount*matrix.RowCount];
Array.Copy(matrix.Values, a, a.Length);
#if CUDA
double[] work = null;
#else
var work = new double[matrix.RowCount];
#endif
Control.LinearAlgebraProvider.LUInverse(a, matrix.RowCount, work);
AssertHelpers.AlmostEqualRelative(a[0], -0.454545454545454, 13);
AssertHelpers.AlmostEqualRelative(a[1], -0.909090909090908, 13);
AssertHelpers.AlmostEqualRelative(a[2], 0.454545454545454, 13);
AssertHelpers.AlmostEqualRelative(a[3], -0.340909090909090, 13);
AssertHelpers.AlmostEqualRelative(a[4], -2.045454545454543, 13);
AssertHelpers.AlmostEqualRelative(a[5], 1.477272727272726, 13);
AssertHelpers.AlmostEqualRelative(a[6], -0.113636363636364, 13);
AssertHelpers.AlmostEqualRelative(a[7], 0.227272727272727, 13);
AssertHelpers.AlmostEqualRelative(a[8], -0.113636363636364, 13);
}
/// <summary>
/// Can compute the inverse of a matrix using LU factorization
/// using a previously factored matrix with a work array.
/// </summary>
[Test]
public void CanComputeLuInverseOnFactoredMatrixWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new double[matrix.RowCount*matrix.RowCount];
Array.Copy(matrix.Values, a, a.Length);
var ipiv = new int[matrix.RowCount];
Control.LinearAlgebraProvider.LUFactor(a, matrix.RowCount, ipiv);
#if CUDA
double[] work = null;
#else
var work = new double[matrix.RowCount];
#endif
Control.LinearAlgebraProvider.LUInverseFactored(a, matrix.RowCount, ipiv, work);
AssertHelpers.AlmostEqualRelative(a[0], -0.454545454545454, 13);
AssertHelpers.AlmostEqualRelative(a[1], -0.909090909090908, 13);
AssertHelpers.AlmostEqualRelative(a[2], 0.454545454545454, 13);
AssertHelpers.AlmostEqualRelative(a[3], -0.340909090909090, 13);
AssertHelpers.AlmostEqualRelative(a[4], -2.045454545454543, 13);
AssertHelpers.AlmostEqualRelative(a[5], 1.477272727272726, 13);
AssertHelpers.AlmostEqualRelative(a[6], -0.113636363636364, 13);
AssertHelpers.AlmostEqualRelative(a[7], 0.227272727272727, 13);
AssertHelpers.AlmostEqualRelative(a[8], -0.113636363636364, 13);
}
/// <summary>
/// Can solve Ax=b using LU factorization.
/// </summary>
@ -686,90 +624,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double
}
}
/// <summary>
/// Can compute QR factorization of a square matrix using a work array.
/// </summary>
[Test]
public void CanComputeQRFactorSquareMatrixWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var r = new double[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, r, r.Length);
var tau = new double[3];
var q = new double[matrix.RowCount*matrix.RowCount];
var work = new double[matrix.ColumnCount*Control.BlockSize];
Control.LinearAlgebraProvider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau, work);
var mq = Matrix<double>.Build.Dense(matrix.RowCount, matrix.RowCount, q);
var mr = Matrix<double>.Build.Dense(matrix.RowCount, matrix.ColumnCount, r).UpperTriangle();
var a = mq*mr;
for (var row = 0; row < matrix.RowCount; row++)
{
for (var col = 0; col < matrix.ColumnCount; col++)
{
AssertHelpers.AlmostEqualRelative(matrix[row, col], a[row, col], 14);
}
}
}
/// <summary>
/// Can compute QR factorization of a tall matrix using a work matrix.
/// </summary>
[Test]
public void CanComputeQRFactorTallMatrixWithWorkArray()
{
var matrix = _matrices["Tall3x2"];
var r = new double[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, r, r.Length);
var tau = new double[3];
var q = new double[matrix.RowCount*matrix.RowCount];
var work = new double[matrix.ColumnCount*Control.BlockSize];
Control.LinearAlgebraProvider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau, work);
var mr = Matrix<double>.Build.Dense(matrix.RowCount, matrix.ColumnCount, r).UpperTriangle();
var mq = Matrix<double>.Build.Dense(matrix.RowCount, matrix.RowCount, q);
var a = mq*mr;
for (var row = 0; row < matrix.RowCount; row++)
{
for (var col = 0; col < matrix.ColumnCount; col++)
{
AssertHelpers.AlmostEqualRelative(matrix[row, col], a[row, col], 14);
}
}
}
/// <summary>
/// Can compute QR factorization of a wide matrix using a work matrix.
/// </summary>
[Test]
public void CanComputeQRFactorWideMatrixWithWorkArray()
{
var matrix = _matrices["Wide2x3"];
var r = new double[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, r, r.Length);
var tau = new double[3];
var q = new double[matrix.RowCount*matrix.RowCount];
var work = new double[matrix.ColumnCount*Control.BlockSize];
Control.LinearAlgebraProvider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau, work);
var mr = Matrix<double>.Build.Dense(matrix.RowCount, matrix.ColumnCount, r).UpperTriangle();
var mq = Matrix<double>.Build.Dense(matrix.RowCount, matrix.RowCount, q);
var a = mq*mr;
for (var row = 0; row < matrix.RowCount; row++)
{
for (var col = 0; col < matrix.ColumnCount; col++)
{
AssertHelpers.AlmostEqualRelative(matrix[row, col], a[row, col], 14);
}
}
}
/// <summary>
/// Can compute thin QR factorization of a square matrix.
/// </summary>
@ -824,61 +678,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double
}
}
/// <summary>
/// Can compute thin QR factorization of a square matrix using a work array.
/// </summary>
[Test]
public void CanComputeThinQRFactorSquareMatrixWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var r = new double[matrix.ColumnCount*matrix.ColumnCount];
var tau = new double[3];
var q = new double[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, q, q.Length);
var work = new double[matrix.ColumnCount*Control.BlockSize];
Control.LinearAlgebraProvider.ThinQRFactor(q, matrix.RowCount, matrix.ColumnCount, r, tau, work);
var mq = Matrix<double>.Build.Dense(matrix.RowCount, matrix.ColumnCount, q);
var mr = Matrix<double>.Build.Dense(matrix.ColumnCount, matrix.ColumnCount, r);
var a = mq*mr;
for (var row = 0; row < matrix.RowCount; row++)
{
for (var col = 0; col < matrix.ColumnCount; col++)
{
AssertHelpers.AlmostEqualRelative(matrix[row, col], a[row, col], 14);
}
}
}
/// <summary>
/// Can compute thin QR factorization of a tall matrix using a work matrix.
/// </summary>
[Test]
public void CanComputeThinQRFactorTallMatrixWithWorkArray()
{
var matrix = _matrices["Tall3x2"];
var r = new double[matrix.ColumnCount*matrix.ColumnCount];
var tau = new double[3];
var q = new double[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, q, q.Length);
var work = new double[matrix.ColumnCount*Control.BlockSize];
Control.LinearAlgebraProvider.ThinQRFactor(q, matrix.RowCount, matrix.ColumnCount, r, tau, work);
var mq = Matrix<double>.Build.Dense(matrix.RowCount, matrix.ColumnCount, q);
var mr = Matrix<double>.Build.Dense(matrix.ColumnCount, matrix.ColumnCount, r);
var a = mq*mr;
for (var row = 0; row < matrix.RowCount; row++)
{
for (var col = 0; col < matrix.ColumnCount; col++)
{
AssertHelpers.AlmostEqualRelative(matrix[row, col], a[row, col], 14);
}
}
}
/// <summary>
/// Can solve Ax=b using QR factorization with a square A matrix.
/// </summary>
@ -931,62 +730,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 13);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a square A matrix
/// using a work array.
/// </summary>
[Test]
public void CanSolveUsingQRSquareMatrixUsingWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new double[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var b = new[] {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
var x = new double[matrix.ColumnCount*2];
var work = new double[matrix.RowCount*Control.BlockSize];
Control.LinearAlgebraProvider.QRSolve(a, matrix.RowCount, matrix.ColumnCount, b, 2, x, work);
NotModified(3, 3, a, matrix);
var mx = Matrix<double>.Build.Dense(matrix.ColumnCount, 2, x);
var mb = matrix*mx;
AssertHelpers.AlmostEqualRelative(mb[0, 0], b[0], 13);
AssertHelpers.AlmostEqualRelative(mb[1, 0], b[1], 13);
AssertHelpers.AlmostEqualRelative(mb[2, 0], b[2], 13);
AssertHelpers.AlmostEqualRelative(mb[0, 1], b[3], 13);
AssertHelpers.AlmostEqualRelative(mb[1, 1], b[4], 13);
AssertHelpers.AlmostEqualRelative(mb[2, 1], b[5], 13);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a tall A matrix
/// using a work array.
/// </summary>
[Test]
public void CanSolveUsingQRTallMatrixUsingWorkArray()
{
var matrix = _matrices["Tall3x2"];
var a = new double[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var b = new[] {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
var x = new double[matrix.ColumnCount*2];
var work = new double[matrix.RowCount*matrix.RowCount];
Control.LinearAlgebraProvider.QRSolve(a, matrix.RowCount, matrix.ColumnCount, b, 2, x, work);
NotModified(3, 2, a, matrix);
var mb = Matrix<double>.Build.Dense(matrix.RowCount, 2, b);
var test = (matrix.Transpose()*matrix).Inverse()*matrix.Transpose()*mb;
AssertHelpers.AlmostEqualRelative(test[0, 0], x[0], 13);
AssertHelpers.AlmostEqualRelative(test[1, 0], x[1], 13);
AssertHelpers.AlmostEqualRelative(test[0, 1], x[2], 13);
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 13);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a square A matrix
/// using a factored A matrix.
@ -1045,66 +788,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 13);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a square A matrix
/// using a factored A matrix with a work array.
/// </summary>
[Test]
public void CanSolveUsingQRSquareMatrixOnFactoredMatrixWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new double[matrix.RowCount*matrix.RowCount];
Array.Copy(matrix.Values, a, a.Length);
var tau = new double[matrix.ColumnCount];
var q = new double[matrix.ColumnCount*matrix.ColumnCount];
var work = new double[2048];
Control.LinearAlgebraProvider.QRFactor(a, matrix.RowCount, matrix.ColumnCount, q, tau, work);
var b = new[] {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
var x = new double[matrix.ColumnCount*2];
Control.LinearAlgebraProvider.QRSolveFactored(q, a, matrix.RowCount, matrix.ColumnCount, tau, b, 2, x, work);
var mx = Matrix<double>.Build.Dense(matrix.ColumnCount, 2, x);
var mb = matrix*mx;
AssertHelpers.AlmostEqualRelative(mb[0, 0], b[0], 13);
AssertHelpers.AlmostEqualRelative(mb[1, 0], b[1], 13);
AssertHelpers.AlmostEqualRelative(mb[2, 0], b[2], 13);
AssertHelpers.AlmostEqualRelative(mb[0, 1], b[3], 13);
AssertHelpers.AlmostEqualRelative(mb[1, 1], b[4], 13);
AssertHelpers.AlmostEqualRelative(mb[2, 1], b[5], 13);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a tall A matrix
/// using a factored A matrix with a work array.
/// </summary>
[Test]
public void CanSolveUsingQRTallMatrixOnFactoredMatrixWithWorkArray()
{
var matrix = _matrices["Tall3x2"];
var a = new double[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var tau = new double[matrix.ColumnCount];
var q = new double[matrix.RowCount*matrix.RowCount];
var work = new double[2048];
Control.LinearAlgebraProvider.QRFactor(a, matrix.RowCount, matrix.ColumnCount, q, tau, work);
var b = new[] {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
var x = new double[matrix.ColumnCount*2];
Control.LinearAlgebraProvider.QRSolveFactored(q, a, matrix.RowCount, matrix.ColumnCount, tau, b, 2, x, work);
var mb = Matrix<double>.Build.Dense(matrix.RowCount, 2, b);
var test = (matrix.Transpose()*matrix).Inverse()*matrix.Transpose()*mb;
AssertHelpers.AlmostEqualRelative(test[0, 0], x[0], 13);
AssertHelpers.AlmostEqualRelative(test[1, 0], x[1], 13);
AssertHelpers.AlmostEqualRelative(test[0, 1], x[2], 13);
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 13);
}
/// <summary>
/// Can solve Ax=b using thin QR factorization with a square A matrix.
/// </summary>
@ -1157,62 +840,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 13);
}
/// <summary>
/// Can solve Ax=b using thin QR factorization with a square A matrix
/// using a work array.
/// </summary>
[Test]
public void CanSolveUsingThinQRSquareMatrixUsingWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new double[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var b = new[] {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
var x = new double[matrix.ColumnCount*2];
var work = new double[matrix.RowCount*matrix.ColumnCount];
Control.LinearAlgebraProvider.QRSolve(a, matrix.RowCount, matrix.ColumnCount, b, 2, x, work, QRMethod.Thin);
NotModified(3, 3, a, matrix);
var mx = Matrix<double>.Build.Dense(matrix.ColumnCount, 2, x);
var mb = matrix*mx;
AssertHelpers.AlmostEqualRelative(mb[0, 0], b[0], 13);
AssertHelpers.AlmostEqualRelative(mb[1, 0], b[1], 13);
AssertHelpers.AlmostEqualRelative(mb[2, 0], b[2], 13);
AssertHelpers.AlmostEqualRelative(mb[0, 1], b[3], 13);
AssertHelpers.AlmostEqualRelative(mb[1, 1], b[4], 13);
AssertHelpers.AlmostEqualRelative(mb[2, 1], b[5], 13);
}
/// <summary>
/// Can solve Ax=b using thin QR factorization with a tall A matrix
/// using a work array.
/// </summary>
[Test]
public void CanSolveUsingThinQRTallMatrixUsingWorkArray()
{
var matrix = _matrices["Tall3x2"];
var a = new double[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var b = new[] {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
var x = new double[matrix.ColumnCount*2];
var work = new double[matrix.RowCount*matrix.ColumnCount];
Control.LinearAlgebraProvider.QRSolve(a, matrix.RowCount, matrix.ColumnCount, b, 2, x, work, QRMethod.Thin);
NotModified(3, 2, a, matrix);
var mb = Matrix<double>.Build.Dense(matrix.RowCount, 2, b);
var test = (matrix.Transpose()*matrix).Inverse()*matrix.Transpose()*mb;
AssertHelpers.AlmostEqualRelative(test[0, 0], x[0], 13);
AssertHelpers.AlmostEqualRelative(test[1, 0], x[1], 13);
AssertHelpers.AlmostEqualRelative(test[0, 1], x[2], 13);
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 13);
}
/// <summary>
/// Can solve Ax=b using thin QR factorization with a square A matrix
/// using a factored A matrix.
@ -1271,66 +898,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 13);
}
/// <summary>
/// Can solve Ax=b using thin QR factorization with a square A matrix
/// using a factored A matrix with a work array.
/// </summary>
[Test]
public void CanSolveUsingThinQRSquareMatrixOnFactoredMatrixWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new double[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var tau = new double[matrix.ColumnCount];
var r = new double[matrix.ColumnCount*matrix.ColumnCount];
var work = new double[2048];
Control.LinearAlgebraProvider.ThinQRFactor(a, matrix.RowCount, matrix.ColumnCount, r, tau, work);
var b = new[] {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
var x = new double[matrix.ColumnCount*2];
Control.LinearAlgebraProvider.QRSolveFactored(a, r, matrix.RowCount, matrix.ColumnCount, tau, b, 2, x, work, QRMethod.Thin);
var mx = Matrix<double>.Build.Dense(matrix.ColumnCount, 2, x);
var mb = matrix*mx;
AssertHelpers.AlmostEqualRelative(mb[0, 0], b[0], 13);
AssertHelpers.AlmostEqualRelative(mb[1, 0], b[1], 13);
AssertHelpers.AlmostEqualRelative(mb[2, 0], b[2], 13);
AssertHelpers.AlmostEqualRelative(mb[0, 1], b[3], 13);
AssertHelpers.AlmostEqualRelative(mb[1, 1], b[4], 13);
AssertHelpers.AlmostEqualRelative(mb[2, 1], b[5], 13);
}
/// <summary>
/// Can solve Ax=b using thin QR factorization with a tall A matrix
/// using a factored A matrix with a work array.
/// </summary>
[Test]
public void CanSolveUsingThinQRTallMatrixOnFactoredMatrixWithWorkArray()
{
var matrix = _matrices["Tall3x2"];
var a = new double[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var tau = new double[matrix.ColumnCount];
var r = new double[matrix.ColumnCount*matrix.ColumnCount];
var work = new double[2048];
Control.LinearAlgebraProvider.ThinQRFactor(a, matrix.RowCount, matrix.ColumnCount, r, tau, work);
var b = new[] {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
var x = new double[matrix.ColumnCount*2];
Control.LinearAlgebraProvider.QRSolveFactored(a, r, matrix.RowCount, matrix.ColumnCount, tau, b, 2, x, work, QRMethod.Thin);
var mb = Matrix<double>.Build.Dense(matrix.RowCount, 2, b);
var test = (matrix.Transpose()*matrix).Inverse()*matrix.Transpose()*mb;
AssertHelpers.AlmostEqualRelative(test[0, 0], x[0], 13);
AssertHelpers.AlmostEqualRelative(test[1, 0], x[1], 13);
AssertHelpers.AlmostEqualRelative(test[0, 1], x[2], 13);
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 13);
}
/// <summary>
/// Can compute the SVD factorization of a square matrix.
/// </summary>
@ -1436,129 +1003,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double
AssertHelpers.AlmostEqualRelative(matrix[1, 2], result[1, 2], 14);
}
/// <summary>
/// Can compute the SVD factorization of a square matrix using
/// a work array.
/// </summary>
[Test]
public void CanComputeSVDFactorizationOfSquareMatrixWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new double[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var s = new double[matrix.RowCount];
var u = new double[matrix.RowCount*matrix.RowCount];
var vt = new double[matrix.ColumnCount*matrix.ColumnCount];
#if CUDA
double[] work = null;
#else
var work = new double[100];
#endif
Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work);
var w = new DenseMatrix(matrix.RowCount, matrix.ColumnCount);
for (var index = 0; index < s.Length; index++)
{
w[index, index] = s[index];
}
var mU = Matrix<double>.Build.Dense(matrix.RowCount, matrix.RowCount, u);
var mV = Matrix<double>.Build.Dense(matrix.ColumnCount, matrix.ColumnCount, vt);
var result = mU*w*mV;
AssertHelpers.AlmostEqualRelative(matrix[0, 0], result[0, 0], 14);
AssertHelpers.AlmostEqualRelative(matrix[1, 0], result[1, 0], 14);
AssertHelpers.AlmostEqualRelative(matrix[2, 0], result[2, 0], 14);
AssertHelpers.AlmostEqualRelative(matrix[0, 1], result[0, 1], 14);
AssertHelpers.AlmostEqualRelative(matrix[1, 1], result[1, 1], 14);
AssertHelpers.AlmostEqualRelative(matrix[2, 1], result[2, 1], 14);
AssertHelpers.AlmostEqualRelative(matrix[0, 2], result[0, 2], 14);
AssertHelpers.AlmostEqualRelative(matrix[1, 2], result[1, 2], 14);
AssertHelpers.AlmostEqualRelative(matrix[2, 2], result[2, 2], 14);
}
/// <summary>
/// Can compute the SVD factorization of a tall matrix using
/// a work array.
/// </summary>
[Test]
public void CanComputeSVDFactorizationOfTallMatrixWithWorkArray()
{
var matrix = _matrices["Tall3x2"];
var a = new double[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var s = new double[matrix.ColumnCount];
var u = new double[matrix.RowCount*matrix.RowCount];
var vt = new double[matrix.ColumnCount*matrix.ColumnCount];
#if CUDA
double[] work = null;
#else
var work = new double[100];
#endif
Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work);
var w = new DenseMatrix(matrix.RowCount, matrix.ColumnCount);
for (var index = 0; index < s.Length; index++)
{
w[index, index] = s[index];
}
var mU = Matrix<double>.Build.Dense(matrix.RowCount, matrix.RowCount, u);
var mV = Matrix<double>.Build.Dense(matrix.ColumnCount, matrix.ColumnCount, vt);
var result = mU*w*mV;
AssertHelpers.AlmostEqualRelative(matrix[0, 0], result[0, 0], 14);
AssertHelpers.AlmostEqualRelative(matrix[1, 0], result[1, 0], 14);
AssertHelpers.AlmostEqualRelative(matrix[2, 0], result[2, 0], 14);
AssertHelpers.AlmostEqualRelative(matrix[0, 1], result[0, 1], 14);
AssertHelpers.AlmostEqualRelative(matrix[1, 1], result[1, 1], 14);
AssertHelpers.AlmostEqualRelative(matrix[2, 1], result[2, 1], 14);
}
/// <summary>
/// Can compute the SVD factorization of a wide matrix using
/// a work array.
/// </summary>
[Test]
public void CanComputeSVDFactorizationOfWideMatrixWithWorkArray()
{
var matrix = _matrices["Wide2x3"];
var a = new double[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var s = new double[matrix.RowCount];
var u = new double[matrix.RowCount*matrix.RowCount];
var vt = new double[matrix.ColumnCount*matrix.ColumnCount];
#if CUDA
double[] work = null;
#else
var work = new double[100];
#endif
Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work);
var w = new DenseMatrix(matrix.RowCount, matrix.ColumnCount);
for (var index = 0; index < s.Length; index++)
{
w[index, index] = s[index];
}
var mU = Matrix<double>.Build.Dense(matrix.RowCount, matrix.RowCount, u);
var mV = Matrix<double>.Build.Dense(matrix.ColumnCount, matrix.ColumnCount, vt);
var result = mU*w*mV;
AssertHelpers.AlmostEqualRelative(matrix[0, 0], result[0, 0], 14);
AssertHelpers.AlmostEqualRelative(matrix[1, 0], result[1, 0], 14);
AssertHelpers.AlmostEqualRelative(matrix[0, 1], result[0, 1], 14);
AssertHelpers.AlmostEqualRelative(matrix[1, 1], result[1, 1], 14);
AssertHelpers.AlmostEqualRelative(matrix[0, 2], result[0, 2], 14);
AssertHelpers.AlmostEqualRelative(matrix[1, 2], result[1, 2], 14);
}
/// <summary>
/// Can solve Ax=b using SVD factorization with a square A matrix.
/// </summary>

558
src/UnitTests/LinearAlgebraProviderTests/Single/LinearAlgebraProviderTests.cs

@ -4,7 +4,7 @@
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 Math.NET
// Copyright (c) 2009-2015 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
@ -435,68 +435,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Single
AssertHelpers.AlmostEqual(a[8], -0.113636363636364, 5);
}
/// <summary>
/// Can compute the inverse of a matrix using LU factorization
/// with a work array.
/// </summary>
[Test]
public void CanComputeLuInverseWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new float[matrix.RowCount*matrix.RowCount];
Array.Copy(matrix.Values, a, a.Length);
#if CUDA
float[] work = null;
#else
var work = new float[matrix.RowCount];
#endif
Control.LinearAlgebraProvider.LUInverse(a, matrix.RowCount, work);
AssertHelpers.AlmostEqual(a[0], -0.454545454545454, 5);
AssertHelpers.AlmostEqual(a[1], -0.909090909090908, 5);
AssertHelpers.AlmostEqual(a[2], 0.454545454545454, 5);
AssertHelpers.AlmostEqual(a[3], -0.340909090909090, 5);
AssertHelpers.AlmostEqual(a[4], -2.045454545454543, 5);
AssertHelpers.AlmostEqual(a[5], 1.477272727272726, 5);
AssertHelpers.AlmostEqual(a[6], -0.113636363636364, 5);
AssertHelpers.AlmostEqual(a[7], 0.227272727272727, 5);
AssertHelpers.AlmostEqual(a[8], -0.113636363636364, 5);
}
/// <summary>
/// Can compute the inverse of a matrix using LU factorization
/// using a previously factored matrix with a work array.
/// </summary>
[Test]
public void CanComputeLuInverseOnFactoredMatrixWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new float[matrix.RowCount*matrix.RowCount];
Array.Copy(matrix.Values, a, a.Length);
var ipiv = new int[matrix.RowCount];
Control.LinearAlgebraProvider.LUFactor(a, matrix.RowCount, ipiv);
#if CUDA
float[] work = null;
#else
var work = new float[matrix.RowCount];
#endif
Control.LinearAlgebraProvider.LUInverseFactored(a, matrix.RowCount, ipiv, work);
AssertHelpers.AlmostEqual(a[0], -0.454545454545454, 5);
AssertHelpers.AlmostEqual(a[1], -0.909090909090908, 5);
AssertHelpers.AlmostEqual(a[2], 0.454545454545454, 5);
AssertHelpers.AlmostEqual(a[3], -0.340909090909090, 5);
AssertHelpers.AlmostEqual(a[4], -2.045454545454543, 5);
AssertHelpers.AlmostEqual(a[5], 1.477272727272726, 5);
AssertHelpers.AlmostEqual(a[6], -0.113636363636364, 5);
AssertHelpers.AlmostEqual(a[7], 0.227272727272727, 5);
AssertHelpers.AlmostEqual(a[8], -0.113636363636364, 5);
}
/// <summary>
/// Can solve Ax=b using LU factorization.
/// </summary>
@ -694,90 +632,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Single
}
}
/// <summary>
/// Can compute QR factorization of a square matrix using a work array.
/// </summary>
[Test]
public void CanComputeQRFactorSquareMatrixWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var r = new float[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, r, r.Length);
var tau = new float[3];
var q = new float[matrix.RowCount*matrix.RowCount];
var work = new float[matrix.ColumnCount*Control.BlockSize];
Control.LinearAlgebraProvider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau, work);
var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q);
var mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r).UpperTriangle();
var a = mq*mr;
for (var row = 0; row < matrix.RowCount; row++)
{
for (var col = 0; col < matrix.ColumnCount; col++)
{
AssertHelpers.AlmostEqualRelative(matrix[row, col], a[row, col], 5);
}
}
}
/// <summary>
/// Can compute QR factorization of a tall matrix using a work matrix.
/// </summary>
[Test]
public void CanComputeQRFactorTallMatrixWithWorkArray()
{
var matrix = _matrices["Tall3x2"];
var r = new float[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, r, r.Length);
var tau = new float[3];
var q = new float[matrix.RowCount*matrix.RowCount];
var work = new float[matrix.ColumnCount*Control.BlockSize];
Control.LinearAlgebraProvider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau, work);
var mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r).UpperTriangle();
var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q);
var a = mq*mr;
for (var row = 0; row < matrix.RowCount; row++)
{
for (var col = 0; col < matrix.ColumnCount; col++)
{
AssertHelpers.AlmostEqualRelative(matrix[row, col], a[row, col], 5);
}
}
}
/// <summary>
/// Can compute QR factorization of a wide matrix using a work matrix.
/// </summary>
[Test]
public void CanComputeQRFactorWideMatrixWithWorkArray()
{
var matrix = _matrices["Wide2x3"];
var r = new float[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, r, r.Length);
var tau = new float[3];
var q = new float[matrix.RowCount*matrix.RowCount];
var work = new float[matrix.ColumnCount*Control.BlockSize];
Control.LinearAlgebraProvider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau, work);
var mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r).UpperTriangle();
var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q);
var a = mq*mr;
for (var row = 0; row < matrix.RowCount; row++)
{
for (var col = 0; col < matrix.ColumnCount; col++)
{
AssertHelpers.AlmostEqualRelative(matrix[row, col], a[row, col], 5);
}
}
}
/// <summary>
/// Can compute thin QR factorization of a square matrix.
/// </summary>
@ -832,61 +686,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Single
}
}
/// <summary>
/// Can compute thin QR factorization of a square matrix using a work array.
/// </summary>
[Test]
public void CanComputeThinQRFactorSquareMatrixWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var r = new float[matrix.ColumnCount*matrix.ColumnCount];
var tau = new float[3];
var q = new float[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, q, q.Length);
var work = new float[matrix.ColumnCount*Control.BlockSize];
Control.LinearAlgebraProvider.ThinQRFactor(q, matrix.RowCount, matrix.ColumnCount, r, tau, work);
var mq = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, q);
var mr = new DenseMatrix(matrix.ColumnCount, matrix.ColumnCount, r);
var a = mq*mr;
for (var row = 0; row < matrix.RowCount; row++)
{
for (var col = 0; col < matrix.ColumnCount; col++)
{
AssertHelpers.AlmostEqualRelative(matrix[row, col], a[row, col], 5);
}
}
}
/// <summary>
/// Can compute thin QR factorization of a tall matrix using a work matrix.
/// </summary>
[Test]
public void CanComputeThinQRFactorTallMatrixWithWorkArray()
{
var matrix = _matrices["Tall3x2"];
var r = new float[matrix.ColumnCount*matrix.ColumnCount];
var tau = new float[3];
var q = new float[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, q, q.Length);
var work = new float[matrix.ColumnCount*Control.BlockSize];
Control.LinearAlgebraProvider.ThinQRFactor(q, matrix.RowCount, matrix.ColumnCount, r, tau, work);
var mq = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, q);
var mr = new DenseMatrix(matrix.ColumnCount, matrix.ColumnCount, r);
var a = mq*mr;
for (var row = 0; row < matrix.RowCount; row++)
{
for (var col = 0; col < matrix.ColumnCount; col++)
{
AssertHelpers.AlmostEqualRelative(matrix[row, col], a[row, col], 5);
}
}
}
/// <summary>
/// Can solve Ax=b using QR factorization with a square A matrix.
/// </summary>
@ -939,62 +738,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Single
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 5);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a square A matrix
/// using a work array.
/// </summary>
[Test]
public void CanSolveUsingQRSquareMatrixUsingWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new float[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var b = new[] {1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f};
var x = new float[matrix.ColumnCount*2];
var work = new float[matrix.RowCount*matrix.RowCount];
Control.LinearAlgebraProvider.QRSolve(a, matrix.RowCount, matrix.ColumnCount, b, 2, x, work);
NotModified(3, 3, a, matrix);
var mx = new DenseMatrix(matrix.ColumnCount, 2, x);
var mb = matrix*mx;
AssertHelpers.AlmostEqualRelative(mb[0, 0], b[0], 5);
AssertHelpers.AlmostEqualRelative(mb[1, 0], b[1], 5);
AssertHelpers.AlmostEqualRelative(mb[2, 0], b[2], 5);
AssertHelpers.AlmostEqualRelative(mb[0, 1], b[3], 5);
AssertHelpers.AlmostEqualRelative(mb[1, 1], b[4], 4);
AssertHelpers.AlmostEqualRelative(mb[2, 1], b[5], 4);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a tall A matrix
/// using a work array.
/// </summary>
[Test]
public void CanSolveUsingQRTallMatrixUsingWorkArray()
{
var matrix = _matrices["Tall3x2"];
var a = new float[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var b = new[] {1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f};
var x = new float[matrix.ColumnCount*2];
var work = new float[matrix.RowCount*matrix.RowCount];
Control.LinearAlgebraProvider.QRSolve(a, matrix.RowCount, matrix.ColumnCount, b, 2, x, work);
NotModified(3, 2, a, matrix);
var mb = new DenseMatrix(matrix.RowCount, 2, b);
var test = (matrix.Transpose()*matrix).Inverse()*matrix.Transpose()*mb;
AssertHelpers.AlmostEqualRelative(test[0, 0], x[0], 5);
AssertHelpers.AlmostEqualRelative(test[1, 0], x[1], 5);
AssertHelpers.AlmostEqualRelative(test[0, 1], x[2], 5);
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 5);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a square A matrix
/// using a factored A matrix.
@ -1053,66 +796,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Single
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 5);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a square A matrix
/// using a factored A matrix with a work array.
/// </summary>
[Test]
public void CanSolveUsingQRSquareMatrixOnFactoredMatrixWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new float[matrix.RowCount*matrix.RowCount];
Array.Copy(matrix.Values, a, a.Length);
var tau = new float[matrix.ColumnCount];
var q = new float[matrix.ColumnCount*matrix.ColumnCount];
var work = new float[2048];
Control.LinearAlgebraProvider.QRFactor(a, matrix.RowCount, matrix.ColumnCount, q, tau, work);
var b = new[] {1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f};
var x = new float[matrix.ColumnCount*2];
Control.LinearAlgebraProvider.QRSolveFactored(q, a, matrix.RowCount, matrix.ColumnCount, tau, b, 2, x, work);
var mx = new DenseMatrix(matrix.ColumnCount, 2, x);
var mb = matrix*mx;
AssertHelpers.AlmostEqualRelative(mb[0, 0], b[0], 5);
AssertHelpers.AlmostEqualRelative(mb[1, 0], b[1], 5);
AssertHelpers.AlmostEqualRelative(mb[2, 0], b[2], 5);
AssertHelpers.AlmostEqualRelative(mb[0, 1], b[3], 5);
AssertHelpers.AlmostEqualRelative(mb[1, 1], b[4], 4);
AssertHelpers.AlmostEqualRelative(mb[2, 1], b[5], 4);
}
/// <summary>
/// Can solve Ax=b using QR factorization with a tall A matrix
/// using a factored A matrix with a work array.
/// </summary>
[Test]
public void CanSolveUsingQRTallMatrixOnFactoredMatrixWithWorkArray()
{
var matrix = _matrices["Tall3x2"];
var a = new float[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var tau = new float[matrix.ColumnCount];
var q = new float[matrix.RowCount*matrix.RowCount];
var work = new float[2048];
Control.LinearAlgebraProvider.QRFactor(a, matrix.RowCount, matrix.ColumnCount, q, tau, work);
var b = new[] {1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f};
var x = new float[matrix.ColumnCount*2];
Control.LinearAlgebraProvider.QRSolveFactored(q, a, matrix.RowCount, matrix.ColumnCount, tau, b, 2, x, work);
var mb = new DenseMatrix(matrix.RowCount, 2, b);
var test = (matrix.Transpose()*matrix).Inverse()*matrix.Transpose()*mb;
AssertHelpers.AlmostEqualRelative(test[0, 0], x[0], 5);
AssertHelpers.AlmostEqualRelative(test[1, 0], x[1], 5);
AssertHelpers.AlmostEqualRelative(test[0, 1], x[2], 5);
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 5);
}
/// <summary>
/// Can solve Ax=b using thin QR factorization with a square A matrix.
/// </summary>
@ -1165,62 +848,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Single
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 5);
}
/// <summary>
/// Can solve Ax=b using thin QR factorization with a square A matrix
/// using a work array.
/// </summary>
[Test]
public void CanSolveUsingThinQRSquareMatrixUsingWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new float[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var b = new[] {1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f};
var x = new float[matrix.ColumnCount*2];
var work = new float[matrix.RowCount*matrix.ColumnCount];
Control.LinearAlgebraProvider.QRSolve(a, matrix.RowCount, matrix.ColumnCount, b, 2, x, work, QRMethod.Thin);
NotModified(3, 3, a, matrix);
var mx = new DenseMatrix(matrix.ColumnCount, 2, x);
var mb = matrix*mx;
AssertHelpers.AlmostEqualRelative(mb[0, 0], b[0], 5);
AssertHelpers.AlmostEqualRelative(mb[1, 0], b[1], 5);
AssertHelpers.AlmostEqualRelative(mb[2, 0], b[2], 5);
AssertHelpers.AlmostEqualRelative(mb[0, 1], b[3], 5);
AssertHelpers.AlmostEqualRelative(mb[1, 1], b[4], 4);
AssertHelpers.AlmostEqualRelative(mb[2, 1], b[5], 4);
}
/// <summary>
/// Can solve Ax=b using thin QR factorization with a tall A matrix
/// using a work array.
/// </summary>
[Test]
public void CanSolveUsingThinQRTallMatrixUsingWorkArray()
{
var matrix = _matrices["Tall3x2"];
var a = new float[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var b = new[] {1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f};
var x = new float[matrix.ColumnCount*2];
var work = new float[matrix.RowCount*matrix.ColumnCount];
Control.LinearAlgebraProvider.QRSolve(a, matrix.RowCount, matrix.ColumnCount, b, 2, x, work, QRMethod.Thin);
NotModified(3, 2, a, matrix);
var mb = new DenseMatrix(matrix.RowCount, 2, b);
var test = (matrix.Transpose()*matrix).Inverse()*matrix.Transpose()*mb;
AssertHelpers.AlmostEqualRelative(test[0, 0], x[0], 5);
AssertHelpers.AlmostEqualRelative(test[1, 0], x[1], 5);
AssertHelpers.AlmostEqualRelative(test[0, 1], x[2], 5);
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 5);
}
/// <summary>
/// Can solve Ax=b using thin QR factorization with a square A matrix
/// using a factored A matrix.
@ -1279,66 +906,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Single
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 5);
}
/// <summary>
/// Can solve Ax=b using thin QR factorization with a square A matrix
/// using a factored A matrix with a work array.
/// </summary>
[Test]
public void CanSolveUsingThinQRSquareMatrixOnFactoredMatrixWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new float[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var tau = new float[matrix.ColumnCount];
var r = new float[matrix.ColumnCount*matrix.ColumnCount];
var work = new float[2048];
Control.LinearAlgebraProvider.ThinQRFactor(a, matrix.RowCount, matrix.ColumnCount, r, tau, work);
var b = new[] {1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f};
var x = new float[matrix.ColumnCount*2];
Control.LinearAlgebraProvider.QRSolveFactored(a, r, matrix.RowCount, matrix.ColumnCount, tau, b, 2, x, work, QRMethod.Thin);
var mx = new DenseMatrix(matrix.ColumnCount, 2, x);
var mb = matrix*mx;
AssertHelpers.AlmostEqualRelative(mb[0, 0], b[0], 5);
AssertHelpers.AlmostEqualRelative(mb[1, 0], b[1], 5);
AssertHelpers.AlmostEqualRelative(mb[2, 0], b[2], 5);
AssertHelpers.AlmostEqualRelative(mb[0, 1], b[3], 5);
AssertHelpers.AlmostEqualRelative(mb[1, 1], b[4], 4);
AssertHelpers.AlmostEqualRelative(mb[2, 1], b[5], 4);
}
/// <summary>
/// Can solve Ax=b using thin QR factorization with a tall A matrix
/// using a factored A matrix with a work array.
/// </summary>
[Test]
public void CanSolveUsingThinQRTallMatrixOnFactoredMatrixWithWorkArray()
{
var matrix = _matrices["Tall3x2"];
var a = new float[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var tau = new float[matrix.ColumnCount];
var r = new float[matrix.ColumnCount*matrix.ColumnCount];
var work = new float[2048];
Control.LinearAlgebraProvider.ThinQRFactor(a, matrix.RowCount, matrix.ColumnCount, r, tau, work);
var b = new[] {1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f};
var x = new float[matrix.ColumnCount*2];
Control.LinearAlgebraProvider.QRSolveFactored(a, r, matrix.RowCount, matrix.ColumnCount, tau, b, 2, x, work, QRMethod.Thin);
var mb = new DenseMatrix(matrix.RowCount, 2, b);
var test = (matrix.Transpose()*matrix).Inverse()*matrix.Transpose()*mb;
AssertHelpers.AlmostEqualRelative(test[0, 0], x[0], 5);
AssertHelpers.AlmostEqualRelative(test[1, 0], x[1], 5);
AssertHelpers.AlmostEqualRelative(test[0, 1], x[2], 5);
AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 5);
}
/// <summary>
/// Can compute the SVD factorization of a square matrix.
/// </summary>
@ -1444,129 +1011,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Single
AssertHelpers.AlmostEqualRelative(matrix[1, 2], result[1, 2], 5);
}
/// <summary>
/// Can compute the SVD factorization of a square matrix using
/// a work array.
/// </summary>
[Test]
public void CanComputeSVDFactorizationOfSquareMatrixWithWorkArray()
{
var matrix = _matrices["Square3x3"];
var a = new float[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var s = new float[matrix.RowCount];
var u = new float[matrix.RowCount*matrix.RowCount];
var vt = new float[matrix.ColumnCount*matrix.ColumnCount];
#if CUDA
float[] work = null;
#else
var work = new float[100];
#endif
Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work);
var w = new DenseMatrix(matrix.RowCount, matrix.ColumnCount);
for (var index = 0; index < s.Length; index++)
{
w[index, index] = s[index];
}
var mU = new DenseMatrix(matrix.RowCount, matrix.RowCount, u);
var mV = new DenseMatrix(matrix.ColumnCount, matrix.ColumnCount, vt);
var result = mU*w*mV;
AssertHelpers.AlmostEqualRelative(matrix[0, 0], result[0, 0], 5);
AssertHelpers.AlmostEqualRelative(matrix[1, 0], result[1, 0], 5);
AssertHelpers.AlmostEqualRelative(matrix[2, 0], result[2, 0], 5);
AssertHelpers.AlmostEqualRelative(matrix[0, 1], result[0, 1], 5);
AssertHelpers.AlmostEqualRelative(matrix[1, 1], result[1, 1], 5);
AssertHelpers.AlmostEqualRelative(matrix[2, 1], result[2, 1], 5);
AssertHelpers.AlmostEqualRelative(matrix[0, 2], result[0, 2], 5);
AssertHelpers.AlmostEqualRelative(matrix[1, 2], result[1, 2], 5);
AssertHelpers.AlmostEqualRelative(matrix[2, 2], result[2, 2], 5);
}
/// <summary>
/// Can compute the SVD factorization of a tall matrix using
/// a work array.
/// </summary>
[Test]
public void CanComputeSVDFactorizationOfTallMatrixWithWorkArray()
{
var matrix = _matrices["Tall3x2"];
var a = new float[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var s = new float[matrix.ColumnCount];
var u = new float[matrix.RowCount*matrix.RowCount];
var vt = new float[matrix.ColumnCount*matrix.ColumnCount];
#if CUDA
float[] work = null;
#else
var work = new float[100];
#endif
Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work);
var w = new DenseMatrix(matrix.RowCount, matrix.ColumnCount);
for (var index = 0; index < s.Length; index++)
{
w[index, index] = s[index];
}
var mU = new DenseMatrix(matrix.RowCount, matrix.RowCount, u);
var mV = new DenseMatrix(matrix.ColumnCount, matrix.ColumnCount, vt);
var result = mU*w*mV;
AssertHelpers.AlmostEqualRelative(matrix[0, 0], result[0, 0], 5);
AssertHelpers.AlmostEqualRelative(matrix[1, 0], result[1, 0], 5);
AssertHelpers.AlmostEqualRelative(matrix[2, 0], result[2, 0], 5);
AssertHelpers.AlmostEqualRelative(matrix[0, 1], result[0, 1], 5);
AssertHelpers.AlmostEqualRelative(matrix[1, 1], result[1, 1], 5);
AssertHelpers.AlmostEqualRelative(matrix[2, 1], result[2, 1], 5);
}
/// <summary>
/// Can compute the SVD factorization of a wide matrix using
/// a work array.
/// </summary>
[Test]
public void CanComputeSVDFactorizationOfWideMatrixWithWorkArray()
{
var matrix = _matrices["Wide2x3"];
var a = new float[matrix.RowCount*matrix.ColumnCount];
Array.Copy(matrix.Values, a, a.Length);
var s = new float[matrix.RowCount];
var u = new float[matrix.RowCount*matrix.RowCount];
var vt = new float[matrix.ColumnCount*matrix.ColumnCount];
#if CUDA
float[] work = null;
#else
var work = new float[100];
#endif
Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work);
var w = new DenseMatrix(matrix.RowCount, matrix.ColumnCount);
for (var index = 0; index < s.Length; index++)
{
w[index, index] = s[index];
}
var mU = new DenseMatrix(matrix.RowCount, matrix.RowCount, u);
var mV = new DenseMatrix(matrix.ColumnCount, matrix.ColumnCount, vt);
var result = mU*w*mV;
AssertHelpers.AlmostEqualRelative(matrix[0, 0], result[0, 0], 5);
AssertHelpers.AlmostEqualRelative(matrix[1, 0], result[1, 0], 5);
AssertHelpers.AlmostEqualRelative(matrix[0, 1], result[0, 1], 5);
AssertHelpers.AlmostEqualRelative(matrix[1, 1], result[1, 1], 5);
AssertHelpers.AlmostEqualRelative(matrix[0, 2], result[0, 2], 5);
AssertHelpers.AlmostEqualRelative(matrix[1, 2], result[1, 2], 5);
}
/// <summary>
/// Can solve Ax=b using SVD factorization with a square A matrix.
/// </summary>

4
src/UnitTests/UnitTests-CUDA.csproj

@ -336,6 +336,10 @@
<Project>{b7cae5f4-a23f-4438-b5be-41226618b695}</Project>
<Name>Numerics</Name>
</ProjectReference>
<ProjectReference Include="..\TestData\TestData.csproj">
<Project>{af3253c9-4db5-45a0-98cf-c105fda9da47}</Project>
<Name>TestData</Name>
</ProjectReference>
</ItemGroup>
<Import Project="$(MSBuildToolsPath)\Microsoft.CSharp.targets" />
<ItemGroup>

15
src/UnitTests/UnitTests-MKL.csproj

@ -130,6 +130,10 @@
<Compile Include="InterpolationTests\LinearSplineTest.cs" />
<Compile Include="InterpolationTests\NevillePolynomialTest.cs" />
<Compile Include="InterpolationTests\StepInterpolationTest.cs" />
<Compile Include="LinearAlgebraProviderTests\Complex32\LinearAlgebraProviderTests.cs" />
<Compile Include="LinearAlgebraProviderTests\Complex\LinearAlgebraProviderTests.cs" />
<Compile Include="LinearAlgebraProviderTests\Double\LinearAlgebraProviderTests.cs" />
<Compile Include="LinearAlgebraProviderTests\Single\LinearAlgebraProviderTests.cs" />
<Compile Include="LinearAlgebraTests\Complex32\DenseMatrixTests.cs" />
<Compile Include="LinearAlgebraTests\Complex32\DenseVectorTests.cs" />
<Compile Include="LinearAlgebraTests\Complex32\DiagonalMatrixTests.cs" />
@ -652,6 +656,10 @@
<Project>{b7cae5f4-a23f-4438-b5be-41226618b695}</Project>
<Name>Numerics</Name>
</ProjectReference>
<ProjectReference Include="..\TestData\TestData.csproj">
<Project>{af3253c9-4db5-45a0-98cf-c105fda9da47}</Project>
<Name>TestData</Name>
</ProjectReference>
</ItemGroup>
<Import Project="$(MSBuildToolsPath)\Microsoft.CSharp.targets" />
<ItemGroup>
@ -664,10 +672,5 @@
<ItemGroup>
<Service Include="{82A7F48D-3B50-4B1E-B82E-3ADA8210C358}" />
</ItemGroup>
<ItemGroup>
<Folder Include="LinearAlgebraProviderTests\Complex32\" />
<Folder Include="LinearAlgebraProviderTests\Complex\" />
<Folder Include="LinearAlgebraProviderTests\Double\" />
<Folder Include="LinearAlgebraProviderTests\Single\" />
</ItemGroup>
<ItemGroup />
</Project>

4
src/UnitTests/UnitTests-OpenBLAS.csproj

@ -336,6 +336,10 @@
<Project>{b7cae5f4-a23f-4438-b5be-41226618b695}</Project>
<Name>Numerics</Name>
</ProjectReference>
<ProjectReference Include="..\TestData\TestData.csproj">
<Project>{af3253c9-4db5-45a0-98cf-c105fda9da47}</Project>
<Name>TestData</Name>
</ProjectReference>
</ItemGroup>
<Import Project="$(MSBuildToolsPath)\Microsoft.CSharp.targets" />
<ItemGroup>

Loading…
Cancel
Save