diff --git a/MathNet.Numerics.NativeProviders.sln b/MathNet.Numerics.NativeProviders.sln index 4e790037..21805ccb 100644 --- a/MathNet.Numerics.NativeProviders.sln +++ b/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 diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex.cs b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex.cs index 31abfbc5..ac729107 100644 --- a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex.cs +++ b/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)); } - /// - /// Computes the inverse of matrix using LU factorization. - /// - /// The N by N matrix to invert. Contains the inverse On exit. - /// The order of the square matrix . - /// Not supported. Should be left null. - /// This is equivalent to the GETRF and GETRI LAPACK routines. - [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)); - } - - /// - /// Computes the inverse of a previously factored matrix. - /// - /// The LU factored N by N matrix. Contains the inverse On exit. - /// The order of the square matrix . - /// The pivot indices of . - /// Not supported. Should be left null. - /// This is equivalent to the GETRI LAPACK routine. - [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)); - } - /// /// Solves A*X=B for X using LU factorization. /// @@ -538,60 +470,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda Solver(SafeNativeMethods.z_cholesky_solve_factored(_solverHandle, orderA, columnsB, a, b)); } - /// - /// Computes the singular value decomposition of A. - /// - /// Compute the singular U and VT vectors or not. - /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The singular values of A in ascending value. - /// If is true, on exit U contains the left - /// singular vectors. - /// If is true, on exit VT contains the transposed - /// right singular vectors. - /// This is equivalent to the GESVD LAPACK routine. - [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt) - { - 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); - } - /// /// Solves A*X=B for X using the singular value decomposition of A. /// @@ -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. /// If is true, on exit VT contains the transposed /// right singular vectors. - /// User work buffers are not supported. Should be null. /// This is equivalent to the GESVD LAPACK routine. [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt, 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)); } } diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex32.cs b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex32.cs index 81679c15..01b8c46e 100644 --- a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex32.cs +++ b/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)); } - /// - /// Computes the inverse of matrix using LU factorization. - /// - /// The N by N matrix to invert. Contains the inverse On exit. - /// The order of the square matrix . - /// Not supported. Should be left null. - /// This is equivalent to the GETRF and GETRI LAPACK routines. - [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)); - } - - /// - /// Computes the inverse of a previously factored matrix. - /// - /// The LU factored N by N matrix. Contains the inverse On exit. - /// The order of the square matrix . - /// The pivot indices of . - /// Not supported. Should be left null. - /// This is equivalent to the GETRI LAPACK routine. - [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)); - } - /// /// Solves A*X=B for X using LU factorization. /// @@ -538,60 +469,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda Solver(SafeNativeMethods.c_cholesky_solve_factored(_solverHandle, orderA, columnsB, a, b)); } - /// - /// Computes the singular value decomposition of A. - /// - /// Compute the singular U and VT vectors or not. - /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The singular values of A in ascending value. - /// If is true, on exit U contains the left - /// singular vectors. - /// If is true, on exit VT contains the transposed - /// right singular vectors. - /// This is equivalent to the GESVD LAPACK routine. - [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt) - { - 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); - } - /// /// Solves A*X=B for X using the singular value decomposition of A. /// @@ -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. /// If is true, on exit VT contains the transposed /// right singular vectors. - /// Not supported. Should be left null. /// This is equivalent to the GESVD LAPACK routine. [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt, 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)); } } diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Double.cs b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Double.cs index 6e73ed46..0f462589 100644 --- a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Double.cs +++ b/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)); } - /// - /// Computes the inverse of matrix using LU factorization. - /// - /// The N by N matrix to invert. Contains the inverse On exit. - /// The order of the square matrix . - /// Not supported. Should be left null. - /// This is equivalent to the GETRF and GETRI LAPACK routines. - [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)); - } - - /// - /// Computes the inverse of a previously factored matrix. - /// - /// The LU factored N by N matrix. Contains the inverse On exit. - /// The order of the square matrix . - /// The pivot indices of . - /// Not supported. Should be left null. - /// This is equivalent to the GETRI LAPACK routine. - [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)); - } - /// /// Solves A*X=B for X using LU factorization. /// @@ -538,60 +469,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda Solver(SafeNativeMethods.d_cholesky_solve_factored(_solverHandle, orderA, columnsB, a, b)); } - /// - /// Computes the singular value decomposition of A. - /// - /// Compute the singular U and VT vectors or not. - /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The singular values of A in ascending value. - /// If is true, on exit U contains the left - /// singular vectors. - /// If is true, on exit VT contains the transposed - /// right singular vectors. - /// This is equivalent to the GESVD LAPACK routine. - [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt) - { - 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); - } - /// /// Solves A*X=B for X using the singular value decomposition of A. /// @@ -650,10 +527,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda /// singular vectors. /// If is true, on exit VT contains the transposed /// right singular vectors. - /// Not supported. Should be left null. /// This is equivalent to the GESVD LAPACK routine. [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt, 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)); } } diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Single.cs b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Single.cs index 0ac4e34f..c68dc1c5 100644 --- a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Single.cs +++ b/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)); } - /// - /// Computes the inverse of matrix using LU factorization. - /// - /// The N by N matrix to invert. Contains the inverse On exit. - /// The order of the square matrix . - /// Not supported. Should be left null. - /// This is equivalent to the GETRF and GETRI LAPACK routines. - [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)); - } - - /// - /// Computes the inverse of a previously factored matrix. - /// - /// The LU factored N by N matrix. Contains the inverse On exit. - /// The order of the square matrix . - /// The pivot indices of . - /// Not supported. This should be left null. - /// This is equivalent to the GETRI LAPACK routine. - [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)); - } - /// /// Solves A*X=B for X using LU factorization. /// @@ -538,60 +469,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda Solver(SafeNativeMethods.s_cholesky_solve_factored(_solverHandle, orderA, columnsB, a, b)); } - /// - /// Computes the singular value decomposition of A. - /// - /// Compute the singular U and VT vectors or not. - /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The singular values of A in ascending value. - /// If is true, on exit U contains the left - /// singular vectors. - /// If is true, on exit VT contains the transposed - /// right singular vectors. - /// This is equivalent to the GESVD LAPACK routine. - [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt) - { - 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); - } - /// /// Solves A*X=B for X using the singular value decomposition of A. /// @@ -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. /// If is true, on exit VT contains the transposed /// right singular vectors. - /// Not supported. Should be left null. /// This is equivalent to the GESVD LAPACK routine. [SecuritySafeCritical] - public override void SingularValueDecomposition(bool computeVectors, float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt, 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)); } } diff --git a/src/Numerics/Providers/LinearAlgebra/ILinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/ILinearAlgebraProvider.cs index 4b49da14..176ba9a3 100644 --- a/src/Numerics/Providers/LinearAlgebra/ILinearAlgebraProvider.cs +++ b/src/Numerics/Providers/LinearAlgebra/ILinearAlgebraProvider.cs @@ -269,29 +269,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// This is equivalent to the GETRI LAPACK routine. void LUInverseFactored(T[] a, int order, int[] ipiv); - /// - /// Computes the inverse of matrix using LU factorization. - /// - /// The N by N matrix to invert. Contains the inverse On exit. - /// The order of the square matrix . - /// 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. - /// This is equivalent to the GETRF and GETRI LAPACK routines. - void LUInverse(T[] a, int order, T[] work); - - /// - /// Computes the inverse of a previously factored matrix. - /// - /// The LU factored N by N matrix. Contains the inverse On exit. - /// The order of the square matrix . - /// The pivot indices of . - /// 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. - /// This is equivalent to the GETRI LAPACK routine. - void LUInverseFactored(T[] a, int order, int[] ipiv, T[] work); - /// /// Solves A*X=B for X using LU factorization. /// @@ -356,23 +333,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// This is similar to the GEQRF and ORGQR LAPACK routines. void QRFactor(T[] a, int rowsA, int columnsA, T[] q, T[] tau); - /// - /// Computes the full QR factorization of A. - /// - /// On entry, it is the M by N A matrix to factor. On exit, - /// it is overwritten with the R matrix of the QR factorization. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// On exit, A M by M matrix that holds the Q matrix of the - /// QR factorization. - /// A min(m,n) vector. On exit, contains additional information - /// to be used by the QR solve routine. - /// 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. - /// This is similar to the GEQRF and ORGQR LAPACK routines. - void QRFactor(T[] a, int rowsA, int columnsA, T[] q, T[] tau, T[] work); - /// /// Computes the thin QR factorization of A where M > N. /// @@ -387,23 +347,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// This is similar to the GEQRF and ORGQR LAPACK routines. void ThinQRFactor(T[] a, int rowsA, int columnsA, T[] r, T[] tau); - /// - /// Computes the thin QR factorization of A where M > N. - /// - /// On entry, it is the M by N A matrix to factor. On exit, - /// it is overwritten with the Q matrix of the QR factorization. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// On exit, A N by N matrix that holds the R matrix of the - /// QR factorization. - /// A min(m,n) vector. On exit, contains additional information - /// to be used by the QR solve routine. - /// 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. - /// This is similar to the GEQRF and ORGQR LAPACK routines. - void ThinQRFactor(T[] a, int rowsA, int columnsA, T[] r, T[] tau, T[] work); - /// /// Solves A*X=B for X using QR factorization of A. /// @@ -417,22 +360,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// Rows must be greater or equal to columns. void QRSolve(T[] a, int rows, int columns, T[] b, int columnsB, T[] x, QRMethod method = QRMethod.Full); - /// - /// Solves A*X=B for X using QR factorization of A. - /// - /// The A matrix. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// On entry the B matrix; on exit the X matrix. - /// The number of columns of B. - /// On exit, the solution matrix. - /// 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. - /// The type of QR factorization to perform. - /// Rows must be greater or equal to columns. - void QRSolve(T[] a, int rows, int columns, T[] b, int columnsB, T[] x, T[] work, QRMethod method = QRMethod.Full); - /// /// Solves A*X=B for X using a previously QR factored matrix. /// @@ -450,26 +377,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The type of QR factorization to perform. void QRSolveFactored(T[] q, T[] r, int rowsA, int columnsA, T[] tau, T[] b, int columnsB, T[] x, QRMethod method = QRMethod.Full); - /// - /// Solves A*X=B for X using a previously QR factored matrix. - /// - /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be - /// null for the native provider. The native provider uses the Q portion stored in the R matrix. - /// The R matrix obtained by calling . - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// Contains additional information on Q. Only used for the native solver - /// and can be null for the managed provider. - /// On entry the B matrix; on exit the X matrix. - /// The number of columns of B. - /// On exit, the solution matrix. - /// 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. - /// Rows must be greater or equal to columns. - /// The type of QR factorization to perform. - void QRSolveFactored(T[] q, T[] r, int rowsA, int columnsA, T[] tau, T[] b, int columnsB, T[] x, T[] work, QRMethod method = QRMethod.Full); - /// /// Computes the singular value decomposition of A. /// @@ -485,23 +392,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// This is equivalent to the GESVD LAPACK routine. void SingularValueDecomposition(bool computeVectors, T[] a, int rowsA, int columnsA, T[] s, T[] u, T[] vt); - /// - /// Computes the singular value decomposition of A. - /// - /// Compute the singular U and VT vectors or not. - /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The singular values of A in ascending value. - /// If is true, on exit U contains the left - /// singular vectors. - /// If is true, on exit VT contains the transposed - /// right singular vectors. - /// The work array. On exit, work[0] contains the optimal work size value. - /// - /// This is equivalent to the GESVD LAPACK routine. - void SingularValueDecomposition(bool computeVectors, T[] a, int rowsA, int columnsA, T[] s, T[] u, T[] vt, T[] work); - /// /// Solves A*X=B for X using the singular value decomposition of A. /// @@ -520,7 +410,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The number of columns in the A matrix. /// The s values returned by . /// The left singular vectors returned by . - /// The right singular vectors returned by . + /// The right singular vectors returned by . /// The B matrix /// The number of columns of B. /// On exit, the solution matrix. diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs index ce777bf1..a6a1fe37 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs +++ b/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); } - /// - /// Computes the inverse of matrix using LU factorization. - /// - /// The N by N matrix to invert. Contains the inverse On exit. - /// The order of the square matrix . - /// 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. - /// This is equivalent to the GETRF and GETRI LAPACK routines. - public virtual void LUInverse(Complex[] a, int order, Complex[] work) - { - LUInverse(a, order); - } - - /// - /// Computes the inverse of a previously factored matrix. - /// - /// The LU factored N by N matrix. Contains the inverse On exit. - /// The order of the square matrix . - /// The pivot indices of . - /// 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. - /// This is equivalent to the GETRI LAPACK routine. - public virtual void LUInverseFactored(Complex[] a, int order, int[] ipiv, Complex[] work) - { - LUInverseFactored(a, order, ipiv); - } - /// /// Solves A*X=B for X using LU factorization. /// @@ -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); - } - - /// - /// Computes the QR factorization of A. - /// - /// On entry, it is the M by N A matrix to factor. On exit, - /// it is overwritten with the R matrix of the QR factorization. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// On exit, A M by M matrix that holds the Q matrix of the - /// QR factorization. - /// A min(m,n) vector. On exit, contains additional information - /// to be used by the QR solve routine. - /// 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. - /// This is similar to the GEQRF and ORGQR LAPACK routines. - 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; } /// @@ -1599,61 +1502,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra } var work = new Complex[rowsA*columnsA]; - ThinQRFactor(a, rowsA, columnsA, r, tau, work); - } - - /// - /// Computes the QR factorization of A where M > N. - /// - /// On entry, it is the M by N A matrix to factor. On exit, - /// it is overwritten with the Q matrix of the QR factorization. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// On exit, A N by N matrix that holds the R matrix of the - /// QR factorization. - /// A min(m,n) vector. On exit, contains additional information - /// to be used by the QR solve routine. - /// 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. - /// This is similar to the GEQRF and ORGQR LAPACK routines. - 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 /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. public virtual void QRSolve(Complex[] a, int rows, int columns, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full) - { - var work = new Complex[rows*columns]; - QRSolve(a, rows, columns, b, columnsB, x, work, method); - } - - /// - /// Solves A*X=B for X using QR factorization of A. - /// - /// The A matrix. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The B matrix. - /// The number of columns of B. - /// On exit, the solution matrix. - /// The 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. - /// The type of QR factorization to perform. - /// Rows must be greater or equal to columns. - public virtual void QRSolve(Complex[] a, int rows, int columns, Complex[] b, int columnsB, Complex[] x, 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; - } - - /// - /// Solves A*X=B for X using a previously QR factored matrix. - /// - /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be - /// null for the native provider. The native provider uses the Q portion stored in the R matrix. - /// The R matrix obtained by calling . - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// Contains additional information on Q. Only used for the native solver - /// and can be null for the managed provider. - /// On entry the B matrix; on exit the X matrix. - /// The number of columns of B. - /// On exit, the solution matrix. - /// 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. - /// The type of QR factorization to perform. - /// Rows must be greater or equal to columns. - 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); } /// @@ -2102,77 +1896,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra } var work = new Complex[rowsA]; - SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work); - } - - /// - /// Computes the singular value decomposition of A. - /// - /// Compute the singular U and VT vectors or not. - /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The singular values of A in ascending value. - /// If is true, on exit U contains the left - /// singular vectors. - /// If is true, on exit VT contains the transposed - /// right singular vectors. - /// The work array. Length should be at least . - /// This is equivalent to the GESVD LAPACK routine. - /// - public virtual void SingularValueDecomposition(bool computeVectors, Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt, 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; } /// @@ -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); } diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs index fcb6a9e3..d68b9142 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs +++ b/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); } - /// - /// Computes the inverse of matrix using LU factorization. - /// - /// The N by N matrix to invert. Contains the inverse On exit. - /// The order of the square matrix . - /// 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. - /// This is equivalent to the GETRF and GETRI LAPACK routines. - public virtual void LUInverse(Complex32[] a, int order, Complex32[] work) - { - LUInverse(a, order); - } - - /// - /// Computes the inverse of a previously factored matrix. - /// - /// The LU factored N by N matrix. Contains the inverse On exit. - /// The order of the square matrix . - /// The pivot indices of . - /// 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. - /// This is equivalent to the GETRI LAPACK routine. - public virtual void LUInverseFactored(Complex32[] a, int order, int[] ipiv, Complex32[] work) - { - LUInverseFactored(a, order, ipiv); - } - /// /// Solves A*X=B for X using LU factorization. /// @@ -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); - } - - /// - /// Computes the QR factorization of A. - /// - /// On entry, it is the M by N A matrix to factor. On exit, - /// it is overwritten with the R matrix of the QR factorization. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// On exit, A M by M matrix that holds the Q matrix of the - /// QR factorization. - /// A min(m,n) vector. On exit, contains additional information - /// to be used by the QR solve routine. - /// 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. - /// This is similar to the GEQRF and ORGQR LAPACK routines. - 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; } /// @@ -1596,62 +1499,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra } var work = new Complex32[rowsA*columnsA]; - ThinQRFactor(a, rowsA, columnsA, r, tau, work); - } - - /// - /// Computes the QR factorization of A where M > N. - /// - /// On entry, it is the M by N A matrix to factor. On exit, - /// it is overwritten with the Q matrix of the QR factorization. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// On exit, A N by N matrix that holds the R matrix of the - /// QR factorization. - /// A min(m,n) vector. On exit, contains additional information - /// to be used by the QR solve routine. - /// 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. - /// This is similar to the GEQRF and ORGQR LAPACK routines. - 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 /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. public virtual void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x, QRMethod method = QRMethod.Full) - { - var work = new Complex32[rows*columns]; - QRSolve(a, rows, columns, b, columnsB, x, work, method); - } - - /// - /// Solves A*X=B for X using QR factorization of A. - /// - /// The A matrix. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The B matrix. - /// The number of columns of B. - /// On exit, the solution matrix. - /// The 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. - /// The type of QR factorization to perform. - /// Rows must be greater or equal to columns. - public virtual void QRSolve(Complex32[] a, int rows, int columns, Complex32[] b, int columnsB, Complex32[] x, 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; - } - - /// - /// Solves A*X=B for X using a previously QR factored matrix. - /// - /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be - /// null for the native provider. The native provider uses the Q portion stored in the R matrix. - /// The R matrix obtained by calling . - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// Contains additional information on Q. Only used for the native solver - /// and can be null for the managed provider. - /// On entry the B matrix; on exit the X matrix. - /// The number of columns of B. - /// On exit, the solution matrix. - /// 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. - /// The type of QR factorization to perform. - /// Rows must be greater or equal to columns. - 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); } /// @@ -2099,77 +1891,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra } var work = new Complex32[rowsA]; - SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work); - } - - /// - /// Computes the singular value decomposition of A. - /// - /// Compute the singular U and VT vectors or not. - /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The singular values of A in ascending value. - /// If is true, on exit U contains the left - /// singular vectors. - /// If is true, on exit VT contains the transposed - /// right singular vectors. - /// The work array. Length should be at least . - /// This is equivalent to the GESVD LAPACK routine. - /// - public virtual void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt, 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; } /// @@ -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); } diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs index c2adea46..d3462126 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs +++ b/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); } - /// - /// Computes the inverse of matrix using LU factorization. - /// - /// The N by N matrix to invert. Contains the inverse On exit. - /// The order of the square matrix . - /// 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. - /// This is equivalent to the GETRF and GETRI LAPACK routines. - public virtual void LUInverse(double[] a, int order, double[] work) - { - LUInverse(a, order); - } - - /// - /// Computes the inverse of a previously factored matrix. - /// - /// The LU factored N by N matrix. Contains the inverse On exit. - /// The order of the square matrix . - /// The pivot indices of . - /// 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. - /// This is equivalent to the GETRI LAPACK routine. - public virtual void LUInverseFactored(double[] a, int order, int[] ipiv, double[] work) - { - LUInverseFactored(a, order, ipiv); - } - /// /// Solves A*X=B for X using LU factorization. /// @@ -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); - } - - /// - /// Computes the QR factorization of A. - /// - /// On entry, it is the M by N A matrix to factor. On exit, - /// it is overwritten with the R matrix of the QR factorization. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// On exit, A M by M matrix that holds the Q matrix of the - /// QR factorization. - /// A min(m,n) vector. On exit, contains additional information - /// to be used by the QR solve routine. - /// 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. - /// This is similar to the GEQRF and ORGQR LAPACK routines. - 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; } /// @@ -1483,62 +1384,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra } var work = new double[rowsA*columnsA]; - ThinQRFactor(a, rowsA, columnsA, r, tau, work); - } - - /// - /// Computes the thin QR factorization of A where M > N. - /// - /// On entry, it is the M by N A matrix to factor. On exit, - /// it is overwritten with the Q matrix of the QR factorization. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// On exit, A N by N matrix that holds the R matrix of the - /// QR factorization. - /// A min(m,n) vector. On exit, contains additional information - /// to be used by the QR solve routine. - /// 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. - /// This is similar to the GEQRF and ORGQR LAPACK routines. - 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 /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. public virtual void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full) - { - var work = new double[rows*columns]; - QRSolve(a, rows, columns, b, columnsB, x, work, method); - } - - /// - /// Solves A*X=B for X using QR factorization of A. - /// - /// The A matrix. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The B matrix. - /// The number of columns of B. - /// On exit, the solution matrix. - /// The 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. - /// The type of QR factorization to perform. - /// Rows must be greater or equal to columns. - public virtual void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x, 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; - } - - /// - /// Solves A*X=B for X using a previously QR factored matrix. - /// - /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be - /// null for the native provider. The native provider uses the Q portion stored in the R matrix. - /// The R matrix obtained by calling . - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// Contains additional information on Q. Only used for the native solver - /// and can be null for the managed provider. - /// On entry the B matrix; on exit the X matrix. - /// The number of columns of B. - /// On exit, the solution matrix. - /// 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. - /// The type of QR factorization to perform. - /// Rows must be greater or equal to columns. - public virtual void QRSolveFactored(double[] q, double[] r, int rowsA, int columnsA, double[] tau, double[] b, int columnsB, double[] x, double[] work, QRMethod method = QRMethod.Full) - { - QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, method); } /// @@ -1986,76 +1776,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra } var work = new double[rowsA]; - SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work); - } - - /// - /// Computes the singular value decomposition of A. - /// - /// Compute the singular U and VT vectors or not. - /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The singular values of A in ascending value. - /// If is true, on exit U contains the left - /// singular vectors. - /// If is true, on exit VT contains the transposed - /// right singular vectors. - /// The work array. Length should be at least . - /// This is equivalent to the GESVD LAPACK routine. - /// - public virtual void SingularValueDecomposition(bool computeVectors, double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt, 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; } /// @@ -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); } diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs index 7ed95c45..05f7989c 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs +++ b/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); } - /// - /// Computes the inverse of matrix using LU factorization. - /// - /// The N by N matrix to invert. Contains the inverse On exit. - /// The order of the square matrix . - /// 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. - /// This is equivalent to the GETRF and GETRI LAPACK routines. - public virtual void LUInverse(float[] a, int order, float[] work) - { - LUInverse(a, order); - } - - /// - /// Computes the inverse of a previously factored matrix. - /// - /// The LU factored N by N matrix. Contains the inverse On exit. - /// The order of the square matrix . - /// The pivot indices of . - /// 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. - /// This is equivalent to the GETRI LAPACK routine. - public virtual void LUInverseFactored(float[] a, int order, int[] ipiv, float[] work) - { - LUInverseFactored(a, order, ipiv); - } - /// /// Solves A*X=B for X using LU factorization. /// @@ -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); - } - - /// - /// Computes the QR factorization of A. - /// - /// On entry, it is the M by N A matrix to factor. On exit, - /// it is overwritten with the R matrix of the QR factorization. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// On exit, A M by M matrix that holds the Q matrix of the - /// QR factorization. - /// A min(m,n) vector. On exit, contains additional information - /// to be used by the QR solve routine. - /// 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. - /// This is similar to the GEQRF and ORGQR LAPACK routines. - 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; } /// @@ -1482,61 +1385,6 @@ namespace MathNet.Numerics.Providers.LinearAlgebra } var work = new float[rowsA*columnsA]; - ThinQRFactor(a, rowsA, columnsA, r, tau, work); - } - - /// - /// Computes the QR factorization of A where M > N. - /// - /// On entry, it is the M by N A matrix to factor. On exit, - /// it is overwritten with the Q matrix of the QR factorization. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// On exit, A N by N matrix that holds the R matrix of the - /// QR factorization. - /// A min(m,n) vector. On exit, contains additional information - /// to be used by the QR solve routine. - /// 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. - /// This is similar to the GEQRF and ORGQR LAPACK routines. - 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 - /// - /// Solves A*X=B for X using QR factorization of A. - /// - /// The A matrix. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The B matrix. - /// The number of columns of B. - /// On exit, the solution matrix. - /// The type of QR factorization to perform. - /// Rows must be greater or equal to columns. - public 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); - } /// /// Solves A*X=B for X using QR factorization of A. @@ -1713,12 +1543,9 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. - /// 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. /// The type of QR factorization to perform. /// Rows must be greater or equal to columns. - public virtual void QRSolve(float[] a, int rows, int columns, float[] b, int columnsB, float[] x, 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; - } - - /// - /// Solves A*X=B for X using a previously QR factored matrix. - /// - /// The Q matrix obtained by QR factor. This is only used for the managed provider and can be - /// null for the native provider. The native provider uses the Q portion stored in the R matrix. - /// The R matrix obtained by calling . - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// Contains additional information on Q. Only used for the native solver - /// and can be null for the managed provider. - /// On entry the B matrix; on exit the X matrix. - /// The number of columns of B. - /// On exit, the solution matrix. - /// 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. - /// The type of QR factorization to perform. - /// Rows must be greater or equal to columns. - public virtual void QRSolveFactored(float[] q, float[] r, int rowsA, int columnsA, float[] tau, float[] b, int columnsB, float[] x, float[] work, QRMethod method = QRMethod.Full) - { - QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, method); } /// @@ -1986,76 +1779,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra } var work = new float[rowsA]; - SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, work); - } - - /// - /// Computes the singular value decomposition of A. - /// - /// Compute the singular U and VT vectors or not. - /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. - /// The number of rows in the A matrix. - /// The number of columns in the A matrix. - /// The singular values of A in ascending value. - /// If is true, on exit U contains the left - /// singular vectors. - /// If is true, on exit VT contains the transposed - /// right singular vectors. - /// The work array. Length should be at least . - /// - 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; } /// @@ -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); } diff --git a/src/UnitTests/LinearAlgebraProviderTests/Complex/LinearAlgebraProviderTests.cs b/src/UnitTests/LinearAlgebraProviderTests/Complex/LinearAlgebraProviderTests.cs index 3998ae08..bc60a10e 100644 --- a/src/UnitTests/LinearAlgebraProviderTests/Complex/LinearAlgebraProviderTests.cs +++ b/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 - /// - /// Can compute the inverse of a matrix using LU factorization - /// with a work array. - /// - [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 - /// - /// Can compute the inverse of a matrix using LU factorization - /// using a previously factored matrix with a work array. - /// - [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 - /// /// Can solve Ax=b using LU factorization. /// @@ -696,92 +630,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex } } -#if ! MKL - /// - /// Can compute QR factorization of a square matrix using a work array. - /// - [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); - } - } - } - - /// - /// Can compute QR factorization of a tall matrix using a work matrix. - /// - [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); - } - } - } - - /// - /// Can compute QR factorization of a wide matrix using a work matrix. - /// - [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 - /// /// Can compute thin QR factorization of a square matrix. /// @@ -836,63 +684,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex } } -#if ! MKL - /// - /// Can compute thin QR factorization of a square matrix using a work array. - /// - [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); - } - } - } - - /// - /// Can compute thin QR factorization of a tall matrix using a work matrix. - /// - [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 - /// /// Can solve Ax=b using QR factorization with a square A matrix. /// @@ -944,94 +735,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 13); } - -#if ! MKL - /// - /// Can solve Ax=b using QR factorization with a square A matrix - /// using a work array. - /// - [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); - } - - /// - /// Can solve Ax=b using QR factorization with a tall A matrix - /// using a work array. - /// - [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); - } - - /// - /// Can solve Ax=b using QR factorization with a square A matrix - /// using a factored A matrix. - /// - [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 /// /// 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 - /// - /// Can solve Ax=b using QR factorization with a square A matrix - /// using a factored A matrix with a work array. - /// - [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); - } - - /// - /// Can solve Ax=b using QR factorization with a tall A matrix - /// using a factored A matrix with a work array. - /// - [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 - /// /// Can solve Ax=b using thin QR factorization with a square A matrix. /// @@ -1174,64 +815,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 13); } -#if ! MKL - /// - /// Can solve Ax=b using thin QR factorization with a square A matrix - /// using a work array. - /// - [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); - } - - /// - /// Can solve Ax=b using thin QR factorization with a tall A matrix - /// using a work array. - /// - [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 - /// /// 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); } - /// - /// Can solve Ax=b using thin QR factorization with a square A matrix - /// using a factored A matrix with a work array. - /// - [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); - } - - /// - /// Can solve Ax=b using thin QR factorization with a tall A matrix - /// using a factored A matrix with a work array. - /// - [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); - } - /// /// Can compute the SVD factorization of a square matrix. /// @@ -1455,129 +978,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex AssertHelpers.AlmostEqualRelative(matrix[1, 2], result[1, 2], 14); } - /// - /// Can compute the SVD factorization of a square matrix using - /// a work array. - /// - [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); - } - - /// - /// Can compute the SVD factorization of a tall matrix using - /// a work array. - /// - [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); - } - - /// - /// Can compute the SVD factorization of a wide matrix using - /// a work array. - /// - [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); - } - /// /// Can solve Ax=b using SVD factorization with a square A matrix. /// diff --git a/src/UnitTests/LinearAlgebraProviderTests/Complex32/LinearAlgebraProviderTests.cs b/src/UnitTests/LinearAlgebraProviderTests/Complex32/LinearAlgebraProviderTests.cs index d0b8a4e3..1b597d22 100644 --- a/src/UnitTests/LinearAlgebraProviderTests/Complex32/LinearAlgebraProviderTests.cs +++ b/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); } - /// - /// Can compute the inverse of a matrix using LU factorization - /// with a work array. - /// - [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); - } - - /// - /// Can compute the inverse of a matrix using LU factorization - /// using a previously factored matrix with a work array. - /// - [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); - } - /// /// Can solve Ax=b using LU factorization. /// @@ -696,90 +634,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex32 } } - /// - /// Can compute QR factorization of a square matrix using a work array. - /// - [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); - } - } - } - - /// - /// Can compute QR factorization of a tall matrix using a work matrix. - /// - [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); - } - } - } - - /// - /// Can compute QR factorization of a wide matrix using a work matrix. - /// - [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); - } - } - } - /// /// Can compute thin QR factorization of a square matrix. /// @@ -834,61 +688,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex32 } } - /// - /// Can compute thin QR factorization of a square matrix using a work array. - /// - [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); - } - } - } - - /// - /// Can compute thin QR factorization of a tall matrix using a work matrix. - /// - [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); - } - } - } - /// /// Can solve Ax=b using QR factorization with a square A matrix. /// @@ -941,62 +740,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex32 AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 5); } - /// - /// Can solve Ax=b using QR factorization with a square A matrix - /// using a work array. - /// - [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); - } - - /// - /// Can solve Ax=b using QR factorization with a tall A matrix - /// using a work array. - /// - [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); - } - /// /// 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); } - /// - /// Can solve Ax=b using QR factorization with a square A matrix - /// using a factored A matrix with a work array. - /// - [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); - } - - /// - /// Can solve Ax=b using QR factorization with a tall A matrix - /// using a factored A matrix with a work array. - /// - [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); - } - /// /// Can solve Ax=b using thin QR factorization with a square A matrix. /// @@ -1167,62 +850,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex32 AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 5); } - /// - /// Can solve Ax=b using thin QR factorization with a square A matrix - /// using a work array. - /// - [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); - } - - /// - /// Can solve Ax=b using thin QR factorization with a tall A matrix - /// using a work array. - /// - [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); - } - /// /// 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); } - /// - /// Can solve Ax=b using thin QR factorization with a square A matrix - /// using a factored A matrix with a work array. - /// - [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); - } - - /// - /// Can solve Ax=b using thin QR factorization with a tall A matrix - /// using a factored A matrix with a work array. - /// - [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); - } - /// /// Can compute the SVD factorization of a square matrix. /// @@ -1446,129 +1013,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex32 AssertHelpers.AlmostEqualRelative(matrix[1, 2], result[1, 2], 4); } - /// - /// Can compute the SVD factorization of a square matrix using - /// a work array. - /// - [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); - } - - /// - /// Can compute the SVD factorization of a tall matrix using - /// a work array. - /// - [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); - } - - /// - /// Can compute the SVD factorization of a wide matrix using - /// a work array. - /// - [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); - } - /// /// Can solve Ax=b using SVD factorization with a square A matrix. /// diff --git a/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs b/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs index 6546bc55..a232431f 100644 --- a/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs +++ b/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); } - /// - /// Can compute the inverse of a matrix using LU factorization - /// with a work array. - /// - [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); - } - - /// - /// Can compute the inverse of a matrix using LU factorization - /// using a previously factored matrix with a work array. - /// - [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); - } - /// /// Can solve Ax=b using LU factorization. /// @@ -686,90 +624,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double } } - /// - /// Can compute QR factorization of a square matrix using a work array. - /// - [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.Build.Dense(matrix.RowCount, matrix.RowCount, q); - var mr = Matrix.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); - } - } - } - - /// - /// Can compute QR factorization of a tall matrix using a work matrix. - /// - [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.Build.Dense(matrix.RowCount, matrix.ColumnCount, r).UpperTriangle(); - var mq = Matrix.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); - } - } - } - - /// - /// Can compute QR factorization of a wide matrix using a work matrix. - /// - [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.Build.Dense(matrix.RowCount, matrix.ColumnCount, r).UpperTriangle(); - var mq = Matrix.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); - } - } - } - /// /// Can compute thin QR factorization of a square matrix. /// @@ -824,61 +678,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double } } - /// - /// Can compute thin QR factorization of a square matrix using a work array. - /// - [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.Build.Dense(matrix.RowCount, matrix.ColumnCount, q); - var mr = Matrix.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); - } - } - } - - /// - /// Can compute thin QR factorization of a tall matrix using a work matrix. - /// - [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.Build.Dense(matrix.RowCount, matrix.ColumnCount, q); - var mr = Matrix.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); - } - } - } - /// /// Can solve Ax=b using QR factorization with a square A matrix. /// @@ -931,62 +730,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 13); } - /// - /// Can solve Ax=b using QR factorization with a square A matrix - /// using a work array. - /// - [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.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); - } - - /// - /// Can solve Ax=b using QR factorization with a tall A matrix - /// using a work array. - /// - [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.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); - } - /// /// 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); } - /// - /// Can solve Ax=b using QR factorization with a square A matrix - /// using a factored A matrix with a work array. - /// - [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.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); - } - - /// - /// Can solve Ax=b using QR factorization with a tall A matrix - /// using a factored A matrix with a work array. - /// - [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.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); - } - /// /// Can solve Ax=b using thin QR factorization with a square A matrix. /// @@ -1157,62 +840,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 13); } - /// - /// Can solve Ax=b using thin QR factorization with a square A matrix - /// using a work array. - /// - [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.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); - } - - /// - /// Can solve Ax=b using thin QR factorization with a tall A matrix - /// using a work array. - /// - [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.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); - } - /// /// 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); } - /// - /// Can solve Ax=b using thin QR factorization with a square A matrix - /// using a factored A matrix with a work array. - /// - [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.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); - } - - /// - /// Can solve Ax=b using thin QR factorization with a tall A matrix - /// using a factored A matrix with a work array. - /// - [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.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); - } - /// /// Can compute the SVD factorization of a square matrix. /// @@ -1436,129 +1003,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double AssertHelpers.AlmostEqualRelative(matrix[1, 2], result[1, 2], 14); } - /// - /// Can compute the SVD factorization of a square matrix using - /// a work array. - /// - [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.Build.Dense(matrix.RowCount, matrix.RowCount, u); - var mV = Matrix.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); - } - - /// - /// Can compute the SVD factorization of a tall matrix using - /// a work array. - /// - [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.Build.Dense(matrix.RowCount, matrix.RowCount, u); - var mV = Matrix.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); - } - - /// - /// Can compute the SVD factorization of a wide matrix using - /// a work array. - /// - [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.Build.Dense(matrix.RowCount, matrix.RowCount, u); - var mV = Matrix.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); - } - /// /// Can solve Ax=b using SVD factorization with a square A matrix. /// diff --git a/src/UnitTests/LinearAlgebraProviderTests/Single/LinearAlgebraProviderTests.cs b/src/UnitTests/LinearAlgebraProviderTests/Single/LinearAlgebraProviderTests.cs index fab970db..a64a6973 100644 --- a/src/UnitTests/LinearAlgebraProviderTests/Single/LinearAlgebraProviderTests.cs +++ b/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); } - /// - /// Can compute the inverse of a matrix using LU factorization - /// with a work array. - /// - [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); - } - - /// - /// Can compute the inverse of a matrix using LU factorization - /// using a previously factored matrix with a work array. - /// - [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); - } - /// /// Can solve Ax=b using LU factorization. /// @@ -694,90 +632,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Single } } - /// - /// Can compute QR factorization of a square matrix using a work array. - /// - [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); - } - } - } - - /// - /// Can compute QR factorization of a tall matrix using a work matrix. - /// - [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); - } - } - } - - /// - /// Can compute QR factorization of a wide matrix using a work matrix. - /// - [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); - } - } - } - /// /// Can compute thin QR factorization of a square matrix. /// @@ -832,61 +686,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Single } } - /// - /// Can compute thin QR factorization of a square matrix using a work array. - /// - [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); - } - } - } - - /// - /// Can compute thin QR factorization of a tall matrix using a work matrix. - /// - [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); - } - } - } - /// /// Can solve Ax=b using QR factorization with a square A matrix. /// @@ -939,62 +738,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Single AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 5); } - /// - /// Can solve Ax=b using QR factorization with a square A matrix - /// using a work array. - /// - [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); - } - - /// - /// Can solve Ax=b using QR factorization with a tall A matrix - /// using a work array. - /// - [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); - } - /// /// 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); } - /// - /// Can solve Ax=b using QR factorization with a square A matrix - /// using a factored A matrix with a work array. - /// - [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); - } - - /// - /// Can solve Ax=b using QR factorization with a tall A matrix - /// using a factored A matrix with a work array. - /// - [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); - } - /// /// Can solve Ax=b using thin QR factorization with a square A matrix. /// @@ -1165,62 +848,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Single AssertHelpers.AlmostEqualRelative(test[1, 1], x[3], 5); } - /// - /// Can solve Ax=b using thin QR factorization with a square A matrix - /// using a work array. - /// - [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); - } - - /// - /// Can solve Ax=b using thin QR factorization with a tall A matrix - /// using a work array. - /// - [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); - } - /// /// 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); } - /// - /// Can solve Ax=b using thin QR factorization with a square A matrix - /// using a factored A matrix with a work array. - /// - [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); - } - - /// - /// Can solve Ax=b using thin QR factorization with a tall A matrix - /// using a factored A matrix with a work array. - /// - [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); - } - /// /// Can compute the SVD factorization of a square matrix. /// @@ -1444,129 +1011,6 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Single AssertHelpers.AlmostEqualRelative(matrix[1, 2], result[1, 2], 5); } - /// - /// Can compute the SVD factorization of a square matrix using - /// a work array. - /// - [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); - } - - /// - /// Can compute the SVD factorization of a tall matrix using - /// a work array. - /// - [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); - } - - /// - /// Can compute the SVD factorization of a wide matrix using - /// a work array. - /// - [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); - } - /// /// Can solve Ax=b using SVD factorization with a square A matrix. /// diff --git a/src/UnitTests/UnitTests-CUDA.csproj b/src/UnitTests/UnitTests-CUDA.csproj index 074df81c..291c5d99 100644 --- a/src/UnitTests/UnitTests-CUDA.csproj +++ b/src/UnitTests/UnitTests-CUDA.csproj @@ -336,6 +336,10 @@ {b7cae5f4-a23f-4438-b5be-41226618b695} Numerics + + {af3253c9-4db5-45a0-98cf-c105fda9da47} + TestData + diff --git a/src/UnitTests/UnitTests-MKL.csproj b/src/UnitTests/UnitTests-MKL.csproj index f0972c67..024b8d9a 100644 --- a/src/UnitTests/UnitTests-MKL.csproj +++ b/src/UnitTests/UnitTests-MKL.csproj @@ -130,6 +130,10 @@ + + + + @@ -652,6 +656,10 @@ {b7cae5f4-a23f-4438-b5be-41226618b695} Numerics + + {af3253c9-4db5-45a0-98cf-c105fda9da47} + TestData + @@ -664,10 +672,5 @@ - - - - - - + \ No newline at end of file diff --git a/src/UnitTests/UnitTests-OpenBLAS.csproj b/src/UnitTests/UnitTests-OpenBLAS.csproj index 94b0ac3e..7cad1e62 100644 --- a/src/UnitTests/UnitTests-OpenBLAS.csproj +++ b/src/UnitTests/UnitTests-OpenBLAS.csproj @@ -336,6 +336,10 @@ {b7cae5f4-a23f-4438-b5be-41226618b695} Numerics + + {af3253c9-4db5-45a0-98cf-c105fda9da47} + TestData +