From 85db0a6884ee539653dc8196b2118c4e3b3da81f Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 26 Oct 2013 11:52:11 +0200 Subject: [PATCH] Managed Provider: cache-optimized matrix L-infinity norm (8-10 times faster) --- .../ManagedLinearAlgebraProvider.Complex.cs | 47 ++++++++------- .../ManagedLinearAlgebraProvider.Complex32.cs | 47 ++++++++------- .../ManagedLinearAlgebraProvider.Double.cs | 57 +++++++++++-------- .../ManagedLinearAlgebraProvider.Single.cs | 49 +++++++++------- 4 files changed, 114 insertions(+), 86 deletions(-) diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs index 1fb2f39d..db336da7 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs @@ -3,9 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// +// // Copyright (c) 2009-2013 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 +14,10 @@ // 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 @@ -200,7 +200,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra } /// - /// Does a point wise add of two arrays z = x + y. This can be used + /// Does a point wise add of two arrays z = x + y. This can be used /// to add vectors or matrices. /// /// The array x. @@ -251,7 +251,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra } /// - /// Does a point wise subtraction of two arrays z = x - y. This can be used + /// Does a point wise subtraction of two arrays z = x - y. This can be used /// to subtract vectors or matrices. /// /// The array x. @@ -431,26 +431,33 @@ namespace MathNet.Numerics.Providers.LinearAlgebra return norm1; case Norm.LargestAbsoluteValue: var normMax = 0d; - for (var i = 0; i < rows; i++) + for (var j = 0; j < columns; j++) { - for (var j = 0; j < columns; j++) + for (var i = 0; i < rows; i++) { normMax = Math.Max(matrix[(j * rows) + i].Magnitude, normMax); } } return normMax; case Norm.InfinityNorm: - var normInf = 0d; - for (var i = 0; i < rows; i++) + var r = new double[rows]; + for (var j = 0; j < columns; j++) { - var s = 0.0; - for (var j = 0; j < columns; j++) + for (var i = 0; i < rows; i++) { - s += matrix[(j*rows) + i].Magnitude; + r[i] += matrix[(j * rows) + i].Magnitude; + } + } + // TODO: reuse + var max = r[0]; + for (int i = 0; i < r.Length; i++) + { + if (r[i] > max) + { + max = r[i]; } - normInf = Math.Max(normInf, s); } - return normInf; + return max; case Norm.FrobeniusNorm: var aat = new Complex[rows*rows]; MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.ConjugateTranspose, 1.0, matrix, rows, columns, matrix, rows, columns, 0.0, aat); @@ -559,7 +566,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra { int m; // The number of rows of matrix op(A) and of the matrix C. int n; // The number of columns of matrix op(B) and of the matrix C. - int k; // The number of columns of matrix op(A) and the rows of the matrix op(B). + int k; // The number of columns of matrix op(A) and the rows of the matrix op(B). // First check some basic requirement on the parameters of the matrix multiplication. if (a == null) @@ -1519,7 +1526,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// 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 + /// On exit, A M by M matrix that holds the Q matrix of the /// QR factorization. /// A min(m,n) vector. On exit, contains additional information /// to be used by the QR solve routine. @@ -1702,7 +1709,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra ComputeQR(work, i, a, i, rowsA, i + 1, columnsA, Control.NumberOfParallelWorkerThreads); } - //copy R + //copy R for (var j = 0; j < columnsA; j++) { var rIndex = j*columnsA; @@ -2546,7 +2553,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra throw new NonConvergenceException(); } - // This section of the program inspects for negligible elements in the s and e arrays, + // 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 @@ -2810,7 +2817,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra if (computeVectors) { - // Finally transpose "v" to get "vt" matrix + // Finally transpose "v" to get "vt" matrix for (i = 0; i < columnsA; i++) { for (j = 0; j < columnsA; j++) diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs index c177e110..cb098939 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs @@ -3,9 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// +// // Copyright (c) 2009-2013 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 +14,10 @@ // 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 @@ -198,7 +198,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra } /// - /// Does a point wise add of two arrays z = x + y. This can be used + /// Does a point wise add of two arrays z = x + y. This can be used /// to add vectors or matrices. /// /// The array x. @@ -249,7 +249,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra } /// - /// Does a point wise subtraction of two arrays z = x - y. This can be used + /// Does a point wise subtraction of two arrays z = x - y. This can be used /// to subtract vectors or matrices. /// /// The array x. @@ -428,26 +428,33 @@ namespace MathNet.Numerics.Providers.LinearAlgebra return norm1; case Norm.LargestAbsoluteValue: var normMax = 0d; - for (var i = 0; i < rows; i++) + for (var j = 0; j < columns; j++) { - for (var j = 0; j < columns; j++) + for (var i = 0; i < rows; i++) { normMax = Math.Max(matrix[(j * rows) + i].Magnitude, normMax); } } return normMax; case Norm.InfinityNorm: - var normInf = 0d; - for (var i = 0; i < rows; i++) + var r = new double[rows]; + for (var j = 0; j < columns; j++) { - var s = 0d; - for (var j = 0; j < columns; j++) + for (var i = 0; i < rows; i++) { - s += matrix[(j*rows) + i].Magnitude; + r[i] += matrix[(j * rows) + i].Magnitude; + } + } + // TODO: reuse + var max = r[0]; + for (int i = 0; i < r.Length; i++) + { + if (r[i] > max) + { + max = r[i]; } - normInf = Math.Max(normInf, s); } - return normInf; + return max; case Norm.FrobeniusNorm: var aat = new Complex32[rows*rows]; MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.ConjugateTranspose, 1.0f, matrix, rows, columns, matrix, rows, columns, 0.0f, aat); @@ -556,7 +563,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra { int m; // The number of rows of matrix op(A) and of the matrix C. int n; // The number of columns of matrix op(B) and of the matrix C. - int k; // The number of columns of matrix op(A) and the rows of the matrix op(B). + int k; // The number of columns of matrix op(A) and the rows of the matrix op(B). // First check some basic requirement on the parameters of the matrix multiplication. if (a == null) @@ -1516,7 +1523,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// 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 + /// On exit, A M by M matrix that holds the Q matrix of the /// QR factorization. /// A min(m,n) vector. On exit, contains additional information /// to be used by the QR solve routine. @@ -1699,7 +1706,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra ComputeQR(work, i, a, i, rowsA, i + 1, columnsA, Control.NumberOfParallelWorkerThreads); } - //copy R + //copy R for (var j = 0; j < columnsA; j++) { var rIndex = j*columnsA; @@ -2543,7 +2550,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra throw new NonConvergenceException(); } - // This section of the program inspects for negligible elements in the s and e arrays, + // 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 @@ -2807,7 +2814,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra if (computeVectors) { - // Finally transpose "v" to get "vt" matrix + // Finally transpose "v" to get "vt" matrix for (i = 0; i < columnsA; i++) { for (j = 0; j < columnsA; j++) diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs index 9271512f..3194284a 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs @@ -3,9 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// +// // Copyright (c) 2009-2013 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 +14,10 @@ // 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 @@ -195,7 +195,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra } /// - /// Does a point wise add of two arrays z = x + y. This can be used + /// Does a point wise add of two arrays z = x + y. This can be used /// to add vectors or matrices. /// /// The array x. @@ -246,7 +246,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra } /// - /// Does a point wise subtraction of two arrays z = x - y. This can be used + /// Does a point wise subtraction of two arrays z = x - y. This can be used /// to subtract vectors or matrices. /// /// The array x. @@ -426,26 +426,33 @@ namespace MathNet.Numerics.Providers.LinearAlgebra return norm1; case Norm.LargestAbsoluteValue: var normMax = 0d; - for (var i = 0; i < rows; i++) + for (var j = 0; j < columns; j++) { - for (var j = 0; j < columns; j++) + for (var i = 0; i < rows; i++) { normMax = Math.Max(Math.Abs(matrix[(j * rows) + i]), normMax); } } return normMax; case Norm.InfinityNorm: - var normInf = 0d; - for (var i = 0; i < rows; i++) - { - var s = 0.0; - for (var j = 0; j < columns; j++) - { - s += Math.Abs(matrix[(j*rows) + i]); - } - normInf = Math.Max(normInf, s); - } - return normInf; + var r = new double[rows]; + for (var j = 0; j < columns; j++) + { + for (var i = 0; i < rows; i++) + { + r[i] += Math.Abs(matrix[(j * rows) + i]); + } + } + // TODO: reuse + var max = r[0]; + for (int i = 0; i < r.Length; i++) + { + if (r[i] > max) + { + max = r[i]; + } + } + return max; case Norm.FrobeniusNorm: var aat = new double[rows*rows]; MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.Transpose, 1.0, matrix, rows, columns, matrix, rows, columns, 0.0, aat); @@ -554,7 +561,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra { int m; // The number of rows of matrix op(A) and of the matrix C. int n; // The number of columns of matrix op(B) and of the matrix C. - int k; // The number of columns of matrix op(A) and the rows of the matrix op(B). + int k; // The number of columns of matrix op(A) and the rows of the matrix op(B). // First check some basic requirement on the parameters of the matrix multiplication. if (a == null) @@ -1406,7 +1413,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// 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 + /// On exit, A M by M matrix that holds the Q matrix of the /// QR factorization. /// A min(m,n) vector. On exit, contains additional information /// to be used by the QR solve routine. @@ -1589,7 +1596,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra ComputeQR(work, i, a, i, rowsA, i + 1, columnsA, Control.NumberOfParallelWorkerThreads); } - //copy R + //copy R for (var j = 0; j < columnsA; j++) { var rIndex = j*columnsA; @@ -2435,7 +2442,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra throw new NonConvergenceException(); } - // This section of the program inspects for negligible elements in the s and e arrays, + // 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 @@ -2701,7 +2708,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra if (computeVectors) { - // Finally transpose "v" to get "vt" matrix + // Finally transpose "v" to get "vt" matrix for (i = 0; i < columnsA; i++) { for (j = 0; j < columnsA; j++) @@ -2722,7 +2729,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra } /// - /// Given the Cartesian coordinates (da, db) of a point p, these function return the parameters da, db, c, and s + /// Given the Cartesian coordinates (da, db) of a point p, these function return the parameters da, db, c, and s /// associated with the Givens rotation that zeros the y-coordinate of the point. /// /// Provides the x-coordinate of the point p. On exit contains the parameter r associated with the Givens rotation diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs index e7f5de64..d472a120 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs @@ -3,9 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// +// // Copyright (c) 2009-2013 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 +14,10 @@ // 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 @@ -195,7 +195,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra } /// - /// Does a point wise add of two arrays z = x + y. This can be used + /// Does a point wise add of two arrays z = x + y. This can be used /// to add vectors or matrices. /// /// The array x. @@ -246,7 +246,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra } /// - /// Does a point wise subtraction of two arrays z = x - y. This can be used + /// Does a point wise subtraction of two arrays z = x - y. This can be used /// to subtract vectors or matrices. /// /// The array x. @@ -426,26 +426,33 @@ namespace MathNet.Numerics.Providers.LinearAlgebra return norm1; case Norm.LargestAbsoluteValue: var normMax = 0d; - for (var i = 0; i < rows; i++) + for (var j = 0; j < columns; j++) { - for (var j = 0; j < columns; j++) + for (var i = 0; i < rows; i++) { normMax = Math.Max(Math.Abs(matrix[(j * rows) + i]), normMax); } } return normMax; case Norm.InfinityNorm: - var normInf = 0d; - for (var i = 0; i < rows; i++) + var r = new double[rows]; + for (var j = 0; j < columns; j++) { - var s = 0d; - for (var j = 0; j < columns; j++) + for (var i = 0; i < rows; i++) { - s += Math.Abs(matrix[(j*rows) + i]); + r[i] += Math.Abs(matrix[(j * rows) + i]); + } + } + // TODO: reuse + var max = r[0]; + for (int i = 0; i < r.Length; i++) + { + if (r[i] > max) + { + max = r[i]; } - normInf = Math.Max(normInf, s); } - return normInf; + return max; case Norm.FrobeniusNorm: var aat = new float[rows*rows]; MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.Transpose, 1.0f, matrix, rows, columns, matrix, rows, columns, 0.0f, aat); @@ -554,7 +561,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra { int m; // The number of rows of matrix op(A) and of the matrix C. int n; // The number of columns of matrix op(B) and of the matrix C. - int k; // The number of columns of matrix op(A) and the rows of the matrix op(B). + int k; // The number of columns of matrix op(A) and the rows of the matrix op(B). // First check some basic requirement on the parameters of the matrix multiplication. if (a == null) @@ -1405,7 +1412,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra /// 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 + /// On exit, A M by M matrix that holds the Q matrix of the /// QR factorization. /// A min(m,n) vector. On exit, contains additional information /// to be used by the QR solve routine. @@ -1588,7 +1595,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra ComputeQR(work, i, a, i, rowsA, i + 1, columnsA, Control.NumberOfParallelWorkerThreads); } - //copy R + //copy R for (var j = 0; j < columnsA; j++) { var rIndex = j*columnsA; @@ -2436,7 +2443,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra throw new NonConvergenceException(); } - // This section of the program inspects for negligible elements in the s and e arrays, + // 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 @@ -2702,7 +2709,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra if (computeVectors) { - // Finally transpose "v" to get "vt" matrix + // Finally transpose "v" to get "vt" matrix for (i = 0; i < columnsA; i++) { for (j = 0; j < columnsA; j++) @@ -2723,7 +2730,7 @@ namespace MathNet.Numerics.Providers.LinearAlgebra } /// - /// Given the Cartesian coordinates (da, db) of a point p, these function return the parameters da, db, c, and s + /// Given the Cartesian coordinates (da, db) of a point p, these function return the parameters da, db, c, and s /// associated with the Givens rotation that zeros the y-coordinate of the point. /// /// Provides the x-coordinate of the point p. On exit contains the parameter r associated with the Givens rotation