diff --git a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs index 4e3b52eb..49ef76c5 100644 --- a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs @@ -1011,7 +1011,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex public override Cholesky Cholesky() { - return new DenseCholesky(this); + return DenseCholesky.Create(this); } public override LU LU() diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/Cholesky.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/Cholesky.cs index 17067796..738b511a 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/Cholesky.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/Cholesky.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// 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 @@ -32,10 +32,11 @@ using MathNet.Numerics.LinearAlgebra.Factorization; namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization { - #if NOSYSNUMERICS +#if NOSYSNUMERICS using Complex = Numerics.Complex; #else using Complex = System.Numerics.Complex; + #endif /// @@ -49,6 +50,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// public abstract class Cholesky : Cholesky { + protected Cholesky(Matrix factor) + : base(factor) + { + } + /// /// Gets the determinant of the matrix for which the Cholesky matrix was computed. /// @@ -60,7 +66,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization for (var j = 0; j < Factor.RowCount; j++) { var d = Factor.At(j, j); - det *= d * d; + det *= d*d; } return det; @@ -77,7 +83,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var det = Complex.Zero; for (var j = 0; j < Factor.RowCount; j++) { - det += 2.0 * Factor.At(j, j).Ln(); + det += 2.0*Factor.At(j, j).Ln(); } return det; diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs index 750cff73..035d0af5 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// 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 @@ -38,6 +38,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization using Numerics; #else using System.Numerics; + #endif /// @@ -49,7 +50,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.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 class DenseCholesky : Cholesky + public sealed class DenseCholesky : Cholesky { /// /// Initializes a new instance of the class. This object will compute the @@ -59,22 +60,22 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// If is null. /// If is not a square matrix. /// If is not positive definite. - public DenseCholesky(DenseMatrix matrix) + public static DenseCholesky Create(DenseMatrix matrix) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare); } // Create a new matrix for the Cholesky factor, then perform factorization (while overwriting). - var factor = (DenseMatrix)matrix.Clone(); + var factor = (DenseMatrix) matrix.Clone(); Control.LinearAlgebraProvider.CholeskyFactor(factor.Values, factor.RowCount); - Factor = factor; + return new DenseCholesky(factor); + } + + DenseCholesky(Matrix factor) + : base(factor) + { } /// @@ -84,18 +85,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// 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"); - } - - // Check for proper dimensions. if (result.RowCount != input.RowCount) { throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); @@ -127,7 +116,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization Array.Copy(dinput.Values, dresult.Values, dinput.Values.Length); // Cholesky solve by overwriting result. - var dfactor = (DenseMatrix)Factor; + var dfactor = (DenseMatrix) Factor; Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Values, dfactor.RowCount, dresult.Values, dresult.ColumnCount); } @@ -138,18 +127,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// The left hand side , x. public override void Solve(Vector input, Vector result) { - // Check for proper arguments. - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - - // Check for proper dimensions. if (input.Count != result.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); @@ -176,7 +153,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization Array.Copy(dinput.Values, dresult.Values, dinput.Values.Length); // Cholesky solve by overwriting result. - var dfactor = (DenseMatrix)Factor; + var dfactor = (DenseMatrix) Factor; Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Values, dfactor.RowCount, dresult.Values, 1); } } diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/UserCholesky.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/UserCholesky.cs index 5d01693a..88a9b7a3 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/UserCholesky.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/UserCholesky.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// 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 @@ -39,6 +39,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization using Numerics; #else using System.Numerics; + #endif /// @@ -50,7 +51,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.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 class UserCholesky : Cholesky + public sealed class UserCholesky : Cholesky { /// /// Initializes a new instance of the class. This object will compute the @@ -60,55 +61,57 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// If is null. /// If is not a square matrix. /// If is not positive definite. - public UserCholesky(Matrix matrix) + public static UserCholesky Create(Matrix matrix) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare); } // Create a new matrix for the Cholesky factor, then perform factorization (while overwriting). - Factor = matrix.Clone(); - var tmpColumn = new Complex[Factor.RowCount]; + var factor = matrix.Clone(); + var tmpColumn = new Complex[factor.RowCount]; // Main loop - along the diagonal - for (var ij = 0; ij < Factor.RowCount; ij++) + for (var ij = 0; ij < factor.RowCount; ij++) { // "Pivot" element - var tmpVal = Factor.At(ij, ij); + var tmpVal = factor.At(ij, ij); if (tmpVal.Real > 0.0) { tmpVal = tmpVal.SquareRoot(); - Factor.At(ij, ij, tmpVal); + factor.At(ij, ij, tmpVal); tmpColumn[ij] = tmpVal; // Calculate multipliers and copy to local column // Current column, below the diagonal - for (var i = ij + 1; i < Factor.RowCount; i++) + for (var i = ij + 1; i < factor.RowCount; i++) { - Factor.At(i, ij, Factor.At(i, ij) / tmpVal); - tmpColumn[i] = Factor.At(i, ij); + factor.At(i, ij, factor.At(i, ij)/tmpVal); + tmpColumn[i] = factor.At(i, ij); } // Remaining columns, below the diagonal - DoCholeskyStep(Factor, Factor.RowCount, ij + 1, Factor.RowCount, tmpColumn, Control.NumberOfParallelWorkerThreads); + DoCholeskyStep(factor, factor.RowCount, ij + 1, factor.RowCount, tmpColumn, Control.NumberOfParallelWorkerThreads); } else { throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite); } - for (var i = ij + 1; i < Factor.RowCount; i++) + for (var i = ij + 1; i < factor.RowCount; i++) { - Factor.At(ij, i, Complex.Zero); + factor.At(ij, i, Complex.Zero); } } + + return new UserCholesky(factor); + } + + UserCholesky(Matrix factor) + : base(factor) + { } /// @@ -120,14 +123,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// Total columns /// Multipliers calculated previously /// Number of available processors - private static void DoCholeskyStep(Matrix data, int rowDim, int firstCol, int colLimit, Complex[] multipliers, int availableCores) + static void DoCholeskyStep(Matrix data, int rowDim, int firstCol, int colLimit, Complex[] multipliers, int availableCores) { var tmpColCount = colLimit - firstCol; if ((availableCores > 1) && (tmpColCount > 200)) { - var tmpSplit = firstCol + (tmpColCount / 3); - var tmpCores = availableCores / 2; + var tmpSplit = firstCol + (tmpColCount/3); + var tmpCores = availableCores/2; CommonParallel.Invoke( () => DoCholeskyStep(data, rowDim, firstCol, tmpSplit, multipliers, tmpCores), @@ -140,7 +143,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization var tmpVal = multipliers[j]; for (var i = j; i < rowDim; i++) { - data.At(i, j, data.At(i, j) - (multipliers[i] * tmpVal.Conjugate())); + data.At(i, j, data.At(i, j) - (multipliers[i]*tmpVal.Conjugate())); } } } @@ -153,17 +156,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// The left hand side , X. public override void Solve(Matrix input, Matrix result) { - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - - // Check for proper dimensions. if (result.RowCount != input.RowCount) { throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); @@ -191,10 +183,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization sum = result.At(i, c); for (var k = i - 1; k >= 0; k--) { - sum -= Factor.At(i, k) * result.At(k, c); + sum -= Factor.At(i, k)*result.At(k, c); } - result.At(i, c, sum / Factor.At(i, i)); + result.At(i, c, sum/Factor.At(i, i)); } // Solve L'*X = Y; @@ -203,10 +195,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization sum = result.At(i, c); for (var k = i + 1; k < order; k++) { - sum -= Factor.At(k, i).Conjugate() * result.At(k, c); + sum -= Factor.At(k, i).Conjugate()*result.At(k, c); } - result.At(i, c, sum / Factor.At(i, i)); + result.At(i, c, sum/Factor.At(i, i)); } } } @@ -218,18 +210,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization /// The left hand side , x. public override void Solve(Vector input, Vector result) { - // Check for proper arguments. - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - - // Check for proper dimensions. if (input.Count != result.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); @@ -250,10 +230,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization sum = result[i]; for (var k = i - 1; k >= 0; k--) { - sum -= Factor.At(i, k) * result[k]; + sum -= Factor.At(i, k)*result[k]; } - result[i] = sum / Factor.At(i, i); + result[i] = sum/Factor.At(i, i); } // Solve L'*X = Y; @@ -262,10 +242,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization sum = result[i]; for (var k = i + 1; k < order; k++) { - sum -= Factor.At(k, i).Conjugate() * result[k]; + sum -= Factor.At(k, i).Conjugate()*result[k]; } - result[i] = sum / Factor.At(i, i); + result[i] = sum/Factor.At(i, i); } } } diff --git a/src/Numerics/LinearAlgebra/Complex/Matrix.cs b/src/Numerics/LinearAlgebra/Complex/Matrix.cs index 6690c7ed..4095c001 100644 --- a/src/Numerics/LinearAlgebra/Complex/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Complex/Matrix.cs @@ -460,7 +460,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex public override Cholesky Cholesky() { - return new UserCholesky(this); + return UserCholesky.Create(this); } public override LU LU() diff --git a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs index 79dcb352..12a7c0bb 100644 --- a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs @@ -1006,7 +1006,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 public override Cholesky Cholesky() { - return new DenseCholesky(this); + return DenseCholesky.Create(this); } public override LU LU() diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/Cholesky.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/Cholesky.cs index 49374cb7..a85a2c26 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/Cholesky.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/Cholesky.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// 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 @@ -45,6 +45,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// public abstract class Cholesky : Cholesky { + protected Cholesky(Matrix factor) + : base(factor) + { + } + /// /// Gets the determinant of the matrix for which the Cholesky matrix was computed. /// @@ -56,7 +61,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization for (var j = 0; j < Factor.RowCount; j++) { var d = Factor.At(j, j); - det *= d * d; + det *= d*d; } return det; @@ -73,7 +78,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var det = Complex32.Zero; for (var j = 0; j < Factor.RowCount; j++) { - det += 2.0f * Factor.At(j, j).NaturalLogarithm(); + det += 2.0f*Factor.At(j, j).NaturalLogarithm(); } return det; diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseCholesky.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseCholesky.cs index 52f9a8cb..28f928fb 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseCholesky.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseCholesky.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// 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 @@ -28,11 +28,12 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System; +using MathNet.Numerics.Properties; + namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { - using System; using Numerics; - using Properties; /// /// A class which encapsulates the functionality of a Cholesky factorization for dense matrices. @@ -43,7 +44,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.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 class DenseCholesky : Cholesky + public sealed class DenseCholesky : Cholesky { /// /// Initializes a new instance of the class. This object will compute the @@ -53,22 +54,22 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// If is null. /// If is not a square matrix. /// If is not positive definite. - public DenseCholesky(DenseMatrix matrix) + public static DenseCholesky Create(DenseMatrix matrix) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare); } // Create a new matrix for the Cholesky factor, then perform factorization (while overwriting). - var factor = (DenseMatrix)matrix.Clone(); + var factor = (DenseMatrix) matrix.Clone(); Control.LinearAlgebraProvider.CholeskyFactor(factor.Values, factor.RowCount); - Factor = factor; + return new DenseCholesky(factor); + } + + DenseCholesky(Matrix factor) + : base(factor) + { } /// @@ -78,18 +79,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// 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"); - } - - // Check for proper dimensions. if (result.RowCount != input.RowCount) { throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); @@ -121,7 +110,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization Array.Copy(dinput.Values, dresult.Values, dinput.Values.Length); // Cholesky solve by overwriting result. - var dfactor = (DenseMatrix)Factor; + var dfactor = (DenseMatrix) Factor; Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Values, dfactor.RowCount, dresult.Values, dresult.ColumnCount); } @@ -132,18 +121,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// The left hand side , x. public override void Solve(Vector input, Vector result) { - // Check for proper arguments. - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - - // Check for proper dimensions. if (input.Count != result.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); @@ -170,7 +147,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization Array.Copy(dinput.Values, dresult.Values, dinput.Values.Length); // Cholesky solve by overwriting result. - var dfactor = (DenseMatrix)Factor; + var dfactor = (DenseMatrix) Factor; Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Values, dfactor.RowCount, dresult.Values, 1); } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserCholesky.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserCholesky.cs index b6b2a229..246844ed 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserCholesky.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserCholesky.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// 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 @@ -28,9 +28,9 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System; using MathNet.Numerics.Properties; using MathNet.Numerics.Threading; -using System; namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { @@ -45,7 +45,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.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 class UserCholesky : Cholesky + public sealed class UserCholesky : Cholesky { /// /// Initializes a new instance of the class. This object will compute the @@ -55,55 +55,57 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// If is null. /// If is not a square matrix. /// If is not positive definite. - public UserCholesky(Matrix matrix) + public static UserCholesky Create(Matrix matrix) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare); } // Create a new matrix for the Cholesky factor, then perform factorization (while overwriting). - Factor = matrix.Clone(); - var tmpColumn = new Complex32[Factor.RowCount]; + var factor = matrix.Clone(); + var tmpColumn = new Complex32[factor.RowCount]; // Main loop - along the diagonal - for (var ij = 0; ij < Factor.RowCount; ij++) + for (var ij = 0; ij < factor.RowCount; ij++) { // "Pivot" element - var tmpVal = Factor.At(ij, ij); + var tmpVal = factor.At(ij, ij); if (tmpVal.Real > 0.0) { tmpVal = tmpVal.SquareRoot(); - Factor.At(ij, ij, tmpVal); + factor.At(ij, ij, tmpVal); tmpColumn[ij] = tmpVal; // Calculate multipliers and copy to local column // Current column, below the diagonal - for (var i = ij + 1; i < Factor.RowCount; i++) + for (var i = ij + 1; i < factor.RowCount; i++) { - Factor.At(i, ij, Factor.At(i, ij) / tmpVal); - tmpColumn[i] = Factor.At(i, ij); + factor.At(i, ij, factor.At(i, ij)/tmpVal); + tmpColumn[i] = factor.At(i, ij); } // Remaining columns, below the diagonal - DoCholeskyStep(Factor, Factor.RowCount, ij + 1, Factor.RowCount, tmpColumn, Control.NumberOfParallelWorkerThreads); + DoCholeskyStep(factor, factor.RowCount, ij + 1, factor.RowCount, tmpColumn, Control.NumberOfParallelWorkerThreads); } else { throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite); } - for (var i = ij + 1; i < Factor.RowCount; i++) + for (var i = ij + 1; i < factor.RowCount; i++) { - Factor.At(ij, i, Complex32.Zero); + factor.At(ij, i, Complex32.Zero); } } + + return new UserCholesky(factor); + } + + UserCholesky(Matrix factor) + : base(factor) + { } /// @@ -115,14 +117,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// Total columns /// Multipliers calculated previously /// Number of available processors - private static void DoCholeskyStep(Matrix data, int rowDim, int firstCol, int colLimit, Complex32[] multipliers, int availableCores) + static void DoCholeskyStep(Matrix data, int rowDim, int firstCol, int colLimit, Complex32[] multipliers, int availableCores) { var tmpColCount = colLimit - firstCol; if ((availableCores > 1) && (tmpColCount > 200)) { - var tmpSplit = firstCol + (tmpColCount / 3); - var tmpCores = availableCores / 2; + var tmpSplit = firstCol + (tmpColCount/3); + var tmpCores = availableCores/2; CommonParallel.Invoke( () => DoCholeskyStep(data, rowDim, firstCol, tmpSplit, multipliers, tmpCores), @@ -135,7 +137,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization var tmpVal = multipliers[j]; for (var i = j; i < rowDim; i++) { - data.At(i, j, data.At(i, j) - (multipliers[i] * tmpVal.Conjugate())); + data.At(i, j, data.At(i, j) - (multipliers[i]*tmpVal.Conjugate())); } } } @@ -148,17 +150,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// The left hand side , X. public override void Solve(Matrix input, Matrix result) { - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - - // Check for proper dimensions. if (result.RowCount != input.RowCount) { throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); @@ -186,10 +177,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization sum = result.At(i, c); for (var k = i - 1; k >= 0; k--) { - sum -= Factor.At(i, k) * result.At(k, c); + sum -= Factor.At(i, k)*result.At(k, c); } - result.At(i, c, sum / Factor.At(i, i)); + result.At(i, c, sum/Factor.At(i, i)); } // Solve L'*X = Y; @@ -198,10 +189,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization sum = result.At(i, c); for (var k = i + 1; k < order; k++) { - sum -= Factor.At(k, i).Conjugate() * result.At(k, c); + sum -= Factor.At(k, i).Conjugate()*result.At(k, c); } - result.At(i, c, sum / Factor.At(i, i)); + result.At(i, c, sum/Factor.At(i, i)); } } } @@ -213,18 +204,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization /// The left hand side , x. public override void Solve(Vector input, Vector result) { - // Check for proper arguments. - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - - // Check for proper dimensions. if (input.Count != result.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); @@ -245,10 +224,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization sum = result[i]; for (var k = i - 1; k >= 0; k--) { - sum -= Factor.At(i, k) * result[k]; + sum -= Factor.At(i, k)*result[k]; } - result[i] = sum / Factor.At(i, i); + result[i] = sum/Factor.At(i, i); } // Solve L'*X = Y; @@ -257,10 +236,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization sum = result[i]; for (var k = i + 1; k < order; k++) { - sum -= Factor.At(k, i).Conjugate() * result[k]; + sum -= Factor.At(k, i).Conjugate()*result[k]; } - result[i] = sum / Factor.At(i, i); + result[i] = sum/Factor.At(i, i); } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Matrix.cs b/src/Numerics/LinearAlgebra/Complex32/Matrix.cs index c4123ee1..0028ac73 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Matrix.cs @@ -455,7 +455,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 public override Cholesky Cholesky() { - return new UserCholesky(this); + return UserCholesky.Create(this); } public override LU LU() diff --git a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs index 635226d8..a628faaa 100644 --- a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs @@ -1037,7 +1037,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double public override Cholesky Cholesky() { - return new DenseCholesky(this); + return DenseCholesky.Create(this); } public override LU LU() diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs b/src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs index 1bb8eef0..63980447 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// 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 @@ -45,6 +45,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// public abstract class Cholesky : Cholesky { + protected Cholesky(Matrix factor) + : base(factor) + { + } + /// /// Gets the determinant of the matrix for which the Cholesky matrix was computed. /// @@ -56,7 +61,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization for (var j = 0; j < Factor.RowCount; j++) { var d = Factor.At(j, j); - det *= d * d; + det *= d*d; } return det; @@ -73,7 +78,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var det = 0.0; for (var j = 0; j < Factor.RowCount; j++) { - det += 2 * Math.Log(Factor.At(j, j)); + det += 2*Math.Log(Factor.At(j, j)); } return det; diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs index 06d587d6..779f5030 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// 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 @@ -28,8 +28,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using MathNet.Numerics.Properties; using System; +using MathNet.Numerics.Properties; namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { @@ -42,7 +42,7 @@ 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 class DenseCholesky : Cholesky + public sealed class DenseCholesky : Cholesky { /// /// Initializes a new instance of the class. This object will compute the @@ -52,22 +52,22 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// If is null. /// If is not a square matrix. /// If is not positive definite. - public DenseCholesky(DenseMatrix matrix) + public static DenseCholesky Create(DenseMatrix matrix) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare); } // Create a new matrix for the Cholesky factor, then perform factorization (while overwriting). - var factor = (DenseMatrix)matrix.Clone(); + var factor = (DenseMatrix) matrix.Clone(); Control.LinearAlgebraProvider.CholeskyFactor(factor.Values, factor.RowCount); - Factor = factor; + return new DenseCholesky(factor); + } + + DenseCholesky(Matrix factor) + : base(factor) + { } /// @@ -77,18 +77,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// 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"); - } - - // Check for proper dimensions. if (result.RowCount != input.RowCount) { throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); @@ -117,10 +105,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization } // Copy the contents of input to result. - Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length * Constants.SizeOfDouble); + Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length*Constants.SizeOfDouble); // Cholesky solve by overwriting result. - var dfactor = (DenseMatrix)Factor; + var dfactor = (DenseMatrix) Factor; Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Values, dfactor.RowCount, dresult.Values, dresult.ColumnCount); } @@ -131,18 +119,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// The left hand side , x. public override void Solve(Vector input, Vector result) { - // Check for proper arguments. - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - - // Check for proper dimensions. if (input.Count != result.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); @@ -166,10 +142,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization } // Copy the contents of input to result. - Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length * Constants.SizeOfDouble); + Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length*Constants.SizeOfDouble); // Cholesky solve by overwriting result. - var dfactor = (DenseMatrix)Factor; + var dfactor = (DenseMatrix) Factor; Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Values, dfactor.RowCount, dresult.Values, 1); } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/UserCholesky.cs b/src/Numerics/LinearAlgebra/Double/Factorization/UserCholesky.cs index 939c2b9a..651a0b45 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/UserCholesky.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/UserCholesky.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// 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 @@ -28,9 +28,9 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System; using MathNet.Numerics.Properties; using MathNet.Numerics.Threading; -using System; namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { @@ -43,7 +43,7 @@ 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 class UserCholesky : Cholesky + public sealed class UserCholesky : Cholesky { /// /// Initializes a new instance of the class. This object will compute the @@ -53,55 +53,57 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// If is null. /// If is not a square matrix. /// If is not positive definite. - public UserCholesky(Matrix matrix) + public static UserCholesky Create(Matrix matrix) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare); } // Create a new matrix for the Cholesky factor, then perform factorization (while overwriting). - Factor = matrix.Clone(); - var tmpColumn = new double[Factor.RowCount]; + var factor = matrix.Clone(); + var tmpColumn = new double[factor.RowCount]; // Main loop - along the diagonal - for (var ij = 0; ij < Factor.RowCount; ij++) + for (var ij = 0; ij < factor.RowCount; ij++) { // "Pivot" element - var tmpVal = Factor.At(ij, ij); + var tmpVal = factor.At(ij, ij); if (tmpVal > 0.0) { tmpVal = Math.Sqrt(tmpVal); - Factor.At(ij, ij, tmpVal); + factor.At(ij, ij, tmpVal); tmpColumn[ij] = tmpVal; // Calculate multipliers and copy to local column // Current column, below the diagonal - for (var i = ij + 1; i < Factor.RowCount; i++) + for (var i = ij + 1; i < factor.RowCount; i++) { - Factor.At(i, ij, Factor.At(i, ij) / tmpVal); - tmpColumn[i] = Factor.At(i, ij); + factor.At(i, ij, factor.At(i, ij)/tmpVal); + tmpColumn[i] = factor.At(i, ij); } // Remaining columns, below the diagonal - DoCholeskyStep(Factor, Factor.RowCount, ij + 1, Factor.RowCount, tmpColumn, Control.NumberOfParallelWorkerThreads); + DoCholeskyStep(factor, factor.RowCount, ij + 1, factor.RowCount, tmpColumn, Control.NumberOfParallelWorkerThreads); } else { throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite); } - for (var i = ij + 1; i < Factor.RowCount; i++) + for (var i = ij + 1; i < factor.RowCount; i++) { - Factor.At(ij, i, 0.0); + factor.At(ij, i, 0.0); } } + + return new UserCholesky(factor); + } + + UserCholesky(Matrix factor) + : base(factor) + { } /// @@ -113,14 +115,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// Total columns /// Multipliers calculated previously /// Number of available processors - private static void DoCholeskyStep(Matrix data, int rowDim, int firstCol, int colLimit, double[] multipliers, int availableCores) + static void DoCholeskyStep(Matrix data, int rowDim, int firstCol, int colLimit, double[] multipliers, int availableCores) { var tmpColCount = colLimit - firstCol; if ((availableCores > 1) && (tmpColCount > 200)) { - var tmpSplit = firstCol + (tmpColCount / 3); - var tmpCores = availableCores / 2; + var tmpSplit = firstCol + (tmpColCount/3); + var tmpCores = availableCores/2; CommonParallel.Invoke( () => DoCholeskyStep(data, rowDim, firstCol, tmpSplit, multipliers, tmpCores), @@ -133,7 +135,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var tmpVal = multipliers[j]; for (var i = j; i < rowDim; i++) { - data.At(i, j, data.At(i, j) - (multipliers[i] * tmpVal)); + data.At(i, j, data.At(i, j) - (multipliers[i]*tmpVal)); } } } @@ -146,17 +148,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// The left hand side , X. public override void Solve(Matrix input, Matrix result) { - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - - // Check for proper dimensions. if (result.RowCount != input.RowCount) { throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); @@ -184,10 +175,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization sum = result.At(i, c); for (var k = i - 1; k >= 0; k--) { - sum -= Factor.At(i, k) * result.At(k, c); + sum -= Factor.At(i, k)*result.At(k, c); } - result.At(i, c, sum / Factor.At(i, i)); + result.At(i, c, sum/Factor.At(i, i)); } // Solve L'*X = Y; @@ -196,10 +187,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization sum = result.At(i, c); for (var k = i + 1; k < order; k++) { - sum -= Factor.At(k, i) * result.At(k, c); + sum -= Factor.At(k, i)*result.At(k, c); } - result.At(i, c, sum / Factor.At(i, i)); + result.At(i, c, sum/Factor.At(i, i)); } } } @@ -211,18 +202,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// The left hand side , x. public override void Solve(Vector input, Vector result) { - // Check for proper arguments. - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - - // Check for proper dimensions. if (input.Count != result.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); @@ -243,10 +222,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization sum = result[i]; for (var k = i - 1; k >= 0; k--) { - sum -= Factor.At(i, k) * result[k]; + sum -= Factor.At(i, k)*result[k]; } - result[i] = sum / Factor.At(i, i); + result[i] = sum/Factor.At(i, i); } // Solve L'*X = Y; @@ -255,10 +234,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization sum = result[i]; for (var k = i + 1; k < order; k++) { - sum -= Factor.At(k, i) * result[k]; + sum -= Factor.At(k, i)*result[k]; } - result[i] = sum / Factor.At(i, i); + result[i] = sum/Factor.At(i, i); } } } diff --git a/src/Numerics/LinearAlgebra/Double/Matrix.cs b/src/Numerics/LinearAlgebra/Double/Matrix.cs index 7440e30c..6d40080c 100644 --- a/src/Numerics/LinearAlgebra/Double/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Double/Matrix.cs @@ -461,7 +461,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double public override Cholesky Cholesky() { - return new UserCholesky(this); + return UserCholesky.Create(this); } public override LU LU() diff --git a/src/Numerics/LinearAlgebra/Factorization/Cholesky.cs b/src/Numerics/LinearAlgebra/Factorization/Cholesky.cs index f45c7fac..07939b4a 100644 --- a/src/Numerics/LinearAlgebra/Factorization/Cholesky.cs +++ b/src/Numerics/LinearAlgebra/Factorization/Cholesky.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// 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 @@ -43,12 +43,17 @@ namespace MathNet.Numerics.LinearAlgebra.Factorization /// /// Supported data types are double, single, , and . public abstract class Cholesky : ISolver - where T : struct, IEquatable, IFormattable + where T : struct, IEquatable, IFormattable { + protected Cholesky(Matrix factor) + { + Factor = factor; + } + /// /// Gets the lower triangular form of the Cholesky matrix. /// - public Matrix Factor { get; protected set; } + public Matrix Factor { get; private set; } /// /// Gets the determinant of the matrix for which the Cholesky matrix was computed. @@ -67,12 +72,6 @@ namespace MathNet.Numerics.LinearAlgebra.Factorization /// The left hand side , X. public virtual Matrix Solve(Matrix input) { - // Check for proper arguments. - if (input == null) - { - throw new ArgumentNullException("input"); - } - var x = input.CreateMatrix(input.RowCount, input.ColumnCount); Solve(input, x); return x; @@ -92,12 +91,6 @@ namespace MathNet.Numerics.LinearAlgebra.Factorization /// 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(input.Count); Solve(input, x); return x; diff --git a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs index c6033e84..b8555c65 100644 --- a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs @@ -1037,7 +1037,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single public override Cholesky Cholesky() { - return new DenseCholesky(this); + return DenseCholesky.Create(this); } public override LU LU() diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/Cholesky.cs b/src/Numerics/LinearAlgebra/Single/Factorization/Cholesky.cs index 312c66b5..00f23409 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/Cholesky.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/Cholesky.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// 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 @@ -45,6 +45,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// public abstract class Cholesky : Cholesky { + protected Cholesky(Matrix factor) + : base(factor) + { + } + /// /// Gets the determinant of the matrix for which the Cholesky matrix was computed. /// @@ -56,7 +61,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization for (var j = 0; j < Factor.RowCount; j++) { var d = Factor.At(j, j); - det *= d * d; + det *= d*d; } return det; @@ -73,7 +78,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var det = 0.0f; for (var j = 0; j < Factor.RowCount; j++) { - det += 2.0f * Convert.ToSingle(Math.Log(Factor.At(j, j))); + det += 2.0f*Convert.ToSingle(Math.Log(Factor.At(j, j))); } return det; diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseCholesky.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseCholesky.cs index 15d9ebce..437e613c 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseCholesky.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/DenseCholesky.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// 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 @@ -28,8 +28,8 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using MathNet.Numerics.Properties; using System; +using MathNet.Numerics.Properties; namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { @@ -42,7 +42,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.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 class DenseCholesky : Cholesky + public sealed class DenseCholesky : Cholesky { /// /// Initializes a new instance of the class. This object will compute the @@ -52,22 +52,22 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// If is null. /// If is not a square matrix. /// If is not positive definite. - public DenseCholesky(DenseMatrix matrix) + public static DenseCholesky Create(DenseMatrix matrix) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare); } // Create a new matrix for the Cholesky factor, then perform factorization (while overwriting). - var factor = (DenseMatrix)matrix.Clone(); + var factor = (DenseMatrix) matrix.Clone(); Control.LinearAlgebraProvider.CholeskyFactor(factor.Values, factor.RowCount); - Factor = factor; + return new DenseCholesky(factor); + } + + DenseCholesky(Matrix factor) + : base(factor) + { } /// @@ -77,18 +77,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// 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"); - } - - // Check for proper dimensions. if (result.RowCount != input.RowCount) { throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); @@ -117,10 +105,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization } // Copy the contents of input to result. - Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length * Constants.SizeOfFloat); + Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length*Constants.SizeOfFloat); // Cholesky solve by overwriting result. - var dfactor = (DenseMatrix)Factor; + var dfactor = (DenseMatrix) Factor; Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Values, dfactor.RowCount, dresult.Values, dresult.ColumnCount); } @@ -131,18 +119,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// The left hand side , x. public override void Solve(Vector input, Vector result) { - // Check for proper arguments. - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - - // Check for proper dimensions. if (input.Count != result.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); @@ -166,10 +142,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization } // Copy the contents of input to result. - Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length * Constants.SizeOfFloat); + Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length*Constants.SizeOfFloat); // Cholesky solve by overwriting result. - var dfactor = (DenseMatrix)Factor; + var dfactor = (DenseMatrix) Factor; Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Values, dfactor.RowCount, dresult.Values, 1); } } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/UserCholesky.cs b/src/Numerics/LinearAlgebra/Single/Factorization/UserCholesky.cs index f6576895..7a6f09bb 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/UserCholesky.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/UserCholesky.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2010 Math.NET +// 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 @@ -28,9 +28,9 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System; using MathNet.Numerics.Properties; using MathNet.Numerics.Threading; -using System; namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { @@ -43,7 +43,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.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 class UserCholesky : Cholesky + public sealed class UserCholesky : Cholesky { /// /// Initializes a new instance of the class. This object will compute the @@ -53,55 +53,57 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// If is null. /// If is not a square matrix. /// If is not positive definite. - public UserCholesky(Matrix matrix) + public static UserCholesky Create(Matrix matrix) { - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare); } // Create a new matrix for the Cholesky factor, then perform factorization (while overwriting). - Factor = matrix.Clone(); - var tmpColumn = new float[Factor.RowCount]; + var factor = matrix.Clone(); + var tmpColumn = new float[factor.RowCount]; // Main loop - along the diagonal - for (var ij = 0; ij < Factor.RowCount; ij++) + for (var ij = 0; ij < factor.RowCount; ij++) { // "Pivot" element - var tmpVal = Factor.At(ij, ij); + var tmpVal = factor.At(ij, ij); if (tmpVal > 0.0) { - tmpVal = (float)Math.Sqrt(tmpVal); - Factor.At(ij, ij, tmpVal); + tmpVal = (float) Math.Sqrt(tmpVal); + factor.At(ij, ij, tmpVal); tmpColumn[ij] = tmpVal; // Calculate multipliers and copy to local column // Current column, below the diagonal - for (var i = ij + 1; i < Factor.RowCount; i++) + for (var i = ij + 1; i < factor.RowCount; i++) { - Factor.At(i, ij, Factor.At(i, ij) / tmpVal); - tmpColumn[i] = Factor.At(i, ij); + factor.At(i, ij, factor.At(i, ij)/tmpVal); + tmpColumn[i] = factor.At(i, ij); } // Remaining columns, below the diagonal - DoCholeskyStep(Factor, Factor.RowCount, ij + 1, Factor.RowCount, tmpColumn, Control.NumberOfParallelWorkerThreads); + DoCholeskyStep(factor, factor.RowCount, ij + 1, factor.RowCount, tmpColumn, Control.NumberOfParallelWorkerThreads); } else { throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite); } - for (var i = ij + 1; i < Factor.RowCount; i++) + for (var i = ij + 1; i < factor.RowCount; i++) { - Factor.At(ij, i, 0.0f); + factor.At(ij, i, 0.0f); } } + + return new UserCholesky(factor); + } + + UserCholesky(Matrix factor) + : base(factor) + { } /// @@ -113,14 +115,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// Total columns /// Multipliers calculated previously /// Number of available processors - private static void DoCholeskyStep(Matrix data, int rowDim, int firstCol, int colLimit, float[] multipliers, int availableCores) + static void DoCholeskyStep(Matrix data, int rowDim, int firstCol, int colLimit, float[] multipliers, int availableCores) { var tmpColCount = colLimit - firstCol; if ((availableCores > 1) && (tmpColCount > 200)) { - var tmpSplit = firstCol + (tmpColCount / 3); - var tmpCores = availableCores / 2; + var tmpSplit = firstCol + (tmpColCount/3); + var tmpCores = availableCores/2; CommonParallel.Invoke( () => DoCholeskyStep(data, rowDim, firstCol, tmpSplit, multipliers, tmpCores), @@ -133,7 +135,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization var tmpVal = multipliers[j]; for (var i = j; i < rowDim; i++) { - data.At(i, j, data.At(i, j) - (multipliers[i] * tmpVal)); + data.At(i, j, data.At(i, j) - (multipliers[i]*tmpVal)); } } } @@ -146,17 +148,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// The left hand side , X. public override void Solve(Matrix input, Matrix result) { - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - - // Check for proper dimensions. if (result.RowCount != input.RowCount) { throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); @@ -184,10 +175,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization sum = result.At(i, c); for (var k = i - 1; k >= 0; k--) { - sum -= Factor.At(i, k) * result.At(k, c); + sum -= Factor.At(i, k)*result.At(k, c); } - result.At(i, c, sum / Factor.At(i, i)); + result.At(i, c, sum/Factor.At(i, i)); } // Solve L'*X = Y; @@ -196,10 +187,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization sum = result.At(i, c); for (var k = i + 1; k < order; k++) { - sum -= Factor.At(k, i) * result.At(k, c); + sum -= Factor.At(k, i)*result.At(k, c); } - result.At(i, c, sum / Factor.At(i, i)); + result.At(i, c, sum/Factor.At(i, i)); } } } @@ -211,18 +202,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization /// The left hand side , x. public override void Solve(Vector input, Vector result) { - // Check for proper arguments. - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - - // Check for proper dimensions. if (input.Count != result.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); @@ -243,10 +222,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization sum = result[i]; for (var k = i - 1; k >= 0; k--) { - sum -= Factor.At(i, k) * result[k]; + sum -= Factor.At(i, k)*result[k]; } - result[i] = sum / Factor.At(i, i); + result[i] = sum/Factor.At(i, i); } // Solve L'*X = Y; @@ -255,10 +234,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization sum = result[i]; for (var k = i + 1; k < order; k++) { - sum -= Factor.At(k, i) * result[k]; + sum -= Factor.At(k, i)*result[k]; } - result[i] = sum / Factor.At(i, i); + result[i] = sum/Factor.At(i, i); } } } diff --git a/src/Numerics/LinearAlgebra/Single/Matrix.cs b/src/Numerics/LinearAlgebra/Single/Matrix.cs index 5486e14f..911220c4 100644 --- a/src/Numerics/LinearAlgebra/Single/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Single/Matrix.cs @@ -461,7 +461,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single public override Cholesky Cholesky() { - return new UserCholesky(this); + return UserCholesky.Create(this); } public override LU LU()