From 2503155bc7ba706ba2829598222952c0f60468cd Mon Sep 17 00:00:00 2001 From: Marcus Cuda Date: Wed, 16 Sep 2009 19:48:05 +0800 Subject: [PATCH] flushed out linear algebra provider interface --- .../LinearAlgebra/ILinearAlgebraProvider.cs | 365 +++++++++++++++++- .../ManagedLinearAlgebraProvider.cs | 144 +++++++ 2 files changed, 508 insertions(+), 1 deletion(-) diff --git a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProvider.cs index 2875df6a..50504a29 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProvider.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProvider.cs @@ -1,4 +1,4 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://mathnet.opensourcedotnet.info // Copyright (c) 2009 Math.NET @@ -22,13 +22,71 @@ // OTHER DEALINGS IN THE SOFTWARE. // +// INITIAL DRAFT MISSING EXCEPTION SPECIFICATIONS namespace MathNet.Numerics.Algorithms.LinearAlgebra { + /// + /// How to transpose a matrix. + /// + public enum Transpose + { + /// + /// Don't transpose a matrix. + /// + DontTranspose = 111, + + /// + /// Transpose a matrix. + /// + Transpose = 112, + + /// + /// Conjugate transpose a complex matrix. + /// + /// If a conjugate transpose is used with a real matrix, then the matrix is just transposed. + ConjugateTranspose = 113 + } + + /// + /// Types of matrix norms. + /// + public enum Norm : byte + { + /// + /// The 1-norm. + /// + OneNorm = (byte)'1', + + /// + /// The Frobenius norm. + /// + FrobeniusNorm = (byte)'f', + + /// + /// The infinity norm. + /// + InfinityNorm = (byte)'i', + + /// + /// The largest absolute value norm. + /// + LargestAbsoluteValue = (byte)'m' + } + /// /// Interface to linear algebra algorithms that work off 1-D arrays. /// public interface ILinearAlgebraProvider { + /// + /// Queries the provider for the optimal, workspace block size + /// for the given routine. + /// + /// Name of the method to query. + /// -1 if the provider cannot compute the workspace size; otherwise + /// the suggested block size. + int QueryWorkspaceBlockSize(string methodName); + /// /// Adds a scaled vector to another: y += alpha*x. /// @@ -45,5 +103,310 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The values to scale. /// This is equivalent to the SCAL BLAS routine. void ScaleArray(double alpha, double[] x); + + /// + /// Computes the dot product of x and y. + /// + /// The vector x. + /// The vector y. + /// The dot product of x and y. + /// This is equivalent to the DOT BLAS routine. + double DotProduct(double[] x, double[] y); + + /// + /// Does a point wise add of two arrays z = x + y. This can be used + /// to add vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the addition. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + void AddArrays(double[] x, double[] y, double[] result); + + /// + /// Does a point wise subtraction of two arrays z = x - y. This can be used + /// to subtract vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the subtraction. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + void SubtractArrays(double[] x, double[] y, double[] result); + + /// + /// Does a point wise multiplication of two arrays z = x * y. This can be used + /// to multiple elements of vectors or matrices. + /// + /// The array x. + /// The array y. + /// The result of the point wise multiplication. + /// There is no equivalent BLAS routine, but many libraries + /// provide optimized (parallel and/or vectorized) versions of this + /// routine. + void PointWiseMultiplyArrays(double[] x, double[] y, double[] result); + + /// + /// Computes the requested of the matrix. + /// + /// The type of norm to compute. + /// The matrix to compute the norm from. + /// + /// The requested of the matrix. + /// + double MatrixNorm(Norm norm, double[] matrix); + + /// + /// Computes the requested of the matrix. + /// + /// The type of norm to compute. + /// The matrix to compute the norm from. + /// The work array. Only used when + /// and needs to be have a length of at least M (number of rows of . + /// + /// The requested of the matrix. + /// + double MatrixNorm(Norm norm, double[] matrix, double[] work); + + /// + /// Multiples two matrices. result = x * y + /// + /// The x matrix. + /// The y matrix. + /// Where to store the result of the multiplication. + /// This is a simplified version of the BLAS GEMM routine with alpha + /// set to 1.0 and beta set to 0.0, and x and y are not transposed. + void MatrixMultiply(double[] x, double[] y, double[] result); + + /// + /// Multiplies two matrices and updates another with the result. c = alpha*op(a)*op(b) + beta*c + /// + /// How to transpose the matrix. + /// How to transpose the matrix. + /// The value to scale matrix. + /// The a matrix. + /// The b matrix + /// The value to scale the matrix. + /// The c matrix. + void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, double[] b, double beta, double[] c); + + /// + /// Computes the LU factorization of A. + /// + /// An m by n matrix. The matrix is overwritten with the + /// the LU factorization On exit. + /// On exit, it contains the pivot indices. The size + /// of the array must be min(m,n). + /// This is equivalent to the GETRF LAPACK routine. + void LUFactor(double[] a, int[] ipiv); + + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// This is equivalent to the GETRF and GETRI LAPACK routines. + void LUInverse(double[] a); + + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The pivot indices of . + /// This is equivalent to the GETRI LAPACK routine. + void LUInverseFactored(double[] a, int[] ipiv); + + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The work array. The array must have a length of at least N, + /// but should be N*blocksize. The blocksize is machine dependent. Use + /// to determine the optimal size of the work array. On exit, work[0] contains the optimal + /// work size value. + /// This is equivalent to the GETRF and GETRI LAPACK routines. + void LUInverse(double[] a, double[] work); + + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// 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. Use + /// to determine the optimal size of the work array. On exit, work[0] contains the optimal + /// work size value. + /// This is equivalent to the GETRI LAPACK routine. + void LUInverseFactored(double[] a, int[] ipiv, double[] work); + + /// + /// Solves A*X=B for X using LU factorization. + /// + /// The number of columns of B. + /// The square matrix A. + /// The B matrix. + /// This is equivalent to the GETRF and GETRS LAPACK routines. + void LUSolve(int columnsOfB, double[] a, double[] b); + + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// The number of columns of B. + /// The factored A matrix. + /// The pivot indices of . + /// The B matrix. + /// This is equivalent to the GETRS LAPACK routine. + void LUSolveFactored(int columnsOfB, double[] a, int ipiv, double[] b); + + /// + /// Solves A*X=B for X using LU factorization. + /// + /// How to transpose the matrix. + /// The number of columns of B. + /// The square matrix A. + /// The B matrix. + /// This is equivalent to the GETRF and GETRS LAPACK routines. + void LUSolve(Transpose transposeA, int columnsOfB, double[] a, double[] b); + + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// How to transpose the matrix. + /// The number of columns of B. + /// The factored A matrix. + /// The pivot indices of . + /// The B matrix. + /// This is equivalent to the GETRS LAPACK routine. + void LUSolveFactored(Transpose transposeA, int columnsOfB, double[] a, int ipiv, double[] b); + + /// + /// Computes the Cholesky factorization of A. + /// + /// A square, positive definite matrix. The matrix is overwritten with the + /// the Cholesky factorization On exit. + /// This is equivalent to the POTRF LAPACK routine. + void CholeskyFactor(double[] a); + + /// + /// Solves A*X=B for X using Cholesky factorization. + /// + /// The number of columns of B. + /// The square, positive definite matrix A. + /// The B matrix. + /// This is equivalent to the POTRF add POTRS LAPACK routines. + void CholeskySolve(int columnsOfB, double[] a, double[] b); + + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// The number of columns of B. + /// The factored A matrix. + /// The B matrix. + /// This is equivalent to the POTRS LAPACK routine. + void CholeskySolveFactored(int columnsOfB, double[] a, double[] b); + + /// + /// Computes the QR factorization of A. + /// + /// On entry, it is the M by N A matrix to factor. On exit, + /// the elements on and above the diagonal of the array + /// contain the min(M,N)-by-N upper trapezoidal matrix R (R is + /// upper triangular if m >= n); the elements below the diagonal, + /// with the array , represent the orthogonal matrix Q as a + /// product of min(m,n) elementary reflectors. + /// On exit, A M by M matrix that holds the Q matrix of the + /// QR factorization. + /// On exit, tau contains information needed by the + /// method. + /// This is equivalent to the GEQRF and ORGQR LAPACK routines. + void QRFactor(double[] a, double[] q, double[] tau); + + /// + /// Computes the QR factorization of A. + /// + /// On entry, it is the M by N A matrix to factor. On exit, + /// the elements on and above the diagonal of the array + /// contain the min(M,N)-by-N upper trapezoidal matrix R (R is + /// upper triangular if m >= n); the elements below the diagonal, + /// with the array , represent the orthogonal matrix Q as a + /// product of min(m,n) elementary reflectors. + /// On exit, A M by M matrix that holds the Q matrix of the + /// QR factorization. + /// On exit, tau contains information needed by the + /// method. + /// The work array. The array must have a length of at least N, + /// but should be N*blocksize. The blocksize is machine dependent. Use + /// to determine the optimal size of the work array. On exit, work[0] contains the optimal + /// work size value. + /// This is equivalent to the GEQRF and ORGQR LAPACK routines. + void QRFactor(double[] a, double[] q, double[] tau, double[] work); + + /// + /// Solves A*X=B for X using a previously QR factored matrix. + /// + /// The number of columns of B. + /// The matrix obtained by calling . + /// The tau vector obtained by calling . + /// The B matrix. + /// On exit, the solution matrix. + /// This is equivalent to the ORMQR LAPACK routine with the TRSM BLAS routine. + void QRSolve(int columnsOfB, double[] a, double[] tau, double[] b, double[] x); + + /// + /// Solves A*X=B for X using a previously QR factored matrix. + /// + /// The number of columns of B. + /// The M by N + /// The tau vector obtained by calling . + /// The B matrix. + /// 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. Use + /// to determine the optimal size of the work array. On exit, work[0] contains the optimal + /// work size value. + /// This is equivalent to the ORMQR LAPACK routine with the TRSM BLAS routine. + void QRSolve(int columnsOfB, double[] a, double[] tau, double[] b, double[] x, double[] 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 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. + void SinguarValueDecomposition(bool computeVectors, double[] a, double[] s, double[] u, double[] 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 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. The array must have a length of at least N, + /// but should be N*blocksize. The blocksize is machine dependent. Use + /// to determine the optimal size of 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, double[] a, double[] s, double[] u, double[] vt, double[] work); + + /// + /// Solves A*X=B for X using a previously SVD decomposed matrix. + /// + /// The s values returned by . + /// The left singular vectors returned by . + /// The right singular vectors returned by . + /// The B matrix. + /// On exit, the solution matrix. + void SvdSolve(double[] s, double[] u, double[] vt, double[] b, double[] x); } } diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs index 258cd7c1..125711b2 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs @@ -85,5 +85,149 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra Parallel.For(0, x.Length, i => x[i] = alpha * x[i]); } + + #region ILinearAlgebraProvider Members + + public int QueryWorkspaceBlockSize(string methodName) + { + throw new NotImplementedException(); + } + + public double DotProduct(double[] x, double[] y) + { + throw new NotImplementedException(); + } + + public void AddArrays(double[] x, double[] y, double[] result) + { + throw new NotImplementedException(); + } + + public void SubtractArrays(double[] x, double[] y, double[] result) + { + throw new NotImplementedException(); + } + + public void PointWiseMultiplyArrays(double[] x, double[] y, double[] result) + { + throw new NotImplementedException(); + } + + public double MatrixNorm(Norm norm, double[] matrix) + { + throw new NotImplementedException(); + } + + public double MatrixNorm(Norm norm, double[] matrix, double[] work) + { + throw new NotImplementedException(); + } + + public void MatrixMultiply(double[] x, double[] y, double[] result) + { + throw new NotImplementedException(); + } + + public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, double[] b, double beta, double[] c) + { + throw new NotImplementedException(); + } + + public void LUFactor(double[] a, int[] ipiv) + { + throw new NotImplementedException(); + } + + public void LUInverse(double[] a) + { + throw new NotImplementedException(); + } + + public void LUInverseFactored(double[] a, int[] ipiv) + { + throw new NotImplementedException(); + } + + public void LUInverse(double[] a, double[] work) + { + throw new NotImplementedException(); + } + + public void LUInverseFactored(double[] a, int[] ipiv, double[] work) + { + throw new NotImplementedException(); + } + + public void LUSolve(int columnsOfB, double[] a, double[] b) + { + throw new NotImplementedException(); + } + + public void LUSolveFactored(int columnsOfB, double[] a, int ipiv, double[] b) + { + throw new NotImplementedException(); + } + + public void LUSolve(Transpose transposeA, int columnsOfB, double[] a, double[] b) + { + throw new NotImplementedException(); + } + + public void LUSolveFactored(Transpose transposeA, int columnsOfB, double[] a, int ipiv, double[] b) + { + throw new NotImplementedException(); + } + + public void CholeskyFactor(double[] a) + { + throw new NotImplementedException(); + } + + public void CholeskySolve(int columnsOfB, double[] a, double[] b) + { + throw new NotImplementedException(); + } + + public void CholeskySolveFactored(int columnsOfB, double[] a, double[] b) + { + throw new NotImplementedException(); + } + + public void QRFactor(double[] a, double[] q, double[] tau) + { + throw new NotImplementedException(); + } + + public void QRFactor(double[] a, double[] q, double[] tau, double[] work) + { + throw new NotImplementedException(); + } + + public void QRSolve(int columnsOfB, double[] a, double[] tau, double[] b, double[] x) + { + throw new NotImplementedException(); + } + + public void QRSolve(int columnsOfB, double[] a, double[] tau, double[] b, double[] x, double[] work) + { + throw new NotImplementedException(); + } + + public void SinguarValueDecomposition(bool computeVectors, double[] a, double[] s, double[] u, double[] vt) + { + throw new NotImplementedException(); + } + + public void SinguarValueDecomposition(bool computeVectors, double[] a, double[] s, double[] u, double[] vt, double[] work) + { + throw new NotImplementedException(); + } + + public void SvdSolve(double[] s, double[] u, double[] vt, double[] b, double[] x) + { + throw new NotImplementedException(); + } + + #endregion } } \ No newline at end of file