diff --git a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs new file mode 100644 index 00000000..b10d4e0c --- /dev/null +++ b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs @@ -0,0 +1,440 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// Copyright (c) 2009 Math.NET +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +// 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. + /// + /// Supported data types are double, single, , and . + public interface ILinearAlgebraProvider where T : struct + { + /// + /// 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. + /// + /// The vector to update. + /// The value to scale by. + /// The vector to add to . + /// This is equivalent to the AXPY BLAS routine. + void AddVectorToScaledVector(T[] y, T alpha, T[] x); + + /// + /// Scales an array. Can be used to scale a vector and a matrix. + /// + /// The scalar. + /// The values to scale. + /// This is equivalent to the SCAL BLAS routine. + void ScaleArray(T alpha, T[] 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. + T DotProduct(T[] x, T[] 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(T[] x, T[] y, T[] 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(T[] x, T[] y, T[] 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(T[] x, T[] y, T[] 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. + /// + T MatrixNorm(Norm norm, T[] 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. + /// + T MatrixNorm(Norm norm, T[] matrix, T[] 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(T[] x, T[] y, T[] 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, T alpha, T[] a, T[] b, T beta, T[] 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(T[] 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(T[] 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(T[] 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(T[] a, T[] 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(T[] a, int[] ipiv, T[] 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, T[] a, T[] 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, T[] a, int ipiv, T[] 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, T[] a, T[] 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, T[] a, int ipiv, T[] b); + + /// + /// Computes the Cholesky factorization of A. + /// + /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the + /// the Cholesky factorization. + /// This is equivalent to the POTRF LAPACK routine. + void CholeskyFactor(T[] 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, T[] a, T[] 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, T[] a, T[] b); + + /// + /// Computes the QR factorization of A. + /// + /// On entry, it is the M by N A matrix to factor. On exit, + /// it is overwritten with the R matrix of the QR factorization. + /// On exit, A M by M matrix that holds the Q matrix of the + /// QR factorization. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + void QRFactor(T[] r, T[] q); + + /// + /// 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. + /// On exit, A M by M matrix that holds the Q matrix of the + /// QR factorization. + /// 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 similar to the GEQRF and ORGQR LAPACK routines. + void QRFactor(T[] r, T[] q, T[] work); + + /// + /// Solves A*X=B for X using QR factorization of A. + /// + /// The number of columns of B. + /// 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. + /// On exit, A M by M matrix that holds the Q matrix of the + /// QR factorization. + /// The B matrix. + /// On exit, the solution matrix. + void QRSolve(int columnsOfB, T[] r, T[] q, T[] b, T[] x); + + /// + /// Solves A*X=B for X using QR factorization of A. + /// + /// The number of columns of B. + /// 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. + /// On exit, A M by M matrix that holds the Q matrix of the + /// QR factorization. + /// 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. + void QRSolve(int columnsOfB, T[] r, T[] q, T[] b, T[] x, T[] work); + + /// + /// Solves A*X=B for X using a previously QR factored matrix. + /// + /// The number of columns of B. + /// The Q matrix obtained by calling . + /// The R matrix obtained by calling . + /// The B matrix. + /// On exit, the solution matrix. + void QRSolveFactored(int columnsOfB, T[] q, T[] r, T[] b, T[] x); + + /// + /// 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, T[] a, 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 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. For real matrices, the work array should be at least + /// Max(3*Min(M, N) + Max(M, N), 5*Min(M,N)). For complex matrices, 2*Min(M, N) + Max(M, N). + /// On exit, work[0] contains the optimal work size value. + /// + /// This is equivalent to the GESVD LAPACK routine. + void SingularValueDecomposition(bool computeVectors, T[] a, T[] s, T[] u, T[] vt, T[] work); + + /// + /// Solves A*X=B for X using the singular value decomposition of A. + /// + /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. + /// The singular values of A in ascending value. + /// On exit U contains the left singular vectors. + /// On exit VT contains the transposed right singular vectors. + /// The B matrix. + /// On exit, the solution matrix. + void SvdSolve(T[] a, T[] s, T[] u, T[] vt, T[] b, T[] x); + + /// + /// Solves A*X=B for X using the singular value decomposition of A. + /// + /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. + /// The singular values of A in ascending value. + /// On exit U contains the left singular vectors. + /// On exit VT contains the transposed right singular vectors. + /// The B matrix. + /// On exit, the solution matrix. + /// The work array. For real matrices, the work array should be at least + /// Max(3*Min(M, N) + Max(M, N), 5*Min(M,N)). For complex matrices, 2*Min(M, N) + Max(M, N). + /// On exit, work[0] contains the optimal work size value. + /// + void SvdSolve(T[] a, T[] s, T[] u, T[] vt, T[] b, T[] x, T[] work); + + /// + /// Solves A*X=B for X using a previously SVD decomposed matrix. + /// + /// The number of columns of B. + /// 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 SvdSolveFactored(int columnsOfB, T[] s, T[] u, T[] vt, T[] b, T[] x); + } +}