flushed out linear algebra provider interface

src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProvider.cs

@@ -1,4 +1,4 @@
-// <copyright file="ILinearAlgebra.cs" company="Math.NET">
+// <copyright file="ILinearAlgebraProvider.cs" company="Math.NET">
 // 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.
 // </copyright>
 
+// INITIAL DRAFT MISSING EXCEPTION SPECIFICATIONS
 namespace MathNet.Numerics.Algorithms.LinearAlgebra
 {
+    /// <summary>
+    /// How to transpose a matrix.
+    /// </summary>
+    public enum Transpose
+    {
+        /// <summary>
+        /// Don't transpose a matrix.
+        /// </summary>
+        DontTranspose = 111,
+
+        /// <summary>
+        /// Transpose a matrix.
+        /// </summary>
+        Transpose = 112,
+
+        /// <summary>
+        /// Conjugate transpose a complex matrix.
+        /// </summary>
+        /// <remarks>If a conjugate transpose is used with a real matrix, then the matrix is just transposed.</remarks>
+        ConjugateTranspose = 113
+    }
+
+    /// <summary>
+    /// Types of matrix norms.
+    /// </summary>
+    public enum Norm : byte
+    {
+        /// <summary>
+        /// The 1-norm.
+        /// </summary>
+        OneNorm = (byte)'1',
+
+        /// <summary>
+        /// The Frobenius norm.
+        /// </summary>
+        FrobeniusNorm = (byte)'f',
+
+        /// <summary>
+        /// The infinity norm.
+        /// </summary>
+        InfinityNorm = (byte)'i',
+
+        /// <summary>
+        /// The largest absolute value norm.
+        /// </summary>
+        LargestAbsoluteValue = (byte)'m'
+    }
+
     /// <summary>
     /// Interface to linear algebra algorithms that work off 1-D arrays.
     /// </summary>
     public interface ILinearAlgebraProvider
     {
+        /// <summary>
+        /// Queries the provider for the optimal, workspace block size
+        /// for the given routine.
+        /// </summary>
+        /// <param name="methodName">Name of the method to query.</param>
+        /// <returns>-1 if the provider cannot compute the workspace size; otherwise
+        /// the suggested block size.</returns>
+        int QueryWorkspaceBlockSize(string methodName);
+        
         /// <summary>
         /// Adds a scaled vector to another: <c>y += alpha*x</c>.
         /// </summary>
@@ -45,5 +103,310 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="x">The values to scale.</param>
         /// <remarks>This is equivalent to the SCAL BLAS routine.</remarks>
         void ScaleArray(double alpha, double[] x);
+
+        /// <summary>
+        /// Computes the dot product of x and y.
+        /// </summary>
+        /// <param name="x">The vector x.</param>
+        /// <param name="y">The vector y.</param>
+        /// <returns>The dot product of x and y.</returns>
+        /// <remarks>This is equivalent to the DOT BLAS routine.</remarks>
+        double DotProduct(double[] x, double[] y);
+
+        /// <summary>
+        /// Does a point wise add of two arrays <c>z = x + y</c>. This can be used 
+        /// to add vectors or matrices.
+        /// </summary>
+        /// <param name="x">The array x.</param>
+        /// <param name="y">The array y.</param>
+        /// <param name="result">The result of the addition.</param>
+        /// <remarks>There is no equivalent BLAS routine, but many libraries
+        /// provide optimized (parallel and/or vectorized) versions of this
+        /// routine.</remarks>
+        void AddArrays(double[] x, double[] y, double[] result);
+
+        /// <summary>
+        /// Does a point wise subtraction of two arrays <c>z = x - y</c>. This can be used 
+        /// to subtract vectors or matrices.
+        /// </summary>
+        /// <param name="x">The array x.</param>
+        /// <param name="y">The array y.</param>
+        /// <param name="result">The result of the subtraction.</param>
+        /// <remarks>There is no equivalent BLAS routine, but many libraries
+        /// provide optimized (parallel and/or vectorized) versions of this
+        /// routine.</remarks>
+        void SubtractArrays(double[] x, double[] y, double[] result);
+
+        /// <summary>
+        /// Does a point wise multiplication of two arrays <c>z = x * y</c>. This can be used
+        /// to multiple elements of vectors or matrices.
+        /// </summary>
+        /// <param name="x">The array x.</param>
+        /// <param name="y">The array y.</param>
+        /// <param name="result">The result of the point wise multiplication.</param>
+        /// <remarks>There is no equivalent BLAS routine, but many libraries
+        /// provide optimized (parallel and/or vectorized) versions of this
+        /// routine.</remarks>
+        void PointWiseMultiplyArrays(double[] x, double[] y, double[] result);
+
+        /// <summary>
+        /// Computes the requested <see cref="Norm"/> of the matrix.
+        /// </summary>
+        /// <param name="norm">The type of norm to compute.</param>
+        /// <param name="matrix">The matrix to compute the norm from.</param>
+        /// <returns>
+        /// The requested <see cref="Norm"/> of the matrix.
+        /// </returns>
+        double MatrixNorm(Norm norm, double[] matrix);
+
+        /// <summary>
+        /// Computes the requested <see cref="Norm"/> of the matrix.
+        /// </summary>
+        /// <param name="norm">The type of norm to compute.</param>
+        /// <param name="matrix">The matrix to compute the norm from.</param>
+        /// <param name="work">The work array. Only used when <see cref="Norm.InfinityNorm"/>
+        /// and needs to be have a length of at least M (number of rows of <paramref name="matrix"/>.</param>
+        /// <returns>
+        /// The requested <see cref="Norm"/> of the matrix.
+        /// </returns>
+        double MatrixNorm(Norm norm, double[] matrix, double[] work);
+
+        /// <summary>
+        /// Multiples two matrices. <c>result = x * y</c>
+        /// </summary>
+        /// <param name="x">The x matrix.</param>
+        /// <param name="y">The y matrix.</param>
+        /// <param name="result">Where to store the result of the multiplication.</param>
+        /// <remarks>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.</remarks>
+        void MatrixMultiply(double[] x, double[] y, double[] result);
+
+        /// <summary>
+        /// Multiplies two matrices and updates another with the result. <c>c = alpha*op(a)*op(b) + beta*c</c>
+        /// </summary>
+        /// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
+        /// <param name="transposeB">How to transpose the <paramref name="b"/> matrix.</param>
+        /// <param name="alpha">The value to scale <paramref name="a"/> matrix.</param>
+        /// <param name="a">The a matrix.</param>
+        /// <param name="b">The b matrix</param>
+        /// <param name="beta">The value to scale the <paramref name="c"/> matrix.</param>
+        /// <param name="c">The c matrix.</param>
+        void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, double[] b, double beta, double[] c);
+
+        /// <summary>
+        /// Computes the LU factorization of A.
+        /// </summary>
+        /// <param name="a">An m by n matrix. The matrix is overwritten with the
+        /// the LU factorization On exit.</param>
+        /// <param name="ipiv">On exit, it contains the pivot indices. The size
+        /// of the array must be min(m,n).</param>
+        /// <remarks>This is equivalent to the GETRF LAPACK routine.</remarks>
+        void LUFactor(double[] a, int[] ipiv);
+
+        /// <summary>
+        /// Computes the inverse of matrix using LU factorization.
+        /// </summary>
+        /// <param name="a">The N by N matrix to invert. Contains the inverse On exit.</param>
+        /// <remarks>This is equivalent to the GETRF and GETRI LAPACK routines.</remarks>
+        void LUInverse(double[] a);
+
+        /// <summary>
+        /// Computes the inverse of a previously factored matrix.
+        /// </summary>
+        /// <param name="a">The LU factored N by N matrix.  Contains the inverse On exit.</param>
+        /// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
+        /// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
+        void LUInverseFactored(double[] a, int[] ipiv);
+
+        /// <summary>
+        /// Computes the inverse of matrix using LU factorization.
+        /// </summary>
+        /// <param name="a">The N by N matrix to invert. Contains the inverse On exit.</param>
+        /// <param name="work">The work array. The array must have a length of at least N,
+        /// but should be N*blocksize. The blocksize is machine dependent. Use <see cref="QueryWorkspaceBlockSize"/>
+        /// to determine the optimal size of the work array. On exit, work[0] contains the optimal
+        /// work size value.</param>
+        /// <remarks>This is equivalent to the GETRF and GETRI LAPACK routines.</remarks>
+        void LUInverse(double[] a, double[] work);
+
+        /// <summary>
+        /// Computes the inverse of a previously factored matrix.
+        /// </summary>
+        /// <param name="a">The LU factored N by N matrix.  Contains the inverse On exit.</param>
+        /// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
+        /// <param name="work">The work array. The array must have a length of at least N,
+        /// but should be N*blocksize. The blocksize is machine dependent. Use <see cref="QueryWorkspaceBlockSize"/>
+        /// to determine the optimal size of the work array. On exit, work[0] contains the optimal
+        /// work size value.</param>
+        /// <remarks>This is equivalent to the GETRI LAPACK routine.</remarks>
+        void LUInverseFactored(double[] a, int[] ipiv, double[] work);
+
+        /// <summary>
+        /// Solves A*X=B for X using LU factorization.
+        /// </summary>
+        /// <param name="columnsOfB">The number of columns of B.</param>
+        /// <param name="a">The square matrix A.</param>
+        /// <param name="b">The B matrix.</param>
+        /// <remarks>This is equivalent to the GETRF and GETRS LAPACK routines.</remarks>
+        void LUSolve(int columnsOfB, double[] a, double[] b);
+
+        /// <summary>
+        /// Solves A*X=B for X using a previously factored A matrix.
+        /// </summary>
+        /// <param name="columnsOfB">The number of columns of B.</param>
+        /// <param name="a">The factored A matrix.</param>
+        /// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
+        /// <param name="b">The B matrix.</param>
+        /// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
+        void LUSolveFactored(int columnsOfB, double[] a, int ipiv, double[] b);
+
+        /// <summary>
+        /// Solves A*X=B for X using LU factorization.
+        /// </summary>
+        /// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
+        /// <param name="columnsOfB">The number of columns of B.</param>
+        /// <param name="a">The square matrix A.</param>
+        /// <param name="b">The B matrix.</param>
+        /// <remarks>This is equivalent to the GETRF and GETRS LAPACK routines.</remarks>
+        void LUSolve(Transpose transposeA, int columnsOfB, double[] a, double[] b);
+
+        /// <summary>
+        /// Solves A*X=B for X using a previously factored A matrix.
+        /// </summary>
+        /// <param name="transposeA">How to transpose the <paramref name="a"/> matrix.</param>
+        /// <param name="columnsOfB">The number of columns of B.</param>
+        /// <param name="a">The factored A matrix.</param>
+        /// <param name="ipiv">The pivot indices of <paramref name="a"/>.</param>
+        /// <param name="b">The B matrix.</param>
+        /// <remarks>This is equivalent to the GETRS LAPACK routine.</remarks>
+        void LUSolveFactored(Transpose transposeA, int columnsOfB, double[] a, int ipiv, double[] b);
+
+        /// <summary>
+        /// Computes the Cholesky factorization of A.
+        /// </summary>
+        /// <param name="a">A square, positive definite matrix. The matrix is overwritten with the
+        /// the Cholesky factorization On exit.</param>
+        /// <remarks>This is equivalent to the POTRF LAPACK routine.</remarks>
+        void CholeskyFactor(double[] a);
+
+        /// <summary>
+        /// Solves A*X=B for X using Cholesky factorization.
+        /// </summary>
+        /// <param name="columnsOfB">The number of columns of B.</param>
+        /// <param name="a">The square, positive definite matrix A.</param>
+        /// <param name="b">The B matrix.</param>
+        /// <remarks>This is equivalent to the POTRF add POTRS LAPACK routines.</remarks>
+        void CholeskySolve(int columnsOfB, double[] a, double[] b);
+
+        /// <summary>
+        /// Solves A*X=B for X using a previously factored A matrix.
+        /// </summary>
+        /// <param name="columnsOfB">The number of columns of B.</param>
+        /// <param name="a">The factored A matrix.</param>
+        /// <param name="b">The B matrix.</param>
+        /// <remarks>This is equivalent to the POTRS LAPACK routine.</remarks>
+        void CholeskySolveFactored(int columnsOfB, double[] a, double[] b);
+
+        /// <summary>
+        /// Computes the QR factorization of A.
+        /// </summary>
+        /// <param name="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 <paramref name="tau"/>, represent the orthogonal matrix Q as a
+        /// product of min(m,n) elementary reflectors. </param>
+        /// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
+        /// QR factorization.</param>
+        /// <param name="tau">On exit, tau contains information needed by the <see cref="QRSolve(int,double[],double[],double[],double[])"/>
+        /// method.</param>
+        /// <remarks>This is equivalent to the GEQRF and ORGQR LAPACK routines.</remarks>
+        void QRFactor(double[] a, double[] q, double[] tau);
+
+        /// <summary>
+        /// Computes the QR factorization of A.
+        /// </summary>
+        /// <param name="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 <paramref name="tau"/>, represent the orthogonal matrix Q as a
+        /// product of min(m,n) elementary reflectors. </param>
+        /// <param name="q">On exit, A M by M matrix that holds the Q matrix of the 
+        /// QR factorization.</param>
+        /// <param name="tau">On exit, tau contains information needed by the <see cref="QRSolve(int,double[],double[],double[],double[])"/>
+        /// method.</param>
+        /// <param name="work">The work array. The array must have a length of at least N,
+        /// but should be N*blocksize. The blocksize is machine dependent. Use <see cref="QueryWorkspaceBlockSize"/>
+        /// to determine the optimal size of the work array. On exit, work[0] contains the optimal
+        /// work size value.</param>
+        /// <remarks>This is equivalent to the GEQRF and ORGQR LAPACK routines.</remarks>
+        void QRFactor(double[] a, double[] q, double[] tau, double[] work);
+
+        /// <summary>
+        /// Solves A*X=B for X using a previously QR factored matrix.
+        /// </summary>
+        /// <param name="columnsOfB">The number of columns of B.</param>
+        /// <param name="a">The matrix obtained by calling <see cref="QRFactor(double[],double[],double[])"/>.</param>
+        /// <param name="tau">The tau vector obtained by calling <see cref="QRFactor(double[],double[],double[])"/>. </param>
+        /// <param name="b">The B matrix.</param>
+        /// <param name="x">On exit, the solution matrix.</param>
+        /// <remarks>This is equivalent to the ORMQR LAPACK routine with the TRSM BLAS routine.</remarks>
+        void QRSolve(int columnsOfB, double[] a, double[] tau, double[] b, double[] x);
+
+        /// <summary>
+        /// Solves A*X=B for X using a previously QR factored matrix.
+        /// </summary>
+        /// <param name="columnsOfB">The number of columns of B.</param>
+        /// <param name="a">The M by N</param>
+        /// <param name="tau">The tau vector obtained by calling <see cref="QRFactor(double[],double[],double[])"/>. </param>
+        /// <param name="b">The B matrix.</param>
+        /// <param name="x">On exit, the solution matrix.</param>
+        /// <param name="work">The work array. The array must have a length of at least N,
+        /// but should be N*blocksize. The blocksize is machine dependent. Use <see cref="QueryWorkspaceBlockSize"/>
+        /// to determine the optimal size of the work array. On exit, work[0] contains the optimal
+        /// work size value.</param>
+        /// <remarks>This is equivalent to the ORMQR LAPACK routine with the TRSM BLAS routine.</remarks>
+        void QRSolve(int columnsOfB, double[] a, double[] tau, double[] b, double[] x, double[] work);
+
+        /// <summary>
+        /// Computes the singular value decomposition of A.
+        /// </summary>
+        /// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
+        /// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
+        /// <param name="s">The singular values of A in ascending value. </param>
+        /// <param name="u">If <paramref name="computeVectors"/> is true, on exit U contains the left
+        /// singular vectors.</param>
+        /// <param name="vt">If <paramref name="computeVectors"/> is true, on exit VT contains the transposed
+        /// right singular vectors.</param>
+        /// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
+        void SinguarValueDecomposition(bool computeVectors, double[] a, double[] s, double[] u, double[] vt);
+
+        /// <summary>
+        /// Computes the singular value decomposition of A.
+        /// </summary>
+        /// <param name="computeVectors">Compute the singular U and VT vectors or not.</param>
+        /// <param name="a">On entry, the M by N matrix to decompose. On exit, A may be overwritten.</param>
+        /// <param name="s">The singular values of A in ascending value. </param>
+        /// <param name="u">If <paramref name="computeVectors"/> is true, on exit U contains the left
+        /// singular vectors.</param>
+        /// <param name="vt">If <paramref name="computeVectors"/> is true, on exit VT contains the transposed
+        /// right singular vectors.</param>
+        /// <param name="work">The work array. The array must have a length of at least N,
+        /// but should be N*blocksize. The blocksize is machine dependent. Use <see cref="QueryWorkspaceBlockSize"/>
+        /// to determine the optimal size of the work array. On exit, work[0] contains the optimal
+        /// work size value.</param>
+        /// <remarks>This is equivalent to the GESVD LAPACK routine.</remarks>
+        void SingularValueDecomposition(bool computeVectors, double[] a, double[] s, double[] u, double[] vt, double[] work);
+
+        /// <summary>
+        /// Solves A*X=B for X using a previously SVD decomposed matrix.
+        /// </summary>
+        /// <param name="s">The s values returned by <see cref="SinguarValueDecomposition(bool,double[],double[],double[],double[])"/>.</param>
+        /// <param name="u">The left singular vectors returned by  <see cref="SinguarValueDecomposition(bool,double[],double[],double[],double[])"/>.</param>
+        /// <param name="vt">The right singular  vectors returned by  <see cref="SinguarValueDecomposition(bool,double[],double[],double[],double[])"/>.</param>
+        /// <param name="b">The B matrix.</param>
+        /// <param name="x">On exit, the solution matrix.</param>
+        void SvdSolve(double[] s, double[] u, double[] vt, double[] b, double[] x);
     }
 }

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
     }
 }