diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs index 3cea55e4..c972edd1 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs @@ -3,9 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// // Copyright (c) 2009-2010 Math.NET -// // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -14,10 +12,8 @@ // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: -// // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. -// // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND @@ -30,6 +26,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra { using System; + using System.Linq; using System.Numerics; using Properties; using Threading; @@ -72,13 +69,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if (alpha == 1.0) { - CommonParallel.For(0, y.Length, - index => { y[index] += x[index]; }); + CommonParallel.For(0, y.Length, index => { y[index] += x[index]; }); } else { - CommonParallel.For(0, y.Length, - index => { y[index] += alpha * x[index]; }); + CommonParallel.For(0, y.Length, index => { y[index] += alpha * x[index]; }); } } @@ -100,8 +95,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra return; } - CommonParallel.For(0, x.Length, - index => { x[index] = alpha * x[index]; }); + CommonParallel.For(0, x.Length, index => { x[index] = alpha * x[index]; }); } /// @@ -198,8 +192,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - CommonParallel.For(0, y.Length, - index => { result[index] = x[index] - y[index]; }); + CommonParallel.For(0, y.Length, index => { result[index] = x[index] - y[index]; }); } /// @@ -234,8 +227,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - CommonParallel.For(0, y.Length, - index => { result[index] = x[index] * y[index]; }); + CommonParallel.For(0, y.Length, index => { result[index] = x[index] * y[index]; }); } public double MatrixNorm(Norm norm, double[] matrix) @@ -325,7 +317,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra // TODO - For small matrices we should get rid of the parallelism because of startup costs. // Perhaps the following implementations would be a good one // http://blog.feradz.com/2009/01/cache-efficient-matrix-multiplication/ - this.MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, 1.0, xdata, xRows, xColumns, ydata, yRows, yColumns, 0.0, result); + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, 1.0, xdata, xRows, xColumns, ydata, yRows, yColumns, 0.0, result); } /// @@ -700,31 +692,32 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra public void LUFactor(double[] data, int order, int[] ipiv) { // Initialize the pivot matrix to the identity permutation. - for (int i = 0; i < order; i++) + for (var i = 0; i < order; i++) { ipiv[i] = i; } - double[] LUcolj = new double[order]; + var LUcolj = new double[order]; // Outer loop. - for (int j = 0; j < order; j++) + for (var j = 0; j < order; j++) { - int indexj = j * order; - int indexjj = indexj + j; - // Make a copy of the j-th column to localize references. - for (int i = 0; i < order; i++) + 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++) { LUcolj[i] = data[indexj + i]; } // Apply previous transformations. - for (int i = 0; i < order; i++) + for (var i = 0; i < order; i++) { // Most of the time is spent in the following dot product. - int kmax = System.Math.Min(i, j); - double s = 0.0; - for (int k = 0; k < kmax; k++) + var kmax = Math.Min(i, j); + var s = 0.0; + for (var k = 0; k < kmax; k++) { s += data[k * order + i] * LUcolj[k]; } @@ -733,32 +726,34 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } // Find pivot and exchange if necessary. - int p = j; - for (int i = j + 1; i < order; i++) + var p = j; + for (var i = j + 1; i < order; i++) { - if (System.Math.Abs(LUcolj[i]) > System.Math.Abs(LUcolj[p])) + if (Math.Abs(LUcolj[i]) > Math.Abs(LUcolj[p])) { p = i; } } + if (p != j) { - for (int k = 0; k < order; k++) + for (var k = 0; k < order; k++) { - int indexk = k * order; - int indexkp = indexk + p; - int indexkj = indexk + j; - double temp = data[indexkp]; + 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 (int i = j + 1; i < order; i++) + for (var i = j + 1; i < order; i++) { data[indexj + i] /= data[indexjj]; } @@ -894,51 +889,171 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentMatrixDimensions); } - if (Object.ReferenceEquals(a, b)) + if (ReferenceEquals(a, b)) { throw new ArgumentException(Resources.ArgumentReferenceDifferent); } CommonParallel.For(0, bColumns, c => + { + var cindex = c * aOrder; + + // Solve L*Y = B; + double 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] * 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. + /// On exit, A M by M matrix that holds the Q matrix of the + /// QR factorization. + /// This is similar to the GEQRF and ORGQR LAPACK routines. + public void QRFactor(double[] r, double[] q) + { + if (r == null) + { + throw new ArgumentNullException("r"); + } + + if (q == null) + { + throw new ArgumentNullException("q"); + } + + // Matrix Q is square (m x m), where "m" is number of rows of initial matrix; + var rowCount = (int)Math.Sqrt(q.Length); + var columnCount = r.Length / rowCount; + + for (var i = 0; i < rowCount; i++) + { + q[i * rowCount + i] = 1.0; + } + + var minmn = Math.Min(rowCount, columnCount); + var u = new double[minmn][]; + + for (var i = 0; i < minmn; i++) + { + u[i] = GenerateColumn(r, rowCount, i, rowCount - 1, i); + UA(u[i], r, rowCount, i, rowCount - 1, i + 1, columnCount - 1); + } + + for (var i = minmn - 1; i >= 0; i--) { - int cindex = c * aOrder; + UA(u[i], q, rowCount, i, rowCount - 1, i, rowCount - 1); + } + } + + public void QRFactor(double[] r, double[] q, double[] work) + { + throw new NotImplementedException(); + } - // Solve L*Y = B; - double sum; - for (int i = 0; i < aOrder; i++) + #region QR Factor Helper functions + + private static void UA(double[] u, double[] A, int rowCount, int rowStart, int rowEnd, int columnStart, int columnEnd) + { + if (rowStart > rowEnd || columnStart > columnEnd) + { + return; + } + + var vector = new double[columnEnd - columnStart + 1]; + for (var j = columnStart; j <= columnEnd; j++) + { + vector[j - columnStart] = 0.0; + } + + for (var i = rowStart; i <= rowEnd; i++) + { + for (var j = columnStart; j <= columnEnd; j++) { - sum = b[cindex + i]; - for (int k = i - 1; k >= 0; k--) - { - sum -= a[k * aOrder + i] * b[cindex + k]; - } - b[cindex + i] = sum / a[i * aOrder + i]; + vector[j - columnStart] = vector[j - columnStart] + u[i - rowStart] * A[j * rowCount + i]; } + } - // Solve L'*X = Y; - for (int i = aOrder - 1; i >= 0; i--) + for (var i = rowStart; i <= rowEnd; i++) + { + for (var j = columnStart; j <= columnEnd; j++) { - sum = b[cindex + i]; - int iindex = i * aOrder; - for (int k = i + 1; k < aOrder; k++) - { - sum -= a[iindex + k] * b[cindex + k]; - } - b[cindex + i] = sum / a[iindex + i]; + A[j * rowCount + i] = A[j * rowCount + i] - u[i - rowStart] * vector[j - columnStart]; } - }); + } } - public void QRFactor(double[] r, double[] q) + private static double[] GenerateColumn(double[] A, int rowCount, int rowStart, int rowEnd, int column) { - throw new NotImplementedException(); - } + var u = new double[rowEnd - rowStart + 1]; - public void QRFactor(double[] r, double[] q, double[] work) - { - throw new NotImplementedException(); + var tmp = column * rowCount; + var index = tmp + rowStart; + + for (var i = rowStart; i <= rowEnd; i++) + { + u[i - rowStart] = A[tmp + i]; + A[tmp + i] = 0.0; + } + + var norm = u.Sum(t => t * t); + norm = Math.Sqrt(norm); + + if (rowStart == rowEnd || norm == 0) + { + A[index] = -u[0]; + u[0] = Math.Sqrt(2.0); + return u; + } + + var scale = 1.0 / norm; + if (u[0] < 0.0) + { + scale *= -1.0; + } + + A[index] = -1.0 / scale; + for (var i = 0; i < u.Length; i++) + { + u[i] *= scale; + } + + u[0] += 1.0; + var s = Math.Sqrt(1.0 / u[0]); + for (var i = 0; i < u.Length; i++) + { + u[i] *= s; + } + + return u; } + #endregion + public void QRSolve(int columnsOfB, double[] r, double[] q, double[] b, double[] x) { throw new NotImplementedException(); @@ -949,9 +1064,89 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new NotImplementedException(); } + /// + /// Solves A*X=B for X using a previously QR factored matrix. + /// + /// The number of columns of B. + /// The Q matrix obtained by calling . + /// The R matrix obtained by calling . + /// The B matrix. + /// On exit, the solution matrix. public void QRSolveFactored(int columnsOfB, double[] q, double[] r, double[] b, double[] x) { - throw new NotImplementedException(); + 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 (b.Length != x.Length) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + + // Matrix Q is square (m x m), where "m" is number of rows of initial matrix; + var rowCount = (int)Math.Sqrt(q.Length); + var columnCount = r.Length / rowCount; + + // Copy B matrix to result, so B data will not be changed + Buffer.BlockCopy(b, 0, x, 0, b.Length * Constants.SizeOfDouble); + + // Compute Y = transpose(Q)*B + var column = new double[rowCount]; + for (var j = 0; j < columnsOfB; j++) + { + var jm = j * rowCount; + for (var k = 0; k < rowCount; k++) + { + column[k] = x[jm + k]; + } + + for (var i = 0; i < rowCount; i++) + { + double s = 0; + var im = i * rowCount; + for (var k = 0; k < rowCount; k++) + { + s += q[im + k] * column[k]; + } + + x[jm + i] = s; + } + } + + // Solve R*X = Y; + for (var k = columnCount - 1; k >= 0; k--) + { + var km = k * rowCount; + for (var j = 0; j < columnsOfB; j++) + { + x[j * rowCount + k] /= r[km + k]; + } + + for (var i = 0; i < k; i++) + { + for (var j = 0; j < columnsOfB; j++) + { + var jm = j * rowCount; + x[jm + i] -= x[jm + k] * r[km + i]; + } + } + } } public void SinguarValueDecomposition(bool computeVectors, double[] a, double[] s, double[] u, double[] vt) @@ -1067,14 +1262,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - var d = 0.0F; - - for (var i = 0; i < y.Length; i++) - { - d += y[i] * x[i]; - } - - return d; + return y.Select((t, i) => t * x[i]).Sum(); } /// @@ -1269,7 +1457,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra // TODO - For small matrices we should get rid of the parallelism because of startup costs. // Perhaps the following implementations would be a good one // http://blog.feradz.com/2009/01/cache-efficient-matrix-multiplication/ - this.MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, 1.0f, x, xRows, xColumns, y, yRows, yColumns, 0.0f, result); + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, 1.0f, x, xRows, xColumns, y, yRows, yColumns, 0.0f, result); } /// @@ -1403,74 +1591,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if ((int)transposeA > 111 && (int)transposeB > 111) { CommonParallel.For(0, aColumns, j => - { - var jIndex = j * cRows; - for (var i = 0; i != bRows; i++) - { - var iIndex = i * aRows; - float s = 0; - for (var l = 0; l != bColumns; l++) - { - s += adata[iIndex + l] * bdata[l * bRows + j]; - } - - c[jIndex + i] = s; - } - }); + { + var jIndex = j * cRows; + for (var i = 0; i != bRows; i++) + { + var iIndex = i * aRows; + float s = 0; + for (var l = 0; l != bColumns; l++) + { + s += adata[iIndex + l] * bdata[l * bRows + j]; + } + + c[jIndex + i] = s; + } + }); } else if ((int)transposeA > 111) { CommonParallel.For(0, bColumns, j => - { - var jcIndex = j * cRows; - var jbIndex = j * bRows; - for (var i = 0; i != aColumns; i++) - { - var iIndex = i * aRows; - float s = 0; - for (var l = 0; l != aRows; l++) - { - s += adata[iIndex + l] * bdata[jbIndex + l]; - } - - c[jcIndex + i] = s; - } - }); + { + var jcIndex = j * cRows; + var jbIndex = j * bRows; + for (var i = 0; i != aColumns; i++) + { + var iIndex = i * aRows; + float s = 0; + for (var l = 0; l != aRows; l++) + { + s += adata[iIndex + l] * bdata[jbIndex + l]; + } + + c[jcIndex + i] = s; + } + }); } else if ((int)transposeB > 111) { CommonParallel.For(0, bRows, j => - { - var jIndex = j * cRows; - for (var i = 0; i != aRows; i++) - { - float s = 0; - for (var l = 0; l != aColumns; l++) - { - s += adata[l * aRows + i] * bdata[l * bRows + j]; - } - - c[jIndex + i] = s; - } - }); + { + var jIndex = j * cRows; + for (var i = 0; i != aRows; i++) + { + float s = 0; + for (var l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[l * bRows + j]; + } + + c[jIndex + i] = s; + } + }); } else { CommonParallel.For(0, bColumns, j => - { - var jcIndex = j * cRows; - var jbIndex = j * bRows; - for (var i = 0; i != aRows; i++) - { - float s = 0; - for (var l = 0; l != aColumns; l++) - { - s += adata[l * aRows + i] * bdata[jbIndex + l]; - } - - c[jcIndex + i] = s; - } - }); + { + var jcIndex = j * cRows; + var jbIndex = j * bRows; + for (var i = 0; i != aRows; i++) + { + float s = 0; + for (var l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[jbIndex + l]; + } + + c[jcIndex + i] = s; + } + }); } } else @@ -1478,74 +1666,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if ((int)transposeA > 111 && (int)transposeB > 111) { CommonParallel.For(0, aColumns, j => - { - var jIndex = j * cRows; - for (var i = 0; i != bRows; i++) - { - var iIndex = i * aRows; - float s = 0; - for (var l = 0; l != bColumns; l++) - { - s += adata[iIndex + l] * bdata[l * bRows + j]; - } - - c[jIndex + i] = c[jIndex + i] * beta + s; - } - }); + { + var jIndex = j * cRows; + for (var i = 0; i != bRows; i++) + { + var iIndex = i * aRows; + float s = 0; + for (var l = 0; l != bColumns; l++) + { + s += adata[iIndex + l] * bdata[l * bRows + j]; + } + + c[jIndex + i] = c[jIndex + i] * beta + s; + } + }); } else if ((int)transposeA > 111) { CommonParallel.For(0, bColumns, j => - { - var jcIndex = j * cRows; - var jbIndex = j * bRows; - for (var i = 0; i != aColumns; i++) - { - var iIndex = i * aRows; - float s = 0; - for (var l = 0; l != aRows; l++) - { - s += adata[iIndex + l] * bdata[jbIndex + l]; - } - - c[jcIndex + i] = s + c[jcIndex + i] * beta; - } - }); + { + var jcIndex = j * cRows; + var jbIndex = j * bRows; + for (var i = 0; i != aColumns; i++) + { + var iIndex = i * aRows; + float s = 0; + for (var l = 0; l != aRows; l++) + { + s += adata[iIndex + l] * bdata[jbIndex + l]; + } + + c[jcIndex + i] = s + c[jcIndex + i] * beta; + } + }); } else if ((int)transposeB > 111) { - CommonParallel.For(0, bRows, j => - { - var jIndex = j * cRows; - for (var i = 0; i != aRows; i++) - { - float s = 0; - for (var l = 0; l != aColumns; l++) - { - s += adata[l * aRows + i] * bdata[l * bRows + j]; - } - - c[jIndex + i] = s + c[jIndex + i] * beta; - } - }); + CommonParallel.For(0, bRows, j => + { + var jIndex = j * cRows; + for (var i = 0; i != aRows; i++) + { + float s = 0; + for (var l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[l * bRows + j]; + } + + c[jIndex + i] = s + c[jIndex + i] * beta; + } + }); } else { CommonParallel.For(0, bColumns, j => - { - var jcIndex = j * cRows; - var jbIndex = j * bRows; - for (var i = 0; i != aRows; i++) - { - float s = 0; - for (var l = 0; l != aColumns; l++) - { - s += adata[l * aRows + i] * bdata[jbIndex + l]; - } - - c[jcIndex + i] = s + c[jcIndex + i] * beta; - } - }); + { + var jcIndex = j * cRows; + var jbIndex = j * bRows; + for (var i = 0; i != aRows; i++) + { + float s = 0; + for (var l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[jbIndex + l]; + } + + c[jcIndex + i] = s + c[jcIndex + i] * beta; + } + }); } } } @@ -1554,74 +1742,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if ((int)transposeA > 111 && (int)transposeB > 111) { CommonParallel.For(0, aColumns, j => - { - var jIndex = j * cRows; - for (var i = 0; i != bRows; i++) - { - var iIndex = i * aRows; - float s = 0; - for (var l = 0; l != bColumns; l++) - { - s += adata[iIndex + l] * bdata[l * bRows + j]; - } - - c[jIndex + i] = c[jIndex + i] * beta + alpha * s; - } - }); + { + var jIndex = j * cRows; + for (var i = 0; i != bRows; i++) + { + var iIndex = i * aRows; + float s = 0; + for (var l = 0; l != bColumns; l++) + { + s += adata[iIndex + l] * bdata[l * bRows + j]; + } + + c[jIndex + i] = c[jIndex + i] * beta + alpha * s; + } + }); } else if ((int)transposeA > 111) { CommonParallel.For(0, bColumns, j => - { - var jcIndex = j * cRows; - var jbIndex = j * bRows; - for (var i = 0; i != aColumns; i++) - { - var iIndex = i * aRows; - float s = 0; - for (var l = 0; l != aRows; l++) - { - s += adata[iIndex + l] * bdata[jbIndex + l]; - } - - c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta; - } - }); + { + var jcIndex = j * cRows; + var jbIndex = j * bRows; + for (var i = 0; i != aColumns; i++) + { + var iIndex = i * aRows; + float s = 0; + for (var l = 0; l != aRows; l++) + { + s += adata[iIndex + l] * bdata[jbIndex + l]; + } + + c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta; + } + }); } else if ((int)transposeB > 111) { CommonParallel.For(0, bRows, j => - { - var jIndex = j * cRows; - for (var i = 0; i != aRows; i++) - { - float s = 0; - for (var l = 0; l != aColumns; l++) - { - s += adata[l * aRows + i] * bdata[l * bRows + j]; - } - - c[jIndex + i] = alpha * s + c[jIndex + i] * beta; - } - }); + { + var jIndex = j * cRows; + for (var i = 0; i != aRows; i++) + { + float s = 0; + for (var l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[l * bRows + j]; + } + + c[jIndex + i] = alpha * s + c[jIndex + i] * beta; + } + }); } else { CommonParallel.For(0, bColumns, j => - { - var jcIndex = j * cRows; - var jbIndex = j * bRows; - for (var i = 0; i != aRows; i++) - { - float s = 0; - for (var l = 0; l != aColumns; l++) - { - s += adata[l * aRows + i] * bdata[jbIndex + l]; - } + { + var jcIndex = j * cRows; + var jbIndex = j * bRows; + for (var i = 0; i != aRows; i++) + { + float s = 0; + for (var l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[jbIndex + l]; + } - c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta; - } - }); + c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta; + } + }); } } } @@ -2095,7 +2283,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra // TODO - For small matrices we should get rid of the parallelism because of startup costs. // Perhaps the following implementations would be a good one // http://blog.feradz.com/2009/01/cache-efficient-matrix-multiplication/ - this.MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex.One, x, xRows, xColumns, y, yRows, yColumns, Complex.Zero, result); + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex.One, x, xRows, xColumns, y, yRows, yColumns, Complex.Zero, result); } /// @@ -2229,74 +2417,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if ((int)transposeA > 111 && (int)transposeB > 111) { CommonParallel.For(0, aColumns, j => - { - var jIndex = j * cRows; - for (var i = 0; i != bRows; i++) - { - var iIndex = i * aRows; - Complex s = 0; - for (var l = 0; l != bColumns; l++) - { - s += adata[iIndex + l] * bdata[l * bRows + j]; - } - - c[jIndex + i] = s; - } - }); + { + var jIndex = j * cRows; + for (var i = 0; i != bRows; i++) + { + var iIndex = i * aRows; + Complex s = 0; + for (var l = 0; l != bColumns; l++) + { + s += adata[iIndex + l] * bdata[l * bRows + j]; + } + + c[jIndex + i] = s; + } + }); } else if ((int)transposeA > 111) { CommonParallel.For(0, bColumns, j => - { - var jcIndex = j * cRows; - var jbIndex = j * bRows; - for (var i = 0; i != aColumns; i++) - { - var iIndex = i * aRows; - Complex s = 0; - for (var l = 0; l != aRows; l++) - { - s += adata[iIndex + l] * bdata[jbIndex + l]; - } - - c[jcIndex + i] = s; - } - }); + { + var jcIndex = j * cRows; + var jbIndex = j * bRows; + for (var i = 0; i != aColumns; i++) + { + var iIndex = i * aRows; + Complex s = 0; + for (var l = 0; l != aRows; l++) + { + s += adata[iIndex + l] * bdata[jbIndex + l]; + } + + c[jcIndex + i] = s; + } + }); } else if ((int)transposeB > 111) { CommonParallel.For(0, bRows, j => - { - var jIndex = j * cRows; - for (var i = 0; i != aRows; i++) - { - Complex s = 0; - for (var l = 0; l != aColumns; l++) - { - s += adata[l * aRows + i] * bdata[l * bRows + j]; - } - - c[jIndex + i] = s; - } - }); + { + var jIndex = j * cRows; + for (var i = 0; i != aRows; i++) + { + Complex s = 0; + for (var l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[l * bRows + j]; + } + + c[jIndex + i] = s; + } + }); } else { CommonParallel.For(0, bColumns, j => - { - var jcIndex = j * cRows; - var jbIndex = j * bRows; - for (var i = 0; i != aRows; i++) - { - Complex s = 0; - for (var l = 0; l != aColumns; l++) - { - s += adata[l * aRows + i] * bdata[jbIndex + l]; - } - - c[jcIndex + i] = s; - } - }); + { + var jcIndex = j * cRows; + var jbIndex = j * bRows; + for (var i = 0; i != aRows; i++) + { + Complex s = 0; + for (var l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[jbIndex + l]; + } + + c[jcIndex + i] = s; + } + }); } } else @@ -2304,74 +2492,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if ((int)transposeA > 111 && (int)transposeB > 111) { CommonParallel.For(0, aColumns, j => - { - var jIndex = j * cRows; - for (var i = 0; i != bRows; i++) - { - var iIndex = i * aRows; - Complex s = 0; - for (var l = 0; l != bColumns; l++) - { - s += adata[iIndex + l] * bdata[l * bRows + j]; - } - - c[jIndex + i] = c[jIndex + i] * beta + s; - } - }); + { + var jIndex = j * cRows; + for (var i = 0; i != bRows; i++) + { + var iIndex = i * aRows; + Complex s = 0; + for (var l = 0; l != bColumns; l++) + { + s += adata[iIndex + l] * bdata[l * bRows + j]; + } + + c[jIndex + i] = c[jIndex + i] * beta + s; + } + }); } else if ((int)transposeA > 111) { CommonParallel.For(0, bColumns, j => - { - var jcIndex = j * cRows; - var jbIndex = j * bRows; - for (var i = 0; i != aColumns; i++) - { - var iIndex = i * aRows; - Complex s = 0; - for (var l = 0; l != aRows; l++) - { - s += adata[iIndex + l] * bdata[jbIndex + l]; - } - - c[jcIndex + i] = s + c[jcIndex + i] * beta; - } - }); + { + var jcIndex = j * cRows; + var jbIndex = j * bRows; + for (var i = 0; i != aColumns; i++) + { + var iIndex = i * aRows; + Complex s = 0; + for (var l = 0; l != aRows; l++) + { + s += adata[iIndex + l] * bdata[jbIndex + l]; + } + + c[jcIndex + i] = s + c[jcIndex + i] * beta; + } + }); } else if ((int)transposeB > 111) { CommonParallel.For(0, bRows, j => - { - var jIndex = j * cRows; - for (var i = 0; i != aRows; i++) - { - Complex s = 0; - for (var l = 0; l != aColumns; l++) - { - s += adata[l * aRows + i] * bdata[l * bRows + j]; - } - - c[jIndex + i] = s + c[jIndex + i] * beta; - } - }); + { + var jIndex = j * cRows; + for (var i = 0; i != aRows; i++) + { + Complex s = 0; + for (var l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[l * bRows + j]; + } + + c[jIndex + i] = s + c[jIndex + i] * beta; + } + }); } else { CommonParallel.For(0, bColumns, j => - { - var jcIndex = j * cRows; - var jbIndex = j * bRows; - for (var i = 0; i != aRows; i++) - { - Complex s = 0; - for (var l = 0; l != aColumns; l++) - { - s += adata[l * aRows + i] * bdata[jbIndex + l]; - } - - c[jcIndex + i] = s + c[jcIndex + i] * beta; - } - }); + { + var jcIndex = j * cRows; + var jbIndex = j * bRows; + for (var i = 0; i != aRows; i++) + { + Complex s = 0; + for (var l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[jbIndex + l]; + } + + c[jcIndex + i] = s + c[jcIndex + i] * beta; + } + }); } } } @@ -2380,74 +2568,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if ((int)transposeA > 111 && (int)transposeB > 111) { CommonParallel.For(0, aColumns, j => - { - var jIndex = j * cRows; - for (var i = 0; i != bRows; i++) - { - var iIndex = i * aRows; - Complex s = 0; - for (var l = 0; l != bColumns; l++) - { - s += adata[iIndex + l] * bdata[l * bRows + j]; - } - - c[jIndex + i] = c[jIndex + i] * beta + alpha * s; - } - }); + { + var jIndex = j * cRows; + for (var i = 0; i != bRows; i++) + { + var iIndex = i * aRows; + Complex s = 0; + for (var l = 0; l != bColumns; l++) + { + s += adata[iIndex + l] * bdata[l * bRows + j]; + } + + c[jIndex + i] = c[jIndex + i] * beta + alpha * s; + } + }); } else if ((int)transposeA > 111) { CommonParallel.For(0, bColumns, j => - { - var jcIndex = j * cRows; - var jbIndex = j * bRows; - for (var i = 0; i != aColumns; i++) - { - var iIndex = i * aRows; - Complex s = 0; - for (var l = 0; l != aRows; l++) - { - s += adata[iIndex + l] * bdata[jbIndex + l]; - } - - c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta; - } - }); + { + var jcIndex = j * cRows; + var jbIndex = j * bRows; + for (var i = 0; i != aColumns; i++) + { + var iIndex = i * aRows; + Complex s = 0; + for (var l = 0; l != aRows; l++) + { + s += adata[iIndex + l] * bdata[jbIndex + l]; + } + + c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta; + } + }); } else if ((int)transposeB > 111) { - CommonParallel.For(0, bRows, j => - { - var jIndex = j * cRows; - for (var i = 0; i != aRows; i++) - { - Complex s = 0; - for (var l = 0; l != aColumns; l++) - { - s += adata[l * aRows + i] * bdata[l * bRows + j]; - } - - c[jIndex + i] = alpha * s + c[jIndex + i] * beta; - } - }); + CommonParallel.For(0, bRows, j => + { + var jIndex = j * cRows; + for (var i = 0; i != aRows; i++) + { + Complex s = 0; + for (var l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[l * bRows + j]; + } + + c[jIndex + i] = alpha * s + c[jIndex + i] * beta; + } + }); } else { CommonParallel.For(0, bColumns, j => - { - var jcIndex = j * cRows; - var jbIndex = j * bRows; - for (var i = 0; i != aRows; i++) - { - Complex s = 0; - for (var l = 0; l != aColumns; l++) - { - s += adata[l * aRows + i] * bdata[jbIndex + l]; - } + { + var jcIndex = j * cRows; + var jbIndex = j * bRows; + for (var i = 0; i != aRows; i++) + { + Complex s = 0; + for (var l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[jbIndex + l]; + } - c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta; - } - }); + c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta; + } + }); } } } @@ -2886,7 +3074,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra // TODO - For small matrices we should get rid of the parallelism because of startup costs. // Perhaps the following implementations would be a good one // http://blog.feradz.com/2009/01/cache-efficient-matrix-multiplication/ - this.MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex32.One, x, xRows, xColumns, y, yRows, yColumns, Complex32.Zero, result); + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex32.One, x, xRows, xColumns, y, yRows, yColumns, Complex32.Zero, result); } /// @@ -3020,74 +3208,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if ((int)transposeA > 111 && (int)transposeB > 111) { CommonParallel.For(0, aColumns, j => - { - var jIndex = j * cRows; - for (var i = 0; i != bRows; i++) - { - var iIndex = i * aRows; - Complex32 s = 0; - for (var l = 0; l != bColumns; l++) - { - s += adata[iIndex + l] * bdata[l * bRows + j]; - } - - c[jIndex + i] = s; - } - }); + { + var jIndex = j * cRows; + for (var i = 0; i != bRows; i++) + { + var iIndex = i * aRows; + Complex32 s = 0; + for (var l = 0; l != bColumns; l++) + { + s += adata[iIndex + l] * bdata[l * bRows + j]; + } + + c[jIndex + i] = s; + } + }); } else if ((int)transposeA > 111) { CommonParallel.For(0, bColumns, j => - { - var jcIndex = j * cRows; - var jbIndex = j * bRows; - for (var i = 0; i != aColumns; i++) - { - var iIndex = i * aRows; - Complex32 s = 0; - for (var l = 0; l != aRows; l++) - { - s += adata[iIndex + l] * bdata[jbIndex + l]; - } - - c[jcIndex + i] = s; - } - }); + { + var jcIndex = j * cRows; + var jbIndex = j * bRows; + for (var i = 0; i != aColumns; i++) + { + var iIndex = i * aRows; + Complex32 s = 0; + for (var l = 0; l != aRows; l++) + { + s += adata[iIndex + l] * bdata[jbIndex + l]; + } + + c[jcIndex + i] = s; + } + }); } else if ((int)transposeB > 111) { CommonParallel.For(0, bRows, j => - { - var jIndex = j * cRows; - for (var i = 0; i != aRows; i++) - { - Complex32 s = 0; - for (var l = 0; l != aColumns; l++) - { - s += adata[l * aRows + i] * bdata[l * bRows + j]; - } - - c[jIndex + i] = s; - } - }); + { + var jIndex = j * cRows; + for (var i = 0; i != aRows; i++) + { + Complex32 s = 0; + for (var l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[l * bRows + j]; + } + + c[jIndex + i] = s; + } + }); } else { CommonParallel.For(0, bColumns, j => - { - var jcIndex = j * cRows; - var jbIndex = j * bRows; - for (var i = 0; i != aRows; i++) - { - Complex32 s = 0; - for (var l = 0; l != aColumns; l++) - { - s += adata[l * aRows + i] * bdata[jbIndex + l]; - } - - c[jcIndex + i] = s; - } - }); + { + var jcIndex = j * cRows; + var jbIndex = j * bRows; + for (var i = 0; i != aRows; i++) + { + Complex32 s = 0; + for (var l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[jbIndex + l]; + } + + c[jcIndex + i] = s; + } + }); } } else @@ -3095,74 +3283,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if ((int)transposeA > 111 && (int)transposeB > 111) { CommonParallel.For(0, aColumns, j => - { - var jIndex = j * cRows; - for (var i = 0; i != bRows; i++) - { - var iIndex = i * aRows; - Complex32 s = 0; - for (var l = 0; l != bColumns; l++) - { - s += adata[iIndex + l] * bdata[l * bRows + j]; - } - - c[jIndex + i] = c[jIndex + i] * beta + s; - } - }); + { + var jIndex = j * cRows; + for (var i = 0; i != bRows; i++) + { + var iIndex = i * aRows; + Complex32 s = 0; + for (var l = 0; l != bColumns; l++) + { + s += adata[iIndex + l] * bdata[l * bRows + j]; + } + + c[jIndex + i] = c[jIndex + i] * beta + s; + } + }); } else if ((int)transposeA > 111) { CommonParallel.For(0, bColumns, j => - { - var jcIndex = j * cRows; - var jbIndex = j * bRows; - for (var i = 0; i != aColumns; i++) - { - var iIndex = i * aRows; - Complex32 s = 0; - for (var l = 0; l != aRows; l++) - { - s += adata[iIndex + l] * bdata[jbIndex + l]; - } - - c[jcIndex + i] = s + c[jcIndex + i] * beta; - } - }); + { + var jcIndex = j * cRows; + var jbIndex = j * bRows; + for (var i = 0; i != aColumns; i++) + { + var iIndex = i * aRows; + Complex32 s = 0; + for (var l = 0; l != aRows; l++) + { + s += adata[iIndex + l] * bdata[jbIndex + l]; + } + + c[jcIndex + i] = s + c[jcIndex + i] * beta; + } + }); } else if ((int)transposeB > 111) { CommonParallel.For(0, bRows, j => - { - var jIndex = j * cRows; - for (var i = 0; i != aRows; i++) - { - Complex32 s = 0; - for (var l = 0; l != aColumns; l++) - { - s += adata[l * aRows + i] * bdata[l * bRows + j]; - } - - c[jIndex + i] = s + c[jIndex + i] * beta; - } - }); + { + var jIndex = j * cRows; + for (var i = 0; i != aRows; i++) + { + Complex32 s = 0; + for (var l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[l * bRows + j]; + } + + c[jIndex + i] = s + c[jIndex + i] * beta; + } + }); } else { CommonParallel.For(0, bColumns, j => - { - var jcIndex = j * cRows; - var jbIndex = j * bRows; - for (var i = 0; i != aRows; i++) - { - Complex32 s = 0; - for (var l = 0; l != aColumns; l++) - { - s += adata[l * aRows + i] * bdata[jbIndex + l]; - } - - c[jcIndex + i] = s + c[jcIndex + i] * beta; - } - }); + { + var jcIndex = j * cRows; + var jbIndex = j * bRows; + for (var i = 0; i != aRows; i++) + { + Complex32 s = 0; + for (var l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[jbIndex + l]; + } + + c[jcIndex + i] = s + c[jcIndex + i] * beta; + } + }); } } } @@ -3171,74 +3359,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra if ((int)transposeA > 111 && (int)transposeB > 111) { CommonParallel.For(0, aColumns, j => - { - var jIndex = j * cRows; - for (var i = 0; i != bRows; i++) - { - var iIndex = i * aRows; - Complex32 s = 0; - for (var l = 0; l != bColumns; l++) - { - s += adata[iIndex + l] * bdata[l * bRows + j]; - } - - c[jIndex + i] = c[jIndex + i] * beta + alpha * s; - } - }); + { + var jIndex = j * cRows; + for (var i = 0; i != bRows; i++) + { + var iIndex = i * aRows; + Complex32 s = 0; + for (var l = 0; l != bColumns; l++) + { + s += adata[iIndex + l] * bdata[l * bRows + j]; + } + + c[jIndex + i] = c[jIndex + i] * beta + alpha * s; + } + }); } else if ((int)transposeA > 111) { CommonParallel.For(0, bColumns, j => - { - var jcIndex = j * cRows; - var jbIndex = j * bRows; - for (var i = 0; i != aColumns; i++) - { - var iIndex = i * aRows; - Complex32 s = 0; - for (var l = 0; l != aRows; l++) - { - s += adata[iIndex + l] * bdata[jbIndex + l]; - } - - c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta; - } - }); + { + var jcIndex = j * cRows; + var jbIndex = j * bRows; + for (var i = 0; i != aColumns; i++) + { + var iIndex = i * aRows; + Complex32 s = 0; + for (var l = 0; l != aRows; l++) + { + s += adata[iIndex + l] * bdata[jbIndex + l]; + } + + c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta; + } + }); } else if ((int)transposeB > 111) { CommonParallel.For(0, bRows, j => - { - var jIndex = j * cRows; - for (var i = 0; i != aRows; i++) - { - Complex32 s = 0; - for (var l = 0; l != aColumns; l++) - { - s += adata[l * aRows + i] * bdata[l * bRows + j]; - } - - c[jIndex + i] = alpha * s + c[jIndex + i] * beta; - } - }); + { + var jIndex = j * cRows; + for (var i = 0; i != aRows; i++) + { + Complex32 s = 0; + for (var l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[l * bRows + j]; + } + + c[jIndex + i] = alpha * s + c[jIndex + i] * beta; + } + }); } else { CommonParallel.For(0, bColumns, j => - { - var jcIndex = j * cRows; - var jbIndex = j * bRows; - for (var i = 0; i != aRows; i++) - { - Complex32 s = 0; - for (var l = 0; l != aColumns; l++) - { - s += adata[l * aRows + i] * bdata[jbIndex + l]; - } + { + var jcIndex = j * cRows; + var jbIndex = j * bRows; + for (var i = 0; i != aRows; i++) + { + Complex32 s = 0; + for (var l = 0; l != aColumns; l++) + { + s += adata[l * aRows + i] * bdata[jbIndex + l]; + } - c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta; - } - }); + c[jcIndex + i] = alpha * s + c[jcIndex + i] * beta; + } + }); } } } @@ -3389,4 +3577,4 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra #endregion } -} \ No newline at end of file +} diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs b/src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs index 5557e403..800e528f 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs @@ -31,7 +31,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { using System; - using Properties; /// /// A class which encapsulates the functionality of a Cholesky factorization. @@ -42,13 +41,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// The computation of the Cholesky factorization is done at construction time. If the matrix is not symmetric /// or positive definite, the constructor will throw an exception. /// - public abstract class Cholesky + public abstract class Cholesky : ISolver { - /// - /// Stores the Cholesky factor. - /// - protected Matrix mFactor; - /// /// Internal method which routes the call to perform the Cholesky factorization to the appropriate class. /// @@ -66,41 +60,44 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization } /// - /// Returns the lower triangular form of the Cholesky matrix. + /// Gets or sets the lower triangular form of the Cholesky matrix. /// public virtual Matrix Factor { - get { return mFactor; } + get; + protected set; } /// - /// The determinant of the matrix for which the Cholesky matrix was computed. + /// Gets the determinant of the matrix for which the Cholesky matrix was computed. /// public virtual double Determinant { get { - double det = 1.0; - for (int j = 0; j < mFactor.RowCount; j++) + var det = 1.0; + for (var j = 0; j < Factor.RowCount; j++) { - det *= (mFactor[j, j] * mFactor[j, j]); + det *= Factor[j, j] * Factor[j, j]; } + return det; } } /// - /// The log determinant of the matrix for which the Cholesky matrix was computed. + /// Gets the log determinant of the matrix for which the Cholesky matrix was computed. /// public virtual double DeterminantLn { get { - double det = 0.0; - for (int j = 0; j < mFactor.RowCount; j++) + var det = 0.0; + for (var j = 0; j < Factor.RowCount; j++) { - det += (2.0 * Math.Log(mFactor[j, j])); + det += 2.0 * Math.Log(Factor[j, j]); } + return det; } } @@ -118,9 +115,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw new ArgumentNullException("input"); } - var X = input.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(input, X); - return X; + var x = input.CreateMatrix(input.RowCount, input.ColumnCount); + Solve(input, x); + return x; } /// diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs index 4d55a7a5..0b1c0881 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs @@ -49,7 +49,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// Cholesky factorization when the constructor is called and cache it's factorization. /// /// The matrix to factor. - /// If is null. + /// If is null. /// If is not a square matrix. /// If is not positive definite. public DenseCholesky(DenseMatrix matrix) @@ -67,7 +67,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization // Create a new matrix for the Cholesky factor, then perform factorization (while overwriting). var factor = (DenseMatrix)matrix.Clone(); Control.LinearAlgebraProvider.CholeskyFactor(factor.Data, factor.RowCount); - mFactor = factor; + Factor = factor; } /// @@ -99,7 +99,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); } - if (input.RowCount != mFactor.RowCount) + if (input.RowCount != Factor.RowCount) { throw new ArgumentException(Resources.ArgumentMatrixDimensions); } @@ -120,7 +120,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization Buffer.BlockCopy(dinput.Data, 0, dresult.Data, 0, dinput.Data.Length * Constants.SizeOfDouble); // Cholesky solve by overwriting result. - var dfactor = mFactor as DenseMatrix; + var dfactor = (DenseMatrix)Factor; Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.RowCount, dresult.ColumnCount); } @@ -148,7 +148,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - if (input.Count != mFactor.RowCount) + if (input.Count != Factor.RowCount) { throw new ArgumentException(Resources.ArgumentMatrixDimensions); } @@ -169,7 +169,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization Buffer.BlockCopy(dinput.Data, 0, dresult.Data, 0, dinput.Data.Length * Constants.SizeOfDouble); // Cholesky solve by overwriting result. - var dfactor = mFactor as DenseMatrix; + var dfactor = (DenseMatrix)Factor; Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.Count, 1); } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs index deea47c8..3283fb06 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs @@ -48,7 +48,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// LU factorization when the constructor is called and cache it's factorization. /// /// The matrix to factor. - /// If is null. + /// If is null. /// If is not a square matrix. public DenseLU(DenseMatrix matrix) { @@ -63,19 +63,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization } // Create an array for the pivot indices. - mPivots = new int[matrix.RowCount]; + Pivots = new int[matrix.RowCount]; // Create a new matrix for the LU factors, then perform factorization (while overwriting). var factors = (DenseMatrix)matrix.Clone(); - Control.LinearAlgebraProvider.LUFactor(factors.Data, factors.RowCount, mPivots); - mFactors = factors; + Control.LinearAlgebraProvider.LUFactor(factors.Data, factors.RowCount, Pivots); + Factors = factors; } /// - /// Solves a system of linear equations, AX = B, with A LU factorized. + /// Solves a system of linear equations, AX = B, with A LU factorized. /// - /// The right hand side , B. - /// The left hand side , X. + /// The right hand side , B. + /// The left hand side , X. public override void Solve(Matrix input, Matrix result) { // Check for proper arguments. @@ -100,7 +100,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); } - if (input.RowCount != mFactors.RowCount) + if (input.RowCount != Factors.RowCount) { throw new ArgumentException(Resources.ArgumentMatrixDimensions); } @@ -121,16 +121,16 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization Buffer.BlockCopy(dinput.Data, 0, dresult.Data, 0, dinput.Data.Length * Constants.SizeOfDouble); // LU solve by overwriting result. - var dfactors = mFactors as DenseMatrix; + var dfactors = (DenseMatrix)Factors; throw new NotImplementedException(); //Control.LinearAlgebraProvider.LUSolveFactored(dfactors.Data, dfactors.RowCount, dresult.Data, dresult.RowCount, dresult.ColumnCount); } /// - /// Solves a system of linear equations, Ax = b, with A LU factorized. + /// Solves a system of linear equations, Ax = b, with A LU factorized. /// - /// The right hand side vector, b. - /// The left hand side , x. + /// The right hand side vector, b. + /// The left hand side , x. public override void Solve(Vector input, Vector result) { // Check for proper arguments. @@ -150,7 +150,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - if (input.Count != mFactors.RowCount) + if (input.Count != Factors.RowCount) { throw new ArgumentException(Resources.ArgumentMatrixDimensions); } @@ -171,7 +171,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization Buffer.BlockCopy(dinput.Data, 0, dresult.Data, 0, dinput.Data.Length * Constants.SizeOfDouble); // LU solve by overwriting result. - var dfactors = mFactors as DenseMatrix; + var dfactors = Factors as DenseMatrix; throw new NotImplementedException(); //Control.LinearAlgebraProvider.LUSolveFactored(dfactors.Data, dfactors.RowCount, dresult.Data, dresult.Count, 1); } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs new file mode 100644 index 00000000..4fd3370e --- /dev/null +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs @@ -0,0 +1,183 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2010 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics.LinearAlgebra.Double.Factorization +{ + using System; + using Properties; + using Threading; + + /// + /// A class which encapsulates the functionality of the QR decomposition. + /// Any real square matrix A may be decomposed as A = QR where Q is an orthogonal matrix + /// (its columns are orthogonal unit vectors meaning QTQ = I) and R is an upper triangular matrix + /// (also called right triangular matrix). + /// + /// + /// The computation of the QR decomposition is done at construction time by Householder transformation. + /// + public class DenseQR : QR + { + /// + /// Initializes a new instance of the class. This object will compute the + /// QR factorization when the constructor is called and cache it's factorization. + /// + /// The matrix to factor. + /// If is null. + public DenseQR(DenseMatrix matrix) + { + if (matrix == null) + { + throw new ArgumentNullException("matrix"); + } + + if (matrix.RowCount < matrix.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + MatrixR = matrix.Clone(); + MatrixQ = new DenseMatrix(matrix.RowCount); + Control.LinearAlgebraProvider.QRFactor(((DenseMatrix)MatrixR).Data, ((DenseMatrix)MatrixQ).Data); + } + + /// + /// Solves a system of linear equations, AX = B, with A QR factorized. + /// + /// The right hand side , B. + /// The left hand side , X. + public override void Solve(Matrix input, Matrix result) + { + // Check for proper arguments. + if (input == null) + { + throw new ArgumentNullException("input"); + } + + if (result == null) + { + throw new ArgumentNullException("result"); + } + + // The solution X should have the same number of columns as B + if (input.ColumnCount != result.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); + } + + // The dimension compatibility conditions for X = A\B require the two matrices A and B to have the same number of rows + if (MatrixR.RowCount != input.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); + } + + // The solution X row dimension is equal to the column dimension of A + if (MatrixR.ColumnCount != result.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); + } + + var dinput = input as DenseMatrix; + if (dinput == null) + { + throw new NotImplementedException("Can only do QR factorization for dense matrices at the moment."); + } + + var dresult = result as DenseMatrix; + if (dresult == null) + { + throw new NotImplementedException("Can only do QR factorization for dense matrices at the moment."); + } + + var solution = new double[dinput.Data.Length]; + Control.LinearAlgebraProvider.QRSolveFactored(input.ColumnCount, ((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, dinput.Data, solution); + CommonParallel.For( + 0, + dresult.RowCount, + row => + { + for (var col = 0; col < dresult.ColumnCount; col++) + { + dresult[row, col] = solution[row + (col * dinput.RowCount)]; + } + }); + } + + /// + /// Solves a system of linear equations, Ax = b, with A QR factorized. + /// + /// The right hand side vector, b. + /// The left hand side , x. + public override void Solve(Vector input, Vector result) + { + if (input == null) + { + throw new ArgumentNullException("input"); + } + + if (result == null) + { + throw new ArgumentNullException("result"); + } + + // Ax=b where A is an m x n matrix + // Check that b is a column vector with m entries + if (MatrixR.RowCount != input.Count) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + + // Check that x is a column vector with n entries + if (MatrixR.ColumnCount != result.Count) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + var dinput = input as DenseVector; + if (dinput == null) + { + throw new NotImplementedException("Can only do QR factorization for dense vectors at the moment."); + } + + var dresult = result as DenseVector; + if (dresult == null) + { + throw new NotImplementedException("Can only do QR factorization for dense vectors at the moment."); + } + + var solution = new double[dinput.Data.Length]; + Control.LinearAlgebraProvider.QRSolveFactored(1, ((DenseMatrix)MatrixQ).Data, ((DenseMatrix)MatrixR).Data, dinput.Data, solution); + CommonParallel.For( + 0, + dresult.Count, + index => { dresult[index] = solution[index]; }); + } + } +} diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/ExtensionMethods.cs b/src/Numerics/LinearAlgebra/Double/Factorization/ExtensionMethods.cs index 6461a553..3eba5289 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/ExtensionMethods.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/ExtensionMethods.cs @@ -30,9 +30,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { - using System; - using Properties; - /// /// Extension methods which return factorizations for the various matrix classes. /// @@ -43,7 +40,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// /// The matrix to factor. /// The Cholesky decomposition object. - public static Factorization.Cholesky Cholesky(this Matrix matrix) + public static Cholesky Cholesky(this Matrix matrix) { return Factorization.Cholesky.Create(matrix); } @@ -57,5 +54,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { return Factorization.LU.Create(matrix); } + + /// + /// Computes the QR decomposition for a matrix. + /// + /// The matrix to factor. + /// The QR decomposition object. + public static QR QR(this Matrix matrix) + { + return Factorization.QR.Create(matrix); + } } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs b/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs index d895d0b0..465330fd 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs @@ -3,9 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// // Copyright (c) 2009-2010 Math.NET -// // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -14,10 +12,8 @@ // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: -// // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. -// // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND @@ -31,7 +27,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { using System; - using Properties; /// /// A class which encapsulates the functionality of an LU factorization. @@ -43,17 +38,25 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// /// The computation of the LU factorization is done at construction time. /// - public abstract class LU + public abstract class LU : ISolver { /// - /// Stores both the L and U factors in the same matrix.. + /// Gets or sets both the L and U factors in the same matrix. /// - protected Matrix mFactors; + protected Matrix Factors + { + get; + set; + } /// - /// Stores the pivot indices of the LU factorization. + /// Gets or sets the pivot indices of the LU factorization. /// - protected int[] mPivots; + protected int[] Pivots + { + get; + set; + } /// /// Internal method which routes the call to perform the LU factorization to the appropriate class. @@ -72,56 +75,64 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization } /// - /// Returns the lower triangular factor. + /// Gets the lower triangular factor. /// public virtual Matrix L { get { - Matrix result = mFactors.LowerTriangle(); - for (int i = 0; i < result.RowCount; i++) + var result = Factors.LowerTriangle(); + for (var i = 0; i < result.RowCount; i++) { result.At(i, i, 1); } + return result; } } /// - /// Returns the upper triangular factor. + /// Gets the upper triangular factor. /// public virtual Matrix U { - get { return mFactors.UpperTriangle(); } + get + { + return Factors.UpperTriangle(); + } } /// - /// Return the permutation applied to LU factorization. + /// Gets the permutation applied to LU factorization. /// public virtual Permutation P { - get { return Permutation.FromInversions(mPivots); } + get + { + return Permutation.FromInversions(Pivots); + } } /// - /// The determinant of the matrix for which the LU factorization was computed. + /// Gets the determinant of the matrix for which the LU factorization was computed. /// public virtual double Determinant { get { - double det = 1.0; - for (int j = 0; j < mFactors.RowCount; j++) + var det = 1.0; + for (var j = 0; j < Factors.RowCount; j++) { - if (mPivots[j] != j) + if (Pivots[j] != j) { - det = -det * mFactors.At(j, j); + det = -det * Factors.At(j, j); } else { - det *= mFactors.At(j, j); + det *= Factors.At(j, j); } } + return det; } } @@ -139,9 +150,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw new ArgumentNullException("input"); } - var X = input.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(input, X); - return X; + var x = input.CreateMatrix(input.RowCount, input.ColumnCount); + Solve(input, x); + return x; } /// diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/QR.cs b/src/Numerics/LinearAlgebra/Double/Factorization/QR.cs new file mode 100644 index 00000000..dab74ed6 --- /dev/null +++ b/src/Numerics/LinearAlgebra/Double/Factorization/QR.cs @@ -0,0 +1,195 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// Copyright (c) 2009-2010 Math.NET +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics.LinearAlgebra.Double.Factorization +{ + using System; + using Properties; + + /// + /// A class which encapsulates the functionality of the QR decomposition. + /// Any real square matrix A (m x n) may be decomposed as A = QR where Q is an orthogonal matrix (m x m) + /// (its columns are orthogonal unit vectors meaning QTQ = I) and R (m x n) is an upper triangular matrix + /// (also called right triangular matrix). + /// + /// + /// The computation of the QR decomposition is done at construction time by Householder transformation. + /// + public abstract class QR : ISolver + { + /// + /// Gets or sets orthogonal Q matrix + /// + protected virtual Matrix MatrixQ + { + get; + set; + } + + /// + /// Gets or sets upper triangular factor R + /// + protected virtual Matrix MatrixR + { + get; + set; + } + + /// + /// Internal method which routes the call to perform the QR factorization to the appropriate class. + /// + /// The matrix to factor. + /// A QR factorization object. + internal static QR Create(Matrix matrix) + { + var dense = matrix as DenseMatrix; + if (dense != null) + { + return new DenseQR(dense); + } + + throw new NotImplementedException(); + } + + /// + /// Gets orthogonal Q matrix + /// + public virtual Matrix Q + { + get + { + return MatrixQ; + } + } + + /// + /// Gets the upper triangular factor R. + /// + public virtual Matrix R + { + get + { + return MatrixR.UpperTriangle(); + } + } + + /// + /// Gets the determinant of the matrix for which the QR matrix was computed. + /// + public virtual double Determinant + { + get + { + if (MatrixR.RowCount != MatrixR.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSquare); + } + + var det = 1.0; + for (var i = 0; i < MatrixR.ColumnCount; i++) + { + det *= MatrixR.At(i, i); + if (Math.Abs(MatrixR.At(i, i)).AlmostEqualInDecimalPlaces(0.0, 15)) + { + return 0; + } + } + + return Math.Abs(det); + } + } + + /// + /// Gets a value indicating whether the matrix is full rank or not. + /// + /// true if the matrix is full rank; otherwise false. + public virtual bool IsFullRank + { + get + { + for (var i = 0; i < MatrixR.ColumnCount; i++) + { + if (Math.Abs(MatrixR.At(i, i)).AlmostEqualInDecimalPlaces(0.0, 15)) + { + return false; + } + } + + return true; + } + } + + /// + /// Solves a system of linear equations, AX = B, with A QR factorized. + /// + /// The right hand side , B. + /// The left hand side , X. + public virtual Matrix Solve(Matrix input) + { + // Check for proper arguments. + if (input == null) + { + throw new ArgumentNullException("input"); + } + + var matrixX = input.CreateMatrix(MatrixR.ColumnCount, input.ColumnCount); + Solve(input, matrixX); + return matrixX; + } + + /// + /// Solves a system of linear equations, AX = B, with A QR factorized. + /// + /// The right hand side , B. + /// The left hand side , X. + public abstract void Solve(Matrix input, Matrix result); + + /// + /// Solves a system of linear equations, Ax = b, with A QR factorized. + /// + /// The right hand side vector, b. + /// The left hand side , x. + public virtual Vector Solve(Vector input) + { + // Check for proper arguments. + if (input == null) + { + throw new ArgumentNullException("input"); + } + + var x = input.CreateVector(MatrixR.ColumnCount); + Solve(input, x); + return x; + } + + /// + /// Solves a system of linear equations, Ax = b, with A QR factorized. + /// + /// The right hand side vector, b. + /// The left hand side , x. + public abstract void Solve(Vector input, Vector result); + } +} diff --git a/src/Numerics/LinearAlgebra/Double/ISolver.cs b/src/Numerics/LinearAlgebra/Double/ISolver.cs new file mode 100644 index 00000000..fcb8658d --- /dev/null +++ b/src/Numerics/LinearAlgebra/Double/ISolver.cs @@ -0,0 +1,64 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2010 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics.LinearAlgebra.Double +{ + /// + /// Classes that solves a system of linear equations, AX = B. + /// + public interface ISolver + { + /// + /// Solves a system of linear equations, AX = B. + /// + /// The right hand side , B. + /// The left hand side , X. + Matrix Solve(Matrix input); + + /// + /// Solves a system of linear equations, AX = B. + /// + /// The right hand side , B. + /// The left hand side , X. + void Solve(Matrix input, Matrix result); + + /// + /// Solves a system of linear equations, Ax = b + /// + /// The right hand side vector, b. + /// The left hand side , x. + Vector Solve(Vector input); + + /// + /// Solves a system of linear equations, Ax = b. + /// + /// The right hand side vector, b. + /// The left hand side , x. + void Solve(Vector input, Vector result); + } +} \ No newline at end of file diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 3b5114f1..39c1e485 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -101,6 +101,9 @@ + + + diff --git a/src/Silverlight/Silverlight.csproj b/src/Silverlight/Silverlight.csproj index 55d984be..174d2513 100644 --- a/src/Silverlight/Silverlight.csproj +++ b/src/Silverlight/Silverlight.csproj @@ -230,12 +230,21 @@ LinearAlgebra\Double\Factorization\DenseLU.cs + + LinearAlgebra\Double\Factorization\DenseQR.cs + LinearAlgebra\Double\Factorization\ExtensionMethods.cs LinearAlgebra\Double\Factorization\LU.cs + + LinearAlgebra\Double\Factorization\QR.cs + + + LinearAlgebra\Double\ISolver.cs + LinearAlgebra\Double\Matrix.Arithmetic.cs diff --git a/src/UnitTests/LinearAlgebraTests/Double/Factorization/QRTests.cs b/src/UnitTests/LinearAlgebraTests/Double/Factorization/QRTests.cs new file mode 100644 index 00000000..51088763 --- /dev/null +++ b/src/UnitTests/LinearAlgebraTests/Double/Factorization/QRTests.cs @@ -0,0 +1,317 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2010 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization +{ + using MbUnit.Framework; + using LinearAlgebra.Double; + using LinearAlgebra.Double.Factorization; + + public class QRTests + { + + [Test] + [ExpectedArgumentNullException] + public void ConstructorNull() + { + new DenseQR(null); + } + + [Test] + [ExpectedArgumentException] + public void WideMatrixThrowsInvalidMatrixOperationException() + { + new DenseQR(new DenseMatrix(3, 4)); + } + + [Test] + [Row(1)] + [Row(10)] + [Row(100)] + public void CanFactorizeIdentity(int order) + { + var I = DenseMatrix.Identity(order); + var factorQR = I.QR(); + + Assert.AreEqual(I.RowCount, factorQR.R.RowCount); + Assert.AreEqual(I.ColumnCount, factorQR.R.ColumnCount); + + for (var i = 0; i < factorQR.R.RowCount; i++) + { + for (var j = 0; j < factorQR.R.ColumnCount; j++) + { + if (i == j) + { + Assert.AreEqual(-1.0, factorQR.R[i, j]); + } + else + { + Assert.AreEqual(0.0, factorQR.R[i, j]); + } + } + } + } + + + [Test] + [Row(1)] + [Row(10)] + [Row(100)] + public void IdentityDeterminantIsOne(int order) + { + var I = DenseMatrix.Identity(order); + var factorQR = I.QR(); + Assert.AreEqual(1.0, factorQR.Determinant); + } + + [Test] + [Row(1,1)] + [Row(2,2)] + [Row(5,5)] + [Row(10,6)] + [Row(50,48)] + [Row(100,98)] + [MultipleAsserts] + public void CanFactorizeRandomMatrix(int row, int column) + { + var matrixA = MatrixLoader.GenerateRandomMatrix(row, column); + var factorQR = matrixA.QR(); + + // Make sure the R has the right dimensions. + Assert.AreEqual(row, factorQR.R.RowCount); + Assert.AreEqual(column, factorQR.R.ColumnCount); + + // Make sure the Q has the right dimensions. + Assert.AreEqual(row, factorQR.Q.RowCount); + Assert.AreEqual(row, factorQR.Q.ColumnCount); + + // Make sure the R factor is upper triangular. + for (var i = 0; i < factorQR.R.RowCount; i++) + { + for (var j = 0; j < factorQR.R.ColumnCount; j++) + { + if (i > j) + { + Assert.AreEqual(0.0, factorQR.R[i, j]); + } + } + } + + // Make sure the Q*R is the original matrix. + var matrixQfromR = factorQR.Q * factorQR.R; + for (int i = 0; i < matrixQfromR.RowCount; i++) + { + for (int j = 0; j < matrixQfromR.ColumnCount; j++) + { + Assert.AreApproximatelyEqual(matrixA[i, j], matrixQfromR[i, j], 1.0e-11); + } + } + } + + [Test] + [Row(1)] + [Row(2)] + [Row(5)] + [Row(10)] + [Row(50)] + [Row(100)] + [MultipleAsserts] + public void CanSolveForRandomVector(int order) + { + var matrixA = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixACopy = matrixA.Clone(); + var factorQR = matrixA.QR(); + + var vectorb = MatrixLoader.GenerateRandomVector(order); + var resultx = factorQR.Solve(vectorb); + + Assert.AreEqual(matrixA.ColumnCount, resultx.Count); + + var bReconstruct = matrixA * resultx; + + // Check the reconstruction. + for (var i = 0; i < order; i++) + { + Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1.0e-11); + } + + // Make sure A didn't change. + for (var i = 0; i < matrixA.RowCount; i++) + { + for (var j = 0; j < matrixA.ColumnCount; j++) + { + Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); + } + } + } + + [Test] + [Row(1)] + [Row(4)] + [Row(8)] + [Row(10)] + [Row(50)] + [Row(100)] + [MultipleAsserts] + public void CanSolveForRandomMatrix(int order) + { + var matrixA = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixACopy = matrixA.Clone(); + var factorQR = matrixA.QR(); + + var matrixB = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixX = factorQR.Solve(matrixB); + + // The solution X row dimension is equal to the column dimension of A + Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); + // The solution X has the same number of columns as B + Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); + + var matrixBReconstruct = matrixA * matrixX; + + // Check the reconstruction. + for (var i = 0; i < matrixB.RowCount; i++) + { + for (var j = 0; j < matrixB.ColumnCount; j++) + { + Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-11); + } + } + + // Make sure A didn't change. + for (var i = 0; i < matrixA.RowCount; i++) + { + for (var j = 0; j < matrixA.ColumnCount; j++) + { + Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); + } + } + } + + [Test] + [Row(1)] + [Row(2)] + [Row(5)] + [Row(10)] + [Row(50)] + [Row(100)] + [MultipleAsserts] + public void CanSolveForRandomVectorWhenResultVectorGiven(int order) + { + var matrixA = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixACopy = matrixA.Clone(); + var factorQR = matrixA.QR(); + var vectorb = MatrixLoader.GenerateRandomVector(order); + var vectorbCopy = vectorb.Clone(); + var resultx = new DenseVector(order); + factorQR.Solve(vectorb,resultx); + + Assert.AreEqual(vectorb.Count, resultx.Count); + + var bReconstruct = matrixA * resultx; + + // Check the reconstruction. + for (var i = 0; i < vectorb.Count; i++) + { + Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1.0e-11); + } + + // Make sure A didn't change. + for (var i = 0; i < matrixA.RowCount; i++) + { + for (var j = 0; j < matrixA.ColumnCount; j++) + { + Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); + } + } + + // Make sure b didn't change. + for (var i = 0; i < vectorb.Count; i++) + { + Assert.AreEqual(vectorbCopy[i], vectorb[i]); + } + } + + [Test] + [Row(1)] + [Row(4)] + [Row(8)] + [Row(10)] + [Row(50)] + [Row(100)] + [MultipleAsserts] + public void CanSolveForRandomMatrixWhenResultMatrixGiven(int order) + { + var matrixA = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixACopy = matrixA.Clone(); + var factorQR = matrixA.QR(); + + var matrixB = MatrixLoader.GenerateRandomMatrix(order, order); + var matrixBCopy = matrixB.Clone(); + + var matrixX = new DenseMatrix(order, order); + factorQR.Solve(matrixB,matrixX); + + // The solution X row dimension is equal to the column dimension of A + Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); + // The solution X has the same number of columns as B + Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); + + var matrixBReconstruct = matrixA * matrixX; + + // Check the reconstruction. + for (var i = 0; i < matrixB.RowCount; i++) + { + for (var j = 0; j < matrixB.ColumnCount; j++) + { + Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-11); + } + } + + // Make sure A didn't change. + for (var i = 0; i < matrixA.RowCount; i++) + { + for (var j = 0; j < matrixA.ColumnCount; j++) + { + Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]); + } + } + + // Make sure B didn't change. + for (var i = 0; i < matrixB.RowCount; i++) + { + for (var j = 0; j < matrixB.ColumnCount; j++) + { + Assert.AreEqual(matrixBCopy[i, j], matrixB[i, j]); + } + } + } + } +} diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 0a1bafa6..67624be3 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -90,6 +90,7 @@ +