diff --git a/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs index c28750a8..41c52f77 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/Atlas/AtlasLinearAlgebraProvider.cs @@ -28,7 +28,7 @@ /* This file is automatically generated - do not modify it. Change NativeLinearAlgebraProvider.include instead. - Last generated on UTC 2010-06-27 13:41:23Z + Last generated on UTC 2010-07-04 13:21:14Z */ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas @@ -235,6 +235,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// Computes the requested of the matrix. /// /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. /// The matrix to compute the norm from. /// /// The requested of the matrix. @@ -248,6 +250,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// Computes the requested of the matrix. /// /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. /// 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 . @@ -452,24 +456,25 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas { throw new ArgumentNullException("a"); } - - if ( order < 1) + + if (order < 1) { - throw new ArgumentException(Properties.Resources.ArgumentMustBePositive, "order"); + throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } - - SafeNativeMethods.d_cholesky_factor(order, a); + + SafeNativeMethods.d_cholesky_factor(order, a); } - + /// /// Solves A*X=B for X using Cholesky factorization. /// /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. - /// The number of columns in the B matrix. - /// This is equivalent to the POTRF add POTRS LAPACK routines. + /// The number of rows in the B matrix. + /// The number of columns in the B matrix. + /// This is equivalent to the POTRF add POTRS LAPACK routines. + /// public void CholeskySolve(double[] a, int aOrder, double[] b, int bRows, int bColumns) { throw new NotImplementedException(); @@ -481,8 +486,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. - /// The number of columns in the B matrix. + /// The number of rows in the B matrix. + /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. public void CholeskySolveFactored(double[] a, int aOrder, double[] b, int bRows, int bColumns) { @@ -493,11 +498,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// 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. + /// 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. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public void QRFactor(double[] r, double[] q) + public void QRFactor(double[] r, int rRows, int rColumns, double[] q) { throw new NotImplementedException(); } @@ -506,13 +513,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// 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 + /// 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. /// 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. - public void QRFactor(double[] r, double[] q, double[] work) + /// This is similar to the GEQRF and ORGQR LAPACK routines. + public void QRFactor(double[] r, int rRows, int rColumns, double[] q, double[] work) { throw new NotImplementedException(); } @@ -520,14 +530,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// 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 + /// 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. /// The B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void QRSolve(int columnsOfB, double[] r, double[] q, double[] b, double[] x) + public void QRSolve(double[] r, int rRows, int rColumns, double[] q, double[] b, int bColumns, double[] x) { throw new NotImplementedException(); } @@ -535,17 +547,19 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// 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 + /// 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. /// 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. - public void QRSolve(int columnsOfB, double[] r, double[] q, double[] b, double[] x, double[] work) + public void QRSolve(double[] r, int rRows, int rColumns, double[] q, double[] b, int bColumns, double[] x, double[] work) { throw new NotImplementedException(); } @@ -553,12 +567,14 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// 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 Q matrix obtained by calling . + /// The R matrix obtained by calling . + /// 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. - public void QRSolveFactored(int columnsOfB, double[] q, double[] r, double[] b, double[] x) + public void QRSolveFactored(double[] q, double[] r, int rRows, int rColumns, double[] b, int bColumns, double[] x) { throw new NotImplementedException(); } @@ -568,13 +584,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// 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. - public void SingularValueDecomposition(bool computeVectors, double[] a, double[] s, double[] u, double[] vt) + public void SingularValueDecomposition(bool computeVectors, double[] a, int aRows, int aColumns, double[] s, double[] u, double[] vt) { throw new NotImplementedException(); } @@ -584,6 +602,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// 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. @@ -593,7 +613,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// 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. - public void SingularValueDecomposition(bool computeVectors, double[] a, double[] s, double[] u, double[] vt, double[] work) + public void SingularValueDecomposition(bool computeVectors, double[] a, int aRows, int aColumns, double[] s, double[] u, double[] vt, double[] work) { throw new NotImplementedException(); } @@ -602,12 +622,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// On exit, the solution matrix. - public void SvdSolve(double[] a, double[] s, double[] u, double[] vt, double[] b, double[] x) + public void SvdSolve(double[] a, int aRows, int aColumns, double[] s, double[] u, double[] vt, double[] b, int bColumns, double[] x) { throw new NotImplementedException(); } @@ -616,15 +639,18 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// 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. - public void SvdSolve(double[] a, double[] s, double[] u, double[] vt, double[] b, double[] x, double[] work) + public void SvdSolve(double[] a, int aRows, int aColumns, double[] s, double[] u, double[] vt, double[] b, int bColumns, double[] x, double[] work) { throw new NotImplementedException(); } @@ -632,13 +658,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// 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 number of rows in the A matrix. + /// 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 B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void SvdSolveFactored(int columnsOfB, double[] s, double[] u, double[] vt, double[] b, double[] x) + public void SvdSolveFactored(int aRows, int aColumns, double[] s, double[] u, double[] vt, double[] b, int bColumns, double[] x) { throw new NotImplementedException(); } @@ -836,6 +864,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// Computes the requested of the matrix. /// /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. /// The matrix to compute the norm from. /// /// The requested of the matrix. @@ -849,6 +879,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// Computes the requested of the matrix. /// /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. /// 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 . @@ -926,6 +958,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas SafeNativeMethods.s_matrix_multiply(transposeA, transposeB, m, n, k, alpha, a, b, beta, c); } + /// /// /// Computes the LUP factorization of A. P*A = L*U. /// @@ -1053,23 +1086,23 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas { throw new ArgumentNullException("a"); } - - if ( order < 1) + + if (order < 1) { - throw new ArgumentException(Properties.Resources.ArgumentMustBePositive, "order"); + throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } - - SafeNativeMethods.s_cholesky_factor(order, a); + + SafeNativeMethods.s_cholesky_factor(order, a); } - + /// /// Solves A*X=B for X using Cholesky factorization. /// /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. - /// The number of columns in the B matrix. + /// The number of rows in the B matrix. + /// The number of columns in the B matrix. /// This is equivalent to the POTRF add POTRS LAPACK routines. public void CholeskySolve(float[] a, int aOrder, float[] b, int bRows, int bColumns) { @@ -1082,8 +1115,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. - /// The number of columns in the B matrix. + /// The number of rows in the B matrix. + /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. public void CholeskySolveFactored(float[] a, int aOrder, float[] b, int bRows, int bColumns) { @@ -1094,11 +1127,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// 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. + /// 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. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public void QRFactor(float[] r, float[] q) + public void QRFactor(float[] r, int rRows, int rColumns, float[] q) { throw new NotImplementedException(); } @@ -1107,13 +1142,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// 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 + /// 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. /// 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. - public void QRFactor(float[] r, float[] q, float[] work) + /// This is similar to the GEQRF and ORGQR LAPACK routines. + public void QRFactor(float[] r, int rRows, int rColumns, float[] q, float[] work) { throw new NotImplementedException(); } @@ -1121,14 +1159,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// 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 + /// 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. /// The B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void QRSolve(int columnsOfB, float[] r, float[] q, float[] b, float[] x) + public void QRSolve(float[] r, int rRows, int rColumns, float[] q, float[] b, int bColumns, float[] x) { throw new NotImplementedException(); } @@ -1136,17 +1176,19 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// 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 + /// 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. /// 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. - public void QRSolve(int columnsOfB, float[] r, float[] q, float[] b, float[] x, float[] work) + public void QRSolve(float[] r, int rRows, int rColumns, float[] q, float[] b, int bColumns, float[] x, float[] work) { throw new NotImplementedException(); } @@ -1154,12 +1196,14 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// 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 Q matrix obtained by calling . + /// The R matrix obtained by calling . + /// 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. - public void QRSolveFactored(int columnsOfB, float[] q, float[] r, float[] b, float[] x) + public void QRSolveFactored(float[] q, float[] r, int rRows, int rColumns, float[] b, int bColumns, float[] x) { throw new NotImplementedException(); } @@ -1169,13 +1213,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// 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. - public void SingularValueDecomposition(bool computeVectors, float[] a, float[] s, float[] u, float[] vt) + public void SingularValueDecomposition(bool computeVectors, float[] a, int aRows, int aColumns, float[] s, float[] u, float[] vt) { throw new NotImplementedException(); } @@ -1185,6 +1231,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// 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. @@ -1194,7 +1242,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// 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. - public void SingularValueDecomposition(bool computeVectors, float[] a, float[] s, float[] u, float[] vt, float[] work) + public void SingularValueDecomposition(bool computeVectors, float[] a, int aRows, int aColumns, float[] s, float[] u, float[] vt, float[] work) { throw new NotImplementedException(); } @@ -1203,12 +1251,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// On exit, the solution matrix. - public void SvdSolve(float[] a, float[] s, float[] u, float[] vt, float[] b, float[] x) + public void SvdSolve(float[] a, int aRows, int aColumns, float[] s, float[] u, float[] vt, float[] b, int bColumns, float[] x) { throw new NotImplementedException(); } @@ -1217,15 +1268,18 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// 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. - public void SvdSolve(float[] a, float[] s, float[] u, float[] vt, float[] b, float[] x, float[] work) + public void SvdSolve(float[] a, int aRows, int aColumns, float[] s, float[] u, float[] vt, float[] b, int bColumns, float[] x, float[] work) { throw new NotImplementedException(); } @@ -1233,13 +1287,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// 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 number of rows in the A matrix. + /// 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 B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void SvdSolveFactored(int columnsOfB, float[] s, float[] u, float[] vt, float[] b, float[] x) + public void SvdSolveFactored(int aRows, int aColumns, float[] s, float[] u, float[] vt, float[] b, int bColumns, float[] x) { throw new NotImplementedException(); } @@ -1437,6 +1493,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// Computes the requested of the matrix. /// /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. /// The matrix to compute the norm from. /// /// The requested of the matrix. @@ -1450,6 +1508,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// Computes the requested of the matrix. /// /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. /// 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 . @@ -1654,23 +1714,23 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas { throw new ArgumentNullException("a"); } - - if ( order < 1) + + if (order < 1) { - throw new ArgumentException(Properties.Resources.ArgumentMustBePositive, "order"); + throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } - - SafeNativeMethods.z_cholesky_factor(order, a); + + SafeNativeMethods.z_cholesky_factor(order, a); } - + /// /// Solves A*X=B for X using Cholesky factorization. /// /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. - /// The number of columns in the B matrix. + /// The number of rows in the B matrix. + /// The number of columns in the B matrix. /// This is equivalent to the POTRF add POTRS LAPACK routines. public void CholeskySolve(Complex[] a, int aOrder, Complex[] b, int bRows, int bColumns) { @@ -1683,8 +1743,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. - /// The number of columns in the B matrix. + /// The number of rows in the B matrix. + /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. public void CholeskySolveFactored(Complex[] a, int aOrder, Complex[] b, int bRows, int bColumns) { @@ -1695,11 +1755,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// 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. + /// 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. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public void QRFactor(Complex[] r, Complex[] q) + public void QRFactor(Complex[] r, int rRows, int rColumns, Complex[] q) { throw new NotImplementedException(); } @@ -1708,13 +1770,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// 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 + /// 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. /// 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. - public void QRFactor(Complex[] r, Complex[] q, Complex[] work) + /// This is similar to the GEQRF and ORGQR LAPACK routines. + public void QRFactor(Complex[] r, int rRows, int rColumns, Complex[] q, Complex[] work) { throw new NotImplementedException(); } @@ -1722,14 +1787,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// 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 + /// 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. /// The B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void QRSolve(int columnsOfB, Complex[] r, Complex[] q, Complex[] b, Complex[] x) + public void QRSolve(Complex[] r, int rRows, int rColumns, Complex[] q, Complex[] b, int bColumns, Complex[] x) { throw new NotImplementedException(); } @@ -1737,17 +1804,19 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// 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 + /// 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. /// 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. - public void QRSolve(int columnsOfB, Complex[] r, Complex[] q, Complex[] b, Complex[] x, Complex[] work) + public void QRSolve(Complex[] r, int rRows, int rColumns, Complex[] q, Complex[] b, int bColumns, Complex[] x, Complex[] work) { throw new NotImplementedException(); } @@ -1755,12 +1824,14 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// 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 Q matrix obtained by calling . + /// The R matrix obtained by calling . + /// 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. - public void QRSolveFactored(int columnsOfB, Complex[] q, Complex[] r, Complex[] b, Complex[] x) + public void QRSolveFactored(Complex[] q, Complex[] r, int rRows, int rColumns, Complex[] b, int bColumns, Complex[] x) { throw new NotImplementedException(); } @@ -1770,13 +1841,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// 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. - public void SingularValueDecomposition(bool computeVectors, Complex[] a, Complex[] s, Complex[] u, Complex[] vt) + public void SingularValueDecomposition(bool computeVectors, Complex[] a, int aRows, int aColumns, Complex[] s, Complex[] u, Complex[] vt) { throw new NotImplementedException(); } @@ -1786,6 +1859,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// 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. @@ -1795,7 +1870,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// 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. - public void SingularValueDecomposition(bool computeVectors, Complex[] a, Complex[] s, Complex[] u, Complex[] vt, Complex[] work) + public void SingularValueDecomposition(bool computeVectors, Complex[] a, int aRows, int aColumns, Complex[] s, Complex[] u, Complex[] vt, Complex[] work) { throw new NotImplementedException(); } @@ -1804,12 +1879,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// On exit, the solution matrix. - public void SvdSolve(Complex[] a, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, Complex[] x) + public void SvdSolve(Complex[] a, int aRows, int aColumns, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, int bColumns, Complex[] x) { throw new NotImplementedException(); } @@ -1818,15 +1896,18 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// 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. - public void SvdSolve(Complex[] a, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, Complex[] x, Complex[] work) + public void SvdSolve(Complex[] a, int aRows, int aColumns, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, int bColumns, Complex[] x, Complex[] work) { throw new NotImplementedException(); } @@ -1834,13 +1915,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// 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 number of rows in the A matrix. + /// 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 B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void SvdSolveFactored(int columnsOfB, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, Complex[] x) + public void SvdSolveFactored(int aRows, int aColumns, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, int bColumns, Complex[] x) { throw new NotImplementedException(); } @@ -2259,23 +2342,23 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas { throw new ArgumentNullException("a"); } - - if ( order < 1) + + if (order < 1) { - throw new ArgumentException(Properties.Resources.ArgumentMustBePositive, "order"); + throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } - - SafeNativeMethods.c_cholesky_factor(order, a); + + SafeNativeMethods.c_cholesky_factor(order, a); } - + /// /// Solves A*X=B for X using Cholesky factorization. /// /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. - /// The number of columns in the B matrix. + /// The number of rows in the B matrix. + /// The number of columns in the B matrix. /// This is equivalent to the POTRF add POTRS LAPACK routines. public void CholeskySolve(Complex32[] a, int aOrder, Complex32[] b, int bRows, int bColumns) { @@ -2288,8 +2371,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. - /// The number of columns in the B matrix. + /// The number of rows in the B matrix. + /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. public void CholeskySolveFactored(Complex32[] a, int aOrder, Complex32[] b, int bRows, int bColumns) { @@ -2300,11 +2383,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// 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. + /// 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. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public void QRFactor(Complex32[] r, Complex32[] q) + public void QRFactor(Complex32[] r, int rRows, int rColumns, Complex32[] q) { throw new NotImplementedException(); } @@ -2313,13 +2398,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// 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 + /// 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. /// 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. - public void QRFactor(Complex32[] r, Complex32[] q, Complex32[] work) + /// This is similar to the GEQRF and ORGQR LAPACK routines. + public void QRFactor(Complex32[] r, int rRows, int rColumns, Complex32[] q, Complex32[] work) { throw new NotImplementedException(); } @@ -2327,14 +2415,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// 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 + /// 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. /// The B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void QRSolve(int columnsOfB, Complex32[] r, Complex32[] q, Complex32[] b, Complex32[] x) + public void QRSolve(Complex32[] r, int rRows, int rColumns, Complex32[] q, Complex32[] b, int bColumns, Complex32[] x) { throw new NotImplementedException(); } @@ -2342,17 +2432,19 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// 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 + /// 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. /// 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. - public void QRSolve(int columnsOfB, Complex32[] r, Complex32[] q, Complex32[] b, Complex32[] x, Complex32[] work) + public void QRSolve(Complex32[] r, int rRows, int rColumns, Complex32[] q, Complex32[] b, int bColumns, Complex32[] x, Complex32[] work) { throw new NotImplementedException(); } @@ -2360,12 +2452,14 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// 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 Q matrix obtained by calling . + /// The R matrix obtained by calling . + /// 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. - public void QRSolveFactored(int columnsOfB, Complex32[] q, Complex32[] r, Complex32[] b, Complex32[] x) + public void QRSolveFactored(Complex32[] q, Complex32[] r, int rRows, int rColumns, Complex32[] b, int bColumns, Complex32[] x) { throw new NotImplementedException(); } @@ -2375,13 +2469,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// 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. - public void SingularValueDecomposition(bool computeVectors, Complex32[] a, Complex32[] s, Complex32[] u, Complex32[] vt) + public void SingularValueDecomposition(bool computeVectors, Complex32[] a, int aRows, int aColumns, Complex32[] s, Complex32[] u, Complex32[] vt) { throw new NotImplementedException(); } @@ -2391,6 +2487,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// 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. @@ -2400,7 +2498,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// 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. - public void SingularValueDecomposition(bool computeVectors, Complex32[] a, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] work) + public void SingularValueDecomposition(bool computeVectors, Complex32[] a, int aRows, int aColumns, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] work) { throw new NotImplementedException(); } @@ -2409,12 +2507,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// On exit, the solution matrix. - public void SvdSolve(Complex32[] a, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, Complex32[] x) + public void SvdSolve(Complex32[] a, int aRows, int aColumns, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, int bColumns, Complex32[] x) { throw new NotImplementedException(); } @@ -2423,15 +2524,18 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// 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. - public void SvdSolve(Complex32[] a, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, Complex32[] x, Complex32[] work) + public void SvdSolve(Complex32[] a, int aRows, int aColumns, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, int bColumns, Complex32[] x, Complex32[] work) { throw new NotImplementedException(); } @@ -2439,13 +2543,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Atlas /// /// 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 number of rows in the A matrix. + /// 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 B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void SvdSolveFactored(int columnsOfB, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, Complex32[] x) + public void SvdSolveFactored(int aRows, int aColumns, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, int bColumns, Complex32[] x) { throw new NotImplementedException(); } diff --git a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs index 01bed7d2..3beb89c7 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs @@ -3,9 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// // Copyright (c) 2009-2010 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 @@ -83,7 +81,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// Interface to linear algebra algorithms that work off 1-D arrays. /// /// Supported data types are double, single, , and . - public interface ILinearAlgebraProvider where T : struct + public interface ILinearAlgebraProvider + where T : struct { /*/// /// Queries the provider for the optimal, workspace block size @@ -93,7 +92,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// -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. /// @@ -210,9 +209,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The number of columns in the matrix. /// The value to scale the matrix. /// The c matrix. - void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, T alpha, T[] a, - int aRows, int aColumns, T[] b, int bRows, int bColumns, T beta, T[] c); - + void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, T alpha, T[] a, int aRows, int aColumns, T[] b, int bRows, int bColumns, T beta, T[] c); + /// /// Computes the LUP factorization of A. P*A = L*U. /// @@ -336,79 +334,93 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// /// 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. /// This is similar to the GEQRF and ORGQR LAPACK routines. - void QRFactor(T[] r, T[] q); + void QRFactor(T[] r, int rRows, int rColumns, 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. + /// 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. /// 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[] r, T[] q, T[] work); + void QRFactor(T[] r, int rRows, int rColumns, 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. + /// 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. /// The B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - void QRSolve(int columnsOfB, T[] r, T[] q, T[] b, T[] x); + void QRSolve(T[] r, int rRows, int rColumns, T[] q, T[] b, int bColumns, 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. + /// 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. /// 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. - void QRSolve(int columnsOfB, T[] r, T[] q, T[] b, T[] x, T[] work); + void QRSolve(T[] r, int rRows, int rColumns, T[] q, T[] b, int bColumns, 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 Q matrix obtained by calling . + /// The R matrix obtained by calling . + /// 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. - void QRSolveFactored(int columnsOfB, T[] q, T[] r, T[] b, T[] x); + void QRSolveFactored(T[] q, T[] r, int rRows, int rColumns, T[] b, int bColumns, 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 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. - void SingularValueDecomposition(bool computeVectors, T[] a, T[] s, T[] u, T[] vt); + void SingularValueDecomposition(bool computeVectors, T[] a, int aRows, int aColumns, 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. @@ -419,43 +431,51 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// 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); + void SingularValueDecomposition(bool computeVectors, T[] a, int aRows, int aColumns, 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// On exit, the solution matrix. - void SvdSolve(T[] a, T[] s, T[] u, T[] vt, T[] b, T[] x); + void SvdSolve(T[] a, int aRows, int aColumns, T[] s, T[] u, T[] vt, T[] b, int bColumns, 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// 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); + void SvdSolve(T[] a, int aRows, int aColumns, T[] s, T[] u, T[] vt, T[] b, int bColumns, 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 number of rows in the A matrix. + /// 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 B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - void SvdSolveFactored(int columnsOfB, T[] s, T[] u, T[] vt, T[] b, T[] x); + void SvdSolveFactored(int aRows, int aColumns, T[] s, T[] u, T[] vt, T[] b, int bColumns, T[] x); } } diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs index 31ca3024..d65226ad 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs @@ -254,6 +254,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra case Norm.FrobeniusNorm: break; } + return ret; } @@ -367,12 +368,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The number of columns in the matrix. /// The value to scale the matrix. /// The c matrix. - public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, - int aRows, int aColumns, double[] b, int bRows, int bColumns, double beta, double[] c) + public void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, int aRows, int aColumns, double[] b, int bRows, int bColumns, double beta, double[] c) { // Choose nonsensical values for the number of rows in c; fill them in depending // on the operations on a and b. - var cRows = -1; + int cRows; // First check some basic requirement on the parameters of the matrix multiplication. if (a == null) @@ -730,7 +730,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra ipiv[i] = i; } - var LUcolj = new double[order]; + var vLUcolj = new double[order]; // Outer loop. for (var j = 0; j < order; j++) @@ -738,10 +738,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra var indexj = j * order; var indexjj = indexj + j; -// Make a copy of the j-th column to localize references. + // Make a copy of the j-th column to localize references. for (var i = 0; i < order; i++) { - LUcolj[i] = data[indexj + i]; + vLUcolj[i] = data[indexj + i]; } // Apply previous transformations. @@ -752,17 +752,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra var s = 0.0; for (var k = 0; k < kmax; k++) { - s += data[k * order + i] * LUcolj[k]; + s += data[k * order + i] * vLUcolj[k]; } - data[indexj + i] = LUcolj[i] -= s; + data[indexj + i] = vLUcolj[i] -= s; } // Find pivot and exchange if necessary. var p = j; for (var i = j + 1; i < order; i++) { - if (Math.Abs(LUcolj[i]) > Math.Abs(LUcolj[p])) + if (Math.Abs(vLUcolj[i]) > Math.Abs(vLUcolj[p])) { p = i; } @@ -899,9 +899,9 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// /// Solves A*X=B for X using a previously factored A matrix. /// - /// The square, positive definite matrix A. Has to be different than . + /// The square, positive definite matrix A. Has to be different than . /// The number of rows and columns in A. - /// The B matrix. Has to be different than . + /// The B matrix. Has to be different than . /// The number of rows in the B matrix. /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. @@ -964,10 +964,325 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// /// 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. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public void QRFactor(double[] r, double[] q) + public void QRFactor(double[] r, int rRows, int rColumns, double[] q) + { + if (r == null) + { + throw new ArgumentNullException("r"); + } + + if (q == null) + { + throw new ArgumentNullException("q"); + } + + if (r.Length != rRows * rColumns) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + } + + if (q.Length != rRows * rRows) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + } + + var work = new double[rRows * rRows]; + QRFactor(r, rRows, rColumns, q, 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. + /// 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 void QRFactor(double[] r, int rRows, int rColumns, double[] q, 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 != rRows * rColumns) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + } + + if (q.Length != rRows * rRows) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + } + + if (work.Length < rRows * rRows) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "work"); + } + + CommonParallel.For(0, rRows, i => q[(i * rRows) + i] = 1.0); + + var minmn = Math.Min(rRows, rColumns); + for (var i = 0; i < minmn; i++) + { + GenerateColumn(work, r, rRows, i, rRows - 1, i); + ComputeQR(work, i, r, rRows, i, rRows - 1, i + 1, rColumns - 1); + } + + for (var i = minmn - 1; i >= 0; i--) + { + ComputeQR(work, i, q, rRows, i, rRows - 1, i, rRows - 1); + } + } + + #region QR Factor Helper functions + + /// + /// Perform calculation of Q or R + /// + /// Work arrat + /// Index of colunn in work array + /// Q or R matrices + /// The number of rows + /// The first row in + /// The last row + /// The first column + /// The last column + private static void ComputeQR(double[] work, int workIndex, double[] a, int rowCount, int rowStart, int rowEnd, int columnStart, int columnEnd) + { + if (rowStart > rowEnd || columnStart > columnEnd) + { + return; + } + + var vector = new double[columnEnd - columnStart + 1]; + for (var i = rowStart; i <= rowEnd; i++) + { + for (var j = columnStart; j <= columnEnd; j++) + { + vector[j - columnStart] += work[(workIndex * rowCount) + i - rowStart] * a[(j * rowCount) + i]; + } + } + + for (var i = rowStart; i <= rowEnd; i++) + { + for (var j = columnStart; j <= columnEnd; j++) + { + a[(j * rowCount) + i] -= work[(workIndex * rowCount) + i - rowStart] * vector[j - columnStart]; + } + } + } + + /// + /// Generate column from initial matrix to work array + /// + /// Work array + /// Initial matrix + /// The number of rows in matrix + /// The firts row + /// The last row + /// Column index + private static void GenerateColumn(double[] work, double[] a, int rowCount, int rowStart, int rowEnd, int column) + { + var tmp = column * rowCount; + var index = tmp + rowStart; + + CommonParallel.For( + rowStart, + rowEnd + 1, + i => + { + var iIndex = tmp + i; + work[iIndex - rowStart] = a[iIndex]; + a[iIndex] = 0.0; + }); + + var norm = 0.0; + for (var i = 0; i < rowEnd - rowStart + 1; ++i) + { + var iIndex = tmp + i; + norm += work[iIndex] * work[iIndex]; + } + + norm = Math.Sqrt(norm); + if (rowStart == rowEnd || norm == 0) + { + a[index] = -work[tmp]; + work[tmp] = Math.Sqrt(2.0); + return; + } + + var scale = 1.0 / norm; + if (work[tmp] < 0.0) + { + scale *= -1.0; + } + + a[index] = -1.0 / scale; + CommonParallel.For(0, rowEnd - rowStart + 1, i => work[tmp + i] *= scale); + work[tmp] += 1.0; + + var s = Math.Sqrt(1.0 / work[tmp]); + CommonParallel.For(0, rowEnd - rowStart + 1, i => work[tmp + i] *= s); + } + + #endregion + + /// + /// Solves A*X=B for X using 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. + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + public void QRSolve(double[] r, int rRows, int rColumns, double[] q, double[] b, int bColumns, double[] x) + { + if (r == null) + { + throw new ArgumentNullException("r"); + } + + if (q == null) + { + throw new ArgumentNullException("q"); + } + + if (b == null) + { + throw new ArgumentNullException("q"); + } + + if (x == null) + { + throw new ArgumentNullException("q"); + } + + if (r.Length != rRows * rColumns) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + } + + if (q.Length != rRows * rRows) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + } + + if (b.Length != rRows * bColumns) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (x.Length != rColumns * bColumns) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + } + + var work = new double[rRows * rRows]; + QRSolve(r, rRows, rColumns, q, b, bColumns, x, work); + } + + /// + /// Solves A*X=B for X using 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. + /// 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. + public void QRSolve(double[] r, int rRows, int rColumns, double[] q, double[] b, int bColumns, double[] x, double[] work) + { + if (r == null) + { + throw new ArgumentNullException("r"); + } + + if (q == null) + { + throw new ArgumentNullException("q"); + } + + if (b == null) + { + throw new ArgumentNullException("q"); + } + + if (x == null) + { + throw new ArgumentNullException("q"); + } + + if (r.Length != rRows * rColumns) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + } + + if (q.Length != rRows * rRows) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + } + + if (b.Length != rRows * bColumns) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (x.Length != rColumns * bColumns) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + } + + if (work.Length < rRows * rRows) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "work"); + } + + QRFactor(r, rRows, rColumns, q, work); + QRSolveFactored(q, r, rRows, rColumns, b, bColumns, x); + } + + /// + /// Solves A*X=B for X using a previously QR factored matrix. + /// + /// The Q matrix obtained by calling . + /// The R matrix obtained by calling . + /// 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. + public void QRSolveFactored(double[] q, double[] r, int rRows, int rColumns, double[] b, int bColumns, double[] x) { if (r == null) { @@ -979,234 +1294,1778 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentNullException("q"); } - // Matrix Q is square (m x m), where "m" is number of rows of initial matrix; - var rowCount = (int)Math.Sqrt(q.Length); - var columnCount = r.Length / rowCount; - - for (var i = 0; i < rowCount; i++) + if (b == null) + { + throw new ArgumentNullException("q"); + } + + if (x == null) + { + throw new ArgumentNullException("q"); + } + + if (r.Length != rRows * rColumns) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + } + + if (q.Length != rRows * rRows) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + } + + if (b.Length != rRows * bColumns) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (x.Length != rColumns * bColumns) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "x"); + } + + var sol = new double[b.Length]; + + // Copy B matrix to "sol", so B data will not be changed + Buffer.BlockCopy(b, 0, sol, 0, b.Length * Constants.SizeOfDouble); + + // Compute Y = transpose(Q)*B + var column = new double[rRows]; + for (var j = 0; j < bColumns; j++) + { + var jm = j * rRows; + CommonParallel.For(0, rRows, k => column[k] = sol[jm + k]); + CommonParallel.For( + 0, + rRows, + i => + { + var im = i * rRows; + sol[jm + i] = CommonParallel.Aggregate(0, rRows, k => q[im + k] * column[k]); + }); + } + + // Solve R*X = Y; + for (var k = rColumns - 1; k >= 0; k--) + { + var km = k * rRows; + for (var j = 0; j < bColumns; j++) + { + sol[(j * rRows) + k] /= r[km + k]; + } + + for (var i = 0; i < k; i++) + { + for (var j = 0; j < bColumns; j++) + { + var jm = j * rRows; + sol[jm + i] -= sol[jm + k] * r[km + i]; + } + } + } + + // Fill result matrix + CommonParallel.For( + 0, + rColumns, + row => + { + for (var col = 0; col < bColumns; col++) + { + x[(col * rColumns) + row] = sol[row + (col * rRows)]; + } + }); + } + + /// + /// 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. + public void SingularValueDecomposition(bool computeVectors, double[] a, int aRows, int aColumns, 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 != aRows * aRows) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); + } + + if (vt.Length != aColumns * aColumns) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); + } + + if (s.Length != Math.Min(aRows, aColumns)) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); + } + + // TODO: Actually "work = new double[aRows]" is acceptable size of work array. I set size proposed in method description + var work = new double[Math.Max((3 * Math.Min(aRows, aColumns)) + Math.Max(aRows, aColumns), 5 * Math.Min(aRows, aColumns))]; + SingularValueDecomposition(computeVectors, a, aRows, aColumns, 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. 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. + public void SingularValueDecomposition(bool computeVectors, double[] a, int aRows, int aColumns, 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 != aRows * aRows) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); + } + + if (vt.Length != aColumns * aColumns) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); + } + + if (s.Length != Math.Min(aRows, aColumns)) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); + } + + if (work.Length == 0) + { + throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work"); + } + + if (work.Length < aRows) + { + // TODO: Actually "work = new double[aRows]" is acceptable size of work array. I set size proposed in method description + work[0] = Math.Max((3 * Math.Min(aRows, aColumns)) + Math.Max(aRows, aColumns), 5 * Math.Min(aRows, aColumns)); + return; + } + + const int Maxiter = 1000; + + var e = new double[aColumns]; + var v = new double[vt.Length]; + var stemp = new double[Math.Min(aRows + 1, aColumns)]; + + int i, j, l, lp1; + + var cs = 0.0; + var sn = 0.0; + double t; + + var ncu = aRows; + + // Reduce matrix to bidiagonal form, storing the diagonal elements + // in "s" and the super-diagonal elements in "e". + var nct = Math.Min(aRows - 1, aColumns); + var nrt = Math.Max(0, Math.Min(aColumns - 2, aRows)); + var lu = Math.Max(nct, nrt); + + for (l = 0; l < lu; l++) + { + lp1 = l + 1; + if (l < nct) + { + // Compute the transformation for the l-th column and + // place the l-th diagonal in vector s[l]. + var l1 = l; + stemp[l] = Math.Sqrt(CommonParallel.Aggregate(l, aRows, i1 => (a[(l1 * aRows) + i1] * a[(l1 * aRows) + i1]))); + + if (stemp[l] != 0.0) + { + if (a[(l * aRows) + l] != 0.0) + { + stemp[l] = Math.Abs(stemp[l]) * (a[(l * aRows) + l] / Math.Abs(a[(l * aRows) + l])); + } + + // A part of column "l" of Matrix A from row "l" to end multiply by 1.0 / s[l] + for (i = l; i < aRows; i++) + { + a[(l * aRows) + i] = a[(l * aRows) + i] * (1.0 / stemp[l]); + } + + a[(l * aRows) + l] = 1.0 + a[(l * aRows) + l]; + } + + stemp[l] = -stemp[l]; + } + + for (j = lp1; j < aColumns; j++) + { + if (l < nct) + { + if (stemp[l] != 0.0) + { + // Apply the transformation. + t = 0.0; + for (i = l; i < aRows; i++) + { + t += a[(j * aRows) + i] * a[(l * aRows) + i]; + } + + t = -t / a[(l * aRows) + l]; + + for (var ii = l; ii < aRows; ii++) + { + a[(j * aRows) + ii] += t * a[(l * aRows) + ii]; + } + } + } + + // Place the l-th row of matrix into "e" for the + // subsequent calculation of the row transformation. + e[j] = a[(j * aRows) + l]; + } + + if (computeVectors && l < nct) + { + // Place the transformation in "u" for subsequent back multiplication. + for (i = l; i < aRows; i++) + { + u[(l * aRows) + i] = a[(l * aRows) + i]; + } + } + + if (l >= nrt) + { + continue; + } + + // Compute the l-th row transformation and place the l-th super-diagonal in e(l). + var enorm = 0.0; + for (i = lp1; i < e.Length; i++) + { + enorm += e[i] * e[i]; + } + + e[l] = Math.Sqrt(enorm); + if (e[l] != 0.0) + { + if (e[lp1] != 0.0) + { + e[l] = Math.Abs(e[l]) * (e[lp1] / Math.Abs(e[lp1])); + } + + // Scale vector "e" from "lp1" by 1.0 / e[l] + for (i = lp1; i < e.Length; i++) + { + e[i] = e[i] * (1.0 / e[l]); + } + + e[lp1] = 1.0 + e[lp1]; + } + + e[l] = -e[l]; + + if (lp1 < aRows && e[l] != 0.0) + { + // Apply the transformation. + for (i = lp1; i < aRows; i++) + { + work[i] = 0.0; + } + + for (j = lp1; j < aColumns; j++) + { + for (var ii = lp1; ii < aRows; ii++) + { + work[ii] += e[j] * a[(j * aRows) + ii]; + } + } + + for (j = lp1; j < aColumns; j++) + { + var ww = -e[j] / e[lp1]; + for (var ii = lp1; ii < aRows; ii++) + { + a[(j * aRows) + ii] += ww * work[ii]; + } + } + } + + if (!computeVectors) + { + continue; + } + + // Place the transformation in v for subsequent back multiplication. + for (i = lp1; i < aColumns; i++) + { + v[(l * aColumns) + i] = e[i]; + } + } + + // Set up the final bidiagonal matrix or order m. + var m = Math.Min(aColumns, aRows + 1); + var nctp1 = nct + 1; + var nrtp1 = nrt + 1; + if (nct < aColumns) + { + stemp[nctp1 - 1] = a[((nctp1 - 1) * aRows) + (nctp1 - 1)]; + } + + if (aRows < m) + { + stemp[m - 1] = 0.0; + } + + if (nrtp1 < m) + { + e[nrtp1 - 1] = a[((m - 1) * aRows) + (nrtp1 - 1)]; + } + + e[m - 1] = 0.0; + + // If required, generate "u". + if (computeVectors) + { + for (j = nctp1 - 1; j < ncu; j++) + { + for (i = 0; i < aRows; i++) + { + u[(j * aRows) + i] = 0.0; + } + + u[(j * aRows) + j] = 1.0; + } + + for (l = nct - 1; l >= 0; l--) + { + if (stemp[l] != 0.0) + { + for (j = l + 1; j < ncu; j++) + { + t = 0.0; + for (i = l; i < aRows; i++) + { + t += u[(j * aRows) + i] * u[(l * aRows) + i]; + } + + t = -t / u[(l * aRows) + l]; + + for (var ii = l; ii < aRows; ii++) + { + u[(j * aRows) + ii] += t * u[(l * aRows) + ii]; + } + } + + // A part of column "l" of matrix A from row "l" to end multiply by -1.0 + for (i = l; i < aRows; i++) + { + u[(l * aRows) + i] = u[(l * aRows) + i] * -1.0; + } + + u[(l * aRows) + l] = 1.0 + u[(l * aRows) + l]; + for (i = 0; i < l; i++) + { + u[(l * aRows) + i] = 0.0; + } + } + else + { + for (i = 0; i < aRows; i++) + { + u[(l * aRows) + i] = 0.0; + } + + u[(l * aRows) + l] = 1.0; + } + } + } + + // If it is required, generate v. + if (computeVectors) + { + for (l = aColumns - 1; l >= 0; l--) + { + lp1 = l + 1; + if (l < nrt) + { + if (e[l] != 0.0) + { + for (j = lp1; j < aColumns; j++) + { + t = 0.0; + for (i = lp1; i < aColumns; i++) + { + t += v[(j * aColumns) + i] * v[(l * aColumns) + i]; + } + + t = -t / v[(l * aColumns) + lp1]; + for (var ii = l; ii < aColumns; ii++) + { + v[(j * aColumns) + ii] += t * v[(l * aColumns) + ii]; + } + } + } + } + + for (i = 0; i < aColumns; i++) + { + v[(l * aColumns) + i] = 0.0; + } + + v[(l * aColumns) + l] = 1.0; + } + } + + // Transform "s" and "e" so that they are double + for (i = 0; i < m; i++) + { + double r; + if (stemp[i] != 0.0) + { + t = stemp[i]; + r = stemp[i] / t; + stemp[i] = t; + if (i < m - 1) + { + e[i] = e[i] / r; + } + + if (computeVectors) + { + // A part of column "i" of matrix U from row 0 to end multiply by r + for (j = 0; j < aRows; j++) + { + u[(i * aRows) + j] = u[(i * aRows) + j] * r; + } + } + } + + // Exit + if (i == m - 1) + { + break; + } + + if (e[i] == 0.0) + { + continue; + } + + t = e[i]; + r = t / e[i]; + e[i] = t; + stemp[i + 1] = stemp[i + 1] * r; + if (!computeVectors) + { + continue; + } + + // A part of column "i+1" of matrix VT from row 0 to end multiply by r + for (j = 0; j < aColumns; j++) + { + v[((i + 1) * aColumns) + j] = v[((i + 1) * aColumns) + j] * r; + } + } + + // Main iteration loop for the singular values. + var mn = m; + var iter = 0; + + while (m > 0) + { + // Quit if all the singular values have been found. + // If too many iterations have been performed throw exception. + if (iter >= Maxiter) + { + throw new ArgumentException(Resources.ConvergenceFailed); + } + + // This section of the program inspects for negligible elements in the s and e arrays, + // on completion the variables kase and l are set as follows: + // kase = 1: if mS[m] and e[l-1] are negligible and l < m + // kase = 2: if mS[l] is negligible and l < m + // kase = 3: if e[l-1] is negligible, l < m, and mS[l, ..., mS[m] are not negligible (qr step). + // kase = 4: if e[m-1] is negligible (convergence). + double ztest; + double test; + for (l = m - 2; l >= 0; l--) + { + test = Math.Abs(stemp[l]) + Math.Abs(stemp[l + 1]); + ztest = test + Math.Abs(e[l]); + if (ztest.AlmostEqualInDecimalPlaces(test, 15)) + { + e[l] = 0.0; + break; + } + } + + int kase; + if (l == m - 2) + { + kase = 4; + } + else + { + int ls; + for (ls = m - 1; ls > l; ls--) + { + test = 0.0; + if (ls != m - 1) + { + test = test + Math.Abs(e[ls]); + } + + if (ls != l + 1) + { + test = test + Math.Abs(e[ls - 1]); + } + + ztest = test + Math.Abs(stemp[ls]); + if (ztest.AlmostEqualInDecimalPlaces(test, 15)) + { + stemp[ls] = 0.0; + break; + } + } + + if (ls == l) + { + kase = 3; + } + else if (ls == m - 1) + { + kase = 1; + } + else + { + kase = 2; + l = ls; + } + } + + l = l + 1; + + // Perform the task indicated by kase. + int k; + double f; + switch (kase) + { + // Deflate negligible s[m]. + case 1: + f = e[m - 2]; + e[m - 2] = 0.0; + double t1; + for (var kk = l; kk < m - 1; kk++) + { + k = m - 2 - kk + l; + t1 = stemp[k]; + + Drotg(ref t1, ref f, ref cs, ref sn); + stemp[k] = t1; + if (k != l) + { + f = -sn * e[k - 1]; + e[k - 1] = cs * e[k - 1]; + } + + if (computeVectors) + { + // Rotate + for (i = 0; i < aColumns; i++) + { + var z = (cs * v[(k * aColumns) + i]) + (sn * v[((m - 1) * aColumns) + i]); + v[((m - 1) * aColumns) + i] = (cs * v[((m - 1) * aColumns) + i]) - (sn * v[(k * aColumns) + i]); + v[(k * aColumns) + i] = z; + } + } + } + + break; + + // Split at negligible s[l]. + case 2: + f = e[l - 1]; + e[l - 1] = 0.0; + for (k = l; k < m; k++) + { + t1 = stemp[k]; + Drotg(ref t1, ref f, ref cs, ref sn); + stemp[k] = t1; + f = -sn * e[k]; + e[k] = cs * e[k]; + if (computeVectors) + { + // Rotate + for (i = 0; i < aRows; i++) + { + var z = (cs * u[(k * aRows) + i]) + (sn * u[((l - 1) * aRows) + i]); + u[((l - 1) * aRows) + i] = (cs * u[((l - 1) * aRows) + i]) - (sn * u[(k * aRows) + i]); + u[(k * aRows) + i] = z; + } + } + } + + break; + + // Perform one qr step. + case 3: + +// calculate the shift. + var scale = 0.0; + scale = Math.Max(scale, Math.Abs(stemp[m - 1])); + scale = Math.Max(scale, Math.Abs(stemp[m - 2])); + scale = Math.Max(scale, Math.Abs(e[m - 2])); + scale = Math.Max(scale, Math.Abs(stemp[l])); + scale = Math.Max(scale, Math.Abs(e[l])); + var sm = stemp[m - 1] / scale; + var smm1 = stemp[m - 2] / scale; + var emm1 = e[m - 2] / scale; + var sl = stemp[l] / scale; + var el = e[l] / scale; + var b = (((smm1 + sm) * (smm1 - sm)) + (emm1 * emm1)) / 2.0; + var c = (sm * emm1) * (sm * emm1); + var shift = 0.0; + if (b != 0.0 || c != 0.0) + { + shift = Math.Sqrt((b * b) + c); + if (b < 0.0) + { + shift = -shift; + } + + shift = c / (b + shift); + } + + f = ((sl + sm) * (sl - sm)) + shift; + var g = sl * el; + + // Chase zeros + for (k = l; k < m - 1; k++) + { + Drotg(ref f, ref g, ref cs, ref sn); + if (k != l) + { + e[k - 1] = f; + } + + f = (cs * stemp[k]) + (sn * e[k]); + e[k] = (cs * e[k]) - (sn * stemp[k]); + g = sn * stemp[k + 1]; + stemp[k + 1] = cs * stemp[k + 1]; + if (computeVectors) + { + for (i = 0; i < aColumns; i++) + { + var z = (cs * v[(k * aColumns) + i]) + (sn * v[((k + 1) * aColumns) + i]); + v[((k + 1) * aColumns) + i] = (cs * v[((k + 1) * aColumns) + i]) - (sn * v[(k * aColumns) + i]); + v[(k * aColumns) + i] = z; + } + } + + Drotg(ref f, ref g, ref cs, ref sn); + stemp[k] = f; + f = (cs * e[k]) + (sn * stemp[k + 1]); + stemp[k + 1] = -(sn * e[k]) + (cs * stemp[k + 1]); + g = sn * e[k + 1]; + e[k + 1] = cs * e[k + 1]; + if (computeVectors && k < aRows) + { + for (i = 0; i < aRows; i++) + { + var z = (cs * u[(k * aRows) + i]) + (sn * u[((k + 1) * aRows) + i]); + u[((k + 1) * aRows) + i] = (cs * u[((k + 1) * aRows) + i]) - (sn * u[(k * aRows) + i]); + u[(k * aRows) + i] = z; + } + } + } + + e[m - 2] = f; + iter = iter + 1; + break; + + // Convergence + case 4: + +// Make the singular value positive + if (stemp[l] < 0.0) + { + stemp[l] = -stemp[l]; + if (computeVectors) + { + // A part of column "l" of matrix VT from row 0 to end multiply by -1 + for (i = 0; i < aColumns; i++) + { + v[(l * aColumns) + i] = v[(l * aColumns) + i] * -1.0; + } + } + } + + // Order the singular value. + while (l != mn - 1) + { + if (stemp[l] >= stemp[l + 1]) + { + break; + } + + t = stemp[l]; + stemp[l] = stemp[l + 1]; + stemp[l + 1] = t; + if (computeVectors && l < aColumns) + { + // Swap columns l, l + 1 + for (i = 0; i < aColumns; i++) + { + var z = v[(l * aColumns) + i]; + v[(l * aColumns) + i] = v[((l + 1) * aColumns) + i]; + v[((l + 1) * aColumns) + i] = z; + } + } + + if (computeVectors && l < aRows) + { + // Swap columns l, l + 1 + for (i = 0; i < aRows; i++) + { + var z = u[(l * aRows) + i]; + u[(l * aRows) + i] = u[((l + 1) * aRows) + i]; + u[((l + 1) * aRows) + i] = z; + } + } + + l = l + 1; + } + + iter = 0; + m = m - 1; + break; + } + } + + if (computeVectors) + { + // Finally transpose "v" to get "vt" matrix + for (i = 0; i < aColumns; i++) + { + for (j = 0; j < aColumns; j++) + { + vt[(j * aColumns) + i] = v[(i * aColumns) + j]; + } + } + } + + // Copy stemp to s with size adjustment. We are using ported copy of linpack's svd code and it uses + // 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(aRows, aColumns) * 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] = aRows; + } + + public void SingularValueDecompositionWithCommonParallel(bool computeVectors, double[] a, int aRows, int aColumns, 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 (u.Length != aRows * aRows) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); + } + + if (vt.Length != aColumns * aColumns) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); + } + + if (s.Length != Math.Min(aRows, aColumns)) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); + } + + if (work == null) + { + throw new ArgumentNullException("work"); + } + + if (work.Length == 0) + { + throw new ArgumentException(Resources.ArgumentSingleDimensionArray); + } + + if (work.Length < aRows) + { + work[0] = Math.Max(3 * Math.Min(aRows, aColumns) + Math.Max(aRows, aColumns), 5 * Math.Min(aRows, aColumns)); + return; + } + + const int Maxiter = 1000; + + var e = new double[aColumns]; + var v = new double[vt.Length]; + var stemp = new double[Math.Min(aRows + 1, aColumns)]; + + int i, j, l, lp1; + + var cs = 0.0; + var sn = 0.0; + double t; + + var ncu = aRows; + + // Reduce matrix to bidiagonal form, storing the diagonal elements + // in "s" and the super-diagonal elements in "e". + var nct = Math.Min(aRows - 1, aColumns); + var nrt = Math.Max(0, Math.Min(aColumns - 2, aRows)); + var lu = Math.Max(nct, nrt); + + for (l = 0; l < lu; l++) + { + lp1 = l + 1; + var lCopy = l; + if (l < nct) + { + // Compute the transformation for the l-th column and + // place the l-th diagonal in vector s[l]. + stemp[l] = Math.Sqrt(CommonParallel.Aggregate(lCopy, aRows, row => (a[(lCopy * aRows) + row] * a[(lCopy * aRows) + row]))); + + if (stemp[l] != 0.0) + { + if (a[(l * aRows) + l] != 0.0) + { + stemp[l] = Math.Abs(stemp[l]) * (a[(l * aRows) + l] / Math.Abs(a[(l * aRows) + l])); + } + + // A part of column "l" of Matrix A from row "l" to end multiply by 1.0 / s[l] + for (i = l; i < aRows; i++) + { + a[(l * aRows) + i] = a[(l * aRows) + i] * (1.0 / stemp[l]); + } + + a[(l * aRows) + l] = 1.0 + a[(l * aRows) + l]; + } + + stemp[l] = -stemp[l]; + } + + for (j = lp1; j < aColumns; j++) + { + var jcopy = j; + + if (l < nct) + { + if (stemp[l] != 0.0) + { + // Apply the transformation. + t = CommonParallel.Aggregate(lCopy, aRows, row => a[(jcopy * aRows) + row] * a[(lCopy * aRows) + row]); + t = -t / a[(l * aRows) + l]; + for (var ii = l; ii < aRows; ii++) + { + a[(j * aRows) + ii] += t * a[(l * aRows) + ii]; + } + } + } + + // Place the l-th row of matrix into "e" for the + // subsequent calculation of the row transformation. + e[j] = a[(j * aRows) + l]; + } + + if (computeVectors && l < nct) + { + // Place the transformation in "u" for subsequent back multiplication. + CommonParallel.For(lCopy, aRows, row => u[(lCopy * aRows) + row] = a[(lCopy * aRows) + row]); + } + + if (l >= nrt) + { + continue; + } + + // Compute the l-th row transformation and place the l-th super-diagonal in e(l). + e[l] = Math.Sqrt(CommonParallel.Aggregate(lp1, e.Length, index => e[index] * e[index])); + if (e[l] != 0.0) + { + if (e[lp1] != 0.0) + { + e[l] = Math.Abs(e[l]) * (e[lp1] / Math.Abs(e[lp1])); + } + + // Scale vector "e" from "lp1" by 1.0 / e[l] + // lp1 > l, so we may use CommonParallel here + CommonParallel.For(lp1, e.Length, index => e[index] = e[index] * (1.0 / e[lCopy])); + e[lp1] = 1.0 + e[lp1]; + } + + e[l] = -e[l]; + + if (lp1 < aRows && e[l] != 0.0) + { + // Apply the transformation. + for (i = lp1; i < aRows; i++) + { + work[i] = 0.0; + } + + for (j = lp1; j < aColumns; j++) + { + for (var ii = lp1; ii < aRows; ii++) + { + work[ii] += e[j] * a[(j * aRows) + ii]; + } + } + + for (j = lp1; j < aColumns; j++) + { + var ww = -e[j] / e[lp1]; + for (var ii = lp1; ii < aRows; ii++) + { + a[(j * aRows) + ii] += ww * work[ii]; + } + } + } + + if (!computeVectors) + { + continue; + } + + // Place the transformation in v for subsequent back multiplication. + CommonParallel.For(lp1, aColumns, index => v[(lCopy * aColumns) + index] = e[index]); + } + + // Set up the final bidiagonal matrix or order m. + var m = Math.Min(aColumns, aRows + 1); + var nctp1 = nct + 1; + var nrtp1 = nrt + 1; + if (nct < aColumns) + { + stemp[nctp1 - 1] = a[(nctp1 - 1) * aRows + (nctp1 - 1)]; + } + + if (aRows < m) + { + stemp[m - 1] = 0.0; + } + + if (nrtp1 < m) + { + e[nrtp1 - 1] = a[(m - 1) * aRows + (nrtp1 - 1)]; + } + + e[m - 1] = 0.0; + + // If required, generate "u". + if (computeVectors) + { + CommonParallel.For(nctp1 - 1, ncu, indexCol => + { + CommonParallel.For(0, aRows, indexRow => u[(indexCol * aRows) + indexRow] = 0.0); + u[(indexCol * aRows) + indexCol] = 1.0; + }); + + for (l = nct - 1; l >= 0; l--) + { + var lCopy = l; + + if (stemp[l] != 0.0) + { + for (j = l + 1; j < ncu; j++) + { + var jCopy = j; + t = CommonParallel.Aggregate(lCopy, aRows, index => u[(jCopy * aRows) + index] * u[(lCopy * aRows) + index]); + t = -t / u[(l * aRows) + l]; + + for (var ii = l; ii < aRows; ii++) + { + u[(j * aRows) + ii] += t * u[(l * aRows) + ii]; + } + } + + // A part of column "l" of matrix A from row "l" to end multiply by -1.0 + CommonParallel.For(lCopy, aRows, index => u[(lCopy * aRows) + index] *= -1.0); + u[(l * aRows) + l] = 1.0 + u[(l * aRows) + l]; + if (lCopy != 0) + { + CommonParallel.For(0, lCopy, index => u[(lCopy * aRows) + index] = 0.0); + } + } + else + { + CommonParallel.For(0, aRows, index => u[(lCopy * aRows) + index] = 0.0); + u[(l * aRows) + l] = 1.0; + } + } + } + + // If it is required, generate v. + if (computeVectors) + { + for (l = aColumns - 1; l >= 0; l--) + { + lp1 = l + 1; + var lCopy = l; + + if (l < nrt) + { + if (e[l] != 0.0) + { + for (j = lp1; j < aColumns; j++) + { + var jCopy = j; + t = CommonParallel.Aggregate(lp1, aColumns, index => v[(jCopy * aColumns) + index] * v[(lCopy * aColumns) + index]); + t = -t / v[(l * aColumns) + lp1]; + for (var ii = l; ii < aColumns; ii++) + { + v[(j * aColumns) + ii] += t * v[(l * aColumns) + ii]; + } + } + } + } + + CommonParallel.For(0, aColumns, index => v[(lCopy * aColumns) + index] = 0.0); + v[(l * aColumns) + l] = 1.0; + } + } + + // Transform "s" and "e" so that they are double + for (i = 0; i < m; i++) + { + double r; + if (stemp[i] != 0.0) + { + t = stemp[i]; + r = stemp[i] / t; + stemp[i] = t; + if (i < m - 1) + { + e[i] = e[i] / r; + } + + if (computeVectors) + { + // A part of column "i" of matrix U from row 0 to end multiply by r + for (j = 0; j < aRows; j++) + { + u[(i * aRows) + j] = u[(i * aRows) + j] * r; + } + } + } + + // Exit + if (i == m - 1) + { + break; + } + + if (e[i] == 0.0) + { + continue; + } + + t = e[i]; + r = t / e[i]; + e[i] = t; + stemp[i + 1] = stemp[i + 1] * r; + if (!computeVectors) + { + continue; + } + + // A part of column "i+1" of matrix VT from row 0 to end multiply by r + for (j = 0; j < aColumns; j++) + { + v[((i + 1) * aColumns) + j] = v[((i + 1) * aColumns) + j] * r; + } + } + + // Main iteration loop for the singular values. + var mn = m; + var iter = 0; + + while (m > 0) + { + // Quit if all the singular values have been found. + // If too many iterations have been performed throw exception. + if (iter >= Maxiter) + { + throw new ArgumentException(Resources.ConvergenceFailed); + } + + // This section of the program inspects for negligible elements in the s and e arrays. + // On completion the variables kase and l are set as follows. + + // kase = 1: if mS[m] and e[l-1] are negligible and l < m + // kase = 2: if mS[l] is negligible and l < m + // kase = 3: if e[l-1] is negligible, l < m, and mS[l, ..., mS[m] are not negligible (qr step). + // kase = 4: if e[m-1] is negligible (convergence). + double ztest; + double test; + for (l = m - 2; l >= 0; l--) + { + test = Math.Abs(stemp[l]) + Math.Abs(stemp[l + 1]); + ztest = test + Math.Abs(e[l]); + if (ztest.AlmostEqualInDecimalPlaces(test, 15)) + { + e[l] = 0.0; + break; + } + } + + int kase; + if (l == m - 2) + { + kase = 4; + } + else + { + int ls; + for (ls = m - 1; ls > l; ls--) + { + test = 0.0; + if (ls != m - 1) + { + test = test + Math.Abs(e[ls]); + } + + if (ls != l + 1) + { + test = test + Math.Abs(e[ls - 1]); + } + + ztest = test + Math.Abs(stemp[ls]); + if (ztest.AlmostEqualInDecimalPlaces(test, 15)) + { + stemp[ls] = 0.0; + break; + } + } + + if (ls == l) + { + kase = 3; + } + else if (ls == m - 1) + { + kase = 1; + } + else + { + kase = 2; + l = ls; + } + } + + l = l + 1; + + // Perform the task indicated by kase. + int k; + double f; + switch (kase) + { + // Deflate negligible s[m]. + case 1: + f = e[m - 2]; + e[m - 2] = 0.0; + double t1; + for (var kk = l; kk < m - 1; kk++) + { + k = m - 2 - kk + l; + t1 = stemp[k]; + + Drotg(ref t1, ref f, ref cs, ref sn); + stemp[k] = t1; + if (k != l) + { + f = -sn * e[k - 1]; + e[k - 1] = cs * e[k - 1]; + } + + if (computeVectors) + { + // Rotate + for (i = 0; i < aColumns; i++) + { + var z = cs * v[(k * aColumns) + i] + sn * v[((m - 1) * aColumns) + i]; + v[((m - 1) * aColumns) + i] = cs * v[((m - 1) * aColumns) + i] - sn * v[(k * aColumns) + i]; + v[(k * aColumns) + i] = z; + } + } + } + + break; + + // Split at negligible s[l]. + case 2: + f = e[l - 1]; + e[l - 1] = 0.0; + for (k = l; k < m; k++) + { + t1 = stemp[k]; + Drotg(ref t1, ref f, ref cs, ref sn); + stemp[k] = t1; + f = -sn * e[k]; + e[k] = cs * e[k]; + if (computeVectors) + { + // Rotate + for (i = 0; i < aRows; i++) + { + var z = cs * u[(k * aRows) + i] + sn * u[((l - 1) * aRows) + i]; + u[((l - 1) * aRows) + i] = cs * u[((l - 1) * aRows) + i] - sn * u[(k * aRows) + i]; + u[(k * aRows) + i] = z; + } + } + } + + break; + + // Perform one qr step. + case 3: + +// calculate the shift. + var scale = 0.0; + scale = Math.Max(scale, Math.Abs(stemp[m - 1])); + scale = Math.Max(scale, Math.Abs(stemp[m - 2])); + scale = Math.Max(scale, Math.Abs(e[m - 2])); + scale = Math.Max(scale, Math.Abs(stemp[l])); + scale = Math.Max(scale, Math.Abs(e[l])); + var sm = stemp[m - 1] / scale; + var smm1 = stemp[m - 2] / scale; + var emm1 = e[m - 2] / scale; + var sl = stemp[l] / scale; + var el = e[l] / scale; + var b = ((smm1 + sm) * (smm1 - sm) + (emm1 * emm1)) / 2.0; + var c = (sm * emm1) * (sm * emm1); + var shift = 0.0; + if (b != 0.0 || c != 0.0) + { + shift = Math.Sqrt((b * b) + c); + if (b < 0.0) + { + shift = -shift; + } + + shift = c / (b + shift); + } + + f = (sl + sm) * (sl - sm) + shift; + var g = sl * el; + + // Chase zeros + for (k = l; k < m - 1; k++) + { + Drotg(ref f, ref g, ref cs, ref sn); + if (k != l) + { + e[k - 1] = f; + } + + f = cs * stemp[k] + sn * e[k]; + e[k] = cs * e[k] - sn * stemp[k]; + g = sn * stemp[k + 1]; + stemp[k + 1] = cs * stemp[k + 1]; + if (computeVectors) + { + for (i = 0; i < aColumns; i++) + { + var z = cs * v[k * aColumns + i] + sn * v[(k + 1) * aColumns + i]; + v[(k + 1) * aColumns + i] = cs * v[(k + 1) * aColumns + i] - sn * v[k * aColumns + i]; + v[k * aColumns + i] = z; + } + } + + Drotg(ref f, ref g, ref cs, ref sn); + stemp[k] = f; + f = cs * e[k] + sn * stemp[k + 1]; + stemp[k + 1] = -sn * e[k] + cs * stemp[k + 1]; + g = sn * e[k + 1]; + e[k + 1] = cs * e[k + 1]; + if (computeVectors && k < aRows) + { + for (i = 0; i < aRows; i++) + { + var z = cs * u[k * aRows + i] + sn * u[(k + 1) * aRows + i]; + u[(k + 1) * aRows + i] = cs * u[(k + 1) * aRows + i] - sn * u[k * aRows + i]; + u[k * aRows + i] = z; + } + } + } + + e[m - 2] = f; + iter = iter + 1; + break; + + // Convergence + case 4: + +// Make the singular value positive + if (stemp[l] < 0.0) + { + stemp[l] = -stemp[l]; + if (computeVectors) + { + // A part of column "l" of matrix VT from row 0 to end multiply by -1 + for (i = 0; i < aColumns; i++) + { + v[(l * aColumns) + i] = v[(l * aColumns) + i] * -1.0; + } + } + } + + // Order the singular value. + while (l != mn - 1) + { + if (stemp[l] >= stemp[l + 1]) + { + break; + } + + t = stemp[l]; + stemp[l] = stemp[l + 1]; + stemp[l + 1] = t; + if (computeVectors && l < aColumns) + { + // Swap columns l, l + 1 + for (i = 0; i < aColumns; i++) + { + var z = v[l * aColumns + i]; + v[l * aColumns + i] = v[(l + 1) * aColumns + i]; + v[(l + 1) * aColumns + i] = z; + } + } + + if (computeVectors && l < aRows) + { + // Swap columns l, l + 1 + for (i = 0; i < aRows; i++) + { + var z = u[l * aRows + i]; + u[l * aRows + i] = u[(l + 1) * aRows + i]; + u[(l + 1) * aRows + i] = z; + } + } + + l = l + 1; + } + + iter = 0; + m = m - 1; + break; + } + } + + if (computeVectors) + { + // Finally transpose "v" to get "vt" matrix + for (i = 0; i < aColumns; i++) + { + for (j = 0; j < aColumns; j++) + { + vt[j * aColumns + i] = v[i * aColumns + j]; + } + } + } + + // Copy stemp to s with size adjustment. We are using ported copy of linpack's svd code and it uses + // 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(aRows, aColumns) * 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] = aRows; + } + + /// + /// Сonstruct givens plane rotation + /// + /// + /// + /// + /// + private static void Drotg(ref double da, ref double db, ref double c, ref double s) + { + // Сonstruct givens plane rotation. + // jack dongarra, linpack, 3/11/78. + double r, z; + + var roe = db; + var absda = Math.Abs(da); + var absdb = Math.Abs(db); + if (absda > absdb) + { + roe = da; + } + + var scale = absda + absdb; + if (scale == 0.0) + { + c = 1.0; + s = 0.0; + r = 0.0; + z = 0.0; + } + else + { + var sda = da / scale; + var sdb = db / scale; + r = scale * Math.Sqrt((sda * sda) + (sdb * sdb)); + if (roe < 0.0) + { + r = -r; + } + + c = da / r; + s = db / r; + z = 1.0; + if (absda > absdb) + { + z = s; + } + + if (absdb >= absda && c != 0.0) + { + z = 1.0 / c; + } + } + + da = r; + db = z; + return; + } + + /// + /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. + /// 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. + /// The number of columns of B. + /// On exit, the solution matrix. + public void SvdSolve(double[] a, int aRows, int aColumns, double[] s, double[] u, double[] vt, double[] b, int bColumns, double[] x) + { + 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 (b == null) + { + throw new ArgumentNullException("b"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (u.Length != aRows * aRows) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); + } + + if (vt.Length != aColumns * aColumns) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); + } + + if (s.Length != Math.Min(aRows, aColumns)) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); + } + + if (b.Length != aRows * bColumns) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (x.Length != aColumns * bColumns) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + // TODO: Actually "work = new double[aRows]" is acceptable size of work array. I set size proposed in method description + var work = new double[Math.Max((3 * Math.Min(aRows, aColumns)) + Math.Max(aRows, aColumns), 5 * Math.Min(aRows, aColumns))]; + SvdSolve(a, aRows, aColumns, s, u, vt, b, bColumns, x, 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 number of rows in the A matrix. + /// The number of columns in the A matrix. + /// 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. + /// The number of columns of B. + /// 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. + public void SvdSolve(double[] a, int aRows, int aColumns, double[] s, double[] u, double[] vt, double[] b, int bColumns, double[] x, double[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (s == null) { - q[i * rowCount + i] = 1.0; + throw new ArgumentNullException("s"); } - var minmn = Math.Min(rowCount, columnCount); - var u = new double[minmn][]; - - for (var i = 0; i < minmn; i++) + if (u == null) { - u[i] = GenerateColumn(r, rowCount, i, rowCount - 1, i); - UA(u[i], r, rowCount, i, rowCount - 1, i + 1, columnCount - 1); + throw new ArgumentNullException("u"); } - for (var i = minmn - 1; i >= 0; i--) + if (vt == null) { - UA(u[i], q, rowCount, i, rowCount - 1, i, rowCount - 1); + throw new ArgumentNullException("vt"); } - } - - public void QRFactor(double[] r, double[] q, double[] work) - { - throw new NotImplementedException(); - } - - #region QR Factor Helper functions - private static void UA(double[] u, double[] A, int rowCount, int rowStart, int rowEnd, int columnStart, int columnEnd) - { - if (rowStart > rowEnd || columnStart > columnEnd) + if (b == null) { - return; + throw new ArgumentNullException("b"); } - var vector = new double[columnEnd - columnStart + 1]; - for (var j = columnStart; j <= columnEnd; j++) + if (x == null) { - vector[j - columnStart] = 0.0; + throw new ArgumentNullException("x"); } - for (var i = rowStart; i <= rowEnd; i++) + if (u.Length != aRows * aRows) { - for (var j = columnStart; j <= columnEnd; j++) - { - vector[j - columnStart] = vector[j - columnStart] + u[i - rowStart] * A[j * rowCount + i]; - } + throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); } - for (var i = rowStart; i <= rowEnd; i++) + if (vt.Length != aColumns * aColumns) { - for (var j = columnStart; j <= columnEnd; j++) - { - A[j * rowCount + i] = A[j * rowCount + i] - u[i - rowStart] * vector[j - columnStart]; - } + throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); } - } - - private static double[] GenerateColumn(double[] A, int rowCount, int rowStart, int rowEnd, int column) - { - var u = new double[rowEnd - rowStart + 1]; - - var tmp = column * rowCount; - var index = tmp + rowStart; - for (var i = rowStart; i <= rowEnd; i++) + if (s.Length != Math.Min(aRows, aColumns)) { - u[i - rowStart] = A[tmp + i]; - A[tmp + i] = 0.0; + throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); } - var norm = u.Sum(t => t * t); - norm = Math.Sqrt(norm); - - if (rowStart == rowEnd || norm == 0) + if (b.Length != aRows * bColumns) { - A[index] = -u[0]; - u[0] = Math.Sqrt(2.0); - return u; + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - var scale = 1.0 / norm; - if (u[0] < 0.0) + if (x.Length != aColumns * bColumns) { - scale *= -1.0; + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - A[index] = -1.0 / scale; - for (var i = 0; i < u.Length; i++) + if (work.Length == 0) { - u[i] *= scale; + throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work"); } - u[0] += 1.0; - var s = Math.Sqrt(1.0 / u[0]); - for (var i = 0; i < u.Length; i++) + if (work.Length < aRows) { - u[i] *= s; + // TODO: Actually "work = new double[aRows]" is acceptable size of work array. I set size proposed in method description + work[0] = Math.Max((3 * Math.Min(aRows, aColumns)) + Math.Max(aRows, aColumns), 5 * Math.Min(aRows, aColumns)); + return; } - return u; - } - - #endregion - - public void QRSolve(int columnsOfB, double[] r, double[] q, double[] b, double[] x) - { - throw new NotImplementedException(); - } - - public void QRSolve(int columnsOfB, double[] r, double[] q, double[] b, double[] x, double[] work) - { - throw new NotImplementedException(); + SingularValueDecomposition(true, a, aRows, aColumns, s, u, vt, work); + SvdSolveFactored(aRows, aColumns, s, u, vt, b, bColumns, x); } /// - /// Solves A*X=B for X using a previously QR factored matrix. + /// Solves A*X=B for X using a previously SVD decomposed matrix. /// - /// The number of columns of B. - /// The Q matrix obtained by calling . - /// The R matrix obtained by calling . + /// The number of rows in the A matrix. + /// 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 B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void QRSolveFactored(int columnsOfB, double[] q, double[] r, double[] b, double[] x) + public void SvdSolveFactored(int aRows, int aColumns, double[] s, double[] u, double[] vt, double[] b, int bColumns, double[] x) { - if (r == null) + if (s == null) { - throw new ArgumentNullException("r"); + throw new ArgumentNullException("s"); } - if (q == null) + if (u == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("u"); + } + + if (vt == null) + { + throw new ArgumentNullException("vt"); } if (b == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("b"); } if (x == null) { - throw new ArgumentNullException("q"); + throw new ArgumentNullException("x"); } - if (b.Length != x.Length) + if (u.Length != aRows * aRows) { - throw new ArgumentException(Resources.ArgumentVectorsSameLength); + throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); } - // Matrix Q is square (m x m), where "m" is number of rows of initial matrix; - var rowCount = (int)Math.Sqrt(q.Length); - var columnCount = r.Length / rowCount; - - // Copy B matrix to result, so B data will not be changed - Buffer.BlockCopy(b, 0, x, 0, b.Length * Constants.SizeOfDouble); + if (vt.Length != aColumns * aColumns) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); + } - // Compute Y = transpose(Q)*B - var column = new double[rowCount]; - for (var j = 0; j < columnsOfB; j++) + if (s.Length != Math.Min(aRows, aColumns)) { - var jm = j * rowCount; - for (var k = 0; k < rowCount; k++) - { - column[k] = x[jm + k]; - } + throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); + } - for (var i = 0; i < rowCount; i++) - { - double s = 0; - var im = i * rowCount; - for (var k = 0; k < rowCount; k++) - { - s += q[im + k] * column[k]; - } + if (b.Length != aRows * bColumns) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } - x[jm + i] = s; - } + if (x.Length != aColumns * bColumns) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); } - // Solve R*X = Y; - for (var k = columnCount - 1; k >= 0; k--) + var mn = Math.Min(aRows, aColumns); + var tmp = new double[aColumns]; + + for (var k = 0; k < bColumns; k++) { - var km = k * rowCount; - for (var j = 0; j < columnsOfB; j++) + for (var j = 0; j < aColumns; j++) { - x[j * rowCount + k] /= r[km + k]; + double value = 0; + if (j < mn) + { + for (var i = 0; i < aRows; i++) + { + value += u[(j * aRows) + i] * b[(k * aRows) + i]; + } + + value /= s[j]; + } + + tmp[j] = value; } - for (var i = 0; i < k; i++) + for (var j = 0; j < aColumns; j++) { - for (var j = 0; j < columnsOfB; j++) + double value = 0; + for (var i = 0; i < aColumns; i++) { - var jm = j * rowCount; - x[jm + i] -= x[jm + k] * r[km + i]; + value += vt[(j * aColumns) + i] * tmp[i]; } + + x[(k * aColumns) + j] = value; } } } - public void SingularValueDecomposition(bool computeVectors, double[] a, double[] s, double[] u, double[] vt) - { - throw new NotImplementedException(); - } - - public void SingularValueDecomposition(bool computeVectors, double[] a, double[] s, double[] u, double[] vt, double[] work) - { - throw new NotImplementedException(); - } - - public void SvdSolve(double[] a, double[] s, double[] u, double[] vt, double[] b, double[] x) - { - throw new NotImplementedException(); - } - - public void SvdSolve(double[] a, double[] s, double[] u, double[] vt, double[] b, double[] x, double[] work) - { - throw new NotImplementedException(); - } - - public void SvdSolveFactored(int columnsOfB, double[] s, double[] u, double[] vt, double[] b, double[] x) - { - throw new NotImplementedException(); - } - #endregion #region ILinearAlgebraProvider Members @@ -1976,52 +3835,179 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new NotImplementedException(); } - public void QRFactor(float[] r, float[] 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. + /// 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. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + public void QRFactor(float[] r, int rRows, int rColumns, float[] q) { throw new NotImplementedException(); } - public void QRFactor(float[] r, float[] q, float[] 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. + /// 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 void QRFactor(float[] r, int rRows, int rColumns, float[] q, float[] work) { throw new NotImplementedException(); } - public void QRSolve(int columnsOfB, float[] r, float[] q, float[] b, float[] x) + /// + /// Solves A*X=B for X using 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. + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + public void QRSolve(float[] r, int rRows, int rColumns, float[] q, float[] b, int bColumns, float[] x) { throw new NotImplementedException(); } - public void QRSolve(int columnsOfB, float[] r, float[] q, float[] b, float[] x, float[] work) + /// + /// Solves A*X=B for X using 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. + /// 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. + public void QRSolve(float[] r, int rRows, int rColumns, float[] q, float[] b, int bColumns, float[] x, float[] work) { throw new NotImplementedException(); } - public void QRSolveFactored(int columnsOfB, float[] q, float[] r, float[] b, float[] x) + /// + /// Solves A*X=B for X using a previously QR factored matrix. + /// + /// The Q matrix obtained by calling . + /// The R matrix obtained by calling . + /// 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. + public void QRSolveFactored(float[] q, float[] r, int rRows, int rColumns, float[] b, int bColumns, float[] x) { throw new NotImplementedException(); } - public void SingularValueDecomposition(bool computeVectors, float[] a, float[] s, float[] u, float[] 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. + /// This is equivalent to the GESVD LAPACK routine. + public void SingularValueDecomposition(bool computeVectors, float[] a, int aRows, int aColumns, float[] s, float[] u, float[] vt) { throw new NotImplementedException(); } - public void SingularValueDecomposition(bool computeVectors, float[] a, float[] s, float[] u, float[] vt, float[] 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. 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. + public void SingularValueDecomposition(bool computeVectors, float[] a, int aRows, int aColumns, float[] s, float[] u, float[] vt, float[] work) { throw new NotImplementedException(); } - public void SvdSolve(float[] a, float[] s, float[] u, float[] vt, float[] b, float[] 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 number of rows in the A matrix. + /// The number of columns in the A matrix. + /// 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. + /// The number of columns of B. + /// On exit, the solution matrix. + public void SvdSolve(float[] a, int aRows, int aColumns, float[] s, float[] u, float[] vt, float[] b, int bColumns, float[] x) { throw new NotImplementedException(); } - public void SvdSolve(float[] a, float[] s, float[] u, float[] vt, float[] b, float[] x, float[] 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 number of rows in the A matrix. + /// The number of columns in the A matrix. + /// 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. + /// The number of columns of B. + /// 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. + public void SvdSolve(float[] a, int aRows, int aColumns, float[] s, float[] u, float[] vt, float[] b, int bColumns, float[] x, float[] work) { throw new NotImplementedException(); } - public void SvdSolveFactored(int columnsOfB, float[] s, float[] u, float[] vt, float[] b, float[] x) + /// + /// Solves A*X=B for X using a previously SVD decomposed matrix. + /// + /// The number of rows in the A matrix. + /// 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 B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + public void SvdSolveFactored(int aRows, int aColumns, float[] s, float[] u, float[] vt, float[] b, int bColumns, float[] x) { throw new NotImplementedException(); } @@ -2767,52 +4753,179 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new NotImplementedException(); } - public void QRFactor(Complex[] r, Complex[] 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. + /// 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. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + public void QRFactor(Complex[] r, int rRows, int rColumns, Complex[] q) { throw new NotImplementedException(); } - public void QRFactor(Complex[] r, Complex[] q, Complex[] 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. + /// 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 void QRFactor(Complex[] r, int rRows, int rColumns, Complex[] q, Complex[] work) { throw new NotImplementedException(); } - public void QRSolve(int columnsOfB, Complex[] r, Complex[] q, Complex[] b, Complex[] x) + /// + /// Solves A*X=B for X using 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. + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + public void QRSolve(Complex[] r, int rRows, int rColumns, Complex[] q, Complex[] b, int bColumns, Complex[] x) { throw new NotImplementedException(); } - public void QRSolve(int columnsOfB, Complex[] r, Complex[] q, Complex[] b, Complex[] x, Complex[] work) + /// + /// Solves A*X=B for X using 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. + /// 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. + public void QRSolve(Complex[] r, int rRows, int rColumns, Complex[] q, Complex[] b, int bColumns, Complex[] x, Complex[] work) { throw new NotImplementedException(); } - public void QRSolveFactored(int columnsOfB, Complex[] q, Complex[] r, Complex[] b, Complex[] x) + /// + /// Solves A*X=B for X using a previously QR factored matrix. + /// + /// The Q matrix obtained by calling . + /// The R matrix obtained by calling . + /// 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. + public void QRSolveFactored(Complex[] q, Complex[] r, int rRows, int rColumns, Complex[] b, int bColumns, Complex[] x) { throw new NotImplementedException(); } - public void SingularValueDecomposition(bool computeVectors, Complex[] a, Complex[] s, Complex[] u, Complex[] 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. + /// This is equivalent to the GESVD LAPACK routine. + public void SingularValueDecomposition(bool computeVectors, Complex[] a, int aRows, int aColumns, Complex[] s, Complex[] u, Complex[] vt) { throw new NotImplementedException(); } - public void SingularValueDecomposition(bool computeVectors, Complex[] a, Complex[] s, Complex[] u, Complex[] vt, Complex[] 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. 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. + public void SingularValueDecomposition(bool computeVectors, Complex[] a, int aRows, int aColumns, Complex[] s, Complex[] u, Complex[] vt, Complex[] work) { throw new NotImplementedException(); } - public void SvdSolve(Complex[] a, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, Complex[] 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 number of rows in the A matrix. + /// The number of columns in the A matrix. + /// 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. + /// The number of columns of B. + /// On exit, the solution matrix. + public void SvdSolve(Complex[] a, int aRows, int aColumns, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, int bColumns, Complex[] x) { throw new NotImplementedException(); } - public void SvdSolve(Complex[] a, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, Complex[] x, Complex[] 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 number of rows in the A matrix. + /// The number of columns in the A matrix. + /// 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. + /// The number of columns of B. + /// 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. + public void SvdSolve(Complex[] a, int aRows, int aColumns, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, int bColumns, Complex[] x, Complex[] work) { throw new NotImplementedException(); } - public void SvdSolveFactored(int columnsOfB, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, Complex[] x) + /// + /// Solves A*X=B for X using a previously SVD decomposed matrix. + /// + /// The number of rows in the A matrix. + /// 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 B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + public void SvdSolveFactored(int aRows, int aColumns, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, int bColumns, Complex[] x) { throw new NotImplementedException(); } @@ -3558,52 +5671,179 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new NotImplementedException(); } - public void QRFactor(Complex32[] r, Complex32[] 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. + /// 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. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + public void QRFactor(Complex32[] r, int rRows, int rColumns, Complex32[] q) { throw new NotImplementedException(); } - public void QRFactor(Complex32[] r, Complex32[] q, Complex32[] 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. + /// 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 void QRFactor(Complex32[] r, int rRows, int rColumns, Complex32[] q, Complex32[] work) { throw new NotImplementedException(); } - public void QRSolve(int columnsOfB, Complex32[] r, Complex32[] q, Complex32[] b, Complex32[] x) + /// + /// Solves A*X=B for X using 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. + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + public void QRSolve(Complex32[] r, int rRows, int rColumns, Complex32[] q, Complex32[] b, int bColumns, Complex32[] x) { throw new NotImplementedException(); } - public void QRSolve(int columnsOfB, Complex32[] r, Complex32[] q, Complex32[] b, Complex32[] x, Complex32[] work) + /// + /// Solves A*X=B for X using 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. + /// 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. + public void QRSolve(Complex32[] r, int rRows, int rColumns, Complex32[] q, Complex32[] b, int bColumns, Complex32[] x, Complex32[] work) { throw new NotImplementedException(); } - public void QRSolveFactored(int columnsOfB, Complex32[] q, Complex32[] r, Complex32[] b, Complex32[] x) + /// + /// Solves A*X=B for X using a previously QR factored matrix. + /// + /// The Q matrix obtained by calling . + /// The R matrix obtained by calling . + /// 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. + public void QRSolveFactored(Complex32[] q, Complex32[] r, int rRows, int rColumns, Complex32[] b, int bColumns, Complex32[] x) { throw new NotImplementedException(); } - public void SingularValueDecomposition(bool computeVectors, Complex32[] a, Complex32[] s, Complex32[] u, Complex32[] 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. + /// This is equivalent to the GESVD LAPACK routine. + public void SingularValueDecomposition(bool computeVectors, Complex32[] a, int aRows, int aColumns, Complex32[] s, Complex32[] u, Complex32[] vt) { throw new NotImplementedException(); } - public void SingularValueDecomposition(bool computeVectors, Complex32[] a, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] 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. 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. + public void SingularValueDecomposition(bool computeVectors, Complex32[] a, int aRows, int aColumns, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] work) { throw new NotImplementedException(); } - public void SvdSolve(Complex32[] a, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, Complex32[] 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 number of rows in the A matrix. + /// The number of columns in the A matrix. + /// 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. + /// The number of columns of B. + /// On exit, the solution matrix. + public void SvdSolve(Complex32[] a, int aRows, int aColumns, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, int bColumns, Complex32[] x) { throw new NotImplementedException(); } - public void SvdSolve(Complex32[] a, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, Complex32[] x, Complex32[] 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 number of rows in the A matrix. + /// The number of columns in the A matrix. + /// 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. + /// The number of columns of B. + /// 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. + public void SvdSolve(Complex32[] a, int aRows, int aColumns, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, int bColumns, Complex32[] x, Complex32[] work) { throw new NotImplementedException(); } - public void SvdSolveFactored(int columnsOfB, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, Complex32[] x) + /// + /// Solves A*X=B for X using a previously SVD decomposed matrix. + /// + /// The number of rows in the A matrix. + /// 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 B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + public void SvdSolveFactored(int aRows, int aColumns, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, int bColumns, Complex32[] x) { throw new NotImplementedException(); } diff --git a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs index 0281e06b..43d26741 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/Mkl/MklLinearAlgebraProvider.cs @@ -28,7 +28,7 @@ /* This file is automatically generated - do not modify it. Change NativeLinearAlgebraProvider.include instead. - Last generated on UTC 2010-06-27 13:41:29Z + Last generated on UTC 2010-07-04 13:21:48Z */ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl @@ -234,6 +234,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// Computes the requested of the matrix. /// /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. /// The matrix to compute the norm from. /// /// The requested of the matrix. @@ -247,6 +249,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// Computes the requested of the matrix. /// /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. /// 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 . @@ -451,24 +455,25 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl { throw new ArgumentNullException("a"); } - - if ( order < 1) + + if (order < 1) { - throw new ArgumentException(Properties.Resources.ArgumentMustBePositive, "order"); + throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } - - SafeNativeMethods.d_cholesky_factor(order, a); + + SafeNativeMethods.d_cholesky_factor(order, a); } - + /// /// Solves A*X=B for X using Cholesky factorization. /// /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. - /// The number of columns in the B matrix. - /// This is equivalent to the POTRF add POTRS LAPACK routines. + /// The number of rows in the B matrix. + /// The number of columns in the B matrix. + /// This is equivalent to the POTRF add POTRS LAPACK routines. + /// public void CholeskySolve(double[] a, int aOrder, double[] b, int bRows, int bColumns) { throw new NotImplementedException(); @@ -480,8 +485,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. - /// The number of columns in the B matrix. + /// The number of rows in the B matrix. + /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. public void CholeskySolveFactored(double[] a, int aOrder, double[] b, int bRows, int bColumns) { @@ -492,11 +497,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// 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. + /// 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. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public void QRFactor(double[] r, double[] q) + public void QRFactor(double[] r, int rRows, int rColumns, double[] q) { throw new NotImplementedException(); } @@ -505,13 +512,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// 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 + /// 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. /// 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. - public void QRFactor(double[] r, double[] q, double[] work) + /// This is similar to the GEQRF and ORGQR LAPACK routines. + public void QRFactor(double[] r, int rRows, int rColumns, double[] q, double[] work) { throw new NotImplementedException(); } @@ -519,14 +529,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// 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 + /// 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. /// The B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void QRSolve(int columnsOfB, double[] r, double[] q, double[] b, double[] x) + public void QRSolve(double[] r, int rRows, int rColumns, double[] q, double[] b, int bColumns, double[] x) { throw new NotImplementedException(); } @@ -534,17 +546,19 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// 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 + /// 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. /// 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. - public void QRSolve(int columnsOfB, double[] r, double[] q, double[] b, double[] x, double[] work) + public void QRSolve(double[] r, int rRows, int rColumns, double[] q, double[] b, int bColumns, double[] x, double[] work) { throw new NotImplementedException(); } @@ -552,12 +566,14 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// 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 Q matrix obtained by calling . + /// The R matrix obtained by calling . + /// 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. - public void QRSolveFactored(int columnsOfB, double[] q, double[] r, double[] b, double[] x) + public void QRSolveFactored(double[] q, double[] r, int rRows, int rColumns, double[] b, int bColumns, double[] x) { throw new NotImplementedException(); } @@ -567,13 +583,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// 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. - public void SingularValueDecomposition(bool computeVectors, double[] a, double[] s, double[] u, double[] vt) + public void SingularValueDecomposition(bool computeVectors, double[] a, int aRows, int aColumns, double[] s, double[] u, double[] vt) { throw new NotImplementedException(); } @@ -583,6 +601,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// 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. @@ -592,7 +612,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// 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. - public void SingularValueDecomposition(bool computeVectors, double[] a, double[] s, double[] u, double[] vt, double[] work) + public void SingularValueDecomposition(bool computeVectors, double[] a, int aRows, int aColumns, double[] s, double[] u, double[] vt, double[] work) { throw new NotImplementedException(); } @@ -601,12 +621,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// On exit, the solution matrix. - public void SvdSolve(double[] a, double[] s, double[] u, double[] vt, double[] b, double[] x) + public void SvdSolve(double[] a, int aRows, int aColumns, double[] s, double[] u, double[] vt, double[] b, int bColumns, double[] x) { throw new NotImplementedException(); } @@ -615,15 +638,18 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// 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. - public void SvdSolve(double[] a, double[] s, double[] u, double[] vt, double[] b, double[] x, double[] work) + public void SvdSolve(double[] a, int aRows, int aColumns, double[] s, double[] u, double[] vt, double[] b, int bColumns, double[] x, double[] work) { throw new NotImplementedException(); } @@ -631,13 +657,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// 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 number of rows in the A matrix. + /// 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 B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void SvdSolveFactored(int columnsOfB, double[] s, double[] u, double[] vt, double[] b, double[] x) + public void SvdSolveFactored(int aRows, int aColumns, double[] s, double[] u, double[] vt, double[] b, int bColumns, double[] x) { throw new NotImplementedException(); } @@ -835,6 +863,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// Computes the requested of the matrix. /// /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. /// The matrix to compute the norm from. /// /// The requested of the matrix. @@ -848,6 +878,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// Computes the requested of the matrix. /// /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. /// 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 . @@ -925,6 +957,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl SafeNativeMethods.s_matrix_multiply(transposeA, transposeB, m, n, k, alpha, a, b, beta, c); } + /// /// /// Computes the LUP factorization of A. P*A = L*U. /// @@ -1052,23 +1085,23 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl { throw new ArgumentNullException("a"); } - - if ( order < 1) + + if (order < 1) { - throw new ArgumentException(Properties.Resources.ArgumentMustBePositive, "order"); + throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } - - SafeNativeMethods.s_cholesky_factor(order, a); + + SafeNativeMethods.s_cholesky_factor(order, a); } - + /// /// Solves A*X=B for X using Cholesky factorization. /// /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. - /// The number of columns in the B matrix. + /// The number of rows in the B matrix. + /// The number of columns in the B matrix. /// This is equivalent to the POTRF add POTRS LAPACK routines. public void CholeskySolve(float[] a, int aOrder, float[] b, int bRows, int bColumns) { @@ -1081,8 +1114,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. - /// The number of columns in the B matrix. + /// The number of rows in the B matrix. + /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. public void CholeskySolveFactored(float[] a, int aOrder, float[] b, int bRows, int bColumns) { @@ -1093,11 +1126,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// 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. + /// 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. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public void QRFactor(float[] r, float[] q) + public void QRFactor(float[] r, int rRows, int rColumns, float[] q) { throw new NotImplementedException(); } @@ -1106,13 +1141,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// 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 + /// 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. /// 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. - public void QRFactor(float[] r, float[] q, float[] work) + /// This is similar to the GEQRF and ORGQR LAPACK routines. + public void QRFactor(float[] r, int rRows, int rColumns, float[] q, float[] work) { throw new NotImplementedException(); } @@ -1120,14 +1158,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// 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 + /// 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. /// The B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void QRSolve(int columnsOfB, float[] r, float[] q, float[] b, float[] x) + public void QRSolve(float[] r, int rRows, int rColumns, float[] q, float[] b, int bColumns, float[] x) { throw new NotImplementedException(); } @@ -1135,17 +1175,19 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// 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 + /// 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. /// 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. - public void QRSolve(int columnsOfB, float[] r, float[] q, float[] b, float[] x, float[] work) + public void QRSolve(float[] r, int rRows, int rColumns, float[] q, float[] b, int bColumns, float[] x, float[] work) { throw new NotImplementedException(); } @@ -1153,12 +1195,14 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// 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 Q matrix obtained by calling . + /// The R matrix obtained by calling . + /// 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. - public void QRSolveFactored(int columnsOfB, float[] q, float[] r, float[] b, float[] x) + public void QRSolveFactored(float[] q, float[] r, int rRows, int rColumns, float[] b, int bColumns, float[] x) { throw new NotImplementedException(); } @@ -1168,13 +1212,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// 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. - public void SingularValueDecomposition(bool computeVectors, float[] a, float[] s, float[] u, float[] vt) + public void SingularValueDecomposition(bool computeVectors, float[] a, int aRows, int aColumns, float[] s, float[] u, float[] vt) { throw new NotImplementedException(); } @@ -1184,6 +1230,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// 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. @@ -1193,7 +1241,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// 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. - public void SingularValueDecomposition(bool computeVectors, float[] a, float[] s, float[] u, float[] vt, float[] work) + public void SingularValueDecomposition(bool computeVectors, float[] a, int aRows, int aColumns, float[] s, float[] u, float[] vt, float[] work) { throw new NotImplementedException(); } @@ -1202,12 +1250,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// On exit, the solution matrix. - public void SvdSolve(float[] a, float[] s, float[] u, float[] vt, float[] b, float[] x) + public void SvdSolve(float[] a, int aRows, int aColumns, float[] s, float[] u, float[] vt, float[] b, int bColumns, float[] x) { throw new NotImplementedException(); } @@ -1216,15 +1267,18 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// 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. - public void SvdSolve(float[] a, float[] s, float[] u, float[] vt, float[] b, float[] x, float[] work) + public void SvdSolve(float[] a, int aRows, int aColumns, float[] s, float[] u, float[] vt, float[] b, int bColumns, float[] x, float[] work) { throw new NotImplementedException(); } @@ -1232,13 +1286,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// 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 number of rows in the A matrix. + /// 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 B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void SvdSolveFactored(int columnsOfB, float[] s, float[] u, float[] vt, float[] b, float[] x) + public void SvdSolveFactored(int aRows, int aColumns, float[] s, float[] u, float[] vt, float[] b, int bColumns, float[] x) { throw new NotImplementedException(); } @@ -1436,6 +1492,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// Computes the requested of the matrix. /// /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. /// The matrix to compute the norm from. /// /// The requested of the matrix. @@ -1449,6 +1507,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// Computes the requested of the matrix. /// /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. /// 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 . @@ -1653,23 +1713,23 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl { throw new ArgumentNullException("a"); } - - if ( order < 1) + + if (order < 1) { - throw new ArgumentException(Properties.Resources.ArgumentMustBePositive, "order"); + throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } - - SafeNativeMethods.z_cholesky_factor(order, a); + + SafeNativeMethods.z_cholesky_factor(order, a); } - + /// /// Solves A*X=B for X using Cholesky factorization. /// /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. - /// The number of columns in the B matrix. + /// The number of rows in the B matrix. + /// The number of columns in the B matrix. /// This is equivalent to the POTRF add POTRS LAPACK routines. public void CholeskySolve(Complex[] a, int aOrder, Complex[] b, int bRows, int bColumns) { @@ -1682,8 +1742,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. - /// The number of columns in the B matrix. + /// The number of rows in the B matrix. + /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. public void CholeskySolveFactored(Complex[] a, int aOrder, Complex[] b, int bRows, int bColumns) { @@ -1694,11 +1754,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// 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. + /// 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. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public void QRFactor(Complex[] r, Complex[] q) + public void QRFactor(Complex[] r, int rRows, int rColumns, Complex[] q) { throw new NotImplementedException(); } @@ -1707,13 +1769,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// 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 + /// 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. /// 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. - public void QRFactor(Complex[] r, Complex[] q, Complex[] work) + /// This is similar to the GEQRF and ORGQR LAPACK routines. + public void QRFactor(Complex[] r, int rRows, int rColumns, Complex[] q, Complex[] work) { throw new NotImplementedException(); } @@ -1721,14 +1786,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// 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 + /// 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. /// The B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void QRSolve(int columnsOfB, Complex[] r, Complex[] q, Complex[] b, Complex[] x) + public void QRSolve(Complex[] r, int rRows, int rColumns, Complex[] q, Complex[] b, int bColumns, Complex[] x) { throw new NotImplementedException(); } @@ -1736,17 +1803,19 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// 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 + /// 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. /// 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. - public void QRSolve(int columnsOfB, Complex[] r, Complex[] q, Complex[] b, Complex[] x, Complex[] work) + public void QRSolve(Complex[] r, int rRows, int rColumns, Complex[] q, Complex[] b, int bColumns, Complex[] x, Complex[] work) { throw new NotImplementedException(); } @@ -1754,12 +1823,14 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// 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 Q matrix obtained by calling . + /// The R matrix obtained by calling . + /// 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. - public void QRSolveFactored(int columnsOfB, Complex[] q, Complex[] r, Complex[] b, Complex[] x) + public void QRSolveFactored(Complex[] q, Complex[] r, int rRows, int rColumns, Complex[] b, int bColumns, Complex[] x) { throw new NotImplementedException(); } @@ -1769,13 +1840,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// 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. - public void SingularValueDecomposition(bool computeVectors, Complex[] a, Complex[] s, Complex[] u, Complex[] vt) + public void SingularValueDecomposition(bool computeVectors, Complex[] a, int aRows, int aColumns, Complex[] s, Complex[] u, Complex[] vt) { throw new NotImplementedException(); } @@ -1785,6 +1858,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// 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. @@ -1794,7 +1869,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// 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. - public void SingularValueDecomposition(bool computeVectors, Complex[] a, Complex[] s, Complex[] u, Complex[] vt, Complex[] work) + public void SingularValueDecomposition(bool computeVectors, Complex[] a, int aRows, int aColumns, Complex[] s, Complex[] u, Complex[] vt, Complex[] work) { throw new NotImplementedException(); } @@ -1803,12 +1878,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// On exit, the solution matrix. - public void SvdSolve(Complex[] a, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, Complex[] x) + public void SvdSolve(Complex[] a, int aRows, int aColumns, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, int bColumns, Complex[] x) { throw new NotImplementedException(); } @@ -1817,15 +1895,18 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// 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. - public void SvdSolve(Complex[] a, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, Complex[] x, Complex[] work) + public void SvdSolve(Complex[] a, int aRows, int aColumns, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, int bColumns, Complex[] x, Complex[] work) { throw new NotImplementedException(); } @@ -1833,13 +1914,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// 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 number of rows in the A matrix. + /// 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 B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void SvdSolveFactored(int columnsOfB, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, Complex[] x) + public void SvdSolveFactored(int aRows, int aColumns, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, int bColumns, Complex[] x) { throw new NotImplementedException(); } @@ -2258,23 +2341,23 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl { throw new ArgumentNullException("a"); } - - if ( order < 1) + + if (order < 1) { - throw new ArgumentException(Properties.Resources.ArgumentMustBePositive, "order"); + throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } - - SafeNativeMethods.c_cholesky_factor(order, a); + + SafeNativeMethods.c_cholesky_factor(order, a); } - + /// /// Solves A*X=B for X using Cholesky factorization. /// /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. - /// The number of columns in the B matrix. + /// The number of rows in the B matrix. + /// The number of columns in the B matrix. /// This is equivalent to the POTRF add POTRS LAPACK routines. public void CholeskySolve(Complex32[] a, int aOrder, Complex32[] b, int bRows, int bColumns) { @@ -2287,8 +2370,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. - /// The number of columns in the B matrix. + /// The number of rows in the B matrix. + /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. public void CholeskySolveFactored(Complex32[] a, int aOrder, Complex32[] b, int bRows, int bColumns) { @@ -2299,11 +2382,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// 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. + /// 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. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public void QRFactor(Complex32[] r, Complex32[] q) + public void QRFactor(Complex32[] r, int rRows, int rColumns, Complex32[] q) { throw new NotImplementedException(); } @@ -2312,13 +2397,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// 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 + /// 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. /// 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. - public void QRFactor(Complex32[] r, Complex32[] q, Complex32[] work) + /// This is similar to the GEQRF and ORGQR LAPACK routines. + public void QRFactor(Complex32[] r, int rRows, int rColumns, Complex32[] q, Complex32[] work) { throw new NotImplementedException(); } @@ -2326,14 +2414,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// 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 + /// 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. /// The B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void QRSolve(int columnsOfB, Complex32[] r, Complex32[] q, Complex32[] b, Complex32[] x) + public void QRSolve(Complex32[] r, int rRows, int rColumns, Complex32[] q, Complex32[] b, int bColumns, Complex32[] x) { throw new NotImplementedException(); } @@ -2341,17 +2431,19 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// 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 + /// 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. /// 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. - public void QRSolve(int columnsOfB, Complex32[] r, Complex32[] q, Complex32[] b, Complex32[] x, Complex32[] work) + public void QRSolve(Complex32[] r, int rRows, int rColumns, Complex32[] q, Complex32[] b, int bColumns, Complex32[] x, Complex32[] work) { throw new NotImplementedException(); } @@ -2359,12 +2451,14 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// 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 Q matrix obtained by calling . + /// The R matrix obtained by calling . + /// 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. - public void QRSolveFactored(int columnsOfB, Complex32[] q, Complex32[] r, Complex32[] b, Complex32[] x) + public void QRSolveFactored(Complex32[] q, Complex32[] r, int rRows, int rColumns, Complex32[] b, int bColumns, Complex32[] x) { throw new NotImplementedException(); } @@ -2374,13 +2468,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// 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. - public void SingularValueDecomposition(bool computeVectors, Complex32[] a, Complex32[] s, Complex32[] u, Complex32[] vt) + public void SingularValueDecomposition(bool computeVectors, Complex32[] a, int aRows, int aColumns, Complex32[] s, Complex32[] u, Complex32[] vt) { throw new NotImplementedException(); } @@ -2390,6 +2486,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// 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. @@ -2399,7 +2497,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// 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. - public void SingularValueDecomposition(bool computeVectors, Complex32[] a, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] work) + public void SingularValueDecomposition(bool computeVectors, Complex32[] a, int aRows, int aColumns, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] work) { throw new NotImplementedException(); } @@ -2408,12 +2506,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// On exit, the solution matrix. - public void SvdSolve(Complex32[] a, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, Complex32[] x) + public void SvdSolve(Complex32[] a, int aRows, int aColumns, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, int bColumns, Complex32[] x) { throw new NotImplementedException(); } @@ -2422,15 +2523,18 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// 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. - public void SvdSolve(Complex32[] a, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, Complex32[] x, Complex32[] work) + public void SvdSolve(Complex32[] a, int aRows, int aColumns, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, int bColumns, Complex32[] x, Complex32[] work) { throw new NotImplementedException(); } @@ -2438,13 +2542,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.Mkl /// /// 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 number of rows in the A matrix. + /// 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 B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void SvdSolveFactored(int columnsOfB, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, Complex32[] x) + public void SvdSolveFactored(int aRows, int aColumns, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, int bColumns, Complex32[] x) { throw new NotImplementedException(); } diff --git a/src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include b/src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include index 0c8bcc52..4c5f734c 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include +++ b/src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include @@ -249,6 +249,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// Computes the requested of the matrix. /// /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. /// The matrix to compute the norm from. /// /// The requested of the matrix. @@ -262,6 +264,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// Computes the requested of the matrix. /// /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. /// 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 . @@ -466,24 +470,25 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> { throw new ArgumentNullException("a"); } - - if ( order < 1) + + if (order < 1) { - throw new ArgumentException(Properties.Resources.ArgumentMustBePositive, "order"); + throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } - - SafeNativeMethods.d_cholesky_factor(order, a); + + SafeNativeMethods.d_cholesky_factor(order, a); } - + /// /// Solves A*X=B for X using Cholesky factorization. /// /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. - /// The number of columns in the B matrix. - /// This is equivalent to the POTRF add POTRS LAPACK routines. + /// The number of rows in the B matrix. + /// The number of columns in the B matrix. + /// This is equivalent to the POTRF add POTRS LAPACK routines. + /// public void CholeskySolve(double[] a, int aOrder, double[] b, int bRows, int bColumns) { throw new NotImplementedException(); @@ -495,8 +500,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. - /// The number of columns in the B matrix. + /// The number of rows in the B matrix. + /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. public void CholeskySolveFactored(double[] a, int aOrder, double[] b, int bRows, int bColumns) { @@ -507,11 +512,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// 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. + /// 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. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public void QRFactor(double[] r, double[] q) + public void QRFactor(double[] r, int rRows, int rColumns, double[] q) { throw new NotImplementedException(); } @@ -520,13 +527,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// 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 + /// 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. /// 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. - public void QRFactor(double[] r, double[] q, double[] work) + /// This is similar to the GEQRF and ORGQR LAPACK routines. + public void QRFactor(double[] r, int rRows, int rColumns, double[] q, double[] work) { throw new NotImplementedException(); } @@ -534,14 +544,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// 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 + /// 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. /// The B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void QRSolve(int columnsOfB, double[] r, double[] q, double[] b, double[] x) + public void QRSolve(double[] r, int rRows, int rColumns, double[] q, double[] b, int bColumns, double[] x) { throw new NotImplementedException(); } @@ -549,17 +561,19 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// 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 + /// 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. /// 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. - public void QRSolve(int columnsOfB, double[] r, double[] q, double[] b, double[] x, double[] work) + public void QRSolve(double[] r, int rRows, int rColumns, double[] q, double[] b, int bColumns, double[] x, double[] work) { throw new NotImplementedException(); } @@ -567,12 +581,14 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// 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 Q matrix obtained by calling . + /// The R matrix obtained by calling . + /// 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. - public void QRSolveFactored(int columnsOfB, double[] q, double[] r, double[] b, double[] x) + public void QRSolveFactored(double[] q, double[] r, int rRows, int rColumns, double[] b, int bColumns, double[] x) { throw new NotImplementedException(); } @@ -582,13 +598,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// 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. - public void SingularValueDecomposition(bool computeVectors, double[] a, double[] s, double[] u, double[] vt) + public void SingularValueDecomposition(bool computeVectors, double[] a, int aRows, int aColumns, double[] s, double[] u, double[] vt) { throw new NotImplementedException(); } @@ -598,6 +616,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// 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. @@ -607,7 +627,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// 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. - public void SingularValueDecomposition(bool computeVectors, double[] a, double[] s, double[] u, double[] vt, double[] work) + public void SingularValueDecomposition(bool computeVectors, double[] a, int aRows, int aColumns, double[] s, double[] u, double[] vt, double[] work) { throw new NotImplementedException(); } @@ -616,12 +636,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// On exit, the solution matrix. - public void SvdSolve(double[] a, double[] s, double[] u, double[] vt, double[] b, double[] x) + public void SvdSolve(double[] a, int aRows, int aColumns, double[] s, double[] u, double[] vt, double[] b, int bColumns, double[] x) { throw new NotImplementedException(); } @@ -630,15 +653,18 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// 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. - public void SvdSolve(double[] a, double[] s, double[] u, double[] vt, double[] b, double[] x, double[] work) + public void SvdSolve(double[] a, int aRows, int aColumns, double[] s, double[] u, double[] vt, double[] b, int bColumns, double[] x, double[] work) { throw new NotImplementedException(); } @@ -646,13 +672,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// 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 number of rows in the A matrix. + /// 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 B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void SvdSolveFactored(int columnsOfB, double[] s, double[] u, double[] vt, double[] b, double[] x) + public void SvdSolveFactored(int aRows, int aColumns, double[] s, double[] u, double[] vt, double[] b, int bColumns, double[] x) { throw new NotImplementedException(); } @@ -862,6 +890,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// Computes the requested of the matrix. /// /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. /// The matrix to compute the norm from. /// /// The requested of the matrix. @@ -875,6 +905,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// Computes the requested of the matrix. /// /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. /// 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 . @@ -952,6 +984,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> SafeNativeMethods.s_matrix_multiply(transposeA, transposeB, m, n, k, alpha, a, b, beta, c); } + /// /// /// Computes the LUP factorization of A. P*A = L*U. /// @@ -1079,23 +1112,23 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> { throw new ArgumentNullException("a"); } - - if ( order < 1) + + if (order < 1) { - throw new ArgumentException(Properties.Resources.ArgumentMustBePositive, "order"); + throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } - - SafeNativeMethods.s_cholesky_factor(order, a); + + SafeNativeMethods.s_cholesky_factor(order, a); } - + /// /// Solves A*X=B for X using Cholesky factorization. /// /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. - /// The number of columns in the B matrix. + /// The number of rows in the B matrix. + /// The number of columns in the B matrix. /// This is equivalent to the POTRF add POTRS LAPACK routines. public void CholeskySolve(float[] a, int aOrder, float[] b, int bRows, int bColumns) { @@ -1108,8 +1141,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. - /// The number of columns in the B matrix. + /// The number of rows in the B matrix. + /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. public void CholeskySolveFactored(float[] a, int aOrder, float[] b, int bRows, int bColumns) { @@ -1120,11 +1153,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// 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. + /// 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. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public void QRFactor(float[] r, float[] q) + public void QRFactor(float[] r, int rRows, int rColumns, float[] q) { throw new NotImplementedException(); } @@ -1133,13 +1168,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// 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 + /// 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. /// 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. - public void QRFactor(float[] r, float[] q, float[] work) + /// This is similar to the GEQRF and ORGQR LAPACK routines. + public void QRFactor(float[] r, int rRows, int rColumns, float[] q, float[] work) { throw new NotImplementedException(); } @@ -1147,14 +1185,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// 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 + /// 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. /// The B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void QRSolve(int columnsOfB, float[] r, float[] q, float[] b, float[] x) + public void QRSolve(float[] r, int rRows, int rColumns, float[] q, float[] b, int bColumns, float[] x) { throw new NotImplementedException(); } @@ -1162,17 +1202,19 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// 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 + /// 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. /// 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. - public void QRSolve(int columnsOfB, float[] r, float[] q, float[] b, float[] x, float[] work) + public void QRSolve(float[] r, int rRows, int rColumns, float[] q, float[] b, int bColumns, float[] x, float[] work) { throw new NotImplementedException(); } @@ -1180,12 +1222,14 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// 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 Q matrix obtained by calling . + /// The R matrix obtained by calling . + /// 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. - public void QRSolveFactored(int columnsOfB, float[] q, float[] r, float[] b, float[] x) + public void QRSolveFactored(float[] q, float[] r, int rRows, int rColumns, float[] b, int bColumns, float[] x) { throw new NotImplementedException(); } @@ -1195,13 +1239,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// 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. - public void SingularValueDecomposition(bool computeVectors, float[] a, float[] s, float[] u, float[] vt) + public void SingularValueDecomposition(bool computeVectors, float[] a, int aRows, int aColumns, float[] s, float[] u, float[] vt) { throw new NotImplementedException(); } @@ -1211,6 +1257,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// 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. @@ -1220,7 +1268,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// 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. - public void SingularValueDecomposition(bool computeVectors, float[] a, float[] s, float[] u, float[] vt, float[] work) + public void SingularValueDecomposition(bool computeVectors, float[] a, int aRows, int aColumns, float[] s, float[] u, float[] vt, float[] work) { throw new NotImplementedException(); } @@ -1229,12 +1277,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// On exit, the solution matrix. - public void SvdSolve(float[] a, float[] s, float[] u, float[] vt, float[] b, float[] x) + public void SvdSolve(float[] a, int aRows, int aColumns, float[] s, float[] u, float[] vt, float[] b, int bColumns, float[] x) { throw new NotImplementedException(); } @@ -1243,15 +1294,18 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// 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. - public void SvdSolve(float[] a, float[] s, float[] u, float[] vt, float[] b, float[] x, float[] work) + public void SvdSolve(float[] a, int aRows, int aColumns, float[] s, float[] u, float[] vt, float[] b, int bColumns, float[] x, float[] work) { throw new NotImplementedException(); } @@ -1259,13 +1313,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// 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 number of rows in the A matrix. + /// 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 B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void SvdSolveFactored(int columnsOfB, float[] s, float[] u, float[] vt, float[] b, float[] x) + public void SvdSolveFactored(int aRows, int aColumns, float[] s, float[] u, float[] vt, float[] b, int bColumns, float[] x) { throw new NotImplementedException(); } @@ -1475,6 +1531,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// Computes the requested of the matrix. /// /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. /// The matrix to compute the norm from. /// /// The requested of the matrix. @@ -1488,6 +1546,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// Computes the requested of the matrix. /// /// The type of norm to compute. + /// The number of rows in the matrix. + /// The number of columns in the matrix. /// 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 . @@ -1692,23 +1752,23 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> { throw new ArgumentNullException("a"); } - - if ( order < 1) + + if (order < 1) { - throw new ArgumentException(Properties.Resources.ArgumentMustBePositive, "order"); + throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } - - SafeNativeMethods.z_cholesky_factor(order, a); + + SafeNativeMethods.z_cholesky_factor(order, a); } - + /// /// Solves A*X=B for X using Cholesky factorization. /// /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. - /// The number of columns in the B matrix. + /// The number of rows in the B matrix. + /// The number of columns in the B matrix. /// This is equivalent to the POTRF add POTRS LAPACK routines. public void CholeskySolve(Complex[] a, int aOrder, Complex[] b, int bRows, int bColumns) { @@ -1721,8 +1781,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. - /// The number of columns in the B matrix. + /// The number of rows in the B matrix. + /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. public void CholeskySolveFactored(Complex[] a, int aOrder, Complex[] b, int bRows, int bColumns) { @@ -1733,11 +1793,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// 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. + /// 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. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public void QRFactor(Complex[] r, Complex[] q) + public void QRFactor(Complex[] r, int rRows, int rColumns, Complex[] q) { throw new NotImplementedException(); } @@ -1746,13 +1808,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// 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 + /// 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. /// 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. - public void QRFactor(Complex[] r, Complex[] q, Complex[] work) + /// This is similar to the GEQRF and ORGQR LAPACK routines. + public void QRFactor(Complex[] r, int rRows, int rColumns, Complex[] q, Complex[] work) { throw new NotImplementedException(); } @@ -1760,14 +1825,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// 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 + /// 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. /// The B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void QRSolve(int columnsOfB, Complex[] r, Complex[] q, Complex[] b, Complex[] x) + public void QRSolve(Complex[] r, int rRows, int rColumns, Complex[] q, Complex[] b, int bColumns, Complex[] x) { throw new NotImplementedException(); } @@ -1775,17 +1842,19 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// 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 + /// 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. /// 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. - public void QRSolve(int columnsOfB, Complex[] r, Complex[] q, Complex[] b, Complex[] x, Complex[] work) + public void QRSolve(Complex[] r, int rRows, int rColumns, Complex[] q, Complex[] b, int bColumns, Complex[] x, Complex[] work) { throw new NotImplementedException(); } @@ -1793,12 +1862,14 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// 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 Q matrix obtained by calling . + /// The R matrix obtained by calling . + /// 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. - public void QRSolveFactored(int columnsOfB, Complex[] q, Complex[] r, Complex[] b, Complex[] x) + public void QRSolveFactored(Complex[] q, Complex[] r, int rRows, int rColumns, Complex[] b, int bColumns, Complex[] x) { throw new NotImplementedException(); } @@ -1808,13 +1879,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// 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. - public void SingularValueDecomposition(bool computeVectors, Complex[] a, Complex[] s, Complex[] u, Complex[] vt) + public void SingularValueDecomposition(bool computeVectors, Complex[] a, int aRows, int aColumns, Complex[] s, Complex[] u, Complex[] vt) { throw new NotImplementedException(); } @@ -1824,6 +1897,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// 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. @@ -1833,7 +1908,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// 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. - public void SingularValueDecomposition(bool computeVectors, Complex[] a, Complex[] s, Complex[] u, Complex[] vt, Complex[] work) + public void SingularValueDecomposition(bool computeVectors, Complex[] a, int aRows, int aColumns, Complex[] s, Complex[] u, Complex[] vt, Complex[] work) { throw new NotImplementedException(); } @@ -1842,12 +1917,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// On exit, the solution matrix. - public void SvdSolve(Complex[] a, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, Complex[] x) + public void SvdSolve(Complex[] a, int aRows, int aColumns, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, int bColumns, Complex[] x) { throw new NotImplementedException(); } @@ -1856,15 +1934,18 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// 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. - public void SvdSolve(Complex[] a, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, Complex[] x, Complex[] work) + public void SvdSolve(Complex[] a, int aRows, int aColumns, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, int bColumns, Complex[] x, Complex[] work) { throw new NotImplementedException(); } @@ -1872,13 +1953,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// 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 number of rows in the A matrix. + /// 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 B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void SvdSolveFactored(int columnsOfB, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, Complex[] x) + public void SvdSolveFactored(int aRows, int aColumns, Complex[] s, Complex[] u, Complex[] vt, Complex[] b, int bColumns, Complex[] x) { throw new NotImplementedException(); } @@ -2309,23 +2392,23 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> { throw new ArgumentNullException("a"); } - - if ( order < 1) + + if (order < 1) { - throw new ArgumentException(Properties.Resources.ArgumentMustBePositive, "order"); + throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); } - - SafeNativeMethods.c_cholesky_factor(order, a); + + SafeNativeMethods.c_cholesky_factor(order, a); } - + /// /// Solves A*X=B for X using Cholesky factorization. /// /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. - /// The number of columns in the B matrix. + /// The number of rows in the B matrix. + /// The number of columns in the B matrix. /// This is equivalent to the POTRF add POTRS LAPACK routines. public void CholeskySolve(Complex32[] a, int aOrder, Complex32[] b, int bRows, int bColumns) { @@ -2338,8 +2421,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// The square, positive definite matrix A. /// The number of rows and columns in A. /// The B matrix. - /// The number of rows in the B matrix. - /// The number of columns in the B matrix. + /// The number of rows in the B matrix. + /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. public void CholeskySolveFactored(Complex32[] a, int aOrder, Complex32[] b, int bRows, int bColumns) { @@ -2350,11 +2433,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// 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. + /// 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. /// This is similar to the GEQRF and ORGQR LAPACK routines. - public void QRFactor(Complex32[] r, Complex32[] q) + public void QRFactor(Complex32[] r, int rRows, int rColumns, Complex32[] q) { throw new NotImplementedException(); } @@ -2363,13 +2448,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// 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 + /// 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. /// 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. - public void QRFactor(Complex32[] r, Complex32[] q, Complex32[] work) + /// This is similar to the GEQRF and ORGQR LAPACK routines. + public void QRFactor(Complex32[] r, int rRows, int rColumns, Complex32[] q, Complex32[] work) { throw new NotImplementedException(); } @@ -2377,14 +2465,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// 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 + /// 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. /// The B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void QRSolve(int columnsOfB, Complex32[] r, Complex32[] q, Complex32[] b, Complex32[] x) + public void QRSolve(Complex32[] r, int rRows, int rColumns, Complex32[] q, Complex32[] b, int bColumns, Complex32[] x) { throw new NotImplementedException(); } @@ -2392,17 +2482,19 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// 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 + /// 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. /// 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. - public void QRSolve(int columnsOfB, Complex32[] r, Complex32[] q, Complex32[] b, Complex32[] x, Complex32[] work) + public void QRSolve(Complex32[] r, int rRows, int rColumns, Complex32[] q, Complex32[] b, int bColumns, Complex32[] x, Complex32[] work) { throw new NotImplementedException(); } @@ -2410,12 +2502,14 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// 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 Q matrix obtained by calling . + /// The R matrix obtained by calling . + /// 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. - public void QRSolveFactored(int columnsOfB, Complex32[] q, Complex32[] r, Complex32[] b, Complex32[] x) + public void QRSolveFactored(Complex32[] q, Complex32[] r, int rRows, int rColumns, Complex32[] b, int bColumns, Complex32[] x) { throw new NotImplementedException(); } @@ -2425,13 +2519,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// 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. - public void SingularValueDecomposition(bool computeVectors, Complex32[] a, Complex32[] s, Complex32[] u, Complex32[] vt) + public void SingularValueDecomposition(bool computeVectors, Complex32[] a, int aRows, int aColumns, Complex32[] s, Complex32[] u, Complex32[] vt) { throw new NotImplementedException(); } @@ -2441,6 +2537,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// 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. @@ -2450,7 +2548,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// 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. - public void SingularValueDecomposition(bool computeVectors, Complex32[] a, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] work) + public void SingularValueDecomposition(bool computeVectors, Complex32[] a, int aRows, int aColumns, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] work) { throw new NotImplementedException(); } @@ -2459,12 +2557,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// On exit, the solution matrix. - public void SvdSolve(Complex32[] a, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, Complex32[] x) + public void SvdSolve(Complex32[] a, int aRows, int aColumns, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, int bColumns, Complex32[] x) { throw new NotImplementedException(); } @@ -2473,15 +2574,18 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// 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 number of rows in the A matrix. + /// The number of columns in the A matrix. /// 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. + /// The number of columns of B. /// 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. - public void SvdSolve(Complex32[] a, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, Complex32[] x, Complex32[] work) + public void SvdSolve(Complex32[] a, int aRows, int aColumns, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, int bColumns, Complex32[] x, Complex32[] work) { throw new NotImplementedException(); } @@ -2489,13 +2593,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> /// /// 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 number of rows in the A matrix. + /// 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 B matrix. + /// The number of columns of B. /// On exit, the solution matrix. - public void SvdSolveFactored(int columnsOfB, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, Complex32[] x) + public void SvdSolveFactored(int aRows, int aColumns, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] b, int bColumns, Complex32[] x) { throw new NotImplementedException(); } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs index 4fd3370e..e41f8755 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs @@ -32,7 +32,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { using System; using Properties; - using Threading; /// /// A class which encapsulates the functionality of the QR decomposition. @@ -65,7 +64,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization MatrixR = matrix.Clone(); MatrixQ = new DenseMatrix(matrix.RowCount); - Control.LinearAlgebraProvider.QRFactor(((DenseMatrix)MatrixR).Data, ((DenseMatrix)MatrixQ).Data); + Control.LinearAlgebraProvider.QRFactor(((DenseMatrix)MatrixR).Data, matrix.RowCount, matrix.ColumnCount, ((DenseMatrix)MatrixQ).Data); } /// @@ -116,18 +115,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw new NotImplementedException("Can only do QR factorization for dense matrices at the moment."); } - var solution = new double[dinput.Data.Length]; - Control.LinearAlgebraProvider.QRSolveFactored(input.ColumnCount, ((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, dinput.Data, solution); - CommonParallel.For( - 0, - dresult.RowCount, - row => - { - for (var col = 0; col < dresult.ColumnCount; col++) - { - dresult[row, col] = solution[row + (col * dinput.RowCount)]; - } - }); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, dinput.Data, input.ColumnCount, dresult.Data); } /// @@ -172,12 +160,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw new NotImplementedException("Can only do QR factorization for dense vectors at the moment."); } - var solution = new double[dinput.Data.Length]; - Control.LinearAlgebraProvider.QRSolveFactored(1, ((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, dinput.Data, solution); - CommonParallel.For( - 0, - dresult.Count, - index => { dresult[index] = solution[index]; }); + Control.LinearAlgebraProvider.QRSolveFactored(((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, MatrixR.RowCount, MatrixR.ColumnCount, dinput.Data, 1, dresult.Data); } } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseSvd.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseSvd.cs new file mode 100644 index 00000000..db81cfac --- /dev/null +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseSvd.cs @@ -0,0 +1,180 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2010 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. +// +namespace MathNet.Numerics.LinearAlgebra.Double.Factorization +{ + using System; + using Properties; + + /// + /// A class which encapsulates the functionality of the singular value decomposition (SVD) for . + /// Suppose M is an m-by-n matrix whose entries are real numbers. + /// Then there exists a factorization of the form M = UΣVT where: + /// - U is an m-by-m unitary matrix; + /// - Σ is m-by-n diagonal matrix with nonnegative real numbers on the diagonal; + /// - VT denotes transpose of V, an n-by-n unitary matrix; + /// Such a factorization is called a singular-value decomposition of M. A common convention is to order the diagonal + /// entries Σ(i,i) in descending order. In this case, the diagonal matrix Σ is uniquely determined + /// by M (though the matrices U and V are not). The diagonal entries of Σ are known as the singular values of M. + /// + /// + /// The computation of the singular value decomposition is done at construction time. + /// + public class DenseSvd : Svd + { + /// + /// Initializes a new instance of the class. This object will compute the + /// the singular value decomposition when the constructor is called and cache it's decomposition. + /// + /// The matrix to factor. + /// Compute the singular U and VT vectors or not. + /// If is null. + /// If SVD algorithm failed to converge with matrix . + public DenseSvd(DenseMatrix matrix, bool computeVectors) + { + if (matrix == null) + { + throw new ArgumentNullException("matrix"); + } + + ComputeVectors = computeVectors; + var nm = Math.Min(matrix.RowCount, matrix.ColumnCount); + VectorS = new DenseVector(nm); + MatrixU = new DenseMatrix(matrix.RowCount); + MatrixVT = new DenseMatrix(matrix.ColumnCount); + Control.LinearAlgebraProvider.SingularValueDecomposition(computeVectors, ((DenseMatrix)matrix.Clone()).Data, matrix.RowCount, matrix.ColumnCount, ((DenseVector)VectorS).Data, ((DenseMatrix)MatrixU).Data, ((DenseMatrix)MatrixVT).Data); + } + + /// + /// Solves a system of linear equations, AX = B, with A SVD factorized. + /// + /// The right hand side , B. + /// The left hand side , X. + public override void Solve(Matrix input, Matrix result) + { + // Check for proper arguments. + if (input == null) + { + throw new ArgumentNullException("input"); + } + + if (result == null) + { + throw new ArgumentNullException("result"); + } + + if (!ComputeVectors) + { + throw new InvalidOperationException(Resources.SingularVectorsNotComputed); + } + + // The solution X should have the same number of columns as B + if (input.ColumnCount != result.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); + } + + // The dimension compatibility conditions for X = A\B require the two matrices A and B to have the same number of rows + if (MatrixU.RowCount != input.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); + } + + // The solution X row dimension is equal to the column dimension of A + if (MatrixVT.ColumnCount != result.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); + } + + var dinput = input as DenseMatrix; + if (dinput == null) + { + throw new NotImplementedException("Can only do SVD factorization for dense matrices at the moment."); + } + + var dresult = result as DenseMatrix; + if (dresult == null) + { + throw new NotImplementedException("Can only do SVD factorization for dense matrices at the moment."); + } + + Control.LinearAlgebraProvider.SvdSolveFactored(MatrixU.RowCount, MatrixVT.ColumnCount, ((DenseVector)VectorS).Data, ((DenseMatrix)MatrixU).Data, ((DenseMatrix)MatrixVT).Data, dinput.Data, input.ColumnCount, dresult.Data); + } + + /// + /// Solves a system of linear equations, Ax = b, with A SVD factorized. + /// + /// The right hand side vector, b. + /// The left hand side , x. + public override void Solve(Vector input, Vector result) + { + if (input == null) + { + throw new ArgumentNullException("input"); + } + + if (result == null) + { + throw new ArgumentNullException("result"); + } + + if (!ComputeVectors) + { + throw new InvalidOperationException(Resources.SingularVectorsNotComputed); + } + + // Ax=b where A is an m x n matrix + // Check that b is a column vector with m entries + if (MatrixU.RowCount != input.Count) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + + // Check that x is a column vector with n entries + if (MatrixVT.ColumnCount != result.Count) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + var dinput = input as DenseVector; + if (dinput == null) + { + throw new NotImplementedException("Can only do QR factorization for dense vectors at the moment."); + } + + var dresult = result as DenseVector; + if (dresult == null) + { + throw new NotImplementedException("Can only do QR factorization for dense vectors at the moment."); + } + + Control.LinearAlgebraProvider.SvdSolveFactored(MatrixU.RowCount, MatrixVT.ColumnCount, ((DenseVector)VectorS).Data, ((DenseMatrix)MatrixU).Data, ((DenseMatrix)MatrixVT).Data, dinput.Data, 1, dresult.Data); + } + } +} diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/ExtensionMethods.cs b/src/Numerics/LinearAlgebra/Double/Factorization/ExtensionMethods.cs index 3eba5289..efd621ea 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/ExtensionMethods.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/ExtensionMethods.cs @@ -64,5 +64,16 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { return Factorization.QR.Create(matrix); } + + /// + /// Computes the SVD decomposition for a matrix. + /// + /// The matrix to factor. + /// Compute the singular U and VT vectors or not. + /// The QR decomposition object. + public static Svd Svd(this Matrix matrix, bool computeVectors) + { + return Factorization.Svd.Create(matrix, computeVectors); + } } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/QR.cs b/src/Numerics/LinearAlgebra/Double/Factorization/QR.cs index f4f216c6..e52326fe 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/QR.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/QR.cs @@ -77,11 +77,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// /// Gets orthogonal Q matrix /// - public virtual Matrix Q + public virtual Matrix Q { get { - return MatrixQ; + return MatrixQ.Clone(); } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/Svd.cs b/src/Numerics/LinearAlgebra/Double/Factorization/Svd.cs new file mode 100644 index 00000000..bb5f9962 --- /dev/null +++ b/src/Numerics/LinearAlgebra/Double/Factorization/Svd.cs @@ -0,0 +1,283 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2010 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. +// + +namespace MathNet.Numerics.LinearAlgebra.Double.Factorization +{ + using System; + using Properties; + + /// + /// A class which encapsulates the functionality of the singular value decomposition (SVD). + /// Suppose M is an m-by-n matrix whose entries are real numbers. + /// Then there exists a factorization of the form M = UΣVT where: + /// - U is an m-by-m unitary matrix; + /// - Σ is m-by-n diagonal matrix with nonnegative real numbers on the diagonal; + /// - VT denotes transpose of V, an n-by-n unitary matrix; + /// Such a factorization is called a singular-value decomposition of M. A common convention is to order the diagonal + /// entries Σ(i,i) in descending order. In this case, the diagonal matrix Σ is uniquely determined + /// by M (though the matrices U and V are not). The diagonal entries of Σ are known as the singular values of M. + /// + /// + /// The computation of the singular value decomposition is done at construction time. + /// + public abstract class Svd : ISolver + { + /// + /// Gets or sets a value indicating whether to compute U and VT matrices during SVD factorization or not + /// + protected bool ComputeVectors + { + get; + set; + } + + /// + /// Gets or sets the singular values (Σ) of matrix in ascending value. + /// + protected Vector VectorS + { + get; + set; + } + + /// + /// Gets or sets left singular vectors (U - m-by-m unitary matrix) + /// + protected Matrix MatrixU + { + get; + set; + } + + /// + /// Gets or sets transpose right singular vectors (transpose of V, an n-by-n unitary matrix + /// + protected Matrix MatrixVT + { + get; + set; + } + + /// + /// Gets the effective numerical matrix rank. + /// + /// The number of non-negligible singular values. + public virtual int Rank + { + get + { + var eps = Math.Pow(2.0, -52.0); + var tol = Math.Max(MatrixU.RowCount, MatrixVT.ColumnCount) * VectorS[0] * eps; + var nm = Math.Min(MatrixU.RowCount, MatrixVT.ColumnCount); + var rank = 0; + for (var h = 0; h < nm; h++) + { + if (VectorS[h] > tol) + { + rank++; + } + } + + return rank; + } + } + + /// + /// Internal method which routes the call to perform the singular value decomposition to the appropriate class. + /// + /// The matrix to factor. + /// Compute the singular U and VT vectors or not. + /// An SVD object. + internal static Svd Create(Matrix matrix, bool computeVectors) + { + var dense = matrix as DenseMatrix; + if (dense != null) + { + return new DenseSvd(dense, computeVectors); + } + + throw new NotImplementedException(); + } + + /// + /// Gets the two norm of the . + /// + /// The 2-norm of the . + public virtual double Norm2 + { + get + { + return VectorS[0]; + } + } + + /// + /// Gets the condition number max(S) / min(S) + /// + /// The condition number. + public virtual double ConditionNumber + { + get + { + var tmp = Math.Min(MatrixU.RowCount, MatrixVT.ColumnCount) - 1; + return VectorS[0] / VectorS[tmp]; + } + } + + /// + /// Gets the determinant of the square matrix for which the SVD was computed. + /// + public virtual double Determinant + { + get + { + if (MatrixU.RowCount != MatrixVT.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSquare); + } + + var det = 1.0; + for (var i = 0; i < VectorS.Count; i++) + { + det *= VectorS[i]; + if (Math.Abs(VectorS[i]).AlmostEqualInDecimalPlaces(0.0, 15)) + { + return 0; + } + } + + return Math.Abs(det); + } + } + + /// Returns the left singular vectors as a . + /// The left singular vectors. The matrix will be null, if computeVectors in the constructor is set to false. + public Matrix U() + { + return ComputeVectors ? MatrixU.Clone() : null; + } + + /// Returns the right singular vectors as a . + /// The right singular vectors. The matrix will be null, if computeVectors in the constructor is set to false. + /// This is the transpose of the V matrix. + public Matrix VT() + { + return ComputeVectors ? MatrixVT.Clone() : null; + } + + /// Returns the singular values as a diagonal . + /// The singular values as a diagonal . + public Matrix W() + { + var rows = MatrixU.RowCount; + var columns = MatrixVT.ColumnCount; + var result = MatrixU.CreateMatrix(rows, columns); + for (var i = 0; i < rows; i++) + { + for (var j = 0; j < columns; j++) + { + if (i == j) + { + result.At(i, i, VectorS[i]); + } + } + } + + return result; + } + + /// Returns the singular values as a . + /// the singular values as a . + public Vector S() + { + return VectorS.Clone(); + } + + /// + /// Solves a system of linear equations, AX = B, with A SVD factorized. + /// + /// The right hand side , B. + /// The left hand side , X. + public virtual Matrix Solve(Matrix input) + { + // Check for proper arguments. + if (input == null) + { + throw new ArgumentNullException("input"); + } + + if (!ComputeVectors) + { + throw new InvalidOperationException(Resources.SingularVectorsNotComputed); + } + + Matrix result = MatrixU.CreateMatrix(MatrixVT.ColumnCount, input.ColumnCount); + Solve(input, result); + return result; + } + + /// + /// Solves a system of linear equations, AX = B, with A SVD factorized. + /// + /// The right hand side , B. + /// The left hand side , X. + public abstract void Solve(Matrix input, Matrix result); + + /// + /// Solves a system of linear equations, Ax = b, with A SVD factorized. + /// + /// The right hand side vector, b. + /// The left hand side , x. + public virtual Vector Solve(Vector input) + { + // Check for proper arguments. + if (input == null) + { + throw new ArgumentNullException("input"); + } + + if (!ComputeVectors) + { + throw new InvalidOperationException(Resources.SingularVectorsNotComputed); + } + + var x = MatrixU.CreateVector(MatrixVT.ColumnCount); + Solve(input, x); + return x; + } + + /// + /// Solves a system of linear equations, Ax = b, with A SVD factorized. + /// + /// The right hand side vector, b. + /// The left hand side , x. + public abstract void Solve(Vector input, Vector result); + } +} diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 39c1e485..7979e567 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -102,8 +102,10 @@ + + diff --git a/src/Numerics/Properties/Resources.Designer.cs b/src/Numerics/Properties/Resources.Designer.cs index ee6eae14..a68693ba 100644 --- a/src/Numerics/Properties/Resources.Designer.cs +++ b/src/Numerics/Properties/Resources.Designer.cs @@ -465,6 +465,15 @@ namespace MathNet.Numerics.Properties { } } + /// + /// Looks up a localized string similar to An algorithm failed to converge.. + /// + internal static string ConvergenceFailed { + get { + return ResourceManager.GetString("ConvergenceFailed", resourceCulture); + } + } + /// /// Looks up a localized string similar to This feature is not implemented yet (but is planned).. /// @@ -591,6 +600,15 @@ namespace MathNet.Numerics.Properties { } } + /// + /// Looks up a localized string similar to The singular vectors were not computed.. + /// + internal static string SingularVectorsNotComputed { + get { + return ResourceManager.GetString("SingularVectorsNotComputed", resourceCulture); + } + } + /// /// Looks up a localized string similar to This special case is not supported yet (but is planned).. /// diff --git a/src/Numerics/Properties/Resources.resx b/src/Numerics/Properties/Resources.resx index 50742f13..3c600f68 100644 --- a/src/Numerics/Properties/Resources.resx +++ b/src/Numerics/Properties/Resources.resx @@ -303,4 +303,10 @@ The integer array does not represent a valid permutation. + + An algorithm failed to converge. + + + The singular vectors were not computed. + \ No newline at end of file diff --git a/src/Silverlight/Silverlight.csproj b/src/Silverlight/Silverlight.csproj index 174d2513..369daafb 100644 --- a/src/Silverlight/Silverlight.csproj +++ b/src/Silverlight/Silverlight.csproj @@ -233,6 +233,9 @@ LinearAlgebra\Double\Factorization\DenseQR.cs + + LinearAlgebra\Double\Factorization\DenseSvd.cs + LinearAlgebra\Double\Factorization\ExtensionMethods.cs @@ -242,6 +245,9 @@ LinearAlgebra\Double\Factorization\QR.cs + + LinearAlgebra\Double\Factorization\Svd.cs + LinearAlgebra\Double\ISolver.cs diff --git a/src/UnitTests/LinearAlgebraTests/Double/Factorization/SvdTests.cs b/src/UnitTests/LinearAlgebraTests/Double/Factorization/SvdTests.cs new file mode 100644 index 00000000..fdc9bb12 --- /dev/null +++ b/src/UnitTests/LinearAlgebraTests/Double/Factorization/SvdTests.cs @@ -0,0 +1,370 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2010 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. +// + +namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization +{ + using System; + using MbUnit.Framework; + using LinearAlgebra.Double; + using LinearAlgebra.Double.Factorization; + + public class SvdTests + { + + [Test] + [ExpectedArgumentNullException] + public void ConstructorNull() + { + new DenseSvd(null, true); + } + + [Test] + [Row(1)] + [Row(10)] + [Row(100)] + public void CanFactorizeIdentity(int order) + { + var I = DenseMatrix.Identity(order); + var factorSvd = I.Svd(true); + + Assert.AreEqual(I.RowCount, factorSvd.U().RowCount); + Assert.AreEqual(I.RowCount, factorSvd.U().ColumnCount); + + Assert.AreEqual(I.ColumnCount, factorSvd.VT().RowCount); + Assert.AreEqual(I.ColumnCount, factorSvd.VT().ColumnCount); + + Assert.AreEqual(I.RowCount, factorSvd.W().RowCount); + Assert.AreEqual(I.ColumnCount, factorSvd.W().ColumnCount); + + for (var i = 0; i < factorSvd.W().RowCount; i++) + { + for (var j = 0; j < factorSvd.W().ColumnCount; j++) + { + Assert.AreEqual(i == j ? 1.0 : 0.0, factorSvd.W()[i, j]); + } + } + } + + [Test] + [Row(1,1)] + [Row(2,2)] + [Row(5,5)] + [Row(10,6)] + [Row(48,52)] + [Row(100,93)] + [MultipleAsserts] + public void CanFactorizeRandomMatrix(int row, int column) + { + var matrixA = MatrixLoader.GenerateRandomMatrix(row, column); + var factorSvd = matrixA.Svd(true); + + // Make sure the U has the right dimensions. + Assert.AreEqual(row, factorSvd.U().RowCount); + Assert.AreEqual(row, factorSvd.U().ColumnCount); + + // Make sure the VT has the right dimensions. + Assert.AreEqual(column, factorSvd.VT().RowCount); + Assert.AreEqual(column, factorSvd.VT().ColumnCount); + + // Make sure the W has the right dimensions. + Assert.AreEqual(row, factorSvd.W().RowCount); + Assert.AreEqual(column, factorSvd.W().ColumnCount); + + // Make sure the U*W*VT is the original matrix. + var matrix = factorSvd.U() * factorSvd.W() * factorSvd.VT(); + for (var i = 0; i < matrix.RowCount; i++) + { + for (var j = 0; j < matrix.ColumnCount; j++) + { + Assert.AreApproximatelyEqual(matrixA[i, j], matrix[i, j], 1.0e-11); + } + } + } + + [Test] + [Row(10, 8)] + [Row(48, 52)] + [Row(100, 93)] + [MultipleAsserts] + public void CheckRankOfNonSquare(int row, int column) + { + var matrixA = MatrixLoader.GenerateRandomMatrix(row, column); + var factorSvd = matrixA.Svd(true); + + var mn = Math.Min(row, column); + Assert.AreEqual(factorSvd.Rank, mn); + } + + [Test] + [Row(1)] + [Row(2)] + [Row(5)] + [Row(9)] + [Row(50)] + [Row(90)] + [MultipleAsserts] + public void CheckRankSquare(int order) + { + var matrixA = MatrixLoader.GenerateRandomMatrix(order, order); + var factorSvd = matrixA.Svd(true); + + if (factorSvd.Determinant != 0) + { + Assert.AreEqual(factorSvd.Rank, order); + } + else + { + Assert.AreEqual(factorSvd.Rank, order - 1); + } + } + + [Test] + [Row(10)] + [Row(50)] + [Row(100)] + [MultipleAsserts] + public void CheckRankOfSquareSingular(int order) + { + var matrixA = new DenseMatrix(order, order); + matrixA[0, 0] = 1; + matrixA[order - 1, order - 1] = 1; + for (var i = 1; i < order - 1; i++) + { + matrixA[i, i - 1] = 1; + matrixA[i, i + 1] = 1; + matrixA[i - 1, i] = 1; + matrixA[i + 1, i] = 1; + } + var factorSvd = matrixA.Svd(true); + + Assert.AreEqual(factorSvd.Determinant, 0); + Assert.AreEqual(factorSvd.Rank, order - 1); + } + + [Test] + [ExpectedException(typeof(InvalidOperationException))] + public void CannotSolveMatrixIfVectorsNotComputed() + { + var matrixA = MatrixLoader.GenerateRandomMatrix(10, 10); + var factorSvd = matrixA.Svd(false); + + var matrixB = MatrixLoader.GenerateRandomMatrix(10, 10); + factorSvd.Solve(matrixB); + } + + [Test] + [ExpectedException(typeof(InvalidOperationException))] + public void CannotSolveVectorIfVectorsNotComputed() + { + var matrixA = MatrixLoader.GenerateRandomMatrix(10, 10); + var factorSvd = matrixA.Svd(false); + + var vectorb = MatrixLoader.GenerateRandomVector(10); + factorSvd.Solve(vectorb); + } + + [Test] + [Row(1, 1)] + [Row(2, 2)] + [Row(5, 5)] + [Row(9, 10)] + [Row(50, 50)] + [Row(90, 100)] + [MultipleAsserts] + public void CanSolveForRandomVector(int row, int column) + { + var matrixA = MatrixLoader.GenerateRandomMatrix(row, column); + var matrixACopy = matrixA.Clone(); + var factorSvd = matrixA.Svd(true); + + var vectorb = MatrixLoader.GenerateRandomVector(row); + var resultx = factorSvd.Solve(vectorb); + + Assert.AreEqual(matrixA.ColumnCount, resultx.Count); + + var bReconstruct = matrixA * resultx; + + // Check the reconstruction. + for (var i = 0; i < vectorb.Count; i++) + { + Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1.0e-11); + } + + // Make sure A didn't change. + for (var i = 0; i < matrixA.RowCount; i++) + { + for (var j = 0; j < matrixA.ColumnCount; j++) + { + Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); + } + } + } + + [Test] + [Row(1, 1)] + [Row(4, 4)] + [Row(7, 8)] + [Row(10, 10)] + [Row(45, 50)] + [Row(80, 100)] + [MultipleAsserts] + public void CanSolveForRandomMatrix(int row, int count) + { + var matrixA = MatrixLoader.GenerateRandomMatrix(row, count); + var matrixACopy = matrixA.Clone(); + var factorSvd = matrixA.Svd(true); + + var matrixB = MatrixLoader.GenerateRandomMatrix(row, count); + var matrixX = factorSvd.Solve(matrixB); + + // The solution X row dimension is equal to the column dimension of A + Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); + // The solution X has the same number of columns as B + Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); + + var matrixBReconstruct = matrixA * matrixX; + + // Check the reconstruction. + for (var i = 0; i < matrixB.RowCount; i++) + { + for (var j = 0; j < matrixB.ColumnCount; j++) + { + Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-11); + } + } + + // Make sure A didn't change. + for (var i = 0; i < matrixA.RowCount; i++) + { + for (var j = 0; j < matrixA.ColumnCount; j++) + { + Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); + } + } + } + + [Test] + [Row(1, 1)] + [Row(2, 2)] + [Row(5, 5)] + [Row(9, 10)] + [Row(50, 50)] + [Row(90, 100)] + [MultipleAsserts] + public void CanSolveForRandomVectorWhenResultVectorGiven(int row, int column) + { + var matrixA = MatrixLoader.GenerateRandomMatrix(row, column); + var matrixACopy = matrixA.Clone(); + var factorSvd = matrixA.Svd(true); + var vectorb = MatrixLoader.GenerateRandomVector(row); + var vectorbCopy = vectorb.Clone(); + var resultx = new DenseVector(column); + factorSvd.Solve(vectorb,resultx); + + var bReconstruct = matrixA * resultx; + + // Check the reconstruction. + for (var i = 0; i < vectorb.Count; i++) + { + Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1.0e-11); + } + + // Make sure A didn't change. + for (var i = 0; i < matrixA.RowCount; i++) + { + for (var j = 0; j < matrixA.ColumnCount; j++) + { + Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); + } + } + + // Make sure b didn't change. + for (var i = 0; i < vectorb.Count; i++) + { + Assert.AreEqual(vectorbCopy[i], vectorb[i]); + } + } + + [Test] + [Row(1, 1)] + [Row(4, 4)] + [Row(7, 8)] + [Row(10, 10)] + [Row(45, 50)] + [Row(80, 100)] + [MultipleAsserts] + public void CanSolveForRandomMatrixWhenResultMatrixGiven(int row, int column) + { + var matrixA = MatrixLoader.GenerateRandomMatrix(row, column); + var matrixACopy = matrixA.Clone(); + var factorSvd = matrixA.Svd(true); + + var matrixB = MatrixLoader.GenerateRandomMatrix(row, column); + var matrixBCopy = matrixB.Clone(); + + var matrixX = new DenseMatrix(column, column); + factorSvd.Solve(matrixB,matrixX); + + // The solution X row dimension is equal to the column dimension of A + Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); + // The solution X has the same number of columns as B + Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); + + var matrixBReconstruct = matrixA * matrixX; + + // Check the reconstruction. + for (var i = 0; i < matrixB.RowCount; i++) + { + for (var j = 0; j < matrixB.ColumnCount; j++) + { + Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-11); + } + } + + // Make sure A didn't change. + for (var i = 0; i < matrixA.RowCount; i++) + { + for (var j = 0; j < matrixA.ColumnCount; j++) + { + Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); + } + } + + // Make sure B didn't change. + for (var i = 0; i < matrixB.RowCount; i++) + { + for (var j = 0; j < matrixB.ColumnCount; j++) + { + Assert.AreEqual(matrixBCopy[i, j], matrixB[i, j]); + } + } + } + } +} diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 67624be3..528eedac 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -90,6 +90,7 @@ +