diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs index 3c10d0f1..ac071b4d 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs @@ -1191,7 +1191,28 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// 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(); + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (aOrder != bRows) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + CholeskyFactor(a, aOrder); + CholeskySolveFactored(a, aOrder, b, bRows, bColumns); } /// @@ -2239,7 +2260,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra break; - // Split at negligible s[l]. + // Split at negligible s[l]. case 2: f = e[l - 1]; e[l - 1] = 0.0; @@ -2267,7 +2288,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra // Perform one qr step. case 3: -// calculate the shift. + // 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])); @@ -2343,7 +2364,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra // Convergence case 4: -// Make the singular value positive + // Make the singular value positive if (stemp[l] < 0.0) { stemp[l] = -stemp[l]; @@ -2421,6 +2442,22 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra work[0] = aRows; } + /// + /// Computes the singular value decomposition of A using CommonParallel where possible. + /// + /// 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 METHOD IS NOT USED. ONLY FOR BENCHMARK PURPOSES. public void SingularValueDecompositionWithCommonParallel(bool computeVectors, double[] a, int aRows, int aColumns, double[] s, double[] u, double[] vt, double[] work) { if (a == null) @@ -5003,7 +5040,94 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// This is equivalent to the GETRF LAPACK routine. public void LUFactor(Complex[] data, int order, int[] ipiv) { - throw new NotImplementedException(); + if (data == null) + { + throw new ArgumentNullException("data"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (data.Length != order * order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "data"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + // Initialize the pivot matrix to the identity permutation. + for (var i = 0; i < order; i++) + { + ipiv[i] = i; + } + + var vLUcolj = new Complex[order]; + + // Outer loop. + for (var j = 0; j < order; j++) + { + var indexj = j * order; + var indexjj = indexj + j; + + // Make a copy of the j-th column to localize references. + for (var i = 0; i < order; i++) + { + vLUcolj[i] = data[indexj + i]; + } + + // Apply previous transformations. + for (var i = 0; i < order; i++) + { + // Most of the time is spent in the following dot product. + var kmax = Math.Min(i, j); + var s = Complex.Zero; + for (var k = 0; k < kmax; k++) + { + s += data[k * order + i] * vLUcolj[k]; + } + + data[indexj + i] = vLUcolj[i] -= s; + } + + // Find pivot and exchange if necessary. + var p = j; + for (var i = j + 1; i < order; i++) + { + if (vLUcolj[i].Magnitude > vLUcolj[p].Magnitude) + { + p = i; + } + } + + if (p != j) + { + for (var k = 0; k < order; k++) + { + var indexk = k * order; + var indexkp = indexk + p; + var indexkj = indexk + j; + var temp = data[indexkp]; + data[indexkp] = data[indexkj]; + data[indexkj] = temp; + } + + ipiv[j] = p; + } + + // Compute multipliers. + if (j < order & data[indexjj] != 0.0) + { + for (var i = j + 1; i < order; i++) + { + data[indexj + i] /= data[indexjj]; + } + } + } } /// @@ -5014,7 +5138,19 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// This is equivalent to the GETRF and GETRI LAPACK routines. public void LUInverse(Complex[] a, int order) { - throw new NotImplementedException(); + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order * order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + var ipiv = new int[order]; + LUFactor(a, order, ipiv); + LUInverseFactored(a, order, ipiv); } /// @@ -5026,7 +5162,34 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// This is equivalent to the GETRI LAPACK routine. public void LUInverseFactored(Complex[] a, int order, int[] ipiv) { - throw new NotImplementedException(); + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order * order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + var inverse = new Complex[a.Length]; + for (var i = 0; i < order; i++) + { + inverse[i + order * i] = Complex.One; + } + + LUSolveFactored(order, a, order, ipiv, inverse); + CommonParallel.For(0, a.Length, index => a[index] = inverse[index]); } /// @@ -5040,7 +5203,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// This is equivalent to the GETRF and GETRI LAPACK routines. public void LUInverse(Complex[] a, int order, Complex[] work) { - throw new NotImplementedException(); + LUInverse(a, order); } /// @@ -5055,191 +5218,1559 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// This is equivalent to the GETRI LAPACK routine. public void LUInverseFactored(Complex[] a, int order, int[] ipiv, Complex[] work) { - throw new NotImplementedException(); + LUInverseFactored(a, order, ipiv); } + /// + /// Solves A*X=B for X using LU factorization. + /// + /// The number of columns of B. + /// The square matrix A. + /// The order of the square matrix . + /// The B matrix. + /// This is equivalent to the GETRF and GETRS LAPACK routines. public void LUSolve(int columnsOfB, Complex[] a, int order, Complex[] b) { - throw new NotImplementedException(); - } - - public void LUSolveFactored(int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b) - { - throw new NotImplementedException(); - } + if (a == null) + { + throw new ArgumentNullException("a"); + } - public void LUSolve(Transpose transposeA, int columnsOfB, Complex[] a, int order, Complex[] b) - { - throw new NotImplementedException(); - } + if (b == null) + { + throw new ArgumentNullException("b"); + } - public void LUSolveFactored(Transpose transposeA, int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b) - { - throw new NotImplementedException(); - } + if (a.Length != order * order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } - /// - /// Computes the Cholesky factorization of A. - /// - /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the - /// the Cholesky factorization. - /// The number of rows or columns in the matrix. - /// This is equivalent to the POTRF LAPACK routine. - public void CholeskyFactor(Complex[] a, int order) - { - throw new NotImplementedException(); - } + if (b.Length != order * columnsOfB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } - /// - /// 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. - public void CholeskySolve(Complex[] a, int aOrder, Complex[] b, int bRows, int bColumns) - { - throw new NotImplementedException(); + var ipiv = new int[order]; + LUFactor(a, order, ipiv); + LUSolveFactored(columnsOfB, a, order, ipiv, b); } /// /// Solves A*X=B for X using a previously factored A matrix. /// - /// The square, positive definite matrix A. - /// The number of rows and columns in A. + /// The number of columns of B. + /// The factored A matrix. + /// The order of the square matrix . + /// The pivot indices of . /// 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) + /// This is equivalent to the GETRS LAPACK routine. + public void LUSolveFactored(int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b) { - throw new NotImplementedException(); - } + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (a.Length != order * order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + if (b.Length != order * columnsOfB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + // Compute the column vector P*B + for (var i = 0; i < ipiv.Length; i++) + { + if (ipiv[i] == i) + { + continue; + } + + var p = ipiv[i]; + for (var j = 0; j < columnsOfB; j++) + { + var indexk = j * order; + var indexkp = indexk + p; + var indexkj = indexk + i; + var temp = b[indexkp]; + b[indexkp] = b[indexkj]; + b[indexkj] = temp; + } + } + + // Solve L*Y = P*B + for (var k = 0; k < order; k++) + { + var korder = k * order; + for (var i = k + 1; i < order; i++) + { + for (var j = 0; j < columnsOfB; j++) + { + var index = j * order; + b[i + index] -= b[k + index] * a[i + korder]; + } + } + } + + // Solve U*X = Y; + for (var k = order - 1; k >= 0; k--) + { + var korder = k + k * order; + for (var j = 0; j < columnsOfB; j++) + { + b[k + j * order] /= a[korder]; + } + + korder = k * order; + for (var i = 0; i < k; i++) + { + for (var j = 0; j < columnsOfB; j++) + { + var index = j * order; + b[i + index] -= b[k + index] * a[i + korder]; + } + } + } - /// - /// 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(); } /// - /// Computes the QR factorization of A. + /// Solves A*X=B for X using LU factorization. /// - /// 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) + /// How to transpose the matrix. + /// The number of columns of B. + /// The square matrix A. + /// The order of the square matrix . + /// The B matrix. + /// This is equivalent to the GETRF and GETRS LAPACK routines. + public void LUSolve(Transpose transposeA, int columnsOfB, Complex[] a, int order, Complex[] b) { - throw new NotImplementedException(); + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (a.Length != order * order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (b.Length != order * columnsOfB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + var ipiv = new int[order]; + LUFactor(a, order, ipiv); + LUSolveFactored(transposeA, columnsOfB, a, order, ipiv, b); } /// - /// Solves A*X=B for X using QR factorization of A. + /// Solves A*X=B for X using a previously factored A matrix. /// - /// 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. + /// How to transpose the matrix. + /// The number of columns of B. + /// The factored A matrix. + /// The order of the square matrix . + /// The pivot indices of . /// 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) + /// This is equivalent to the GETRS LAPACK routine. + public void LUSolveFactored(Transpose transposeA, int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b) { - throw new NotImplementedException(); + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (a.Length != order * order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + if (b.Length != order * columnsOfB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (transposeA == Transpose.Transpose) + { + var aT = new Complex[a.Length]; + for (var i = 0; i < order; i++) + { + for (var j = 0; j < order; j++) + { + aT[(j * order) + i] = a[(i * order) + j]; + } + } + + LUSolveFactored(columnsOfB, aT, order, ipiv, b); + } + else if (transposeA == Transpose.ConjugateTranspose) + { + var acT = new Complex[a.Length]; + for (var i = 0; i < order; i++) + { + for (var j = 0; j < order; j++) + { + acT[(j * order) + i] = a[(i * order) + j].Conjugate(); + } + } + + LUSolveFactored(columnsOfB, acT, order, ipiv, b); + } + else + { + LUSolveFactored(columnsOfB, a, order, ipiv, b); + } } /// - /// Solves A*X=B for X using QR factorization of A. + /// Computes the Cholesky 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) + /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the + /// the Cholesky factorization. + /// The number of rows or columns in the matrix. + /// This is equivalent to the POTRF LAPACK routine. + public void CholeskyFactor(Complex[] a, int order) { - throw new NotImplementedException(); + if (a == null) + { + throw new ArgumentNullException("a"); + } + + for (var i = 0; i < order; i++) + { + var d = Complex.Zero; + int index; + for (var j = 0; j < i; j++) + { + var s = Complex.Zero; + for (var k = 0; k < j; k++) + { + s += a[k * order + i] * a[k * order + j].Conjugate(); + } + + var tmp = j * order; + index = tmp + i; + a[index] = s = (a[index] - s) / a[tmp + j]; + d += s * s.Conjugate(); + } + + index = i * order + i; + d = a[index] - d; + if (d.Real <= 0.0) + { + throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite); + } + + a[index] = d.SquareRoot(); + for (var k = i + 1; k < order; k++) + { + a[k * order + i] = 0.0; + } + } } /// - /// Solves A*X=B for X using a previously QR factored matrix. + /// Solves A*X=B for X using Cholesky factorization. /// - /// 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 square, positive definite matrix A. + /// The number of rows and columns in A. /// 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) + /// 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) { - throw new NotImplementedException(); + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (aOrder != bRows) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + CholeskyFactor(a, aOrder); + CholeskySolveFactored(a, aOrder, b, bRows, bColumns); + } + + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// 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 POTRS LAPACK routine. + public void CholeskySolveFactored(Complex[] a, int aOrder, Complex[] b, int bRows, int bColumns) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (aOrder != bRows) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + CommonParallel.For(0, bColumns, c => + { + var cindex = c * aOrder; + + // Solve L*Y = B; + Complex sum; + for (var i = 0; i < aOrder; i++) + { + sum = b[cindex + i]; + for (var k = i - 1; k >= 0; k--) + { + sum -= a[k * aOrder + i] * b[cindex + k]; + } + + b[cindex + i] = sum / a[i * aOrder + i]; + } + + // Solve L'*X = Y; + for (var i = aOrder - 1; i >= 0; i--) + { + sum = b[cindex + i]; + var iindex = i * aOrder; + for (var k = i + 1; k < aOrder; k++) + { + sum -= a[iindex + k].Conjugate() * b[cindex + k]; + } + + b[cindex + i] = sum / a[iindex + i]; + } + }); + } + + /// + /// 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) + { + 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 Complex[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(Complex[] r, int rRows, int rColumns, Complex[] q, Complex[] work) + { + if (r == null) + { + throw new ArgumentNullException("r"); + } + + if (q == null) + { + throw new ArgumentNullException("q"); + } + + if (work == null) + { + throw new ArgumentNullException("q"); + } + + if (r.Length != rRows * rColumns) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "r"); + } + + if (q.Length != rRows * rRows) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "q"); + } + + if (work.Length < rRows * rRows) + { + work[0] = rRows * rRows; + throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + } + + CommonParallel.For(0, rRows, i => q[(i * rRows) + i] = Complex.One); + + 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); + } + + work[0] = rRows * rRows; + } + + #region QR Factor Helper functions + + /// + /// Perform calculation of Q or R + /// + /// Work array + /// 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(Complex[] work, int workIndex, Complex[] a, int rowCount, int rowStart, int rowEnd, int columnStart, int columnEnd) + { + if (rowStart > rowEnd || columnStart > columnEnd) + { + return; + } + + var vector = new Complex[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].Conjugate() * 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(Complex[] work, Complex[] 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] = Complex.Zero; + }); + + var norm = Complex.Zero; + for (var i = 0; i < rowEnd - rowStart + 1; ++i) + { + var iIndex = tmp + i; + norm += work[iIndex].Magnitude * work[iIndex].Magnitude; + } + + norm = norm.SquareRoot(); + if (rowStart == rowEnd || norm.Magnitude == 0) + { + a[index] = -work[tmp]; + work[tmp] = new Complex(2.0, 0).SquareRoot(); + return; + } + + if (work[tmp].Magnitude != 0.0) + { + norm = norm.Magnitude * (work[tmp] / work[tmp].Magnitude); + } + + a[index] = -norm; + CommonParallel.For(0, rowEnd - rowStart + 1, i => work[tmp + i] /= norm); + work[tmp] += 1.0; + + var s = (1.0 / work[tmp]).SquareRoot(); + CommonParallel.For(0, rowEnd - rowStart + 1, i => work[tmp + i] = work[tmp + i].Conjugate() * s); } - /// - /// 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(); - } + #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(Complex[] r, int rRows, int rColumns, Complex[] q, Complex[] b, int bColumns, Complex[] 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 Complex[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(Complex[] r, int rRows, int rColumns, Complex[] q, Complex[] b, int bColumns, Complex[] x, Complex[] 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) + { + work[0] = rRows * rRows; + throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + } + + QRFactor(r, rRows, rColumns, q, work); + QRSolveFactored(q, r, rRows, rColumns, b, bColumns, x); + + work[0] = rRows * rRows; + } + + /// + /// 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) + { + 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 sol = new Complex[b.Length]; + + // Copy B matrix to "sol", so B data will not be changed + CommonParallel.For(0, b.Length, index => sol[index] = b[index]); + + // Compute Y = transpose(Q)*B + var column = new Complex[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].Conjugate() * 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, Complex[] a, int aRows, int aColumns, Complex[] s, Complex[] u, Complex[] vt) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (s == null) + { + throw new ArgumentNullException("s"); + } + + if (u == null) + { + throw new ArgumentNullException("u"); + } + + if (vt == null) + { + throw new ArgumentNullException("vt"); + } + + if (u.Length != 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 Complex[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, Complex[] a, int aRows, int aColumns, Complex[] s, Complex[] u, Complex[] vt, Complex[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (s == null) + { + throw new ArgumentNullException("s"); + } + + if (u == null) + { + throw new ArgumentNullException("u"); + } + + if (vt == null) + { + throw new ArgumentNullException("vt"); + } + + if (work == null) + { + throw new ArgumentNullException("work"); + } + + if (u.Length != 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) + { + work[0] = aRows; + throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + } + + const int Maxiter = 1000; + + var e = new Complex[aColumns]; + var v = new Complex[vt.Length]; + var stemp = new Complex[Math.Min(aRows + 1, aColumns)]; + + int i, j, l, lp1; + + var cs = 0.0; + var sn = 0.0; + Complex 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 sum = 0.0; + for (i = l; i < aRows; i++) + { + sum += a[(l * aRows) + i].Magnitude * a[(l * aRows) + i].Magnitude; + } + + stemp[l] = Math.Sqrt(sum); + if (stemp[l] != 0.0) + { + if (a[(l * aRows) + l] != 0.0) + { + stemp[l] = stemp[l].Magnitude * (a[(l * aRows) + l] / a[(l * aRows) + l].Magnitude); + } + + // 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[(l * aRows) + i].Conjugate() * a[(j * 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].Conjugate(); + } + + 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].Magnitude * e[i].Magnitude; + } + + e[l] = Math.Sqrt(enorm); + if (e[l] != 0.0) + { + if (e[lp1] != 0.0) + { + e[l] = e[l].Magnitude * (e[lp1] / e[lp1].Magnitude); + } + + // 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].Conjugate(); + + 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]).Conjugate(); + 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[(l * aRows) + i].Conjugate() * u[(j * 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[(l * aColumns) + i].Conjugate() * v[(j * 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++) + { + Complex r; + if (stemp[i] != 0.0) + { + t = stemp[i].Magnitude; + 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].Magnitude; + 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 = stemp[l].Magnitude + stemp[l + 1].Magnitude; + ztest = test + e[l].Magnitude; + 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 + e[ls].Magnitude; + } + + if (ls != l + 1) + { + test = test + e[ls - 1].Magnitude; + } + + ztest = test + stemp[ls].Magnitude; + 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].Real; + e[m - 2] = 0.0; + double t1; + for (var kk = l; kk < m - 1; kk++) + { + k = m - 2 - kk + l; + t1 = stemp[k].Real; + Drotg(ref t1, ref f, ref cs, ref sn); + stemp[k] = t1; + if (k != l) + { + f = -sn * e[k - 1].Real; + 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].Real; + e[l - 1] = 0.0; + for (k = l; k < m; k++) + { + t1 = stemp[k].Real; + Drotg(ref t1, ref f, ref cs, ref sn); + stemp[k] = t1; + f = -sn * e[k].Real; + 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, stemp[m - 1].Magnitude); + scale = Math.Max(scale, stemp[m - 2].Magnitude); + scale = Math.Max(scale, e[m - 2].Magnitude); + scale = Math.Max(scale, stemp[l].Magnitude); + scale = Math.Max(scale, e[l].Magnitude); + var sm = stemp[m - 1].Real / scale; + var smm1 = stemp[m - 2].Real / scale; + var emm1 = e[m - 2].Real / scale; + var sl = stemp[l].Real / scale; + var el = e[l].Real / 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].Real) + (sn * e[k].Real); + e[k] = (cs * e[k]) - (sn * stemp[k]); + g = sn * stemp[k + 1].Real; + 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].Real) + (sn * stemp[k + 1].Real); + stemp[k + 1] = -(sn * e[k]) + (cs * stemp[k + 1]); + g = sn * e[k + 1].Real; + 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].Real < 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].Real >= stemp[l + 1].Real) + { + 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].Conjugate(); + } + } + } - /// - /// 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(); + // 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. + CommonParallel.For(0, Math.Min(aRows, aColumns), index => s[index] = stemp[index]); + + // 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; } /// @@ -5256,7 +6787,64 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// 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(); + 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 Complex[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); } /// @@ -5276,7 +6864,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// 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(); + 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"); + } + + if (work.Length == 0) + { + throw new ArgumentException(Resources.ArgumentSingleDimensionArray, "work"); + } + + if (work.Length < aRows) + { + work[0] = aRows; + throw new ArgumentException(Resources.WorkArrayTooSmall, "work"); + } + + SingularValueDecomposition(true, a, aRows, aColumns, s, u, vt, work); + SvdSolveFactored(aRows, aColumns, s, u, vt, b, bColumns, x); } /// @@ -5292,7 +6947,88 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// 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(); + 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"); + } + + var mn = Math.Min(aRows, aColumns); + var tmp = new Complex[aColumns]; + + for (var k = 0; k < bColumns; k++) + { + for (var j = 0; j < aColumns; j++) + { + var value = Complex.Zero; + if (j < mn) + { + for (var i = 0; i < aRows; i++) + { + value += u[(j * aRows) + i].Conjugate() * b[(k * aRows) + i]; + } + + value /= s[j]; + } + + tmp[j] = value; + } + + for (var j = 0; j < aColumns; j++) + { + var value = Complex.Zero; + for (var i = 0; i < aColumns; i++) + { + value += vt[(j * aColumns) + i].Conjugate() * tmp[i]; + } + + x[(k * aColumns) + j] = value; + } + } } #endregion diff --git a/src/Numerics/LinearAlgebra/Double/SparseVector.cs b/src/Numerics/LinearAlgebra/Double/SparseVector.cs index 238cbe09..6a123451 100644 --- a/src/Numerics/LinearAlgebra/Double/SparseVector.cs +++ b/src/Numerics/LinearAlgebra/Double/SparseVector.cs @@ -768,12 +768,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double NonZerosCount = NonZerosCount }; - Buffer.BlockCopy(_nonZeroIndices, 0, result._nonZeroIndices, 0, NonZerosCount * Constants.SizeOfInt); - - CommonParallel.For( - 0, - NonZerosCount, - index => result._nonZeroValues[index] = -_nonZeroValues[index]); + if (NonZerosCount != 0) + { + CommonParallel.For( + 0, + NonZerosCount, + index => result._nonZeroValues[index] = -_nonZeroValues[index]); + Buffer.BlockCopy(_nonZeroIndices, 0, result._nonZeroIndices, 0, NonZerosCount * Constants.SizeOfInt); + } return result; } @@ -921,7 +923,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double { if (NonZerosCount == 0) { -// No non-zero elements. Return 0 + // No non-zero elements. Return 0 return 0; } @@ -1101,7 +1103,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double } var copy = new SparseVector(Count); - for (int i = 0; i < _nonZeroIndices.Length; i++) + for (var i = 0; i < _nonZeroIndices.Length; i++) { var d = _nonZeroValues[i] * other[_nonZeroIndices[i]]; if (d != 0.0) diff --git a/src/UnitTests/LinearAlgebraTests/Double/DenseVectorTests.cs b/src/UnitTests/LinearAlgebraTests/Double/DenseVectorTests.cs index 46d28e97..4669946e 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/DenseVectorTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/DenseVectorTests.cs @@ -53,8 +53,8 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanCreateDenseVectorFromArray() { - var data = new double[this._data.Length]; - Array.Copy(this._data, data, this._data.Length); + var data = new double[Data.Length]; + Array.Copy(Data, data, Data.Length); var vector = new DenseVector(data); Assert.AreSame(data, vector.Data); @@ -71,12 +71,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanCreateDenseVectorFromAnotherDenseVector() { - var vector = new DenseVector(this._data); + var vector = new DenseVector(Data); var other = new DenseVector(vector); - Assert.AreNotSame(vector, other); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { Assert.AreEqual(vector[i], other[i]); } @@ -86,11 +85,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanCreateDenseVectorFromAnotherVector() { - var vector = (Vector)new DenseVector(this._data); + var vector = (Vector)new DenseVector(Data); var other = new DenseVector(vector); Assert.AreNotSame(vector, other); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { Assert.AreEqual(vector[i], other[i]); } @@ -100,10 +99,10 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanCreateDenseVectorFromUserDefinedVector() { - var vector = new UserDefinedVector(this._data); + var vector = new UserDefinedVector(Data); var other = new DenseVector(vector); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { Assert.AreEqual(vector[i], other[i]); } @@ -131,7 +130,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanConvertDenseVectorToArray() { - var vector = new DenseVector(this._data); + var vector = new DenseVector(Data); var array = (double[])vector; Assert.IsInstanceOfType(typeof(double[]), array); Assert.AreSame(vector.Data, array); @@ -151,9 +150,9 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Test] public void CanCallUnaryPlusOperatorOnDenseVector() { - var vector = new DenseVector(this._data); + var vector = new DenseVector(Data); var other = +vector; - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { Assert.AreEqual(vector[i], other[i]); } @@ -163,26 +162,26 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanAddTwoDenseVectorsUsingOperator() { - var vector = new DenseVector(this._data); - var other = new DenseVector(this._data); + var vector = new DenseVector(Data); + var other = new DenseVector(Data); var result = vector + other; - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i], vector[i], "Making sure the original vector wasn't modified."); - Assert.AreEqual(this._data[i], other[i], "Making sure the original vector wasn't modified."); - Assert.AreEqual(this._data[i] * 2.0, result[i]); + Assert.AreEqual(Data[i], vector[i], "Making sure the original vector wasn't modified."); + Assert.AreEqual(Data[i], other[i], "Making sure the original vector wasn't modified."); + Assert.AreEqual(Data[i] * 2.0, result[i]); } } [Test] public void CanCallUnaryNegationOperatorOnDenseVector() { - var vector = new DenseVector(this._data); + var vector = new DenseVector(Data); var other = -vector; - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(-this._data[i], other[i]); + Assert.AreEqual(-Data[i], other[i]); } } @@ -190,14 +189,14 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanSubtractTwoDenseVectorsUsingOperator() { - var vector = new DenseVector(this._data); - var other = new DenseVector(this._data); + var vector = new DenseVector(Data); + var other = new DenseVector(Data); var result = vector - other; - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i], vector[i], "Making sure the original vector wasn't modified."); - Assert.AreEqual(this._data[i], other[i], "Making sure the original vector wasn't modified."); + Assert.AreEqual(Data[i], vector[i], "Making sure the original vector wasn't modified."); + Assert.AreEqual(Data[i], other[i], "Making sure the original vector wasn't modified."); Assert.AreEqual(0.0, result[i]); } } @@ -206,32 +205,32 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanMultiplyDenseVectorByScalarUsingOperators() { - var vector = new DenseVector(this._data); + var vector = new DenseVector(Data); vector = vector * 2.0; - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i] * 2.0, vector[i]); + Assert.AreEqual(Data[i] * 2.0, vector[i]); } vector = vector * 1.0; - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i] * 2.0, vector[i]); + Assert.AreEqual(Data[i] * 2.0, vector[i]); } - vector = new DenseVector(this._data); + vector = new DenseVector(Data); vector = 2.0 * vector; - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i] * 2.0, vector[i]); + Assert.AreEqual(Data[i] * 2.0, vector[i]); } vector = 1.0 * vector; - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i] * 2.0, vector[i]); + Assert.AreEqual(Data[i] * 2.0, vector[i]); } } @@ -239,26 +238,26 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanDivideDenseVectorByScalarUsingOperators() { - var vector = new DenseVector(this._data); + var vector = new DenseVector(Data); vector = vector / 2.0; - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i] / 2.0, vector[i]); + Assert.AreEqual(Data[i] / 2.0, vector[i]); } vector = vector / 1.0; - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i] / 2.0, vector[i]); + Assert.AreEqual(Data[i] / 2.0, vector[i]); } } [Test] public void CanCalculateOuterProductForDenseVector() { - var vector1 = this.CreateVector(this._data); - var vector2 = this.CreateVector(this._data); + var vector1 = CreateVector(Data); + var vector2 = CreateVector(Data); Matrix m = Vector.OuterProduct(vector1, vector2); for (var i = 0; i < vector1.Count; i++) { @@ -274,7 +273,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double public void OuterProducForDenseVectortWithFirstParameterNullShouldThrowException() { DenseVector vector1 = null; - var vector2 = this.CreateVector(this._data); + var vector2 = CreateVector(Data); Vector.OuterProduct(vector1, vector2); } @@ -282,7 +281,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [ExpectedArgumentNullException] public void OuterProductForDenseVectorWithSecondParameterNullShouldThrowException() { - var vector1 = this.CreateVector(this._data); + var vector1 = CreateVector(Data); DenseVector vector2 = null; Vector.OuterProduct(vector1, vector2); } diff --git a/src/UnitTests/LinearAlgebraTests/Double/SparseVectorTest.cs b/src/UnitTests/LinearAlgebraTests/Double/SparseVectorTest.cs index c08baaaa..ce5f5c1f 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/SparseVectorTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/SparseVectorTest.cs @@ -57,8 +57,8 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanCreateSparseVectorFromArray() { - var data = new double[_data.Length]; - System.Array.Copy(_data, data, _data.Length); + var data = new double[Data.Length]; + Array.Copy(Data, data, Data.Length); var vector = new SparseVector(data); for (var i = 0; i < data.Length; i++) @@ -71,12 +71,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanCreateSparseVectorFromAnotherSparseVector() { - var vector = new SparseVector(_data); + var vector = new SparseVector(Data); var other = new SparseVector(vector); - Assert.AreNotSame(vector, other); - for (var i = 0; i < _data.Length; i++) + for (var i = 0; i < Data.Length; i++) { Assert.AreEqual(vector[i], other[i]); } @@ -86,11 +85,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanCreateSparseVectorFromAnotherVector() { - var vector = (Vector)new SparseVector(_data); + var vector = (Vector)new SparseVector(Data); var other = new SparseVector(vector); Assert.AreNotSame(vector, other); - for (var i = 0; i < _data.Length; i++) + for (var i = 0; i < Data.Length; i++) { Assert.AreEqual(vector[i], other[i]); } @@ -100,10 +99,10 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanCreateSparseVectorFromUserDefinedVector() { - var vector = new UserDefinedVector(_data); + var vector = new UserDefinedVector(Data); var other = new SparseVector(vector); - for (var i = 0; i < _data.Length; i++) + for (var i = 0; i < Data.Length; i++) { Assert.AreEqual(vector[i], other[i]); } @@ -131,7 +130,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanConvertSparseVectorToArray() { - var vector = new SparseVector(_data); + var vector = new SparseVector(Data); var array = vector.ToArray(); Assert.IsInstanceOfType(typeof(double[]), array); Assert.AreElementsEqual(vector, array); @@ -150,9 +149,9 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Test] public void CanCallUnaryPlusOperatorOnSparseVector() { - var vector = new SparseVector(_data); + var vector = new SparseVector(Data); var other = +vector; - for (var i = 0; i < _data.Length; i++) + for (var i = 0; i < Data.Length; i++) { Assert.AreEqual(vector[i], other[i]); } @@ -162,26 +161,26 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanAddTwoSparseVectorsUsingOperator() { - var vector = new SparseVector(_data); - var other = new SparseVector(_data); + var vector = new SparseVector(Data); + var other = new SparseVector(Data); var result = vector + other; - for (var i = 0; i < _data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(_data[i], vector[i], "Making sure the original vector wasn't modified."); - Assert.AreEqual(_data[i], other[i], "Making sure the original vector wasn't modified."); - Assert.AreEqual(_data[i] * 2.0, result[i]); + Assert.AreEqual(Data[i], vector[i], "Making sure the original vector wasn't modified."); + Assert.AreEqual(Data[i], other[i], "Making sure the original vector wasn't modified."); + Assert.AreEqual(Data[i] * 2.0, result[i]); } } [Test] public void CanCallUnaryNegationOperatorOnSparseVector() { - var vector = new SparseVector(_data); + var vector = new SparseVector(Data); var other = -vector; - for (var i = 0; i < _data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(-_data[i], other[i]); + Assert.AreEqual(-Data[i], other[i]); } } @@ -189,14 +188,14 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanSubtractTwoSparseVectorsUsingOperator() { - var vector = new SparseVector(_data); - var other = new SparseVector(_data); + var vector = new SparseVector(Data); + var other = new SparseVector(Data); var result = vector - other; - for (var i = 0; i < _data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(_data[i], vector[i], "Making sure the original vector wasn't modified."); - Assert.AreEqual(_data[i], other[i], "Making sure the original vector wasn't modified."); + Assert.AreEqual(Data[i], vector[i], "Making sure the original vector wasn't modified."); + Assert.AreEqual(Data[i], other[i], "Making sure the original vector wasn't modified."); Assert.AreEqual(0.0, result[i]); } } @@ -205,32 +204,32 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanMultiplySparseVectorByScalarUsingOperators() { - var vector = new SparseVector(_data); + var vector = new SparseVector(Data); vector = vector * 2.0; - for (var i = 0; i < _data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(_data[i] * 2.0, vector[i]); + Assert.AreEqual(Data[i] * 2.0, vector[i]); } vector = vector * 1.0; - for (var i = 0; i < _data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(_data[i] * 2.0, vector[i]); + Assert.AreEqual(Data[i] * 2.0, vector[i]); } - vector = new SparseVector(_data); + vector = new SparseVector(Data); vector = 2.0 * vector; - for (var i = 0; i < _data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(_data[i] * 2.0, vector[i]); + Assert.AreEqual(Data[i] * 2.0, vector[i]); } vector = 1.0 * vector; - for (var i = 0; i < _data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(_data[i] * 2.0, vector[i]); + Assert.AreEqual(Data[i] * 2.0, vector[i]); } } @@ -238,26 +237,26 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanDivideSparseVectorByScalarUsingOperators() { - var vector = new SparseVector(_data); + var vector = new SparseVector(Data); vector = vector / 2.0; - for (var i = 0; i < _data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(_data[i] / 2.0, vector[i]); + Assert.AreEqual(Data[i] / 2.0, vector[i]); } vector = vector / 1.0; - for (var i = 0; i < _data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(_data[i] / 2.0, vector[i]); + Assert.AreEqual(Data[i] / 2.0, vector[i]); } } [Test] public void CanCalculateOuterProductForSparseVector() { - var vector1 = this.CreateVector(this._data); - var vector2 = this.CreateVector(this._data); + var vector1 = CreateVector(Data); + var vector2 = CreateVector(Data); Matrix m = Vector.OuterProduct(vector1, vector2); for (var i = 0; i < vector1.Count; i++) { @@ -273,7 +272,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double public void OuterProducForSparseVectortWithFirstParameterNullShouldThrowException() { SparseVector vector1 = null; - var vector2 = this.CreateVector(this._data); + var vector2 = CreateVector(Data); Vector.OuterProduct(vector1, vector2); } @@ -281,7 +280,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [ExpectedArgumentNullException] public void OuterProductForSparseVectorWithSecondParameterNullShouldThrowException() { - var vector1 = this.CreateVector(this._data); + var vector1 = CreateVector(Data); SparseVector vector2 = null; Vector.OuterProduct(vector1, vector2); } @@ -291,11 +290,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanCreateSparseVectorFromDenseVector() { - var vector = (Vector)new DenseVector(_data); + var vector = (Vector)new DenseVector(Data); var other = new SparseVector(vector); Assert.AreNotSame(vector, other); - for (var i = 0; i < _data.Length; i++) + for (var i = 0; i < Data.Length; i++) { Assert.AreEqual(vector[i], other[i]); } @@ -398,7 +397,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double public void PointwiseMultiplySparseVector() { var zeroArray = new[] { 0.0, 1.0, 0.0, 1.0, 0.0 }; - var vector1 = new SparseVector(this._data); + var vector1 = new SparseVector(Data); var vector2 = new SparseVector(zeroArray); var result = new SparseVector(vector1.Count); @@ -406,7 +405,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double for (var i = 0; i < vector1.Count; i++) { - Assert.AreEqual(this._data[i] * zeroArray[i], result[i]); + Assert.AreEqual(Data[i] * zeroArray[i], result[i]); } Assert.AreEqual(2, result.NonZerosCount); } diff --git a/src/UnitTests/LinearAlgebraTests/Double/VectorTests.Arithmetic.cs b/src/UnitTests/LinearAlgebraTests/Double/VectorTests.Arithmetic.cs index 588bdb36..21108807 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/VectorTests.Arithmetic.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/VectorTests.Arithmetic.cs @@ -35,10 +35,10 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Test] public void CanCallPlus() { - var vector = this.CreateVector(this._data); + var vector = CreateVector(Data); var other = vector.Plus(); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { Assert.AreEqual(vector[i], other[i]); } @@ -55,10 +55,10 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Test] public void CanCallUnaryPlusOperator() { - var vector = this.CreateVector(this._data); + var vector = CreateVector(Data); var other = +vector; - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { Assert.AreEqual(vector[i], other[i]); } @@ -68,18 +68,18 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanAddScalarToVector() { - var copy = this.CreateVector(this._data); + var copy = CreateVector(Data); var vector = copy.Add(2.0); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i] + 2.0, vector[i]); + Assert.AreEqual(Data[i] + 2.0, vector[i]); } vector.Add(0.0); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i] + 2.0, vector[i]); + Assert.AreEqual(Data[i] + 2.0, vector[i]); } } @@ -87,67 +87,67 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanAddScalarToVectorUsingResultVector() { - var vector = this.CreateVector(this._data); - var result = this.CreateVector(this._data.Length); + var vector = CreateVector(Data); + var result = CreateVector(Data.Length); vector.Add(2.0, result); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i], vector[i], "Making sure the original vector wasn't modified."); - Assert.AreEqual(this._data[i] + 2.0, result[i]); + Assert.AreEqual(Data[i], vector[i], "Making sure the original vector wasn't modified."); + Assert.AreEqual(Data[i] + 2.0, result[i]); } vector.Add(0.0, result); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i], result[i]); + Assert.AreEqual(Data[i], result[i]); } } [Test] public void ThrowsArgumentNullExceptionWhenAddingScalarWithNullResultVector() { - var vector = this.CreateVector(this._data.Length); + var vector = CreateVector(Data.Length); Assert.Throws(() => vector.Add(0.0, null)); } [Test] public void ThrowsArgumentExceptionWhenAddingScalarWithWrongSizeResultVector() { - var vector = this.CreateVector(this._data.Length); - var result = this.CreateVector(this._data.Length + 1); + var vector = CreateVector(Data.Length); + var result = CreateVector(Data.Length + 1); Assert.Throws(() => vector.Add(0.0, result)); } [Test] public void ThrowsArgumentNullExceptionWhenAddingTwoVectorsAndOneIsNull() { - var vector = this.CreateVector(this._data); + var vector = CreateVector(Data); Assert.Throws(() => vector.Add(null)); } [Test] public void ThrowsArgumentExceptionWhenAddingTwoVectorsOfDifferingSize() { - var vector = this.CreateVector(this._data.Length); - var other = this.CreateVector(this._data.Length + 1); + var vector = CreateVector(Data.Length); + var other = CreateVector(Data.Length + 1); Assert.Throws(() => vector.Add(other)); } [Test] public void ThrowsArgumentNullExceptionWhenAddingTwoVectorsAndResultIsNull() { - var vector = this.CreateVector(this._data.Length); - var other = this.CreateVector(this._data.Length + 1); + var vector = CreateVector(Data.Length); + var other = CreateVector(Data.Length + 1); Assert.Throws(() => vector.Add(other, null)); } [Test] public void ThrowsArgumentExceptionWhenAddingTwoVectorsAndResultIsDifferentSize() { - var vector = this.CreateVector(this._data.Length); - var other = this.CreateVector(this._data.Length); - var result = this.CreateVector(this._data.Length + 1); + var vector = CreateVector(Data.Length); + var other = CreateVector(Data.Length); + var result = CreateVector(Data.Length + 1); Assert.Throws(() => vector.Add(other, result)); } @@ -155,7 +155,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double public void AdditionOperatorThrowsArgumentNullExpectionIfAVectorIsNull() { Vector a = null; - var b = this.CreateVector(this._data.Length); + var b = CreateVector(Data.Length); Assert.Throws(() => a += b); a = b; @@ -166,21 +166,21 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Test] public void AdditionOperatorThrowsArgumentExpectionIfVectorsAreDifferentSize() { - var a = this.CreateVector(this._data.Length); - var b = this.CreateVector(this._data.Length + 1); + var a = CreateVector(Data.Length); + var b = CreateVector(Data.Length + 1); Assert.Throws(() => a += b); } [Test] public void CanAddTwoVectors() { - var copy = this.CreateVector(this._data); - var other = this.CreateVector(this._data); + var copy = CreateVector(Data); + var other = CreateVector(Data); var vector = copy.Add(other); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i] * 2.0, vector[i]); + Assert.AreEqual(Data[i] * 2.0, vector[i]); } } @@ -188,16 +188,16 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanAddTwoVectorsUsingResultVector() { - var vector = this.CreateVector(this._data); - var other = this.CreateVector(this._data); - var result = this.CreateVector(this._data.Length); + var vector = CreateVector(Data); + var other = CreateVector(Data); + var result = CreateVector(Data.Length); vector.Add(other, result); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i], vector[i], "Making sure the original vector wasn't modified."); - Assert.AreEqual(this._data[i], other[i], "Making sure the original vector wasn't modified."); - Assert.AreEqual(this._data[i] * 2.0, result[i]); + Assert.AreEqual(Data[i], vector[i], "Making sure the original vector wasn't modified."); + Assert.AreEqual(Data[i], other[i], "Making sure the original vector wasn't modified."); + Assert.AreEqual(Data[i] * 2.0, result[i]); } } @@ -205,27 +205,27 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanAddTwoVectorsUsingOperator() { - var vector = this.CreateVector(this._data); - var other = this.CreateVector(this._data); + var vector = CreateVector(Data); + var other = CreateVector(Data); var result = vector + other; - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i], vector[i], "Making sure the original vector wasn't modified."); - Assert.AreEqual(this._data[i], other[i], "Making sure the original vector wasn't modified."); - Assert.AreEqual(this._data[i] * 2.0, result[i]); + Assert.AreEqual(Data[i], vector[i], "Making sure the original vector wasn't modified."); + Assert.AreEqual(Data[i], other[i], "Making sure the original vector wasn't modified."); + Assert.AreEqual(Data[i] * 2.0, result[i]); } } [Test] public void CanAddVectorToItself() { - var copy = this.CreateVector(this._data); + var copy = CreateVector(Data); var vector = copy.Add(copy); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i] * 2.0, vector[i]); + Assert.AreEqual(Data[i] * 2.0, vector[i]); } } @@ -233,14 +233,14 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanAddVectorToItselfUsingResultVector() { - var vector = this.CreateVector(this._data); - var result = this.CreateVector(this._data.Length); + var vector = CreateVector(Data); + var result = CreateVector(Data.Length); vector.Add(vector, result); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i], vector[i], "Making sure the original vector wasn't modified."); - Assert.AreEqual(this._data[i] * 2.0, result[i]); + Assert.AreEqual(Data[i], vector[i], "Making sure the original vector wasn't modified."); + Assert.AreEqual(Data[i] * 2.0, result[i]); } } @@ -248,25 +248,25 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanAddTwoVectorsUsingItselfAsResultVector() { - var vector = this.CreateVector(this._data); - var other = this.CreateVector(this._data); + var vector = CreateVector(Data); + var other = CreateVector(Data); vector.Add(other, vector); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i], other[i], "Making sure the original vector wasn't modified."); - Assert.AreEqual(this._data[i] * 2.0, vector[i]); + Assert.AreEqual(Data[i], other[i], "Making sure the original vector wasn't modified."); + Assert.AreEqual(Data[i] * 2.0, vector[i]); } } [Test] public void CanCallNegate() { - var vector = this.CreateVector(this._data); + var vector = CreateVector(Data); var other = vector.Negate(); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(-this._data[i], other[i]); + Assert.AreEqual(-Data[i], other[i]); } } @@ -281,11 +281,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Test] public void CanCallUnaryNegationOperator() { - var vector = this.CreateVector(this._data); + var vector = CreateVector(Data); var other = -vector; - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(-this._data[i], other[i]); + Assert.AreEqual(-Data[i], other[i]); } } @@ -293,18 +293,18 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanSubtractScalarFromVector() { - var copy = this.CreateVector(this._data); + var copy = CreateVector(Data); var vector = copy.Subtract(2.0); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i] - 2.0, vector[i]); + Assert.AreEqual(Data[i] - 2.0, vector[i]); } vector.Subtract(0.0); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i] - 2.0, vector[i]); + Assert.AreEqual(Data[i] - 2.0, vector[i]); } } @@ -312,67 +312,67 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanSubtractScalarFromVectorUsingResultVector() { - var vector = this.CreateVector(this._data); - var result = this.CreateVector(this._data.Length); + var vector = CreateVector(Data); + var result = CreateVector(Data.Length); vector.Subtract(2.0, result); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i], vector[i], "Making sure the original vector wasn't modified."); - Assert.AreEqual(this._data[i] - 2.0, result[i]); + Assert.AreEqual(Data[i], vector[i], "Making sure the original vector wasn't modified."); + Assert.AreEqual(Data[i] - 2.0, result[i]); } vector.Subtract(0.0, result); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i], result[i]); + Assert.AreEqual(Data[i], result[i]); } } [Test] public void ThrowsArgumentNullExceptionWhenSubtractingScalarWithNullResultVector() { - var vector = this.CreateVector(this._data.Length); + var vector = CreateVector(Data.Length); Assert.Throws(() => vector.Subtract(0.0, null)); } [Test] public void ThrowsArgumentExceptionWhenSubtractingScalarWithWrongSizeResultVector() { - var vector = this.CreateVector(this._data.Length); - var result = this.CreateVector(this._data.Length + 1); + var vector = CreateVector(Data.Length); + var result = CreateVector(Data.Length + 1); Assert.Throws(() => vector.Subtract(0.0, result)); } [Test] public void ThrowsArgumentNullExceptionWhenSubtractingTwoVectorsAndOneIsNull() { - var vector = this.CreateVector(this._data); + var vector = CreateVector(Data); Assert.Throws(() => vector.Subtract(null)); } [Test] public void ThrowsArgumentExceptionWhenSubtractingTwoVectorsOfDifferingSize() { - var vector = this.CreateVector(this._data.Length); - var other = this.CreateVector(this._data.Length + 1); + var vector = CreateVector(Data.Length); + var other = CreateVector(Data.Length + 1); Assert.Throws(() => vector.Subtract(other)); } [Test] public void ThrowsArgumentNullExceptionWhenSubtractingTwoVectorsAndResultIsNull() { - var vector = this.CreateVector(this._data.Length); - var other = this.CreateVector(this._data.Length + 1); + var vector = CreateVector(Data.Length); + var other = CreateVector(Data.Length + 1); Assert.Throws(() => vector.Subtract(other, null)); } [Test] public void ThrowsArgumentExceptionWhenSubtractingTwoVectorsAndResultIsDifferentSize() { - var vector = this.CreateVector(this._data.Length); - var other = this.CreateVector(this._data.Length); - var result = this.CreateVector(this._data.Length + 1); + var vector = CreateVector(Data.Length); + var other = CreateVector(Data.Length); + var result = CreateVector(Data.Length + 1); Assert.Throws(() => vector.Subtract(other, result)); } @@ -380,7 +380,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double public void SubtractionOperatorThrowsArgumentNullExpectionIfAVectorIsNull() { Vector a = null; - var b = this.CreateVector(this._data.Length); + var b = CreateVector(Data.Length); Assert.Throws(() => a -= b); a = b; @@ -391,19 +391,19 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Test] public void SubtractionOperatorThrowsArgumentExpectionIfVectorsAreDifferentSize() { - var a = this.CreateVector(this._data.Length); - var b = this.CreateVector(this._data.Length + 1); + var a = CreateVector(Data.Length); + var b = CreateVector(Data.Length + 1); Assert.Throws(() => a -= b); } [Test] public void CanSubtractTwoVectors() { - var copy = this.CreateVector(this._data); - var other = this.CreateVector(this._data); + var copy = CreateVector(Data); + var other = CreateVector(Data); var vector = copy.Subtract(other); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { Assert.AreEqual(0.0, vector[i]); } @@ -413,15 +413,15 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanSubtractTwoVectorsUsingResultVector() { - var vector = this.CreateVector(this._data); - var other = this.CreateVector(this._data); - var result = this.CreateVector(this._data.Length); + var vector = CreateVector(Data); + var other = CreateVector(Data); + var result = CreateVector(Data.Length); vector.Subtract(other, result); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i], vector[i], "Making sure the original vector wasn't modified."); - Assert.AreEqual(this._data[i], other[i], "Making sure the original vector wasn't modified."); + Assert.AreEqual(Data[i], vector[i], "Making sure the original vector wasn't modified."); + Assert.AreEqual(Data[i], other[i], "Making sure the original vector wasn't modified."); Assert.AreEqual(0.0, result[i]); } } @@ -430,14 +430,14 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanSubtractTwoVectorsUsingOperator() { - var vector = this.CreateVector(this._data); - var other = this.CreateVector(this._data); + var vector = CreateVector(Data); + var other = CreateVector(Data); var result = vector - other; - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i], vector[i], "Making sure the original vector wasn't modified."); - Assert.AreEqual(this._data[i], other[i], "Making sure the original vector wasn't modified."); + Assert.AreEqual(Data[i], vector[i], "Making sure the original vector wasn't modified."); + Assert.AreEqual(Data[i], other[i], "Making sure the original vector wasn't modified."); Assert.AreEqual(0.0, result[i]); } } @@ -445,10 +445,10 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Test] public void CanSubtractVectorFromItself() { - var copy = this.CreateVector(this._data); + var copy = CreateVector(Data); var vector = copy.Subtract(copy); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { Assert.AreEqual(0.0, vector[i]); } @@ -458,13 +458,13 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanSubtractVectorFromItselfUsingResultVector() { - var vector = this.CreateVector(this._data); - var result = this.CreateVector(this._data.Length); + var vector = CreateVector(Data); + var result = CreateVector(Data.Length); vector.Subtract(vector, result); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i], vector[i], "Making sure the original vector wasn't modified."); + Assert.AreEqual(Data[i], vector[i], "Making sure the original vector wasn't modified."); Assert.AreEqual(0.0, result[i]); } } @@ -473,13 +473,13 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanSubtractTwoVectorsUsingItselfAsResultVector() { - var vector = this.CreateVector(this._data); - var other = this.CreateVector(this._data); + var vector = CreateVector(Data); + var other = CreateVector(Data); vector.Subtract(other, vector); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i], other[i], "Making sure the original vector wasn't modified."); + Assert.AreEqual(Data[i], other[i], "Making sure the original vector wasn't modified."); Assert.AreEqual(0.0, vector[i]); } } @@ -488,18 +488,18 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanDivideVectorByScalar() { - var copy = this.CreateVector(this._data); + var copy = CreateVector(Data); var vector = copy.Divide(2.0); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i] / 2.0, vector[i]); + Assert.AreEqual(Data[i] / 2.0, vector[i]); } vector.Divide(1.0); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i] / 2.0, vector[i]); + Assert.AreEqual(Data[i] / 2.0, vector[i]); } } @@ -507,20 +507,20 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanDivideVectorByScalarUsingResultVector() { - var vector = this.CreateVector(this._data); - var result = this.CreateVector(this._data.Length); + var vector = CreateVector(Data); + var result = CreateVector(Data.Length); vector.Divide(2.0, result); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i], vector[i], "Making sure the original vector wasn't modified."); - Assert.AreEqual(this._data[i] / 2.0, result[i]); + Assert.AreEqual(Data[i], vector[i], "Making sure the original vector wasn't modified."); + Assert.AreEqual(Data[i] / 2.0, result[i]); } vector.Divide(1.0, result); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i], result[i]); + Assert.AreEqual(Data[i], result[i]); } } @@ -528,18 +528,18 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanMultiplyVectorByScalar() { - var copy = this.CreateVector(this._data); + var copy = CreateVector(Data); var vector = copy.Multiply(2.0); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i] * 2.0, vector[i]); + Assert.AreEqual(Data[i] * 2.0, vector[i]); } vector.Multiply(1.0); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i] * 2.0, vector[i]); + Assert.AreEqual(Data[i] * 2.0, vector[i]); } } @@ -547,50 +547,50 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanMultiplyVectorByScalarUsingResultVector() { - var vector = this.CreateVector(this._data); - var result = this.CreateVector(this._data.Length); + var vector = CreateVector(Data); + var result = CreateVector(Data.Length); vector.Multiply(2.0, result); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i], vector[i], "Making sure the original vector wasn't modified."); - Assert.AreEqual(this._data[i] * 2.0, result[i]); + Assert.AreEqual(Data[i], vector[i], "Making sure the original vector wasn't modified."); + Assert.AreEqual(Data[i] * 2.0, result[i]); } vector.Multiply(1.0, result); - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i], result[i]); + Assert.AreEqual(Data[i], result[i]); } } [Test] public void ThrowsArgumentNullExceptionWhenMultiplyingScalarWithNullResultVector() { - var vector = this.CreateVector(this._data.Length); + var vector = CreateVector(Data.Length); Assert.Throws(() => vector.Multiply(1.0, null)); } [Test] public void ThrowsArgumentNullExceptionWhenDividingScalarWithNullResultVector() { - var vector = this.CreateVector(this._data.Length); + var vector = CreateVector(Data.Length); Assert.Throws(() => vector.Divide(1.0, null)); } [Test] public void ThrowsArgumentExceptionWhenMultiplyingScalarWithWrongSizeResultVector() { - var vector = this.CreateVector(this._data.Length); - var result = this.CreateVector(this._data.Length + 1); + var vector = CreateVector(Data.Length); + var result = CreateVector(Data.Length + 1); Assert.Throws(() => vector.Multiply(0.0, result)); } [Test] public void ThrowsArgumentExceptionWhenDividingScalarWithWrongSizeResultVector() { - var vector = this.CreateVector(this._data.Length); - var result = this.CreateVector(this._data.Length + 1); + var vector = CreateVector(Data.Length); + var result = CreateVector(Data.Length + 1); Assert.Throws(() => vector.Divide(0.0, result)); } @@ -598,32 +598,32 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanMultiplyVectorByScalarUsingOperators() { - var vector = this.CreateVector(this._data); + var vector = CreateVector(Data); vector = vector * 2.0; - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i] * 2.0, vector[i]); + Assert.AreEqual(Data[i] * 2.0, vector[i]); } vector = vector * 1.0; - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i] * 2.0, vector[i]); + Assert.AreEqual(Data[i] * 2.0, vector[i]); } - vector = this.CreateVector(this._data); + vector = CreateVector(Data); vector = 2.0 * vector; - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i] * 2.0, vector[i]); + Assert.AreEqual(Data[i] * 2.0, vector[i]); } vector = 1.0 * vector; - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i] * 2.0, vector[i]); + Assert.AreEqual(Data[i] * 2.0, vector[i]); } } @@ -631,18 +631,18 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanDivideVectorByScalarUsingOperators() { - var vector = this.CreateVector(this._data); + var vector = CreateVector(Data); vector = vector / 2.0; - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i] / 2.0, vector[i]); + Assert.AreEqual(Data[i] / 2.0, vector[i]); } vector = vector / 1.0; - for (var i = 0; i < this._data.Length; i++) + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(this._data[i] / 2.0, vector[i]); + Assert.AreEqual(Data[i] / 2.0, vector[i]); } } @@ -666,8 +666,8 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Test] public void CanDotProduct() { - var dataA = this.CreateVector(this._data); - var dataB = this.CreateVector(this._data); + var dataA = CreateVector(Data); + var dataB = CreateVector(Data); Assert.AreEqual(55.0, dataA.DotProduct(dataB)); } @@ -676,7 +676,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [ExpectedArgumentNullException] public void DotProductThrowsExceptionWhenArgumentIsNull() { - var dataA = this.CreateVector(this._data); + var dataA = CreateVector(Data); Vector dataB = null; dataA.DotProduct(dataB); @@ -686,8 +686,8 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [ExpectedArgumentException] public void DotProductThrowsExceptionWhenArgumentHasDifferentSize() { - var dataA = this.CreateVector(this._data); - var dataB = this.CreateVector(new double[] { 1, 2, 3, 4, 5, 6 }); + var dataA = CreateVector(Data); + var dataB = CreateVector(new double[] { 1, 2, 3, 4, 5, 6 }); dataA.DotProduct(dataB); } @@ -695,8 +695,8 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Test] public void CanDotProductUsingOperator() { - var dataA = this.CreateVector(this._data); - var dataB = this.CreateVector(this._data); + var dataA = CreateVector(Data); + var dataB = CreateVector(Data); Assert.AreEqual(55.0, dataA * dataB); } @@ -705,7 +705,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [ExpectedArgumentNullException] public void OperatorDotProductThrowsExceptionWhenLeftArgumentIsNull() { - var dataA = this.CreateVector(this._data); + var dataA = CreateVector(Data); Vector dataB = null; var d = dataA * dataB; @@ -716,7 +716,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double public void OperatorDotProductThrowsExceptionWhenRightArgumentIsNull() { Vector dataA = null; - var dataB = this.CreateVector(this._data); + var dataB = CreateVector(Data); var d = dataA * dataB; } @@ -725,8 +725,8 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [ExpectedArgumentException] public void OperatorDotProductThrowsExceptionWhenArgumentHasDifferentSize() { - var dataA = this.CreateVector(this._data); - var dataB = this.CreateVector(new double[] { 1, 2, 3, 4, 5, 6 }); + var dataA = CreateVector(Data); + var dataB = CreateVector(new double[] { 1, 2, 3, 4, 5, 6 }); var d = dataA * dataB; } @@ -734,25 +734,25 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Test] public void PointwiseMultiply() { - var vector1 = this.CreateVector(this._data); + var vector1 = CreateVector(Data); var vector2 = vector1.Clone(); var result = vector1.PointwiseMultiply(vector2); for (var i = 0; i < vector1.Count; i++) { - Assert.AreEqual(this._data[i] * this._data[i], result[i]); + Assert.AreEqual(Data[i] * Data[i], result[i]); } } [Test] public void PointwiseMultiplyUsingResultVector() { - var vector1 = this.CreateVector(this._data); + var vector1 = CreateVector(Data); var vector2 = vector1.Clone(); var result = CreateVector(vector1.Count); vector1.PointwiseMultiply(vector2, result); for (var i = 0; i < vector1.Count; i++) { - Assert.AreEqual(this._data[i] * this._data[i], result[i]); + Assert.AreEqual(Data[i] * Data[i], result[i]); } } @@ -760,7 +760,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [ExpectedArgumentNullException] public void PointwiseMultiplyWithOtherNullShouldThrowException() { - var vector1 = this.CreateVector(this._data); + var vector1 = CreateVector(Data); Vector vector2 = null; var result = CreateVector(vector1.Count); vector1.PointwiseMultiply(vector2, result); @@ -770,7 +770,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [ExpectedArgumentNullException] public void PointwiseMultiplyWithResultNullShouldThrowException() { - var vector1 = this.CreateVector(this._data); + var vector1 = CreateVector(Data); var vector2 = vector1.Clone(); Vector result = null; vector1.PointwiseMultiply(vector2, result); @@ -780,34 +780,34 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [ExpectedArgumentException] public void PointwiseMultiplyWithInvalidResultLengthShouldThrowException() { - var vector1 = this.CreateVector(this._data); + var vector1 = CreateVector(Data); var vector2 = vector1.Clone(); - var result = this.CreateVector(vector1.Count + 1); + var result = CreateVector(vector1.Count + 1); vector1.PointwiseMultiply(vector2, result); } [Test] public void PointWiseDivide() { - var vector1 = this.CreateVector(this._data); + var vector1 = CreateVector(Data); var vector2 = vector1.Clone(); var result = vector1.PointwiseDivide(vector2); for (var i = 0; i < vector1.Count; i++) { - Assert.AreEqual(this._data[i] / this._data[i], result[i]); + Assert.AreEqual(Data[i] / Data[i], result[i]); } } [Test] public void PointWiseDivideUsingResultVector() { - var vector1 = this.CreateVector(this._data); + var vector1 = CreateVector(Data); var vector2 = vector1.Clone(); var result = CreateVector(vector1.Count); vector1.PointwiseDivide(vector2, result); for (var i = 0; i < vector1.Count; i++) { - Assert.AreEqual(this._data[i] / this._data[i], result[i]); + Assert.AreEqual(Data[i] / Data[i], result[i]); } } @@ -815,7 +815,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [ExpectedArgumentNullException] public void PointwiseDivideWithOtherNullShouldThrowException() { - var vector1 = this.CreateVector(this._data); + var vector1 = CreateVector(Data); Vector vector2 = null; var result = CreateVector(vector1.Count); vector1.PointwiseDivide(vector2, result); @@ -825,7 +825,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [ExpectedArgumentNullException] public void PointwiseDivideWithResultNullShouldThrowException() { - var vector1 = this.CreateVector(this._data); + var vector1 = CreateVector(Data); var vector2 = vector1.Clone(); Vector result = null; vector1.PointwiseDivide(vector2, result); @@ -835,17 +835,17 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [ExpectedArgumentException] public void PointwiseDivideWithInvalidResultLengthShouldThrowException() { - var vector1 = this.CreateVector(this._data); + var vector1 = CreateVector(Data); var vector2 = vector1.Clone(); - var result = this.CreateVector(vector1.Count + 1); + var result = CreateVector(vector1.Count + 1); vector1.PointwiseDivide(vector2, result); } [Test] public void CanCalculateOuterProduct() { - var vector1 = this.CreateVector(this._data); - var vector2 = this.CreateVector(this._data); + var vector1 = CreateVector(Data); + var vector2 = CreateVector(Data); Matrix m = Vector.OuterProduct(vector1, vector2); for (var i = 0; i < vector1.Count; i++) { @@ -861,7 +861,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double public void OuterProductWithFirstParameterNullShouldThrowException() { Vector vector1 = null; - var vector2 = this.CreateVector(this._data); + var vector2 = CreateVector(Data); Vector.OuterProduct(vector1, vector2); } @@ -869,7 +869,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [ExpectedArgumentNullException] public void OutercProductWithSecondParameterNullShouldThrowException() { - var vector1 = this.CreateVector(this._data); + var vector1 = CreateVector(Data); Vector vector2 = null; Vector.OuterProduct(vector1, vector2); } @@ -877,8 +877,8 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Test] public void CanCalculateTensorMultiply() { - var vector1 = this.CreateVector(this._data); - var vector2 = this.CreateVector(this._data); + var vector1 = CreateVector(Data); + var vector2 = CreateVector(Data); var m = vector1.TensorMultiply(vector2); for (var i = 0; i < vector1.Count; i++) { @@ -893,7 +893,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [ExpectedArgumentNullException] public void TensorMultiplyWithNullParameterNullShouldThrowException() { - var vector1 = this.CreateVector(this._data); + var vector1 = CreateVector(Data); Vector vector2 = null; vector1.TensorMultiply(vector2); } diff --git a/src/UnitTests/LinearAlgebraTests/Double/VectorTests.Norm.cs b/src/UnitTests/LinearAlgebraTests/Double/VectorTests.Norm.cs index 64d1b9cc..4cddb71c 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/VectorTests.Norm.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/VectorTests.Norm.cs @@ -38,21 +38,21 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Test] public void CanComputeNorm() { - var vector = CreateVector(_data); + var vector = CreateVector(Data); AssertHelpers.AlmostEqual(7.416198487095663, vector.Norm(2), 15); } [Test] public void CanComputeNorm1() { - var vector = CreateVector(_data); + var vector = CreateVector(Data); AssertHelpers.AlmostEqual(15.0, vector.Norm(1), 15); } [Test] public void CanComputeSquareNorm() { - var vector = CreateVector(_data); + var vector = CreateVector(Data); AssertHelpers.AlmostEqual(55.0, vector.Norm(2) * vector.Norm(2), 15); } @@ -63,14 +63,14 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Row(10, 5.0540557845353753)] public void CanComputeNormP(int p, double expected) { - var vector = CreateVector(_data); + var vector = CreateVector(Data); AssertHelpers.AlmostEqual(expected, vector.Norm(p), 15); } [Test] public void CanComputeNormInfinity() { - var vector = CreateVector(_data); + var vector = CreateVector(Data); AssertHelpers.AlmostEqual(5.0, vector.Norm(Double.PositiveInfinity), 15); } @@ -78,7 +78,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanNormalizeVector() { - var vector = CreateVector(_data); + var vector = CreateVector(Data); var result = vector.Normalize(2); AssertHelpers.AlmostEqual(0.134839972492648, result[0], 14); AssertHelpers.AlmostEqual(0.269679944985297, result[1], 14); diff --git a/src/UnitTests/LinearAlgebraTests/Double/VectorTests.cs b/src/UnitTests/LinearAlgebraTests/Double/VectorTests.cs index a340cc1f..dad18d75 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/VectorTests.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/VectorTests.cs @@ -37,18 +37,18 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [TestFixture] public abstract partial class VectorTests { - protected readonly double[] _data = { 1, 2, 3, 4, 5 }; + protected readonly double[] Data = { 1, 2, 3, 4, 5 }; [Test] [MultipleAsserts] public void CanCloneVector() { - var vector = this.CreateVector(this._data); + var vector = CreateVector(Data); var clone = vector.Clone(); Assert.AreNotSame(vector, clone); Assert.AreEqual(vector.Count, clone.Count); - for (var index = 0; index < this._data.Length; index++) + for (var index = 0; index < Data.Length; index++) { Assert.AreEqual(vector[index], clone[index]); } @@ -58,12 +58,12 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanCloneVectorUsingICloneable() { - var vector = this.CreateVector(this._data); + var vector = CreateVector(Data); var clone = (Vector)((ICloneable)vector).Clone(); Assert.AreNotSame(vector, clone); Assert.AreEqual(vector.Count, clone.Count); - for (var index = 0; index < this._data.Length; index++) + for (var index = 0; index < Data.Length; index++) { Assert.AreEqual(vector[index], clone[index]); } @@ -72,7 +72,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Test] public void CanConvertVectorToString() { - var vector = this.CreateVector(this._data); + var vector = CreateVector(Data); var str = vector.ToString(); var sep = CultureInfo.CurrentCulture.TextInfo.ListSeparator; Assert.AreEqual(string.Format("1{0}2{0}3{0}4{0}5", sep), str); @@ -82,8 +82,8 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanCopyPartialVectorToAnother() { - var vector = this.CreateVector(this._data); - var other = this.CreateVector(this._data.Length); + var vector = CreateVector(Data); + var other = CreateVector(Data.Length); vector.CopyTo(other, 2, 2, 2); @@ -98,7 +98,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanCopyPartialVectorToSelf() { - var vector = this.CreateVector(this._data); + var vector = CreateVector(Data); vector.CopyTo(vector, 0, 2, 2); Assert.AreEqual(1.0, vector[0]); @@ -112,12 +112,12 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanCopyVectorToAnother() { - var vector = this.CreateVector(this._data); - var other = this.CreateVector(this._data.Length); + var vector = CreateVector(Data); + var other = CreateVector(Data.Length); vector.CopyTo(other); - for (var index = 0; index < this._data.Length; index++) + for (var index = 0; index < Data.Length; index++) { Assert.AreEqual(vector[index], other[index]); } @@ -132,7 +132,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Test] public void CanCreateVector() { - var expected = this.CreateVector(5); + var expected = CreateVector(5); var actual = expected.CreateVector(5); Assert.AreEqual(expected.GetType(), actual.GetType(), "vectors are same type."); } @@ -140,19 +140,19 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Test] public void CanEnumerateOverVector() { - var vector = this.CreateVector(this._data); - Assert.AreElementsEqual(this._data, vector); + var vector = CreateVector(Data); + Assert.AreElementsEqual(Data, vector); } [Test] [MultipleAsserts] public void CanEnumerateOverVectorUsingIEnumerable() { - var enumerable = (IEnumerable)this.CreateVector(this._data); + var enumerable = (IEnumerable)CreateVector(Data); var index = 0; foreach (var element in enumerable) { - Assert.AreEqual(this._data[index++], (double)element); + Assert.AreEqual(Data[index++], (double)element); } } @@ -160,9 +160,9 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanEquateVectors() { - var vector1 = this.CreateVector(this._data); - var vector2 = this.CreateVector(this._data); - var vector3 = this.CreateVector(4); + var vector1 = CreateVector(Data); + var vector2 = CreateVector(Data); + var vector3 = CreateVector(4); Assert.IsTrue(vector1.Equals(vector1)); Assert.IsTrue(vector1.Equals(vector2)); Assert.IsFalse(vector1.Equals(vector3)); @@ -173,46 +173,46 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void ThrowsArgumentExceptionIfSizeIsNotPositive() { - Assert.Throws(() => this.CreateVector(-1)); - Assert.Throws(() => this.CreateVector(0)); + Assert.Throws(() => CreateVector(-1)); + Assert.Throws(() => CreateVector(0)); } [Test] public void TestingForEqualityWithNonVectorReturnsFalse() { - var vector = this.CreateVector(this._data); + var vector = CreateVector(Data); Assert.IsFalse(vector.Equals(2)); } [Test] public void CanTestForEqualityUsingObjectEquals() { - var vector1 = this.CreateVector(this._data); - var vector2 = this.CreateVector(this._data); + var vector1 = CreateVector(Data); + var vector2 = CreateVector(Data); Assert.IsTrue(vector1.Equals((object)vector2)); } [Test] public void VectorGetHashCode() { - var vector = this.CreateVector(new double[] { 1, 2, 3, 4 }); + var vector = CreateVector(new double[] { 1, 2, 3, 4 }); Assert.AreEqual(2145910784, vector.GetHashCode()); } [Test] public void GetIndexedEnumerator() { - var vector = this.CreateVector(this._data); + var vector = CreateVector(Data); foreach (var pair in vector.GetIndexedEnumerator()) { - Assert.AreEqual(this._data[pair.Key], pair.Value); + Assert.AreEqual(Data[pair.Key], pair.Value); } } [Test] public void CanConvertVectorToArray() { - var vector = this.CreateVector(this._data); + var vector = CreateVector(Data); var array = vector.ToArray(); Assert.AreElementsEqual(vector, array); } @@ -221,7 +221,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanConvertVectorToColumnMatrix() { - var vector = this.CreateVector(this._data); + var vector = CreateVector(Data); var matrix = vector.ToColumnMatrix(); Assert.AreEqual(vector.Count, matrix.RowCount); @@ -237,7 +237,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [MultipleAsserts] public void CanConvertVectorToRowMatrix() { - var vector = this.CreateVector(this._data); + var vector = CreateVector(Data); var matrix = vector.ToRowMatrix(); Assert.AreEqual(vector.Count, matrix.ColumnCount); @@ -252,11 +252,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Test] public void CanSetValues() { - var vector = this.CreateVector(this._data); - vector.SetValues(this._data); - for (var i = 0; i < this._data.Length; i++) + var vector = CreateVector(Data); + vector.SetValues(Data); + for (var i = 0; i < Data.Length; i++) { - Assert.AreEqual(vector[i], this._data[i]); + Assert.AreEqual(vector[i], Data[i]); } } @@ -269,7 +269,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Row(1, -10, ExpectedException = typeof(ArgumentOutOfRangeException))] public void CanCalculateSubVector(int index, int length) { - var vector = this.CreateVector(this._data); + var vector = CreateVector(Data); var sub = vector.SubVector(index, length); Assert.AreEqual(length, sub.Count); for (var i = 0; i < length; i++) @@ -281,81 +281,81 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Test] public void CanFindAbsoluteMinimumIndex() { - var source = this.CreateVector(this._data); - var expected = 0; + var source = CreateVector(Data); + const int Expected = 0; var actual = source.AbsoluteMinimumIndex(); - Assert.AreEqual(expected, actual); + Assert.AreEqual(Expected, actual); } [Test] public void CanFindAbsoluteMinimum() { - var source = this.CreateVector(this._data); - double expected = 1; + var source = CreateVector(Data); + const double Expected = 1; var actual = source.AbsoluteMinimum(); - Assert.AreEqual(expected, actual); + Assert.AreEqual(Expected, actual); } [Test] public void CanFindAbsoluteMaximumIndex() { - var source = this.CreateVector(this._data); - var expected = 4; + var source = CreateVector(Data); + const int Expected = 4; var actual = source.AbsoluteMaximumIndex(); - Assert.AreEqual(expected, actual); + Assert.AreEqual(Expected, actual); } [Test] public void CanFindAbsoluteMaximum() { - var source = this.CreateVector(this._data); - double expected = 5; + var source = CreateVector(Data); + const double Expected = 5; var actual = source.AbsoluteMaximum(); - Assert.AreEqual(expected, actual); + Assert.AreEqual(Expected, actual); } [Test] public void CanFindMaximumIndex() { - var vector = this.CreateVector(this._data); + var vector = CreateVector(Data); - var expected = 4; + const int Expected = 4; var actual = vector.MaximumIndex(); - Assert.AreEqual(expected, actual); + Assert.AreEqual(Expected, actual); } [Test] public void CanFindMaximum() { - var vector = this.CreateVector(this._data); + var vector = CreateVector(Data); - double expected = 5; + const double Expected = 5; var actual = vector.Maximum(); - Assert.AreEqual(expected, actual); + Assert.AreEqual(Expected, actual); } [Test] public void CanFindMinimumIndex() { - var vector = this.CreateVector(this._data); + var vector = CreateVector(Data); - var expected = 0; + const int Expected = 0; var actual = vector.MinimumIndex(); - Assert.AreEqual(expected, actual); + Assert.AreEqual(Expected, actual); } [Test] public void CanFindMinimum() { - var vector = this.CreateVector(this._data); + var vector = CreateVector(Data); - double expected = 1; + const double Expected = 1; var actual = vector.Minimum(); - Assert.AreEqual(expected, actual); + Assert.AreEqual(Expected, actual); } [Test] @@ -367,7 +367,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Row(1, -10, ExpectedException = typeof(ArgumentOutOfRangeException))] public void CanGetSubVector(int index, int length) { - var vector = this.CreateVector(this._data); + var vector = CreateVector(Data); var sub = vector.SubVector(index, length); Assert.AreEqual(length, sub.Count); for (var i = 0; i < length; i++) @@ -382,8 +382,8 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double double[] testData = { -20, -10, 10, 20, 30, }; var vector = CreateVector(testData); var actual = vector.Sum(); - double expected = 30; - Assert.AreEqual(expected, actual); + const double Expected = 30; + Assert.AreEqual(Expected, actual); } [Test] @@ -392,15 +392,15 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double double[] testData = { -20, -10, 10, 20, 30, }; var vector = CreateVector(testData); var actual = vector.SumMagnitudes(); - double expected = 90; - Assert.AreEqual(expected, actual); + const double Expected = 90; + Assert.AreEqual(Expected, actual); } [Test] [ExpectedArgumentNullException] public void CanSetValuesWithNullParameterShouldThrowException() { - var vector = this.CreateVector(this._data); + var vector = CreateVector(Data); vector.SetValues(null); } @@ -408,8 +408,8 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [ExpectedArgumentException] public void CanSetValuesWithNonEqualDataLengthShouldThrowException() { - var vector = this.CreateVector(this._data.Length + 2); - vector.SetValues(this._data); + var vector = CreateVector(Data.Length + 2); + vector.SetValues(Data); } @@ -417,8 +417,8 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [ExpectedArgumentException] public void RandomWithNumberOfElementsLessThanZeroShouldThrowException() { - var vector = this.CreateVector(4); - vector = vector.Random(-2, new ContinuousUniform()); + var vector = CreateVector(4); + vector.Random(-2, new ContinuousUniform()); } [Test]