diff --git a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs index 8bf94015..d42ac4cf 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs @@ -210,14 +210,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra int aRows, int aColumns, T[] b, int bRows, int bColumns, T beta, T[] c); /// - /// Computes the LU factorization of A. + /// Computes the LUP factorization of A. P*A = L*U. /// - /// An m by n matrix. The matrix is overwritten with the - /// the LU factorization On exit. - /// On exit, it contains the pivot indices. The size - /// of the array must be min(m,n). + /// An by matrix. The matrix is overwritten with the + /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0 + /// for the L factor). The upper triangular factor U is stored on and above the diagonal of . + /// The order of the square matrix . + /// On exit, it contains the pivot indices. The size of the array must be . /// This is equivalent to the GETRF LAPACK routine. - void LUFactor(T[] a, int[] ipiv); + void LUFactor(T[] a, int aOrder, int[] ipiv); /// /// Computes the inverse of matrix using LU factorization. diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs index 260d38d1..a49eea9a 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs @@ -688,7 +688,16 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } } - public void LUFactor(double[] a, int[] ipiv) + /// + /// Computes the LUP factorization of A. P*A = L*U. + /// + /// An by matrix. The matrix is overwritten with the + /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0 + /// for the L factor). The upper triangular factor U is stored on and above the diagonal of . + /// The order of the square matrix . + /// On exit, it contains the pivot indices. The size of the array must be . + /// This is equivalent to the GETRF LAPACK routine. + public void LUFactor(double[] a, int aOrder, int[] ipiv) { throw new NotImplementedException(); } diff --git a/src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include b/src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include index f1a2c689..9ee36977 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include +++ b/src/Numerics/Algorithms/LinearAlgebra/NativeAlgebraProvider.include @@ -333,16 +333,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> SafeNativeMethods.d_matrix_multiply(transposeA, transposeB, m, n, k, alpha, a, b, beta, c); } - + /// - /// Computes the LU factorization of A. + /// Computes the LUP factorization of A. P*A = L*U. /// - /// An m by n matrix. The matrix is overwritten with the - /// the LU factorization On exit. - /// On exit, it contains the pivot indices. The size - /// of the array must be min(m,n). + /// An by matrix. The matrix is overwritten with the + /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0 + /// for the L factor). The upper triangular factor U is stored on and above the diagonal of . + /// The order of the square matrix . + /// On exit, it contains the pivot indices. The size of the array must be . /// This is equivalent to the GETRF LAPACK routine. - public void LUFactor(double[] a, int[] ipiv) + public void LUFactor(double[] a, int aOrder, int[] ipiv) { throw new NotImplementedException(); } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs b/src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs index 02f8839b..5557e403 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/Cholesky.cs @@ -110,7 +110,18 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// /// The right hand side , B. /// The left hand side , X. - public abstract Matrix Solve(Matrix input); + 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; + } /// /// Solves a system of linear equations, AX = B, with A Cholesky factorized. @@ -124,7 +135,18 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization /// /// The right hand side vector, b. /// The left hand side , x. - public abstract Vector Solve(Vector input); + 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; + } /// /// Solves a system of linear equations, Ax = b, with A Cholesky factorized. diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs index 67ad3bb5..4d55a7a5 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs @@ -70,24 +70,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization mFactor = factor; } - /// - /// Solves a system of linear equations, AX = B, with A Cholesky factorized. - /// - /// The right hand side , B. - /// The left hand side , X. - public override Matrix Solve(Matrix input) - { - // Check for proper arguments. - if (input == null) - { - throw new ArgumentNullException("input"); - } - - var X = new DenseMatrix(input.RowCount, input.ColumnCount); - Solve(input, X); - return X; - } - /// /// Solves a system of linear equations, AX = B, with A Cholesky factorized. /// @@ -142,24 +124,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization Control.LinearAlgebraProvider.CholeskySolveFactored(dfactor.Data, dfactor.RowCount, dresult.Data, dresult.RowCount, dresult.ColumnCount); } - /// - /// Solves a system of linear equations, Ax = b, with A Cholesky factorized. - /// - /// The right hand side vector, b. - /// The left hand side , x. - public override Vector Solve(Vector input) - { - // Check for proper arguments. - if (input == null) - { - throw new ArgumentNullException("input"); - } - - var x = new DenseVector(input.Count); - Solve(input, x); - return x; - } - /// /// Solves a system of linear equations, Ax = b, with A Cholesky factorized. /// diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs new file mode 100644 index 00000000..deea47c8 --- /dev/null +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs @@ -0,0 +1,179 @@ +// +// 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 an LU factorization. + /// For a matrix A, the LU factorization is a pair of lower triangular matrix L and + /// upper triangular matrix U so that A = L*U. + /// + /// + /// The computation of the LU factorization is done at construction time. + /// + public class DenseLU : LU + { + /// + /// Initializes a new instance of the class. This object will compute the + /// LU factorization when the constructor is called and cache it's factorization. + /// + /// The matrix to factor. + /// If is null. + /// If is not a square matrix. + public DenseLU(DenseMatrix matrix) + { + if (matrix == null) + { + throw new ArgumentNullException("matrix"); + } + + if (matrix.RowCount != matrix.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSquare); + } + + // Create an array for the pivot indices. + mPivots = 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; + } + + /// + /// Solves a system of linear equations, AX = B, with A LU 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"); + } + + // Check for proper dimensions. + if (result.RowCount != input.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSameRowDimension); + } + + if (result.ColumnCount != input.ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension); + } + + if (input.RowCount != mFactors.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + var dinput = input as DenseMatrix; + if (dinput == null) + { + throw new NotImplementedException("Can only do LU factorization for dense matrices at the moment."); + } + + var dresult = result as DenseMatrix; + if (dresult == null) + { + throw new NotImplementedException("Can only do LU factorization for dense matrices at the moment."); + } + + // Copy the contents of input to result. + Buffer.BlockCopy(dinput.Data, 0, dresult.Data, 0, dinput.Data.Length * Constants.SizeOfDouble); + + // LU solve by overwriting result. + var dfactors = mFactors as DenseMatrix; + 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. + /// + /// The right hand side vector, b. + /// 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); + } + + if (input.Count != mFactors.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + var dinput = input as DenseVector; + if (dinput == null) + { + throw new NotImplementedException("Can only do LU factorization for dense vectors at the moment."); + } + + var dresult = result as DenseVector; + if (dresult == null) + { + throw new NotImplementedException("Can only do LU factorization for dense vectors at the moment."); + } + + // Copy the contents of input to result. + Buffer.BlockCopy(dinput.Data, 0, dresult.Data, 0, dinput.Data.Length * Constants.SizeOfDouble); + + // LU solve by overwriting result. + var dfactors = mFactors 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/LU.cs b/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs new file mode 100644 index 00000000..de7a25c2 --- /dev/null +++ b/src/Numerics/LinearAlgebra/Double/Factorization/LU.cs @@ -0,0 +1,171 @@ +// +// 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 an LU factorization. + /// For a matrix A, the LU factorization is a pair of lower triangular matrix L and + /// upper triangular matrix U so that A = L*U. + /// In the Math.Net implementation we also store a set of pivot elements for increased + /// numerical stability. The pivot elements encode a permutation matrix P such that P*A = L*U. + /// + /// + /// The computation of the LU factorization is done at construction time. + /// + public abstract class LU + { + /// + /// Stores both the L and U factors in the same matrix.. + /// + protected Matrix mFactors; + + /// + /// Stores the pivot indices of the LU factorization. + /// + protected int[] mPivots; + + /// + /// Internal method which routes the call to perform the LU factorization to the appropriate class. + /// + /// The matrix to factor. + /// An LU factorization object. + internal static LU Create(Matrix matrix) + { + var dense = matrix as DenseMatrix; + if (dense != null) + { + return new DenseLU(dense); + } + + throw new NotImplementedException(); + } + + /// + /// Returns the lower triangular factor. + /// + public virtual Matrix L + { + get + { + Matrix result = mFactors.GetLowerTriangle(); + for (int i = 0; i < result.RowCount; i++) + { + result.At(i, i, 1); + } + return result; + } + } + + /// + /// Returns the upper triangular factor. + /// + public virtual Matrix U + { + get { return mFactors.GetUpperTriangle(); } + } + + /// + /// 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++) + { + if (mPivots[j] != j) + { + det = -det * mFactors.At(j, j); + } + else + { + det *= mFactors.At(j, j); + } + } + return det; + } + } + + /// + /// Solves a system of linear equations, AX = B, with A LU 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 X = input.CreateMatrix(input.RowCount, input.ColumnCount); + Solve(input, X); + return X; + } + + /// + /// Solves a system of linear equations, AX = B, with A LU 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 LU 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(input.Count); + Solve(input, x); + return x; + } + + /// + /// Solves a system of linear equations, Ax = b, with A LU 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/Matrix.cs b/src/Numerics/LinearAlgebra/Double/Matrix.cs index d36cd2fa..c5da5b90 100644 --- a/src/Numerics/LinearAlgebra/Double/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Double/Matrix.cs @@ -33,6 +33,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double using System; using System.Text; using Properties; + using MathNet.Numerics.Threading; /// /// Defines the base class for Matrix classes. @@ -456,6 +457,108 @@ namespace MathNet.Numerics.LinearAlgebra.Double } } + /// + /// Returns a new matrix containing the lower triangle of this matrix. + /// + /// The lower triangle of this matrix. + public virtual Matrix GetLowerTriangle() + { + Matrix ret = CreateMatrix(RowCount, ColumnCount); + CommonParallel.For(0, ColumnCount, j => + { + for (int i = j; i < RowCount; i++) + { + ret.At(i, j, At(i, j)); + } + }); + return ret; + } + + /// + /// Puts the lower triangle of this matrix into the result matrix. + /// + /// Where to store the lower triangle. + /// If is . + /// If the result matrix's dimensions are not the same as this matrix. + public virtual void GetLowerTriangle(Matrix result) + { + if (result == null) + { + throw new ArgumentNullException("result"); + } + + if (result.RowCount != RowCount || result.ColumnCount != ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions, "result"); + } + + CommonParallel.For(0, ColumnCount, j => + { + for (int i = 0; i < RowCount; i++) + { + if (i >= j) + { + result.At(i, j, At(i, j)); + } + else + { + result.At(i, j, 0); + } + } + }); + } + + /// + /// Returns a new matrix containing the upper triangle of this matrix. + /// + /// The upper triangle of this matrix. + public virtual Matrix GetUpperTriangle() + { + Matrix ret = CreateMatrix(RowCount, ColumnCount); + CommonParallel.For(0, ColumnCount, j => + { + for (int i = 0; i <= j; i++) + { + ret.At(i, j, At(i, j)); + } + }); + return ret; + } + + /// + /// Puts the upper triangle of this matrix into the result matrix. + /// + /// Where to store the lower triangle. + /// If is . + /// If the result matrix's dimensions are not the same as this matrix. + public virtual void GetUpperTriangle(Matrix result) + { + if (result == null) + { + throw new ArgumentNullException("result"); + } + + if (result.RowCount != RowCount || result.ColumnCount != ColumnCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions, "result"); + } + + CommonParallel.For(0, ColumnCount, j => + { + for (int i = 0; i < RowCount; i++) + { + if (i <= j) + { + result.At(i, j, At(i, j)); + } + else + { + result.At(i, j, 0); + } + } + }); + } + #region Implemented Interfaces #if !SILVERLIGHT diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index ccac2022..fe5ca4f5 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -145,6 +145,8 @@ + +